New physics in rare B decays after Moriond 2021

The anomalies in rare B decays endure. We present results of an updated global analysis that takes into account the latest experimental input – in particular the recent results on \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}RK and BR\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(B_s \rightarrow \mu ^+\mu ^-)$$\end{document}(Bs→μ+μ-) – and that qualitatively improves the treatment of theory uncertainties. Fit results are presented for the Wilson coefficients of four-fermion contact interactions. We find that muon specific Wilson coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_9 \simeq -0.73$$\end{document}C9≃-0.73 or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_9 = -C_{10} \simeq -0.39$$\end{document}C9=-C10≃-0.39 continue to give an excellent description of the data. If only theoretically clean observables are considered, muon specific \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{10} \simeq 0.60$$\end{document}C10≃0.60 or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_9=-C_{10} \simeq -0.35$$\end{document}C9=-C10≃-0.35 improve over the Standard Model by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{\Delta \chi ^2} \simeq 4.7\sigma $$\end{document}Δχ2≃4.7σ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{\Delta \chi ^2} \simeq 4.6\sigma $$\end{document}Δχ2≃4.6σ, respectively. In various new physics scenarios we provide predictions for lepton flavor universality observables and CP asymmetries that can be tested with more data. We update our previous combination of ATLAS, CMS, and LHCb data on BR\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(B_s \rightarrow \mu ^+\mu ^-)$$\end{document}(Bs→μ+μ-) and BR\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(B^0\rightarrow \mu ^+\mu ^-)$$\end{document}(B0→μ+μ-) taking into account the full two-dimensional non-Gaussian experimental likelihoods.


Introduction
Since several years there exist persistent discrepancies between the Standard Model (SM) predictions and the experimental results for rare decays of B mesons based on the neutral current b → s transitions. Discrepancies are seen in the branching fractions of the rare decays B → K μ + μ − , B → K * μ + μ − , B s → φμ + μ − , and B s → μ + μ − , in the angular distribution of B → K * μ + μ − and in lepton flavor universality (LFU) ratios. Of particular interest are the hints for LFU violation that have been observed by LHCb in the following ratios of branching fractions a e-mail: waltmann@ucsc.edu b e-mail: stangl@itp.unibe.ch (corresponding author) While the SM predictions for most absolute branching fractions and also the angular observables are potentially subject to large hadronic uncertainties, the LFU ratios R K and R K * can be predicted with high accuracy. Significant deviations in these observables would thus constitute clean indirect evidence for new physics. Also the absolute branching ratio of the purely leptonic decay B s → μ + μ − can be considered as theoretically clean. Non-perturbative physics enters through a single hadronic parameter, the B s meson decay constant, which is know with high precision from lattice QCD calculations. Intriguingly, the simplest new physics scenarios that address the theoretically clean hints for LFU violation simultaneously explain also the other discrepancies. Parameterizing the new physics in terms of four fermion contact interactions, global fits of rare B decay data find consistently very strong preference for new physics in the form of the operator 1 2 NP (sγ α P L b)(μγ α μ) or 1 2 NP (sγ α P L b)(μγ α P L μ) with a generic new physics scale of NP 35 TeV (for recent work see [1][2][3][4][5][6][7][8][9][10]).
Very recently, the LHCb collaboration presented updated results for two theoretically clean observables that have previously shown tensions with the SM predictions: the LFU ratio R K and the branching ratio BR(B s → μ + μ − ). Using the full run 2 data set the value for R K is [11] R K = 0.846 +0.042 −0.039 +0.013 −0.012 , for 1.1 GeV 2 < q 2 < 6 GeV 2 , where the first uncertainty is statistical and the second one systematic, and q 2 is the di-muon invariant mass squared.
The new result has exactly the same central value as the previous result R K = 0.846 +0.060 −0.054 +0.016 −0.014 [12], while the statistical uncertainty has been reduced by approximately 30%, commensurate with the increased statistics. Consequently, the tension between the experimental measurement and the SM prediction, which is unity to an excellent approximation, has increased from previously 2.5σ to now 3.1σ .
The branching ratio of the B s → μ + μ − decay measured with the full run 2 data is found to be [13,14] BR(B s → μ + μ − ) = 3.09 +0.46 −0.43 +0.15 −0.11 × 10 −9 , where the first uncertainty is statistical and the second one systematic. This result by itself has a precision close to the previous world average BR(B s → μ + μ − ) = (2.69 +0.37 −0.35 ) × 10 −9 [15] that was based on results from ATLAS, CMS, and LHCb [16][17][18]. Compared to the previous measurement by LHCb, BR(B s → μ + μ − ) = (3.0 ± 0.6 +0. 3 −0.2 ) × 10 −9 [16], the new update finds nearly the same central value. While the LHCb result is compatible with the SM prediction within 1σ , the previous world average was below the SM prediction by more than 2σ . Here, we provide an update of the world average of the B s → μ + μ − branching ratio and the correlated B 0 → μ + μ − branching ratio, taking into account the new LHCb results. A Gaussian approximation to our combined two-dimensional likelihood is given by with an error correlation coefficient ρ = −0.27. We find a one-dimensional pull with the SM predictions of 2.3σ . Details on how the combination and the discrepancy with the SM are obtained are given in the Appendix A.
The main goal of this paper is to interpret the impact of the new experimental results in a model independent way, using the well established effective Hamiltonian approach. We parameterize new physics contributions by Wilson coefficients of dimension 6 interactions evaluated at the renormalization scale μ = 4.8 GeV We consider the following set of semi-leptonic operators We do not consider semi-leptonic tensor operators, because they are not generated at dimension 6 in the Standard Model Effective Field Theory (SMEFT). Similarly, in the case of the scalar operators, we will impose the following relations among the corresponding Wilson coefficients C bs S = −C bs P and C bs S = C bs P , as they hold at dimension 6 in the SMEFT [19]. We also do not consider semi-tauonic operators or 4-quark operators, as they affect the observables we consider only at the loop level [20,21].
A critical aspect of global fits is the treatment of theory uncertainties. In our previous studies [6,[22][23][24] we have evaluated theory uncertainties and their correlations for the Wilson coefficients fixed to their SM values. This is typically a good approximation as long as the best fit results are in the vicinity of the SM point. Possible exceptions are observables that have negligible uncertainties in the SM but not in the presence of new physics contributions. A prominent example of such observables are the LFU ratios R K and R K * . While the experimental uncertainties still dominate for R K and R K * , the precision of the new R K result in Eq. (2) is strong motivation to improve our treatment of theory uncertainties. In this paper we incorporate the new physics dependence of the theory uncertainties for the first time in our fit.
The paper is organized as follows: in Sect. 2 we discuss in detail the improved treatment of theory uncertainties and illustrate the size of the effect in the case of LFU observables and CP asymmetries in the presence of new physics. In Sect. 3 we collect the results of the updated global fit. We consider scenarios with one real Wilson coefficient at a time, scenarios with two real Wilson coefficients as well as scenarios with complex Wilson coefficients. Sect. 4 contains new physics predictions for a number of LFU observables and CP asymmetries that can be tested with future data. We conclude in Sect. 5. Our combination of the experimental results on the B s → μ + μ − branching ratio is described in Appendix A.

Improved treatment of theory uncertainties
Our global fits are based on a χ 2 function that depends on the Wilson coefficients in the effective Hamiltonian and that takes into account both the theoretical and experimental uncertainties in terms of covariance matrices, exp and th In the above expression, the O exp are the measured central values of the observables of interest and O th are the corresponding theory predictions that have dependence on the considered set of Wilson coefficients C i . The covariance matrix th includes uncertainties from parametric input, in partic-ular CKM matrix elements and form factor parameters, as well as from non-factorisable power corrections. Our treatment of the non-factorisable power corrections follows [22] and is summarized in Appendix B.1. In previous global fits, we made the assumption that the theoretical uncertainties are well described by the covariance matrix th determined with the SM values for the Wilson coefficients and neglected possible dependence of th on the new physics. This has the advantage that the time consuming evaluation of th has to be performed only once. We have developed a computationally efficient method to determine the new physics dependence of th . The procedure is summarized in the following and described in detail in Appendix B.2.
As the rare B decay amplitudes are linear functions of the Wilson coefficients it is possible to express the branching ratios as second order polynomials in the Wilson coefficients. The coefficients of the polynomials are independent of new physics and their correlated uncertainties can be described by a covariance matrix that needs to be determined only once. The covariance matrix of the branching ratios can then be expressed in a straight forward way in terms of the covariance matrix of the polynomial coefficients and the Wilson coefficients.
The CP averaged angular observables S i , the CP asymmetries A i , and the LFU ratios can be written in terms of ratios of second order polynomials, while the P i observables involve also irrational functions. In those cases we obtain an approximation of the covariance matrix for the observables by expanding the functions to second order in the Wilson coefficients and then following the same procedure as for the branching ratios. We find that this procedure gives reliable estimates as long as the absolute values of the new physics Wilson coefficients are somewhat smaller than the corresponding relevant SM coefficients. In principle, the accuracy of the approximation could be systematically improved by expanding to higher orders.
The new error treatment is particularly relevant for quantities that are predicted with very high precision in the SM but that have non-negligible uncertainties in the presence of new physics. In that case, the corresponding entries in the theoretical covariance matrix evaluated in the SM and the ones in the presence of new physics may differ significantly. The most important examples can be grouped into three categories: (i) lepton flavor universality tests, (ii) CP asymmetries, (iii) observables that vanish in the absence of right-handed currents. In most cases, the current experimental uncertainties of these observables are considerably larger than the theory uncertainties both in the SM as well as in viable new physics scenarios and the impact of the theory uncertainties in the global fit is moderate. However, with the expected improvement in experimental sensitivity, the theoretical uncertainties will become more and more important and their new physics dependence needs to be taken into account.
Among the lepton universality tests, ratios of branching ratios, like R K and R K * , are known with high precision in the SM, with uncertainties of around 1% [25,26]. In the presence of new physics, however, the uncertainties can be several percent. On the experimental side, the most precisely known quantity is R K , with an uncertainty of ∼ 4% [11], c.f. Eq. (2). After run 3 of the LHC, with ∼ 25 fb −1 of integrated luminosity collected by LHCb, one expects an experimental uncertainty of R K (R K * ) of ∼ 2.5% (2.8%) [27] assuming that systematic uncertainties can be controlled. The precision might reach ∼ 1% with 300 fb −1 . This clearly shows the need to consistently take into account the theory uncertainties including their new physics dependence. Other lepton universality tests, like the differences of angular observables [28] (denoted by Q i in [29,30]), have currently sizeable experimental uncertainties [30] and do not play a major role in global fits, yet. However, given the expected future experimental precision of a few percent [27] it becomes desirable to have a robust treatment of their theory uncertainties as well.
In Figs. 1 and 2 we illustrate the above points with a few examples. The plots in Fig. 1 show the theory predictions for R K and R K * (in the q 2 bin from 1.1 GeV 2 to 6 GeV 2 ) in the presence of new physics parameterized by various Wilson coefficients. As is well known, the Wilson coefficients with left-handed quark currents (C 9 and C 10 ) lead to a correlated effect in R K and R K * , while for right-handed quark currents (C 9 and C 10 ) one finds an anti-correlation [31]. For C 9 = −C 10 one has to an excellent approximation R K R K * . The various colored bands show the theoretical uncertainties at the 1σ and 2σ level. Circle, square, and diamond markers correspond to Wilson coefficient magnitudes of 0.5, 1.0, and 1.5. Colored markers correspond to positive, white markers to negative values. While the uncertainties are negligible close to the SM point, they become sizeable away from it. For comparison, we also show the current experimental results with 1σ uncertainties [11,32], as well as the expected uncertainties after run 3, assuming the same central value.
Similarly, the plot in Fig. 2 shows the theory predictions for D P 4 and D P 5 (in the q 2 bin from 1 GeV 2 to 6 GeV 2 ) in the presence of a few combinations of non standard Wilson coefficients. Also here we observe that the theory uncertainties can be sizable away from the SM point. As the current experimental uncertainties are still large [30], we show as comparison the expected experimental uncertainties with the full Belle II data-set which we expect to be around 5%, 1  With regards to CP violation, we note that results on CP asymmetries in B → K * μ + μ − are available from LHCb with 3fb −1 of run 1 data [33]. The most interesting asymmetries are A 7 , A 8 , and A 9 as they are not suppressed by small strong phases and therefore could in principle be O (1) in the presence of CP violating new physics [34] (Interesting CP asymmetries in B → K μ + μ − have been recently discussed in [35]). In the SM, they are strongly Cabibbo suppressed, [36]. The available experimental results are all compatible with zero with uncertainties of approximately 5% [33] both at low q 2 ∈ (1 GeV 2 , 6 GeV 2 ) and at high q 2 ∈ (15 GeV 2 , 19 GeV 2 ). Scaling with √ N , we expect sensitivities with the run 2 data set of approximately 2%−3% and ultimate sensitivities of below 1% with 300fb −1 .
In Fig. 3 we show the theory predictions for the B → K * μ + μ − CP asymmetries A 7 and A 8 (in the q 2 bin from 1.1 GeV 2 to 6 GeV 2 ) in the presence of imaginary parts of Wilson coefficients. Similarly to the LFU observables discussed above, also here we observe non-negligible theory uncertainties away from the SM point. For comparison, we also show the current experimental results with 1σ uncertainties [33], as well as uncertainties of 1%, assuming the same central value.

The updated global fit
In comparison to our previous fit in [6], we improve the treatment of the theory uncertainties as described in the previous section and we include a series of new experimental results: • The update of the B 0 → K * 0 μ + μ − angular analysis with 2016 data from LHCb [37]. The P 5 anomaly persists in this recent update, with a slightly reduced significance compared to the run 1 results [33]. Included in our fit are the angular observables F L , P 1 , P 2 , P 3 , P 4 , P 5 , P 6 , and P 8 in all available q 2 bins below 6 GeV 2 and the one large q 2 bin above the narrow charmonium resonances. • The new B ± → K * ± μ + μ − angular analysis [38]. While the experimental uncertainties of the B ± → K * ± μ + μ − angular analysis are still sizeable, deviations from SM predictions are observed that are broadly showing the same pattern as in the B 0 → K * 0 μ + μ − angular analysis. • The latest results on B s → μ + μ − from CMS [18] and the very recent result from LHCb [13,14]. We combine these results with the ATLAS result [17], as described in Appendix A. Compared to the previous LHC combination [15], our combination has a slightly larger central value and a slightly reduced relative uncertainty. • The recent update of R K [11]. The new result has exactly the same central value but reduced uncertainty compared to the previous result [12], increasing the tension with the SM from 2.5σ to 3.1σ . • The latest results from LHCb and CMS on the effective B s → μ + μ − lifetime, τ eff = (2.07 ± 0.29 ± 0.03) × 10 −12 s [13,14] and τ eff = (1.70 +0.61 −0.44 ) × 10 −12 s [18] (see [16] for the previous LHCb result). Precision measurements of τ eff can lead to non-trivial constraints on new physics in the form of the scalar Wilson coefficients C ( ) S,P [39,40].
• The recent update of the B s → φμ + μ − branching ratio [41] that confirms the previously seen tension [42] with the SM prediction.
Our numerical code is based on the Python package flavio [43], which provides all the theory predictions including their uncertainties and correlations. We use the full set of b → s observables and measurements as implemented in the Python package smelli v2.3.1 [44,45], which builds upon flavio v2.3.0. We plan to implement our new error treatment (cf. Sect. 2) in future versions of flavio and smelli.

One parameter scenarios
We start by considering simple one parameter new physics scenarios, switching on one real new physics Wilson coefficient at a time. We consider several fits, including certain subsets of observables. In Table 1 we report the best fit values for the Wilson coefficients as well as the 1σ best-fit regions and the "pull" in σ , defined as the χ 2 between the best fit point and the χ 2 of the SM.
In the column "b → sμμ" in Table 1, we focus on the b → sμμ observables that include the differential branching ratios of B → K μ + μ − , B → K * μ + μ − , B s → φμ + μ − , and b → μ + μ − as well as all available CP averaged angular observables in these decays. Note that these observables are subject to potentially large hadronic uncertainties. While existing calculations indicate that long distance effects are well within the assumed uncertainties [46], it cannot be fully excluded that such effects are unexpectedly large. As the considered decay modes do neither involve electrons nor are sensitive to scalar operators, only results for vector and axial-vector muonic Wilson coefficients are shown. Consistent with previous findings, we observe that a negative C bsμμ 9 −0.75 or the left-handed muon combination C bsμμ 9 = −C bsμμ 10 −0.53, are strongly preferred by the fit. For those values of the Wilson coefficients the agreement between theory and data is improved by more than 3σ compared to the SM In the column "LFU, B s → μμ" in Table 1, we consider the neutral current LFU observables (R K ( * ) , D P 4,5 ) and BR(B s → μ + μ − ) only, including in particular the new Table 1 Best-fit values with corresponding 1σ ranges as well as pulls in sigma between the best-fit point and the SM point for scenarios with NP in a single real Wilson coefficient. Column "b → sμμ": fit including only the b → sμμ observables (branching ratios and angular observables). Column "LFU, B s → μμ": fit including only the neutral current LFU observables (R K ( * ) , D P 4,5 ) and = −C bsμμ 10 −0.35, which have a pull of 4.7σ and 4.6σ , respectively. These scenarios do not only address the anomalies in R K and R K * , but also the slightly reduced branching ratio of B s → μ + μ − . The coefficients C bsμμ 9 , C bsee 9 , and C bsee 10 can explain the R K and R K * data, but do not affect the B s → μ + μ − decay. Their pulls are therefore a bit lower, around 4σ . The scalar Wilson coefficients show a slight (∼ 2σ ) preference for negative values, that lead to a suppression of the B s → μμ branching ratio in accordance with the data. Note that we include the effect of the scalar Wilson coefficients only in the B s → μ + μ − decay. In the parameter space allowed by B s → μ + μ − , the scalar Wilson coefficients have negligible impact on all the other b → sμμ transitions.
Finally, in the the column "all rare B decays" in Table 1 we show the results of the global fit. Included are the b → sμμ observables, the LFU observables, and the B s → μ + μ − branching ratio. 2 The largest pulls of 5.6σ and 5.2σ are found for C −0.73, respectively. As expected, the pulls for the electronic Wilson coefficients are very similar to the values in the "LFU, 2 Note that in previous fits [6] we had also included F = 2 observables that are correlated to the B s → μ + μ − branching ratio and the various b → sμμ branching ratios, mainly through their dependence on common CKM input. Adding F = 2 observables in the fit further increases the pulls slightly. B s → μμ" column. We observe a small change in the preferred values for the scalar Wilson coefficients, which is due to the correlations of the theory uncertainties of BR(B s → μ + μ − ) and the b → sμμ observables.
To illustrate the impact of our improved treatment of theory uncertainties, we compare in Table 2  scenario, in which the pull from the b → sμμ observables is somewhat reduced once the new physics dependence of the theory errors is taken into account. We expect the effect to become much more pronounced with more precise data.

Two parameter scenarios
Next, we discuss scenarios where two Wilson coefficients are turned on simultaneously. In Fig. 4 we show the best fit regions in the C bsμμ 9 vs. C bsμμ 10 plane. The plot on the left focuses on the constraints from the LFU ratios R K and R K * . The R K constraint before the update [11] is shown by the dashed contours. As the measured R K > R K * the best fit range prefers a sizable positive C bsμμ 10 . The plot on the right shows the result of the global fit. The B s → μ + μ − branching ratio prefers a modest positive C bsμμ 10 , while the b → sμμ observables mainly prefer a negative C bsμμ 9 . Overall, the best fit point corresponds to (C    [11,13,41] In Fig. 5 we show the viable parameter space of a couple of other Wilson coefficient pairs, that were found to give good fits in the past. The plot on the left shows the C bsμμ 9 vs. C bsμμ 9 plane, while the plot on the right shows the C univ. ). The best fit points are given by (C ) (−0.32, −0.34) and correspond to pulls of 5.0σ and 5.4σ , respectively. The scenario on the left gives an excellent fit of R K and R K * , but the slightly reduced B s → μ + μ − branching ratio remains unexplained. The scenario on the right can resolve the tension in BR(B s → μ + μ − ), but leaves a tension between R K and R K * . Note that C univ. 9 could in principle be mimicked by a hadronic effect. A lepton flavor universal C univ. 9 of the pre-ferred size can also be generated through renormalization group running from semi-tauonic operators that are motivated by the R D ( * ) anomalies [21] or from four-quark operators [6].
As clearly seen in the plots of Figs. 4 and 5, the branching ratio of B s → μ + μ − plays an important role in constraining the Wilson coefficient C 10 . It is well known that B s → μ + μ − is also very sensitive to new physics in the scalar Wilson coefficients (see e.g. [40]). In Fig. 6 we show the constraints in the Wilson coefficient plane C (left) and C univ.   (0, 0). The dashed lines show the constraints before the recent update [13] SM prediction and the experimental world average is clearly reflected in the plot. With the recent BR(B s → μ + μ − ) update, the preferred region in the Wilson coefficient space moved slightly towards the SM point. We observe that the measurements of the effective B s → μ + μ − lifetime already have some impact on the allowed parameter space of the scalar Wilson coefficients. The region of parameter space that corresponds to a mass eigenstate rate asymmetry A = −1 is excluded at the 1σ level. Note that the latest LHCb result for BR(B s → μ + μ − ) assumes the SM value A = +1. Due to the lifetime dependence of the acceptance, the experimentally determined BR(B s → μ + μ − ) is larger by approximately 5% or 11% for A = 0 or −1, respectively [13]. A similar effect is observed in the ATLAS and CMS analyses [15]. We do not attempt to model this effect in our fit of the scalar Wilson coefficients. In the region that is currently slightly disfavored by the measured effective B s → μ + μ − lifetime, we expect a few percent shift of the best fit band.

Generic scenarios
We also consider more generic scenarios with more then two Wilson coefficients. In particular, we consider a four parameter scenario including the muon-specific semi-leptonic Wilson coefficients C , and C bsee 10 . In the four parameter scenario we perform two fits: (1) a fit including only the b → sμμ observables (branching ratios and CP averaged angular observables) and (2) the global fit of all rare B decay data, including the LFU observables and B s → μ + μ − . In both cases we identify the best fit point in Wilson coefficient space and approximate the likelihood function in its vicinity by a multivariate Gaussian. The Table 3 Best fit values, uncertainties, and correlation matrix of the four-parameter fit to the Wilson coefficients C   with large significance. The corresponding central value is close to the result found in the one-parameter fit to C bsμμ 9 discussed in Sect. 3.1. We find sizable correlations among the Wilson coefficients. One of the main contributors to the correlations is the new precise measurement of R K , as can be seen in the two parameter scenarios shown in Figs. 4 and 5. We find a positive correlation between C bsμμ 9 and C bsμμ 10 that increases when R K is included, as expected from Fig. 4. In the four parameter +0.14 ± 0.23 and slightly negative correlation in [49]. We checked that excluding the high-q 2 bins from our fit (as done in [49]) improves the agreement with [49] to some extent, but differences remain.
We find similar results in the six parameter scenario. The parameters of the multivariate Gaussian that approximates the likelihood function in the vicinity of the best fit point of the global fit is reported in Table 5. The results for the muon specific Wilson coefficients are very similar to the four parameter fit discussed above. New physics effects in the electron-specific Wilson coefficients C bsee 9 and C bsee 10 are compatible with zero. The uncertainties of C bsee 9 and C bsee 10 are large and highly correlated. Shown separately are the constraints from LFU observables, CP con-serving b → sμμ observables, the B → K * μ + μ − CP asymmetries, and the global fit. The dashed lines show the constraints before the recent updates [11,41]

Complex Wilson coefficients
In the presence of new physics, the contributions to the flavor changing Wilson coefficients can generically be CP violating. While the observables that show tensions with SM predictions are CP conserving, it is interesting to investigate the impact that imaginary parts of Wilson coefficients have on the fit, and to which extent imaginary parts are constrained by existing data (see also [50] for a recent study that considers complex Wilson coefficients).
In Fig. 7 we show constraints in the planes of complex  and Im . The superscripts on the observables indicate the q 2 range in GeV 2 observables, the B → K * μ + μ − CP asymmetries from [33], and the global fit.
In the case of C bsμμ 9 , the experimental data does not lead to relevant constraints on the imaginary part of the Wilson coefficient, yet. In fact the strongest constraint on Im(C bsμμ 9 ) arises due to the fact that a sizeable imaginary part universally enhances the b → sμμ rates. We observe that the other scenarios Im(C bsμμ 9 ), Im(C bsμμ 10 ), and Im(C ) =Im(C bsμμ 10 ) are already being constrained by the experimental data on the CP asymmetries. Still, the current measurements do leave room for imaginary parts that are at least as large as the corresponding real parts. All imaginary parts are compatible with zero at the 2σ level. The best fit points of the real part of the Wilson coefficients are very close to the values that we obtain setting the imaginary parts to zero.

Predictions for LFU observables and CP asymmetries
As discussed in the previous section, several new physics Wilson coefficients (or combinations of Wilson coefficients) can significantly improve the agreement between data and theory predictions. The various best fit points show comparable pulls, and it is therefore interesting to identify predictions that allow us to distinguish the new physics scenarios.
We consider six different two parameter new physics scenarios: (i) Re C . In each of these cases, we sample the likelihood of the Wilson coefficients and show in Table 6 the predictions for several observables.
The first set of rows shows the predictions for the LFU ratios R K , R K * , and R φ both at low q 2 and at high q 2 . Overall, the predictions are fairly similar in all the considered new physics scenarios. Given the precise measurement of R K at low q 2 that enters the global fits, all scenarios reproduce the measurement of 0.85 at the 1σ level. The predicted values for all other LFU ratios are similar in all scenarios (i) -(vi).
The central values are all expected between 0.8 and 0.9. This is in particularly true for R K * where the current experimental result is considerably lower.
The second set of rows shows predictions for LFU differences of B → K * μ + μ − angular observables: D P 5 , D P 4 , and D A FB . Here we find significant differences in the various scenarios. In particular, precise measurements of D P 5 will allow to narrow down new physics scenarios.
The last set of rows shows predictions for the B → K * μ + μ − CP asymmetries A 7 and A 8 . The CP asymmetries remain close to zero (i.e. SM-like) in the scenarios (iv)-(vi) as they do not contain any new sources of CP violation. In scenarios (i)-(iii), A 7 and A 8 can be non-zero. Interestingly, an imaginary part of C bsμμ 9 leads to an effect in A 8 , while an imaginary part of C bsμμ 10 leads to an effect in A 7 . The predicted ranges for A 7 and A 8 can already be probed with run 2 data.
In Figs. 8, 9, and 10, we show the most distinctive cases in graphical form. The plots of Fig. 8 contain the predictions for the LFU ratios in scenarios (i), (iii), and (iv). The new physics predictions are compared to the SM predictions (with uncertainties from [25]) and the current experimental results [11,32]. Similarly, the plots of Fig. 9 show predictions and experimental results [30] for the LFU differences in scenarios (i), (ii), and (iii). The uncertainties of the SM predictions are illustrated with ±0.01. Finally, the plots of Fig. 10 show the CP Asymmetries in the scenarios with imag- Fig. 8 Predictions for the LFU ratios R K , R K * , and R φ in three new physics scenarios and the SM. For comparison the current measurements from LHCb [12,32] are shown as well Fig. 9 Predictions for the LFU differences D P 5 , D P 4 , and D AFB in three new physics scenarios and the SM. For comparison the current measurements from Belle [30] are shown as well Fig. 10 Predictions for the CP asymmetries A 7 and A 8 in three new physics scenarios and the SM. For comparison the current measurements from LHCb [33] are shown as well inary parts (i), (ii), and (iii). The tiny SM uncertainties are neglected and the experimental results are taken from [33]. The plots clearly show the discrimination power of the different observables.

Conclusions
With the recent updates of R K and BR(B s → μ + μ − ) by LHCb, the case for new physics in rare B decays has been further strengthened. Our improved global fit shows very strong preference for the muon specific Wilson coefficients from the experimental results on the B → K * μ + μ − CP asymmetries.
Finally, we give new physics predictions for a large set of observables including LFU ratios, LFU differences of CP averaged B → K * μ + μ − observables, and B → K * μ + μ − CP asymmetries. Future more precise measurements of these observables will allow us to distinguish between different new physics scenarios. Note Added Another model independent interpretation of the new results can be found in [49]. First interpretations in new physics models have been presented in [51,52].

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: There is no associated data. We will implement our new error treatment in future versions of the publicly available software packages flavio and smelli. The numerical results and plots can then be reproduced using these tools.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .
Since the B 0 and B s have a similar mass the measurements of the B 0 → μ + μ − and B s → μ + μ − branching ratios are correlated and the experimental results are given by two-dimensional likelihoods. We combine them assuming the likelihoods of different experiments are uncorrelated. The individual likelihoods are shown as thin lines in Fig. 11 while our combination is shown as thick solid red line. We also determine a Gaussian approximation (shown as thick dashed red line) and compare the experimental results to the SM predictions.
The two-dimensional Gaussian approximation is given by with an error correlation of ρ = −0.27.
Comparing the SM predictions with the two dimensional experimental likelihood we get the following onedimensional pulls 3 : • if both branching ratios are SM-like, 2.3σ , 4 • if B s → μ + μ − is SM-like and B 0 → μ + μ − profiled over, 1.9σ , Given its prominent role in constraining new physics in b → sμμ transitions, it is of great interest to have confidence regions for the B s → μ + μ − branching ratio itself, fixing B 0 → μ + μ − either to its SM central value or profiling over Fig. 11 Likelihood contours in the plane of BR(B 0 → μ + μ − ) and BR(B s → μ + μ − ) from the individual ATLAS, CMS, and LHCb measurements (thin contours), our combination (thick solid contours), and the Gaussian approximation (thick dashed contours). Also shown are the SM predictions and their 1σ correlated uncertainties it. Using our two-dimensional likelihood, we find For B 0 → μ + μ − we get analogously B Appendix: Details on theory uncertainties

B.1 Parameterization of Non-Factorizable Effects
We parameterize the non-factorizable effects in the decay amplitudes of semileptonic rare B decays following [22,54]. For B → K decays, the Wilson coefficient C eff 9 (q 2 ) is modified in the following way where low q 2 and high q 2 refers to di-lepton invariant masses below and above the narrow charmonium resonances, respectively. The central values of the complex parameters a K , b K , and c K are set to zero and the 1σ uncertainties enclose the effects considered in [55][56][57] Re(a K ) = 0.0 ± 0.08, Re(b K ) = 0.0 ± 0.03, Im(a K ) = 0.0 ± 0.08, Im(b K ) = 0.0 ± 0.03, We use the same ranges for B + → K + and B 0 → K 0 decays and assume that the corresponding coefficients are correlated by +99% due to iso-spin symmetry. For B → K * and B s → φ decays we use the following parameterization where the replacement of C eff 7 is performed only in the λ = 0, − helicity amplitudes, and the replacement of C 7 only in the λ = + amplitude. Furthermore, we have in all the helicity amplitudes. We use the following values for the hadronic parameters Re(a + ) = 0.0 ± 0.004, Re(b + ) = 0.0 ± 0.005, Im(a + ) = 0.0 ± 0.004, Im(b + ) = 0.0 ± 0.005, Re(a − ) = 0.0 ± 0.015, Re(b − ) = 0.0 ± 0.01, Im(a − ) = 0.0 ± 0.015, Im(b − ) = 0.0 ± 0.01, Re(a 0 ) = 0.0 ± 0.12, Re(b 0 ) = 0.0 ± 0.05, Im(a 0 ) = 0.0 ± 0.12, Im(b 0 ) = 0.0 ± 0.05, The same ranges of the parameters are considered for B 0 → K * 0 , B + → K * + , and B s → φ decays. A +99% correlation is assumed between the B 0 → K * 0 and B + → K * + coefficients (due to iso-spin), and a +90% correlation between the coefficients for the B s → φ decay and the B → K * decays (due to SU (3) symmetry). The above treatment of the non-factorizable effects is implemented in flavio since version 1.0.

B.2 Implementation of the New Physics Dependence
The decay amplitudes of rare semileptonic b hadron decays are linear functions of the Wilson coefficients. Thus, in the presence of new physics, the angular coefficients in the differential decay rates are second order polynomials in the new physics Wilson coefficients. Any observable O k in rare semileptonic decays that we consider can therefore be written as a function of second order polynomials p i O k = f k ( p 1 , p 2 , . . . , p n ).
For example, binned branching ratios are given directly in terms of a single second order polynomial, f k ( p 1 ) = p 1 . The CP averaged angular observables S i , the CP asymmetries A i , and the LFU ratios are ratios of two second order polynomials f k ( p 1 , p 2 ) = p 1 / p 2 . The angular observable P 5 has the form f k ( p 1 , p 2 ) = p 1 / √ p 2 (1 − p 2 ), and so on. The polynomials can be written in terms of a vector product factors of are introduced to track the order in the Wilson coefficients and they will be set to = 1 in the end. The vectors p i = (a i , b T i , c T i ) T in (32) are independent of the new physics. They depend on the considered observable and are given in terms of known input parameters. For any set of observables we can determine the covariance matrix p for the corresponding set of vectors p i . If N polynomials and M Wilson coefficients are involved, p is a N (1 + M + M 2 ) × N (1 + M + M 2 ) matrix. 6 We infer p by varying the input parameters within uncertainties, assuming Gaussian distributions.
For branching ratios, the functions f k are the identity, the observables depend linearly on the p i , and the number of polynomials, N , is equal to the number of observables. In this case, the N × N theory covariance matrix th that enters the χ 2 function (11) can simply be written as (see e.g. [58]) This th contains the exact dependence on the new physics Wilson coefficients. If the new physics Wilson coefficients are set to zero, it reduces to the theory covariance matrix in the SM. Expressing th as above has the big advantage that the new physics dependence is given analytically and the time consuming numerical determination of p has to be performed only once. In cases where the functions f k are non-trivial, the th with the exact new physics dependence can not be found in a simple analytical way from p . However, one can still find an analytic approximation in the limit of small new physics. If the new physics Wilson coefficients are small compared to the SM values, we can expand the functions f k in and write them as polynomials f k ( p 1 , p 2 , . . . , p n ) = p k = p k · V + O( 3 ) The coefficients of these polynomials are given by f k (a 1 , a 2 , . . . , a n ), b k = g i k b i ,  6 In practice, the size of the covariance matrix p can be slightly reduced by using the fact that only M(M + 1)/2 out of the M 2 entries Footnote 6 continued in D = vec( C ⊗ C) are independent and that usually some of the components of the p i are exactly zero.
As above, it is straight forward to determine the covariance matrix p of the vectors p k . Since all approximated observables are linear in p k , we find analogously to (33) where N is the number of polynomials p k , which equals the number of observables. The approximation can be improved systematically by expanding the functions f k in (34) to higher order in . In that case, the vector V has to be extended to include higher powers of the Wilson coefficients. As the observables are still linear in the coefficients p k , (37) continues to hold at any fixed order of the expansion. Note, however, that the size of the covariance matrix p grows rapidly with the order of the expansion.