Gluino-mediated electroweak penguin with flavor-violating trilinear couplings

In light of a discrepancy of the direct $CP$ violation in $K\to\pi\pi$ decays, $\varepsilon'/\varepsilon_K$, we investigate gluino contributions to the electroweak penguin, where flavor violations are induced by squark trilinear couplings. Top-Yukawa contributions to $\Delta S = 2$ observables are taken into account, and vacuum stability conditions are evaluated in detail. It is found that this scenario can explain the discrepancy of $\varepsilon'/\varepsilon_K$ for the squark mass smaller than 5.6 TeV. We also show that the gluino contributions can amplify $\mathcal{B}(K \to \pi \nu \overline{\nu})$, $\mathcal{B}(K_S \to \mu^+ \mu^-)_{\rm eff}$ and $\Delta A_{\rm CP}(b\to s\gamma)$. Such large effects could be measured in future experiments.


Introduction
Flavor-changing neutral current (FCNC) processes are sensitive probes for new physics beyond the Standard Model (SM), since in the SM there are no FCNC processes at the tree level, and they are suppressed further by the Glashow-Iliopoulos-Maiani (GIM) mechanism. One of the FCNC observables, the CP -violating ratio ε /ε K in neutral kaon decays into two pions, has attracted attentions recently because of a discrepancy between the experimental data and the theoretical predictions based on the first lattice calculation of the hadronic parameters B (1/2) 6 and B (3/2) 8 by the RBC-UKQCD collaboration [1][2][3][4]. 1 The next-toleading order (NLO) prediction for ε /ε K has been calculated in ref. [9], and it has been confirmed by an improved calculation in ref. [10]. The latter result is given by ε /ε K SM = (1.06 ± 5.07) × 10 −4 , (1.1)

JHEP04(2018)019
which deviates from the experimental data [11][12][13][14] Re ε /ε K exp = (16.6 ± 2.3) × 10 −4 , (1.2) at the 2.8 σ level. The theoretical result which is much smaller than the data is supported by analyses in the large-N c dual QCD approach [15,16]. Note that improvements of the lattice calculation and independent confirmations of the result by other lattice collaborations are highly important to establish the presence of new physics in ε /ε K .
In this paper, we study ε /ε K in the minimal supersymmetric standard model (MSSM) with introducing large off-diagonal entries in the trilinear couplings of the down-type squarks to the Higgs boson. The off-diagonal couplings generate gluino contributions to the flavor-changing Z penguin which affects ε /ε K via the I = 2 amplitude. Although such a scenario has been studied in ref. [17], top-Yukawa contributions to ∆F = 2 observables have not been taken into account. In the scenario, ε K receives those contributions from the Z penguin through the renormalization group (RG) running from the new physics scale to the electroweak (EW) scale, and through the matching onto the low-energy FCNC operators at the EW scale [18,19]. They can be comparable in size to ordinary gluino box contributions. Moreover, since the LHC experiment is pushing up the lower bounds on the squark and gluino masses [20,21], the situation changes: larger trilinear couplings are required to explain the ε /ε K discrepancy.
The large off-diagonal trilinear couplings also affect other FCNC observables. We consider constraints on the couplings as well as on other MSSM parameters from the branching ratios of K L → µ + µ − ,B → X s γ andB → X d γ in addition to ε K . Furthermore, such large trilinear couplings can make the EW vacuum unstable. Although the vacuum instability was overlooked in ref. [17], we investigate the vacuum (meta-)stability condition in detail and show that the constraint is significant. In ref. [22], the vacuum condition has been studied in another scenario with large off-diagonal trilinear couplings of the up-type squarks, which bring chargino contributions to the Z penguin. An alternative scenario for the explanation of the ε /ε K discrepancy in the MSSM has been proposed in refs. [23,24].
The discrepancy in ε /ε K requires large CP -violating phases in the off-diagonal trilinear couplings. They also contribute to the branching ratios of K + → π + νν and K L → π 0 νν, the effective branching ratio of K S → µ + µ − [25,26] and the CP asymmetry difference ∆A CP (b → sγ). We investigate SUSY effects on these observables in our scenario, and examine if the effects can be observed at current and/or near-future experiments. This paper is organized as follows. In section 2 we summarize the effective Lagrangian together with the RG equations and the one-loop matching conditions that are relevant to our analysis. Top-Yukawa contributions are also explained. In section 3 we present the gluino contributions associated with the Z penguin. In section 4 we explain how each FCNC observable receive gluino contributions. In section 5 we discuss the constraints from the vacuum stability condition. In section 6 we present our numerical analysis. Our conclusions are drawn in section 7.

JHEP04(2018)019 2 Effective Lagrangian and top-Yukawa contributions
In this paper, we study flavor-changing processes via the gluino one-loop contributions and the Z-boson exchanges. The latter is described by higher dimensional operators in the SM effective field theory (SMEFT), where the gauge invariance is guaranteed. The effective Lagrangian is defined as where the first term in the right-hand side is the SM Lagrangian, and the second one is composed by higher dimensional operators [27]. In particular, those relevant to the ∆F = 1 Z-boson penguin are given by Here, q is the (left-handed) SU(2) quark doublets and d is the (right-handed) down-type quark singlets with quark-flavor indices, i, j, and an SU(2) index, a. The Higgs doublet carries a hypercharge +1/2, and thus, has a vacuum-expectation value (VEV), H = (0, v/ √ 2) T , with v 246 GeV after the EW symmetry breaking (EWSB). The covariant derivative is defined for the Higgs doublet as and On the other hand, ∆F = 2 processes are described by the following four-Fermi operators,

JHEP04(2018)019
where Y t is the top-quark Yukawa coupling. It is noticed that there are no O(α s ) corrections at the one-loop level. The operators also contribute to ∆S = 2 four-quark operators as HQ ] 12 + . . . , (2.14) 12 . In the first leading logarithm approximation, the Wilson coefficients after the RG running from Λ to µ (Λ > µ) are estimated as Irrelevant operator mixings and higher-order corrections during the evolutions are neglected. In particular, C QQ and C (1) QD are generated by C HQ and C HD . After the EWSB, O HQ and O HD are matched to the flavor-changing Z couplings through the expansion, where the terms irrelevant for the matching onto the ∆S = 2 operators are omitted.
The operators also contribute to ∆F = 2 observables through the effective Hamiltonian, (2.18) where the effective operators are with color indices α, β. In this paper, chirality-flipped operators and their Wilson coefficients are denoted with a prime. At the tree level, the SMEFT operators are matched at the weak scale to these operators as [31] [  where N c = 3 is the number of colors. In addition, these low-energy ∆F = 2 operators are generated by the ∆F = 1 ones in the SMEFT through the one-loop matchings at the weak scale [31]. The conditions for C HQ and C HD at the scale µ W are approximated as [18,19] [ These results are gauge-independent. The loop functions are defined as Here, we discarded box contributions which are suppressed by CKM factors or by m 2 c,u /m 2 W in the ∆S = 2 case (see ref. [19]

SUSY contributions
Here, md L(R)i is the left-(right-) handed squark soft mass for the i-th generation, mg is the gluino mass, and T D is the scalar trilinear coupling of the down-type squarks. In this paper, the SUSY Les Houches Accord (SLHA) notation [32,33] is used, and flavor violations are discussed in the basis where the Yukawa matrix of the downtype quark is diagonalized. The Wilson coefficients are set at the SUSY scale. 2 The loop function is defined as In the limit of y, z → x, it becomes Other SUSY contributions are explained in the next section. Note that, in literature, e.g., ref. [34], it has been argued that gluino-mediated contributions to EW penguin are suppressed compared to the other penguins, by assuming that the gluino contributions to the EW penguin are proportional to those to the photon penguin. However, this is not the case in our scenario, where the SU(2)×U(1) symmetry is broken by large scalar trilinear couplings. Such couplings can generate the Z penguin significantly via double-mass insertion contributions, as was pointed out in ref. [35] and explicitly shown in this section.

ε /ε K
The direct CP violation of the K → ππ decays, ε /ε K , includes the SM and SUSY Zpenguin contributions,

JHEP04(2018)019
The latter contribution is approximated to be (cf., ref. [36] where the Wilson coefficients are estimated at the Z-boson mass scale, µ = m Z . By using lattice simulations [2][3][4], B (3/2) 8 (m c ) = 0.76 ± 0.05 is obtained [9,44]. Here, ε K in the denominator is evaluated by the experimental value. The right-handed contribution is amplified by c 2 W /s 2 W 3.33 compared to the left-handed one. Currently, the SM prediction deviates from the experimental result at the 2.8 σ level. In this paper, the discrepancy of ε /ε K is required to be explained within the 1 σ range, where ref. [10] is used for the SM prediction at the NLO level.

ε K
Both the SM and SUSY affect to the indirect CP violation of the neutral kaon system, is composed by gluino box diagrams as well as C HQ and C HD . In our scenario, although the gluino box contributions are sizable, their dominant contributions arise as dimension-ten operators in the SMEFT. In order to include them in our formalism, we separately calculate them in the broken phase, where the Higgs VEV is involved. 4 At the one-loop level, they are obtained as [45] [ at the SUSY scale (µ SUSY ) with generation indices i = j and x r = m 2 dr /m 2 g , where R d ri for r = 1, 2, . . . , 6 is the squark rotation matrix defined in the SLHA notation [32,33]. C 1,2,3 are obtained by flipping the chirality of R d( * ) ri in C 1,2,3 . The loop functions are defined as , (4.10) .
From µ SUSY to the hadronic scale, we solve the RG equations at the NLO level [46] and use the hadronic matrix elements in ref. [47].
Additionally, [C 1 ] ij and [C 5 ] ij receive the top-Yukawa contributions depending on C HQ and C HD as at the Z-boson mass scale. These results are derived as follows: the Wilson coefficients are evolved by solving the RG equations with the beta function (2.14) in the first leading logarithm approximation (2.15), and then, matched onto the low-scale operators at the weak scale (2.24)-(2.26). Also, the one-loop matchings, (2.27) and (2.28), are taken into account to include the additional contributions of C HQ and C HD at the weak scale (see ref. [19]). 5 Equivalently, the same results are reproduced by substituting µ W → µ SUSY in eqs. (2.27) and (2.28). This is because the logarithmic scale dependence of the one-loop matching conditions has the same origin as the one-loop beta functions (see ref. [18]). It is also noticed that, in eq. (4.12), the logarithmic dependence of µ SUSY cancels out because of [C The SM value is estimated to be where the input SM parameters are found in ref. [48] (cf., ref. [49]). Especially, the Wolfenstein parameters are determined by the angle-only fit [50], and |V cb | obtained from inclusive semileptonic B decays (B → X c −ν ) [51] is used. 6 We use lattice results for the ξ 0 parameter [1], which parametrizes the absorptive part of long-distance effects, and refrain from 5 The results are independent of the matching scale µW by including the one-loop matching conditions.
Consequently, the logarithmic function becomes ln(µSUSY/mW ). 6 Recently, there are debates about systematic uncertainties of the exclusive determinations of

JHEP04(2018)019
relying on the experimental result of ε /ε K , because we consider SUSY contributions to ε /ε K . On the other hand, the experimental result is (cf., ref. [14]) Therefore, the SUSY contributions are required to be within the range, at the 2 σ level. 7

K → πνν
The Z-penguin contributions induce the decays, K + → π + νν and K L → π 0 νν. They are expressed as [36,44] 8 , and the charm contribution gives P c (X) = (9.39±0.31)×10 −4 /λ 4 +(0.04± 0.02). In terms of C HQ and C HD , X eff is approximated to be (cf., ref. [18]) Im X eff = 2.12 × 10 −4 + 5.62 × 10 6 GeV 2 Im C H+ , (4.20) where the first terms in the right-hand sides are the SM contributions in each equation, and The Wilson coefficients are estimated at the Z-boson mass scale. The SM predictions are known to be [18] B(K + → π + νν) SM = (8.5 ± 0.5) × 10 −11 , (4.22) while the experimental results are [55,56] These experimental values will be improved in the near future. The NA62 experiment at CERN has already started the physics run and aims to measure B(K + → π + νν) with a precision of 10% relative to the SM prediction [57]. The KOTO experiment at J-PARC aims to measure B(K L → π 0 νν) around the SM sensitivity by 2021 [58,59]. JHEP04(2018)019 The decay rate of K L → µ + µ − , which is a CP -conserving process, is sensitive to a real component of the flavor-changing Z couplings. There are large theoretical uncertainties from a long-distance (LD) contribution. In addition, an unknown sign of A (K L → γγ) conceals a relative sign between the LD and a short-distance (SD) amplitudes. One can, therefore, estimate only the SD branching ratio, which is expressed as [36,60,61] Here, Y eff is approximately given as (cf., ref. [18]) where the first term in the right-hand side is the SM contribution, and The Wilson coefficients are estimated at the Z-boson mass scale.
The SM value is obtained as [18] B It is challenging to extract the SD contribution from the experimental value. An upper bound is estimated as [62] B(K L → µ + µ − ) exp SD < 2.5 × 10 −9 . (4.30) Since the constraint is much weaker than the SM uncertainties, we simply impose a bound, The decay, K S → µ + µ − , proceeds via LD CP -conserving P-wave and SD CP -violating Swave processes. Since the decay rate is dominated by the former, whose uncertainty is large, the sensitivity to the imaginary component of the flavor-changing Z couplings is diminished [62][63][64]. Interestingly, the SD contribution is enhanced through an interference between the K L and K S states in the neutral kaon beam [26]. The effective branching ratio of K S → µ + µ − after including the interference is expressed as (cf., ref. [26]) where a dilution factor D is an initial asymmetry between the numbers of K 0 and K 0 ,

JHEP04(2018)019
In the right-hand side, the branching ratio is approximated to be where the first and second terms in the right-hand side come from the LD and SD contributions, respectively. Here, the Wilson coefficients are estimated at the Z-boson mass scale. On the other hand, the interference contribution is given as The Wilson coefficients are estimated at the Z-boson mass scale. The unknown relative sign between the LD and SD contributions in K L → µ + µ − gives two different predictions of B (K S → µ + µ − ) int , which are expressed by η A , (see refs. [26,65]) . (4.36) Here, scalar operator contributions are discarded in the above formulae: they can be significant especially when tan β is large and m A is small [25]. The SM prediction depends on D and η A , which are determined by experiments. For D = 0, it is obtained as [26,62,63]  The experimental sensitivity is expected to reach B(K S → µ + µ − ) = O(10 −11 ) by the end of the LHCb Run-2, and the Run-3 project is aiming to achieve the sensitivity as precise as the SM level [67].
where the effective operators are defined as where e > 0 and g 3 > 0, and the covariant derivatives for the quark and squark follow the same sign convention as eq. (2.5). At the one-loop level, the gluino contributions are obtained as where x r = m 2 dr /m 2 g , and the loop functions are defined to be Also, C 7γ and C 8g are obtained by flipping the chirality of R d( * ) ri in C 7γ and C 8g , respectively. In the analysis, an approximation formula in ref. [68]

JHEP04(2018)019
for E γ > 1.6 GeV. In the analysis, the theoretical prediction including the SM and SUSY contributions is required to be consistent with the experimental result at the 2σ level. CP violations of b → d i γ are sensitive to the imaginary parts of flavor-violating scalar trilinear couplings. Long-distance effects tend to spoil the sensitivity [74]. This could be resolved by taking a difference of the CP asymmetries [74], where the right-handed contributions are taken into account [75]. The hadronic parameter Λ 78 introduces an uncertainty to the analysis and is estimated to be 12 MeV < Λ 78 < 190 MeV [68]. We take an average value, Λ 78 = 89 MeV, in the analysis. The Wilson coefficients include both the SM and SUSY contributions, which are evaluated at the scale µ b = 2 GeV. The SM prediction is expected to be much suppressed, ∆A CP (b → sγ) SM ≈ 0 [74]. On the other hand, the experimental result is [76] ∆A CP (b → sγ) exp = (5.0 ± 3.9 stat ± 1.5 syst )% (4.53) from the BaBar experiment. The Belle experiment also published a result on ∆A CP (B → K * γ) [77], The asymmetry of the inclusive decay is expected to be comparable to that of the exclusive mode [78]. Both results are consistent with a null asymmetry difference. Since the uncertainties are large, the SUSY parameters will not be constrained in the region of our interest. In future, the uncertainty is projected to achieve 0.37% for ∆A CP (b → sγ) at Belle II with 50 ab −1 [79]. 9

Vacuum stability
The Wilson coefficients in eqs. (3.1)-(3.3) are enhanced by large off-diagonal trilinear couplings, (T D ) i3 and (T D ) 3i (i = 1, 2). Such large trilinear couplings tend to generate dangerous charge and color breaking (CCB) global minima in the scalar potential [80]. Hence, they are limited by the vacuum (meta-)stability condition: the lifetime of the EW vacuum must be longer than the age of the Universe. In this section, we will investigate the vacuum stability conditions of (T D ) i3 and (T D ) 3i . The vacuum decay rate per unit volume is represented by Γ/V = A exp (−S E ), where S E is the Euclidean action of the bounce solution [81]. CosmoTransition 2.0.2 [82] is used to estimate S E at the semiclassical level. The prefactor A cannot be determined unless radiative corrections are taken into account [83,84]. We adopt an order-of-magnitude 9 Although the experimental uncertainty of the direct CP asymmetry ACP(b → sγ) is also projected to be sub-percent level [79], long-distance contributions as well as hadronic uncertainties spoil the SM prediction [74].

JHEP04(2018)019
estimation, A ∼ (100 GeV) 4 . By requiring (Γ/V ) 1/4 to be smaller than the current Hubble parameter, the lifetime of the EW vacuum becomes longer than the age of the Universe. The condition corresponds to S E 400. In this paper, thermal effects and radiative corrections to the vacuum transitions are discarded.
The bounce solution and S E are determined by the scalar potential. The potential relevant for the vacuum decay generated by (T D ) 13 and/or (T D ) 31 is where the coefficients are Let us first consider the vacuum stability condition when only (T D ) 13 is large. The scalar potential is simplified to be When m A ∼ mQ ,1 ∼ mD ,3 , CCB vacua appear around a h d -d L -b R plane. In figure 2, the solid lines show upper bounds on |(T D ) 13 | for tan β = 5, 10, 30, and 50. We assumed m A = mQ ,1 = mD ,3 . It is shown that the upper bounds are proportional to mQ. Also, the results depend on tan β slightly. This is because the scalar potential is stabilized by a quartic coupling y 2 bb 2  When m A is larger than mQ ,1 ∼ mD ,3 , the position of the CCB vacuum approaches to a H-d L -b R plane, where H includes the SM-like Higgs boson, H = h SM + v. In figure 3, the m A dependence of the upper bound is shown. Here, tan β = 5 and 30 are taken. We found that the vacuum stability condition is relaxed for large m A .
In the decoupling limit of the heavy Higgs bosons (m 2 A m 2 Z , α → β −π/2), the scalar potential can be expressed by H,d L , andb R as The upper bounds on |(T D ) 13 | are shown by the dashed lines in figure 2. 10 Again, they are proportional to mQ. In contrast to the case of m A ∼ mQ, the result is almost proportional to tan β. This is understood by cos β associated to (T D ) 13 . A fitting formula of the vacuum stability condition in the large m A limit with mQ ,1 = mD ,3 ≡ mQ is derived as where the phase of (T D ) i3 is taken into account. This formula works well for mQ > 1 TeV. Let us next turn on (T D ) 23 in addition to (T D ) 13 . The scalar trilinear term becomes Here, (T D ) 13,23 are taken to be real by rephasing the scalar fields. By mixingd L ands L , one can obtain whered L =d L cos θ −s L sin θ ands L =d L sin θ +s L cos θ with tan θ = (T D ) 23 / (T D ) 13 . When m 2Q ,1 = m 2Q ,2 ≡ m 2Q , the scalar potential ofd L is obtained from that ofd L by as well asd L →d L . Therefore, the vacuum stability condition (5.7) is extended to be where the phases of (T D ) 13,23 are taken into account appropriately. The formula is valid when mQ ≡ mQ ,1 = mQ ,2 = mD ,3 > 1 TeV and m A is decoupled. 11 When only (T D ) 31 is large, the potential becomes By repeating the above procedure, one can obtain quantitatively the same fitting formula for (T D ) 3i as eq. (5.10), where mQ ≡ mQ ,3 = mD ,1 = mD ,2 > 1 TeV and m A is decoupled. 10 In this scalar potential, the SM-like Higgs boson is lighter than 125 GeV. The vacuum stability condition can be evaluated naively by adding top-stop radiative corrections, H sin 4 βH 4 /8, [85-88] to eq. (5.6) in order to achieve the 125 GeV SM-like Higgs boson at the EW vacuum. We found that eq. (5.7) is barely changed. Dedicated studies are needed to fully include the radiative corrections (see ref. [84]). 11 We have validated the formula (5.10) explicitly by analyzing the bounce action of the scalar potential of H,dL,sL, andbR.

JHEP04(2018)019 6 Numerical analysis
In this section, we study gluino contributions to ε /ε K via the Z penguin. They are enhanced by large scalar trilinear couplings as shown in section 3. Since (T D ) 13,23,31,32 are complex variables, there are 8 degrees of freedom. For simplicity, we restrict the parameter space such that two of (T D ) 13,23,31,32 are real. When (T D ) 23,32 are real, we checked that wide parameter regions to explain the discrepancy of ε /ε K are tightly excluded by B (B → X d,s γ). Therefore, we consider the cases when (T D ) 13,31 are real. The scalar trilinear coupling are parameterized as where α i , β i and γ i are real parameters. Then, one obtains (see section 3) 3) The L variables contribute to the left-handed Wilson coefficients, and the R variables to the right-handed ones. In order to evaluate the observables, we scan the whole parameter region of α i , β i , and γ i where the vacuum stability conditions are satisfied. 12 When β L γ L > 0 and β R γ R > 0, the SUSY contribution to ε /ε K is maximized, because the left-handed contribution, C HQ , constructively interferes with the right-handed one, C HD . In this case, B(K L → π 0 νν) cannot exceed the SM prediction, because positive β L γ L and β R γ R tends to decrease the branching ratio, as can be seen from eq. (4.20). We consider this case in section 6.1. In contrast, ε /ε K cannot be accommodated with the result (4.3) for β L γ L < 0 and β R γ R < 0. When either β L γ L or β R γ R is negative, the discrepancy of ε /ε K can also be explained. Because the right-handed contribution to ε /ε K is larger than the left-handed one, β R γ R > 0 is favored to amplify ε /ε K . At the same time, B(K L → π 0 νν) can be enhanced and may exceed the SM value. Hence, we consider the case when β L γ L < 0 and β R γ R > 0 in section 6.2.
Before proceeding to the analysis, let us summarize assumptions on model parameters. Since the vacuum stability condition is relaxed by large m A , the heavy Higgs bosons are supposed to be decoupled. The squark masses are set to be degenerate, mQ ≡ mQ ,1 = mQ ,2 = mQ ,3 = mD ,1 = mD ,2 = mD ,3 , for simplicity. The Higgsino mass parameter is also equal to mQ, though dependences of the observables on it are weak. We take tan β = 5, though the following results are insensitive to the choice, because the observables as well as the vacuum stability condition depend on it dominantly in a combination of T D cos β.
In figure 4, the maximal values of the SUSY contributions to ε /ε K are shown for β L γ L > 0 and β R γ R > 0 as a function of mQ. There is a peak structure for each line. In smaller squark mass regions, the maximal value is determined by B(B → X d γ). Defining the squark JHEP04(2018)019 Figure 4. The maximal gluino contributions to ε /ε K as a function of mQ. The parameters are γ R /β R = γ L /β L = 1 and mg/mQ = 1 on the black line. In the left plot, γ R /β R = γ L /β L = 0.6, 0.8, 1.2 with mg/mQ = 1 from left to right of the red lines. In the right plot, mg/mQ = 1.8, 1.4, 0.8 with γ R /β R = γ L /β L = 1 from left to right of the green lines.
mixing parameter, δ D = (T D ) ij v cos β/m 2Q , the SUSY contributions to ε /ε K depend on it where mg ∼ mQ is supposed. Thus, the maximal value of ε /ε K increases as mQ becomes larger. In larger squark mass regions, the maximal value is determined by ε K , B(B → X s γ) and the vacuum stability condition as well as B(B → X d γ). In particular, the gluino box contribution to ε K depends on δ D as ∼ δ 4 D /m 2Q , whereas the SUSY contributions via C HQ and C HD are not suppressed by mQ, i.e., behaves as ∼ λ t δ 2 D /m 2 Z . When mQ is small, the latter contribution can be canceled enough by the former one. However, as mQ increases, the cancellation becomes weaker in the parameter region allowed by the other constraints. Hence, the bounds on the trilinear couplings become severer to satisfy the constraint of ε K . Consequently, the maximal value of ε /ε K decreases.
In the figures, γ i /β i or mg/mQ is also varied. On the black line, γ R /β R = γ L /β L = 1 and mg/mQ = 1 are chosen. In the left plot, γ R /β R = γ L /β L = 0.6, 0.8, 1.2 with mg/mQ = 1 from left to right of the red lines. On the other hand, mg/mQ = 1.8, 1.4, 0.8 with γ R /β R = γ L /β L = 1 from left to right of the green lines in the right plot. The maximum value increases when γ i /β i is small and mg/mQ is large. Also, it is found that the current discrepancy of ε /ε K can be explained if the squark mass is smaller than 5.6 TeV.
We study other observables with keeping the SUSY contribution to ε /ε K sizable for β L γ L < 0 and β R γ R > 0. The SUSY parameters are determined to achieve (ε /ε K ) SUSY = 10.0 × 10 −4 , where the current discrepancy between the experimental and SM values is explained at the 1σ level.  In figure 5, B(K L → π 0 νν) is maximized for given mQ. One finds a peak structure for each line. On the left side of the peak, the parameters are constrained by B(B → X d γ). If the soft masses are too small, ε /ε K cannot be large sufficiently. On the right side, the constraints from ε K and B(B → X s γ) become relevant. When SUSY particles are very heavy, the SUSY contribution to ε K via C HQ and C HD cannot be canceled enough by that via the gluino box contribution in the parameter region allowed by the other constraints.
One can see that B(K L → π 0 νν) can be larger than the SM value. This result is contrasted with the case when β L γ L > 0 and β R γ R > 0.
In the figures, γ i /β i or mg/mQ is also varied. On the black line, γ R /β R = −γ L /β L = 1 and mg/mQ = 1 are chosen. In the left plot, γ R /β R = −γ L /β L = 0.6, 0.8, 1.2 with mg/mQ = 1 from left to right of the red lines. On the other hand, mg/mQ = 1.8, 1.4, 0.8 with γ R /β R = −γ L /β L = 1 from left to right of the green lines in the right plot. In both plots, the peak positions depend on the setup. The maximum value increases when |γ i /β i | is small and/or mg/mQ is large. It is found that B(K L → π 0 νν) can be about 1.5 times larger than the SM prediction. Such a branching ratio could be discovered in future KOTO experiment.
Next, B(K + → π + νν) is maximized for given mQ in figure 6. The branching ratio depends on C HQ and C HD similarly to the case of B(K L → π 0 νν). Hence, it can be larger than the SM prediction when either β L γ L or β R γ R is negative. The real component of C HQ and C HD contributes to the ratio, which is different from the case of B(K L → π 0 νν) and ε /ε K . Consequently, the peak structure in figure 5 disappears. The maximal value tends to decrease as mQ increases. They are enhanced when |γ i /β i | is small and mg/mQ is large. The maximal value can be about 1.6-1.7 times larger than the SM prediction. The deviation could be measured in the current NA62 experiment. ��(� + →π + νν)/�� Figure 6. The maximum value of B(K + → π + νν) normalized by the SM prediction as a function of mQ. Here, (ε /ε K ) SUSY = 10.0 × 10 −4 is fixed. The parameters are γ R /β R = −γ L /β L = 1 and mg/mQ = 1 on the black line. In the left plot, γ R /β R = −γ L /β L = 0.6, 0.8, 1.2 with mg/mQ = 1 from left to right of the red lines. In the right plot, mg/mQ = 1.8, 1.4, 0.8 with γ R /β R = −γ L /β L = 1 from left to right of the green lines.    Let us also mention about the CP -violating observable, ∆A CP (b → sγ). In the analysis, since the CP -violating phases arise in (T D ) 23 and (T D ) 32 , the asymmetry can be sizable. In figure 7, the maximum value of ∆A CP (b → sγ) is shown as a function of mQ.
Here, (ε /ε K ) SUSY = 10.0 × 10 −4 is fixed. On the black line, γ R /β R = −γ L /β L = 1 and mg/mQ = 1 are chosen. In the left plot, the trilinear coupling is varied as γ R /β R = −γ L /β L = 0.6, 0.8, 1.2 with mg/mQ = 1 from left to right of the red lines. In the right plot, the gluino mass is set as mg/mQ = 1.8, 1.4, 0.8 with γ R /β R = −γ L /β L = 1 from left to right of the green lines. It is found that the asymmetry is enhanced especially when |γ i /β i | is small, because smaller ratios lead to larger (T D ) 23 and (T D ) 32 to achieve (ε /ε K ) SUSY = 10.0 × 10 −4 . Also, when mg/mQ is small, the asymmetry becomes large. The CP asymmetry can be as large as 14% for γ R /β R = −γ L /β L = 0.6. We also find that ∆A CP (b → sγ) is likely to be positive when it is enhanced in our scenario. Such an asymmetry seems to be large enough to be measured at Belle II with 50 ab −1 . 13 Finally, we study the SUSY contribution to K S → µ + µ − as a function of mQ. They are enhanced when the sign of the left-handed contribution is opposite to that of the righthanded one. Such a setup is realized in this subsection. In figure 8, the effective branching ratio of K S → µ + µ − is shown. Here, the dilution factor D = 1 and the relative sign η A = −1 are chosen as a reference case. 14 Since the interference term is almost independent of a real component of C H− in the parameter regions of our interest, B (K S → µ + µ − ) eff is determined once (ε /ε K ) SUSY and B(K L → π 0 νν) are given. Therefore, in figure 8, we 13 Although a part of the parameter regions seems to be constrained by the current experimental result (4.53), the theoretical uncertainty is large, and thus, we have not employed this limit.
14 In the case of D = 0, we find that the branching ratio B(KS → µ + µ − ) in eq. (4.34) is not deviated from the SM value (4.37) sizably.

JHEP04(2018)019
take the same α i , β i and γ i as those in figure 5, which maximize B(K L → π 0 νν). It is found that B (K S → µ + µ − ) eff is enhanced especially when |γ i /β i | is small. The effective branching ratio can be 1.9 × 10 −11 , which is larger than the SM prediction (4.38). Such a branching ratio might be measured by the end of the LHCb Run-2, and it is large enough to be detected at the LHCb Run-3 [67].

Conclusions
In this paper, we studied CP violations in the neutral kaon decay in the MSSM scenario where non-minimal flavor mixings and CP -violating phases reside in the trilinear scalar couplings of the down-type squarks. We calculated SUSY contributions that are induced by one-loop diagrams involving gluino and squarks, and evaluated their effects on flavor observables. We took the top-Yukawa contributions to ∆S = 2 observables into account. Considering constraints from the vacuum stability and the measurements of ε K , B(K L → µ + µ − ), B(B → X s γ) and B(B → X d γ), we searched for the allowed parameter regions of the trilinear coupling parameters and investigated possible effects on ε /ε K , B(K L → π 0 νν), B(K + → π + νν), B(K S → µ + µ − ) eff and ∆A CP (b → s γ).
We found that the difference between the measured value and the SM prediction of ε /ε K can be explained by the gluino-mediated Z-penguin contribution to the s → d transition amplitude for the squark mass smaller than 5.6 TeV. In addition, B(K L → π 0 νν) and B(K + → π + νν) can be enhanced by about 50 % and 70 % of the SM values, respectively. It is also shown that B(K S → µ + µ − ) eff and ∆A CP (b → s γ) are significantly enhanced.
The deviations from the SM predictions of these observables can be probed in nearfuture experiments such as KOTO, NA62, LHCb and Belle II. Since the pattern of the deviations is closely related to the structure of the trilinear coupling matrix in the model, the measurements would provide us with important clues to explore flavor structures in physics beyond the SM.