Global fits to b ->s ll data and signs for lepton non-universality

There are some slight tensions with the SM predictions within the latest LHCb measurements. Besides the known anomaly in one angular observable of the rare decay B ->K* mu+ mu-, another small discrepancy recently occurred. The ratio R_K = BR(B+ ->K+ mu+ mu-) / BR(B+ ->K+ e+ e-) in the low-q^2 region has been measured by LHCb showing a 2.6 sigma deviation from the SM prediction. In contrast to the anomaly in the rare decay B ->K* mu+ mu- which is affected by power corrections, the ratio R_K is theoretically rather clean. We analyse all the b ->s ll data with global fits and in particular explore the possibility of breaking of lepton universality. Possible cross-checks with an analysis of the inclusive B ->X_s l+ l- decay are also explored.


Introduction
There are some small tensions with the SM predictions within the latest LHCb measurements: The first measurement of new angular observables in the exclusive decay B → K * µ + µ − has shown a kind of anomaly [1]. Due to the large hadronic uncertainties it is not clear whether this anomaly is a first sign for new physics beyond the SM, or a consequence of underestimated hadronic power corrections [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]; but of course, it could just turn out being a statistical fluctuation. So the LHCb analysis based on the 3 fb −1 dataset is eagerly awaited to clarify the situation.
Besides this known anomaly in the angular analysis of the rare decay B → K * µ + µ − , another small discrepancy recently occurred. The ratio R K = BR(B + → K + µ + µ − )/ BR(B + → K + e + e − ) in the low-q 2 region has been measured by LHCb using the full 3 fb −1 of data, showing a 2.6σ deviation from the SM prediction [17]. In contrast to the anomaly in the rare decay B → K * µ + µ − which is affected by unknown power corrections, the ratio R K is theoretically rather clean.
The R K discrepancy has been addressed in a few recent studies [18][19][20][21][22]. In Ref. [19], two leptoquark models have been proposed to explain the discrepancy. In Ref. [21], R-parity violating supersymmetry is used to explain the R K anomaly together with lepton non-universal effects in the W R search observed by CMS. The authors of Ref. [20] have studied the R K result by performing a Bayesian statistical fit to the Wilson coefficients considering the data on some b → s transitions.
In this paper, we analyse the latest LHCb data within various global fits applying a frequentist statistical approach using the available data on all the relevant |∆B| = |∆S| = 1 processes in order to explore the R K and B → K * µ + µ − anomalies.
Moreover, we discuss the inclusive decay mode B → X s + − which is an important theoretically clean mode of the indirect search for new physics via flavour observables [23,24]; especially it allows for a non-trivial cross-check of the recent LHCb data on the exclusive mode [13].
This paper is organised as follows. In section 2 we describe our model independent analysis and highlight the main new theoretical and experimental inputs of the present study. In section 3 we present our results for the global fits for different sets of operators. The comparison with the inclusive decay mode is provided in section 4 and our conclusions are given in section 5.

Input of model-independent analysis
Compared to the analysis in Ref. [13], we have the following changes within the experimental and theoretical inputs: • We now also consider the recent LHCb measurements of B → K + − and in particular the observable R K [25]. Our theory analysis in the low-q 2 region is based on the method of QCD factorisation (QCDf). The factorisable and non-factorisable order one α s corrections are known from Refs. [26][27][28][29]. In our numerical analysis we have considered the aforementioned corrections following [30]. We have chosen the factorisation scheme according to [29] where one of the three independent B → K form factors, f + , is fixed by f + ≡ ξ K to all orders in perturbation theory. This f + form factor is taken from LCSR calculations of Ref. [15]. For the high-q 2 region we have followed Ref. [31], where instead of using LCSR extrapolated form factors and applying Isgur-Wise relation, we have used the three full form factors from lattice calculations [32] which significantly reduces the form factor uncertainties.
• For the B → K * µ + µ − angular observables, we use the lattice results for the form factors of Ref. [33] in the high-q 2 region thereby getting a significant reduction in the theoretical uncertainties. The theoretical uncertainties due to power corrections are estimated following Refs. [34,35]. In order to be conservative in our theoretical error estimation, we  have doubled the power correction budget in the overall error compared to the theoretical predictions in Ref. [36].
• The three loop QCD corrections for C 10 [37] as well as the NLO electroweak corrections [38] are included.
• We also consider the exclusive electron B → K * e + e − decay for the region where q 2 takes the whole [0.1, • For the inclusive modes B → X s + − we consider the electron and muon final states separately and for the experimental values we use the BaBar results [40].
• We have updated all the numerical values for the input parameters [41,42]. Some of the main input parameters which have been updated compared to [13] are given in Table 1.
• Compared to Refs. [13,46], where there is an 11% theoretical uncertainty for the SM prediction of BR(B s → µ + µ − ), we now have an overall error of about 7.5%. This reduction is due to the improved precision of the lattice results for the decay constant f Bs [42] as well as inclusion of higher order QCD and EW corrections 2 .
The list of all the observables that we have considered in this work together with the SM values and the experimental results are given in Tables 2 and 3. The theoretical predictions of all the observables are computed using the program SuperIso [53,54]. For more details on the theoretical framework we refer to Ref. [13,55]. We perform a model independent χ 2 analysis using all the observables given in Tables 2 and 3. For B → K * 0 µ + µ − angular observables we consider the experimental correlations as described in Ref. [13]. As one of the main objectives of this study is to investigate lepton universality and assess the effect of R K , we use the ∆χ 2 approach which is more suitable to find the preferred directions in the new physics parameter space. But we have to check first that the χ 2 method signals the overall consistency of the fit.
For each Wilson coefficient, one defines δC i = C NP i − C SM i . We consider separately the electron and muon semileptonic Wilson coefficients C e 9,10 and C µ 9,10 respectively, which are equal in the SM or in models with lepton universality.

Results
We first make a global fit to all observables considering new physics contributions to two Wilson coefficients only at a time. We then extend the study by considering new physics contributions to four Wilson coefficients and highlight the limitations of considering arbitrarily only a subset 1 Although it seems that using the experimental results of LHCb [39] ((3.1 ± 0.9) × 10 −7 ) available for the [0.0009, 1.0] GeV 2 bin and its corresponding SM prediction ((2.43 ± 0.57) × 10 −7 ) [2] one could get an observable with less (theoretical + experimental) error, but we use the full range of q 2 from BaBar since the [0.0009, 1.0] GeV 2 bin is dominated by C eff 7 . 2 See Ref. [47] for the various sources of uncertainties.

SM prediction Measurement
Observable SM prediction Measurement  [1,48], here the experimental errors have been symmetrised by taking the largest side error and whenever there have been several sources of uncertainty the total error has been obtained by adding them in quadrature.
Observable SM prediction Measurement 13.6 ± 2.0 9.5 ± 1.7 [52] 10 0.24 ± 0.07 0.60 ± 0.31 [40]  In the following we use all the observables given in Tables 2 and 3 in the global fits. It was shown in Ref. [18,19] that the current data on R K cannot be explained by tensor operators only. Moreover, the bounds from B s → µ + µ − disfavour also the possibility of scalar and pseudoscalar operators accounting for R K . However, a fine-tuned solution at the 2σ level remains possible if one assumes large electron contribution and one accepts cancellations in order to fulfil the B s → µ + µ − constraints. Therefore we consider here new physics contributions to the vector and axial vector operators allowing for flavour non-universality.

Fit results for two operators
We first consider new physics effects in O 9 and O 10 and make a χ 2 fit by scanning over δC 9 and δC 10 . The minimum χ 2 here is 52, with 52 degrees of freedom. The best fit point has therefore a correct value with respect to the goodness-of-fit. The result is presented in Fig. 1 where δC 9 and δC 10 are normalised to their SM values. We notice that deviations of the order of 40% compared to the SM predictions are allowed for C 9 and C 10 as a result of the global fit to all observables. However, the SM value (corresponding to δC 9 = δC 10 = 0) is slightly disfavoured by 2.3σ which represents a small tension for C 9 .
We then analyse further the effect of C 9 by considering separately the electron and muon contributions, and make a χ 2 fit by scanning over δC µ 9 , δC e 9 . We obtain a minimum χ 2 of 44, with 52 degrees of freedom, which shows that the fit is better than in the previous case.   , which shows that the SM value for the muon coefficient is disfavoured at more than 3σ. Yet, the universality condition, δC µ 9 = δC e 9 , is still barely allowed, at the border of the 2σ level, meaning that non-universality improves the fit. This shows clearly that in this fit the tension with the SM originates from the muon contribution.
Finally we consider the possibility of chirality flipped operator O 9 and scan over δC 9 , δC 9 . The minimum χ 2 is 52, with 52 degrees of freedom. The results are shown in Fig. 3. The tension with the SM is still present in this set but only for C 9 and not for C 9 .

Fit results for four operators
We expand here the study in the previous section by considering four operators in the fits, since there is a priori no reason that new physics should affect only two operators. We consider three possible sets including chirality flipped operators, and electron and muon contributions separately which are described in the following subsections.
3.2.1 Fit results for {C 9 , C 9 , C 10 , C 10 } The first set of four operators that we consider is an extension of {O 9 , O 10 } presented in section 3.1 by adding also the chirality flipped operators {O 9 , O 10 }. For this, we scan over δC 9 , δC 9 , δC 10 , δC 10 . In Fig. 4 we show the results of the global fit. We see that in this scenario the agreement with the SM is better but only at the 2σ level. Negative contributions to δC 9 are more favoured, whereas all the other Wilson coefficients can keep their SM values. To compare with the result of section 3.1, we superimpose the 1 and 2σ contours from the fit to two operators only in Fig. 4. This shows that considering only two operators leads to more restrictive results and some new physics contributions could be overlooked.
The set {C 9 , C 9 , C 10 , C 10 } is an extension of both the {C 9 , C 10 } and {C 9 , C 9 } sets, it is hence instructive to compare the χ 2 values of the best fit-points. The best fit point of {C 9 , C 9 , C 10 , C 10 } has a value of 51 (for 50 degrees of freedom), which is very similar to the values for the two other sets. Therefore, we can conclude that in the lepton universal scenario, adding the C 10 coefficients or the primed coefficients does not improve the fit. In the lower left plot the contours corresponding to 1 and 2σ regions from the fit to (C 9 , C 9 ) is superimposed.
In order to compare the results with the scan for two operators {O µ 9 , O e 9 } only, we overlay the contours corresponding to the 1 and 2σ fit result for {C µ 9 , C e 9 } in Fig. 5 in the plane (δC µ 9 , δC e 9 ). The comparison with the results obtained in the {C µ 9 , C e 9 } case clearly shows that considering only the modification of two Wilson coefficients leads to much more restrictive results.
We see that the SM is disfavoured at the 2σ level, even if the agreement is now improved compared to the two operator case. In particular, it is now possible to have simultaneously δC e 9 = 0 and δC µ 9 = 0 at the 2σ level. Yet there is still tension in the muon sector for C 9 . However, we emphasize also that sizeable new physics contributions to the other three operators are allowed at the 1σ level as it is visible in the plots of Fig. 5.
The set {C µ 9 , C µ 9 , C e 9 , C e 9 } is a direct extension of {C µ 9 , C e 9 }. Its best fit point has a χ 2 of 42, therefore adding the primed coefficients only slightly improves the fit. If we now compare {C µ 9 , C µ 9 , C e 9 , C e 9 } to the {C 9 , C 9 } set (which can be obtained by setting δC µ 9 = δC e 9 and δC µ 9 = δC e 9 ), we notice an improvement of about 2.6σ, which shows that considering nonuniversal lepton couplings improves the fit. The rigorous statement is the following: Assuming that the four operator scenario with {C µ 9 , C µ 9 , C e 9 , C e 9 } is correct, then the one with the two operators with {C 9 , C 9 } is ruled out at 2.6σ. This allows at least for a qualitative comparison of the various four operator fits. Figure 5: Global fit results for C µ 9 , C µ 9 , C e 9 , C e 9 . The red and black contours in the upper left plot correspond to the 1 and 2σ regions respectively, for the (C e 9 , C µ 9 ) fit presented in section 3.1.

3.2.3
Fit results for {C µ 9 , C e 9 , C µ 10 , C e 10 } We consider finally the set with both O 9 and O 10 but allow for lepton non-universality by considering separately the electron and muon contributions, neglecting the chirality flipped operators. We therefore perform a scan over δC µ 9 , δC e 9 , δC µ 10 , δC e 10 . The results are displayed in Fig. 6. Figure 6: Global fit results for C µ 9 , C e 9 , C µ 10 , C e 10 . The red and black contours correspond to the 1 and 2σ regions respectively of the (C e 9 , C µ 9 ) fit presented in section 3.1.
As in the other cases, having all the Wilson coefficients at their SM values is disfavoured at the 2σ level, yet allowing for a negative contribution to C µ 9 sorts this problem. However, as the plots illustrate, we note again that also large new physics contributions to C e 9 , C µ 10 and C e 10 are allowed within the 1σ level.
The comparison with the results obtained in the {C µ 9 , C e 9 } case clearly shows that considering only the modification of two Wilson coefficients leads to much restrictive results and overlooks other viable possibilities for new physics contributions.
The set {C µ 9 , C e 9 , C µ 10 , C e 10 } is a direct extension of {C µ 9 , C e 9 }. The best fit point in that scenario has a χ 2 of 43. So the addition of the C µ,e 10 coefficients only slightly improves the fit. If we compare it to the {C 9 , C 10 } set, we see an improvement of 2.5σ, again favouring the non-universal extension against the universal one.

Fit results assuming left-handed leptons
Finally, we make the assumption that we have left-handed leptons only, which represents an attractive option in model building beyond the SM. In this case one finds the following relations between the Wilson coefficients: Here we introduced the quantities C i XY where i is the flavour index, X denotes the chirality of the quark current and Y of the lepton one. 3 Figure 7: Global fit results for C µµ LL , C ee LL , C µµ RL , C ee RL .
We perform a fit for C µµ LL , C ee LL , C µµ RL , C ee RL . The results are shown in Fig. 7. The fit shows a large 2σ region and three separate 1σ areas. The coefficient C µµ LL deviates slightly from its SM value at the 2σ level, while the other coefficients are compatible with their SM values. Both C µµ LL and C µµ RL can vary by only 30% away from the SM value, while the electron operators can have a larger deviation. The minimum χ 2 of the fit is 41.
We also made a fit with two coefficients C µµ LL , C ee LL , setting the RL operators to zero. The result is displayed in Fig. 8. As seen in section 3.2, the resulting fit regions are much smaller, which shows that it is important to consider the possibility of RL new physics operators as well. In this case, the minimum χ 2 is 46, showing that assuming the four operator set is the correct one, the two operator set is disfavoured by about 1.5σ. 3 For completeness, we note that assuming right-handed leptons only, one gets the following relations:

Comparison of exclusive and inclusive b → s observables
The inclusive mode B → X s + − can only be measured at e + e − machines and is theoretically cleaner than the exclusive modes [56,57]. The theoretical accuracy in the low-q 2 region is of the order of 10% [58]. But the branching fraction has been measured by Belle and BaBar using the sum-of-exclusive technique only. The latest published measurement of Belle [59] is based on a sample of 152 × 10 6 BB events only, which corresponds to less than 30% of the dataset available at the end of the Belle experiment. BaBar has just recently presented an analysis based on the whole dataset using a sample of 471 × 10 6 BB events [40] overwriting the previous measurement from 2004 based on 89 × 10 6 BB events [60].
In order to compare different sets of observables, we consider the operators O 7 , O 9 and O 10 which are the most relevant operators for the semileptonic B decays. Again we use the ∆χ 2 fit method to obtain the exclusion plots of the Wilson coefficients. 4 First we make the global fit using only B → X s + − branching ratios in the low-q 2 and high-q 2 regions with = µ, e. We note that there is no sign of lepton non-universality in the published data. We find that the χ 2 fit using the two branching ratios with = e only is not very good; there is no compatibility at 68% C.L. , so for the ∆χ 2 -metrology we restrict ourselves to the case = µ. We compare this with two other fits using exclusive observables, one with BR(B → K 0 µ + µ − ), BR(B + → K + µ + µ − ), R K , and one with all B → K * µ + µ − observables. In Fig. 9, we illustrate the results of the ∆χ 2 fit for the relevant Wilson coefficients. The upper row shows the fit based on the exclusive observables BR(B → K 0 µ + µ − ), BR(B + → K + µ + µ − ), R K ; the middle row shows fit based on B → K * + − observables; and the lower row the one based on the measurements of the inclusive (B → X s µ + µ − ) branching ratio in the low-and high-q 2 regions. It is remarkable that all three sets of exclusion plots are nicely compatible with each other. This is a non-trivial consistency check.
At the moment, the measurements of the B → K * + − are the most powerful ones. However, the final word of Belle is still pending and, moreover, there will be a Super-B factory Belle-II with a final integrated luminosity of 50 ab −1 [61]. There is a recent analysis [62] of the expected total uncertainty on the partial decay width and the forward-backward asymmetry in several bins of dilepton mass-squared for the fully inclusive B → X s + − decays assuming a 50 ab −1 total integrated luminosity. Based on some reasonable assumptions (for details see Ref. [13]), one finds a relative fractional uncertainty of 2.9% (4.1%) for the branching fraction in the low-(high-) q 2 region and a total absolute uncertainty of 0.050 in the low-q 2 bin 1 (1 < q 2 < 3.5 GeV 2 ), 0.054 in the low-q 2 bin 2 (3.5 < q 2 < 6 GeV 2 ) and 0.058 in the high-q 2 interval (q 2 > 14.4 GeV 2 ) for the normalised A F B . Hence, the inclusive mode will lead to very strong constraints on the Wilson coefficients and to a more significant cross-check of the new physics hypothesis. Figure 9: Fit results for the new physics contributions to C 7 , C 9 and C 10 , using only BR(B → K 0 µ + µ − ), BR(B + → K + µ + µ − ), R K (first row), using only B → K * µ + µ − observables (second row), and using the current measurements of BR(B → X s µ + µ − ) at low-and high-q 2 (third row).
We illustrate the usefulness of these future measurements of the inclusive mode at Belle-II in the following way: We make a model independent fit for the coefficients C 7 , C 8 , C 9 , C 10 and C l (for notation see Ref. [13]). In addition to the observables given in Tables 2 and 3, we consider the inclusive branching ratio of B → X s γ as well as the isospin asymmetry in B → K * γ decay which are relevant to constrain C 7 and C 8 . Based on our model-independent analysis we predict the branching ratios at low-and high-q 2 . In Fig. 10, we show the 1, 2, and 3σ ranges for these observables. In addition, we add the future measurements at Belle-II assuming the best fit solution of our model-independent analysis as central value. These measurements are indicated by the black error bars. They should be compared with the theoretical SM predictions given by the red (grey) error bars. Fig. 10 indicates that the future measurement of the inclusive branching ratios separates nicely the SM prediction and the model-independent best fit point. Moreover, the future measurement of the forward-backward asymmetry at Belle-II will allow to separate the potential new physics measurement from the SM prediction in a significant way as shown in Fig. 11. Figure 10: 1, 2 and 3σ predictions for the branching ratio at low-and high-q 2 within the modelindependent analysis. Future measurement at the high-luminosity Belle-II Super-B-Factory assuming the best-fit point of the model-independent analysis as central value (black) and the SM predictions (red/grey). Figure 11: 1, 2 and 3σ predictions for the unnormalised forward-backward asymmetry in bin 1 (1 < q 2 < 3.5 GeV 2 ) and in bin 2 (3.5 < q 2 < 6 GeV 2 ) within the model-independent analysis. Future measurement at the high-luminosity Belle-II Super-B-Factory assuming the best-fit point of the model-independent analysis as central value (black) and the SM predictions (red/grey).

Conclusions
We present here for the first time global fits to the complete b → s dataset available from B factories and from LHC, in particular addressing the two observed tensions in the angular analysis of the exclusive decay B → K * µ + µ − and the ratio R K . We perform several global fits using different sets of two or four vector and axial vector operators allowing for non-universality.
Comparing the 4-operator and 2-operator fit results we have shown that while considering 2-operator fits can be illustrative for where in the parameter space new physics could be found it can be too restrictive and maybe even misleading since a large area of new physics parameter space might be unjustifiably overlooked.
Considering the full set of 8 semileptonic operators {C µ 9 , C µ 9 , C e 9 , C e 9 , C µ 10 , C µ 10 , C e 10 , C e 10 } and comparing the minimum χ 2 of the different subsets can lead to a determination of the most relevant operators. This comparison reveals that the sets with lepton non-universality, namely {C µ 9 , C µ 9 , C e 9 , C e 9 } and {C µ 9 , C e 9 , C µ 10 , C e 10 } give the best fit to the data. Our global fits show that simultaneous agreement of all the Wilson coefficients with the Standard Model is ruled out by 2σ at least. This problem can be resolved if the new physics contribution to C µ 9 is negative. However, we emphasize that sizeable new physics contributions to the other operators are allowed at the 1σ level.
If these tensions are not resolved in the near future, we have demonstrated that the future measurements of the inclusive b → s observables by Belle II will allow for a powerful crosscheck. The present data on inclusive and exclusive decays are nicely compatible with each other and there is no sign of lepton non-universality in the published data on the inclusive mode.