Implications of new evidence for lepton-universality violation in $b\to s\ell^+\ell^-$ decays

Motivated by renewed evidence for new physics in $b \to s\ell\ell$ transitions in the form of LHCb's new measurements of theoretically clean lepton-universality ratios and the purely leptonic $B_s\to\mu^+\mu^-$ decay, we quantify the combined level of discrepancy with the Standard Model and fit values of short-distance Wilson coefficients. A combination of the clean observables $R_K$, $R_{K^*}$, and $B_s\to \mu\mu$ alone results in a discrepancy with the Standard Model at $4.0\sigma$, up from $3.5\sigma$ in 2017. One-parameter scenarios with purely left-handed or with purely axial coupling to muons fit the data well and result in a $\sim 5 \sigma$ pull from the Standard Model. In a two-parameter fit of %$C_9$ and $C_{10}$, new-physics contributions with both vector and axial-vector couplings to muons the allowed region is much more restricted than in 2017, principally due to the much more precise result on $B_s \to \mu^+ \mu^-$, which probes the axial coupling to muons.Including angular observables data restricts the allowed region further.A by-product of our analysis is an updated average of $\text{BR}(B_s \to \mu^+ \mu^-) = (2.8\pm 0.3) \times 10^{-9}$.


INTRODUCTION
Flavor physics played a central role in the development of the Standard Model (SM) and could well spearhead the discovery of new physics (NP) beyond the SM (BSM). In fact, although the vast majority of particle-physics data is consistent with the predictions of the SM, a conspicuous series of discrepancies has appeared in rare flavor-changing processes mediated by quark-level b → s transitions. These are suppressed by the "GIM mechanism" in the SM and are, therefore, potentially sensitive to very high-energy NP scales [1]. A perennial question in this context is how to distinguish longdistance strong-interaction effects from genuine new physics. Several years ago, following LHCb's first measurement of the lepton-universality violating ratio R K * , we demonstrated [2] the power of using observables which are almost entirely free from hadronic uncertainties to provide a high-significance rejection of the SM, and its potential to narrow down the chiral structure of the BSM effect. In particular, we pointed out the importance of the B s → µ + µ − decay and lepton-flavorviolating ratios of forward-backward asymmetries in lifting a degeneracy between axial and vectorial couplings to leptons. Motivated by LHCb's updates to the ratio R K and of BR(B s → µ + µ − ) we revisit this set of decays in the present work.

OBSERVABLES
The measurement of several rare b → s decays yields results in tension with the SM expectations implying the presence of new interactions breaking lepton universality (see Refs. [1,3] for recent reviews). Among them stand out a subset of observables with theoretical uncertainties at or below the percent level [4,5]. Our main observables of interest com-prise the lepton-universality ratios R K = Γ(B → Kµµ)/Γ(B → Kee) [4] and R K * = Γ(B → K * µµ)/Γ(B → K * ee), and the purely leptonic decay B s → µ + µ − .
In particular, the LHCb collaboration has just reported the most precise measurement of R K in the q 2 -bin [1.1, 6] GeV 2 using the full run 1 and 2 datasets [6], where the first uncertainty is statistical and the second systematic. This result deviates from the SM predictions (see Table I in Ref. [2]) 1 R SM K = 1.0004 +0.0008 −0.0007 , with a significance of 3.1σ. Compared to the first LHCb measurement reported in 2014 [9], the tension with respect to the SM has significantly increased. At the same time, LHCb has published new results for the branching faction of B s → µ + µ − , obtained with the same full dataset [10] and also known to about 1% accuracy in the SM. This result, along with other recent measurements done by ATLAS [11] and CMS [12], indicate a decay rate lower than the SM prediction. This set of key inputs is summarized in Table I. In place of the asymmetric errors on R K and R K * published by the experiments, we conservatively employ a symmetric error equal to the upper, larger error (combining statistical and systematic in quadrature), in line with the treatment in Ref. [2]. In 2020 LHCb also reported a new measurement of the CP-averaged angular observables of the decay B 0 → K * 0 µ + µ − [17] and of its isospin partner, B + → K * + µ + µ − [18]. The new data seems to confirm the previous measurements pointing to possible tensions with the SM . However, and contrary to the lepton-universality ratios and B s → µµ, the SM predictions for the B → K * µµ angular observables suffer from significant hadronic uncertainties which hinder a clear interpretation of the discrepancies in terms of NP [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58].
In this work we combine the experimental data focusing on the clean observables as in Ref. [2] and carry out global fits of the Wilson coefficients (or short-distance coefficients) of the low-energy b → s effective Lagrangian to the data. We find that the data on clean observables is at variance with the SM at a level of 4.0σ. We also find that one-parameter scenarios with purely left-handed or axial currents provide a good description of the data, excluding the SM point in each case at close to 5σ. As discussed abundantly in the literature, such new lepton-universality-violating (LUV) interactions can arise at tree or loop level from new mediators such as neutral vector bosons (Z ) or leptoquarks (see Ref. [59] which includes a review of NP interpretations).
An important aspect to note is that the three measurements of B s → µ + µ − cannot be naively averaged together, as a result of correlations with B d → µ + µ − . We therefore construct a two-dimensional joint likelihood from the published measurements [10][11][12]. In doing so, we assume a correlation coefficient of −0.5 for ATLAS, which reproduces the results reported in Ref. [11], and neglect correlations in the LHCb measurement. The resulting combination is represented in Fig. 1.
with χ 2 min = 3.72 (5 d.o.f.). 2 As with the existing combination [60], the central value of the average is lower than the average of the three individual central values.
We combine the experimental measurements and the SM prediction of the B s → µ + µ − branching fraction in the ratio obtaining R = 0.78(9) by using the most up to date theoretical prediction of Ref. [5]. We end this section by noting that only the recent LHCb result implements a newer (and larger by ∼ 6%) measurement of the ratio of hadronization fractions f s / f d . However, including the corresponding increase in the branching fractions of B s → µ + µ − measured by ATLAS and CMS (keeping the correlation with B d → µ + µ − ) leads to a very small increase (of about ∼ 3%) in the average in Eq. (4) that will be neglected in this work.

THEORETICAL APPROACH
The low-energy effective Hamiltonian for semileptonic b → s + − processes at the scale µ ∼ m b in the SM is written as [61] where G F is the Fermi constant, and λ ps = V pb V * ps is a combination of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements with p = u, c. The short-distance contributions, stemming from scales above µ ∼ m b , are matched to a set of Wilson coefficients C i . The O p 1,2 , O 3−6 and O 8 are the "currentcurrent", "QCD-penguin" and "chromomagnetic" operators, [11], CMS [12], and LHCb [10], compared to the SM prediction (red square). Contours of the combination correspond to 1σ, 2σ and 3σ, and those of each experiment to just 3σ.
respectively. The explicit forms for these operators can be found in Ref. [61]. The remaining operators O 7 , O 9 , and O 10 from electromagnetic penguin-, electroweak penguin-, and box-loop diagrams are defined as follows: In the presence of NP, one has nine more operators, i.e., O 7 and O 9,10 , the opposite-quark-chirality counterparts of O 7 and O 9,10 , plus four scalar and two tensor operators [62]. However, if the NP effect enters in the couplings to muons, only O 9 and O 10 can explain the R K ( * ) data [62]. The tensor-operator contributions to b → s decays are suppressed if the NP scale is heavier than the electroweak scale. Other operators cannot induce LUV, or are tightly constrained by the B q → decays. In the following analysis, we assume that all the Wilson coefficients are real and that the presence of NP only appears in the b → sµ + µ − sector. The last assumption is justified by the tension in the data in For reliable predictions of observables, it is of vital importance to estimate theoretical uncertainties. They mainly stem from nonperturbative contributions including form factors F(q 2 )'s and "nonfactorizable" terms h λ (q 2 )'s [40]. At low q 2 the nonperturbative contributions can be addressed in the heavy quark and large-energy limits: they can be expanded as [51], where F ∞ (q 2 ) and h ∞ λ (q 2 ) can be calculated in light-cone sum rules [43,63] and within the QCD factorization approach [40], and the rest are power-correction terms. Among them, r λ (q 2 ) is dominated by the long-distance charm contributions involving the "current-current" operators O c 1 and O c 2 [51]. We parametrize the charm loop contributions, r c λ (q 2 ), by where A λ and B λ are dimensionless constants. Therefore, the overall uncertainties at low q 2 arise from the leading-power terms, power corrections of form factors and charm loop contributions, characterized by 27 nuisance parameters. For more details and ranges of values taken for these, see Refs. [2,51].
In the high-q 2 region, form factors have been calculated in lattice QCD [48] whereas the nonfactorizable contribution can be computed with an operator product expansion [41,44]. We omit the analysis of this region as none of the LUV ratios have been measured there yet. We use the frequentist statistical approach to quantify the compatibility between the experimental data and the theoretical predictions. We define the χ 2 function as where χ 2 exp ( C, y) includes the correlations reported by the experiments and χ 2 th ( y) is a theoretical component. The theoretical predictions for the observables O th are functions of Wilson coefficients C and nuisance hadronic parameters y. We choose two models for χ 2 th ( y); one in which y i follows a normal distribution (that we call "Gaussian") and another (that we call "R-fit") where it is restricted to a range, see Ref. [51], with a flat distribution. Here, we assume that these nuisance parameters are uncorrelated [51].
In order to obtain best-fit values in a particular scenario, we can construct a profile χ 2 depending only on certain Wilson coefficients with the remaining Wilson coefficients set to their SM values.
Here, χ 2 is minimized by varying the nuisance parameters.
In our statistical analysis, we adopt the widely used p-value and Pull SM to denote how well the experimental data can be described and how significant is the deviation from the SM. Results of fits are reported for δC i ≡ C i − C SM I .

THE THEORETICALLY CLEAN FIT
We first restrict ourselves to the analysis of the theoretically clean observables R K , R K * and B s → µ + µ − . The relevant data  Contour plots at the 1σ, 3σ and 5σ confidence levels in the (δC µ 9 , δC µ 10 ) plane for the Gaussian form of χ 2 th (regions in light red and orange). We also show the 1σ and 3σ contours for the individual constraints (R K ( * ) and B s → µµ) and, for the same fit from 2017 [2] (dashed in red).
Minimizing χ 2 over all SM and theoretical nuisance parameters, one obtains χ 2 min,SM = 31.32 (χ 2 min,SM = 30.54) and p SM = 5.4 × 10 −5 (p SM = 7.6 × 10 −5 ) using the Gaussian (R-fit) form of the χ 2 th . We emphasize that the sole role of the minimization here is to implement the tiny theoretical uncertainties. In arriving at this p-value and in the rest of this paper, we treat our BR(B s → µ + µ − ) average as a single measurement, following common practice and in line with Ref. [2]. In other words, with 7 d.o.f. the clean data is at variance with the null hypothesis (Standard Model) at a level of 4.0σ, up from 3.5σ in our previous work [2]. Were we to treat the 6 B q → µ + µ − measurements as separate inputs, we would instead obtain χ 2 min,SM = 35.04 and p SM = 3.5σ (12 d.o.f.); the reduction in significance is unsurprising given that we now include (nonanomalous) data on b → d transitions. We next fit one-and two-parameter BSM scenarios. Only δC µ 9 and δC µ 10 can describe a deficit in both R K and R K * (Ref. [2] and Fig. 1 there). Therefore we analyze the data in the clean observables by fitting only these two Wilson coefficients. Figure 2 shows the constraints imposed in the (δC µ 9 , δC µ 10 ) plane and using the Gaussian model for χ 2 th (see below for the Rfit results). In Table III we also show the numerical results for this fit, as well as fits involving a single Wilson coefficient, for the Gaussian approach. We also fit the left-and right-handed combinations (related to the couplings to muons) of Wilson coefficients δC µ L = (δC µ 9 − δC µ 10 )/2 and δC µ R = (δC µ 9 + δC µ 10 )/2. We observe that one-parameter scenarios with purely lefthanded coupling δC µ L or purely axial coupling δC µ 10 describe the data well (p > 10%). Compared to either scenario, the SM (identified as the δC i = 0 point) is excluded at a confidence level of close to 5σ. On the other hand, and in contrast with our previous analysis from 2017 [2], the pure-C 9 scenario is in tension with the data at 2.3σ although it still compares favorably with the data compared to the SM. Finally, as we will see below, using the R-fit version of χ 2 th and increasing the theoretical uncertainties in the predictions of B → K ( * ) only produce very small changes in the results of the fit to clean observables and do not change our conclusions.

THE GLOBAL FIT
For the sake of completeness we also perform a global fit including all the measurements of angular observables reported by the LHCb, ATLAS, and CMS experiments in the low-q 2 region. As mentioned above, these observables are afflicted by larger theoretical uncertainties compared to LUV ratios and B s → µµ. However, it is important to analyze how the conclusions change when including these data within a modelindependent framework for the theoretical uncertainties such as ours.  More specifically, compared to our 2017 analysis [2], we replace the CP-averaged angular observables for the B 0 → K * 0 µ + µ − , the ratio R for the B s → µ + µ − decay, and R K in the bin [1.1, 6.0] GeV 2 with the latest measurements by the LHCb, CMS and ATLAS experiments [11,12,17,60,64]. In addition, we also include 32 new measurements of F L , P 1 , P 2 , P 3 , P 4 , P 5 , P 6 and P 8 in four low bins (q 2 ≤ 6 GeV 2 ) for the B + → K * + µ + µ − decay [18] as well as three Belle data R K and R K * in Table I. As a result the total number of data fitted becomes 94. 3 In this fit strategy, we obtain a χ 2 min,SM = 126.88 and p SM = 0.01 with 94 d.o.f.. Compared to the global fit results in Ref. [2], the updated fit results in Table III and Fig. 3 show that the confidence level of the exclusion of the SM point increases by 1.3σ for the δC µ 10 scenario and by 0.7σ for the δC µ L scenario, but only by 0.2σ and 0.4σ in the δC µ 9 and (δC µ 9 , δC µ 10 ) scenarios, respectively. Interestingly, these updated fits constrain better these Wilson coefficients and exclude positive δC µ 9 and negative δC µ 10 at more than the 3σ confidence level. We also perform a four-dimensional global fit with the Gaussian χ 2 th including C µ 9 and C µ 10 . The resulting Wilson coefficients from the fit are, with the correlation matrix, and where χ 2 min = 96.88 for 90 d.o.f., corresponding to a pvalue of 0.29 and a Pull SM = 4.57.

IMPACT OF THEORETICAL UNCERTAINTIES
Finally, we briefly investigate the robustness of the fits with respect to the hadronic uncertainties. We do so by comparing the results obtained above with those obtained by using the Rfit model for χ 2 th , with nominal hadronic uncertainties in B → K ( * ) or multiplied by a factor 2 and 3. The relevant results are shown in Tables IV, V and VI. In Fig. 4 we also show the new results in the (δC µ 9 , δC µ 10 ) plane overlaid with the ones obtained with the Gaussian model of χ 2 th . The treatment of uncertainties has a significant impact on the global fit, especially in the parameter ranges obtained for the Wilson coefficient δC µ 9 . As discussed in Refs. [47,51] this is due to the fact that a shift to C 9 in the amplitude is indistinguishable from a nonfactorizable charm contribution or a shift to a certain combination of B-decay form factors. Therefore, increasing the ranges allowed for these parameters in a framework such as R-fit tends to reduce the significance of a NP effect in C 9 .
This effect is clearly seen in Fig. 4 where the contours in the global fit approach those of the clean fit when increasing the errors and the tension of the data with the SM becomes dominated by the LUV ratios and B s → µµ. In contrast, the results and conclusions derived from the clean fit to these latter observables are robust with respect to the same variation of hadronic uncertainties. This is illustrated in Fig. 5.

SUMMARY AND OUTLOOK
In conclusion, we have presented a statistical analysis of recent data on LUV ratios and B s → µµ using the low-energy b → s effective Lagrangian. We find that the data on these clean observables disagree with the SM at a level of 4.0σ. Scenarios with pure left-handed or axial currents provide a good description of the data, and each of them excludes the SM point at ∼ 5σ confidence level. Therefore, our results reinforce the NP interpretation of the anomalies in the b → s transitions, which could correspond to the contribution of the tree-level exchange of a leptoquark or Z with a mass Λ ∼ 30 TeV and ∼ O(1) couplings to the SM. Further data on LUV observables from LHCb and Belle II should be able to clarify soon this tantalizing possibility.  In this section, we show that the latest measurements of the branching fraction ratios R K 0 s and R K * + , in combination with the already known ratios R K + , R K * 0 , and the branching fraction B s → µµ, point at a discrepancy with the Standard Model at 4.2σ. One-parameter scenarios, C µ 10 and C µ L , fit the data well and result in a pull from the SM crossing the 5.0σ threshold. The two-parameter fit of C 9 and C 10 yields also a pull of 5.02σ. On the other hand, the one-parameter scenario of C 9 alone still has a pull of 4.5σ, up from the previous 4.1σ (see Table II).
Very recently, the LHCb Collaboration reported measurements of two new lepton-universality ratios and R K * + = Γ(B + → K * + µ + µ − )/Γ(B + → K * + e + e − ) in the q 2 ranges [1.1, 6.0] GeV 2 and [0.045, 6.0] GeV 2 , respectively, using proton-proton collision data corresponding to an integrated luminosity of 9 fb −1 [73]. They represent the first observation of the B 0 → K 0 s e + e − and B + → K * + e + e − decays. The two ratios are where the first error is statistical and the second is systematic. These results show again tension with respect to the SM predictions (see Table I) with a significance of 1.5σ and 1.4σ, respectively. In the following analysis, we conservatively employ a symmetric error equal to the larger error. It should be noted that these ratios are the isospin partners of R K + and R K * 0 , and therefore should receive the same NP contributions, if they exist. As a result, it is of utmost importance to update the clean fit performed in the main text and to check whether these new measurements increase or decrease the significance of the tension with the SM.
Note that compared to the clean fit performed in the main text, the total number of fitted data becomes 9 after adding the two new LHCb measurements. Setting the Wilson coefficients to their SM values, we obtain a χ 2 min,SM = 36.50, corresponding to a p-value of 3.23×10 −5 or a 4.2σ deviation (for 9 d.o.f.). As a result, with the new data, the LU ratios are in discrepancy with the SM predictions at a level of 4.2σ, up from the 4.0σ in March 2021.
In Table VII and Fig. 6, we show the results of four fits which are obtained by allowing for lepton-specific contributions δC µ 9 , δC µ 10 , δC µ L and a two-dimensional scenario (δC µ 9 , δC µ 10 ). We find that the significance of all the NP scenarios except for δC µ 9 are more than 5σ. Such a result, obtained by considering only clean observables, seems to point unambiguously to the presence of new physics. We conclude that although the newly observed ratios are the isospin partners of the existing ones, they indeed have increased the significance of the tension with the SM by about 0.5σ compared with the fits without them in March 2021.  Fig. 2 but for the new clean fit to 9 data including the latest LHCb measurements [73] vs. the clean fit to 7 data shown in Fig. 2 (Gaussian in 2021).