Implications for New Physics in $b\to s \mu\mu$ transitions after recent measurements by Belle and LHCb

We present a Bayesian analysis of the implications for new physics in semileptonic $b\to s$ transitions after including new measurements of $R_K$ at LHCb and new determinations of $R_{K^*}$ and $R_{K^{*+}}$ at Belle. We perform global fits with 2, 4, and 8 input Wilson coefficients, plus one CKM nuisance parameter to take into account uncertainties that are not factorizable. We infer the 68% and 95.4% credibility regions of the marginalized posterior probability density for all scenarios and perform comparisons of models in pairs by calculating the Bayes factor given a common data set. We then proceed to analyzing a few well known BSM models that can provide a high energy framework for the EFT analysis. These include the exchange of a heavy $Z^{'}$ boson in models with heavy vector-like fermions and a scalar field, and a model with scalar leptoquarks. We provide predictions for the BSM couplings and expected mass values.


Introduction
The LHCb Collaboration has recently presented a new measurement of the observable R K , the ratio of the branching fraction of B-meson decay into a kaon and muons, over the decay to electrons, from the combined analyses of the Run 1 and partially of Run 2 data set [1], which reads R K = 0.846 +0.060+0.016 −0.054−0.014 .
This new measurement of R K is compatible with the Standard Model (SM) prediction at 2.5 σ significance. At the same time, the Belle Collaboration has presented new results for the observable R K * in B 0 -meson decays, as well as the first measurement of its counterpart R K * + in B + decays [2]. These new results, listed in Table 1, are consistent with the SM at 1 σ, mainly due to large experimental uncertainties. The rare decays of B mesons are known to provide fertile testing ground for physics beyond the Standard Model (BSM), as in the SM they are highly suppressed by the smallness of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements and by helicity. While for many observables an anomalous determination does not necessarily imply the presence of new physics, as the QCD uncertainties can be sizable, ratios like R K and R K * provide fairly  61 −0.36 ± 0.10 1.40 +1.99 −0.68 ± 0.11 Table 1: R K * and R K * + results from Belle. clean probes, with parametric uncertainties that cancel out to high precision. Additionally, a deviation from the SM in these observables would imply a violation of lepton-flavor universality (LFUV), a purely BSM phenomenon. The rare decays of B-meson involving b → sll transitions have attracted a lot of attention for the search of New Physics (NP) beyond the SM over the last few years. The update of the LHCb results had been long awaited. The Run 1 determinations of R K and R K * [3,4] both featured a 2-3 σ deficit with respect to the SM and were part of a broader set of anomalous measurements in rare semileptonic B decays obtained at LHC and Belle [5,6,7,8,9,10], which involved b → s transitions and muons in the final state, and which in global statistical analyses [11,12,13,14,15,16] had been shown to favor strongly the presence of a few new physics operators over the SM, with a significance that according to some studies reached the ∼ 6 σ level.
In light of the recent measurements of R K and R K * by LHCb and Belle, in this paper we perform a global Bayesian analysis of the BSM effects appearing in the combination of the new data with previous results from the LHC and Belle. In this sense the paper follows the spirit of Run 1 global fits [11,12,13,14,15,16,17,18], and of earlier Bayesian analyses of radiative B-mesons decays [19,20]. The effect of incorporating new data has been also analyzed in Refs. [21,22,23,24,25].
We first consider model-independent fits to the global set of flavor observables, within the framework of the electroweak effective field theory (EFT), assuming as input the Wilson coefficients of four-fermion vector operators that were shown to be able to accommodate the observed data in Run 1 better than the SM. We perform global fits with 2, 4, and 8 independent input parameters, plus one nuisance parameter, the V cb element of the CKM matrix, which takes into account uncertainties that are not factorizable. For each scan we infer the 68.0% and 95.4% credibility regions of the marginalized posterior probability density function (pdf). We then compare models in pairs and provide the Bayes factor for a given set of data. Additionally, we make contact with frequentist analyses by providing for each scan the best-fit point, its pull from the SM, and an estimate of the goodness of fit.
We then proceed to analyzing a few well-known models that can provide a high energy framework for the EFT analysis. Depending on the case and on which set of observables is included in the global fit, BSM interpretations have involved the exchange of a heavy Z boson (see [26,27,28,29,30,31] for early studies) or of a leptoquark (early studies include [32,33,34,35,36,37,38,39]), both with non-universal and flavor-violating couplings to leptons and quarks. Moving one step further, such couplings can be generated assuming heavy vector-like (VL) fermions that can mix with the SM fermions [40,41]. We perform a few global scans for some of the Z models with VL fermions, and for some of the leptoquark models that have been shown to be consistent with the flavor anomalies.
The paper is organized as follows. In Sec. 2, we describe the methodology which is used in our analysis. In Sec. 3, we perform a model-independent global fit to a set of b → s observables. Predictions for several extensions of the SM that can accommodate the observed anomalies are presented in Sec. 4. Finally, we summarize our findings in Sec. 5.

Fit methodology
For each model described by a set of input parameters we map out the regions of the parameter space that are in best agreement with all relevant experimental constraints. To this end we use Bayesian statistics, whose main features we briefly summarize here.
In the Bayesian approach, for a theory described by some parameters m, experimental observables ξ(m) can be compared with data d and a pdf p(m|d), of the model parameters m, can be calculated through Bayes' Theorem. This reads where the likelihood p(d|ξ(m)) ≡ L(m) gives the probability density for obtaining d from a measurement of ξ given a specific value of m, and the prior π(m) parametrizes assumptions about the theory prior to performing the measurement. The evidence, p(d) ≡ Z, is a function of the data that depends globally on the model's parameter space. As long as one considers only one model the evidence is a normalization constant, but it serves as a comparative measure for different models or scenarios.
Bayes' theorem provides an efficient and natural procedure for drawing inferences on a subset of r variables in the parameter space, ψ i=1,..,r ⊂ m. One just needs to marginalize, or integrate, the posterior pdf over the remaining parameters, where n denotes the dimension of the full parameter space. In this work we will be interested in drawing the 68% (1 σ) and 95.4% (2 σ) marginalized 2-dimensional credible regions of the posterior pdf for each model under consideration. We will also compare in pairs different models fitting to the same data set, to determine which one is favored by the data distribution. We do this by computing the Bayes factor, defined as the ratio of evidences for two arbitrary models, M 1 and M 2 : We estimate the significance of Bayes factors according to Jeffrey's scale [42,43]. The central object in our statistical analysis is the likelihood function, constructed using the following prescription. Given the set m of input parameters, which can be, depending on the case, Wilson coefficients, particle masses, coupling constants, or other, the likelihood function is where O th gives a vector of theoretical predictions of the observables of interest and O exp is the vector of the experimental measurements of those observables. We have taken into account the available experimental correlation which is encoded in the matrix C exp . The experimental correlation is available in angular observables for B → K * µµ [6] and B s → φµµ [7]. The theoretical correlation is given by the matrix C th , which is computed using flavio [44]. The theoretical uncertainties, including possible correlations, are estimated by the standard deviation of values of observables, calculated by taking N random values of all input parameters (form factors, bag parameters, decays constants, masses of the particles) distributed according to their probability distribution [44]. In this procedure, the precision depends on the number of random points and we take N = 2000 random points which gives a precision of about 2%. The V cb element of the CKM matrix is treated as a real nuisance parameter. We scan it together with the models' input parameters, following a Gaussian distribution around its central Particle Data Group (PDG) value [45], and adopting PDG uncertainties.

Effective field theory analysis
In the model-independent approach we adopt the weak EFT framework. The effective Hamiltonian for the b → sll transition can be written as: where G F is the Fermi constant and V tb , V ts are elements of the CKM matrix. In Eq. (5) the short-distance physics is encoded in the Wilson coefficients C ( )l i after integrating out the heavy degrees of freedom, whereas the long-distance physics is described by the four-fermion dimension-six interaction operators O ( )l i , invariant under the SU(3) c ×U(1) em gauge group. In this study we will assume the presence of new physics in the following semi-leptonic operators: where the lepton l can be an electron or a muon. We restrict ourselves to the analysis of CP-conserving new physics effects, so that the Wilson coefficient are assumed to be real. We do not consider here new physics in scalar and pseudoscalar operators, O ( ) S and O ( ) P , as they are severely constrained by the B s → µ + µ − measurement [46,47]. Similarly, the electromagnetic dipole operator O ( ) 7 is tightly constrained by radiative decays [48]. The remaining dimension-six operators, chromomagnetic dipole operators, and four-quark operators, at the leading order can only contribute to the semi-leptonic decays through the mixing into semi-leptonic operators. All other BSM contributions enter at a higher order, so that it is safe to consider them as negligible for the purposes of this analysis.
The Wilson coefficients defined in Eq. (5) contain both the SM and new physics (NP) contributions, which can be written as, e.g., New physics contributions to the primed operators can in principle be significant and as such can be considered a smoking gun of BSM phenomena. We perform in this work 6 separate EFT scan fits to the data, each with a different combination of Wilson coefficients as input parameters. We summarize their input ranges and prior distributions in the first 6 lines of Table 2 (here and in what follows we drop the superscript "NP" from the Wilson coefficients' names, but we always take as parameters the new physics contribution). In each determination we also simultaneously scan over the CKM matrix element V cb , which we treat as a nuisance parameter for the Bayesian analysis. We have checked with several preliminary scans that the latter is the only CKM matrix element that can substantially interfere with NP effects due to its large uncertainty, as was pointed out, e.g., in Ref. [11].
This section is dedicated to the discussion of the EFT fits. In Sec. 4, we will instead investigate the implications of the new data for a few popular BSM models, involving a new gauge boson Z or a leptoquark, which are able to provide the favored parameter space regions for the Wilson coefficients. The input parameters, ranges and priors of two of those models, which will be described in the next section, occupy the remaining two lines of Table 2.

Discussion of results
All the observables are calculated with flavio [44], according to the procedure outlined in Sec. 2. In order to efficiently scan the multidimensional parameter space we used MultiNest Parameter Range Prior v.2.7 [49] and pyMultiNest [50] for sampling. The 68% (1 σ) and 95.4% (2 σ) credible regions of the marginalized posterior pdf are computed and plotted with the public tool Superplot [51]. In Fig. 1(a) we show the 1 σ and 2 σ credible regions of the posterior pdf for the model in the first row of Table 2, parametrized by C µ 9 , C µ 10 and the nuisance parameter. The posterior is compared to the one obtained with the previous data, pre LHCb Run 2. The overall value of R K , higher than in the previous determination, has the effect of bringing the 2 σ region closer to the axes origin. The modification of the posterior pdf is not large, but visible. One encounters a less substantial modification of the pdf in the case of the scan parametrized by C µ 9 , C µ 9 (second row of Table 2), for which the credible regions are shown in Fig. 1 The overall effect appears to be a very slight detachment of the 2 σ region from the C µ 9 = 0 axis. A fit to the new data with 4 input NP parameters, C µ 9 , C µ 10 , C µ 9 , C µ 10 (third row of Table 2) shows that the introduction of C µ 10 as a free parameter leads to an interesting interference with C µ 9 . The latter can be made thus comfortably consistent with zero, at the price of introducing a substantial negative value of C µ 10 . This is presented in Fig. 2. In Fig. 2(a) we show a comparison between the marginalized pdf in the (C µ 9 , C µ 10 ) plane for the scan with 2 input NP parameters (first row of Table 2), and the one with 4 NP parameters (third row of Table 2). Larger negative values of C µ 9 are favored by the data with 4 parameters. In Fig. 2(b) we show an equivalent comparison in the (C µ 9 , C µ 9 ) plane, between the marginalized posterior pdf of the 2 parameter scan (second row of Table 2), versus the 4 parameter scan (third row of Table 2). An ample region of the parameter space with C µ 9 ≤ 0 appears, due to the introduction of the parameter C µ 10 . The correlation with C µ 9 is explicitly shown in Fig. 2(c).
Best-fit point 2σ region 1σ region (b) Figure 1: (a) In green, the 1 σ (dark) and 2 σ (light) credible regions of the posterior pdf for the scan in the input parameter C µ 9 , C µ 10 (first row of Table 2), marginalized over the nuisance parameter. The red star marks the position of the best-fit point. The gray solid (dashed) line shows the 1 σ (2 σ) credible region of the pdf corresponding to the data pre-LHCb Run 2. The associated best-fit point is also shown in gray. (b) Same as (a) for the scan parametrized by C µ 9 , C µ 9 (second row of Table 2).
The 2-dimensional regions of the posterior pdf undergo less dramatic modifications if one scans a different set of 4 input parameters: C µ 9 , C µ 10 , C e 9 , C e 10 , see the fourth row of Table 2, or C µ 9 , C µ 9 , C e 9 , C e 9 , see the fifth row of Table 2. We perform the comparison between the relative marginalized 2-dimensional regions of these different models in Fig. 3(a) and Fig. 3(b). The details of the plots are explained in the caption. Note, that the Wilson coefficients of the electron sector, whose pdf's are presented in Fig. 3(c) and Fig. 3(d) remain consistent at 2 σ with zero, implying that the global data set can be easily explained by the presence of NP in the muon sector only.
We finally present in shades of blue the marginalized pdf of the 8-parameter scan introduced in the sixth row of Table 2 in the most relevant planes, (C µ 9 , C µ 10 ) in Fig. 4(a), and (C µ 9 , C µ 9 ) in Fig. 4(b). The posterior regions are compared to the 1 σ and 2 σ regions of the 2-parameters scans in the same plane, presented in shades of green. As one can see, the figures do not differ significantly from Fig. 2(a) and Fig. 2(b), as can be expected given the limited impact the Wilson coefficient of the electron sector bring to the fit.
We summarize the main characteristics of the six fits analyzed in this section in Table 3. The main Bayesian quantity is the negative logarithm of the evidence Z, featured in the second column. We use Jeffrey scale [42,43] to quickly assess the Bayes factor, defined in Sec. 2, which will point to which model is favored by the data. In general, models that are characterized by a smaller number of input parameters tend to fare significantly better than those with a larger input set, as the latter are penalized by volume effects. Unless, of course, these volume effects are counterbalanced by a significantly higher likelihood function.
Best-fit (4 pars.) 2σ region 1σ region Best-fit (4 pars.) 2σ region 1σ region Best-fit (2 pars.) 2σ region 1σ region Best-fit (4 pars.) 2σ region 1σ region (c) Figure 2: (a) In green, the 1 σ (dark) and 2 σ (light) credible regions of the posterior pdf for the scan in the input parameter C µ 9 , C µ 10 (first row of Table 2), compared with the marginalized 2-dimensional regions in the same parameters for the scan with C µ 9 , C µ 10 , C µ 9 , C µ 10 all floating (third row of Table 2), which are shown in brown (1 σ) and orange (2 σ). The red stars mark the position of the best-fit points. (b) A similar comparison of the posterior pdf for the scan in C µ 9 , C µ 9 (shades of green) and the one with C µ 9 , C µ 10 , C µ 9 , C µ 10 all floating (orange/brown). (c) The marginalized 2-dimensional credible regions in C µ 9 , C µ 10 for the scan with C µ 9 , C µ 10 , C µ 9 , C µ 10 all floating (third row of Table 2).
In the specific cases considered here we see that, for example, Best-fit (4 pars.) 2σ region 1σ region Best-fit (4 pars.) 2σ region 1σ region Best-fit (2 pars.) 2σ region 1σ region Best-fit (4 pars.) 2σ region 1σ region (d) Figure 3: (a) In green, the 1 σ (dark) and 2 σ (light) credible regions of the posterior pdf for the scan in the input parameter C µ 9 , C µ 10 (first row of Table 2), compared with the marginalized 2-dimensional regions in the same parameters for the scan with C µ 9 , C µ 10 , C e 9 , C e 10 all floating (fourth row of Table 2), which are shown in brown (1 σ) and orange (2 σ). The red stars mark the position of the best-fit points. (b) A similar comparison of the posterior pdf for the scan in C µ 9 , C µ 9 (second row of Table 2, in shades of green) and the one with C µ 9 , C µ 9 , C e 9 , C e 9 all floating (fifth row of Table 2, in orange/brown). (c) The marginalized 2dimensional credible regions in C e 9 , C e 10 for the scan with C µ 9 , C µ 10 , C e 9 , C e 10 all floating. (d) The marginalized 2-dimensional credible regions in C e 9 , C e 9 for the scan with C µ 9 , C µ 9 , C e 9 , C e 9 all floating.
Best-fit (8 pars.) 2σ region 1σ region Best-fit (2 pars.) 2σ region 1σ region (b) Figure 4: (a) In green, the 1 σ (dark) and 2 σ (light) credible regions of the posterior pdf for the scan in the input parameter C µ 9 , C µ 10 (first row of Table 2), compared with the marginalized 2-dimensional regions in the same parameters for the scan with all 8 NP parameters floating (sixth row of Table 2), which are shown in shades of blue. The stars mark the position of the best-fit points. (b) A similar comparison of the posterior pdf for the scan in C µ 9 , C µ 9 (second row of Table 2, shades of green) and the one with 8 parameters floating (sixth row of Table 2, shades of blue). Z C µ 9 ,C µ 9 Z C µ 9 ,C µ 10 ,C µ 9 ,C µ 10 = 1.5 (Barely worth mentioning) Z C µ 9 ,C µ 10 ,C µ 9 ,C µ 10 Z C µ 9 ,C µ 10 ,C e 9 ,C e 10 = 221 (Very strong) Z C µ 9 ,C µ 10 ,C µ 9 ,C µ 10 Z C µ 9 ,C µ 9 ,C e 9 ,C e 9 = 109 , (Strong) (13) which brings us to the conclusion that the models favored by the data are the ones in the second and third rows of Table 2.
In the third and fourth column of Table 3 we make contact with frequentist approaches by presenting the pull of the best-fit point from the SM, the minimum chi-squared χ 2 TOT , the minimum chi-squared per degree of freedom, and the relative chi-squared of the muon observables, χ 2 µ , electron observables, χ 2 e , and LFUV observables, χ 2 R K and χ 2 R * K , of the six EFT scans analyzed here. We calculate the minimum chi-squared per degree of freedom very roughly, neglecting all correlations, as an indicative measure of the relative goodness of fit: where the ±1 is placed as a reminder of the nuisance parameter. The full list of constraints is collected in Appendix A.     Finally, again to favor comparison with frequentist analyses, we show in Table 4 the numerical value of the Wilson coefficients and of the most important observables at the best-fit points.

Simple models for b → sll anomalies
In the previous section we determined the preferred 1σ and 2σ ranges for the NP Wilson coefficients relevant to explaining b → s anomalies. In the following we will discuss several simple BSM scenarios that are known to naturally lead to the desired EFT operator structure due to an exchange of a BSM boson with flavor-violating couplings to the SM quarks and leptons. Since we confirmed in our model-independent fit (see Table 3) that the presence of the electron Wilson coefficients does not improve the goodness of the fit, we will only consider models in which the BSM boson couples exclusively to muons.

Heavy Z
As a first example, we discuss the SM extended by an additional vector boson, commonly denoted as Z . The most generic Lagrangian, parametrizing non-universal couplings of Z to the b-s current and the muons reads The relevant Wilson coefficients are then given by where ∆ µµ 9 ≡ (∆ µµ R + ∆ µµ L )/2, ∆ µµ 10 ≡ (∆ µµ R − ∆ µµ L )/2, m Z is the mass of the Z boson, and is the typical effective scale of the new physics. If the heavy Z is the gauge boson of a new U(1) X gauge group, its couplings to the gauge eigenstates must be family-universal, and an additional structure is required to generate ∆ sb L and ∆ sb R . Thus, in this work we also consider the impact of the new LHCb and Belle data on the masses and couplings of a few simplified but UV complete models.
In Eqs. (18)- (20) and in the following text we label SM fields by lower-case letters and BSM matter with capital ones. Besides Z , we also add to the SM a scalar singlet field S to spontaneously break the U(1) X symmetry and one family of VL quark pairs Q, Q and D, D to create the flavorchanging couplings ∆ bs L,R [55]: Q : (3, 2, 1/6, −1) Q : (3, 2, −1/6, 1) , With the above choice of the U (1) X charges, the Lagrangian features new Yukawa couplings λ Q,i and λ D,i that mix the SM and BSM fields, as well as VL mass terms M Q,D : where i = 1, 2, 3 label the SM generations. After spontaneously breaking U(1) X by a new scalar vev v S , and rotating the Lagrangian (24) to the mass basis, one obtains the flavor-generating couplings of the form and ∆ µµ 9 = g X ∆ µµ 10 = 0 , where g X is the gauge coupling of the U(1) X group. By recalling that v S ≡ m Z /g X one finally obtains while C µ 10,NP = C µ 10,NP = 0. Without loss of generality one can assume that the couplings of the second and third generations are unified, denoted as λ Q,D . Therefore, the model can be parametrized in terms of only 3 free parameters: m Z /g X , M Q /λ Q , and M D /λ D . Their scanning ranges imposed in the global fit are shown in Table 2. We scan m Z /g X for values larger than 500 GeV to evade the strong bound from neutrino trident production [56].
Recall, finally, that any scenario with a non-universal ∆ bs L,R coupling is subject to the strong constraints from B s mixing [11,57]: R BB ≤ 0.014. In our VL model the latter can be expressed in terms of the Wilson coeffiecints as [26]: where S 0 ≈ 2.3 is a loop factor. For this reason, we impose an upper bound on the prior range of m Z /g X ≤ 5 TeV.

Model 2.
Another realization of the L µ − L τ model we consider is an extension of the SM characterized by one pair of VL quark doublets Q, Q , to generate the flavor-violating coupling of the Z in the quark sector, ∆ bs L , and one pair of VL U(1) X neutral leptons E, E [58,59], which have to be SU(2) singlets for reasons that will be clear below. Thus, Q : (3, 2, 1/6, −1) Q : (3, 2, −1/6, 1) , The Yukawa and VL mass part of the Lagrangian relative to the new leptons reads Assuming that in Eq. (33) only second-generation Yukawa couplings are nonzero one gets Note that if we had chosen a VL lepton SU(2) doublet instead of a singlet we would have obtained C µ 9,NP and C µ 10,NP of the same sign, which is disfavored by the data. Again we parametrize this model in terms of 3 free parameters: m Z /g X , M Q /λ Q , and M E /λ E,2 , where λ Q represent equal couplings to the second and third generation quarks. Their scanning ranges imposed in the global fit are shown in Table 2. We also apply the bounds from B s mixing.
We present in Fig. 5(a) the marginalized 2-dimensional posterior pdf in the (m Z /g X , M Q /λ Q ) plane for Model 2. The VL mass range is determined by the 2 σ range in C µ 9,NP and lies around a 20-30 TeV scale for a coupling λ Q of order unity. Within 2 σ probability, the m Z /g X mass is limited to values below 5 TeV, as a result of the B s mixing constraint, which directly depends on v 2 S . Note that the 2-dimensional posterior pdf in the (m Z /g X , M Q /λ Q ) plane for Model 1 looks very similar to Fig. 5(a), as a direct consequence of the bounds on C µ 9,NP and R BB . In both Model 1 and Model 2, the second VL mass is unbounded from above at the 2 σ level, as depicted in Fig. 5(b), where we show scatter plots of the χ 2 distributions of Model 1 (red points) and Model 2 (blue points) versus the VL mass rescaled by the Yukawa coupling. This is a consequence of the fact that C µ 9,NP in Model 1 and, especially C µ 10,NP in Model 2, are consistent with zero at the 2 σ level. On the other hand, the values of M D /λ D and M E /λ E,2 emerging at the best-fit point are very different for the two cases.
We summarize the main characteristics of the scans of Model 1 and Model 2 in Table 5. In analogy to what we observed in the EFT analysis, the Bayes factor favors Model 1 over Model 2 by 4.9:1, which reads "positive" evidence on Jeffrey's scale.

A model with U(1) X charged VL leptons
We finally consider an alternative to the L µ − L τ model, obtained if one charges the VL leptons under the U(1) X symmetry, and leaves the SM leptons uncharged, see, e.g., [60]. Model 3. We add to the SM the following particle content The lepton Yukawa and VL mass part of the Lagrangian reads After rotating to the quark and lepton mass bases one finds This model can be parametrized in terms of 3 parameters: m Z /g X , M Q /λ Q , and the hierarchical , defined such that M L /λ L,2 = M Q /λ Q .
The 2 σ region of the profile likelihood reads for this 1-dimensional model We apply this bound to Model 3, together with the bound from B s mixing. The favored 2 σ regions are shown in Fig. 5(c), with different color code for different . The severe bound on R BB limits this model to steep hierarchies between VL quark and lepton masses.

Leptoquarks
A second, well known class of models that can easily generate NP contributions to the Wilson coefficients of the EFT are leptoquarks. Much work has been done in the past few years on the phenomenology of leptoquarks in relation to the flavor anomalies, for some early references, see, e.g., [32,33,34,35,36,37,38,39]. Here, we limit ourselves to the analysis of but one of these models, the scalar SU(2) triplet S 3 , which can generate a C µ 9 = −C µ 10 contribution at the tree level, like Model 3 of the previous subsection. We introduce the scalar leptoquark The Lagrangian acquires a Yukawa term where and where the electric charge is indicated in the subscript. After rotating to the mass basis one writes where indices ij span the SM generations, and Y T ij =Ŷ T ik (V † CKM ) kj . Note that couplings of the type QS † 3 Q, which are very dangerous for proton decay, are allowed in the SM, so that in UV complete models one should make sure they are forbidden by an additional symmetry.
By matching to the EFT one finds where in this and what follows we have assumed that the mass of the triplet states is the same, m S 3 . The constraint from the 1-dimensional EFT at 2 σ is given in Eq. (41). This leads to If one starts withŶ bµŶ * sµ = 0, the CKM matrix generates additional Yukawa couplings, so, that possible complementary constraints come from B → K ( * ) νν, b → cµ −ν decays, t → c µ + µ − , and t → c νν. The most dangerous constraint is possibly given by B → K ( * ) νν decay. The bound can be expressed as [61,33] where Br (SM) = (4.0 ± 0.5) × 10 −6 , Br (90% CL) = 1.6 × 10 −5 , and C SM L = −6.38 ± 0.06. We get the limit which does not constrain the parameter space emerging in Eq. (47).

Summary and conclusions
In this paper we have presented a global Bayesian analysis of the new physics effects on effective operators of semileptonic b → s transitions after the very recent measurement of R K at LHCb and new results for the observable R K * in B 0 -meson decays, as well as the first measurement of its counterpart R K * + in B + decays at Belle. We have performed global fits with 2, 4, and 8 Wilson coefficients as inputs, plus one CKM nuisance parameter to take into account uncertainties that are not factorizable with the NP effects. From the fits, we then inferred the 68% and 95.4% credibility regions of the marginalized posterior probability density for all models.
The new measurement of R K is closer in central value to the SM prediction than the Run 1 determination, but the much improved precision of the new data keeps it at 2.5σ from the SM. As a result the high-probability region of the posterior pdf in the NP Wilson coefficients C µ 9 and C µ 10 shifts slightly towards the zero value with respect to the scans with the Run 1 determination of R K , but the overall pull remains quite large, at the level of 4 − 5 σ, quite independently of the number of scanned input coefficients.
We have confirmed previous observations that the impact of the Wilson coefficients of the electron sector on the data is negligible with respect to the muon sector. Moreover, a pair-like comparison of the Bayes factors of different models has allowed us to determine that the two scans characterized by the inputs C µ 9 , C µ 9 , and C µ 9 , C µ 10 , C µ 9 , C µ 10 are strongly or even very strongly favored by the data, with respect to all other combinations, even if frequentist measures of the goodness of fit like the minimum chi-squared per degree of freedom show a less pronounced preference.
Finally, we have also analyzed a few well-known BSM models that can provide a high energy framework for the EFT analysis. These include the exchange of a heavy Z gauge boson in models with heavy vector-like fermions and a scalar field whose vev breaks spontaneously the new symmetry, and a model with scalar leptoquarks. But despite the introduction of new constraints that are specific to the model-dependent analysis, when it comes to determining which hypotheses are strongly favored by the data, the Bayes factors mirror the results of the EFT fits, i.e., models that can generate the C µ 9 and C µ 9 Wilson coefficients after integrating out heavy degrees of freedom are preferred with respect to other combinations.

A List of observables used in the global analysis
In this appendix we provide a tabularized list of all the observables included in our global analysis as components of the likelihood function. For each of them we show the experimental measurement and the SM prediction derived with flavio, which includes the theoretical error obtained by calculating the spread of values for a given observable, when a set of input parameters (form factors, bag parameters, decays constants, masses of the particles) were randomly generated for 2000 times. In the last column we also present a deviation of the measurement from the SM prediction that quantifies the significance of a potential anomaly.