$\varepsilon'/\varepsilon$ in 331 Models

Motivated by the recent findings that the ratio $\varepsilon'/\varepsilon$ in the Standard Model (SM) appears to be significantly below the data we investigate whether the necessary enhancement of this ratio can be obtained in 331 models, in which new physics (NP) contributions to $\varepsilon'/\varepsilon$ and other flavour observables are dominated by tree-level exchanges of a $Z^\prime$. NP contributions to $\varepsilon'/\varepsilon$ in these models are governed by the electroweak operator $Q_8$ for which the hadronic matrix element is rather well known so that our analysis of NP contributions is subject to much smaller uncertainties than within the SM. We consider seven 331 models selected in our earlier analysis on the basis of electroweak precision data. Imposing the constraints from $\Delta F=2$ transitions we find that only three of these models can provide a significant positive shift in $\varepsilon'/\varepsilon$ up to $6\times 10^{-4}$ for $M_{Z^\prime}=3$ TeV. Two of them allow simultaneously a supression of ${\cal B}(B_{s}\to \mu^+\mu^-)$ by $20\%$, bringing the theory closer to the data without any significant impact on the Wilson coefficient $C_9$. The third one provides also the shift $\Delta C_9=-0.6$, softening the anomalies in $B\to K^*\mu^+\mu^-$, without any significant impact on $B_{s}\to \mu^+\mu^-$. NP effects in rare $K$ decays and in $B\to K(K^*)\nu\bar\nu$ turn out to be small. The flavour structure of 331 models implies that even for $M_{Z^\prime}=30$ Tev a shift of $\varepsilon'/\varepsilon$ up to $8\times 10^{-4}$ and a significant shift in $\varepsilon_K$ can be obtained, while the effects in other flavour observables are small.


Introduction
The recent findings that the ratio ε /ε predicted by the Standard Model (SM) appears to be significantly below the experimental data [1] poses a natural question what kind of new physics (NP) could be responsible for this new anomaly. In the present paper we will address this question in 331 models based on the gauge group SU (3) C × SU (3) L × U (1) X [2,3]. These models display several appealing features. The first one is that the requirement of asymptotic freedom of QCD together with that of anomaly cancelation constrains the number of generations to be necessarily equal to the number of colours, providing an explanation for the existence of three generations. Moreover, under the action of SU (3) L two quark generations should transform as triplets, one as an antitriplet. Adopting the choice that the third generation is the one transforming as an antitriplet, this different treatment could be at the origin of the large top mass. It should be recalled that some of the generators of the group are connected by the relation Q = T 3 + βT 8 + X where Q is the electric charge, T 3 and T 8 are two of the SU (3) generators and X the generator of U (1) X . β is a parameter that defines a specific variant of the model.
Several new particles are present in these models, their features depending on the chosen variant. However, in all the variants a new neutral gauge boson Z exists that can mediate tree level flavour changing neutral currents (FCNC) in the quark sector.
In the framework of 331 models, the ratio ε /ε has been studied in our earlier analysis [4]. Here we update it and improve it.
Recent analyses addressing the implications of new value of ε /ε within the SM for NP in the Littlest Higgs model with T-parity and simplified Z and Z models can be found in [5] and [6], respectively.
The present status of ε /ε in the SM has been reviewed recently in [1], where references to rich literature can be found. After the new results for the hadronic matrix elements of QCD penguin and electroweak penguin (V − A) ⊗ (V + A) operators from RBC-UKQCD lattice collaboration [7,8] and the extraction of the corresponding matrix elements of penguin (V − A) ⊗ (V − A) operators from the CP-conserving K → ππ amplitudes one finds [1] ε /ε = (1.9 ± 4.5) × 10 −4 .
This result differs with 2.9 σ significance from the experimental world average from NA48 [9] and KTeV [10,11] collaborations, suggesting evidence for NP in K decays. But even discarding the lattice results and using instead newly derived upper bounds on the matrix elements of the dominant penguin operators from large N approach [12], one finds at most [1] (ε /ε) SM = (8.6 ± 3.2) × 10 −4 , still 2 σ below the experimental data.
The dominant uncertainty in the SM prediction for ε /ε originates from the partial cancellation between QCD penguin contributions and electroweak penguin contributions that depend sensitively on the parameters B , respectively. QCD penguins give a positive contribution while electroweak penguins a negative one. Fortunately a new insight in the values of these parameters has been obtained recently through the results from the RBC-UKQCD collaboration on the relevant hadronic matrix elements of the operators Q 6 [8] and Q 8 [7] and upper bounds on both B (1/2) 6 and B (3/2) 8 which can be derived from large N approach [12].
The results from the RBC-UKQCD collaboration imply the following values for B and the bounds from large N approach read [12] B While one finds in this approach B is less precise but there is a strong indication that B (1/2) 6 < B (3/2) 8 in agreement with (4). For further details, see [12].
In 331 models we have with the ∆(ε /ε) resulting from tree-level Z and Z exchanges. In this paper we will concentrate exclusively on this shift in ε /ε, which as we will see has significantly smaller theoretical uncertainties than the SM part. Indeed, as demonstrated by us in [4], the shift in ε /ε in question is in 331 models governed by the electroweak (V − A) × (V + A) penguin operator with only small contributions of other operators that we will neglect in what follows.
As the relevant non-perturbative parameter B is much better known than B (1/2) 6 , the non-perturbative uncertainty in NP contributions in (6) is significantly smaller than in the SM term. Moreover, except for the value of B (3/2) 8 , the NP contributions are fully independent from the SM one. Consequently we can fully concentrate on the these new contributions and investigate which 331 models can bring theory closer to data. It will also be interesting to see what this implies for other flavour observables, in particular branching ratios for B s,d → µ + µ − , B → K(K * )νν, K + → π + νν, K L → π 0 νν and the Wilson coefficient C 9 that enters the discussion of B → K * µ + µ − anomalies.
In this context our detailed analyses of FCNC processes in 331 models in [14,15] will turn out to be useful. References to earlier analyses of flavour physics in 331 models can be found there and in [16,17]. But the ratio ε /ε has been analyzed in 331 models only in [4]. However, in that paper values of B (1/2) 6 as high as 1.25 and thus violating the bounds in (5) have been considered. As seen in Fig. 14 of that paper in this case the resulting ε /ε in the SM can even be larger than the data so that dependently on the chosen value of B (1/2) 6 both enhancements and suppressions of ε /ε in a given model were required to fit the data in (2) with different implications for K L → π 0 νν. With the new results in (4) and (5) the situation changed drastically and one needs a significant enhancement relative to the SM prediction.
The main new aspects of the present paper relative to [4] are as follows: • We update our analysis of ∆(ε /ε) by taking new result on B into account. We also include in the discussion second fermion representation (F 2 ) which was not considered by us in the context of ε /ε previously.
• After the constraints from ∆F = 2 transitions have been taken into account the size of the possible enhancement of ε /ε depends on a given model and in certain models it is too small to be relevant. Such models are then disfavoured.
• Further selection of the models is provided through the correlation of ε /ε with other flavour observables, in particular the decays B s → µ + µ − , B → K * µ + µ − , K L → π 0 νν, K + → π + νν and B → K(K * )νν. While a definite selection is not possible at present, because the data is not sufficiently precise, it will be possible in the coming years.
In [4] we have considered several 331 models corresponding to four different values of β, three values of tanβ related to Z − Z mixing and two fermion representations F 1 and F 2 . 24 models in total. Among them 7 have been favoured by electroweak precision tests and we will concentrate our analysis on them. The important result of the present paper is that the requirement of significant enhancement of ε /ε reduces the number of favourite models to 3.
Our paper is organized as follows. In Section 2 we present the general formula for ε /ε in 331 models which is now valid for both fermion representations considered in [4]. In Section 3 we first briefly introduce the 7 favourite 331 models that have been selected in [4] on the basis of electroweak precision data. Subsequently, imposing the constraints from ∆F = 2 transitions in K and B s,d systems, we find that only three models are of interest for ε /ε. Subsequently we demonstrate that the correlations of ∆(ε /ε) with rare decays in these three models can provide further selection between them when the data on flavour observables considered by us improves. In Section 4 we consider the case of a Z outside the reach of the LHC. The particular flavour structure of 331 models implies that for M Z ≥ 10 TeV, NP effects in rare B s,d decays, K + → π + νν and K L → π 0 νν are very small, while the ones in K and ε /ε can be significant even for M Z = 30 TeV. We conclude in Section 5. Except for the formulae for ε /ε, all other expressions for observables considered by us can be found in [4,14,15] and we will not repeat them here. But in Table 3 we give all relevant input parameters which occassionaly differ from the ones used in [4,14,15].

Preliminaries
In 331 models ε /ε receives the dominant new contribution from tree-level Z exchanges but through Z − Z mixing, analyzed in detail in [4], contributions from tree-level Z exchanges in certain models and for certain values of new parameters cannot be neglected. We begin with general expressions valid for both Z and Z contributions which will allow us to recall all the relevant parameters of the 331 models. Subsequently we will specify these expressions to Z and Z cases which differ only by the value of the Wilson coefficient of the Q 8 operator at the low renormalization scale at which the relevant hadronic matrix element of Q 8 is calculated.
The basic expression for V = (Z , Z) contribution to ε /ε is given by The factor ω differs by 10% from the corresponding factor ω + in [4] as discussed in [1].
In evaluating (8) we use, as in the case of the SM, the experimental values for ReA 2 and ε K : The amplitude A 2 (V ) is dominated by the contribution of the Q 8 operator and is given by with C 8 (m c , V ) given below. The hadronic matrix element is given by The normalization of our amplitudes is such that Q 8 (m c ) 2 is by a factor of 3 2 smaller than the one in [1]. Correspondingly the value of ReA 2 in (9) is by this factor smaller than in [1]. The choice of the scale µ = m c is convenient as it is used in analytic formulae for ε /ε in [1]. In our numerical calculations we will use B New sources of flavour and CP violation in 331 models are parametrized by new mixing parameters and phasess 13 ,s 23 , withs 13 ands 23 positive definite and 0 ≤ δ 1,2 ≤ 2π. They can be constrained by flavour observables as demonstrated in detail in [14].
Noticeably, constraints deriving from ∆F = 2 observables do not depend on the choice of the fermion representation. We recall here that, as already mentioned, the choice of the transformation properties of the fermions under the gauge group of 331 models, and in particular under the action of SU (3) L , is not unique. Following [4] we denote by F 1 the fermion representation in which the first two generations of quarks transform as triplets under SU (3) L while the third one as well as leptons transform as antitriplets. On the other hand, F 2 corresponds to the case in which the choice of triplets and antitriplets is reversed. Right-handed fermions are always singlets.

Z Contribution
For the fermion representation F 1 we find with 1.35 being renormalization group factor calculated for M Z = 3 TeV in [4].
We have indicated that the values of g 2 and s 2 W should be evaluated at M Z with s 2 W = sin 2 θ W = 0.249 and g 2 (M Z ) = 0.633 corresponding to M Z = 3 TeV.
The coupling ∆ sd L (Z ) is given in terms of the parameters in (14) as follows The formulae (15)- (17) are valid for the fermion representation F 1 . For a given β, the formulae for the fermion representation F 2 are obtained by reversing the sign in front of β. We note that ∆ sd L (Z ) is independent of the fermion representation as f (β) depends only on β 2 .
Combining all these formulae we find with the upper sign for F 1 and the lower for F 2 We observe that the contribution of Z to ε /ε is invariant under the transformation This invariance is in fact valid for other flavour observables in the absence of Z − Z mixing. But as pointed out in [4] in the presence of this mixing it is broken as we will see soon.

Z Contribution
In the case of tree-level Z contribution also the operator Q 8 dominates but its Wilson coefficient is given first by [20] Here g 2 = g 2 (M Z ) = 0.652 is the SU (2) L gauge coupling and the factor 0.76 is the outcome of the RG evolution.
In the 331 models the flavour violating couplings of Z are generated through Z − Z mixing. They are given by [4] with describing the Z −Z mixing, where s 2 W = 0.23126. It should be stressed that this mixing is independent of the fermion representation. Here As the vacuum expectation values of the Higgs triplets ρ and η are responsible for the masses of up-quarks and down-quarks, respectively, we expressed the parameter a in terms of the usual tanβ where we introduced a bar to distinguish the usual angle β from the parameter β in 331 models. See [4] for further details.
We thus find Combining these formulae allows to derive a simple relation where The upper sign in the expression (26) is valid for fermion representation F 1 , the lower for F 2 . The values of R ε are listed in Table 1 for various values of β, tanβ and the representations F 1 and F 2 . Evidently Z dominates NP contributions to ε /ε implying that Z − Z mixing effects are small in this ratio. The two exceptions are the case of β = −1/ √ 3 and tanβ = 5 and the case of β = 1/ √ 3 and tanβ = 0.2 for which Z contribution reaches 50% of the Z one. However, as seen in Table 2 both models are not among favourits and the largest Z − Z effect of 25% among the chosen models in that table is found in M6.
It should be noted that whereas the Z contribution to ε /ε for the representation F 2 differs from the one for F 1 by sign, the contribution of Z to ε /ε is independent of the fermion representation. This disparity breaks the invariance in (19). Analogous feature is observed in several other flavour observables.

Final Formula
The final expression for the shift of ε /ε in 331 models is given by In the next section we will investigate, which 331 models can provide a significant shift ∆(ε /ε) for M Z = 3 TeV and how this shift is correlated with other flavour observables.

Favourite 331 Models
The favourite 331 models selected in [4] on the basis of their perfomance in electroweak precision tests are listed in the notation of that paper in Table 2. In addition to the fermion representation and the values of β and tanβ in a given model we indicate how in that model NP effects in the branching ratio for B s → µ + µ − are correlated with the ones in C 9 . For B(B s → µ + µ − ) the signs ± denote the enhancement and suppression of it with respect to its SM value, respectively. C 9 in the SM is positive and ± also here denote the enhancements and suppressions with respect to its SM value, respectively. In M6 and M11 NP contributions to C 9 are fully negligible. These correlations are shown in Fig. 15 of [4]. Before entering the discussion of ε /ε let us recall that the present data favour simultaneous suppressions of B(B s → µ + µ − ) and C 9 . From Fig. 15 in [4] and Table 2 we reach the following conclusions.
• Qualitatively models M3, M14 and M16 can provide simultaneous suppression of B(B s → µ + µ − ) and a negative shift ReC NP 9 but the suppression of B(B s → µ + µ − ) is not significant.
• For softening the B d → K * µ + µ − anomaly the most interesting is the model M16.
If the anomaly in question remains but decreases with time also models M3 and M14 would be of interest.
• The remaining four models, in fact the four top models on our list of favourites in (28) below, as far as electroweak precision tests are concerned, do not provide any explanation of B d → K * µ + µ − anomaly but are interesting for B s → µ + µ − decay. These are M6, M8, M9 and M11, the first two with F 1 and the last two with F 2 In the last column we list the values of sin(δ 2 − δ 1 ) for which the maximal positive shifts of ε /ε in a given model can be obtained.
fermion representation. It turns out that the strongest suppression of the rate for B s → µ + µ − can be achieved in M8 and M9. In fact these two models are the two leaders on the list of favourites in (28). But in these models C 9 is enhanced and not suppressed as presently observed in the data. The suppression of the B s → µ + µ − rate is smaller in M6 and M11 but there the shift in C 9 can be neglected.
We conclude that when the data for B(B s → µ + µ − ) and C 9 improve we will be able to reduce the number of favourite models. But if both will be significantly suppressed none of the models considered here will be able to describe the data. In fact model M2 with F 1 , β = −2/ √ 3 and tanβ = 5 could in principle do this work here but it is disfavoured through electroweak precision tests.
Concerning these tests the ranking is given as follows M9, M8, M6, M11, M3, M16, M14, (favoured) with the first five performing better than the SM while the last two basically as the SM. The models with odd index I correspond to tanβ = 1.0 and the ones with even one to tanβ = 5.0. None of the models with tanβ = 0.2 made this list implying reduced impact of Z − Z mixing on ε /ε and small NP effects in decays with neutrinos in the final state.

Predictions for ε /ε in Favourite Models
After the recollection of the correlations among B physics observables in the seven models in questions we are in the position to investigate which of these models allows for significant enhancement of ε /ε. To this end we set the CKM parameters to This choice is in the ballpark of present best values for these three parameters but is also motivated by the fact that NP contributions to ε K in 331 models are rather small for M Z of a few TeV and SM should perform well in this case. Indeed for this choice of CKM parameters we find |ε K | SM = 2.14 × 10 −3 , (∆M K ) SM = 0.467 · 10 −2 ps −1 (30) and |ε K | in the SM only 4% below the data. Due to the presence of long distance effects in ∆M K also this value is compatible with the data. Moreover, the resulting Imλ t = 1.42 × 10 −4 is very close to the central value Imλ t = 1.40 × 10 −4 used in [1]. While our choice of CKM parameters is irrelevant for the shift in ε /ε it matters in the predictions for NP contributions to rare K and B decays due to the their intereference with SM contributions. Next, as in [14], we perform a simplified analysis of ∆M d,s , S ψK S and S ψφ in order to identify oases in the space of four parameters in (14) for which these four observables are consistent with experiment. To this end we use the formulae for ∆F = 2 observables in [4,14] and set all input parameters listed in Table 3 at their central values. But in order to take partially hadronic and experimental uncertainties into account we require the 331 models to reproduce the data for ∆M s,d within ±10%(±5%) and the data on S ψK S and S ψφ within experimental 2σ. As seen in Table 3 the present uncertainties in hadronic parameters relevant for ∆M s,d are larger than 10% but we anticipate progress in the coming years. The accuracy of ±5% should be achieved at the end of this decade.
Specifically, our search is governed by the following allowed ranges: 16.0 (16.9)/ps ≤ ∆M s ≤ 19.5 (18.7)/ps, −0.055 ≤ S ψφ ≤ 0.085, where the values in parentheses correspond to decreased uncertainty. For the central parameters we find in the SM (∆M s ) SM = 18.45/ps, (∆M d ) SM = 0.558/ps, S SM ψφ = 0.037, S SM ψK S = 0.688 . (33) In the case of ε K the status of hadronic parameters is better than for ∆M s,d but the CKM uncertainties are larger and the result depends on whether inclusive or exclusive determinations of |V ub | and |V cb | are used. As we keep these parameters fixed we include this uncertainty by choosing the allowed range for |ε K | below to be roughly the range one would get in the SM by varying |V ub | and |V cb | in their ranges known from tree-level determinations.
The uncertainties in ∆M K are very large both due to the presence of long distance effects and large uncertainty in η cc . We could in principle ignore this constraint but as we will see in the next section it plays a role for M Z above 30 TeV not allowing for large shifts in ε /ε in 331 models for such high values of M Z . In fact as we will explain in the next section it is ∆M K and not ε K which is most constraining the maximal values of ε /ε. But at the LHC scales and even at M Z = 10 TeV the ∆M K constraint is irrelevant. Only for scales above M Z = 30 TeV it starts to play an important role bounding the maximal values of ε /ε. Once the knowledge of long distance effects improves and the error on η cc decreases it will be possible to improve our analysis in this part.
We will then impose the ranges  Table 3: Values of the experimental and theoretical quantities used as input parameters as of June 2015. For future updates see PDG [18], FLAG [19] and HFAG [24].
The search for the oases in question is simplified by the fact that the pair (∆M s , S ψφ ) depends only on (s 23 , δ 2 ), while the pair (∆M d , S ψK S ) only on (s 13 , δ 1 ). The result of this search is similar to the one found in Figs. 5 and 6 in [14] but the oases differ in details because of slight changes in input parameters and the reduced allowed range on S ψφ by about a factor of three.
Having determined the ranges for the parameters (14) we can calculate all the remaining flavour observables of interest. In particular we can eliminate those models listed in Table 2 which are not capable of providing a shift in ε /ε larger than say 4 × 10 −4 . To this end we show in Fig. 1 this shift as a function of ε K for models M8, M9 and M16 and in Fig. 2 this shift for the remaining models. On the basis of these results we observe the following: • Only models M8, M9 and M16 are of interest to us as far as ε /ε is concerned and in what follows we will concentrate our numerical analysis on these three models.
• Interestingly, as mentioned above, the strongest suppression of the rate for B s → µ + µ − can be achieved in M8 and M9 although they have presently difficulties with the LHCb anomalies. Using the formulae in [4] this can be expressed in terms of the relations between the coefficients C NP  [31,32].
• On the other hand M16 is the most interesting model for softening the B d → K * µ + µ − anomaly but cannot help by much in suppressing B s → µ + µ − . One finds in this case which is much closer to one of the favourite solutions in which NP resides dominantly in the coefficient C 9 .
• Thus already on the basis of B physics observables we should be able to distinguish between (M8,M9) and M16. But the common feature of the three models is that they provide a bigger shift in ε /ε when the SM value of ε K is below the data and a positive shift in ε K is required. This is in particular seen in the case of lighter colours describing decreased uncertainties in ∆M s,d .
• Most importantly positive shifts in ε /ε in the ballpark of 6 × 10 −4 are possible in these three models, but they are somewhat reduced when the allowed range for ∆M s,d is reduced. Such shifts could be in principle sufficient to describe the data for ε /ε.
The question then arises whether the models M8 and M9 while enhancing ε /ε can simultaneously suppress the rate for B s → µ + µ − and whether M16 while enhancing ε /ε can simultaneously suppress ReC 9 . Moreover correlation of ∆(ε /ε) with the branching ratios for K + → π + νν and K L → π 0 νν is of great interest in view of NA62 and KOPIO experiments.
We observe: • Models M8 and M9 have similar pattern of deviations from the SM with NP effects only relevant in ε /ε and B s → µ + µ − and in fact positive shift of ε /ε up to 6 × 10 −4 and simultaneous suppression of the rate for B s → µ + µ − up to (15 − 20)% are possible in both models. NP effects in C 9 , K + → π + νν and K L → π 0 νν are very small.
• In M16 NP effects are only relevant in ε /ε, C 9 and K L → π 0 νν. A positive shift of ε /ε up to 5 × 10 −4 and simultaneous negative shift of C 9 up to −0.55 is possible in this model. But then the rate for K L → π 0 νν is predicted to be suppressed by 15% and the one for K + → π + νν by roughly 5%.
We conclude therefore that the distinction between M8 and M9 will be very difficult on the basis of observables considered by us but the distinction between these models and M16 should be possible with improved measurements of B s → µ + µ − and improved determination of NP contribution to C 9 which is expected in the flavour precision era. But if NA62 collaboration finds the rate for K + → π + νν significantly above its SM value all these models will fail in describing the data. The same applies to the models in Fig. 2.

Z Outside the Reach of the LHC
We will next investigate how the pattern of NP effects changes if M Z is above 5 TeV and out of the reach of the LHC. As already pointed out in [14], with increased M Z NP effects in ε K and ∆M K increase relative to the ones in ∆M s,d in 331 models due to particular structure of flavour violating couplings in these models. While the coupling ∆ sd L (Z ) in (17) is proportional to the products 13s23 , the corresponding couplings relevant for B s,d system are given by ∆ bs L (Z ) = each involving only one small parameters ij . ∆M d and the CP asymmetry S ψK S specify the allowed ranges for s 13 and δ 1 , while ∆M s and the CP asymmetry S ψφ specify the allowed ranges for s 23 and δ 2 . In this manner for a given M Z the allowed ranges of the four parameters entering the K meson system are determined. But, can be further constrained by ε K and in particular by ∆M K for sufficiently large M Z as explained below. We refer to the plots in [14].
In order to proceed we would like to point out that with increasing M Z the RG analysis leading to (18) has to be improved modifying this formula to with the upper sign for F 1 and the lower for F 2 . The parameter r ε takes into account additional RG evolution above µ = M Z = 3 TeV into account and reaches r ε = 1.45 for M Z = 100 TeV. This could turn out to be useful in models in which the ∆F = 2 constraints could be eliminated, for instance in the presence of other operators. But in 331 models this is not possible and as we will see for M Z ≥ 50 TeV NP effects in ε /ε are suppressed in all 331 models. In Table 4    In order to analyze NP effects beyond the LHC scales we recall the formulae for the shifts due to NP in ∆F = 2 observables.
In the K meson system we have where r ε describes RG effects above M Z = 3 TeV. These effects are much smaller than in the case of ε /ε and in fact suppress slightly NP contribution to ε K , ∆M K and also ∆M s,d . But even for M Z = 100 TeV this factor amounts to r ε ≈ 0.95 in the NP contributions to these observables and for M Z ≤ 50 TeV this effect can be fully neglected. But we keep this factor in formulae below for the future in case various uncertainties decrease. From (39) and (40) we find with − for F 1 and + for F 2 . From (39) and (41) we have on the other hand 1.60 × 10 4 ε ε Z (43) with + for F 1 and − for F 2 . The SM contribution to ∆M K is subject to much larger hadronic uncertainties than is the case of ε K and it is harder to find out what is the room left for NP contributions. In fact we do not even know the required sign of this contribution. But as we will see soon this constraint could become relevant with improved theory for high values of M Z and we will impose the constraint in (34). For the present discussion we neglect the effects of Z − Z mixing which will be included in the numerics.
The formula (42) represents a correlation between ∆M s,d , S ψK S and S ψφ which for a given β determine all parameters on its r.h.s and the ratio of ε /ε and (∆ε K ) Z . Even if M Z does not enter explicitly this expression, the allowed values fors 13 ands 23 depend on it. Similar comments apply to (43).
The relation between the Z effects in ∆F = 2 master functions S i are related as follows [14] ∆S Indeed NP contributions to ∆M s,d are proportional to ∆S s,d and in ε K to ∆S K . This relation follows then from the fact that Presently the strongest constraints on the parameterss ij come from ∆F = 2 processes. With increasing value of M Z the maximal values ofs ij allowed by these constraints increase with increasing M Z . This in turn has impact on the M Z dependence of maximal values of NP contributions to ∆F = 1 observables that in addition to explicit M Z dependence through Z propagator depend sensitively ons ij . Now if the B 0 s,d −B 0 s,d mixing constraints dominate, which turns out still the case for M Z ≤ 30 TeV, one has where we neglect RG effects. In this case with increasing M Z : • Maximal NP effects in ε /ε increase slowly with RG effects represented by r ε .
• Maximal NP effects in B s,d decays decrease like 1/M Z when their intereference with SM contribution dominates the modifications in the branching ratios.
• Maximal NP effects in K + → π + νν and K L → π 0 νν are independent of M Z .
• Maximal NP effects in ε K and ∆M K increase quadratically with M Z up to the point at whichs 13 ands 23 reach maximal values allowed by the unitarity of the new mixing matrix. But this point is never reached as in particular ∆M K constraint becomes important with increasing M Z much earlier. As maximal NP contributions to ε K and ∆M K allowed by B 0 s,d −B 0 s,d mixing constraints increase fast with increasing M Z these two observables will dominate the allowed ranges fors ij at sufficiently high value of M Z and the pattern of M Z dependences changes. Assuming for simplicity that the maximal values ofs 13 ands 23 have the same M Z dependence we have this time at fixed δ 2 − δ 1 In this case with increasing M Z : • Maximal NP effects in ε /ε decrease up to RG effects represented by r ε as 1/M Z .
• Maximal NP effects in B s,d decays decrease like 1/M 1.5 Z when the intereference with SM contribution dominates the modifications in the branching ratios.
• Maximal NP effects in ∆M s,d decrease as 1/M Z .
A closer inspection of formulae (39), (42) and (43) shows that it is the ∆M K constraint that is most important. Indeed, in order to have NP in ε /ε to be significant, we need sin(δ 2 − δ 1 ) ≈ ±1 with the sign dependent on the model considered as listed in the last column in Table 2. This is allowed by B 0 s,d −B 0 s,d mixing constraints. But then as seen in (42) the shift in ε K can be kept small in the presence of a significant shift in ε /ε by having cos(δ 2 − δ 1 ) very small. However, in the case of ∆M K this is not possible as in this limit (43) reduces to with − for F 1 and and + for F 2 . From the signs of β and sin(δ 2 − δ 1 ) in Table 2 we find therefore that for large values of M Z significant enhancement of ε /ε in 331 models implies uniquely suppression of ∆M K for all 331 models considered and for sufficiently large values of M Z this suppression will be too large to agree with data and this in turn will imply suppression of ε /ε. This suppression of ∆M K is not accidental and is valid in any Z model in which flavour violating couplings of Z to quarks are dominantly imaginary as one can easily derive from (41). For sin(δ 2 − δ 1 ) ≈ ±1 as required in 331 models to get large shift in ε /ε the relevant couplings must indeed be dominantly imaginary but in general Z models this could not necessarily be the case and also enhancements of ∆M K could be possible. When ∆M s,d , ε K and ∆M K constraints are equally important the pattern is more involved but these dependences indicate what we should roughly expect. The main message from this analysis is that with increasing M Z the importance of NP effects in K meson system is likely to increase relative to the one in B s,d systems. But one should be cautioned that this depends also on other parameters and on the size of departures of SM predictions for various observables from the data.
Therefore, a detailed quantitative analysis of this pattern will only be possible when Figure 8: Correlations of ∆(ε /ε) with ε K and K L → π 0 νν in M16 at M Z = 50 TeV . Colour coding as in Fig 1. the room left for NP in the quantities in question will be better known. But, the message is clear: possible tensions in ε /ε and ε K can be removed in 331 models for values of M Z beyond the LHC easier than in rare B s,d decays.
As an example we show in Fig. 6 ∆(ε /ε) versus ε K for the favourite models M8, M9 and M16 at M Z = 10 TeV. We observe in accordance with our arguments that at this mass the maximal effects in ε /ε found at M Z = 3 TeV are still possible and the range for possible values of ε K is significant increased. NP contributions to ∆M K at these scales at M Z = 10 TeV are still at most of ±4% and the ∆M K constraint begins to play a role only for M Z ≥ 30 TeV.
In Fig. 7 we show the correlation of ∆(ε /ε) with B s → µ + µ − for M Z = 10 TeV in models M8 and M9 and the correlations of ∆(ε /ε) with C 9 and K L → π 0 νν in M16. To this end we imposed the constraints in (31), (32) and (34). As expected from M Z dependences discussed above, we observe that when the B 0 s,d −B 0 s,d mixing constraints still dominate, NP effects in B s → µ + µ − and C 9 are significantly decreased while they remain practically unchanged in the case of ε /ε and K L → π 0 νν.
With further increase of M Z NP effects in B s → µ + µ − and C 9 further decrease but the maximal effects in K L → π 0 νν are unchanged. In ε /ε they even increase due to the increase of r ε so that for M Z ≈ 30 TeV the shift in ε /ε can reach approximately 8 × 10 −4 in all three models. But for higher values of M Z the ∆M K constraint becomes important and NP effects in both ε /ε and K L → π 0 νν are suppressed relative to the region M Z ≤ 30 TeV as expected from our discussion of the M Z dependence. We illustrate this in the case of M16 in Fig. 8 for M Z = 50 TeV. The suppression is slightly stronger in the case of K L → π 0 νν because in the case of ε /ε it is compansated by roughly 10% by the increase of r ε . The result for the first correlation in M8 and M9 is similar but NP effects in K L → π 0 νν are tiny in these models.

Summary
In this paper we have updated and improved our analysis of ε /ε in 331 models presented in [4]. The new analysis has been motivated by the new results for this ratio from [1,7,8,12] which show that ε /ε within the SM is significantly below the data.
Considering first seven 331 models selected by us in [4] by electroweak precision tests and requiring a shift ∆(ε /ε) ≥ 4.0×10 −4 by NP, we reduced the number of 331 models to three: M8, M9 and M16 in the terminology of [4]. All three can provide for M Z = 3 TeV a shift in ε /ε of (5 − 6) · 10 −3 and this could in principle be sufficient to bring the theory to agree with data if B (1/2) 6 increases towards its upper bound in the future. Moreover: • Models M8 and M9 can simultaneously suppress B s → µ + µ − but do not offer the explanation of the suppression of the Wilson coefficient C 9 in B → K * µ + µ − (the so-called LHCb anomaly).
• On the contrary M16 offers an explanation of this anomaly simultaneously enhancing ε /ε but does not provide suppression of B s → µ + µ − which could be required when the data improves and the inclusive value of |V cb | will be favoured.
• NP effects in K + → π + νν, K L → π 0 νν and B → K(K * )νν are small which can be regarded as prediction of these models to be confronted in the future with NA62, KOPIO and Belle II results.
• Interestingly for values of M Z well above the LHC scales our favourite 331 models can still successfully face the ε /ε anomaly and also possible tensions in ε K can easier be removed than for M Z in the reach of the LHC. This is clearly seen in Fig. 6 obtained for M Z = 10 TeV and similar behaviour is found for M Z up to 30 TeV with maximal NP contribution to ε /ε increased by RG effects up to 8 × 10 −4 . On the other hand, as seen in Fig. 7, the effects in B s → µ + µ − and C 9 are found above M Z = 10 TeV to be very small. For M Z = 50 TeV also effects in ε /ε become too small to be able to explain the ε /ε anomaly.
The possibility of accessing masses of M Z far beyond the LHC reach in 331 models with the help of ε /ε and ε K is very appealing but one should keep in mind that the future of 331 models will crucially depend on the improved theory for ε /ε and ∆F = 2 observables and improved data on rare B s,d and K decays as we have stressed at various places of this writing, in particular when presenting numerous plots.