On the anomalies in the latest LHCb data

Depending on the assumptions on the power corrections to the exclusive b ->s l+ l- decays, the latest data of the LHCb collaboration - based on the 3 fb^-1 data set and on two different experimental analysis methods - still shows some tensions with the SM predictions. We present a detailed analysis of the theoretical inputs and various global fits to all the available b ->s l+ l- data. This constitutes the first global analysis of the new data of the LHCb collaboration based on the hypothesis that these tensions can be at least partially explained by new physics contributions. In our model-independent analysis we present one-, two-, four-, and also five-dimensional global fits in the space of Wilson coefficients to all available b ->s l+ l- data. We also compare the two different experimental LHCb analyses of the angular observables in B ->K* mu+ mu-. We explicitly analyse the dependence of our results on the assumptions about power corrections, but also on the errors present in the form factor calculations. Moreover, based on our new global fits we present predictions for ratios of observables which may show a sign of lepton non-universality. Their measurements would crosscheck the LHCb result on the ratio R_K = BR(B+ ->K+ mu+ mu-) / BR(B+ ->K+ e+ e-) in the low-q^2 region which deviates from the SM prediction by 2.6 sigma.


Introduction
The LHCb collaboration has recently presented the angular analysis of the B 0 → K * 0 µ + µ − decay with the 3 fb −1 data set. They use two analysis methods. The observables are determined using an unbinned maximum likelihood fit and by the principal angular moments [1]. In addition, a new analysis on the angular observables in B s → φµ + µ − has been presented [2].
These new analyses of the LHCb collaboration have been eagerly awaited in view of the previous LHCb analysis of the B 0 → K * 0 µ + µ − based on the 1 fb −1 data set [3]. The LHCb collaboration had announced a local discrepancy of 3.7σ from the Standard Model (SM) predictions in one bin for one of the angular observables [3]. There had been also more, yet smaller tensions with the SM predictions in other observables. This announcement was followed by a large number of theoretical analyses showing that, due to the large hadronic uncertainties in exclusive modes, it is not clear at all whether this anomaly is a first sign for new physics beyond the SM or a consequence of underestimated hadronic power corrections or just a statistical fluctuation [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].
In the recent analysis based on the 3 fb −1 data set the LHCb collaboration now announced a 3.4σ tension with predictions based on the SM within a global fit to the complete set of CPaveraged observables [1]. They point out that this tension could be explained by contributions from physics beyond the SM or by unexpectedly large hadronic effects that are underestimated in the SM predictions.
Regarding the latter option it is important to note that there is a significant difference in the theoretical accuracy of the inclusive and exclusive b → s + − decays in the low-q 2 region. In the inclusive case, there is a theoretical description of power corrections; they can be calculated or at least estimated within the theoretical approach (for reviews see [23][24][25]). 1 In contrast, in the exclusive case there is no theoretical description of power corrections existing within the theoretical framework of QCD factorisation and SCET which is the standard theoretical framework for these exclusive decay modes in the low-q 2 region. Thus, power corrections can only be guesstimated. This issue makes it rather difficult or even impossible to separate new physics effects from such potentially large hadronic power corrections within these exclusive angular observables. So these tensions might stay unexplained until Belle II will clarify the situation by measuring the corresponding inclusive b → s observables as was demonstrated in Ref. [16,31,32].
Thus, it is also obvious that the significance of the tension depends in principle on the precise guesstimate of the unknown power corrections within the SM prediction. Because the two sets of SM predictions -LHCb compares with in their first and in their latest analysis -use rather different guesstimates [13,33] the quoted standard deviations in both analyses cannot be compared directly. 2 This situation motivated a recent theory analysis in which the unknown power corrections were just fitted to the data [34] using an ansatz with 18 additional real parameters in the fit. 1 In the inclusive case one can show that if only the leading operator of the electroweak hamiltonian is considered, one is led to a local operator product expansion (OPE). In this case, the leading hadronic power corrections in the decayB → Xs + − scale with 1/m 2 b and 1/m 3 b only and have already been analysed [26]. A systematic and careful analysis of hadronic power corrections including all relevant operators has been performed in the case of the decaȳ B → Xsγ [27]. Such linear power corrections can be analysed within soft-collinear effective theory (SCET). This analysis goes beyond the local OPE. An additional uncertainty of ±5% has been identified. The analysis in the case ofB → Xs + − is work in progress. There is no reason to expect any large deviation from theB → Xsγ result. Nonfactorisable power corrections that scale with 1/m 2 c have first been considered in Ref. [28], but can be now included in the systematic analysis of hadronic power corrections and be calculated quite analogously to those in the decayB → Xsγ [27]. Moreover, in the KS approach [29,30] one absorbs factorisable long-distance charm rescattering effects (in which theB → Xscc transition can be factorised into the product ofsb and cc colour-singlet currents) into the matrix element of the leading semileptonic operator O9. Following the inclusion of nonperturbative corrections scaling with 1/m 2 c , the KS approach avoids double-counting. 2 In this sense the significance of the tension with the SM has not been really reduced within the new 3 fb −1 measurement compared to the one based on the 1 fb −1 data set.
However, this fit to the data needs very large power corrections in the critical bins. If one compares the fitted theory predictions in this analysis with the leading contributions based on QCD factorisation then one finds 20% but also 50% or even larger power corrections relative to the leading contribution. The existence of such large hadronic corrections cannot be ruled out in principle. The authors of Ref. [34] come to the well-known conclusion that it is difficult to deduce the existence of new physics effects unambiguously from the measurements of exclusive b → s + − observables.
In this respect, another tension in the LHCb data of b → s + − gains importance. The ratio R K = BR(B + → K + µ + µ − )/ BR(B + → K + e + e − ) in the low-q 2 region had been measured by LHCb using the full 3 fb −1 of data, showing a 2.6σ deviation from the SM prediction [35]. This discrepancy has been addressed in many studies [32,. It is often claimed that the electromagnetic corrections might not have been fully taken into account in this measurement. Thus one might wonder whether this sign of lepton non-universality could be traced back to logarithmically enhanced QED corrections. These corrections were calculated in the inclusive case in Ref. [61]. However, LHCb uses the PHOTOS Monte Carlo to eliminate the impact of collinear photon emissions from the final state electrons. Therefore, such corrections do not seem to apply to the ratio R K , especially if one considers the agreement between the PHOTOS results and the analytical calculations in the inclusive case in Ref. [61]. Nevertheless, it would be an interesting check to correct for photon radiation using data-driven methods that do not rely on PHOTOS.
In contrast to the anomaly in the rare decay B → K * µ + µ − which is affected by power corrections, the ratio R K is theoretically rather clean and its tension with the SM cannot be explained by power corrections. But independent of this difference, both tensions might be healed by new physics in the semi-leptonic operator contribution C 9 . Therefore the measurement of other ratios which could show lepton non-universality would be a very important crosscheck of the present R K result but also of the anomalies in the angular observables in B → K * µ + µ − . Thus, we present predictions for such ratios based on the global fits in this paper.
Finally, we mention that there is an alternative approach to calculate the nonlocal charmloop effects. In Ref. [19] a local OPE near the light cone for the gluon emission from the c-quark loop at q 2 m 2 c is used to derive a nonlocal effective quark-antiquark gluon operator. Then the LCSR approach leads to an effective resummation of the soft-gluon part valid for q 2 m 2 c only. But a prediction of the charm loops effects up to the open charm threshold is achieved by a phenomenological model of the charmonium resonances via a dispersion relation. Because not all contributions were included in the dispersion relation this prediction for larger q 2 should not be regarded as the final result. An analogous complete calculation of the soft-gluon nonfactorisable contributions in the case of the decay of B 0 → K 0 + − [20] leads to a moderate effect. The corresponding predictions of the branching ratios for the various bins are given in Appendix A next to our predictions based on QCD factorisation. The two sets of predictions are compatible with each other at the 1σ level. Based on QCD factorisation there are two strategies to calculate the decay amplitudes: the so-called soft form factor (soft FF) approach (see [62,63]) and the full form factor (full FF) approach (see [64]). Both methods have been implemented in SuperIso v3.5 [65,66] which is used for all the calculations presented in this work. We discuss in detail the advantages and the disadvantages of these two approaches and critically analyse the guesstimate of power corrections which are used in the literature (Section 2). Among the input parameters and experimental data used in our study, we mainly discuss the theoretical error estimation and the form factor calculations via the light cone sum rule method and also define the statistical method used for the fitting analysis of the b → s data in Section 3. In contrast to many analyses in the literature, we consider one-and two-but also four-and five-dimensional global fits (in the space of Wilson coefficients) within our model-independent analysis of the present data. We analyse the dependence of our results on the theoretical approach used, as well as the assumptions about power corrections, and also on the errors present in the form factor calculations. We also investigate the effect of the S 5 and the R K measurements on the global fit results. Moreover, we compare the Wilson coefficient fits when using the two different experimental LHCb analysis for B → K * µ + µ − angular observables which this collaboration presented recently [1] (Section 4). Our results represent the first global fit analysis with the latest LHCb results for B → K * µ + µ − angular observables assuming that the tensions can be at least partially explained by new physics contributions. Finally, we use our global fit results to present predictions for other ratios of observables which may indicate lepton non-universality. As mentioned above, these measurements will be important crosschecks of the anomalies discussed in this paper (Section 5). We give our conclusions in Section 6 and additional analyses and more details on our theoretical predictions in the appendices.

Soft vs. full form factor approach
The theoretical description of the B → K ( * ) µ + µ − or B s → φµ + µ − decays in the low-q 2 region is based on the QCD-improved factorisation (QCDf) approach and its field-theoretical formulation of Soft-Collinear Effective Theory (SCET) [67,68]. The combined limit of a heavy b quark and an energetic meson M like K * , K, or φ leads to the schematic form of the decay amplitude which is valid to leading order in Λ QCD /m b and to all orders in α s . Thus the decay amplitude factorises into process-independent non-perturbative quantities like B → M soft form factors (ξ) and light-cone distribution amplitudes (LCDAs) of the heavy and light mesons (Φ) and perturbatively calculable quantities (C, T ). The key issue of this factorisation formula is that there are non-factorisable contributions (second term) beyond the ones described by form factors (first term). However, at order Λ QCD /m b the factorisation formula breaks down. This means that the power corrections cannot be calculated within QCDf in general. As already emphasised in the introduction, this is the main drawback of this approach which makes it difficult to identify new physics effects. There are power corrections to both terms in the factorisation formulae, which are called factorisable and non-factorisable power corrections respectively. The large energy limit allows for simplifications [69]. The seven independent full QCD B → K * form factors V, A 0,1,2 and T 1,2,3 , for example, are reduced to two universal soft form factors ξ ⊥ and ξ in this limit. These relations can also be written via a factorisation formula in a schematic way: where D and T F are perturbatively calculable functions. There are non-factorisable contributions (second term), but also power corrections (third term). Based on QCDf there are two strategies to calculate the decay amplitudes: the so-called soft form factor (soft FF) approach (see [62,63]) and the full form factor (full FF) approach (see [64]).
Using the former strategy, one takes advantage of the factorisation formula (2.1) and the universal soft form factors. The various meson spin amplitudes at leading order in Λ QCD /m b and α s turn out to be linear in the soft form factors ξ ⊥, and also in the short-distance Wilson coefficients. As was explicitly shown in Refs. [63,70], these simplifications allow to design a set of optimised observables, in which any soft form factor dependence (and its corresponding uncertainty) cancels out for all low dilepton mass squared q 2 at leading order in α s and Λ QCD /m b . An optimised set of independent observables was constructed in Refs. [33,71], in which almost all observables are free from hadronic uncertainties which are related to the form factors.
However, there are additional hadronic uncertainties in the factorisation formula (2.1), namely the unknown factorisable and non-factorisable power corrections. They are not cal-culable and can only be guesstimated through dimensional arguments within the QCDf/SCET approach [63,70].
As Eq. (2.2) indicates, the factorisable power corrections can be avoided by using the full QCD form factors. Thus, within this full FF approach one faces only unknown power corrections to the non-factorisable contributions; this means one has to guesstimate only power corrections to leading contributions from 4-quark operators O 1−6 and the chromomagnetic operator O 8 . 3 It is important to note that the operators O 7 , O 9 , and O 10 do not induce non-factorisable contributions. 4 As it was argued in Ref. [72], the correlations between QCD form factors in the large energy limit of the meson can be established also in this full FF approach using the equation of motions. We will discuss this issue further in Section 3.1.
In a global analysis both strategies should lead to similar results if all correlations are taken into account. If one focuses on a single observable, the advantage of the optimised observables on the theory side is unfortunately diminished by the large errors on the experimental side [1].
In this work we have used the full FF approach for calculating the B → K ( * )¯ and B s → φ¯ observables, as the main approach. For the uncertainty regarding the SM predictions only the non-factorisable power corrections are relevant in the full FF approach. We have guesstimated the effect due to these corrections by considering either 5, 10 or 20% and also 60% errors as described in Section 3.2.
As a crosscheck we use also the soft FF approach for an independent global fit with the same input data. In this case we make a guesstimate of the factorisable and the non-factorisable power corrections.
For the sake of completeness, we note that in the high-q 2 region a local operator product expansion can be used. One finds small power corrections below 5% [73,74]. Duality violation effects have been estimated within a model of resonances and have been determined to be again below 5%. But new resonances are found in this region (see e.g. Ref. [21]), so an update of this calculation would be desirable. The different theory framework in the high-q 2 region allows for a non-trivial crosscheck, especially of any new physics hypothesis in the low-q 2 region.

Input of model-independent analysis
We follow the methodology used in Ref. [32], but consider the following important improvement and updates within the experimental and theoretical inputs: • Unless otherwise specified, we have used the full form factor approach for the SM predictions of the B → K ( * )¯ and B s → φ¯ observables, assuming 10% error for the power corrections (see Sections 2 and 3.2). Crosschecks are also performed considering the soft form factor approach and different assumptions for the power corrections.
• We consider the LHCb measurements on the branching ratio as well as the angular observables of the B s → φ¯ decay. The theoretical treatment of the B s → φ¯ decay is similar to the B → K * ¯ decay, with the requisite replacements of the masses and hadronic parameters as well as the necessary changes resulting from the spectator effects. The required modifications can be found in [68]. Since the B s → φ¯ decay is not self-tagging, unlike 3 We note that there are also additional leading (calculable) non-factorisable corrections due to the rewriting of the soft form factors into full form factors according to Eq. (2.2). 4 In Ref. [18], the central value of factorisable power corrections are determined by fitting the soft form factors to the relevant full form factors using Eq. (2.2), but the corresponding uncertainties should still be guesstimated. From the principle point of view this procedure does not allow for any advantages in respect to the full FF form factor approach where one uses the full QCD form factors directly with the uncertainties of the light cone sum rule (LCSR) calculation. However, in view of the fact that there is only one independent LCSR calculation including the correlations [72], the authors of Ref. [18] prefer to use the full form factor calculation of Ref. [19] with much larger uncertainties. But no correlations are given in this calculation, so only this hybrid approach is possible.  the B → K * ¯ decay, the time integrated untagged average over theB 0 s and B 0 s decay distributions, including the effects ofB 0 s − B 0 s mixing, should be calculated [80]. Details about how to include this effect can be found in Refs. [80][81][82]. We have implemented this effect following Ref. [82].
• For the B → K * and B s → φ form factors we have used the combined fit results of LCSR calculations [72,83] and lattice computations [84,85] given in Ref. [72]. The combined fit results are applicable for the entire kinematic range of q 2 (see Section 3.1).
• For the B → K form factors we have used the combined fit results of LCSR calculations [86,87] and lattice computations [88] given in Ref. [89]. The combined fit results are applicable for the entire kinematic range of q 2 .
• For the branching ratios of B 0(+) → K 0(+) µ + µ − decays, the [1. • We have calculated the theoretical errors and correlations for the B → K ( * ) µ + µ − and B s → φµ + µ − decays with various assumptions for the power corrections as described in Section 3.2.
• For the experimental values of B → K * µ + µ − angular observables, we have considered both the LHCb results determined by the maximum likelihood fit method or the results for the finer bins determined by the method of moments [1] for comparison. Unless otherwise specified, for our analysis we have considered the experimental results obtained by the method of moments and used the recent results of the LHCb collaboration based on 3 fb −1 of data [1] for the S i observables. Alternatively, the experimental results of the optimised observables (P i ) can be used. The latter are determined by reparametrising the S i observables. The F L (1 − F L ) dependence of the optimised observables results in large and asymmetrical experimental errors in the bins where F L is close to zero or one, specially for the results which have been obtained using the method of moments (see Table 9 of Ref. [1] or Table 8 in Appendix A of the present paper).
• Some of the relevant updated input parameters of this work have been given in Table 1.
• The full list of observables that we have used as well as their SM predictions and experimental measurements can be found in Appendix A.

Form factors
As mentioned above, for the B → K * and B s → φ form factors we have used the combined fit results of LCSR and lattice calculations [72] which are applicable for the entire kinematic range of q 2 . One should consider the following aspects: • The form factor uncertainties are drastically reduced in the BSZ parametrisation compared to the ones in Refs. [19,90]. The main reason is very simple. In the latter reference the authors used another LCSR approach in which the role of the B meson and the light meson is exchanged [91][92][93]. Our knowledge about the B meson distribution amplitude is rather restricted compared to the one about the light-meson amplitude. This simple fact leads to much larger errors in this alternative LCSR approach somehow by construction. The BSZ parametrisation goes back to the original calculation given in Ref. [83]. An independent check of the calculation is missing. In view of this, we vary the errors given in LCSR calculation of Ref. [72] in our global fits in order to test the impact of the error estimation in this analysis.
• In Ref. [72] strong correlations between the QCD form factors are derived. In addition to the standard ones due to physical input parameters, the authors argue that the equation of motions imply the additional correlations between sum rule specific input parameters like the threshold parameters. This might be problematic because the latter parameters are intrinsic parameters of each separate LCSR calculations. However, these additional correlations are introduced in a two-step procedure: First, separate LCSR calculations are done without these additional correlations between the sum rule specific parameters in order to verify that the intrinsic threshold parameters of the LCSR of specific tensor and vector form factors have to be equal up to some variations and that also ratios of the corresponding form factors are equal to one. In a second step these correlations are implemented in the direct LCSR calculation of all QCD form factors. These additional correlations have analogous implications as the form factor relations in the large energy limit [69]. We will also test the impact of these additional correlations on our final theoretical predictions.
An explicit numerical comparison between the form factors calculated in Ref. [19] (KMPW) and in Ref. [72] (BSZ) is given in Appendix B. Moreover, a comparison for some of the angular observables when using these two different sets of form factors, as well as when employing the two different theoretical approaches (soft FF and full FF) have been given in Appendix C.

Error estimation
We have used a Monte Carlo program taking 200,000 random points for the error calculation: • All the input parameters are varied with Gaussian distributions in their 1σ range. ), respectively.
• The form factor parameters are varied using distributions reproducing the correlation matrices.
• In the soft FF approach for B → K * ¯ in the low-q 2 region, the factorisable and nonfactorisable power corrections have been considered collectively at the transversity spin amplitude level. The amplitudes are varied according to ranges which we refer to as 10, 20 and 30% error for the power corrections, respectively. For the high-q 2 region the multiplicity factor that we use is 3) for the 10, 20 and 30% errors, respectively.
• For B s → φ¯ and B → K¯ error estimation we have used a procedure similar to the B → K * ¯ decay.
• We compute the theoretical errors and correlations between all the observables which are given in the appendices and we combine them with the experimental correlation matrices to obtain the χ 2 .
We finally note that fitting a polynomial ansatz for the non-factorisable power corrections as it has been done in Ref. [34] and varying randomly the coefficients of the polynomial for an estimation of the power correction error are two different approaches. However, if the nonfactorisable power corrections needed in the fit and the power correction errors used in the FF approach are of the same magnitude, one expects a similar compatibility with the SM predictions. As mentioned in the introduction, the fit to the data presented in Ref. [34] needs rather large non-factorisable power corrections in the critical bins. A comparison of the fitted theory predictions with the leading contributions based on QCD factorisation finds that the power corrections needed in the fit within the critical bin from 4 GeV 2 to 6 GeV 2 are at the level of observables for example +23% in S 3 , -42% in S 5 , and -44% in S 4 relative to the leading contribution. On the other hand, a 5%, 10%, or 20% estimation of the non-factorisable power correction error described above (following Ref. [18]) leads to an error of maximal 1.5%, 3% or 6% at the observable level in S 3 , S 4 and S 5 respectively. Nevertheless the 60% error estimation in this framework leads to 17%-20% error in these three observables. And a 150% estimation of the power correction error is needed to reproduce errors up to 50% at the observable level which are comparable in size with the corrections needed in the fit presented in Ref. [34]. As a consequence, one should not expect a large impact of the 5%, 10%, or 20% power correction error in our global analysis.

Statistical methods
We use the absolute χ 2 method to verify the goodness of the fit in the first step. In a second step we consider the difference of the χ 2 with the minimum χ 2 to obtain the allowed regions for the Wilson coefficients (∆χ 2 ). We always scan over a specific number of Wilson coefficients δC i at the µ b scale and include all available correlations.
We also directly obtain the allowed regions from the absolute χ 2 . This procedure leads to larger allowed regions for the Wilson coefficients with respect to the use of the ∆χ 2 . One reason for this is that some of the observables are less sensitive to some Wilson coefficients, while they contribute in a democratic way to the number of degrees of freedom. However, the statistical meaning of the two dimensional contours within the absolute χ 2 is that for a point in the 1σ interval allowed region, there is at least one solution with the corresponding values of the Wilson coefficients that has a χ 2 probability corresponding to less than one Gaussian standard deviation with respect to the full set of measurements (see Ref. [16] for more details).
Within our analysis below we will arrive at the case that the absolute χ 2 fit leads to no compatibility at 68% C.L. and the 95% C.L. regions are small. In this case, the ∆χ 2 metrology does not make much sense. This indicates the possible role of the absolute χ 2 fit as a check of the goodness of the fit.

Observable dependency on Wilson coefficients
The tensions between the experimental measurement for B → K * µ + µ − observables and their SM predictions with 10% non-factorisable power corrections can be explained by modified Wilson and there are those which are dependent to only certain Wilson coefficients such as S 3 which is much more sensitive to primed Wilson coefficients, most specifically to C 7 in q 2 3 GeV 2 range [94] and less so to C 10 in the q 2 3 GeV 2 region. In addition to the B → K * µ + µ − observables, the other b → s transitions are also dependent on the Wilson coefficients. Of course, in order to get the best agreement with data one should find the best value for the Wilson coefficient(s) through a fitting with a method such as the method of least squares.

Global fit results assuming new physics in one operator only
Considering that new physics effects only appear in one operator, we make a χ 2 fit by scanning over a single Wilson coefficient while keeping the other Wilson coefficients to their SM values. We here use the full FF approach with a 10% power correction error. In the case of lepton flavour universality (C µ i = C e i ), an 18% reduction to C 9 gives the most probable scenario with a χ 2 min of 123.8. Assuming this scenario to be the correct description of the b → s data, the SM value for C 9 (corresponding to δC 9 = 0) is in 3.0σ tension with the best fit value (Pull SM ). On the other hand, considering contributions from electrons and muons to be different, the most probable scenario is for C µ 9 to receive a 21% reduction compared to the SM value for C 9 . The χ 2 min of this scenario is 115.5, which is significantly reduced compared to the one in which lepton universality is assumed. Here the SM value is in 4.2σ tension with the best fit value of C µ 9 . The best fit values of the different single Wilson coefficient fits as well as their 68 and 95% confidence level regions are given in Table 2.  Table 2: Best fit values and the corresponding 68 and 95% confidence level regions in the one operator global fit to the b → s data as described in the text. In the last two rows the χ 2 fits are done when considering lepton non-universality.

Global two operator fit
In this section we have considered new physics effects to appear in two operators by varying {C 9 , C 10 }, {C 9 , C 9 } and {C µ 9 , C e 9 } separately.

Fit results for {C 9 , C 10 }
We have first obtained the best fit value and the corresponding 1 and 2σ allowed regions when varying the {C 9 , C 10 } set of Wilson coefficients, where for the theoretical predictions of the observables we have used both the full FF and the soft FF approaches in Fig. 1 and Fig. 2, respectively. In the full FF approach, the allowed regions for the {C 9 , C 10 } set are similar when considering 5, 10 and 20% error for the power corrections which indicates that up to 20%, the power corrections are sub-dominant compared to the other theoretical uncertainties. As derived in Section 3.2, the 5, 10 and 20% error for the power corrections at the amplitude level leads to a 6% error at maximum at the observable level only. This is not the case for the soft FF approach since the power corrections affect both the factorisable as well as non-factorisable corrections, while in the full FF approach only the nonfactorisable part is affected. However, in both methods, considering the power corrections to Figure 3: Global fit results for C 9 , C 10 by using full form factors, with ∆χ 2 method and 10% power correction error. On the left plot, the 2σ contours when removing the form factor correlations, as well as when doubling (quadrupling) the form factor errors are shown with solid and dashed (dotted) lines, respectively. On the right plot, the solid lines correspond to the 1 and 2σ contours when considering 60% power correction error. Figure 4: Global fit results for C 9 , C 10 by using full form factors and 10% power correction errors, with absolute χ 2 . For B → K + µ + µ − only the low-and high-q 2 bins are used in the plot on the left while all the small bins are used in the plot on the right. be up to 20%, the SM is disfavoured at more than 2σ. For example, assuming a 10% power correction error within the full FF method leads to a SM pull of 2.6σ, (meaning that the SM value is in 2.6σ tension with the best fit values of C 9 and C 10 ). For comparison we also show the global fits to the Wilson coefficients based on the absolute χ 2 method. As expected (see Section 3.3) they lead to weaker bounds.
In the left plot of Fig. 3 the 1 and 2σ allowed regions are demonstrated, when removing the form factor correlations or doubling as well as quadrupling the form factor errors. The size of the form factor errors has a crucial role in constraining the allowed region; doubling (quadrupling) the error decreases the tension from 2.6σ to 2.1σ (1.4σ). But removing the form factor correlations does not have a significant impact. This is due to the small uncertainties of the BSZ form factors. If the quadruple form factor errors were considered, then the correlations would play a more important role. Assuming a 60% power correction error in the global fit has not a big impact either as the right plot of Fig. 3 shows. As discussed above, such a guesstimate of the non-factorisable power correction leads to errors of 20% at the observable level. That the best fit point gets slightly moved away from the SM point is a consequence of how such a guesstimate is implemented at the amplitude level, namely within terms without dependences Figure 5: Global fit results for C 9 , C 9 , with ∆χ 2 method and by considering 10% power correction error. On the left the full form factor method has been used, where the solid line corresponds to 5% power correction errors and the dashed line to 20% power correction errors. On the right the soft form factor method has been used, where the solid line corresponds to 20% power correction errors and the dashed line to 30% power correction errors. The solid and dashed lines are at 2σ. Figure 6: Global fit results for C 9 , C 9 by using full form factors, with ∆χ 2 method and 10% power correction errors. On the left plot, the 2σ contours when removing the form factor correlations, as well as when doubling (quadrupling) the form factor errors are shown with solid and dashed (dotted) lines, respectively. On the right plot, the solid lines correspond to the 1 and 2σ contours when considering 60% power correction error. on the Wilson coefficients C 7 , C 9 , and C 10 (see Section 3.2).
Finally, we analyse the global fit to C 9 and C 10 with 10% power correction using the absolute χ 2 method when we do not use only one global low-q 2 and one global high-q 2 , but also smaller bins of the observable B → K + µ + µ − . The right plot in Fig. 4 shows the latter option. We see that such a fit does not lead to any compatibility at 68% C.L. which clearly indicates that the fit is not very good, when the data with smaller bins are used. This reflects the well-known observation that the large resonance structure in the high-q 2 region of this observable does not allow for a theoretical description of this observable with smaller binning.

4.3.2
Fit results for {C 9 , C 9 } In Fig. 5 we give the 1 and 2σ allowed regions for the fit to the {C 9 , C 9 } set within both the full FF as well as soft FF approaches. Assuming 10% power correction for the SM predictions Figure 7: Global fit results for C e 9 , C µ 9 by using full FF approach and considering 10% power correction when employing the ∆χ 2 (left) and absolute χ 2 (right) methods. The solid and dashed lines correspond to the allowed regions at 2σ when considering 5 and 20% power corrections, respectively. Figure 8: Global fit results for C e 9 , C µ 9 by using soft FF approach and considering 10% power correction when employing the ∆χ 2 (left) and absolute χ 2 (right) methods. The solid and dashed lines correspond to the allowed regions at 2σ when considering 20 and 30% power corrections, respectively.
of B → K ( * ) (φ)¯ decays, the SM value has a slight tension of more than 2σ with the best fit point (2.7σ (2.1σ) with χ 2 = 123.1 (118.9) for the full (soft) FF approach), where the tension is mostly in C 9 . Compared to the {C 9 , C 10 } fit, the very similar χ 2 values indicate that there is no preference between the fits when considering new physics effects in C 10 or C 9 . The effects of doubling and quadrupling the form factor errors and of a 60% power correction error as shown in Fig. 6 are also similar to the {C 9 , C 10 } fit.

Fit results for {C e
9 , C µ 9 } In Figs. 7,8, the allowed regions of the fit at 1 and 2σ are presented when considering different contributions for muons and electrons in C 9 . Using the full (soft) FF approach with 10% power correction, the χ 2 of the best fit point is 114.6 (109.8), which indicates a considerable improvement of the fit compared to the {C 9 , C 10 } and {C 9 , C 9 } fits. The SM value which respects lepton universality has a tension of 3.9σ (3.6σ) with the best fit point where most of the tension appears in C µ 9 . Within the full FF approach, the Pull SM is reduced from 3.9σ to 3.1σ only by quadrupling the form factor error within the full FF approach (see Fig. 9).

Role of R K and S 5
While the experimental measurement of S 5 and its tension with the SM prediction seems to be the main reason for a best fit point of C 9 about 20% less than its SM value, the S 5 turns out to be not the only observable which drives δC 9 /C SM 9 to negative values. In Fig. 10 we have given the two operator fits when removing the data on S 5 while keeping all the other b → s data. It can be seen that while the tension of C SM 9 and best fit point value of C 9 is slightly reduced in the various two operator fits, still the tension exists at more than 2σ. This feature indicates that S 5 is not the only observable that drives δC 9 to negative values and other observables play a similar role. However, the situation is different in the case of the observable R K . In Fig. 11 from the lower right figure it can be seen that R K is the main measurement resulting in the best fit values for C µ 9 and C e 9 which are in more than 2σ tension with lepton-universality (and with 3.9σ tension with the SM). Removing the data on R K from the fit, lepton-universality can be restored at slightly larger than 1σ. So at present R K is the only observable within the b → s + − transitions which shows some sign of lepton non-universality.

Fit results considering only
The experimental measurements of the B → K * µ + µ − angular observables have been obtained using the method of moments as well as the most likelihood method. While the former gives less precise results, it is more robust compared to the latter one, specially for low signal yields. Figure 12: Global fit results with the ∆χ 2 method, using only the B → K * µ + µ − data in the full FF approach with 10% power correction. The coloured regions correspond to the best fit value using the method of moment results at 1 and 2σ. The solid lines depict the 1 and 2σ allowed regions using the likelihood method results for the B → K * µ + µ − observables.
The new physics analysis clearly depends on the method which has been used to obtain the B → K * µ + µ − experimental measurements. To compare the best fit points when using the two different results we have done two operator fits using only the data on B → K * µ + µ − observables. In Fig. 12, the coloured areas show the allowed regions at 1 and 2σ when using the method of moment results while the solid and dashed lines correspond to the 1 and 2σ contours by using the method of most likelihood results. We see that when using the method of moment results, the agreement of the fit with the SM is better. E.g., in the C 9 , C 10 plane the SM has a pull of 1.7σ with the best fit point when using the method of moments results, compared to 3.3σ for the most likelihood results. The smaller tension between the SM value and the best fit point within the method of moments is partly due to small shifts in the central values and mostly due to the larger experimental errors associated with the results of this method.

Global four operator fits
There is in principle no reason to assume that new physics shows up only in one or two Wilson coefficients. Hence, it is of importance to check the global agreement of the experimental data when allowing new physics contributions to several Wilson coefficients at the time. In this section we consider four operator fits, where in all the cases the full FF approach with 10% power correction is used.

Global fit results for four operators {C e
9 , C µ 9 , C e 9 , C µ 9 } Assuming new physics to appear in the O 9 and O 9 operators with different contributions for the electron and muon sectors, we have fitted the {C e 9 , C µ 9 , C e 9 , C µ 9 } set of Wilson coefficients to the b → s data. In this four operator fit the SM has 3.4σ tension with the best fit point which has a χ 2 of 113.3. In Fig. 13, a projection of the allowed regions at 1 and 2σ on different Wilson coefficient planes are shown. In the upper right plot of Fig. 13, the projection of the 1 and 2σ regions in the {C e 9 , C µ 9 } plane is given, where each point within the coloured region indicates that there exist at least some values of C e 9 and C µ 9 for which the corresponding C e 9 and C µ 9 value give a χ 2 that is within the 2σ region. In this plot, besides the projected four operator fit, we have overlaid the 1 and 2σ contours of the two operator fit of {C e 9 , C µ 9 } for comparison. The comparison of the results shows that considering only the modification of two Wilson coefficients leads to much more restrictive results. And while in the latter, lepton universality is in more than 2σ tension with the best fit point, in the four operator fit lepton universality in C 9 is respected even within the 1σ region. However, in this case the lepton non-universality is appearing in C 9 which is not shown in the projected plane.

Global fit results for four operators {C e
9 , C µ 9 , C e 10 , C µ 10 } We consider here the four operator fit assuming lepton non-universality in O 9 and O 10 . The best fit point of the {C e 9 , C µ 9 , C e 10 , C µ 10 } fit has a χ 2 = 114, which indicates that there is no improvement in the fit when replacing the operator O 9 with O 10 , even when considering different contributions for muons and electrons (the same result was also seen in the two operator fit with lepton universality). In the {C e 9 , C µ 9 , C e 10 , C µ 10 } fit the SM value has a pull of 3.4σ with respect to the best fit point, similar to the {C e 9 , C µ 9 , C e 9 , C µ 9 } fit. Some of the possible two dimensional projections of the four operator fit are presented in Fig. 14. In the upper left plot, the 1 and 2σ contours of the {C e 9 , C µ 9 } two operator fit has also been shown. A comparison between the allowed regions in the two and four operator fits shows that considering four operator fits considerably relaxes the constraints on the Wilson coefficients leaving room for more diverse new physics contributions which are otherwise overlooked. 4.6.3 Global fit results for four operators {C 9 , C 10 , C 9 , C 10 } In Fig. 15, the projection of the {C 9 , C 10 , C 9 , C 10 } fit on different 2-dimensional planes are demonstrated. This four operator fit has a best fit point with χ 2 = 121.6 which indicates that the experimental measurements are better described assuming lepton non-universality as in the two previous subsections. In this case the SM value of the Wilson coefficients has a pull of 2.3σ with the best fit point. Including the primed operators with respect to the two operator fit for {C 9 , C 10 } (with χ 2 = 123.7) does not improve the fit 5 (see also the upper left plot of Fig. 15). The two-operator fits are overlaid again in the projection plot of the four-operator fit. The comparison shows that the bounds based on the two-operator fits are always stronger by construction.

Global fit results in MFV
In this section we show the impact of the b → s data within the framework of minimal flavour violation (MFV), see e.g. [95][96][97][98][99] and [100] for a recent review. There are different definitions for the MFV framework. We follow the canonical one which is based on a symmetry principle introduced in Ref. [97], which implies that in a MFV model all flavour-violating interactions can be traced back to the well-known structure of the Yukawa couplings and all CP-violating interactions are due to the physical phase in the CKM matrix.
The specific hierarchy of the quark masses and CKM matrix elements implies that only a small number of operators are relevant in the MFV framework compared to a general modelindependent analysis. Here, besides the operators present in the SM, especially O 7 , O 8 , O 9 , and O 10 , we also have to consider the scalar-density operator with right-handed b-quarks (O l = e 2 /16π 2 × (s L b R )(¯ R L )) (see Ref. [99] for more details).
The MFV hypothesis represents an important benchmark in the sense that any measurement which is inconsistent with the general constraints and relations induced by the MFV hypothesis unambiguously indicates the existence of new flavour structures. Thus, any incompatibility of the b → s measurements with the MFV hypothesis would imply that new flavour structures are needed to explain the data.
We have made a global fit within the MFV framework using the five operators listed above. In addition to the observables used in the previous fits, we include BR(B → X s γ) and the isospin asymmetry of B → K * γ, which are sensitive to the O 7 and O 8 operators. The best fit point has a χ 2 of 123.5 for 137 degrees of freedom which represents a good fit. The five operator fit within the MFV framework shows compatibility with the MFV hypothesis. In Fig. 16, the resulting bounds on the Wilson coefficients are shown.

R K and predictions for other ratios
In the introduction, we have identified the observable R K as a possible key observable to clarify also the origin of the anomalies in the LHCb data. The ratio R K = BR(B + → K + µ + µ − )/ BR(B + → K + e + e − ) in the low-q 2 region had been measured by LHCb using the full 3 fb −1 of data, showing a 2.6σ deviation from the SM prediction [35]. This might be a sign for lepton non-universality. In contrast to the anomalies in the angular observables in B → K * µ + µ − , the ratio R K is theoretically rather clean, in particular it is unaffected by power corrections and also the electromagnetic corrections are under control. Because both tensions might be explained by new physics in the Wilson coefficients C 9 , other ratios of observables which may indicate lepton non-universality will be important crosschecks of the anomalies discussed in this paper. In Section 4.4 we found that at present the R K ratio is the only driving force for lepton non-universality.
In Ref. [89], the authors predict the central values for such ratios based on a specific choice of the Wilson coefficients, in particular on the assumption that the electron modes are SM like. In contrast, we make our predictions for such ratios based on the global fit considering two Wilson coefficients C µ 9 and C e 9 . We find that in most cases the SM point is outside the 2σ region of our indirect predictions reflecting the present tension in R K (see Table 3). Moreover, the ratio in the case of the A F B looks most promising from the theoretical point of view.
Finally we note that we have shown previously [16,31,32], that the present data on inclusive and exclusive decays are compatible with each other and there is no sign of lepton non-universality in the published data on the inclusive mode.

Conclusions
The LHCb collaboration has recently presented new data on exclusive b → s + − penguin decays [1,2]. These data were eagerly awaited because of some tensions with the SM in the angular observables of B → K * µ + µ − . In view of these new data we have stressed again that there is a significant difference in the theoretical accuracy of the inclusive and exclusive b → s + − decays in the low-q 2 region. The theoretical description of power corrections exists in the inclusive case so that they can be at least estimated within the theoretical approach. On the contrary, no theoretical description of power corrections exists in the exclusive case in the framework of QCD factorisation and SCET, and power corrections can only be guesstimated. This issue makes it rather difficult or even impossible to separate new physics effects from such potentially large hadronic power corrections within these exclusive angular observables. Therefore, these tensions might stay unexplained until Belle II will clarify the situation by measuring the corresponding inclusive b → s + − observables as we have demonstrated previously [16,31,32].
The present situation motivated a recent theory analysis in which the unknown power corrections were just fitted to the data using an ansatz with 18 additional real parameters in the fit [34]. However, we have shown that this fit to the data needs very large power corrections from 20% up to 50%, and even larger, in the critical bins of the angular observables. The existence of such large hadronic corrections cannot be ruled out in principle, but they somehow question the validity of the QCD factorisation approach for such observables.
We have analysed how the tensions in the present LHCb data depend on the input parameters like form factor calculations and corresponding correlations and have shown that they are rather insensitive to these inputs. Only quadrupling the form factor error makes a real difference. We have also found that the standard guesstimate of non-factorisable power corrections from 5% to 20% at the amplitude level has no real impact on the theoretical predictions. As we have shown such variations lead to a 6% error at the observable level for the three observables S 3 , S 4 and S 5 . Only variations significantly larger than 60% -corresponding to errors larger that 20% at the observable level -have a real impact. For example we have shown that for the three observables S 3 , S 4 and S 5 , a 150% of the power correction error is needed in the critical bins to reproduce errors up to 50% at the observable level which are comparable in size with the corrections needed in the aforementioned recent fit to the SM.
In addition, we have explicitly shown that within a new physics analysis the observable S 5 is not the only observable that drives the new physics Wilson coefficient δC 9 to negative values, but also other observables play a similar role.
We have noted that the other tensions in the LHCb data of b → s may play a crucial role in the near future, namely the observable R K with a 2.6σ deviation from the SM prediction [35] what might signal a first sign of lepton non-universality. In contrast to the anomaly in the rare decay B → K * µ + µ − which is affected by power corrections, the ratio R K is theoretically rather clean and its tension with the SM cannot be explained by power corrections. Neither some missing electromagnetic corrections can serve as explanation for the discrepancy from the SM.
It might be just an accidental coincidence that the tensions in R K and in the angular observables can simultaneously be resolved by a negative new physics contribution to the Wilson coefficient of the semileptonic operator O 9 . If not, then the measurement of analogous ratios which may show signs of lepton non-universality will be also an important crosscheck of the anomalies in the angular observables and might even resolve the puzzle. Therefore, we have presented predictions for such ratios based on our global fits.

A SM predictions and experimental values
The experimental values and SM predictions of the observables considered in this work have been given in Tables 4-10. The angular observables in Tables 7-9 0.60 ± 0.31 [105]  B → Kµ + µ − differential branching ratios Bin (GeV 2 ) SM prediction Measurement 26.6 ± 5.7 12.2 +5.9 12.1 ± 0.4 ± 0.6 Table 5: SM predictions and experimental values for the differential branching ratio of the B → Kµ + µ − decay. The uncertainties of the experimental values [106] are (from left to right) statistical and systematic.
B → K * µ + µ − differential branching ratios Bin (GeV 2 ) SM prediction Measurement 10 9 × dBR/dq 2 (B + → K * + µ + µ − ) 46.7 ± 6.0 11.6 9.1 −7.6 ± 0.8  Table 6: SM predictions and experimental values for the differential branching ratio of the B → K * µ + µ − decay. The experimental error of the B + → K * + µ + µ − decay [106] are (from left to right) statistical and systematic. The experimental error of the B 0 → K * 0 µ + µ − decay [107] have been added in quadrature, taking the largest side error in case of non-symmetrical uncertainties.      Table 11: Comparison of the SM predictions for BR(B 0 → K 0 µ + µ − ) from SuperIso using the full FF approach and assuming 10% power correction with the result of Table 4 of Ref. [20].  Figure 17: LCSR results for the B → K * form factors. In the upper row the KMPW results [19] are shown while in the lower row the BSZ results [72] are presented. The relation among T 23 and T 2,3 as well as the relation between A 12 and A 1,2 are defined in Ref. [84].

B Form factors
The methods used for obtaining form factor results depend on the recoil energy of the outgoing light meson. At high-q 2 which corresponds to the low recoil region, the B → K * and B s → φ form factor results are available from unquenched lattice calculations [84,85], while for the low-q 2 region the form factors can be taken from LCSR calculations which are available from Ref. [19] (KMPW) as well as from the Ref. [72] (BSZ). The theoretical errors associated with the KMPW form factors are larger than the ones from BSZ. The larger error of the KMPW form factor is mostly by construction and due to the different choices of distribution amplitudes where they employ the B-meson distribution for which less information is available compared to the K * meson which has been used for the BSZ form factors. Moreover, for the BSZ form factors, interpolation with lattice data and additional correlations between LCSR intrinsic parameters have been used which also have some effect in reducing the theoretical uncertainty.
The form factor uncertainties are correlated through hadronic inputs as well as kinematic relations at the endpoint q 2 =0. As claimed in Ref. [72] there are additional correlations due to intrinsic LCSR parameters. While in the soft FF approach due the relations among form factors at high recoil energy the number of independent form factors reduces from seven to two form factors, and most of the latter correlations are included by construction, in the full FF approaches analogous implications can be derived directly within the LCSR results via these additional correlations mentioned above.
Unfortunately the form factor correlations have not been given for the KMPW results, however, they have been provided in Ref. [72] for the BSZ form factors.
The LCSR results for the seven independent B → K * form factors including their theoretical uncertainties are shown in Fig. 17, where the KMPW form factors which are applicable only at The blue crosses represent the LHCb measurements [1]. The solid and dotted blue lines correspond to SM predictions using BSZ form factors within the full form factor and soft form factor approaches, respectively. The solid and dotted red lines correspond to SM predictions using KMPW form factors within the full form factor and soft form factor approaches, respectively. low q 2 have been extrapolated to high-q 2 as well. For the BSZ form factors we have presented the fit results of Ref. [72] which are applicable for both the low-and high-q 2 regions.
C SM predictions of B → K * µ + µ − observables with different theoretical approaches and form factors In order to compare with the results of Ref. [22], we reproduced the SM predictions using KMPW form factors within the soft FF approach and we also added the fit results of Ref. [18] which are meant to take into account the missing 1/m b factorisable corrections. These predictions coincided nicely with central values quoted in Ref. [22] up to deviations in observables which are very small or in the bins where there is a zero-crossing (the difference being mostly due to slight disagreements in the choice of SM Wilson coefficient values). A similar comparison was done with the results of Ref. [72], when using the BSZ form factors in the full FF approach. Again the results were in good agreement, and only in the bins with zero crossings and for observables having very small values, there were some slight differences.
In Figs. 18,19, we have presented SM predictions of the central values for the most relevant B → K * µ + µ − observables, employing the soft FF and full FF approaches both when using KMPW and BSZ form factors (Appendix B). In these figures, for the soft FF approach we have not included the 1/m b power corrections which have been estimated through fitting with ad hoc functions in Ref. [18] (for the KMPW form factors). In Figs. 18,19 it can be seen that while there are good agreement between the different approaches and the different form factor choices, the significance of the tension between central value of the SM predictions and the experimental data depends on the particular choice of the theoretical approach as well as which set of form factors are used. E.g., both for the SM prediction of S 5 in Fig 18 and    . The black crosses correspond to the LHCb measurements where dBR/dq 2 is from the 1 fb −1 of data [107] and the angular observables are from the 3 fb −1 of data [1]. The blue bands correspond to the binned SM predictions with their relevant uncertainties.