Continuing search for new physics in $b \to s \mu \mu$ decays: two operators at a time

The anomalies in the measurements of observables involving $b \to s \mu\mu$ decays, namely $R_K$, $R_{K^*}$, $P_5^{\prime}$, and $B_s^\phi$, may be addressed by adding lepton-universality-violating new physics contributions to the effective operators ${\cal O}_9, {\cal O}_{10}, {\cal O}^\prime_9, {\cal O}^\prime_{10}$. We analyze all the scenarios where the new physics contributes to a pair of these operators at a time. We perform a global fit to all relevant data in the $b \to s$ sector to estimate the corresponding new Wilson coefficients, $C_9^{\rm NP}, C_{10}^{\rm NP}, C_9^\prime, C_{10}^\prime$. In the light of the new data on $R_K$ and $R_{K^*}$ presented in Moriond 2019, we find that the scenarios with new physics contributions to the ($C_9^{\rm NP}$, $C_9^\prime$) or ($C_9^{\rm NP}$, $C_{10}^\prime$) pair remain the most favored ones. On the other hand, though the competing scenario ($C_9^{\rm NP}$, $C_{10}^{\rm NP}$) remains attractive, its advantage above the SM reduces significantly due to the tension that emerges between the $R_K$ and $R_{K^*}$ measurements with the new data. The movement of the $R_K$ measurement towards unity would also result in the re-emergence of the one-parameter scenario $C_9^{\rm NP} = -C_9^\prime$.


Introduction
The Standard Model (SM) of particle physics cannot be the ultimate theory of fundamental interactions of nature. The necessity for new physics (NP) beyond SM is indicated from multiple directions, such as the neutrino masses, baryon asymmetry in the universe, dark matter, etc. Flavor physics is one of the most incisive probe of such NP, since new particles with masses beyond the reach of current experiments can contribute to low-energy processes through quantum corrections. These NP effects may be measurable at dedicated flavor experiments like LHCb [1] and Belle-II [2], as well as at multipurpose experiments like ATLAS [3] and CMS [4]. Deviations from the SM predictions, observed in the measurements of processes sensitive to such effects, can provide indirect indications of heavy particles or new interactions. These NP effects may be quantified in a model-agnostic way, using the language of effective field theory, by introducing additional operators to the SM effective Hamiltonian governing the relevant processes.
Over the last few years, the rare decays of B mesons, in particular the decays induced by the quark level transition b → s + − ( = e, µ) have already provided some such tantalizing hints of NP.
• The R K anomaly: The LHCb collaboration, in 2014, reported the measurement of the ratio R K ≡ Γ(B + → K + µ + µ − )/Γ(B + → K + e + e − ) in the "low q 2 " range (1.0 GeV 2 ≤ q 2 ≤ 6.0 GeV 2 ), where q 2 is the invariant mass-squared of the dilepton [5]. This measurement deviates from the SM value of 1 [6,7] by 2.6 σ, and is an indication of lepton flavor universality (LFU) violation. This measurement was recently updated in Moriond 2019, including the Run-II data and an update of the Run-I analysis. The measurement of R K from the Run-II data is reported to be R K (Run-II)= 0.928 +0.089+0.020 −0.076−0.017 , while the combined measurement from both the runs is R K (new)= 0.846 +0.060+0.016 −0.054−0.014 [8]. Clearly the central value of R K is moving towards unity, however the discrepancy with SM has remained ≈ 2.5σ.
• The R K * anomaly: The LFU violation in b → s µ + µ − sector was further corroborated by the measurement of the related quantity R K * ≡ Γ(B 0 → K * 0 µ + µ − )/Γ(B 0 → K * 0 e + e − ) in April 2017. The ratio R K * was measured in the low-q 2 (0.045 GeV 2 ≤ q 2 ≤ 1.1 GeV 2 ), as well as in the central-q 2 (1.1 GeV 2 ≤ q 2 ≤ 6.0 GeV 2 ) bin [9]. These measurements differ from the SM predictions of R K * 1 [6,7] by ≈ 2.4σ each. The Belle collaboration has presented their first measurements of R K * in B 0 decays, and the world's first measurement of R K * in B + decays, in Moriond 2019 [10]. These measurements, in multiple q 2 bins, have comparatively large uncertainties, and hence the anomaly in R K * still stands at ≈ 2.4σ level.
The SM predictions of R K and R K * are theoretically clean [6,7], therefore the deviations of these measurements from the SM are clear indications of NP. On the other hand, the calculations of P 5 and B φ s involve form factor uncertainties and undetermined power corrections [22][23][24][25], so by themselves these two anomalies cannot be considered as unambiguous signals of NP. However, since all these four observables are in the same (b → s + − ) sector, simultaneous anomalies observed in them should be taken seriously and addressed within the same framework. While the R K and R K * anomalies could be due to NP in b → sµ + µ − and/or b → se + e − decays [26][27][28], the discrepancies in P 5 and B φ s can be attributed to the presence of new physics only in b → s µ + µ − . Hence it would be natural to account for all of these anomalies by assuming new physics only in the b → sµ + µ − sector, which naturally breaks the LFU. We follow this assumption throughout this work.
We analyze the above four anomalies within the framework of effective field theory, with the aim of gauging the effects of new operators with different Lorentz structures that may contribute to b → sµµ processes. While the possible Lorentz structures are vector (V), axial vector (A), scalar (S), pseudo-scalar (P), and tensor (T), the last three are heavily constrained from the measurements of B s → µµ and b → sγ [29][30][31]. Hence in our analysis, we consider NP in the form of V and A operators only. Among possible operators, O 9 = (sγ µ P L b) (μγ µ µ) and O 10 = (sγ µ P L b) (μγ µ γ 5 µ) already exist in the SM effective Hamiltonian, however their Wilson coefficients (WCs) may be modified due to NP. There are also two chirality-flipped operators, O 9 = (sγ µ P R b) (μγ µ µ) and O 10 = (sγ µ P R b) (μγ µ γ 5 µ), which do not exist in the SM but may be provided by NP. We represent the WCs of these operators by C 9 , C 10 , C 9 and C 10 , respectively. The NP contribution to C 9 and C 10 are denoted by C NP 9 and C NP 10 , respectively, i.e. C 9 = C SM 9 + C NP 9 and C 10 = C SM 10 + C NP 10 . After the advent of the R K * result in 2017, several analyses were performed with an aim of identifying the Lorentz structure of possible NP [27,[32][33][34][35][36][37][38][39]. Most of these analyses showed that these anomalies, except the low-q 2 bin R K * measurement, may be explained by using a combination of C NP 9 , C NP 10 , C 9 , and C 10 . The explanation of the R K * (low-q 2 ) anomaly would need the introduction of a tensor operator [31]. On the other hand, a tensor operator alone cannot help in resolving the other anomalies considered in this paper. The resolution of the R K (low-q 2 ) anomaly is therefore decoupled from that of the others, and we do not dwell on that in this paper.
The most parsimoneous solutions to the anomalies would be the "1D" scenarios, where only one new WC contributes, or the values of two new WCs are related, so that there is only one extra parameter. The scenarios with only-C NP 9 , C NP 9 = −C NP 10 , or C NP 9 = −C 9 fit the data much better than the SM [27], though the last one seems to be disfavored since it predicts R K ≈ 1 [38]. The above 1D scenarios can indeed be generated in several proposed new physics models that contribute to b → s µ + µ − at the tree level. For example, Z models with gauge couplings to leptons can generate the only-C NP 9 scenario [40]. Some leptoquark models [41][42][43][44][45][46], and Z models with loop-induced couplings or with heavy vector-like fermions [47][48][49], can give rise to C NP 9 = −C NP 10 scenarios. In Z models with vector-like fermions and L µ − L τ symmetry, the C NP 9 = −C 9 scenario may be generated [50]. The "2D" scenarios, where NP contributes to two of the WCs, would be expected to give much better fits to the data than the SM or the 1D fits. The scenarios contributing to the pairs (C NP 9 , C NP 10 ), (C NP 9 , C 9 ) and (C NP 9 , C 10 ) have been shown to be able to account for all the above anomalies, except the low-q 2 bin R K * measurement, to a reasonable extent [27]. Out of these scenarios, the (C NP 9 , C 9 ) may be generated in Z models with couplings to leptons through the L µ − L τ portal [50]. The relative importance of these different 2D scenarios needs to be freshly analyzed in the light of the updated R K and R K * results.
In this paper, we analyze all the 2D scenarios, i.e. where NP contributes to two WCs at a time in an uncorrelated manner, with the inclusion of the 2019 Moriond update of the R K and R K * data. We perform a global fit to the anomalies as well as to the related data on observables that involve b → sµµ transitions and would be affected by the same WCs. Since all the observables we consider are CP-conserving, we restrict the WCs to be real. We also consider the fate of the 1D scenarios, which naturally emerge as subsets of the relevant 2D scenarios. We focus on pointing out any changes in the fits to the different scenarios due to the 2019 update. We also interpret these changes in terms of analytic approximations to R K and R K * in various scenarios.
The plan of the paper is as follows. In Sec. 2, we discuss the methodology adopted in our analyses. In Sec. 3, we provide the results of our fits and discuss various 2D scenarios and their 1D sub-scenarios. Finally, we summarize and conclude in Sec. 4, with a comparison among different scenarios.

Methodology
We represent the effective Hamiltonian for the decay b → sµµ in the presence of new physics V and A operators by where the SM effective Hamiltonian is Here G F is the Fermi constant and V ij are the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements. The Wilson coefficients C i of the four-fermi operators O i encode the short-distance contributions to the Hamiltonian, where the scale-dependence is implicit, i.e. 6,8) contribute to these processes through the modifications C 7 (µ) → C eff 7 (µ, q 2 ) and C 9 (µ) → C eff 9 (µ, q 2 ), where q 2 is the invariant mass-squared of the final state muon pair. The NP effective Hamiltonian is The NP effects are thus encoded in the Wilson coefficients C NP 9 , C NP 10 , C 9 and C 10 . While NP can in principle contribute to all the above four WCs, we focus on those scenarios where only two of these coefficients are nonzero. While this restriction is somewhat arbitrary at this stage, it is possible that symmetries of the NP at high scales can naturally make some of these coefficients vanish. The scenarios we consider may provide clearer insights on the role of NP Lorentz structures, due to the smaller number of parameters involved. We consider all six possible pairs of these coefficients, viz. (C NP 9 , C NP 10 ), (C NP 9 , C 9 ), (C NP 9 , C 10 ), (C NP 10 , C 9 ), (C NP 10 , C 10 ) and (C 9 , C 10 ). This analysis is also naturally applicable to the scenarios where only one of these coefficients is nonzero, or the two are linearly related, as considered in [51][52][53].
Note that all these observables are CP-conserving, as a result we do not expect to be sensitive to the complex nature of the new WCs. We therefore take C NP 9 , C NP 10 , C 9 and C 10 to be real for the sake of this article. We perform a two-dimensional (2D) χ 2 fit using the CERN minimization code MINUIT [62]. The χ 2 function is defined as Here O th (C i , C j ) are the theoretical predictions of the N=116 (122) observables before (after) the Moriond 2019 update used in the fit, while O exp are the experimental measurements. The N × N total covariance matrix C is obtained by adding the individual theoretical and experimental covariance matrices. The values of O th (C i , C j ) and the theoretical covariance matrix are calculated using flavio [63]. The correlations among O exp are included for the angular observables in B → K ( * ) µ + µ − [14] and B s → φµ + µ − [20]. For the other observables, we add the statistical and systematic errors in quadrature. Wherever the errors are asymmetric, we use the conservative approach of using the larger error on both sides of the central value.
We denote the value of χ 2 in the SM by χ 2 SM , and the best-fit value in the presence of NP by χ 2 bf . Clearly the addition of two degrees of freedom provided by the two new WCs decreases the χ 2 , and hence χ 2 SM > χ 2 bf . We define ∆χ 2 ≡ χ 2 SM − χ 2 bf for each pair of WCs, which would enable us to quantify the extent to which a particular combination of WCs is able to provide a better fit to the data. For convenience of notation, we denote the value of ∆χ 2 before (after) the 2019 update as ∆χ 2 old (∆χ 2 new ).

Results and discussions
We present the results of our 2D fits in the form of contour plots in the parameter space of the two relevant WCs, as shown in Fig. 1. The six plots correspond to the six scenarios with nonzero NP contributions to (C NP 9 , C NP 10 ), (C NP 9 , C 9 ), (C NP 9 , C 10 ), (C NP 10 , C 9 ), (C NP 10 , C 10 ) and (C 9 , C 10 ), respectively. In all plots, SM corresponds to the point (0, 0).
In the figure, we show the 1σ regions allowed from the measurements of (i) the ratio R K * (central bin: 1.0 GeV 2 < q 2 < 6.0 GeV 2 ), (ii) the average of the angular observable P 5 (4.0 GeV 2 < q 2 < 6.0 GeV 2 ) from the ATLAS and LHCb experiments [16], and (iii) the branching ratio B(B s → φ µ + µ − ), with bands of blue, pink, and green color, respectively. The  1σ allowed region of R K from the 2014 data [5] and the updated 2019 data [8] are shown by light and dark yellow bands, respectively. The overlaps (or lack of them) of these bands contain information about the consistency (or tension) among different anomalies. Note that none of these scenarios is able to account for the measured value of R K * in the low-q 2 bin within 2σ. So the band corresponding to this measurement is not shown in the plots, though  Figure 1: The 1σ allowed bands for R K (1.0 GeV 2 < q 2 < 6.0 GeV 2 ) before and after 2019 update, R K * in the central bin (1.0 GeV 2 < q 2 < 6.0 GeV 2 ), P 5 (4.0 GeV 2 < q 2 < 6.0 GeV 2 ) from ATLAS and LHCb, and B φ s ≡ B(B s → φµ + µ − ) in the range (1.0 GeV 2 < q 2 < 6.0 GeV 2 ), for the six 2D scenarios. The 1σ and 2σ allowed regions from the global fit using data before (after) the 2019 R K update are shown by dashed (solid) contours. Specific 1D sub-scenarios that give a good fit to the data are also shown. it contributes to the global fit. Also, the CMS results on P 5 [59] are not shown in the bands since they correspond to a different q 2 -range. The new R K * result from Belle [10] are also not shown, since they currently have large uncertainties. These results are, however, included in the global fit.
Superimposed on the above bands are the 1σ and 2σ contours, shown in brown and red, respectively, corresponding to the global fit to all 116 (122) observables, before (after) the Moriond 2019 update. The contours corresponding to the data before (after) the update have dashed (solid) boundaries. A comparison of these two sets of contours gives us an indication of how the preferred parameter space in the particular NP scenario has changed due to the 2019 update. The superposition of these contours on the 1σ bands of key individual measurements above allows us to check whether the best-fit region is indeed able to account for all the anomalies.
Some of the plots also indicate the lines corresponding to selected scenarios with linear relations between the two WCs which give good fits to the data. While the viability of these 1D sub-scenarios may be judged qualitatively from the figures, Table 1 lists the best-fit values of parameters, along with the ∆χ 2 old , ∆χ 2 new , and 1σ allowed regions for them. Below we list some important observations that may be made for the six scenarios. Since the measurements of R K and R K * are theoretically clean, and are expected to dominate the fits, we also try to understand the impact of new R K and R K * measurements by using analytic approximations for R K and R K * (central-q 2 ) in the presence of the corresponding NP. Henceforth in this section, we shall refer to R K * (central-q 2 ) simply as R K * for the sake of brevity.

The (C NP
9 , C NP 10 ) scenario This scenario improves the global fit significantly as compared to the SM, however ∆χ 2 new ≈ 42 has decreased substantially from its older value of ∆χ 2 old ≈ 47. This is partly an effect of the new R K measurement having moved closer to the SM prediction. The new measurements have also increased the tension of the global best fit with all the four individual anomalies marginally. This scenario still stands as one of the favored ones to account for these anomalies. The 1D sub-scenarios C NP 9 = −C NP 10 and C NP 10 = 0 also continue to improve the global fit, however the extent of improvement has reduced to ∆χ 2 new ≈ 37 with the new data, compared to ∆χ 2 old ≈ 41 from earlier. The relatively sharp decrease (compared to the other scenarios) in the value of ∆χ 2 after the Moriond 2019 update may be understood from the approximate functional forms [64] R K = R K * ≈ 1 + 0.24 (C NP 9 − C NP 10 ) . It can be seen that the values of R K and R K * are forced to be approximately equal in this scenario. While this was indeed the case before the update, after the update one has R K ≈ 0.85 and R K * ≈ 0.69. Thus, a tension has emerged in the measurements of these two quantities, thereby decreasing the overall goodness of fit.

3.2
The (C NP 9 , C 9 ) scenario This scenario already provided a slightly better fit to the data than the (C NP 9 , C NP 10 ) scenario, even before the 2019 update. With the update, ∆χ 2 old ≈ 50 for this scenario has stayed almost the same at ∆χ 2 new ≈ 49, indicating that it is still able to explain most of the data much better than the SM. Indeed, the fit is still consistent with R K * and P 5 , while its agreement with R K has improved with the new data. The 1D sub-scenario C NP 9 = −C 9 also has continued to provide a good fit to the data (∆χ 2 new ≈ 40), however earlier it was considered to be disfavored as it predicted R K ≈ 1 [38]. The updated data, however, has moved R K closer to unity. If this trend continues, this scenario could re-emerge as a favored NP solution.
In the (C NP 9 , C 9 ) scenario, the choices for C NP 9 and C 9 can allow R K and R K * to vary independently: No significant tension is therefore created because of the updated value of R K . The increase in the central value of R K after the update has only shifted the best fit point in the (C NP 9 , C 9 ) plane to higher values of C NP 9 and C 9 . More importantly, the increase in the R K measurement has directly decreased the value of the combination C NP 9 + C 9 , making the 1D sub-scenario C NP 9 = −C 9 more viable.

The (C NP
9 , C 10 ) scenario This scenario was the one with the largest ∆χ 2 old ≈ 53 among all the 2D global fits before the update, and stays so (∆χ 2 new ≈ 53) even with the update. It can accommodate R K and R K * anomalies within 1σ, and is quite close to the 1σ allowed regions for P 5 and B φ s . Note that the possible 1D sub-scenarios, C NP 9 = 0, C 10 = 0 or C NP 9 = −C 10 , do not improve the SM fit significantly, while C NP 9 = C 10 improves it by ∆χ 2 new ≈ 27. As far as the dependence of R K and R K * on the NP parameters is concerned, this scenario is similar to the previous one: While both these scenarios perform equally well in accounting for R K , R K * , and P 5 , the (C NP 9 , C 10 ) scenario can accommodate B φ s values closer to its measurement, and hence has a slightly better ∆χ 2 than (C NP 9 , C 9 ). The updated R K measurement shifts the best fit point to higher C NP 9 and lower C 10 .

The (C NP
10 , C 9 ) scenario This scenario offers a moderate improvement over the SM, with ∆χ 2 new ≈ 32. The best fit for this scenario continues to be able to account for the R K and R K * anomalies to within 1σ, however it cannot explain P 5 even within 2σ. The 1D sub-scenarios C NP 10 = C 9 and C 9 = 0 offer some improvement (∆χ 2 new ≈ 22) over the SM, however C NP 10 = −C 9 can only allow ∆χ 2 new ≈ 9. With only-C NP 10 , one can get ∆χ 2 new ≈ 27, which is only a moderate improvement.
The approximate functional forms of R K and R K * in this scenario are Since C 9 contributes to R K and R K * with opposite signs, in order to have both R K and R K * values less than unity, one would need a large value of C NP 10 . However, such a large value of C NP 10 is disfavoured by B s → µ + µ − measurement, which is close to its SM prediction. As a result, the improvement above SM is not significant in this scenario.

The (C NP
10 , C 10 ) scenario This scenario offers a moderate improvement over the SM, with ∆χ 2 new ≈ 29. The best fit for this scenario continues to be able to account for the R K and R K * anomalies to within 1σ, however it cannot explain P 5 even within 2σ. The 1D sub-scenarios, C NP 10 = −C 10 and C 10 = 0 offer some improvement (∆χ 2 new ≈ 16 and ∆χ 2 new ≈ 27, respectively) over the SM, however C NP 10 = C 10 can only allow ∆χ 2 new ≈ 8. The reason for only a moderate improvement in the goodness of fit over the SM is similar to the one in the previous scenario. Here, Thus C 10 contributes to R K and R K * with opposite signs, forcing C NP 10 to have unreasonably large values.
3.6 The (C 9 , C 10 ) scenario This scenario is not able to offer any significant improvement over the SM: both ∆χ 2 old and ∆χ 2 new are less than 1. As can be seen from the figure, the pairs of measurements (R K , P 5 ) and (R K * , B φ s ) pull the best fit point in almost opposite directions, thus keeping it close to the SM, without offering any solution to the anomalies. These opposite pulls are mainly the result of R K and R K * measurements. We have R K ≈ 1 + 0.24 (C 9 − C 10 ) , R K * ≈ 1 + 0.17 (−C 9 + C 10 ) . (3.6) In the presence of only these two new WCs, the values of R K and R K * are forced in opposite directions from unity. As long as the measured values of R K and R K * are both less than unity, the allowed values of C 9 and C 10 will stay small and cannot contribute to resolving both the anomalies simultaneously. The global fit will therefore stay poor.

Summary and conclusions
In this paper, we have explored whether pairs of new vector or axial vector effective operators would allow us to explain the anomalies observed in b → s decays, namely R K , R K * , P 5 , and B φ s . We have analyzed all the six pairwise combinations of the NP Wilson coefficients C NP 9 , C NP 10 , C 9 , C 10 that may contribute to the resolutions of these anomalies. We have performed global fits to data available before and after the Moriond 2019 update of R K and R K * , in order to obtain the favored values of the relevant WCs in these six scenarios. Our 2D global fits lead to the following observations: • The two scenarios (C NP 9 , C 9 ) and (C NP 9 , C 10 ) continue to offer significantly better fits to the data as compared to the SM (∆χ 2 new > 49), even with the 2019 update to the data. Both of these best fits can account for R K , R K * anomalies within 1σ, and P 5 , B φ s anomalies within 2σ.
• The scenario (C NP 9 , C NP 10 ), which used to give a significantly better fit (∆χ 2 old ≈ 47) than the SM before the 2019 update, cannot offer as good an improvement (∆χ 2 new ≈ 42) over the SM after the update. Indeed it is the only 2D scenario whose ∆χ 2 has undergone such a sharp decrease after the update, compared to the other ones. The scenario is still viable, though the tensions with individual experiments have increased with the update. The root cause of this may be traced to the approximately identical functional dependence of R K and R K * to the two WCs, C NP 9 and C NP 10 , in this scenario.
• The scenarios (C NP 10 , C 9 ) and (C NP 10 , C 10 ) continue to offer only moderate improvements (∆χ 2 new ≈ 25 − 30) over the SM. The worst scenario for explaining the anomalies turns out to be (C 9 , C 10 ). The best fit for this scenario is very close to the SM, and does not help in the simultaneous explanation of the anomalies.
Many features of the above global fits, and the changes in these fits after the R K and R K * update, may be understood in terms of the effect of new WCs on R K and R K * using analytic approximations. Note that the anomaly in the low-q 2 bin of R K * cannot be explained by any of these 2D fits, as has been pointed out earlier.
These 2D fits also allow us to explore their 1D sub-scenarios where only one new WC is nonzero, or where the two new WCs are linearly related. Such scenarios may be interesting not only from the point of view of smaller number of parameters, but also because such relations may prevent unwelcome effective operators from getting generated. The following 1D sub-scenarios offer significant improvements above the SM: • The C NP 9 = −C 9 scenario can give ∆χ 2 new ≈ 40. While this was still the case before the update, it was not considered to be a favored scenario since it predicted R K ≈ 1, in conflict with the older data. The update has moved R K in the direction of unity, and has made this scenario more attractive.
In our analysis, we have taken the data-driven approach and considered the addition of only a single, or a couple of, NP operators. While these would appear to be the most economical solutions in the language of effective field theory, they may not be always so from the point of view of constructing a high scale theory. While reducing the high scale theory to a low scale effective theory, the desired new effective operator(s) may be necessarily accompanied by other additional effective operators with different Lorentz structures. Putting the coefficients of these effective operators to zero is a possible way out, however the stability of such a scenario needs to be guaranteed by a symmetry at the high scale, or the scenario would involve some fine tuning of parameters. Here we take the approach that having a good fit in a 2D scenario guarantees an equally good (if not better) fit in the space with more than two NP parameters. The favored scenarios that have emerged with the updated data could help in narrowing down possible NP models and guiding constructions of models beyond the current paradigm.