331 models facing new b ->s mu^+ mu^- data

We investigate how the 331 models face new data on B_{s,d}->mu^+mu^- and B_d->K^*mu^+mu^- taking into account present constraints from Delta F=2 observables, low energy precision measurements, LEP-II and the LHC data. In these models new flavour violating interactions of ordinary quarks and leptons are dominated by a new heavy Z^prime gauge boson with the strength of the relevant couplings governed by four new parameters in the quark sector and the parameter beta. We study the implications of these models for beta=+-n/sqrt{3} with n=1,2,3. The case beta=-sqrt{3} leading to Landau singularities for M_{Z'}approx 4Tev can be ruled out when the present constraints on Z^prime couplings, in particular from LEP-II, are taken into account. For n=1,2 interesting results are found for M_{Z^prime}<4Tev with largest NP effects for beta<0 in B_d->K^*mu^+mu^- and the ones in B_{s,d}->mu^+mu^- for beta>0. As Re(C_9^NP) can reach the values -0.8 and -0.4 for n=2 and n=1, respectively the B_d->K^*mu^+mu^- anomalies can be softened with the size depending on Delta M_s/(Delta M_s)_SM and the CP-asymmetry S_{psiphi}. A correlation between Re(C_9^NP) and BR(B_s->mu^+mu^-), identified for beta<0, implies for negative Re(C_9^NP) uniquely suppression of BR(B_s->mu^+mu^-) relative to its SM value which is favoured by the data. In turn also S_{psiphi}K^*mu^+mu^- anomalies in the future data and confirmation of the suppression of BR(B_s->mu^+mu^-) relative to its SM value would favour beta=1/sqrt{3} and M_{Z'}approx 3Tev. Assuming lepton universality, we find an upper bound |C_9^NP|<1.1(1.4) from LEP-II data for all Z' models with only left-handed flavour violating couplings to quarks when NP contributions to Delta M_s at the level of 10%(15%) are allowed.


Introduction
The great expectations to find New Physics (NP) at the LHC did not materialize until now. In particular the order of magnitude enhancements of the branching ratio for B s → µ + µ − decay over its Standard Model (SM) value, possible in supersymmetric models and models with tree-level heavy neutral scalar and pseudoscalar exchanges, are presently ruled out. This is also the case of O(1) values of the CP-asymmetry S ψφ which could also be accommodated in these models. A recent review can be found in [1].
While for the models in question these new flavour data are a big disappointment, for other models like the ones with constrained minimal flavour violation (CMFV), 331 models [2] and Littlest Higgs Model with T-parity [3] they brought a relief as in these models NP effects were naturally predicted to be small. On the other hand the most recent data from LHCb and CMS bring new challenges for the latter models: • The LHCb and CMS collaborations presented new results on B s,d → µ + µ − [4][5][6].
While the branching ratio for B s → µ + µ − as stated above turns out to be rather close to the SM prediction, although a bit lower than the latter, the central value for the one of B d → µ + µ − is by a factor of 3.5 higher than its SM value.
• LHCb collaboration reported new results on angular observables in B d → K * µ + µ − that show departures from SM expectations [7,8]. Moreover, new data on the observable F L , consistent with LHCb value in [7] have been presented by CMS [9].
In particular the anomalies in B d → K * µ + µ − triggered two sophisticated analyses [10,11] with the goal to understand the data and to indicate what type of NP could be responsible for these departures from the SM. Subsequently several other analyses of these data have been presented in [12][13][14][15][16][17] and very recently in [18].
The outcome of these efforts can be summarized briefly as follows. There seems to be a consensus among different groups that NP definitely affects the Wilson coefficient C 9 [10,11,13,17,18] with the value of the shift in C 9 depending on the analysis considered: There is also a consensus that small negative NP contributions to the the Wilson coefficient C 7γ could together with C NP 9 provide the explanation of the data [10,11]. On the other hand as seen in the analyses in [11,13,17], a particularly successful scenario is the one with participation of right-handed currents C NP 9 < 0, C 9 > 0, However a very recent analysis in [18] challenges this solution favouring the one with NP contributions dominantly represented by C NP 9 ≈ −1.5 with much smaller NP contributions to the remaining Wilson coefficients, in particular C 9 .
For the models presented here sorting out these differences is important as in these models C 9 = 0 and as demonstrated in [2,19] NP contributions to C 7γ are totally negligible. Thus C NP 1 Introduction 2 Assuming that indeed NP is at work here, one of the physical mechanisms behind these deviations that seems to emerge from these studies is the presence of tree-level Z exchanges. In [19] we have presented an anatomy of Z contributions to flavour changing neutral current processes (FCNC) identifying various correlations between various observables characteristic for this NP scenario. Recently we have analyzed how this scenario faces the new data listed above [13] including the correlation with the values of C Bq = ∆M q /(∆M q ) SM , S ψφ and S ψK S which should be precisely determined in this decade.
The dominant role in [13] was played by the so-called LHS scenario in which the flavour violating couplings of Z to quarks were purely left-handed. While, in agreement with [11] and recently with [17] it has been found that the presence of right-handed couplings leading to a non-vanishing C 9 gives a better description of the data than the LHS scenario, it is clear that in view of theoretical and experimental uncertainties the LHS scenario remains as a viable alternative.
The nice virtue of the LHS scenario is that for certain choices of the Z couplings the model resembles the structure of CMFV or models with U (2) 3 flavour symmetry. Moreover as no new operators beyond those present in the SM are present, the nonperturbative uncertainties are the same as in the SM, still allowing for non-MFV contributions beyond those present in U (2) 3 models. In particular the stringent CMFV relation between ∆M s,d and B(B s,d → µ + µ − ) [25] valid in the simplest U (2) 3 models is violated in the LHS scenario as analyzed in detail in [13]. Another virtue of the LHS scenario is the paucity of its parameters that enter all flavour observables in a given meson system which should be contrasted with most NP scenarios outside the MFV framework. Indeed, if we concentrate on B 0 s −B 0 s mixing, b → sµ + µ − and b → sνν observables, for a given mass M Z there are only four new parameters to our disposal: the three couplings (our normalizations of couplings are given in Section 2) of which the first one is generally complex and the other two real. The couplings ∆ µμ A,V (Z ) are defined in (19) and due to SU (2) L symmetry implying in LHS ∆ νν L (Z ) = ∆ µμ L (Z ) one also has Extending these considerations to B d and K meson systems brings in four additional parameters, the complex couplings: Thus in this general LHS scenario we deal with eight new parameters. Further reduction of parameters is only possible in a concrete dynamical model. In this context an interesting class of dynamical models representing LHS scenario are the 331 models based on the gauge group SU (3) C × SU (3) L × U (1) X [26,27]. A detailed analysis of FCNC processes in one of these models has been presented by us in [2]. Selection of earlier analyses of various aspects of these models related to our paper can be found in [28][29][30][31][32][33][34][35][36][37].
The nice feature of these models is a small number of free parameters which is lower than present in the general LHS scenario considered in [13,19]. This allows to find certain correlations between different meson systems which is not possible in the general 1 Introduction 3 case. Indeed the strength of the relevant Z couplings to down-quarks is governed by two mixing parameters, two CP-violating phases and the parameter β which defines 1 a given 331 model and determines the charges of new heavy fermions and gauge bosons. Thus for a given M Z and β there are only four free parameters to our disposal. In particular for a given β, the couplings of Z to leptons are fixed. As evident from the general analysis of LHS scenario in [13], knowing the latter couplings simplifies the analysis significantly, increasing simultaneously the predictive power of the theory.
In [2] the relevant couplings have been presented for arbitrary β but the detailed FCNC analysis has been only performed for β = 1/ √ 3. While this model provides interesting results for B s,d → µ + µ − , it fails in the case of anomalies in B d → K * µ + µ − because in this model the coupling ∆ µμ V (Z ) and consequently the Wilson coefficient C NP 9 turn out to be very small. It has been pointed out recently in [14] that for β = − √ 3 a very different picture arises. Indeed in this case ∆ µμ V (Z ) is much larger than for β = 1/ √ 3 so that the B d → K * µ + µ − anomaly can be in principle successfully addressed. Simultaneously the coupling ∆ µμ A (Z ) turns out to be small so that NP contributions to B s → µ + µ − are small in agreement with the data. Moreover aligning the new mixing matrix V L with the CKM matrix, the authors end up with a very simple model in which the only new parameter relevant for their analysis is M Z and the negative sign of C NP 9 required by the B d → K * µ + µ − anomaly is uniquely predicted.
Unfortunately this model has several problems, in particular in the MFV limit considered in [14], which in our view eliminates it as a valid description of the present flavour data. As discussed in Appendix C these problems originate in the known fact that the 331 models with β = ± √ 3 imply a Landau singularity for sin 2 θ W = 0.25 and this value is reached through the renormalization group evolution of the SM couplings for M Z typically around 4 TeV, scales not much higher than the present lower bounds on M Z .
Yet, the observation of the authors of [14] that negative values of β could provide solution to B d → K * µ + µ − anomalies motivates us to generalize our phenomenological analysis of 331 model in [2] from β = 1/ √ 3 to arbitrary values of β, both positive and negative, for which the Landau singularity in question is avoided up to the very high scales, even as high as GUTs scales. This generalization is in fact straight forward as in [2] we have provided formulae for the Z couplings to quarks and leptons for arbitrary β 2 and the expressions for various flavour observables as functions of β can be directly obtained from the formulae of that paper. In this context we will concentrate our analysis on the cases β = ±n/ √ 3 with n = 1, 2 choosing M Z = 3 TeV in order to satisfy existing bounds from flavour conserving observables. A simple scaling law allows then to obtain predictions for other values of M Z .
However, in contrast to our numerical analysis in [2] which assumed certain fixed values of B Bs F Bs and B B d F B d we will investigate in the spirit of our recent paper [13] how our results depend on which should be precisely determined in this decade. 1 Up to the choice of the representation of the gauge group according to which the fermions should transform, as will be better clarified in the next section.
2 The 331 models and their couplings 4 Our paper is organized as follows. In Section 2 we review very briefly the basic aspects of 331 models, recalling their free parameters and the general formulae for the couplings of Z to quarks and leptons for arbitrary β. We also present a table with the values of flavour diagonal couplings of Z to quarks and leptons for n = 1, 2, 3 which should facilitate other researchers to test precisely these models in processes not considered by us. In Section 3 we collect formulae for various Wilson coefficients and one-loop master functions in terms of the couplings of Section 2. This will allow us to identify certain properties and correlations between various observables that will be explicitly seen in our numerical analysis. As all relevant formulae for various branching ratios and other observables have been presented in [2] we recall in Section 4 only crucial observables and their status in the SM and experiment. The strategy for our analysis is presented in Section 5 and its execution in Section 6. In Section 7 we present predictions for low energy precision observables which could provide additional tests of the models considered and analyse also the bounds from LEP-II. We also comment on the bounds on M Z from the LHC. We summarize the main results of our paper in Section 8. Some useful information can also be found in three appendices.

The 331 models and their couplings 2.1 The 331 models
The name 331 encompasses a class of models based on the gauge group SU (3) C × SU (3) L ×U (1) X [26,27], that is at first spontaneously broken to the SM group SU (3) C × SU (2) L × U (1) Y and then undergoes the spontaneous symmetry breaking to SU (3) C × U (1) Q . The extension of the gauge group with respect to SM leads to interesting consequences. The first one is that the requirement of anomaly cancellation together with that of asymptotic freedom of QCD implies that the number of generations must necessarily be equal to the number of colours, hence giving an explanation for the existence of three generations. Furthermore, quark generations should transform differently under the action of SU (3) L . In particular, two quark generations should transform as triplets, one as an antitriplet. Choosing the latter to be the third generation, this different treatment could be at the origin of the large top quark mass. This choice imposes that the leptons should transform as antitriplets. However, one could choose a different scenario in which the role of triplets and antitriplets is exchanged, provided that the number of triplets equals that of antitriplets, in order to fulfil the anomaly cancellation requirement. Therefore, different versions of the model are obtained according to the way one fixes the fermion representations. The fermion representations for specific 331 models analyzed in our paper are described in detail in [2].
A fundamental relation holds among some of the generators of the group: where Q indicates the electric charge, T 3 and T 8 are two of the SU (3) generators and X is the generator of U (1) X . β is a key parameter that defines a specific variant of the model. The 331 models comprise several new particles. There are new gauge bosons Y and V and new heavy fermions, all with electric charges depending on β. Also the Higgs system is extended.
2 The 331 models and their couplings 5 As analyzed in detail in [31] and stated in that paper β can be arbitrary. Yet due to the fact that in 331 models where u is the vacuum expectation value related to the first symmetry breaking it is evident that only values of β satisfying are allowed. With the known value of s 2 W this means that |β| ≤ √ 3 (10) and in fact the only explicit models analyzed in the literature are the ones with β = ±1/ √ 3 and β = ± √ 3. But only for β = ±1/ √ 3 one can avoid the presence of exotic charges both in the fermion and gauge boson sectors. If one considers only then for n = 1 there are singly charged Y ± bosons and neutral ones V 0 (V 0 ), while for n = 3 one finds instead two new singly charged bosons V ± and two doubly charged ones Y ±± . For n = 2 exotic charges ±1/2 and ±3/2 for gauge bosons are found. From Table 1 in [2] we also find that while for n = 1 no exotic charges for heavy fermions are present, for n = 2 heavy quarks carry exotic electric charges ±5/6 and ±7/6 while heavy leptons ±1/2 and ±3/2. Discovering such fermions at the LHC would be a spectacular event. We refer to [2] for further details. In principle β could be a continuous variable satisfying (10) but in the present paper we will only consider the cases n = 1, 2, 3. Most importantly for our paper for all β a new neutral gauge boson Z is present. This represents a very appealing feature, since Z mediates tree level flavour changing neutral currents (FCNC) in the quark sector and could be responsible for the recent anomalies as indicated by their recent extensive analyses.
As in the SM, quark mass eigenstates are defined upon rotation of flavour eigenstates through two unitary matrices U L (for up-type quarks) and V L (for down-type quarks). The relation holds in analogy with the SM case. However, while in the SM V CKM appears only in charged current interactions and the two rotation matrices never appear individually, this is not the case in this model and both U L and V L can generate tree-level FCNCs mediated by Z in the up-quark and down-quark sector, respectively. But these two matrices have to satisfy the relation (12). A useful parametrization for V L which we have used in [2] is 2 The 331 models and their couplings 6 This matrix implies through (12) new sources of flavour violation in the up-sector. However, when U L = 1 as used in [14] V L = V CKM and we deal with a particular simple CMFV model.
With this parametrization, the Z couplings to quarks, for the three meson systems, depend only on four new parameters (explicit formulae are given in [2]): withs 13 ands 23 being positive definite and δ i in the range [0, 2π]. Therefore for fixed M Z and β, the Z contributions to all processes analyzed by us depend only on these parameters implying very strong correlations between NP effects to various observables. Indeed, as seen in (13) the B d system involves only the parameterss 13 and δ 1 while the B s system depends ons 23 and δ 2 . Moreover, stringent correlations between observables in B d,s sectors and in the kaon sector are found since kaon physics depends ons 13 ,s 23 and δ 2 − δ 1 . A very constraining feature of this models is that the diagonal couplings of Z to quarks and leptons are fixed for a given β, except for a weak dependence on M Z due to running of sin 2 θ W provided β differs significantly from ± √ 3.

The couplings
We will now recall those couplings for arbitrary β that are relevant for our paper. The expressions for other couplings and masses of new gauge bosons and fermions as well as expressions for their electric charges that depend on β can be found in [2]. Central for our analysis is the function where the positivity of this function results from the reality of M Z as stressed above.
The following properties should be noted: • For β ≈ √ 3 there is a Landau singularity for s 2 W = 0.25. As at M W one has s 2 W ≈ 0.23 (with exact number depending on its definition considered) and renormalization group evolution of weak couplings increases s 2 W with increasing scale, While we will specifically consider only the cases β = ±n/ √ 3 with n = 1, 2, 3 we list here the formulae for the relevant couplings for arbitrary real β = √ 3 satisfying (10). The case β = √ 3 is considered separately in Appendix A. The important point which we would like to make here is that the couplings of Z to quarks and leptons have to be evaluated at the scale µ at which Z is integrated out, that is at µ = O(M Z ) and not at M W . For n = 1 this difference is irrelevant. For n = 2 it plays a role if acceptable precision is required and it is crucial for n = 3. The values of couplings listed by us and in Appendix A correspond to µ = M Z with the latter specified below. Denoting the elements of the matrix V L in (13) by v ij , the relevant couplings for quarks are then given as follows: The diagonal couplings are valid for the first two generations of quarks neglecting the very small non-diagonal contributions in the matrices V L and U L . For the third generation there is an additional term which can be found in [2]. For leptons we have In Fig. 1 we show ∆ νν L (Z ), ∆ µμ V (Z ) and ∆ µμ A (Z ) as functions of β for s 2 W = 0.249 and g = 0.633 corresponding to M Z = 3 TeV. We observe the following features:  • ∆ µμ V (Z ) can have both signs and for fixed |β| its magnitude is larger for β < 0. In fact the models with negative β are then favoured by the B d → K * µ + µ − anomalies. As noticed in [14] in this case in the limit of CMFV one automatically obtains C NP 9 < 0 as required by experiment. If there are new sources of flavour violation represented by the matrix V L then the oasis in the space of new parameters has to be chosen for which in the case of negative β one still has Re(C NP 9 ) < 0. This will in turn have consequences for other processes as we will see below.

Various 331 models
It is instructive to list the values of the resulting couplings for the models with n = 1, 2, 3 in (11). We do this for flavour diagonal couplings in Table 1 for s 2 W = 0.249 and g = 0.633 valid at M Z = 3 TeV except for n = 3, where we use M Z = 2 TeV in order not to be too close to the Landau singularity. In Appendix A we give explicit formulae for these couplings in terms of sin 2 θ W as well as expressions for flavour violating couplings. Here and in Appendix A we give also Z-couplings.
3 Master formulae for one-loop functions and Wilson coefficients

New physics contributions
We collect here for completeness the corrections to SM one-loop functions and relevant Wilson coefficients as functions of the couplings listed in the previous section.
In the case of ∆F = 2 transitions governed by the function S we have where 3 Master formulae for one-loop functions and Wilson coefficients andr is a QCD factor calculated in [2]. One findsr = 0.965,r = 0.953 andr = 0.925 for M Z = 2, 3, 10 TeV, respectively. It should be remarked that g 2 SM and sin 2 θ W appearing outside the Z couplings, like in (26) and (27) below, should be evaluated at M Z with input values given in Table 3 as they are just related to the overall normalization of Wilson coefficients and rescale relative to SM contributions.
For decays B q → µ + µ − with q = d, s governed by the function Y one has and for Similarly for b → qνν transitions governed by the function X one finds and for K + → π + νν and K L → π 0 νν The corrections from NP to the Wilson coefficients C 9 and C 10 relevant for b → sµ + µ − 3 Master formulae for one-loop functions and Wilson coefficients 10 transitions and used in the recent literature are given as follows 3

Correlations
These formulae imply certain relations that are useful for the subsequent sections. First of all we have the ratio which involves only leptonic couplings and depends only on β. This ratio is given in Table 2 for different values of β and s 2 W = 0.249 except for β = ± √ 3 where we use s 2 W = 0.246. We observe that for β < 0 these two coefficients are predicted to have opposite signs independently of the Z couplings to quarks and as C SM 10 and C SM 9 have also opposite signs (see (53)). Re(C NP 9 ) and NP contributions to B(B s → µ + µ − ) are correlated with each other. This means that Re(C NP 9 ) < 0 required by B d → K * µ + µ − data implies for β < 0 uniquely suppression of B(B s → µ + µ − ) relative to its SM value which is favoured by the data. On the other hand for β = 1/ √ 3 the ratio in (28) is tiny and for β > 1/ √ 3 it is positive implying that NP contributions to B(B s → µ + µ − ) and Re(C NP 9 ) are anti-correlated with each other. Consequently in this case Re(C NP 9 ) < 0 required by B d → K * µ + µ − anomaly implies the enhancement of B(B s → µ + µ − ) which is presently not supported by the data but this could change in the future. We will see all this explicitly in Section 6.
A complementary relation valid in any LHS model that this time does not depend on the lepton couplings is [13] Im(C NP 9 ) Re(C NP 9 ) = Im(C NP 10 ) Re(C NP 10 ) Two important points should be noticed here. These two ratios have to be equal to each other. Moreover they are the same in the two oases resulting from ∆F = 2 constraint. Next the ratios express the relative importance of Z contributions within a given meson system to decays with muons and neutrinos in the final state. While investigating the numbers for R 2 in Table 2 we should recall that in the SM Y ≈ 1.0 while X ≈ 1.5 which means that it is easier to make an impact on decays to muons. While from this ratio we cannot conclude whether a given branching ratio is enhanced or suppressed as the quark couplings cancel in this ratio, the message is clear:  • For β > 0 NP effects in decays to muons governed by axial-vector couplings are much larger than in decays to neutrinos, whereas the opposite is true for β < 0. Therefore in the latter case which is chosen by the B d → K * µ + µ − anomaly we can expect also measurable effects in decays to neutrinos.
• Very importantly in all models, Important are also the relations between the Z contributions to ∆F = 1 (X and Y functions) and ∆F = 2 (S functions) observables. We have where a d = 1 and a s = −1.
We also have ∆X(K) As presently the constraints on 331 models are dominated by ∆F = 2 transitions we observe that for a given allowed size of ∆S(B q ), NP effects in the one-loop functions in question are proportional to 1/M Z and this dependence is transfered to branching ratios in view of the fact that the dominant NP contributions are present there as interference between SM and NP contributions. That these effects are only suppressed like 1/M Z and not like 1/M 2 Z is the consequence of the increase with M Z of the allowed values of the couplings ∆ ij L (Z ) extracted from ∆F = 2 observables, the point already stressed in [2]. In summary, denoting by ∆O NP (M 3 Master formulae for one-loop functions and Wilson coefficients 12 independently of β and C Bq . However the size of NP effect will depend on these two parameters as seen already in the case of β in Table 2 and we will see this more explicitly in Section 6. While this scaling law would apply in the case of the absence of correlations between B q and K systems also to K decays, in the 331 models the situation is different as we will now demonstrate. Indeed in these models there is a correlation between the Z effects in ∆F = 2 master functions in different meson systems where Here and in following equations we set |V tb | = 1 andc 13 =c 23 = 1 if necessary. As the present data and lattice results imply |∆S(B q )| < 0.25 and R 4 < 0.7 in all models, NP effects in ε K are typically below 10%, which is welcome as with input parameters in Table 3 ε K within the SM is in good agreement with the data. Combining then relations (30), (31), (33) and (35) we find with the values of R 5 and R 6 given in Table 2. We observe that ∆X(K) and ∆Y (K) do not depend on M Z when the parameters in V L are constrained through B 0 s,d −B 0 s,d mixings. This fact has already been noticed in [2] but these explicit relations are new. For the models considered in detail by us R 5 < 1.3 and as |∆S(B q )| < 0.25 we find that |∆X(K)| ≤ 0.05 which implies a correction to K + → π + νν and K L → π 0 νν of at most 10% at the level of the branching ratios.
As far as the Wilson coefficients C NP 9 and C NP 10 are concerned we have two important relations which we have written in a form suitable for the analysis in Section 6. We recall that S SM = S 0 (x t ) = 2.31. Both relations are valid for B s and B d systems as indicated on the l.h.s of these equations and the Wilson coefficients on the r.h.s should be appropriately adjusted to the case considered. The virtue of these relation is their independence of the new parameters in (15) so that for a given β the size of C NP 9 and C NP 10 allowed by the ∆F = 2 constraints can be found. In particular in the case of a real C NP 9 and C NP 10 , corresponding to V L = V CKM , ∆S(B q ) and ∆M q will be enhanced which is only allowed if the SM values of ∆M q will turn out to be below the data. If this will not be the case the only solution is to misalign V L and V CKM which results in complex C NP 9 and C NP 10 and consequently novel CP-violating effects.

.1 ∆F = 2 observables
The B 0 s −B 0 s observables are fully described in 331 models by the function where S 0 (x t ) is the real one-loop SM box function and the additional generally complex term has been given in (20). The two observables of interest, ∆M s and S ψφ are then given by and with β s −1 • . In the case of B 0 d system the corresponding formulae are obtained from (41) and (42) by replacing s by d. Moreover (43) is replaced by With the input for |V ub | and |V cb | in Table 3 and γ = 68 • there is a good agreement of the SM with data on S ψK S and ε K . In the SM one has 4 For the central values of the parameters in Table 3 there is a good agreement with the very accurate data [39]: that are not yet included in the FLAG average, the central value of ∆M s would go down to 18.2/ps. Concerning S ψφ and S ψK S we have with the second value known already for some time [39] and the first one being the most recent average from HFAG [39] close to the earlier result from the LHCb [41]. The first value is consistent with the SM expectation of 0.04. This is also the case of S ψK S for the values of |V ub | and |V cb | used by us.
The two Wilson coefficients that receive NP contributions in 331 models are C 9 and C 10 . We decompose them into the SM and NP contributions 5 : where NP contributions have been given in (26) and (27) and the SM contributions are given as follows with all the entries given in [13,19] except for η eff which is new and given below. We have then In the case of B s → µ + µ − decay one has [42][43][44] where with the later value being the latest world average [39]. The bar indicates that ∆Γ s effects have been taken into account. In the SM and CMFV A µµ ∆Γ = 1 but in the 331 models it is slightly smaller and we take this effect into account. Generally as shown in [44] A µµ ∆Γ can serve to test NP models as it can be determined in time-dependent measurements [42,43]. Of interest is also the CP asymmetry which has been studied in detail in [19,44] in the context of Z models. In the case of B d → µ + µ − decay the formulae given above apply with s replaced by d and y d ≈ 0. Explicit formulae for B d → µ + µ − can be found in [19].
Concerning the the status of the branching ratios for B s,d → µ + µ − decays we have The SM values are based on [45] in which NLO corrections of electroweak origin [46] and QCD corrections up to NNLO [47] have been taken into account. These values are rather close to the ones presented previously by us [44,48] but the inclusion of these new higher order corrections that were missing until now reduced significantly various scale uncertainties so that non-parametric uncertainties in both branching ratios are below 2%. The experimental data are the most recent averages of the results from LHCb and CMS [4][5][6].
The calculations performed in [46,47] are very involved and in analogy to the QCD factors, like η B and η 1−3 in ∆F = 2 processes, we find it useful to include all QCD and electroweak corrections into η eff introduced in (52) that without these corrections would be equal to unity. Inspecting the analytic formulae in [45] one finds then The small departure of η eff from unity was already anticipated in [48,49] but only the calculations in [45][46][47] could put these expectations and conjectures on firm footing. Indeed, in order to end up with such a simple result it was crucial to perform such involved calculations as these small corrections are only valid for particular definitions of the top-quark mass and of other electroweak parameters involved. In particular one has to use in Y 0 (x t ) the MS-renormalized top-quark mass m t (m t ) with respect to QCD but on-shell with respect to electroweak interactions. This means m t (m t ) = 163.5 GeV as calculated in [45]. Moreover, in using (61) to calculate observables like branching ratios it is important to have the same normalization of effective Hamiltonian as in the latter paper. There this normalization is expressed in terms of G F and M W only. Needless to say one can also use directly the formulae in [45].
In the present paper we follow the normalization of effective Hamiltonian in [50] which uses G F , α(M Z ) and sin 2 θ W and in order to be consistent with the calculation in [45] our η eff = 0.991 with m t (m t ) unchanged. Interestingly also in the case of K + → π + νν and K L → π 0 νν the analog of η eff , multiplying this time X 0 (x t ), is found with the normalizations of effective Hamiltonian in [50] and definition of m t as given above to be within 1% from unity [51].
In the case of B → K * µ + µ − we will concentrate our discussion on the Wilson coefficient C NP 9 which can be extracted from the angular observables, in particular F L , S 5 and A 8 , in which within the 331 models NP contributions enter exclusively through this coefficient. On the other hand Im(C NP 10 ) governs the CP-asymmetry A 7 . Useful approximate expressions for these angular observables at low q 2 in terms of C NP 9 , C NP 10 and other Wilson coefficients have been provided in [11]. Specified to 331 models they are given as follows S 5 ≈ −0.14 − 0.09 Re C NP 9 .
A 7 ≈ 0.07 Im C NP 10 , 5 Strategy for numerical analysis 16 Note that NP contributions to S 4 and A 9 vanish in 331 models due to the absence of right-handed currents in these models.
Eliminating Re C NP 9 from these expressions in favour of S 5 one finds [13] which shows analytically the point made in [10,11] that NP effects in F L and S 5 are anti-correlated as observed in the data. Indeed, the recent B → K * µ + µ − anomalies imply the following ranges for C NP 9 [10,11] As C SM 9 ≈ 4.1 at µ b = 4.8 GeV, these are very significant suppressions of this coefficient. We note that C 9 remains real as in the SM but the data do not yet preclude a significant imaginary part for this coefficient. The details behind these two results that differ by a factor of two is discussed in [11]. In fact inspecting Figs. 3 and 4 of the latter paper one sees that if the constraints from A FB and B → Kµ + µ − were not taken into account C NP only as suggested originally in [10]. Similarly a very recent comprehensive Bayesian analysis of the authors of [52,53] in [15] finds that although SM works well, if one wants to interpret the data in extensions of the SM then NP scenarios with dominant NP effect in C 9 are favoured although the inclusion of chirality-flipped operators in agreement with [11] would help to reproduce the data. This is also confirmed in [13,17]. However, as we remarked at the beginning of our paper, a very recent analysis in [18] challenges the solution with significant righthanded currents and we are looking forward to the consensus on this point in the future. References to earlier papers on B → K * µ + µ − by all these authors can be found in [10,11,53] and [1].

Strategy for numerical analysis
In our numerical analysis we will follow our recent strategy applied to general LHS models in [13] with the following significant simplification in the case of 331 models. The leptonic couplings of Z are fixed for a given β and this allows us to avoid a rather involved numerical analysis that in [13] had as a goal finding the optimal values of these couplings. Even if β is not fixed and varying it changes the leptonic couplings in question, the ∆ µμ V , ∆ µμ A and ∆ νν L couplings are correlated with each other and finding one day their optimal values in 331 models will also select the optimal value of β fixing the electric charges of new heavy gauge bosons and fermions. Of course also quark couplings will play a prominent role in this analysis, although even if they depend on β, their correlation with leptonic couplings is washed out by the new parameters in (15).
Clearly NP contributions in any extension of the SM are constrained by ∆F = 2 processes which presently are subject to theoretical and experimental uncertainties. However, it is to be expected that in the flavour precision era ahead of us, which will include both advances in experiment and theory, in particular lattice calculations, it will G F = 1.16637(1) × 10 −5 GeV −2 [54]   be possible to decide with high precision whether ∆M s and ∆M d within the SM agree or disagree with the data. For instance already the need for enhancements or suppressions of these observables would be an important information. Similar comments apply to S ψφ and S ψK S as well as to the branching ratios B(B s,d → µ + µ − ) and angular observables in B d → K * µ + µ − . In particular correlations and anti-correlations between suppressions and enhancements allow to distinguish between various NP models as can be illustrated with the DNA charts proposed in [1]. In order to monitor this progress in the context of the 331 models we will consider similarly to [13] the following five bins for C Bs and C B d in (6) C Bs = 0.90 ± 0.01 (yellow), 0.96 ± 0.01 (green), 1.00 ± 0.01 (red), C Bs = 1.04 ± 0.01 (blue), 1.10 ± 0.01 (purple) (70) and similarly for C B d . This strategy avoids variations over non-perturbative parameters like F Bs B Bs and can be executed here because in 331 models these ratios have a very simple form and thanks to the presence of a single operator do not involve any non-perturbative uncertainties. Of course in order to find out the experimental values of these ratios one 6 Numerical analysis 18 has to handle these uncertainties but this is precisely what we want to monitor in the coming years. The most recent update from Utfit collaboration reads C Bs = 1.08 ± 0.09, However, it should be stressed that such values are sensitive to the CKM input and in fact as seen in (45) and (46) with the central values of CKM parameters in Table 3 we would rather expect the central values of C Bs and C B d to be below unity. In order to have full picture we will not use the values in (72) but rather investigate how the results depend on the bins in (70). Concerning S ψφ and S ψK S we will vary them in the ranges − 0.14 ≤ S ψφ ≤ 0.09, 0.639 ≤ S ψK S ≤ 0.719 (73) corresponding to 1σ and 2σ ranges around the central experimental values for S ψφ and S ψK S , respectively. Finally, in order to be sure that the lower bounds from LEP-II and LHC on M Z are satisfied we will present the results for β = ±1/ √ 3 and β = ±2/ √ 3 for M Z = 3 TeV. We will return to this issue in Section 7. The scaling law in (34) allows to translate our results for observables in B s and B d decays into results for other choices of M Z . As we have shown in (37) and (38) the M Z dependence cancels out in ∆X(K) and ∆Y (K).
6 Numerical analysis 6.1 CMFV case It will be instructive to begin our numerical analysis with a particular case, considered in [14], in which In this case the CP-asymmetries S ψφ and S ψK S equal the SM ones and the Wilson coefficients C NP 9 and C NP 10 remain real as in the SM. Moreover, having only two new variables to our disposal, β and M Z , we find very concrete predictions and a number of correlations.
In presenting our results in this section we choose the following colour coding for β: The cases of β = ± √ 3 will be considered separately. In Fig. 2  for different values of β, varying M Z in the range 2 − 5 TeV. We observe: • As already pointed out in [13] and known from CMFV scenario C Bs and C B d are bound to be above unity but this enhancement for the values of β in (75) is not as severe as in the β = − √ 3 case considered in [14]. It should be noted that the sign of β does not matter in these plots and the red and blue lines shown there are equivalent to yellow and green lines, respectively.
6 Numerical analysis

19
• The case of β = ± √ 3 is shown separately in Fig. 4 for fixed s 2 W = 0.246, corresponding to M Z = 2 TeV only as an illustration. Only values of M Z away from singularity are shown and to bring C Bs and C B d down to the acceptable values M Z has to be increased well above M Z = 4 TeV at which Landau singularity is present, that is beyond the range of validity of this model. The authors of [14] working with s 2 W = 0.231 could not see these large enhancements of C Bs and C B d . One can also check that for β = − √ 3 and M Z < 3.5 TeV, in order to stay away from the singularity, |C NP 9 | is much larger than indicated by the data. Clearly, as suggested in [66,67] one could improve this situation by shifting the singularity above 4 TeV through addition of other matter but then the model is a different one and one would have to investigate what impact this additional matter has for observables considered here. To summarize, in this scenario for quark couplings the case β = −2/ √ 3 is performing best as due to large value of ∆ µμ V (Z ) it allows to obtain C NP 9 ≈ −1.0 for C Bs ≈ 1.2. Yet a value C Bs ≈ 1.2 could be problematic when the data improve. For β = −1/ √ 3 which does not introduce exotic charges the required values of C Bs to get sufficiently negative C NP 9 are even larger and the positive values of β are excluded by the required sign of the latter coefficient.
Thus on the whole the idea of the authors of [14] to consider negative values of β was a good one but their choice β = − √ 3 is excluded on the basis of the constraints on C Bs and C B d when the correct values of sin θ 2 W at M Z are used. Moreover LEP-II data on leptonic Z couplings exclude this case as we will see in the next section.
It should finally be noted that NP physics effects in Figs. 2 and 3 appear to be significantly larger than found by us in [2]. One reason for this are different values of β considered here but the primary reason is that the constraints from ∆M s and ∆M d require the matrix V L to be even more hierarchical than V CKM and V L = V CKM in 331 models appears to be problematic as we have just seen. The case V L = V CKM is much more successful as we will demonstrate now.

Preliminaries
Assuming next that the matrix V L in (13) differs from the CKM matrix one has to find first the ranges of parameters (oases) for which a given 331 model agrees with the 6 Numerical analysis  ∆F = 2 data. The outcome of this search, for a given C Bs (or C B d below), are two oases in the space (s 23 , δ 2 ) for B s -system and in the space (s 13 , δ 1 ) for the B d -system. We will not show these oases as they have structure similar to the ones shown in [2,19]. We only recall that the two oases differ by a 180 • shift in the phases δ 1,2 which implies flips of signs of NP contributions to various ∆F = 1 observables. As the ∆F = 2 observables do not change under this shift of phases, the correlation of a given ∆F = 2 observable like S ψφ and ∆F = 1 observable like B(B s → µ + µ − ) in one oasis is changed to the anti-correlation in another oasis and vice versa. Measuring these two observables one can then determine the favoured oasis and subsequently make predictions for other observables.
We will next investigate what happens for the four different values of β considered by us and how the correlations between observables depend on the value of C Bs using the colour coding in (70). We recall that all results for β = ±1/ √ 3 and β = ±2/ √ 3 are obtained for M Z = 3 TeV in order to be sure that the LHC lower bound on M Z is satisfied, although as we will discuss in the next section, for β = ±1/ √ 3 also values M Z ≈ 2.5 TeV are consistent with these bounds. The results for ∆F = 1 observables for M Z = 3 TeV can be obtained by using the scaling law in (34). Then NP effects in β = ±1/ √ 3 could still be by a factor 1.2 larger than shown in the plots below. Concerning the case of β = − √ 3, the problem with too a large C Bs can now be avoided by properly choosing V L but the other problems of this scenario, mentioned at the beginning of our paper and listed in Appendix C, cannot be avoided in this manner and we will not discuss this case any more.

Results
In Fig. 5 we show B(B s → µ + µ − ) versus Re(C NP 9 ) for the four models considered. These four plots exhibit the structure identified through the ratio R 1 in (28) for which numerical values have been given in Table 2. In particular we observe the following features 6 : • For a given C Bs = 1, one can always find an oasis in which Re(C NP 9 ) is negative softening significantly the B d → K * µ + µ − anomalies. However while for β = −2/ • For β < 0 the branching ratio B(B s → µ + µ − ) remains SM-like although in accordance with the relation (28) it is suppressed relative to its SM value for negative Re(C NP 9 ). For the largest values of C Bs (purple and blue lines) this suppression can reach for most negative values of Re(C NP 9 ) 4% and 7% for β = −2/ √ 3 and β = −1/ √ 3, respectively. The slope of the strict correlation between these two observables depends on β. This correlation is presently supported by the data for both observables even if the effects in B(B s → µ + µ − ) are small.
• Looking at these four plots simultaneously we note that going from negative to positive values of β the correlation line moves counter clock-wise with the center of the clock placed at the SM value. This of course means that with increasing beta the correlation B(B s → µ + µ − ) versus Re(C NP 9 ) observed for β < 0 changes into anti-correlation for β > 0, which is rather pronounced in the case of β = 2/ √ 3. Consequently the suppression of B(B s → µ + µ − ) implies positive values of Re(C NP 9 ) which is not what we want to understand B d → K * µ + µ − data. We also note that for β > 0 Re(C NP 9 ) remains small but the effects in B(B s → µ + µ − ) can be larger than for β < 0. These scenarios would be the favorite ones if the B d → K * µ + µ − anomalies decreased or disappeared in the future while the experimental branching ratio B(B s → µ + µ − ) turned out to be indeed by 20% suppressed below its SM value as present central experimental and SM values seem to indicate. In this case the model with β = 1/ √ 3 would be the winner.
This pattern of effects for negative and positive values of β is also seen in Figs. 6 and 7 where we show the correlations between S ψφ versus B(B s → µ + µ − ) and S ψφ versus Re(C NP 9 ), respectively. In particular we find that for models with β < 0 for most negative values of Re(C NP 9 ) and smallest values of B(B s → µ + µ − ) the negative values of S ψφ are favoured. But as the values of S ψφ are rather sensitive for a given value of C Bs to the value of B(B s → µ + µ − ), away from the lower bound on this branching ratio also positive values of S ψφ are allowed. This is in particular the case for largest values of C Bs .
Bearing this ambiguity in mind, we identify therefore for a given C Bs a triple correlation between Re(C NP 9 ), B(B s → µ + µ − ) and S ψφ that is an important test of this model. Interestingly the requirement of a most negative Re(C NP 9 ) shifts automatically the other two observables closer to the data.
While the departure of S ψφ from its SM value is already a clear signal of new sources of CP-violation in ∆F = 2 transitions, non-vanishing imaginary parts of C 9 and C 10 are signals of such new effects in ∆F = 1 transitions. In Figs. 8 and 9 we show the correlations Im(C NP 9 ) versus Re(C NP 9 ) and Im(C NP 10 ) versus Re(C NP 10 ), respectively. The fact that the pattern in both figures for a given β is the same, even if the size of NP effects differs, is related to the relation (29).
We again observe that for β < 0 NP effects are mainly seen in Im(C NP 9 ) while for β > 0 in Im(C NP 10 ). In particular for β = 1/ √ 3 NP effects in C 9 practically vanish which is a good test of this model. Dependently on the values of C Bs and |β|, the CP-asymmetry A 8 could reach (2 − 3)% and the asymmetry A 7 even (3 − 4)% for β < 0 and β > 0, respectively.
Finally in Figs. 10 and 11 we show the correlations Re(C NP 10 ) versus Re(C NP 9 ) and Im(C NP 10 ) versus Im(C NP 9 ) for the four 331 models considered by us. These results follow from (28).
New sources of CP-violation can also be tested in B s → µ + µ − through the asymmetry S µµ defined in (58) and studied in detail in [19,44] in the context of general Z models. In Fig. 12 we show the correlation of S s µµ versus B(B s → µ + µ − ) in 331 model considered. As expected the effects in the models with β > 0 are larger than for β < 0. Similar to the case of S ψφ the required suppression of B(B s → µ + µ − ) favours negative values of S s µµ in all models. As stressed in [11] the Wilson coefficient C NP 9 by itself has difficulty in removing completely the anomalies in B d → K * µ + µ − due to the constraint from B d → Kµ + µ − . We have seen that even without this constraint the values of Re C NP 9 have to be larger than −0.8 but this could turn out to be sufficient to reproduce the data when they improve. Still it is of interest to have a closer look at B d → Kµ + µ − within the four 331 6 Numerical analysis To this end the approximate formula for the branching ratio confined to large q 2 region by the authors of [11] is very useful. Lattice calculations of the relevant form factors are making significant progress here [68,69] and the importance of this decay will increase in the future. Neglecting the interference between NP contributions the formula of [11] reduces in 331 models to where we have used (28) to express Re(C NP 10 ) in terms of Re(C NP 9 ). The error on the first SM term is estimated to be 10% [68,69]. This should be compared with the LHCb result Using (76) we show in Fig. 13 the correlation between B(B d → Kµ + µ − ) [14.2,22] and Re(C NP 9 ) for the four 331 models in question. We observe that the pattern of the correlations is similar to the ones in Fig 5 which originates in the fact that B(B d → Kµ + µ − ) [14.2,22] is strongly correlated with B(B s → µ + µ − ) within LHS models as already shown in [13] for a general LHS model. Moreover, as R 1 is fixed in a given model and its values have been collected in Table 2  • Our results for Re(C NP 9 ) are fully in accordance with the present data on B(B d → Kµ + µ − ) [14.2,22] .
• On the basis of Figs. 5 and 13 there is a triple correlation between B(B d → Kµ + µ − ) [14.2,22] , Re(C NP 9 ) and B(B s → µ + µ − ) which constitutes an important test for the models in question. We indicate this correlation in Fig. 13 by showing when the latter branching ratio is suppressed (black) or enhanced (yellow) with respect to its SM value in accordance with the colour coding in DNA-charts of [1].

Non-CMFV case (B d -system)
We have seen in the case of the MFV limit that B(B d → µ + µ − ) is predicted to be suppressed relative to its SM value when Re(C NP 9 ) is negative. This moves the theory away from the central value of the experimental branching ratio. However, in the non-MFV case we can choose the particular oasis in the space (s 13 , δ 1 ) in which B(B d → µ + µ − ) is enhanced.
In show B(B d → µ + µ − ) versus S ψK S for β = 2/ √ 3. We observe that now enhancement of B(B d → µ + µ − ) can reach 20% over its SM value but still far away from the central experimental value. For β = −1/ √ 3 and β = 1/ √ 3 NP effects turn out to be larger and smaller relative to β = ∓2/ √ 3 respectively, as one could deduce from the values of the axial-vector couplings.
Next in Fig. 15 we show B(B d → µ + µ − ) versus B(B s → µ + µ − ) for the four models considered with the colour coding for β given in (75). We also show the CMFV line. As the uncertainty in the latter line should be reduced to a few percent in this decade, this plot could turn out to be useful for testing and distinguishing the four 331 models.

Preliminaries
Finally, we turn our discussion to decays with neutrinos in the final state. We recall that for given β, C B d , C Bs and the chosen oases in B d and B s systems the corresponding oasis including its size is fixed so that definite predictions for b → sνν transition, K + → π + νν and K L → π 0 νν can be made.
The inspection of the correlations presented in Section 4 teaches us about the follow-6 Numerical analysis 28 Figure 11: Im(C NP 10 ) versus Im(C NP 9 ) for β = ±1/ √ 3 and β = ±2/ √ 3 setting M Z = 3 TeV and different values of C Bs with their colour coding in (70).
ing facts: • NP effects in ε K are small but this is not a problem as with our nominal values of |V ub |, |V cb | and γ SM value of ε K agrees well with the data.
• For β > 0 NP effects in these decays are found to be small but are larger in the cases with β < 0 where Z couplings to neutrinos are largest.
• Similarly NP effects in K + → π + νν and K L → π 0 νν are small as we have already expected on the basis of the relation (37).

The b → sνν transitions
In the absence of right-handed currents one finds [70] and ∆X(B s ) given in (24). The QCD factor η X = 0.994 [71]. In this case the NLO electroweak corrections are of the order of one per mil [51] when similarly to our discussion of η eff in the context of B s,d → µ + µ − decays one uses the normalization of effective Hamiltonian in [50] and the top quark mass is evaluated in the MS scheme for QCD and on-shell with respect to electroweak interactions. Thus accidentally η eff that includes both QCD and electroweak corrections turns out in this scheme to be practically the same for K → πνν and B s,d → µ + µ − decays.
The equality of these three ratios is an important test of any LHS scenario. The violation of them would imply the presence of right-handed couplings at work [70,72,73]. In the context of Z models this is clearly seen in Fig. 20 of [19].
The SU (2) L relation in (4) satisfied in any LHS model, therefore also in the 331 models presented by us, implies a correlation between R νν , B(B s → µ + µ − ) and C NP 9 as shown for a general LHS model in Fig. 9 of [13].
In Fig. 16 we show one of these ratios versus B(B s → µ + µ − ) for the models considered. We observe that in all models considered we have an anti-correlation between these two observables. But the predicted NP effects in all models are rather small. The same conclusion has been reached for general LHS models in [11,13].
effects are rather small. What is interesting are the SM values in the middle of both plots that are enhanced over the usual values quoted as a consequence of inclusive value of |V cb | used by us.

Low and high energy constraints 7.1 Low energy precision observables
Low energy precision observables provide additional bounds on the parameters of the models considered, in particular on the allowed range of M Z as investigated recently in the context of β = − √ 3 model in [14]. We want to add that in concrete models studied here the signs of deviations from SM predictions for these observables are fixed providing additional tests beyond the lower bounds on M Z . In what follows we will present the predictions for three such observables, considered also in [14], separately in each model from which the lower bounds on M Z follow.
We begin with the effect due to a Z gauge boson on the weak charge of a nucleus consisting of Z protons and N neutrons calculated in [74]. In translating this result into our notation one should note that the vector and axial-vector couplings f V,A defined 7 Low and high energy constraints for the four models considered in the paper. The colour coding for β is given in (75). The straight line represents CMFV.
in [74] are not equal to our couplings ∆ V,A (Z ) but are related through 7 Low and high energy constraints  We find then (∆ eē A (Z ) = ∆ µμ A (Z )) 7 Low and high energy constraints  Table 4: Prediction for various observables for different β setting M Z = 3 TeV. Only lower bounds on M Z above 1 TeV resulting from present constraints on these observables are shown.
which has an additional overall factor of −1/4 relative to the corresponding expression in [14] where f V,A = ∆ V,A have been used 7 . We have then Similarly for the effective shift in the weak charge of electron that can be studied in Møller scattering we find For the violation of the first row CKM unitarity expressed through one has for M Z M W [12][13][14]75] In Table 4 we show predictions for these shifts in four models considered by us and in each case the lower bound on M Z that follows from present experimental bounds. In the first case we use, as in [14], Cesium nucleus with Z = 55 and N = 78. We observe that the 90% CL experimental bounds [55,75,76] |∆Q Cs W | ≤ 0.6, |∆Q e W | ≤ 0.016, are well satisfied and the lower bounds on M Z are significantly below the values used by us. We indicated by dashes lower bounds on M Z below 1 TeV. In order to obtain these bounds we neglected running of sin 2 θ W from 3 TeV down to these bounds. Including it would further weaken these bound but this effect is minor.
7 Low and high energy constraints 34

LEP-II constraints
Recently the final analysis of LEP-II data by the LEP electroweak working group appeared in [77] which allows us to check whether the values for M Z for the six 331 models considered by us are consistent with these data. The data relevant for us correspond to the range of center of mass energy 189 GeV ≤ √ s ≤ 207 GeV. In our numerical calculations we will set √ s = 200 GeV. The fundamental for this analysis is the formula (3.8) in this paper [78] In this formula η ij = ±1 or η ij = 0. The different signs of η ij allow to distinguish between constructive (+) and destructive (−) interference between the SM and NP contribution. Λ ± is the scale of the contact interaction which can be related to M Z after proper rescaling of η ij . The lower bounds on Λ ± presented in table 3.15 of [77] apply to certain choices of η ij that are defined in table 3.14 of that paper. In the models considered by us there is the overall minus sign due to Z propagator relative to the SM contribution which we include in the definition of η ij so that with we obtain the relation As we know the signs of η ij in each model we know in each case whether the bound on Λ + or Λ − should be used. In Tables 5-7 we list the values of the couplings η ij for the six models considered by us together with the corresponding values for η ij (Z) for which the minus sign in (88) should be omitted as the energies involved at LEP-II √ s > M Z . The case of β = − √ 3 is easy to test in the case of e + e − → µ + µ − as in this case we deal with the model V V − of [77]. We find then the lower bound for M Z of 11 TeV, well above the validity of this model. We would like to emphasize that this bound is quoted here only as an illustration. As discussed in Appendix C the coupling α X at scales above 1 TeV is too large to trust perturbation theory and calculating only tree diagrams misrepresents the real situation. Whether a non-perturbative dynamics would cure this model remains to be seen.
Before turning to explicit four models analyzed by us let us note that in the LL − , RR − and V V − models for couplings in [77], which correspond to ∆ ll R (Z ) = 0, ∆ ll L (Z ) = 0 and ∆ ll L (Z ) = ∆ ll R (Z ) , respectively, the combination of (89) and (39) allows to derive upper bounds on |C NP 9 | that go beyond the 331 models and apply to LHS scenario for Z generally. Indeed for the case e + e − → µ + µ − we obtain from (89) the bound with a=1 for LL − and RR − and a = 1/2 for V V − . From The last factor becomes unity for a 10% contribution from NP to ∆M s and consequently in this case the maximal by LEP-II allowed values for |C NP 9 | read: 0.91, 0.96 and 1.10, for LL − , RR − and V V − , respectively. The latter case is the one considered in [10] and also has similar structure to β = − √ 3 model without specification of actual values of the muon couplings.
We conclude that for a 10% shift in S it is impossible in these models to obtain C NP 9 = −1.5 as found in [10]. Only for effects S in the ballpark of 20% could such large negative values of C NP 9 be obtained. While these results look similar to the ones shown in Fig. 3, they are more general as they do not assume CMFV and 331 models at work and moreover take into account LEP-II data. Needless to say these LEP-II bounds can be significantly weakened by breaking lepton universality in Z couplings and suppressing Z couplings to electrons relative to the muon ones.
As far as the bound on |C NP 10 | is concerned the bounds obtained for LL − and RR − apply also to this coefficient with ∆ V replaced by ∆ A . For V V − this coefficient vanishes. But for the case AA − in [77] that corresponds to ∆ ll L (Z ) = −∆ ll R (Z ), we find the analogue of the ratio K to be 1.89 TeV and slightly weaker bound than for |C NP 9 | in the V V − case. Thus LEP-II bounds on |C NP 10 | are weaker than the bounds presently available from B s → µ + µ − .
For the remaining models considered by us a complication arises due to the fact that the values of η ij in the simple models studied in [77] and listed in table 3.14 of that paper do not correspond to our models in which generally all combinations of L and R contribute.
However, even without new global fits in these models, which would be beyond the scope of our paper we have checked by using the Tables 5-7, the formulae in Appendix B and the table 3.15 of [77] that the four models considered in detail by us satisfy all LEP-II bounds. In fact our findings are as follows: • For the cases n = −1, 1, 2 the lower bounds on M Z are significantly below 2 TeV, typically close to 1 TeV.
• For β = −2 √ 3 the lower bound on M Z is below 2 TeV but its precise value would require a more sophisticated analysis. In any case it appears that the LHC bound of approximately 3 TeV in this model is stronger that LEP-II bounds.
Concerning the LHC bounds on M Z from ATLAS and in particular CMS [79], the authors of [14] using MAdGraph5 and CTEQ611 parton distribution functions derived for the β = − √ 3 model a 95% CL bound of M Z > 3.9 TeV. As these bounds are based on the Drell-Yan process and are dominated by Z couplings to up-quarks and muons that in 331 models equal to those of electrons, the values of η ij in Table 7 can give us a hint what happens in the models considered by us.
As the relevant η ij in the models considered in detail by us are much lower than the ones in the β = − √ 3 model, the lower bounds on M Z in these models must be significantly lower than 3.9 TeV. On the other hand the couplings in the β = ±2/ √ 3 models are comparable, even if slightly larger than the ones of Z boson. Therefore, we expect that lower bound on M Z could be slightly larger than the one reported by CMS (M Z > 2.9 TeV) and our choice of M Z = 3.0 TeV could be consistent with LHC  Table 5: Values of 10 × η ij for different β relevant for e + e − → + − using sin 2 θ W = 0.249 for β = ±1/ √ 3 and β = ±2/ √ 3 and sin 2 θ W = 0.246 for β = √ 3. sin 2 θ W = 0.231 for Z-couplings.  bounds. Yet in order to find it out a dedicated analysis would be necessary 9 . As far as the LHC bounds for β = ±1/ √ 3 are concerned the analysis in [81] indicates that in these models one could still have M Z ≈ 2.5 TeV. This would allow to enhance NP effects in all ∆F = 1 observables in these models by roughly a factor of 1.2 A complementary lower bound M Z ≥ 1 TeV for 331 models with β = −1/ √ 3 was derived in [82] using dark matter data.
As we have provided all information on the couplings necessary to perform such an analysis in the β = ±2/ √ 3 and β = ±1/ √ 3 models, collider experimentalists and phenomenologists having the relevant codes could derive precise lower bounds on M Z in the models in question. If these bounds turn out in the future to be stronger or weaker than M Z = 3.0 TeV our scaling law in (34) will allow us to translate all results presented in our paper into the new ones.
• The 331 models analyzed by us do not account for the B d → K * µ + µ − anomalies if the latter require Re(C NP 9 ) ≤ −1.3 as indicated by the model independent analysis in [10]. On the other hand, these models could be in accordance with the outcome of the analyses in [11,13,15,17] provided the required size of C 9 in some of these papers will decrease with time (see also [18] where the impact of a NP contribution to C 9 on these anomalies is discussed).
• Going beyond 331 models and assuming lepton universality we find an upper bound |C NP 9 | ≤ 1.1(1.4) from LEP-II data for all Z models within LHS scenario, when NP contributions to ∆M s at the level of 10%(15%) are allowed. We conclude therefore that it is unlikely that values like Re(C NP 9 ) = −1.5 can be accommodated in Z models of LHS type when lepton universality is assumed. As the 331 models not analyzed by us belong to this class of models, this finding applies to them as well.
• The central experimental value of B(B d → µ + µ − ) from LHCb and CMS cannot be reproduced in the 331 models, although an enhancement by 20% over its SM value is possible. A general LHS scenario can do much better as demonstrated in [13]. But then the universality in lepton couplings has to be broken to satisfy LEP-II constraints and the diagonal Z couplings to quarks must be smaller than in 331 models considered by us to avoid the bounds on M Z from LHC.
In more detail our findings are as follows: • Analyzing the models with β = ±1/ √ 3 and β = ±2/ √ 3 we find that for β > 0 measurable NP effects are allowed in B s,d → µ + µ − , sufficient to suppress B(B s → µ + µ − ) down to its central experimental value. On the other hand as mentioned above B(B d → µ + µ − ) even if reaching values 20% above SM result, is still well below the experimental central value. Thus we expect that the experimental value of B(B d → µ + µ − ) must go down if these models should stay alive. For β > 0, the B d → K * µ + µ − anomaly cannot be explained and in fact the anti-correlation between B(B s → µ + µ − ) and Re(C NP 9 ) predicted in this case is not in accordance with the present data. On the other hand in the case of the absence of B d → K * µ + µ − anomalies in the future data and confirmation of the suppression of B(B s → µ + µ − ) relative to its SM value the model with β = 1/ √ 3 and M Z ≈ 3 TeV would be favoured.
• Presently, more interesting appear models with β < 0 where NP effects in B(B s → µ + µ − ) and Re(C NP 9 ) bring the theory closer to the data. Moreover we identified a triple correlation between Re(C NP 9 ), B(B s → µ + µ − ) and S ψφ that for Re(C NP 9 ) < −0.5 required by B d → K * µ + µ − anomalies implies uniquely suppression of B(B s → 8 Summary and conclusions 38 µ + µ − ) relative to its SM value which is favoured by the data. In turn also S ψφ < S SM ψφ is favoured with S ψφ having dominantly opposite sign to S SM ψφ and closer to its central experimental value. Figs. 5-7 show these correlations in explicit terms.
• Another important triple correlation is the one between Re(C NP 9 ), B(B s → µ + µ − ) and B d → Kµ + µ − . It can be found in Fig. 13.
• Our study of b → sνν transitions, K + → π + νν and K L → µ + µ − shows that NP effects in these decays in the models considered are typically below 10% at the level of the branching ratios. NP effects in K L → π 0 νν can reach 20%.
• We have demonstrated how the effects found by us are correlated with the departures of C Bs and C B d from unity. As the latter departures depend sensitively on the precision of lattice non-perturbative calculations, the future of 331 models does not only depend on experimental progress but also on progress of latter calculations.
• As a by-product we have presented bounds on 331 models from low energy precision experiments and provided enough information on the couplings of Z to quarks and leptons that a sophisticated analyses of LEP-II observables and of LHC constraints could be performed in the future.
• Finally, the model with β = − √ 3 can be ruled out on the basis of the data for various observables, in particular the final results from LEP-II. But even if renormalization group effects in sin 2 θ W are not taken into account, the resulting lower bounds on M Z are higher than the upper bounds implied by the Landau singularity. On the other hand the model with β = √ 3 does not predict significant departures from the SM.
Whether the models with β = −1/ √ 3 and β = −2/ √ 3 or with β = 1/ √ 3 and β = 2/ √ 3 will be favoured by the data will depend on the future of the experimental results for B s,d → µ + µ − , B d → K * (K)µ + µ − and future values of C Bq . The numerous plots presented in our paper should allow to monitor these developments. Most importantly, the values of M Z considered in our paper are sufficiently low that this new gauge boson could be discovered in the next run of the LHC and its properties could even be studied at a future ILC [80].

A Expressions for couplings in various 331 models
In obtaining the results below we use for β = ±1/ √ 3 and β = ±2/ √ 3 the values sin 2 θ W = 0.249 and g = 0.633 corresponding to M Z = 3 TeV. We stress that for these models the dependence of the couplings on M Z for 1 TeV ≤ M Z ≤ 5 TeV, unless they are very small, is basically negligible assuring the scaling law (34).
For both signs we have and Next for β = 2/ √ 3 we have and for β = −2/ √ 3 ∆ νν L (Z ) = These results confirm the ones seen in Fig. 1. For completeness we also list the formulae for β = ± √ 3 in order to demonstrate that for β = √ 3 the couplings are too small to provide relevant NP effects, while for β = − √ 3 they are too large to be consistent with the flavour data and LEP-II bounds for M Z < 4 TeV, for which this model is valid because of the Landau singularities in question. In order to stay away from this singularity we give the values of couplings for M Z = 2 TeV, that is for sin 2 θ W = 0.246 and g = 0.636.
For both signs we have In the case of β = √ 3 the formulae for leptonic couplings are modified [2]. We have then and for β = − √ 3

SM couplings of Z
For comparison we give the couplings of Z boson that we evaluate with g = 0.652 and sin 2 θ W = 0.23116 as valid at M Z . The non-diagonal couplings vanish at tree-level and the diagonal ones are given as follows: ∆ µμ A (Z) = g 2c W = 0.372 (102g)

B LEP-II constraints
We will list here formulae which we used to verify that the four 331 models investigated by us satisfy LEP-II constraints on M Z . To this end we generalized the usual SM expressions to include Z contribution. In this context we found the presentation in the book of Burgess and Moore [83] for M Z typically around 4 TeV [28,84] 10 . Therefore these models as they stand, even if V L = V CKM , can only be valid for M Z < 4 TeV. Although in principle some new dynamics entering around these scales could shift the Landau singularity to higher scales, in particular supersymmetry [66,67], one should realize that even at µ = 80 GeV the coupling would be as large as α X ≈ 0.6 that is much larger than all couplings of the SM. At the relevant scales of order few TeV α X ≥ 2.5 implying that perturbative calculations cannot be trusted even in the presence of a large M Z . This is not the problem for other four models discussed by us, where at µ = 3 TeV, the coupling α X equals approximately 0.07 and 0.11 for β = ±1/ √ 3 and β = ±2/ √ 3, respectively. The related problems are as follows • Noting that the masses of the new charged gauge bosons V and Y are related within 1% accuracy to M Z through we find for |β| = √ 3 and s 2 W = 0.24 − 0.25, valid for M Z in the ballpark of a few TeV, the masses of other heavy gauge bosons M V = M Y ≤ M Z /5. This is basically ruled out by the LHC for M Z ≤ 4 TeV. However, a dedicated study would be necessary in order to put this statement on the firm footing. This is not a problem for |β| = 1/ • With the matrix V L equal to the CKM matrix we find that even for values of M Z = (5 − 7) TeV as considered in [14] the mass differences ∆M s and ∆M d are enhanced at least by a factor of two (C B s,d ≈ 2) relative to the SM values. In our view it is unlikely that the future lattice values of B Bs F Bs and B B d F B d would change so much to allow for a satisfactory description of the data for ∆M s,d in this model. While choosing V L = V CKM would remove this problem, this does not help because of the last difficulty.
• It turns out that the size of predicted coupling ∆ µμ V (Z ) in this model implies through LEP-II data a lower bound on M Z of order of 10 TeV when RG effects in sin 2 θ W are taken into account 11 . This value is outside the validity of the model unless complicated new dynamics is introduced at scales of few TeV.