Global analysis of b → sℓℓ anomalies

We present a detailed discussion of the current theoretical and experimental situation of the anomaly in the angular distribution of B → K*(→ Kπ)μ+μ−, observed at LHCb in the 1 fb−1 dataset and recently confirmed by the 3 fb−1 dataset. The impact of this data and other recent measurements on b → sℓ+ℓ− transitions (ℓ = e, μ) is considered. We review the observables of interest, focusing on their theoretical uncertainties and their sensitivity to New Physics, based on an analysis employing the QCD factorisation approach including several sources of hadronic uncertainties (form factors, power corrections, charm-loop effects). We perform fits to New Physics contributions including experimental and theoretical correlations. The solution that we proposed in 2013 to solve the B → K*μ+μ− anomaly, with a contribution C9NP≃−1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal{C}}_9^{\mathrm{NP}}\simeq -1 $$\end{document}, is confirmed and reinforced. A wider range of New-Physics scenarios with high significances (between 4 and 5 σ) emerges from the fit, some of them being particularly relevant for model building. More data is needed to discriminate among them conclusively. The inclusion of b → se+e− observables increases the significance of the favoured scenarios under the hypothesis of New Physics breaking lepton flavour universality. Several tests illustrate the robustness of our conclusions.


Introduction
Flavour-Changing Neutral Currents (FCNC) have been prominent tools in high-energy physics in the search for new degrees of freedom, due to their quantum sensitivity to energies much higher than the external particles involved. In the current context where the LHC has discovered a scalar boson completing the Standard Model (SM) picture but no additional particles that would go beyond this framework, FCNC can be instrumental in order to determine where to look for New Physics (NP). One particularly interesting instance of FCNC is provided by b → s and b → sγ transitions, which can be probed through various decay channels, currently studied in detail at the LHCb, CMS and ATLAS experiments. In addition, in some kinematic configurations it is possible to build observables with a very limited sensitivity to hadronic uncertainties, and thus enhancing the discovery potential of these decays for NP, based on the use of effective field theories adapted to the problem at hand. Finally, it is possible to analyse all these decays using a model-independent approach, namely the effective Hamiltonian [1,2] where heavy degrees of freedom have been integrated out in short-distance Wilson coefficients C i , leaving only a set of operators O i describing the physics at long distances: (up to small corrections proportional to V ub V * us in the SM). In the following, the factorisation scale for the Wilson coefficients is µ b = 4.8 GeV. We focus our attention on the operators

JHEP06(2016)092
section 3. In section 4 we discuss a set of scenarios with large NP contributions to one or two Wilson coefficients, confirming that a negative contribution to C 9 yields a significant improvement compared to the SM. We discuss which of these scenarios are able to reduce the anomalies observed in b → s transitions. By performing a global fit to all six Wilson coefficients simultaneously, we show that the most economic scenarios do indeed capture the main patterns suggested by the data. In this case we provide, in addition, confidencelevel regions for all Wilson coefficients when all of them are allowed to deviate from their SM values simultaneously. We also consider scenarios with violation of lepton-flavour universality, and describe tests of the robustness of the fits presented. In section 5, we provide tests of the various sources of hadronic uncertainties that could affect our results (choice of form factors, power corrections, long-distance charm corrections). We present our conclusions in section 6. Appendices A and B are devoted to tables presenting our predictions for the SM as well as the best-fit point for NP in C 9 only. In appendix C, the confidence regions for less favoured, but theoretically interesting, scenarios are shown. Appendix D describes how various changes in the analysis affect its outcome for the scenario with NP in C 9 only. Appendix E contains further details on power corrections to B s → φ and B → K form factors. Appendix F gathers basic features of Z models relevant for the b → s anomalies.

General approach
In the effective Hamiltonian approach and in the SM (the extension to NP operators is straightforward), the B → K * µµ transversity amplitudes can be written in a compact way as A ∝ C 7 2im b q 2 q ρ K * |sσ ρµ (1 + γ 5 )b|B + C 9 K * |sγ µ (1 − γ 5 )b|B + H µ ū γ µ v +C 10 K * |sγ µ (1 − γ 5 )b|B ū γ µ γ 5 v , (2.1) with H µ ∝ i d 4 x e iq·x K * |T [cγ µ c]H c |B , (2.2) where H c denotes the part of the weak effective Hamiltonian involving four-quark operators with two charm fields. For simplicity, we have neglected contributions from CKMsuppressed terms here (they are included in our numerical evaluations). One can see from eq. (2.1) the existence of two different kinds of contributions: local ones yielding form factors (seven for B → K * ) and non-local ones (involving cc loops propagating). The former can be determined using non-perturbative methods (light-cone sum rules, lattice), whereas the latter must be estimated using 1/m b expansion (QCD factorisation, OPE), with different tools depending on the kinematic regime considered (large-or low-K * recoil). We will illustrate these points in the large-recoil region where the strongest deviations have been observed between SM predictions and data. A first step in the evaluation of the amplitudes comes from the contributions due to O 7,9,10 , involving seven form factors. In the large-recoil region there are basically two approaches:

JHEP06(2016)092
• "Improved QCD Factorisation (QCDF) approach": in this framework [7] the largerecoil symmetries between form factors are used to implement the dominant correlations among them. This general approach is easy to cross-check and to implement for any form factor parametrisation (e.g. for the light-cone sum rules parametrisations [17,20,39]). The symmetries allow the 7 form factors to be written in terms of only two so-called soft form factors ξ ⊥, [40]: To this soft-form factor representation one should add (perturbatively computable) hard-gluon O(α s ) corrections as well as (non-perturbative) O(Λ/m b ) corrections [41]. The soft form factors can be computed in a specific parametrisation. The basis of optimized observables P i is usually taken in this approach [5][6][7][42][43][44]. We follow ref. [23] where we considered all symmetry-breaking corrections to the relations in eq. (2.3). Our predictions take into account factorizable α s -corrections computed within QCDF [41,45,46], as well as factorizable power corrections. We will consider most of the time the full form factors of ref. [17], but for completeness we will also compare some of our results with the results using the form factors in ref. [20].
• "Full Form Factor approach": here a specific set of full form factors determined from light-cone sum rules [20,39] is used. Factorizable α s and factorizable power corrections are automatically included with correlations associated to this particular parametrisation. Other corrections to the amplitudes (non-factorisable pieces, see below) have to be included and/or estimated exactly as in the previous approach. This approach has been employed in refs. [10,15,47].
Both approaches are useful and complementary, should converge and give comparable results and error sizes, as long as the correlations among the form factors are dominated by the large-recoil relations. It is interesting to notice that the relevant form factors for the transversity amplitudes are not those defined in the usual transversity basis (V, A i , T i ) but rather the helicity form factors [21,48] being linear combinations of the usual transversity ones. It is therefore important to determine properly the correlations among the usual form factors in order to determine correctly the transversity amplitudes. The first approach allows one to restore correlations that are expected among the various form factors, even when these correlations were not given initially. The second one requires one to compute the complete set of form factors and to achieve a very good control of the applied theoretical method in order to determine a meaningful correlation matrix. Of course, both methods can be used to compute both types of observables P i and S i , and they are expected to yield similar results. We will discuss this point further in section 5.
Once the issue of the form factors has been settled, one can proceed with the determination of the amplitudes involving not only the form factors but also non-local cc loop contributions. QCD factorisation [41,45,46] yields an expression of the amplitudes in terms of JHEP06(2016)092 soft form factors, α s -and power corrections, which can be further split into factorisable and non-factorisable contributions (stemming or not from the expression of full form factors in terms of soft form factors). The factorisable power corrections have already been considered at the level of the form factors, whereas the non-factorisable ones still have to be addressed. First we take the three hadronic form factors T i (q 2 ) that parametrise the matrix element K * γ * |H eff |B [41], and we single out the hadronic contribution that is not related to the radiative Wilson coefficients (obtained taking the limit T had i = T i | C 7 ( ) →0 ). We multiply each of these amplitudes serving as a normalisation with a complex q 2 -dependent factor [23] T had We define our central values as the ones with r i (s) ≡ 0, and estimate the uncertainties from non-factorizable power corrections by varying r a,b,c i ∈ [0, 0.1] and φ a,b,c i ∈ [−π, π] independently, corresponding to a ∼ 10% correction with an arbitrary phase.
Part of the cc-loop contributions have been already included in the non-factorizable contributions (hard-gluon exchange). The remaining long-distance contributions from cc loops are still under debate. For these contributions we will rely on the partial computation in ref. [17]. It is important to remark that the soft-gluon contribution of ref. [17] coming from 4-quark and penguin operators induces a positive contribution to C eff 9 whose effect is to enhance the anomaly. Since we are interested only in the long-distance contribution δC LD 9 (q 2 ), we subtract the perturbative LO part and include the shift due to a different reference value for m c . Ref. [23] provides more details on this procedure. We introduce two different parametrisations, corresponding to the contribution to transverse amplitudes and to the longitudinal amplitude (which does not exhibit a pole at q 2 = 0) setting s 0 = 1 GeV 2 . We tune the parameters in order to cover the results obtained in section 7 of ref. [17] in the q 2 -region between 1 and 9 GeV 2 , where results for the three transversity amplitudes (denoted M 1 , M 2 and M 3 ) have been derived. 1 We get a ⊥ , a || = 9.25 ± 2.25 ,  where all parameters will be taken as uncorrelated. The resulting functions δC LD,(⊥,||) 9 (q 2 ) and δC LD,0 9 (q 2 ) are shown in figure 1. In order to be conservative, and in particular given the discussion on the sign of this contribution, we use the result of ref. [17] as an order of magnitude estimate, performing the following shift in each pair of transversity amplitudes with three independent parameters s i = 0 ± 1 (we recall that we include the perturbative cc contribution in C eff 9 and that the direct inclusion of the result from ref. [17] would correspond to choosing s i = 1).
For the low-recoil region [49][50][51], one can perform a similar analysis based on Operator Product Expansion and Heavy-Quark Effective Theory, or using directly form factors provided by lattice QCD simulations. In the following, we will use the latter approach for the computation of the observables at low recoil. In this region, one has also to deal with resonances such as those observed by LHCb in the data of the partner channel B + → K + µ + µ − . This observation prevents one from taking small bins afflicted by the resonance structures. In ref. [52] a quantitative estimate of duality violation is given. Unavoidably, one needs to use a model for this estimate, still the result is that the low recoil bin, integrated over a large energy range, gets a duality-violation impact of a few percent at the level of the branching ratio (estimated to 5% in ref. [53] or 2% in ref. [52]). It remains to be determined if this estimate also applies for angular observables in B → K * µµ. Moreover, the exact definition of the ends of the single large bin has some impact on the analysis in the framework of the effective Hamiltonian [54]. In order to take into account such effect of duality violation for angular observables and the sensitivity to the position of the ends of the bin, we add a contribution of O(10%) (with an arbitrary phase) to the term proportional to C eff 9 for each transversity amplitude. We notice that for all exclusive processes at low recoil, we include the NNLL corrections for b → s processes as described in ref. [55].

Optimised basis of observables: definition, properties and impact of data
The structure of the amplitudes at large recoil led to the construction of the optimised observables P i and P CP i [5][6][7][42][43][44] that exhibit a sensitivity to the soft form factors JHEP06(2016)092 suppressed by α s or Λ/m b . The observables that we consider can be found in appendix A, including the branching ratio, its longitudinal fraction F L and the optimised observables P i . As discussed in ref. [5,6], the optimised observables P i together with two additional (form factor dependent) observables exhaust the information provided by the angular coefficients. 2 These optimised observables have been measured by LHCb: the latest results incorporating the full 3 fb −1 of data collected during LHC run I can be found in ref. [31], which includes the results for the CP-averaged coefficients S i introduced in ref. [47], as well as the corresponding correlation matrices.
We should stress at this point that our definition of some optimised observables P i and CP-averaged angular coefficients S i differs from that adopted by the LHCb collaboration, due to two different issues. First, our convention for the angles to define theB →K * kinematics (identical to ref. [47]) differs from the LHCb choice. Refs. [58,59] provided the angular coefficients J i in terms of the transversity amplitudes using the LHCb convention. Comparing with the expressions in ref. [47], one can confirm that the two conventions can be related using This induces different signs in both conventions when the angular coefficients J i (and their CP-averaged versions S i ) are expressed in terms of transversity amplitudes, leading to the identification S LHCb 4,6c,6s,7,9 = −S 4,6c,6s,7,9 , (2.13) the other coefficients S i being identical in both conventions. Second, our definition of the optimised observables P i in terms of the angular coefficients J i is different from the definition used by the LHCb collaboration [60]. This induces further sign and normalisation differences when expressing P i in terms of transversity amplitudes, finally leading to 3 14) The presence of discrepancies with respect to the SM in the LHCb measurements at 1 fb −1 and 3 fb −1 can be interpreted as a sign of additional contributions to some of the Wilson coefficients. It is thus interesting to study the sensitivity of the P i observables to such shifts, see table 1. One can see interesting patterns, and in particular the global preference for a negative contribution to C 9 , as already observed with previous data [8] and in other frameworks [10,15,16]. We will now discuss the features of each of the P i observables in more detail, as well as the status of LHCb data for these quantities. The results given here are based on the final results provided in ref. [31]. We will focus on the results obtained using the maximum likelihood approach, and we will not consider the results obtained using the amplitude method discussed recently in ref. [61]. [15,19]

JHEP06(2016)092
(2.15) and P 8 . These observables are predicted tiny in the case of real NP contributions and are measured compatible with zero, so that this update of dictionary has no actual consequences on the results of the fit in ref. [8]. 4 In this definition and in the following ones in this section, it should be understood that each term is combined with the corresponding CP-conjugated term and the two leptonic chiralities are included (for instance, . In addition, we will ignore various factors of βµ ≡ 1 − 4m 2 µ /q 2 , which are important for the observables at very low q 2 . For precise definitions see [5][6][7], where also the bin-integrated observables are given. Evidently, we use the exact expressions in all the numerical results throughout the paper. Figure 2. Data (blue crosses) and SM prediction (red boxes) for P 1 , P 4 . The sources of uncertainties (added in quadrature) are shown as boxes in the following order from the center towards the outside: parametric, form factors, factorisable corrections, non-factorisable corrections, charm loop. P 1 is particularly well suited to detect the presence of right-handed currents. The lefthanded structure of the SM implies that a b quark in the helicity state −1/2 would produce an s quark in the same helicity state (neglecting the s quark mass), combined with the spectator quark to generate a K * meson in an helicity state −1 or 0, but not +1. The suppression of H +1 = (A + A ⊥ )/ √ 2 0 implies A ⊥ −A and consequently P SM 1 0. In an completely analogous manner, a b quark in the helicity state +1/2 leads to

JHEP06(2016)092
Deviations from this prediction would signal contributions from a new right-handed structure.
As seen in figure 2, all bins are consistent with the SM, however with very large error bars, so that no robust conclusion can be extracted from this observable with present data.. In table 1 we present the impact on P 1 [0.1,0.98] , P 1 [6,8] and P 1 [15,19] of shifting one of 10 at a time. This is useful to see the relative size of the impact and if a corresponding NP contribution improves or not the agreement with data. Only significant improvements towards data are indicated. As expected, shifting Wilson coefficients for the SM operators does not induce any sizeable change. On the other hand, P 1 exhibits a relatively large sensitivity to right-handed operators. In particular should be noted a high sensitivity to contributions to C 7 in the first bin [33] as compared to other coefficients and also to other bins.

P 4
The next observable that we would like to discuss is .
(2. 16) In conjunction with P 5 , P 4 establishes bounds on P 1 and enters consistency relations [62]. In particular, the bound is very efficient in two bins: [6,8] and low recoil. The preference of data for P 4 ≥ 1 in the [6,8] bin requires P 1 ≤ 0, in agreement with the 2015 LHCb data. Strictly speaking, this bound holds among the q 2 dependent observables, but it should also apply when the functions are only slowly varying (or almost constant) for the binned observables. As an illustration of the usefulness as a test on data of the bounds provided by eq. (2.17) we have checked which value would imply for P 1 the measured values of P 4 and P 5 at low recoil. Taking central experimental values for this illustrative example we find that P 1 should be roughly in the range −0.54 ≤ P 1 ≤ −0.44 which is the right ball park as compared to the central measured value As can be seen in figure 2, P 4 exhibits a perfect agreement with the SM in all bins, still with very large error bars. For completeness we provide also the bins [6,8] and [15,19] in table 1 to make manifest the lower sensitivity of this observable to shifts of Wilson coefficients (particularly at low recoil) as compared to other observables, a fact that should not downgrade its status to a mere "control" observable.

P 2
The definition is [5,7] . (2.18) This observable is the optimised and clean version of the forward-backward asymmetry, as illustrated in figure 3 where the difference in the size of the uncertainties is obvious. It was originally called A (re) T = 2P 2 and proposed in ref. [65]. P 2 measures a particular correlation between A FB and F L that is independent of form factors at LO, and combined with either A FB or F L shows a higher NP sensitivity than the pair {A FB , F L } itself.
The observable P 2 contains some important pieces of information, such as the position of its zero q 2 0 (identical to the zero of A FB ), the position of its maximum q 2 1 , and its maximum value P 2 (q 2 1 ). To leading order and assuming no contribution from right-handed currents, i.e. C i = 0, they are given by: and where for the position of the maximum we have neglected a term of O(Im(C eff 9 ) 2 ) following ref. [57]. These expressions illustrate that a NP contribution to C 9 and C 7 would shift both the zero and the maximum, but with a different magnitude. Moreover, the maximum can be also shifted by a contribution to C 10 . The NLO prediction in the SM for these quantities are: with P 2 (q 2 NLO 1 ) = 0.501 ± 0.004. In refs. [57,65], a NP contribution to C 7,9,10 was shown to shift the position of the maximum but not the value of its maximum that is fixed at P 2 (q 2 1 ) = 1/2. On the other hand, NP contributions to the chirally flipped operators would reduce the maximum below 1/2, even if not by a large amount. Unfortunately, a fluctuation of the F L [2. 5,4] bin has induced a large experimental error in the corresponding bin of P 2 . This will be cured with more data and a finer binning. Table 1 shows the sensitivity to shifts of Wilson coefficients for the two interesting [6,8] and low-recoil bins. It is clear the low sensitivity to NP of this observable at low-recoil, where the largest shift is only of +0.06. Indeed this is consistent with the perfect agreement of this observable with SM at low-recoil. Concerning the large-recoil bin, it is interesting to notice that the shifts of the Wilson coefficients pushing P 2 [6,8] towards the data also shifts P 2 [2,5,4] in the right direction (assuming that data is above the SM prediction), while all chirally flipped coefficients (positive or negative) always shift down this observable in this bin but by a relatively small amount. Finally, P 2 offers different consistency checks based on the relation [62] This relation is very useful to check the internal consistency of experimental or theoretical results for the observables. 5 Figure 4. Data (blue crosses) and SM prediction (red boxes) for P 5 . Same conventions as in figure 2.

JHEP06(2016)092
A first example is given by the bin [6,8] (or even [4,6]). Setting P 2 = − (with > 0) one immediately obtains from the previous equation Using central values for illustration and taking P 2 [6,8] [6,8] and P 5 [6,8] , as well as in the bin [4,6] (in agreement with LHCb data). A second example comes from the zero of eq. (2.22). In ref. [62], it was shown that the following relation should be fulfilled at the position q 2 0 of the zero of P 2 (or A F B ):

P 5
This observable is defined as [6,7] . (2.24) One can provide an interpretation of P 4 and P 5 based on the expression in terms of the two-dimensional complex transversity vectors n ⊥, ,0 (see ref. [5] for the definition of these JHEP06(2016)092 vectors defined in a basis of transversity amplitudes with left-and right-handed structure for the dimuons). If we assume for simplicity that the transversity amplitudes are real, these two observables can be understood as the "cosine" of the relative angle between the parallel (respectively perpendicular) transversity vector and the longitudinal one It is interesting to translate these expressions in the helicity basis by introducing two vectors based on the helicity h = −1 components of the K * : n (a) In the absence of right-handed currents (H +1 0), these observables correspond to the projection of the longitudinal helicity vector on one of the two negative helicity states, namely (2.26) Given the dominance of the left-handed part of the amplitude, this explains that P 4 and P 5 exhibit q 2 -dependences that are almost the reflection of each other with respect to the axis q 2 = 0. Of course, this discussion is only qualitative and the details on the role of the righthanded amplitude n a,b − are fundamental to assess the sensitivity of these two observables to semileptonic coefficients. P 5 exhibits the largest deviation with respect to the SM prediction in some bins, as seen in figure 4, corresponding to the so-called anomaly [8]. An illustrative exercise consists in determining how this observable can receive a large impact while keeping P 4 near the SM value (in agreement with data). 6 A numerical analysis allows one to identify two mechanisms to enforce a suppression of P 5 with respect to P 4 . The first mechanism relies on lifting the suppression of the right-handed amplitudes with respect to the lefthanded amplitudes and to profit from the relative minus sign between the two terms in the numerator of P 5 versus the plus sign in P 4 . The suppression of the right-handed amplitudes is due to the C SM 9 ∼ −C SM 10 cancellation, altered if the NP contribution to the Wilson coefficients does not follow the same direction. 7 The second mechanism is much more simple and relies on introducing a new physics contribution that suppresses A L ⊥ without affecting all other amplitudes.
In table 1 we show the sensitivity to shifts of Wilson coefficients for the [6,8] and lowrecoil bins. One can notice the large sensitivity of P 5 [6,8] to a change of only C NP 9 as compared to P 4 [6,8] in agreement with the data. Similar results are found for P 5 [4,6] albeit with a different importance. At low recoil, P 5 [15,19] exhibits a better sensitivity to NP than other observables in this region (though less than in the large-recoil region). This observable is already at 1 σ consistent with SM at low-recoil, but the shifts in Wilson coefficients improving the agreement with data at large recoil go into the opposite direction at low recoil. 6 In table 1, one can notice the large impact of a variation of C9 in P 5 compared to the negligible impact on P 4 in the bin [6,8]. 7 This can be easily seen using the large-recoil expression of the amplitudes. The numerator of P 2 4 contains a term proportional to C 2 10 that dominates and screens the partial cancellation between the C9 and C7 terms. There is no such C 2 10 term surviving in the numerator of P 5 , so that the partial cancellation between C9 and C7 suppresses P 5 with respect to P 4 . Figure 5. Data (blue crosses) and SM prediction (red boxes) for P 3 (top), P 6 (bottom left), P 8 (bottom right). Same conventions as in figure 2.

JHEP06(2016)092
2.2.5 P 3 , P 6 and P 8 These observables are defined as [6,7] , (2.27) and . (2.28) They are mainly sensitive to phases, either strong or weak, in the SM or beyond. Present data is compatible with the SM with huge error bars, including also a local fluctuation of around 2 σ in one bin of P 6 that will plausibly disappear with more data. This set of observables also are required to fulfill bounds like which is a natural extension of the bounds discussed in ref. [62]. Let us mention that a more direct way to test the presence of new weak phases is the measurement of the P CP i observables [7]. The still limited statistics of LHCb data requires taking the limit of massless leptons for the determination of angular observables. The impact of this assumption is completely negligible in all bins except for the lowest bin [0.1,0.98]. Once included in the computation, the lepton mass yields a sizeable effect, pushing the SM prediction in the direction of data for P 2 , P 4,5 and F L . Indeed, the first terms of the distribution at LHCb are given by which is modified once lepton masses are considered [57] 1 whereF T,L and F L,T are detailed in ref. [56]. 8 All our observables are thus written and computed in terms of the longitudinal and transverse polarisation fractions F L,T (2.32) However, LHCb measures F L from the expression eq. (2.30) without lepton masses, where the dominant term is the cos 2 θ K term. This means that the experimental analysis actually extractsF L , whereF (2.33) The difference between F L andF L has a negligible impact in all bins except for the bin [0.1,0.98]. We have recomputed the first bin of P 2 , P 4,5 usingF L instead of F L and imposing the LHCb conditionF T = 1 −F L . For these observables, the central value for the SM prediction is shifted towards the data The bin [6,8] Some recent analyses of B → K * µµ data [15,16] have discarded the [6,8] bin because of the proximity of the J/ψ resonance. It is obviously possible to perform analyses without this bin, as some judgement must be exerted to decide which observables are sufficiently well controlled to be included in the fit. However, we want to emphasise the role played by this bin in our analysis. The smooth behaviour of P 5 up to bin [6,8] does not support claims of extremely large charm-loop contributions inducing a positive contribution to C 9 which would affect mainly bins above 6 GeV 2 [19]. A direct comparison of the relative positions of P 5 [4,6] and P 5 [6,8] observables supports a global deviation with respect to SM predictions over a large q 2 range, rather than an effect localised near the J/ψ resonance that would push up P 5 [6,8] with respect to P 5 [4,6] . Indeed, current data exhibits a pattern opposite to what was proposed in ref. [19] (see the plot for P 5 in figure 12 of ref. [19]). Of course, this cannot be considered as a proof that there are no effects coming from charm resonances, but it supports the concept of a limited impact which does not reach the size advocated in ref. [19].
On the other hand, this bin exhibits a significant discrepancy from SM expectations in P 5 and impacts our analysis. As discussed in section 2.1, we include in our predictions an estimate of the impact of charm resonances, but we also perform cross-checks concerning the role of this bin in section 4.4.

Other observables involved in the fit
Here we discuss a large set of observables that we include in the fit organized in two sets, the first one involving muons and photons in the final state and the second one involving electrons.

b → sµµ and b → sγ observables
This class of observables corresponds to exclusive and inclusive processes where either a real photon or a pair of muons is produced. It includes the decay B → K * µµ discussed at length in the previous sections, but also many other modes of interest.

B s → φµµ
The main difference between this mode and the decay B → K * µµ originates from the fact that B s → φµµ is not self-tagging, i.e. the final state does not contain information on whether the initial meson was a B s or aB s . In the absence of flavour tagging, only a subset of angular observables can be easily measured at a hadron collider, some of them corresponding to CP-averaged angular coefficients (J 1s,1c,2s,2c,3,4,7 ) and some to CP-violating ones (J 5,6s,6c, 8,9 ). Moreover, B s -B s mixing can interfere with direct decay providing additional contributions to the amplitude. This issue was addressed in detail in ref. [66], where it was shown that additional observables could be measured through a time-dependent analysis of the angular coefficients (in particular, promising optimised observables Q 8 and Q 9 ). Furthermore, the measurement of time-integrated angular coefficients in a hadronic JHEP06(2016)092 environment yields O(∆Γ s /Γ s ) corrections to the analogous B + → K * + expressions in terms of transversity amplitudes (related to interference between mixing and decay).
One of the guidelines in our analysis is to try to test the sensitivity of the results on different choices of form factor parametrisations and thus on the specific details and assumptions of a particular form factor computation. Therefore we compare whenever possible the predictions obtained with our default form factor parametrisation to those obtained with other choices, e.g. in the case of B → Kµµ and B → K * µµ results based on KMPW [17] (B-meson LCSR) to results based on BSZ [20] (light-meson LCSR). On the other hand, for the case of B s → φµµ, only two form factor determinations were available at low-q 2 (BZ [39] and BSZ [20]) following rather similar approaches with the latter being an update of the former one.
For this reason and given the importance of this mode, we implemented an alternative approach, based on the B-meson LCSR computation discussed in ref. [67] (corresponding to the same type of method as in KMPW [17]). Unfortunately, ref. [67] does not provide the complete set of form factors necessary for a calculation of the B s → φµµ amplitudes in the full-form factor approach, but the available subset is sufficient to construct the two soft form factors. These are extracted from the full form factors V , A 1 and A 2 in ref. [67] using the value of decay constants, masses and hadronic inputs (we use the same threshold parameter as for K * and the Borel parameter is set to M 2 = 1.0 GeV 2 ). The results obtained for ξ ⊥ and ξ are plotted in figure 6 where they are compared to the corresponding functions from BZ and BSZ. Only central values are shown, illustrating the excellent agreement between the parametrisation using ref. [67] and the BSZ parametrisation up to 5 GeV 2 , and a small deviation (below 8%) in the 5 to 8 GeV 2 region. Considering the very good agreement with the independent computation in ref. [67], we feel confident to use the complete information available for the BSZ parametrisation to implement our soft form factor approach for B s → φµµ.
We thus compute the relevant B s → φµµ observables with the same approach as for B → K * µµ, applied to the form factors from ref. [20] as our default. The O(∆Γ s /Γ s ) corrections to these observables are included using the expressions given in ref. [66], assuming all Wilson coefficients to be real. We use a similar approach for power corrections and du-JHEP06(2016)092 ality violation effects as in the case of B → K * µµ, without assuming any correlation even though SU(3) symmetry is expected to hold approximately. Similarly, for long-distance cc contributions, we use the same estimates for δC cc 9 as in B → K * µµ, 9 but we do not correlate the coefficients s i , a, b, c with those appearing for B → K * µµ.
On the experimental side, LHCb has recently updated the measurement of this mode, providing the branching ratio, its longitudinal fraction F L as well as several CP-averaged angular observables S 3,4,7 which can be recast into optimised observables P 1 , P 4 , P 6 using the correlation matrix provided in ref. [68]. We have checked the linear propagation of errors used to obtain the optimised observables, by converting the B → K * µµ S i observables into optimised P i observables using the same procedure and checking that they agree very well with the values of P i quoted by the LHCb collaboration. Due to differences in the convention in kinematics and the definition of observables, the same dictionary has to be used as in the B → K * case to relate our definitions of angular and optimised observables to those from LHCb articles.

B → Kµµ
In addition to the differential branching ratio, the angular distribution for B → Kµµ features two further observables, the forward-backward asymmetry A F B and the coefficient F H [69]. LHCb data does not suggest any deviation from SM expectations in these two quantities which are sensitive to the presence of scalar/pseudoscalar and tensor operators. Since we do not consider such NP operators, we will only examine the B → Kµµ branching ratio.
The theoretical description of the decay B → Kµµ with the scalar K meson in the final state is considerably simpler than the one of the decay B → K * µµ with the vector K * meson, even though similar conceptual issues are involved. In the large-recoil region, we apply QCD factorisation [41] to the form factors in ref. [17], taking them as uncorrelated. The large-recoil symmetries allows us to reduce the three form factors f + , f 0 , f T to a single soft-form factor. We use the most common scheme [41] where the soft form factor is identified with f + , which dominates the computation of the branching ratio (contributions involving the scalar form factor f 0 are suppressed by the lepton mass, and the tensor form factor f T arises only in the presence of scalar or tensor operators). The dominance of the form factor f + renders correlations with the other two f 0 , f T less important, and therefore the gain of the implementation of correlations via the soft form factor approach is less significant for B → Kµµ than for the vector mode (we checked that using the full form factors for B → Kµµ yields indeed very similar results). Long-distance charm-loop corrections are neglected here, as they are expected to have very little impact on branching ratios [17].
At low recoil, we use the lattice determination from ref. [70]. In this region, the question on the size of duality-violation effects arises as in the case of B → K * . Again we consider a single large bin covering this region and we implement an O(10%) correction (with an arbitrary phase) to the term proportional to C 9 for this bin.

JHEP06(2016)092
There are other important b → s penguin modes sensitive to magnetic and dimuonic operators. We consider the branching ratios B(B → X s γ) Eγ >1.6GeV and B(B → X s µ + µ − ) [1,6] . In both cases, the SM prediction [71][72][73] has gained some recent improvement, with a better control of higher QCD orders for B → X s γ [34][35][36] and the inclusion of logarithmically enhanced electromagnetic corrections for B → X s µ + µ − [37]. This has induced a shift of the SM prediction, both for the central value and the uncertainty. We update the SM predictions entering the relevant formulas for these observables in ref. [3], but we do not modify the part depending on the NP coefficients C i (with NP being constrained to small effects, the inclusion of higher-order effects in this part can be neglected, considering the accuracy aimed at).

B s → µµ
The CMS and LHCb correlations have both measured the branching ratio for B s → µµ, and provided an average of the two measurements [74]. The SM theoretical prediction has been improved significantly over the past year, including NNLO QCD corrections and NLO electroweak corrections, inducing a change in the central value and the uncertainty [75][76][77]. We follow the same approach as for inclusive decays and modify the relevant formulas for these observables in ref. [3] by updating the SM predictions, but without changing the part depending on the NP coefficients C i .
We follow the discussion in ref. [3] for B → K * γ in order to constrain significantly C 7 (and C 7 ). The observables included in our analysis are the isospin asymmetry A I (B → K * γ) and the B → K * γ time-dependent CP asymmetry S K * γ .

Λ b → Λµµ
Another example of a b → sµµ transition is the baryonic mode Λ b → Λµµ, for which the branching ratio and several angular observables have been measured by the LHCb collaboration [78]. Due to limitations of the current theoretical description of this decay (restricted to naive factorisation, with only a limited knowledge of form factors) [79][80][81], we prefer not to include these results in our fits. We note, however, that the current measurement of the differential branching ratio tends to lie below its SM prediction using lattice QCD inputs [78,80].

JHEP06(2016)092
Instead of including directly the ratio R K together with Br(B → Kµµ) in the fit, we use the two branching ratios Br(B → Kµµ) and Br(B → Kee), keeping track of all theoretical correlations among them. Note that in this way we do not lose the information concerning the cancellation of hadronic uncertainties as it would occur in the observable R K because this cancellation is implicitly encoded in the correlations among the two branching ratios. On the experimental side, R K is significantly correlated with Br(B → Kµµ) (in a way not quantified yet), whereas Br(B → Kee) may only have part of the (sub-dominant) systematic uncertainties correlated with Br(B → Kµµ). It seems thus safer to include Br(B → Kee) in the global fit (rather R K ) to avoid a double counting of (correlated) deviations.
Another source of information on b → see is provided by B → K * ee at very low invariant squared masses q 2 of the electron pair, close to the photon pole. An angular analysis [33] provides four observables F L , A T , A re T , A im T (or equivalently F L , P 1,2,3 ), which can be included in the fit to constrain the Wilson coefficients, in particular C 7 and C 7 due to the proximity to the photon pole. 10 Finally, we do not include information on B → X s ee, as this decay provides already little information in the muon case.
In the generic NP models, the effective Hamiltonian involves different effective bs couplings for different lepton species ( = µ, e), so that one should distinguish the Wilson Hence it is not possible to include the above data in a model-independent fit to Wilson coefficients C i , unless an additional hypothesis concerning the value of the C e i or their relationship to the C i is made. Therefore we will not include this set of data in our reference fit described above and in appendix A, but we will consider it in combined µ+e fits, assuming that NP is either absent from C i e , or that it enters flavour-universally in C i e and C i µ .

General framework
We start with a global analysis of the data, in scenarios with potential (real) NP contributions to the Wilson coefficients C 7,9,10,7 ,9 ,10 . 11 We split the SM and NP contributions at 7,9,10 = −0.29, 4.07, −4.31). Our reference fit is obtained using • the observables for b → sγ discussed in section 3.1, • the form factors in ref. [17], apart from B s → φ form factors [20], 10 We will not provide predictions for the branching ratio BR(B → K * ee) at very low recoil: the photon pole magnifies the uncertainty coming from the form factor T1(0), which is very large due to our choice of input for V (0) with large uncertainties. Contrary to the case of angular observables, our estimate for this branching ratio at very large recoil is thus affected by large uncertainties (though in the right ball park of other estimates [21]). 11 We will not consider imaginary contributions to Wilson coefficients and we do not include CP-violating observables in our fits.
For our experimental inputs, we include only LHCb data for the exclusive modes considered here [4,9,24,31,33,68,98], as they dominate the current analysis of the anomalies and allow for a consistent inclusion of correlations. Inclusive modes and b → sγ inputs are taken from the HFAG review [99] and BR(B s → µµ) from the current CMS and LHCb combination [74]. In case of asymmetric error bars, we symmetrise by taking the largest of the two uncertainties quoted, without modifying the central value.
We have to include the experimental and theoretical correlations between the different observables (and bins) for B → K * µµ and B → Kµµ. The experimental correlations are available for B → K * µµ [31], B s → φµµ [68] and B → Kµµ [100]. For B → K * µµ, the correlations are given for both S i and P i observables, whereas they are given only for S i observables for B s → φµµ. We have performed a linear propagation of errors in the latter case in order to obtain the correlations among P i observables (we checked the validity of this procedure by reproducing the correlations among P i observables in B → K * µµ quoted in ref. [31] starting from the information on the S i observables given in the same reference).
For theoretical correlations, we have produced a correlation matrix by performing a propagation of error. This is achieved by varying all input parameters following a Gaussian distribution including known correlations, and determining the resulting distribution of the observables of interest. This is particularly necessary for the form factors: we include correlations between parameters from the lattice QCD computation at low recoil in ref. [101,102]. We treat all parameters as uncorrelated at large recoil in the case of ref. [17], whereas we include the available correlations when we use ref. [20]. We stress that even the uncorrelated scan of parameters (like power corrections) induces correlations among the observables (for instance branching ratios at large recoil) because the latter have a correlated functional dependence on these parameters. The large error bars in ref. [17] for B → K * µµ may lead to excursions in parameter space that distort the distribution of the P i observables and yield significant non-Gaussianities. These non-Gaussianities are avoided by scanning over the input parameters after scaling down all uncertainties by a global factor ρ, producing the correlation matrix for the P i observables, and multiplying all its entries by ρ 2 . The resulting covariance matrix is an accurate representation of the uncertainties and correlations for the P i observables in the vicinity of the central values of the input parameters, as long as it is possible to propagate errors in a linearised way. This matrix encodes all the relevant information concerning uncertainties and correlations among observables, with all uncertainties effectively added in quadrature (we explicitly checked that the results are independent on the exact numerical choice of the rescaling factor ρ, and in practice ρ = 3 is sufficient). The other sets of form factors yield Gaussian distributions for the B s → φµµ and B → Kµµ observables, because of the smaller uncertainty ranges.
Finally, we construct a single covariance matrix as the sum of the experimental (C exp ij ) and the theoretical one (C th ij ), and we use it to build the usual χ 2 function corresponding JHEP06(2016)092 to observables with correlated Gaussian distributions: 12 Once the χ 2 function is computed, it remains to exploit the information that it carries. Following standard frequentist analysis, a first piece of information is provided by the global minimum χ 2 min , which provides an indication of the goodness-of-fit. It can be expressed as a p-value assessing the agreement between the measurements and the scenario tested, given as the probability for a χ 2 -distributed random variable with the corresponding number of degrees of freedom (number of data points minus number of free parameters) to reach a higher value than the one obtained from the data.
If the fit is good enough, one can move on and perform the metrology of the n parameters (NP in Wilson coefficients) by considering the test statistic ∆χ 2 (C i ) ≡ χ 2 (C i ) − χ 2 min . Assuming that this quantity is distributed as a χ 2 random variable with n degrees of freedom, the k-sigma confidence region is obtained as ∆χ 2 (C i ) ≤ ξ(k, n), where ξ(k, n) is the value at which the χ 2 (n)-cumulative distribution function reaches the probability P k σ associated to k sigmas. In practice, ξ(k, 1) = {1, 4, 9}, ξ(k, 2) = {2.3, 6.18, 11.83} and ξ(k, 6) = {5.89, 11.31, 18.21} for k = {1, 2, 3}, corresponding to P kσ = {68.3, 95.4, 99.7}% defined as the probability for a Gaussian random variable to be measured within n standard deviations from the mean. In addition, the pull of the SM is the p-value corresponding to ∆χ 2 (C k = 0), i.e., the probability described above and converted in units of sigma, in the case of a χ 2 (n)-distributed random variable.
When we compare scenarios with different number of parameters, some care is thus needed both for the goodness-of-fit (p-value for χ 2 min ) and the metrology (pull of the SM). For instance, we note that a fit to two parameters (C NP i , C NP j ) may contain the hypothesis C NP i = C NP j = 0 within the 2 σ region, while the corresponding fit to the single parameter C NP i (with C NP j = 0 fixed) might not. In general, p-values and pulls tend to decrease when adding more parameters, unless the added parameters are essential in improving the agreement with data. Having more free parameters in a fit typically reduces the significance of the SM pull and decreases the p-value for χ 2 min if these parameters are not relevant and do not affect the χ 2 function.

One-dimensional fits to Wilson coefficients
First of all, the SM itself does not yield a particularly good fit when considering all the b → sµµ and b → sγ data, with χ 2 min = 110 for N dof = 96, corresponding to a p-value of 12 The theoretical correlation matrices are obtained for the observables in the context of the SM computation. In the following, we will assume that the theory covariance matrix has only a mild dependence on the values of the Wilson coefficients, and we will keep its SM value in the construction of our χ 2 test statistics [15]. We have checked that for the scenarios considered in this paper this assumption holds, by calculating the covariance matrix at the best-fit point and comparing the outcome of the fit with the one using the SM covariance matrix.   16%. We then include NP and start by considering 1D scenarios where only one of the Wilson coefficients is let free to receive NP contributions. The corresponding p-values and pulls for the SM hypothesis gathered in table 2 show clearly that a scenario with NP in C 9 is the most favoured by far. A scenario with NP in C 10 and C 9 is also preferred compared to the pure SM case, but to a lesser extent.

JHEP06(2016)092
It is also interesting to test some scenarios where NP enters in a correlated way in two Wilson coefficients. This occurs in particular in models preserving SU(2) L invariance in the lepton sector [103], or models assuming a vector or axial preference for quark couplings [26][27][28][29]. From table 2, the most favoured scenario corresponds to C NP 9 = −C NP 9 , which could for instance be generated by a Z boson with axial quark-flavour changing and vector muon couplings. This scenario yields a large pull due to the fact that it leads to an excellent agreement with the angular observables at low recoil; however, it has no impact on B → K branching ratios, so that R K remains unexplained. The scenario C NP 9 = −C NP 10 preserving the SU(2) L symmetry can also be considered as interesting. One should however  Table 3. Pulls obtained by allowing successively NP in two Wilson coefficients: for the C j column, the second row gives the pull of the SM hypothesis in the case where C j is let free to vary, whereas the C i row yields the pull of the hypothesis C i = C SM i in the scenario where C i and C j are let free to vary.
be careful not to overinterpret these results: any scenario allowing for NP in C 9 yields a large pull, and the modification of the other Wilson coefficients might slightly improve or worsen the agreement between predictions and measurements, but only with limited impact. We confirm our previous result of 2013 [8] with the 3 fb −1 dataset, namely that C 9 plays a central role in the interpretation of the anomalies, and it is the main Wilson coefficient unavoidably present in any scenario with a pull above 4 sigmas. We find that this Wilson coefficient receives typically a negative contribution of order 25% with respect to the SM. More details on the impact of various experimental inputs and theoretical hypotheses can be found in appendix D.

Two-dimensional fits to Wilson coefficients
It is also interesting to proceed as in ref. [8] and consider nested scenarios where NP is added to one Wilson coefficient after the other, starting from the SM hypothesis. In a given scenario (where some Wilson coefficients C j 1 ,...j N receive NP and the others do not), the improvement obtained by allowing one more Wilson coefficient C i to receive NP contributions can be quantified by computing the pull of the C i = C SM i hypothesis. This allows us to determine the NP scenarios which manage best to reproduce data. From the results in table 3, the most favoured scenarios correspond to (C NP 9 , C NP 9 ) and (C NP 9 , C NP 10 ). This is supported by the actual 2D fits, with results shown in table 4, which also indicates that (C NP 9 , C NP 10 ) is interesting to consider. Other scenarios are also interesting where constraints are used to relate the various NP contributions, for instance C NP 9 = −C NP 9 , C NP 10 = C NP 10 , as well as C NP 9 = −C NP 9 , C NP 10 = −C NP 10 . In figures 7 and 8, we show the 3 σ regions corresponding to the constraints coming from branching ratios and angular observables, and from individual decay channels (respectively) for 4 favoured scenarios. Each constraint is built by considering one of the above subsets and adding the inputs from b → sγ and inclusive decays. Both branching ratios and angular observables favour a negative value of C 9 . As far as channels are concerned, the  discrepancy with the Standard Model is triggered by B → K * µµ and by B s → φµµ (to a lesser extent). Both scenarios with NP in (C 9 , C 9 ) or (C 9 , C 10 ) favour non-zero contributions for both Wilson coefficients, whereas the two scenarios C NP 9 = −C NP 9 , C NP 10 = C NP 10 and C NP 9 = −C NP 9 , C NP 10 = −C NP 10 favour NP in C NP 9 = −C NP 9 mainly (even though contributions to C 10 and C 10 are allowed).
We emphasise that not all those scenarios have an interpretation in terms of a Z which was first proposed by three of us in ref. [8], and was discussed in more detail in refs. [10,15,[26][27][28][29]. Indeed, an interpretation within a Z context would reduce the subset of 2D constrained scenarios to the set of scenarios that fulfills C NP 9 × C NP 10 = C NP 9 × C NP 10 (see appendix F). Notice that this constraint is fulfilled by the scenarios with NP contribution JHEP06(2016)092 Figure 7. For 4 favoured scenarios, we show the 3 σ regions allowed by branching ratios only (dashed green), by angular observables only (long-dashed blue) and by considering both (red, with 1,2,3 σ contours, corresponding to 68.3%, 95.5% and 99.7% confidence levels). Each constraint corresponding to a subset of data includes also the inclusive and b → sγ data. only in C 9 or (C 9 , C 9 ) since both sides of the equation vanish trivially. On the other hand, if one wants to switch on NP in all four coefficients and preserve some simple pattern among them, there are four options that may agree with a Z interpretation: • (C NP 9 = −C NP 9 , C NP 10 = −C NP 10 ), with a large pull for the b → sµµ reference fit, but giving R K = 1 by construction, • (C NP 9 = C NP 10 , C NP 9 = C NP 10 ), disfavoured by the data on B s → µµ, which prefer a SM value for C 10 , leading to a tension with the value of C NP 9 needed for B → K * µµ JHEP06(2016)092 = C NP 9 , C NP 10 = C NP 10 ) which could be interesting candidates but get lower pulls (2.0 and 3.9 σ respectively).
We see therefore that Z scenarios could alleviate part of the discrepancies observed in b → sµµ data, but with only one or two Wilson coefficients receiving NP contributions, corresponding to Z models with definite parity/chirality in its coupling to muons/quarks.
Another important criterion of choice among scenarios comes from considering the main anomalies, namely, P 5 (B → K * µµ), R K and BR(B s → φµµ), and how they are weakened or strengthened in each scenario. As can be seen from appendix A, besides the large JHEP06(2016)092 Impact R K P 5 [4,6], [6,8] B Bs→φµµ  B low recoil   High  I,II  I,VI  VI  III,IV,VI   III  II,IV  II,III,IV,V  I,II IV,VI III,V I V Low V Table 5. Relative impact of each scenario on the anomalies for R K , P 5 , B Bs→φµµ and on the low-recoil bins of the different branching fractions.
deviations of order 2.5 to 3 σ in different observables P 5 , R K and B(B s → φµµ) (that we called generically anomalies), there are also a large set of smaller deviations (many of them at low recoil) that can push in different or similar directions. In appendix B, we illustrate how observables are affected in the presence of NP by providing the predictions and the pulls for the observables at the best-fit point for NP in C 9 only. In table 5 we compare the best fit points for 1D and 2D scenarios involving C 9( ),10( ) , 13 leading to a pull above 4.4 σ: I : C NP 9 , II : (C NP 9 , C NP 10 ), III : (C NP 9 , C NP 10 ), IV : (C NP 9 , C NP 9 ), V : (C NP 9 = −C NP 9 , C NP 10 = −C NP 10 ), VI : (C NP 9 = −C NP 9 , C NP 10 = C NP 10 ) .
We classify these scenarios according to how well they can fix a given anomaly or tension at their best-fit point (reducing it below 1 σ level awards the first position, failing to improve leads to the last position). Some scenarios are unable to improve on certain anomalies: for instance, R K which depends on the combination (C NP 9 + C NP 9 ) − (C NP 10 + C NP 10 ) cannot be explained by a scenario of the type V. In other cases, the observables obtain contributions of opposite signs from the different NP contributions. This is the case for instance for scenarios with C NP 9 = −C NP 9 where a negative C NP 9 goes in the right direction to alleviate the tension in P 5 whereas a positive C NP 9 goes in the wrong direction. However the impact of C NP 9 is only 25% of C SM 9 , so a small positive contribution or the contribution from other coefficients like a positive but small C NP 10 remains a viable possibility to explain the discrepancy in P 5 . The result of this classification is that the scenario V is clearly disfavoured compared to the others which fare almost equally well, with a mild preference for I, II and VI. Only the scenarios II and III improve on the tiny deviation of data with respect to the SM for B s → µµ.
One can compare these results with the recent analysis in ref. [15], which relied on a different approach (full form factor analysis, based on a different set of form factors with correlations [20], use of CP-averaged angular coefficients for the B → V angular analysis). We see that similar 1D scenarios are preferred with a contribution to C 9 alone, C 10 alone to a lesser extent, as well as C NP 9 = −C NP 10 . For 2D scenarios, the best-fit points for (C NP 9 , C NP 9 ) and (C NP 9 , C NP 10 ) are also similar.

Six-dimensional fits to Wilson coefficients
Even though we have seen that the anomalies can be described allowing for additional contributions in two Wilson coefficients, it is interesting to consider the most general scenario 13 We do not consider C NP 9 = −C NP 9 which has a large pull, but is not able to solve the discrepancy in RK . with contributions to all six coefficients. In this case, the best fit point is C (C NP 9 ) are favoured strongly (mildly), whether the other coefficients may vanish at 1 σ but may also accommodate small C NP i in their fairly large confidence intervals. The SM pull is 3.6 σ, which is lower than the pulls for successful 1D and 2D-scenarios, since this scenario allows for more degrees of freedom which are not all relevant to explain the anomalies.

Fits considering lepton flavour (non-) universality
As stated in section 3.2, several measurements have been performed for b → see and can be included in our analysis, as long as we assume some relationship between the Wilson coefficients in the electron and muon sectors C i e and C i µ . In the following, we add to our reference fit the data described in section 3.2, and assume that NP enters b → see and b → sµµ the same way (NP Lepton Flavour Universality [LFU]), that it enters in a different way (NP LFU Violation), or even that there is no NP in the b → see (Maximal NP LFU Violation).
Even in the case of Maximal NP LFU Violation, adding b → see data on the fit may have an impact through the additional constraints that b → see data sets on hadronic inputs (in particular form factors). The main input here is BR(B → Kee), which has a very strong theoretical correlation with BR(B → Kµµ) and thus amounts to including the constraint from R K . Tables 7 and 8 (with b → see data) can be compared with tables 2 and 4 (without it). Since the discrepancy in R K is mainly driven by the disagreement between the SM predictions and the measurements for BR(B → Kµµ), it is not surprising that the 1D scenarios modifying C 9 µ see their significance increase, as well as the p-value associated with the fit (apart from C NP 9 = −C NP 9 which remains unchanged). In particular, scenarios with contribution to C NP 9 only and C NP 9 = −C NP 10 have a large SM pull and a decent p-value. A similar situation occurs for the favoured 2D hypotheses.
It is also interesting to predict R K , R K * and R φ for different scenarios in the intermediate region [1,6] Table 7. Best-fit point, confidence intervals, pulls for the SM hypothesis and p-value for different 1D NP scenarios, including b → see data but assuming NP only in b → sµµ.
the most efficient way to get R K in agreement with the current LHCb values, with values of R K * and R φ around 0.85. Other scenarios yield larger values of R K and smaller for R K * and R φ , apart from the scenario C NP 9 = −C NP 10 which leads to R K , R K * , R φ all around 0.7. The increase in the uncertainties for our predictions for R K * and R φ in NP scenarios comes from the fact that a part of the effects proportional to the lepton mass come from the angular coefficient J 1s which involves 4m 2 /s multiplied by Re( . This term is small in the SM where C 9 −C 10 and thus A R ⊥,|| 0, but in presence of NP not following the same SU(2) L relationship, this contribution increases, with an uncertainty coming mainly from the form factors. We illustrate the sensitivity to the choice of form factors for B → K * where we provide the results using the form factors of ref. [17], compared to ref. [20] (in brackets). The larger uncertainties in the former case come mainly from the normalisation of the form factors. Moreover, one may notice that R K * and R φ are almost identical when using the form factors of ref. [20]: these ratios are driven by the ratios F (0)/V (0) with F = A 1 , A 2 , T 1 , T 2 are almost identical for B → K * and B s → φ in ref. [20].   If NP Lepton Flavour Universality Violation is assumed, NP may enter both b → see and b → sµµ decays though potentially with different values. We show the corresponding constraints in figure 9 for two different scenarios, namely (C NP 9µ , C NP 9e ) and (C NP 9µ = −C NP 10µ , C NP 9e = −C NP 10e ). For each scenario, we see that there is no clear indication of a NP contribution in the electron sector, whereas one has clearly a non-vanishing contribution for the muon sector, with a deviation from the Lepton Flavour Universality line, in global agreement with ref. [15] but with a lower significance.  Table 9. Predictions for R K , R K * , R φ at the best fit point of different scenarios of interest, assuming that NP enters only in the muon sector, and using the inputs of our reference fit, in particular the KMPW form factors in ref. [17] for B → K and B → K * , and ref. [20] for B s → φ.

JHEP06(2016)092
In the case of B → K * , we also indicate in brackets the predictions using the form factors in ref. [20].

Fits to magnetic operators at very low q 2
Traditionally, the main constraints on C 7 , C 7 have been provided by b → sγ observables, both inclusive and exclusive (see e.g. ref. [6]). Recent measurement of b → see observables at very low q 2 provides an alternative source for such constraints, as the photon pole enhances the relative importance of C 7,7 with respect to C 9 ( ) ,10 ( ) . In order to compare the constraining power of both sets of observables separately, and to gauge the impact of including b → s modes in the fit, we have performed separate fits to C NP 7 , C NP 7 in two different scenarios: a) all other Wilson coefficients have their SM value, and b) C NP 9µ = −1.1 and the other coefficients have their SM value (a solution preferred globally by the data, as shown above).
In figure 10 we show the resulting fits. The constraints from b → see observables alone (shown in green) are milder than the b → sγ ones (shown in blue) but the two set of constraints are largely complementary, leading to much tighter constraints once combined (dashed contours). As expected, all these constraints are independent of the value of C NP 9µ . The result of the global fit including all observables (b → sγ, b → see and b → sµµ) is also JHEP06(2016)092 Figure 9. For two scenarios where NP occurs in the two Wilson coefficients C 9µ and C 9e , we show the 1,2,3 σ regions obtained using only BR(B + → K + µ + µ − ) and BR(B + → K + e + e − ) for bins in the [1,6] region (dashed green), and 1,2,3 σ regions using all data from the reference fit and b → see data (solid red). The two NP scenarios correspond to: (C NP 9µ , C NP 9e ) (left) and The diagonal line corresponds to the limit of Lepton Flavour Universality. Same conventions for the constraints as in figure 7. shown (red contours). The constraints are then (slightly) tighter, as b → sµµ observables also constrain magnetic operators, with a clear dependence on C NP 9µ . As C NP 9µ is varied from zero to -1.1 the overall compatibility among the different sets of observables improves.

Role of low-and large-recoil regions in the fit
The issues related to the first and last bins of the large-recoil region were already discussed in section 2.3. One may wonder to which extent our results depend on the inclusion of these bins, in particular the [6][7][8] bin where part of the discrepancies with the SM arises. In section 2.3, we also recalled a different issue, the size of duality-violating effects, Figure 12. For the scenario where NP occurs in the two Wilson coefficients C 7 and C 9 , we compare the situation from the analysis in figure 1 of ref. [8] (on the left) and the current situation (on the right). On the right, we show the 3 σ regions allowed by large-recoil only (dashed green), by bins in the [1][2][3][4][5][6] range (long-dashed blue), by low recoil (dot-dashed purple) and by considering all data (red, with 1,2,3 σ contours). Same conventions for the constraints as in figure 7.
affecting the low-recoil bin. Even though some estimates indicate that they should not affect branching ratios significantly, we are not aware of a similar discussion for angular observables which are an important part of the reference fit. We illustrate the role played by the different bins by considering fits with only the low-recoil region, the large-recoil region, or the bins in the [1,6] GeV 2 range in figure 11. It should be noticed that low recoil favours the same range of NP contributions as the large-recoil bins, but in a milder way. In addition, the [1,6] region provides similar constraints as the whole large-recoil range, implying that our results for the different NP scenarios hold even considering ranges for the dilepton invariant mass where charm contributions are expected to be less relevant. Figure 12 illustrates a similar analysis for the (C 7 , C 9 ) scenario, which updates figure 1 in ref. [8]. There is an overall similarity, with a best-fit point requiring almost no NP contributions to C 7 . We stress that the right-hand plot involves a larger set of experimental measurements and a more complete understanding of the sources of theoretical uncertainties on the right. In addition, "only [1,6] bins" refers to observables in the single bin [1,6] only on the 2013 plot (on the left), but to those taken in any of the (smaller) bins inside the [1,6] range on the 2015 plot (on the right).

Tests of SM theoretical uncertainties
The previous studies show the robustness of the results when only part of the experimental information is included in the fit. On the other hand, since the main discrepancies in the previous fits come from exclusive b → sµµ transitions (B → K * µµ, B s → φµµ and B → Kµµ), one ought to consider the sources of systematics entering the SM the-JHEP06(2016)092 Figure 13. For 4 favoured scenarios, we show the 3 σ regions allowed by S i angular observables for B → K * µµ and B s → φµµ only (dashed green), by P i angular observables for B → K * µµ and B s → φµµ only (long-dashed blue), and by considering all data with P i angular observables (red, with 1,2,3 σ contours). Same conventions for the constraints as in figure 7. oretical predictions carefully, namely: form factor uncertainties, power corrections and long-distance corrections due to cc loops. We will consider these different sources of uncertainties in the following.

Role of the form factors
Predictions for B → K * µµ observables depend on seven hadronic form factors whose calculation via non-perturbative methods like light-cone sum rules (LCSR) suffers from relatively large uncertainties (typically ∼ 20 − 50%). It is thus natural to raise the question whether an underestimation of the form factor uncertainties could be the origin of the JHEP06(2016)092 Figure 14. For 4 favoured scenarios, we show the 3 σ regions allowed using form factors in ref. [20] in the full form factor approach (long-dashed blue) compared to our reference fit with the soft form factor approach (red, with 1,2,3 σ contours). Same conventions for the constraints as in figure 7.
observed anomaly [21]. There are two different issues to be distinguished, namely on one hand the overall size of the form factor uncertainties, and on the other hand the correlations among the errors of the different form factors.

Overall size of uncertainties
Let us first stress that the overall size of the form factor uncertainties has a minor impact on global fits, and in the case of clean observables P ( ) i even on the predictions for individual observables. The reason is that assuming a precise knowledge of the correlations among the form factors, they cancel at leading order in the construction of the observables P   1,2,3 σ contours). We also show the 3 σ region obtained from the fit to power-correction-insensitive observables (mostly low recoil), which would correspond to the limiting fit with completely arbitrary power corrections. Same conventions for the constraints as in figure 7.

JHEP06(2016)092
reducing the impact of their errors to a next-to-leading-order effect O(α s , Λ/m B ). For the observables S i this effect only occurs in a global fit where the correlation between different observables effectively reduces the sensitivity to the form factors, while individual S i observables display a form factor dependence at leading order. Note that the size of the form factor errors entering our analysis is much more conservative than what is typically assumed in other analyses [15,22,38], as we are taking form factors from ref. [17] where particularly large errors are assigned. In ref. [22] the error of the normalisation of the soft form factor ξ ⊥ (0) = 0.31 ± 0.04 is determined by considering the spread of the central values of various different non-perturbative form factor calculations like light-cone sum rules [17,39] and Dyson-Schwinger equations [104]. This has to be compared with our value ξ ⊥ (0) = 0.31 +0.20 −0.10 that has an error band exceeding by far the one in ref. [22], implying that it covers the form factor values that would be obtained by the other methods [39,104].
We performed various tests on the sensitivity of our results to the choice of form factors. First, we checked the dependence on the choice of form factors for the observables that are most sensitive to the form factors, namely the branching ratios, in the Standard Model case. To this end we have compared our prediction for BR(B → K * µµ) using the B-meson LCSR determination (KMPW [17]) with other predictions available in the literature based on a different form factor determination (BSZ [20]). We found a good agreement at the 1 σ level for the different bins we compared, while for the total BR(B → K * µµ) the agreement is stronger (below the 1 σ level). In the case of B → Kµµ we observe a systematic difference in the branching ratio at the order of 30% compared to ref. [15], which entirely stems from the difference between the set of form factors chosen (KMPW versus BSZ) and illustrates the sensitivity of these observables to the set of form factors considered.
To demonstrate the limited role of the size of the form factor uncertainties in a global analysis, one can trade the optimised angular observables P i for the CP-averaged angular observables S i [47] which are known to be more sensitive to form factor inputs [7]. The comparison presented in figure 13 shows that the outcome of the fit is very similar in both cases, which is owed to the correlations among the seven form factors restored via the approximate large-recoil symmetries (see below) and reducing the sensitivity to the overall size of uncertainties. We observe a systematic albeit small increase in significance of around 0.3 σ when P i observables are used compared to S i observables.

Correlations
The correlations among the different form factors can in principle be extracted from the corresponding calculation, as it has been done for example in ref. [20]. On the other hand, the dominant correlations can also be assessed from first principles relying on symmetry relations fulfilled by the form factors at low q 2 . While this second approach is more general and avoids any dependence on the details of a particular non-perturbative calculation, it provides the correlations only up to symmetry-breaking corrections of the order Λ/m b (factorisable power corrections). In our analysis we explicitly introduce these symmetry-breaking corrections by hand and assign to them an error of the order of 10% of the respective full form factor, corresponding to a 100% error of the factorisable power corrections. We could confirm that this assumption of 10% power corrections is a realistic JHEP06(2016)092 estimate by determining the central values for the power corrections from a fit to a particular non-perturbative calculation: it has been done in ref. [23] for the B → K * form factors in refs. [39] and [17] and in appendix E for B s → φ and B → K .
We illustrate the compatibility of the two approaches at the level of the global fit analysis in figure 14. We compare the results of our reference fit, performed applying the soft-form factor approach based on the large-recoil symmetries described in section 2.1 and using mainly the B-meson LCSR results of ref. [17], with the full-form factor approach applied to the light-meson LCSR results of ref. [20] (including correlations, similarly to ref. [10,15]). We see that both sets of results are very similar, even though in the soft-form factor approach we started from a set of form factors with larger uncertainties and no knowledge of correlations. This highlights the advantages of the soft-form factor approach to restore correlations among form factors. Not surprisingly, the full-form factor approach based on ref. [20] is more constraining than our soft-form factor approach based on the results of ref. [17], which exhibits much larger uncertainties for the form factor parameters.
For the reasons mentioned above, our SM predictions as well as our fit results are in good agreement with ref. [15]. It is thus surprising that the authors of ref. [22] find much larger errors from factorisable power corrections. This necessarily implies that they must implicitly have introduced a much stronger breaking of the large-recoil symmetry relations, in contradiction to expectations from dimensional arguments as well as to the explicit results for the particular LCSR calculation [20]. In other words, the results in ref. [22] taken at face value should imply that the recent LCSR estimates performed in ref. [20] are not correct.
One may wonder how big the large-recoil symmetry breaking effects (i.e. the factorisable power corrections) should be in order to produce a similar pattern of deviations as observed in the data. In order to study this, we performed the test of assuming twice or four times larger power corrections (corresponding to 20% or 40% of the corresponding full form factors). The results in figure 15 show that the factorisable power corrections only play a minor role in the uncertainties and the outcome for our reference fit. When power corrections are increased from 10% to 40%, the fit is more and more driven by observables with no sensitivity to power corrections, such as low-recoil observables. Indeed, one can see that the shape of the 3 σ regions in figure 15 evolves into the low-recoil regions shown in figure 11 as the size of power corrections increases, which are reproduced in figure 15 to ease the comparison.
If one wants to solve the anomalies exhibited in b → sµµ processes through power corrections, it is important not to focus on one single observable, like P 5 , alone but on the full set. Since power corrections enter many observables, trying to adjust them to fix one observable may generate a problem in another one. The effect of correlations is illustrated in figure 16, inspired by figure 5 of ref. [22]. For comparison we work here in the soft-form factor scheme employed in ref. [22] with the soft form factors defined from the full form factors T 1 and A 0 . We also switch to the helicity basis used in ref. [22] where for example JHEP06(2016)092 Figure 16. Power corrections a V − and a V + needed to obtain agreement between SM predictions and experiment at 1 σ, considering different observables. This plot illustrates that a V ± can indeed be used to obtain agreement between SM prediction and experiment in one observable, but correlations hinder a similar agreement when a larger set of observables is considered. the helicity form factors V ± are defined in terms of the transversity form factors V and A 1 as For the constant terms a V + and a V − in the series of power corrections for the form factors V ± this implies Figure 16 shows that power corrections can explain the data for P 5 [4,6] within the 1 σ range if they occur at the level of 20% for both a V + and a V − in the combination In a scheme where the soft form factor ξ ⊥ is defined from the full form factor V [23], such power corrections are absorbed into ξ ⊥ . Figure 16 displays also the power corrections needed for the observables P 5 [6,8] , P 2 [4,6] and P 1 [4,6] . Comparing the region for P 5 [4,6] and P 1 [4,6] one notices that a solution for P 5 [4,6] through large power corrections moves P 1 [4,6] away from the measured value. An explanation of all three observables within the SM in terms of power corrections requires to reach the limit of the 20% region. Only marginal agreement is obtained then, and once P 5 [6,8] is added to the list, no common solution is found even for power corrections beyond 20%.
This situation seems to be in contradiction with figure 5 of ref. [22]. Note, however, that the (a V + , a V − ) profile shown there corresponds to a scenario where all other power JHEP06(2016)092 correction parameters have been fixed in such a way to describe best the experimental data, without specifying their presumably quite large values. In fact, already the point a V + = a V − = 0 is in agreement with data nearly at the 1 σ level, even though the power correction parameters a V + and a V − are the most relevant ones for the observable P 5 . It is further irritating that, while it is claimed that power correction parameters are scanned only in a range of ±10% of the soft-form factor value, for the plot a region covering |a V ± | ≤ 0.2 has been chosen, corresponding to power corrections of order ±66%.

Role of long-distance charm corrections
Another frequent attempt to explain the B → K * µµ anomaly consists in assuming a very large charm-loop contribution (or non-factorizable power correction). It is not difficult to imagine that with a sufficiently general q 2 -dependent parametrisation one might easily fit any data [63]. However, one must check that the parametrisation itself and the resulting fit respect all known properties of the related charm-loop correlator, as well as its behaviour at large recoil. In the end, the assumption that the charm contribution is responsible for the anomalies leads to two kinds of predictions: first, those arising from correlations that might survive among observables under the most general parametrisation of the correlator (which would give little information on C 9 ), and second, the prediction that R K = 1 due to SM lepton-flavour universality. Whatever explanation might be first assumed concerning b → sµµ transitions, one has still to invoke a NP contribution to explain R LHCb K 0.75, most plausibly in the form of a non-SM contribution to C 9µ . But once such NP contribution has been introduced, the other b → s anomalies are reduced and there is no actual need for abnormally large non-perturbative hadronic effects on top of the NP contribution.
We would like to stress that explaining some of the anomalies through such large charm contributions leads to further difficulties. First of all, these explanations predict an enhanced effect when one gets closer to the J/ψ peak. They typically lead to an increase for the observable P 5 in the [6,8] region with respect to [4,6] (see figure12 of ref. [19] and 'prediction column' of table 6 of ref. [63] for P 5 ). But current data does not seem to follow this behaviour: an increase in statistics in this particular region will be very important to settle this point definitely. Another important issue comes from the comparison between lowand large-recoil regions: the charm effects advocated in refs. [19,63] to explain the current data at large recoil within the SM do not provide any clue about the deviations observed in B → K * µµ and B s → φµµ branching ratios at low recoil, whereas a single short-distance contribution to C 9 is able to accommodate the deviations in both regions simultaneously.
A confirmation of the deviation measured in R K with a higher significance, as well as the measurement of other observables exhibiting lepton-flavour universality violation (see e.g. ref. [105]) would strongly disfavour solutions involving non-perturbative charm-loop effects such as the ones proposed in refs. [19,63]. Conversely, a clear evidence for a q 2dependent effect, or the need for different contributions in different transversity amplitudes in B → K * µµ, would lead to prefer non-perturbative QCD effects rather than New Physics. However, there is no evidence for such a dependence on q 2 or transversity in the present data. A further discussion of this issue can also be found in ref. [106].

Increasing the size of the charm contributions
Long-distance charm corrections have been subject to many recent discussions, with different estimates [17][18][19]. We recalled in section 2.1 that we use the work of ref. [17] as an estimate of this effect to be added on top of the perturbative contribution, but without assuming a specific sign for this contribution. In our reference fit, for each transversity amplitude of B → K * µµ and B s → φµµ we multiply this contribution by s i = 0 ± 1 (hence six uncorrelated parameters). We present in figure 17 the corresponding results if we take contributions twice or four times larger. Increasing the size of the charm contributions reduces the significance of the deviations from Standard Model, but the discrepancy remains above 3 σ for the various scenarios considered even if the long-distance cc contribution is multiplied by 4 compared to our reference fit.

Distinguishing new physics from charm contribution in C 9
Another way of checking the robustness of our approach with respect to charm consists in determining if the fit to the data favours an additional q 2 -dependent contribution to C 9 . In that case, this would be a clear indication that some long-distance contribution has been underestimated in our analysis, as NP contributions cannot have any such dependence.
We have performed fits to the same data as in the reference fit, but limited to particular q 2 -ranges, in order to check the stability of the value of C 9 needed in different bins. We can perform this fit under different hypotheses: for instance, one can leave only C NP 9 , or assume that C NP 9 = −C NP 9 , or that C NP 9 = −C NP 10 . An underestimated hadronic contribution from charm loop would correspond to a q 2 -dependent contribution to C 9 only, i.e., the first case. In the two other cases, the need for a q 2 -dependent contribution might indicate a problem of consistency in the fit that could not be understood only through a hadronic contribution. Figure 18 shows no need for a q 2 -dependent contribution in these three situations. 14 In the case of the last scenario C NP 9 = −C NP 10 (bottom plot in figure 18), the fit tries to accommodate both B s → µµ (constraining C NP 10 to remain small) and B → K ( * ) µµ and B s → φµµ observables (favouring a significant contribution to C NP 9 ). The fit exhibits thus a tension between the two types of constraints. In order to assess this effect, we have performed a fit without B s → µµ. The result, indicated with dotted lines, favours lower values of C NP 9 = −C NP 10 , within a narrower range, without spoiling the good agreement between the global and bin-by-bin analyses.
As an alternative test, we added three q 2 -dependent contributions to C SM 9 of the form C had 9,p (s) = A p + B p × s for p = K, K * , φ. We assumed that the same contribution entered the three transversity amplitudes identically for B → K * µµ (we assumed the same in the case of B s → φµµ). A 6D fit to the real parameters A K,K * ,φ and B K,K * ,φ in the large-recoil region showed a clear preference for A K * and A φ negative and different from zero, a mild preference for A K negative and different from zero, whereas B K,K * ,φ remained unconstrained, confirming that the current fit needs a negative contribution to C 9 in order to explain the data, but that it does not exhibit a preference for a q 2 -dependent contribution.

JHEP06(2016)092
theoretical improvements concerning various sources of uncertainties (form factors, power corrections, charm contribution), we have updated and considerably extended the analysis of ref. [8] performed by three of us in 2013, using the recently published LHCb results based on the full 3 fb −1 dataset from LHC run I.
We confirm the previous result of ref. [8], namely that C 9 plays a central role in explaining the anomaly: a negative NP contribution to this Wilson coefficient (typically of order 25% with respect to the SM value) is unavoidably present in any scenario with a pull above 4 σ. Other coefficients play a secondary role but might lead to an increase in the significance. In this sense we found several scenarios with one or two free parameters that exhibit a pull of more than 4 σ compared to the SM hypothesis. Oneparameter scenarios with that property are C NP 9 , C NP 9 = −C NP 10 and C NP 9 = −C NP 9 , and two-parameter scenarios are (C NP 7 , C NP 9 ), (C NP 9 , C NP 10 ), (C NP 9 , C NP 7 ), (C NP 9 , C NP 9 ), (C NP 9 , C NP 10 ) and (C NP 9 = −C NP 9 , C NP 10 = C NP 10 ), (C NP 9 = −C NP 9 , C NP 10 = −C NP 10 ). We have performed a global fit to all six Wilson coefficients simultaneously, and found that all the coefficients are consistent with their SM values at the level of 1-2 σ, except for C 9 , in line with the results from the more economical scenarios mentioned above.
We have also briefly discussed the situation in the context of models violating leptonflavour universality, by allowing NP contributions of different sizes in the electron and muon sectors. While the data requires a NP contribution in the muon sector to explain the anomalies, it does not show preferences for a contribution in the electron sector (and thus more generally disfavours a lepton-flavour universal NP contribution). If one restricts NP to the muon sector, some of the above scenarios see their significance increase, with a SM pull very close (or equal) to 5 σ in some instances.
We have performed several checks to test the robustness of our results. We have compared different possible choices for the analysis: QCD factorisation with soft-form factors versus the computation with full form factors, different choices for the set of LCSR form factors taken as input, optimised observables P i versus CP-averages S i , different choices for the binning. We find a very good agreement between the various approaches. In particular, we have checked that the details of the form factor computation are not very significant for the optimised observables.
Non-perturbative effects from power corrections and from long-distance cc contributions can only be estimated. We have studied the effect of increasing the size of these contributions, without finding a large impact on the overall picture presented above. Moreover, the above-mentioned hadronic effects (in particular the cc contributions) are expected to exhibit a q 2 -dependence which allows them to be distinguished from a q 2 -independent NP effect. We have studied a possible q 2 -dependence in a twofold way: on one hand, performing a bin-by-bin analysis, on the other hand introducing a separate linear q 2 -dependence in C 9 for the fit to B → Kµµ, B → K * µµ and B s → φµµ. In both cases, we found no conclusive evidence for a q 2 -dependence.
One should notice that our results are in good agreement with those obtained in ref. [15], even though the applied methods differ in many central points: different sets of form factors, a different approach to the computation (soft form factors versus full form factors), different angular observables and different estimates of hadronic uncertainties JHEP06(2016)092 R K P 5 [4,6], [6,8]   (power corrections, charm contribution). While our method is to a large extent independent of the modelling of non-perturbative effects but has to rely on an estimation of subleading contributions based on dimensional arguments, the analysis in ref. [15] is based on (and limited to) a particular non-perturbative LCSR calculation. Strengths and weaknesses of the two approaches are of complementary nature, and the comparison of the obtained results is thus a useful cross-check of the different hypotheses on which the two analyses rely.
While the observables S i become competitive to the P i in a global fit, where their LO form factor dependence gets cured thanks to correlations, the P i exhibit a much larger sensitivity to NP on the level of the individual observables as they are shielded to a large extent from hadronic uncertainties. Whereas for example the observable P 5 can be predicted in the SM with a precision of ∼ 10%, basically independently of the underlying form factor parametrisation, predictions for S 5 can develop uncertainties up to ∼ 40% depending on the form factors used as input. This feature makes the experimental measurement of the observables P i indispensable in the search for NP where it will be essential to find apart from global tensions in combined fits also some clear-cut discrepancies in individual observables.
The results we obtained from our fits are particularly encouraging as they show that at the level of the Wilson coefficients several NP scenarios provide a consistent explanation of the deviations observed in b → s transitions. On the other hand, the most favoured scenarios are difficult to generate in terms of simple NP models (such as a heavy Z boson, or leptoquarks). Obviously this might change over time with new experimental results. In this respect, we found it interesting to summarise in table 10 how a given NP contribution to a Wilson coefficient would affect the different anomalies. As expected, only C NP 9 < 0 is able to provide a consistent explanation for all of them. There is also certain preference for C NP 10 and C NP 10 to be positive in order to explain two out of three anomalies, and negative for C NP 9 and C NP 9 . However, whereas the best-fit point of the 1D and 2D (unconstrained)

JHEP06(2016)092
scenarios with NP in C 9 and C 10 (see tables 2 and 4) indeed shows a preference for a negative (respectively, positive) value in agreement with table 10, the best-fit point for C NP 9 and C NP 10 prefers a positive (respectively, negative) value in contradiction with table 10. This suggests that the chirally-flipped operators are not particularly useful to solve the anomalies but they are quite efficient (especially C 9 ) in fixing small deviations in various bins, summing up to an overall large significance. Such a situation arises in the scenario C NP 9 = −C NP 9 which fixes neither R K nor the anomaly in P 5 , but still manages to yield a very large pull with respect to the SM hypothesis.
In summary, C NP 9 < 0 is very much favoured, providing a consistent picture for the anomalies in agreement with the results of global fits. A contribution C NP 10 > 0 comes in second place, while the situation with respect to NP contributions to the chirally-flipped operators is less clear. Obviously, this guess work is completely tributary to the current experimental situation. Updates of these measurements, and in particular B → K * µµ observables with a finer binning, will prove particularly important to provide a more definite answer concerning the origin the anomalies observed in b → s transitions. channels (Marseille), where part of this work was discussed. We thank Tobias Huber for sharing with us the results of ref. [37] prior to publication. We thank Roman Zwicky for bringing ref. [58] to our attention while we were completing this work. JV is funded by the DFG within research unit FOR 1873 (QFET), and acknowledges financial support from CNRS. SDG, JM and JV acknowledge financial support from FPA2014-61478-EXP. L.H. has been supported by FPA2011-25948 and the grant 2014 SGR 1450, and in part by the Centro de Excelencia Severo Ochoa.

A SM predictions
The prediction column corresponds to the Standard Model case.  = −C NP 10 ) (lower right), we show the 3 σ regions allowed by branching ratios only (dashed green), by angular observables only (long-dashed blue) and by considering both (red, with 1,2,3 σ contours). Same conventions for the constraints as in figure 7.

C Confidence regions for selected 2D new physics scenarios
In figure 19, we provide the confidence regions of interest for two-dimensional scenarios less favoured from the point of view of the fit, but which might be of interest for model building, namely contributions to (C NP 9 , C NP 10 ), (C NP 7 , C NP 9 ), (C NP 9 = −C NP 10 , C NP 9 = C NP 10 ) and (C NP 9 = −C NP 10 , C NP 9 = −C NP 10 ). F (q 2 ) = F soft (q 2 ) + ∆F αs (q 2 ) + ∆F Λ (q 2 ). (E.1)

JHEP06(2016)092
The soft component F soft is a linear combination of two soft form factors ξ ⊥ and ξ for M vector, and proportional to a single soft form factor ξ P for M pseudoscalar. The decomposition eq. (E.1) is not unique: depending on the exact definition of the soft form factors ξ i , a part of ∆F αs and ∆F Λ can be reabsorbed into the soft contribution F soft . This introduces a scheme dependence for ∆F αs and ∆F Λ which has been discussed in detail for the B → K * form factors in ref. [23]. While QCD corrections ∆F αs can be calculated employing QCD factorisation, the power corrections ∆F Λ cannot be computed directly and in general, they must be estimated on dimensional grounds. However, one can perform explicit computations of the full form factors F (q 2 ) (say, from light-cone sum rules) in order to extract ∆F Λ (q 2 ) through a fit. In the case of the B → K * form factors, this determination has been performed in ref. [23]. The parametersâ F ,b F ,ĉ F arising in the parametrisation can be found in tables 1 and 2 of that paper for two different choices of scheme for ξ ⊥ and ξ (considering either the LCSR calculation in ref. [17] or that in ref. [39] as inputs). In table 13, we give the corresponding results for B s → φ and B → K form factors using LCSR input from refs. [20] and [17], respectively. We follow scheme 1 in ref. [23] and define the soft form factors as for B s → φ, and as ξ P (q 2 ) = f + (q 2 ), (E.4) for B → K. We further quantify the relative size of power corrections for the various form factors though the ratio at q 2 = 0, 4, 8 GeV 2 . From dimensional arguments one expects r(q 2 ) = O(Λ/m B ) 10%. The results in table 13 show that the LCSR form factors from refs. [20] and [17] indeed comply with this expectation, except for the form factor A 2 where larger power corrections occur.

JHEP06(2016)092
Bs → φâFbFĉF r(0 GeV 2 ) r(4 GeV 2 ) r(8  Table 13. Fit results for the power-correction parameters to the B s → φ and B → K form factors, choosing a scheme with the soft form factors (ξ ⊥ , ξ ) defined from V and the difference of A 1 and A 2 in the case of B s → φ, and with ξ P defined from f + in the case of B → K. The corresponding LCSR input has been taken from ref. [20] for B s → φ and from ref. [17] for B → K. Furthermore, the relative size r(q 2 ) with which the power corrections contribute to the full form factors is shown for q 2 = 0, 4, 8 GeV 2 .
In our SM predictions as well as in the NP fits, we use the results from table 13 as central values for the parameters a F , b F , c F , to which we assign error ranges of the order of 10% × F . Comparing with r(q 2 ), we see that this corresponds to the assumption of O(100%) uncertainties for the coefficients a F , b F , c F . Since our error estimate is based only on dimensional arguments, it is independent of the detail of the particular LCSR calculation. On the other hand, taking into account correlations among the LCSR form factors, it is also possible to determine the uncertainties of a F , b F , c F from a particular set of LCSR input, which will be detailed in an upcoming publication [106].

JHEP06(2016)092
These couplings obey the relationship C NP 9 × C NP 10 = C NP 10 × C NP 9 . (F.5) A Z model can therefore belong to the following categories: • NP only in the following pairs, with a priori arbitrary contributions, (C 9 , C 10 ) , (C 9 , C 9 ) , (C 10 , C 10 ) , (C 9 , C 10 ) , (F. 6) each case corresponding to the vanishing of some of the couplings ∆ sb L,R , ∆ µµ V,A . These models have a definite chirality for quark-flavour changing coupling currents and/or a definite parity for the couplings to muons.
• NP enters all four semileptonic coefficients with the following relationships Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.