Some results on lepton flavour universality violation

Motivated by recent experimental measurements on flavour physics, in tension with Standard Model predictions, we perform an updated analysis of new physics violating lepton flavour universality, by using the effective Lagrangian approach and in the Z′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z^{'}$$\end{document} and S3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_3$$\end{document} leptoquark models. We explicitly analyze the impact of considering complex Wilson coefficients in the analysis of B-anomalies, by performing a global fit of RK\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{K}$$\end{document} and RK∗0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{K^{*0}}$$\end{document} observables, together with ΔMs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta M_s$$\end{document} and ACPmix\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_{CP}^{\mathrm{mix}}$$\end{document}. The inclusion of complex couplings provides a slightly improved global fit, and a marginally improved ΔMs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta M_s$$\end{document} prediction.


Introduction
At present, many interesting measurements on flavour physics are performed at the LHC [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Some relevant flavour transition processes in order to constrain new physics at the LHC are the leptonic, semi-leptonic, baryonic and radiative exclusive decays. Some of these decays allow us to build optimized observables, as ratios of these decays, that are theoretically clean observables and whose measurements are in tension with Standard Model (SM) predictions. One example is the case of observables in b → sll transitions. Recently, the LHCb collaboration observed a deviation from the SM predictions in the neutral-current b → s transition [1,2,[5][6][7]11,13], hinting at lepton flavour universality violation effects. The results for ratios of branching ratios involving different lepton flavours are given by [2,11,13], (1 GeV 2 < q 2 < 6 GeV 2 ) = 0.660 +0.110 −0.070 ± 0.024 (0.045 GeV 2 < q 2 < 1.1 GeV 2 ) = 0.685 +0.113 −0.069 ± 0.047 (1.1 GeV 2 < q 2 < 6 GeV 2 ) (1) where the first uncertainty is statistical and the second one comes from systematic effects. In the SM R K = R K * 0 = 1 with theoretical uncertainties of the order of 1% [15,16], as a consequence of Lepton Flavour Universality. The compatibility of the above results with respect to the SM predictions is of 2.6 σ deviation in the first case and for R K * 0 , in the low q 2 di-lepton invariant mass region is of about 2.3 standard deviations; being in the central−q 2 of 2.4 σ . A discrepancy of about 3 σ is found when the measurements of R K and R K * 0 are combined [17]. Anomalous deviations were also observed in the angular distributions of the decay rate of B → K * μ + μ − , being the most significant discrepancy for the P 5 observable [1,6]. The Belle Collaboration has also reported a discrepancy in angular observables consistent with LHCb results [18]. In addition, ATLAS and CMS collaborations have presented their updated results for the angular parameters of the B meson decay, B 0 → K 0 μ + μ − [19][20][21][22]. A great theoretical effort has been devoted to the understanding of these deviations, see for example [15,17,[23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] and references therein. From the theoretical side, the ratios R K and R K * 0 are very clean observables; essentially free of hadronic uncertainties that cancel in the ratios [15]. The experimental data has been used to constrain new physics (NP) models. One useful way to analyze the effects of NP in these observables and to quantify the possible deviations from the SM predictions is through the effective Hamiltonian approach, allowing us a model-independent analysis of new physics effects. In addition, one can compute this effec-tive Hamiltonian in the context of specific NP models. It has been shown that Z and leptoquark models could explain the R K , R K * 0 deviations.
On the other hand, NP models are also severely constrained by other flavour observables, for example in B smixing. Recently an updated computation for the B s mesons mass difference in the SM has been presented [43][44][45][46][47], which shows a deviation with the experimental result [47,48]: such that ΔM SM s > ΔM exp s at about 2 σ . This fact imposes additional constraints over the NP parameter space. Therefore, a global fit is mandatory when considering all updated flavour observables. A negative contribution to ΔM s is needed to reconcile it with the experimental result, in the context of some NP models (like Z or leptoquarks) it implies complex Wilson coefficients in the effective Hamiltonian of R K , R K * 0 [47] (see also below). To the best of our knowledge, most previous works have used only real Wilson coefficients in global fits of R K and R K * 0 observables together with ΔM s , an exception being Ref. [34]. An effect of introducing complex couplings is the generation of C P asymmetries. The mixing-induced C P asymmetry in the B-sector can be , experimentally it is measured to be [48]: In the SM it is given by A mix [47,49,50], with β s = 0.01852 ± 0.00032 [51] we obtain A mix C P SM = −0.03703 ± 0.00064, which is consistent with the experimental result (3) at the ∼ 0.5 σ level.
Reference [34] performed fits for the B-decay physics observables using complex Wilson coefficients, in the model independent and model dependent approaches. The analysis of Ref. [34] performs fits for the B-decay observables using complex couplings, without including the ΔM s or A mix C P observables, then Ref. [34] proceeds to provide predictions to C P-violation observables. Reference [34] only includes ΔM s and A mix C P in the Z -model fit. Our results agree with the ones of Ref. [34] wherever comparable.
The aim of the present work is to investigate the effects of complex Wilson coefficients in the global analyses of NP in B-meson anomalies. We assume a model independent effective Hamiltonian approach and we study the region of NP parameter space compatible with the experimental data, by considering the dependence of the results on the assumptions of imaginary and/or complex Wilson coefficients. We compare our results with the case of considering only real Wilson coefficients. A brief summary of the NP contributions to the effective Lagrangian relevant for b → s transitions and B s -mixing is presented in Sect. 2, where we also recall the need to consider complex Wilson coefficients in the analysis. In Sect. 3 we discuss the effects of having imaginary or complex Wilson coefficients on R K observables. The impact of these complex Wilson coefficients in the analysis of Bmeson anomalies in two specific models, Z and leptoquarks, is included in Sect. 4. We consider a global fit of R K and R K * 0 observables, together with ΔM s and C P-violation observable A mix C P in this analysis. Finally, conclusions are given in Sect. 5.

Effective Hamiltonians and new physics models
The effective Lagrangian for b → s transitions is conventionally given by [52], The Wilson coefficients have contributions from the SM and NP, In the present work we analyze the NP contributions C ( ) NP i . In most of our analysis we will consider the left-handed Wilson coefficients C NP i , the righ-handed Wilson coefficients C NP i are treated briefly in the model-independent approach of Sect. 3 (see Table 1 below). The NP contributions to B smixing are described by the effective Lagrangian [52]: where C L L bs is a Wilson coefficient. In order to study the allowed NP parameter space we follow the same procedure as given in [47], comparing the experimental measurement of the mass difference with the prediction in the SM and NP. Therefore, the effects can be parametrized as [47], where R loop SM = 1.3397 × 10 −3 [47]. The NP prediction to the C P-asymmetry A mix C P is given by [47,49,50] where β s and R the NP Wilson coefficient C L L bs (7) must be negative (C L L bs < 0). In a generic effective Hamiltonian approach, each Wilson coefficient is independent, and setting C L L bs < 0 has no effect on C NP μ 9 , C NP μ 10 , etc. However, explicit NP models give predictions on the Wilson coefficients which introduce correlations among them. We will concentrate on two specific models that have been proposed to solve the semi-leptonic B s -decay anomalies: Z and leptoquarks.
We start with the Z model that contains a Z boson with mass M Z and whose extra NP operators can involve different chiralities. The part of the effective Lagrangian relevant for b → sμ + μ − transitions and B s -mixing is given by [47], where d i and α denote down-quark and charged-lepton mass eigenstates, and λ Q and λ L are hermitian matrices in flavour space. When matching the above equation with Eqs. (4) and (7) one obtains the expressions for the Wilson coefficients at the tree level [47], and where η L L (M Z ) > 0 encodes the running down to the bottom mass scale.
From (12) it is clear that to obtain a negative C L L bs one needs an imaginary number inside the square (λ Q 23 /(V tb V * ts ) ∈ I), but this is the same factor that appears in C NP μ 9 = −C NP μ 10 in (11). λ L 22 ∈ R, since λ is an hermitic matrix, then it follows that C NP μ 9,10 would be imaginary (C NP μ 9,10 ∈ I). Of course, a purely imaginary coupling (or Wilson coefficient) is just a particular and extreme case of having a generic complex coupling. Once one abandons the restriction of considering real couplings it seems more natural to consider the most generic case of complex couplings. There is, however, a motivation to try also the extreme case of imaginary couplings: an imaginary λ Q 23 /(V tb V * ts ) provides a real C L L bs (12), which in turn provides no additional contributions to the C Pasymmetry A mix C P (9), so imaginary couplings might provide a way of improving the predictions on ΔM s without introducing unwanted C P-asymmetries. Now we focus on leptoquark models. Specifically, we consider the scalar leptoquark S 3 ∼ (3, 3, 1/3). The quantum number in brackets indicate colour, weak and hypercharge representation, respectively. The interaction Lagrangian reads [47] where σ a are the Pauli matrices, ε = iσ 2 , and Q i and L α are the left-handed quark and lepton doublets. In this case, the contribution to the Wilson coefficients C NP μ 9,10 arises at the tree level and is given by [47], For C L L bs the contribution appears at the one loop level and can be written as [47,53]: where α = 1, 2, 3 is a lepton family index. Again, in order to obtain C L L bs < 0 in (15), the couplings must ts ) ∈ I, then the expression in Eq. (14) suggests C NP μ 9,10 ∈ I. Of course, the expression (15) is a sum over all generations, so it is possible to set up a model with y QL 32 y QL * 22 /(V tb V * ts ) ∈ R, and to have a cancellation such that the sum in Eq. (15) is imaginary, but this would be a highly fine-tuned scenario. If the sum in (15) has an imaginary part, it would be most natural if all its addends have some imaginary part.
Here we have shown two examples of new physics models which justify the choice of imaginary (or complex) values for the Wilson coefficients C NP μ 9,10 . In the next section we take  Fig. 1 Values of R K and R K * 0 with a imaginary and b real Wilson coefficients an effective Hamiltonian approach and explore whether an imaginary or complex NP Wilson coefficients can accommodate the experimental R K deviations.

Imaginary Wilson coefficients and R K observables
Several groups have analyzed the predictions for the ratios (1) based on different global fits [17,29,30,[34][35][36][37][38][39][40], extracting possible NP contributions or constraining it. As it is well known, an excellent fit to the experimental data is obtained when C NP 9 = −C NP 10 ; corresponding to left-handed lepton currents. By considering this relation, we investigate the effects of having imaginary Wilson coefficients on R K observables. For the numerical evaluation we use inputs values as given in [54]. The SM input parameters most relevant for our computation are: (11), note that the product V tb V * ts , which appears in the computation of Wilson coefficients in NP models (11), (12), (14), (15) is approximately a negative real number (V tb V * ts −0.04). Figure 1 shows the values of the ratios R K and R K * 0 , in their respective q 2 ranges, when both Wilson coefficients C NP μ 9 and C NP μ 10 are imaginary (Fig. 1a) and when they are real (Fig. 1b), by assuming that C . This behaviour was already pointed out in Ref. [26], where it is shown that the interference of purely imaginary Wilson with the SM vanishes, and therefore they can not provide negative contributions to R K , R K * 0 (see also below). In contrast, as shown in the right panel, values of R K ( * 0) ∼ 0.7 (as in the experimental measurements) are possible when the Wilson coefficients are real.
We have done a global fit by including the ratios R K and R K * 0 , and the angular observables P 4 and P 5 [6,19,21,22]. 1 Results are shown in Fig. 2. The allowed regions for imaginary values of C NP μ 9 and C NP μ 10 when fitting to measurements of a series of b → sμ + μ − observables are presented in Fig. 2a, by assuming all other Wilson coefficients to be SMlike. The numerical analysis has been done by using the open source code flavio 0.28 [55], which computes the χ 2 function with each (C NP μ 9 , C NP μ 10 ) pair. The χ 2 difference is evaluated with respect to the SM point, Δχ 2 SM = χ 2 SM −χ 2 min . Then, the pull in σ is defined as Δχ 2 SM , in the case of only one Wilson coefficient, and for the two-dimensional case it can be evaluated by using the inverse cumulative distribution function of a χ 2 distribution having two degrees of freedom; for instance, Δχ 2 = 2.29 for 1 σ . The darker red shaded regions in Fig. 2 correspond to the points with Δχ 2 = χ 2 − χ 2 min ≤ 2.29, that is, they are less than 1 σ away from the best fit point, whereas the lighter red shaded regions correspond to Δχ 2 ≤ 6.18 (≡ 2 σ ). The crosses mark the position of the best fit points. In Fig. 2a the χ 2 function has a broad flat region centered around the origin, with two nearly symmetric minima found Reference [26] showed that imaginary Wilson coefficients do not interfere with the SM amplitude, an therefore imaginary C NP μ 9,10 can not decrease the prediction for R K , R K * 0 . This is numerically shown in the above analysis, where imaginary Wilson coefficients C NP μ 9,10 are not able to reduce significantly the prediction for R K , R K * 0 . To further investigate this question we have analytically computed a numerical approximation to R K * 0 as a function of C NP μ 9 , C NP μ 10 in the region 1.1 ≤ q 2 ≤ 6.0 GeV 2 . After integration and some approximations regarding the scalar products of final state momenta, we obtain: − 0.3013 Re C NP e 10 + · · · · · · +0.0212|C NP μ 9 We have checked that this approximation reproduces the flavio-computed value of R K * 0 to better than 4% in a large region of the parameter space. Now, if we assume that NP does not affect the electron channel (C NP e 9 = C NP e 10 = 0), it is clear that to obtain R K * 0 < R SM K * 0 one needs to introduce C can not bring the prediction of R K * 0 closer to the experimental value. In addition, this expression tells us that the better option to reduce the prediction of R K * 0 is using a real negative C NP μ 9 , and a real positive C NP μ 10 . This is actually the result that we have obtained in our numerical analysis. Figure 1b shows that, for real Wilson coefficients, the lowest prediction for R K * 0 is obtained for C NP μ 9 = −C NP μ 10 < 0, and Fig. 2b shows that the best fit is obtained for negative C NP μ 9 and positive C NP μ 10 . Figure 1a shows that, in general, imaginary Wilson coefficients give positive contributions to R K , R K * 0 , in accordance with Eq. (17). Of course, the full expression is richer than Eq. (17), and we expect some deviations, Fig. 2a shows that the best fit point is not the SM (C NP μ 9 = C NP μ 10 = 0), but the best fit regions are centered around it, and the SM pull with respect the best fit points is small.
We conclude that, actually, a NP explanation for R K , R K * 0 requires that C  , C NP μ 10 should be general complex numbers. Following this reasoning we have performed a global fit to the semi-leptonic decay observables R K , R K * 0 , P 4 and P 5 using generic complex Wilson coefficients allowing only one free Wilson coefficient at a time. Table 1 shows the best fit values, pulls (defined as Δχ 2 SM ) and χ 2 min /d.o.f., for scenarios with NP in one individual complex Wilson coefficient, and Table 2 shows the prediction for R K , R K * 0 for the corresponding central values of each fit, together with the 1 σ uncertainties. The primed Wilson coefficients are also included. We found that the best fit of R K and R K * 0 and the angular distributions is obtained for C Reference [34] also provides fits for complex generic Wilson coefficients. Their scenario I corresponds to our first line in Table 1, our best fit value agrees with their result (C NP μ 9 = (−1.1 ± 0.2) + (0 ± 0.9 i)), within the large uncertainties they give for the imaginary part, but we obtain larger pulls (5.6 σ vs. 4.2 σ of Ref. [34]). Their scenario II corresponds to our third line in Table 1  Choosing complex Wilson coefficients also implies additional constraints from C P-violating observables. This fact has not been considered in the previous analysis. In the next section we study the consequences of having these coefficients in the analysis of B-meson anomalies on some NP models and we consider a global fit of both the ratios R K and R K * 0 and the angular observables P 4 and P 5 , and also the C P-mixing asymmetry.

B s -mixing and NP models
Several NP models that are able to explain the lepton flavour universality violation effects are constrained by other flavour observables like B s -mixing. In particular the parameter space of Z and leptoquark models are severely constrained by the present experimental results of ΔM s [47]. Besides, as already mentioned, additional constraints emerge from C Pviolating observables when considering complex couplings. Reference [47] argues that nearly imaginary Wilson coefficients could explain the discrepancies with the ΔM s experimental measurement, but a global fit of R K and R K * 0 observables, together with ΔM s and C P-violation observable A mix C P in B s → J/ψφ decays should be performed. In the next subsections we investigate these issues for the case of Z and leptoquark models.

Z fit
From now on, a global fit of R K and R K * 0 observables, ΔM s and the C P-violation observable A mix C P is included in our analysis. Figure 3 shows the fits on the Z mass M Z and the imaginary coupling λ Q 23 (setting λ L 22 = 1) imposed by b → sμ + μ − decays and B s -mixing. The red lines (dotted, dash-dotted) correspond to the fit using only semi-leptonic B-meson decays, i.e. b → sμ + μ − as in Fig. 2 plus the branching ratios BR(B s → μ + μ − ) and BR(B 0 → μ + μ − ). The best fit region is the one between the curves; dotted lines: Δχ 2 = 1, dash-dotted lines: Δχ 2 = 4. Blue lines (solid, dashed) correspond to the fit to B s -mixing observables ΔM s and A mix C P . The best fit region is the one between the lines; solid lines Δχ 2 = 1, dashed lines Δχ 2 = 4, there are two regions with Δχ 2 < 1, but between them Δχ 2 is always smaller than 4. The green regions are the combined global fit: dark region Δχ 2 ≤ 1, medium Δχ 2 ≤ 4 and light Δχ 2 ≤ 9.
The best fit for the b → sμ + μ − observables in the region under study is M Z = 11 TeV, λ Q 23 = 0.015 i, with a tiny  (12), with a ΔM s prediction close to the experimental value (2), while the contributions to |C NP μ 9,10 | decrease as M −1 Z (11). Since imaginary couplings worsen the R K , R K * 0 prediction, the larger M Z provides better predictions for them, bringing them closer to the SM value. The best fit Δχ 2 SM grows very slowly with growing allowed M Z . Table 3 summarizes the best fit values for λ Q 23 and M Z , and corresponding pulls, to R K and R K * 0 observables, ΔM s and A mix C P ; considering real, imaginary and complex Wilson coefficients. Results for the above observables in each scenario are included in this table. It is clear that R K and R K * 0 observables prefer real Wilson coefficients, as expected. For real couplings the description is better than the SM, with a pull of 5.39 σ but it does not improve the prediction for ΔM s . Contrary, to improve the prediction for ΔM s imaginary couplings are required in the Z model, however the pull with respect the SM is small, and it has a large χ 2 min /d.o.f.. When allowing generic complex couplings (third column in Table 3) we find that the best fit point is close to the best fit point using only real couplings (first column in Table 3), and the pull with respect the  SM improves slightly (5.43 σ versus 5.39 σ ), and the predictions for the observables are also close to the pure real couplings case, showing a slight improvement in the prediction for ΔM s . Figure 4 shows the best fit regions in the complex λ Q 23 plane for the best fit mass value M Z = 1.08 TeV ( Table 3). The red region shows the 2-dimensional 1 and 2-σ allowed values (Δχ 2 = 2.29, 6.18) including only the b → sμ + μ − observables, the blue region shows the 1 and 2-σ allowed values including only ΔM s , and the green region show the 1 and 2-σ allowed values including only A mix C P , the violet region shows the combined fit. Here we see the tension between the b → sμ + μ − and ΔM s fits. b → sμ + μ − selects a region around the real axis of the coupling, whereas ΔM s selects regions away from it. There are two small intersection regions for the 1-σ allowed values of both fits. The A mix C P fit selects one of these regions, and breaks the degeneracy. Actually, the b → sμ + μ − fit selects fixed values of C  Fig. 4) around the real axis will grow as M 2 Z , but, at the same time, the allowed region will move away from the imaginary axis as M 2 Z . On the other hand, the fit on ΔM s selects fixed values of C L L bs , Eq. (12), since C L L bs ∼ (λ Q 23 ) 2 /M 2 Z , for fixed C L L bs the 1-σ unfavored region around the origin (light blue region in Fig. 4) will grow as λ Q 23 ∼ M Z . As M Z grows, the red region moves away from the origin as M 2 Z , but the blue region expands only as M Z , so that at some M Z value their 1-σ regions do not longer intersect. This is the reason why we obtain a relatively low M Z in the fits of Table 3.
Reference [34] provides also a fit for the Z model, using a fixed M Z = 1 TeV, this value is close to our best fit value of Table 3. For λ L 22 = 1 they obtain the best fit coupling λ Q 23 = (−0.8 ± 0.3) × 10 −3 + (−0.4 ± 3.1) × 10 −3 i with a pull of 4.0 σ . Our best fit values agree with them within uncertainties. Note that we do not provide uncertainties for the best fit values, the reason being that the parameters are not independent, the 2-dimensional best fit regions in Fig. 4 are not ellipses, and the best fit points are not on the center of the figures, so that giving a central value with 1-dimensional uncertainties overestimates the uncertainty and leads to confusion about the meaning and position of the best fit point.
We conclude that, in the framework of Z models, R K − R K * 0 observables are better described than in the SM, with a pull > ∼ 5.39 σ for M Z 1 − 1.3 TeV, and a coupling with a real part Re(λ − 0.0021 improves slightly the fit, as well as the ΔM s prediction.

Leptoquark fit
The leptoquark model has three independent couplings contributing to ΔM s (15). For the global fits we will assume that the dominant coupling is the muon coupling y QL 32 y QL * 22 , which is the one contributing to R K , R K * 0 (14). The fits on the S 3 leptoquark mass M S 3 and the imaginary coupling imposed by b → sμ + μ − decays and B s -mixing are presented in Fig. 5. The observables used in the respective fits are the same as in Fig. 3. The red lines (dotted, dash-dotted) correspond to the fit using only semi-leptonic B-meson decays, i.e. b → sμ + μ − plus the branching ratios BR(B s → μ + μ − ) and BR(B 0 → μ + μ − ), the best fit region is the one between the curves; dotted lines: Δχ 2 = 1, dashdotted lines: Δχ 2 = 4. Blue lines (solid, dashed) correspond to the fit to B s -mixing observables ΔM s and A mix C P . The best fit region is the one between the lines; solid lines Δχ 2 = 1, dashed lines Δχ 2 = 4, there are two regions with Δχ 2 < 1, but between them Δχ 2 is always smaller than 4. The green regions are the combined global fit: dark region Δχ 2 ≤ 1, medium Δχ 2 ≤ 4 and light Δχ 2 ≤ 9. In the b → sμ + μ − fit the best fit parameters for imaginary couplings is y 16. Larger M S 3 masses provide similar values for the best fit couplings, and observable predictions, and the pulls improve slowly. The situation is similar than in the Z case: by allowing larger M S 3 masses the best fit coupling reaches an asymptotic straight line, where the contribution to ΔM s is constant (15), whereas the contribution to |C  Table 4 shows the best fit parameters for the leptoquark model considered in this work, corresponding pulls, predictions to the observables R K , R K * 0 , ΔM s and A mix C P and χ 2 min /d.o.f., considering real, imaginary and complex Wilson coefficients. Table 4 shows that only imaginary couplings do not improve the results, they cannot explain the R K ( * ) anomaly. However, when complex couplings are considered, we found a better global fit of R K , R K * 0 observables, the best global fit parameters emerge at M S 3 = 4.1 TeV and y  Figure 6 shows the best fit regions in the complex y QL 32 y QL * 22 plane, for the best fit mass parameter M S 3 = 4.1 TeV, Table 4. The meaning of each region is as in Fig. 4. In this model there is no intersection between the 1-σ best fit regions of the b → sμ + μ − and the ΔM s fits. Here we also find the tension between the b → sμ + μ − and ΔM s observables, and the different evolution of the best fit regions with the leptoquark mass M S 3 . The ΔM s fit moves the best fit point away from the real axis, and the A mix C P fit selects of the of the signs for the imaginary part, however the global best fit region lies outside the 1-σ region for ΔM s , and the ΔM s prediction does not improve with respect the SM.
Reference [34] also provides a fit for the leptoquark scenario, our model corresponds to their Δ 1/3 [S3] model. Reference [34] performs a fit fixing the leptoquark mass to M S 3 = 1 TeV, and they obtain a two nearly degenerate minimums with positive and negative imaginary parts. The reason for that is that they do not include the A mix C P observable in the fit. Since the C  Table 4, and is inside the best fit region of Fig. 6 (14) are no longer correlated, it would be possible to choose: a purely real coupling to muons, such that it fulfils the first column of Table 4; a vanishing coupling for electrons, such that it does not contribute to R K , R K * 0 ; and a complex coupling for taus, such that y  Table 4. Of course, this would be a quite strange arrangement for leptoquark couplings! Another option would be to take an specific model construction for the relations among the leptoquark couplings, and make a global fit on these parameters. This analysis is beyond the scope of the present work.

Conclusions
In this work, we have updated the analysis of New Physics violating lepton flavour universality, by using the effective Lagrangian approach and also in the Z and leptoquark models. By considering generic complex Wilson coefficients we found that purely imaginary coefficients do not improve significantly B-meson physics observable predictions, whereas complex coefficients (Table 1) do improve the predictions, with a slightly improved pull than using only real coefficients [36]. We have analyzed the impact of considering complex Wilson coefficients in the analysis of B-meson anomalies in two specific models: Z and leptoquarks, and we have presented a global fit of R K and R K * 0 observables, together with ΔM s and C P-violation observable A mix C P when these complex couplings are included in the analysis. We confirm that real Wilson coefficients cannot explain the B s -mixing anomaly; but also only imaginary Wilson coefficients cannot explain the R K , R K * 0 anomaly. Contrary, complex couplings offer a slightly better global fit. For complex couplings the predictions for R K , R K * 0 and ΔM s are similar than for real couplings (Tables 3, 4). For Z models the best fit in both cases is obtained for M Z 1-1.3 TeV, a negative real part of the coupling Re(λ Q 23 ) −0.002, with possibly a similar imaginary coupling part Im(λ Q 23 ) − 0.0021. For leptoquark models the situation is similar, with a best fit mass of M S 3 = 4-5 TeV and a coupling with a positive real part y QL 32 y QL * 22 0.03-0.04, the presence of a similar imaginary part does not improve significantly the fit. One can obtain better fits in the leptoquark models by relaxing the assumption on the leptoquark couplings, or providing specific models for leptoquark couplings, this analysis is beyond the scope of the present work. In summary, new physics Z or leptoquark models with complex couplings provide a slightly improved global fit to B-meson physics observables as compared with models with real couplings.

Note added
After the completion of this work, Ref. [56] appeared also analysing the presence of complex couplings in the Bsystem. Our results agree with Ref. [56] wherever comparable.
(2017SGR754) and CPAN (CSD2007-00042). J.G. thanks the warm hospitality of the Universidad de Zaragoza during the completion of this work.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and no experimental data has been listed.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .