On flavourful Easter eggs for New Physics hunger and lepton flavour universality violation

Within the standard approach of effective field theory of weak interactions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta B = 1$$\end{document}ΔB=1 transitions, we look for possibly unexpected subtle New Physics effects, here dubbed “flavourful Easter eggs”. We perform a Bayesian global fit using the publicly available HEPfit package, taking into account state-of-the-art experimental information concerning these processes, including the suggestive measurements from LHCb of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{K}$$\end{document}RK and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{K^{*}}$$\end{document}RK∗, the latter available only very recently. We parametrise New Physics contributions to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b \rightarrow s$$\end{document}b→s transitions in terms of shifts of Wilson coefficients of the electromagnetic dipole and semileptonic operators, assuming CP-conserving effects, but allowing in general for violation of lepton flavour universality. We show how optimistic/conservative hadronic estimates can impact quantitatively the size of New Physics extracted from the fit. With a conservative approach to hadronic uncertainties we find nonzero New Physics contributions to Wilson coefficients at the level of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 3\sigma $$\end{document}∼3σ, depending on the model chosen. Furthermore, given the interplay between hadronic contributions and New Physics effects in the leptonic vector current, a scenario with nonstandard leptonic axial currents is comparable to the more widely advocated one with New Physics in the leptonic vector current.


Introduction
Easter eggs nowadays also refer to inside jokes and/or secret messages usually hidden e.g. in computer gaming and hi-tech software. In this work, we take advantage of this terminology to motivate the search for New Physics Beyond the Standard Model in the radiative and in the (semi)leptonic channels of rare B meson decays.
In the decades that have followed the original formulation of flavour mixing [1], the flavour structure of the SM has been experimentally tested and well established. The tremendous progress of the experimental facilities has probed the flavour of the SM to an exquisite level of precision [2], along with the substantial effort on the part of the theoretical community to go well beyond leading order computations [3]. From this perspective of "precision tests", radiative and (semi)leptonic ΔB = 1 processes, related at the partonic level to b → sγ, s transitions, occupy a special place in probing the SM and its possible extensions in terms of New Physics (NP) models [4,5].
Firstly, these rare B meson decays belong to the class of flavour-changing neutral current (FCNC) processes, which are well known to be sensitive probes of Physics Beyond the Standard Model (BSM): in fact -within the SM -the flavour structure of the theory allows FCNC to arise only at loop level, as a consequence of the GIM mechanism [6]. This allows for significant room for heavy new degrees of freedom to sizeably contribute to these rare processes.
Secondly, from the experimental side, the study of rare B meson decays offers us some of the most precise measurements amongst the |ΔF| = 1 processes. For instance, the measurement of the inclusive branching fraction of B → X s γ is currently performed with a relative uncertainty of a few percent [7][8][9], while the study of an exclusive mode such as B → K * allows for a detailed analysis of the angular distribution of the four final state particles, yielding rich experimental information in terms of angular functions of the dilepton invariant mass, with full kinematic coverage of the latter [10] and -starting from Ref. [11] -also with available experimental correlations among the angular observables.
In B Physics, the recent years have been characterised by the emergence of a striking pattern of anomalies in multiple independent studies of some of these rare b → s transitions [12]. Of particular importance, the measurement of the P 5 angular observable [13][14][15][16] stands out from all the other ones related to the angular distribution of B → K * μμ ; first realised by the LHCb collaboration [17,18] and later on also by the Belle collaboration [19], the experimental analysis of P 5 in the large recoil region of the decay points to a deviation of about 3σ with respect to the SM prediction presented in Ref. [20]. The latter, however, suffers from possible hadronic uncertainties which are sometimes even hard to guesstimate [21][22][23][24], and this observation has been at the origin of a quite vivid debate in the recent literature about the size of (possibly) known and (yet) unknown QCD power corrections to the amplitude of this process in the infinite mass limit [25][26][27][28]. To corroborate even more the cumbersome picture of the "P 5 anomaly", two new independent measurements of this angular observable (among others) have been recently released by the ATLAS [29] and CMS [30] collaborations, showing, respectively, an appreciable increase and reduction of the tension between data and the SM prediction in Ref. [20], as reported by these experiments.
For the sake of completeness, one should also remark that other smaller tensions have been around, concerning the measurement of differential branching fractions of B → K μμ [31,32] and B s → φμμ [33]. It is worth noting that, while for the latter mode an explanation in terms of hadronic physics may easily be conceivable, the theoretical computation of the former seems to be under control [34].
Quite surprisingly, a possible smoking gun for NP in rare B meson decays already came out in 2014, when the LHCb collaboration presented for the first time the measurement of the ratio of branching fractions [35]: where the subscript refers to the dilepton mass (denoted hereafter q 2 ) range going from 1 to 6 GeV 2 . This experimental value shows a deviation of about 2.6σ with respect to the standard theoretical prediction. Indeed, the SM value of R K in the bin provided by the LHCb collaboration is expected to be equal to unity beyond the percent level of accuracy [36,37]. In fact, contrary to observables such as P 5 , it is important to stress that R K may be, in general, regarded as insensitive to QCD effects [36]. From the model building point of view, R K can certainly be considered as quite informative, hinting at a UV completion of the SM where lepton flavour universality violation (LFUV) takes place in the flavour-violating couplings of new heavy degrees of freedom, e.g. leptoquarks and/or Z gauge bosons . Most importantly, the tantalising correlation of this signature of LFUV with the P 5 anomaly, suggested by several global analyses [4,[68][69][70][71][72] has triggered different proposals of measurements of such effect in the angular analysis of the K * channel [73,74]. Interestingly enough, an analysis from the Belle collaboration aiming at separating the leptonic flavours in B → K * [75] shows a consistent ∼ 2.6σ deviation from the SM prediction reported in Ref. [20] in the dimuon leptonic final state only. This is compatible with previous experimental findings related only to the mode with muonic final states.
Sitting on similar theoretical grounds to R K , another intriguing ratio of B decay branching fractions can be measured in the K * channel: These measurements for the low-q 2 bin and the central-q 2 one have just been presented by the LHCb collaboration [76], pointing again to a discrepancy of about 2σ with respect to the expected SM prediction -again equal to 1 to a very good accuracy for the central-q 2 bin and close to 0.9 for the low-q 2 one -and yielding more than a 3σ deviation when naively combined with the measurement of R K . Note that with higher degree of braveness (or, depending on the taste of the reader, of unconsciousness), the disagreement of the SM with precision B physics may reach the exciting level of 5σ when one naively combines together the single significances coming from R K ,K * ratios, P 5 measurements and the minor deviations observed in the other exclusive branching fractions.
Given the excitement of these days for all the above hints of a possible NP discovery in rare B meson decays, in this work we take our first steps towards a positive attitude in the search of a definite BSM pattern aimed at addressing these B anomalies. We perform our study in a model-independent fashion, within the framework of effective field theories for weak interactions [77][78][79]. In particular, in Sect. 2 we define the setup characterising the whole global analysis, presenting six different benchmark scenarios for NP, together with a discussion as regards two different approaches in the estimate of the hadronic uncertainties that can affect quantitatively our final results. In Sect. 3, we list all the experimental measurements we use to construct the likelihood in our fit, and we discuss in detail our most important findings. The latter are effectively depicted in Figs. 1, 2, 3, 4, 5, 6, and collected in Tables 2, 3, 4, 5 in Appendix A. In Sect. 4 we summarise our conclusions.

Theoretical framework of the analysis
In this section we present the effective field theory framework at the basis of this work and introduce the benchmark scenarios we focus on for our study of NP effects in rare B decays. We then illustrate the two distinct broad classes of assumptions that characterise our global analysis: the case where we take an optimistic attitude towards the estimate of hadronic uncertainty plaguing the amplitude of both B → K * /γ and B s → φ /γ channels, and a second one where we aim at providing a more conservative approach. All the results in Sect. 3.2 will be classified under these two different setups.

New Physics benchmarks for ΔB = 1
Integrating out the heavy degrees of freedom, the resulting effective Hamiltonian of weak interactions for b → sγ, s transitions involves the following set of dimension six operators within the SM [80]: where = e, μ, p = u, c and we have neglected the chirally suppressed SM dipoles. The ΔB = 1 effective Hamiltonian can be cast in full generality in the form of a combination of two distinct parts: where, within the SM, the hadronic term involves the first seven operators in Eq. (5): while the second piece includes the electromagnetic dipole and semileptonic operators: with λ i corresponding to the CKM combination V ib V * is for i = u, c, t and where C i=1,...,10 are the Wilson coefficients (WCs) encoding the short-distance physics of the theory. All the SM WCs in this work are evolved from the mass scale of the W boson down to μ b = 4.8 GeV, using state-of-the-art perturbative QCD and QED calculations for the matching conditions [81][82][83] and the anomalous dimension matrices [83][84][85][86] relevant for the processes considered in this analysis.
While a general UV completion of the SM may enter in the effective couplings present in both pieces of Eq. (5), general NP effects in b → sγ, s can be phenomenologically parametrised as shifts of the Wilson coefficients of the electromagnetic and semileptonic operators at the typical scale of the processes, μ b . In particular, the most general basis for NP effects in radiative and (semi)leptonic B decays can be enlarged by the presence of scalar, pseudo-scalar and tensorial semileptonic operators, together with right-handed quark currents as the analogue of Q 7γ , Q 9V , Q 10 A SM operators [21,87]. In this work, motivated by previous interesting findings concerning LFUV [69][70][71] and the measurement of R K and R K * , we focus on the contributions of NP appearing as shifts of the SM WCs related to the electromagnetic dipole and semileptonic operators with left-handed quark currents only. A comprehensive analysis with different chiral structures as well as a more general effective theory framework will be presented elsewhere [88]. Furthermore, we restrict ourselves to CP-conserving effects, taking NP WCs to be real.
For NP in semileptonic operators we discriminate between couplings to muon and electron fields both in the axial and vector leptonic currents. We characterise our phenomenological analysis for NP through six different benchmark scenarios, studying the impact of combinations of the following NP WCs: (I) C NP 9,μ and C (VI) C NP 7 , C NP 9,μ , C NP 9,e , C NP 10,μ and C NP 10,e varied simultaneously in the respective ranges defined above, i.e. a NP scenario described by five different parameters.
We remark that while benchmarks (I) and (II) have been already studied in the literature, none of the other cases has been analysed so far. In particular, NP scenarios (III) and (IV) allow us to study, for the first time, the interesting impact of a NP radiative dipole operator in combination with vector-like and axial-like LFUV effects generated by NP. Most interestingly, scenario (V) allows us to explore the correlation C NP 9 = −C NP 10 , possibly hinting at a SU (2) L preserving BSM theory. As an additional interesting case to explore, we eventually generalise to simultaneously nonvanishing C NP 7 , C NP 9,μ , C NP 9,e , C NP 10,μ and C NP 10,e in case (VI). We wish to stress that all of the six benchmarks defined above will be studied for the first time under two different approaches in the estimate of QCD hadronic power corrections, as presented in next section.

Treatment of the hadronic uncertainties
In our previous work [24,27,89], we went into considerable detail on the treatment of hadronic contributions in the angular analysis of B → K * . Our approach there was to study how large these contributions can be assuming that the LHCb data on branching fractions and angular distributions of these decay modes could be described within the SM. For that purpose we considered four scenarios for the hadronic contributions, with increasing theoretical input from the phenomenological analysis presented in Ref. [90]. The underlying functional form that we used for the hadronic contribution was given by where we fitted for the complex, helicity dependent, coefficients h (i) λ , (i = 0, 1, 2) and (λ = 0, +, −) using the data and the phenomenological model in [90]. Since h 0 enters the decay amplitude with an additional factor of q 2 with respect to h ± , we drop h (2) 0 in our analysis. In this work we proceed to study the possible existence of NP contributions in semileptonic and radiative b → s decays which requires a re-evaluation of the hadronic uncertainties. For the sake of simplicity, to address hadronic contributions we use the same functional parameterisation as given in Eq. (8). However, we limit ourselves to only two hadronic models. The first, corresponding to the most widely used assumption, relies completely on the phenomenological model in [90] below q 2 < 4m 2 c . The second is a more conservative approach, where we impose the latter only in the large recoil region at q 2 ≤ 1 GeV 2 while letting the data drive the hadronic contributions in the higher invariant mass region. We will refer to the first approach as phenomenological model driven (PMD) and the second as phenomenologically and data driven (PDD). In our fit we vary the h i λ parameters over generous ranges. More detailed discussion of these can be found in [24,27].
In the present analysis we also need to address modes that were not considered in our previous work, namely B → K , B s → φ and B s → φγ . The decay B → K has been studied in detail in [34], where the authors show that the hadronic uncertainties are smaller than in B → K * . A comparison of the LCSR estimate of the soft gluon contribution and the QCDF estimate of the hard gluon contribution reveals that the soft gluon exchange is subdominant with respect to QCDF hard gluon exchange. Therefore, although in principle the same concerns on the soft gluon contribution we raised for B → K * apply in this case, in practice the overall effect of soft gluons can be reasonably neglected. In our computation we therefore only include hard gluon exchange computed using the QCDF formalism in Ref. [91].
The long distance contributions for B s → φ and B s → φγ follow a similar theoretical derivation to those for B → K * and B → K * γ , respectively, barring the fact that the spectator quark in the former is different from that in the latter. No theoretical estimates of power corrections to the infinite mass limit are available for the B s → φ /γ decays and one has to rely on the ones for the B → K * /γ decays to get a handle on the long distance contributions. The spectator quark effects can come through the hard spectator scattering involving matrix elements of Q 2 , P 6 and Q 8g computable in QCD factorisation [91] which we include in our computation. However, we do not include the subleading, and numerically small, QCDF power corrections to spectator scattering involving Q 8g [92][93][94] and contributions to weak spectator scattering involving Q 8g beyond QCDF computed in LCSR [95][96][97]. The effect of the difference in all these spectator contributions is expected to be low firstly because they are numerically small and, secondly, because the effect is proportional to the small flavour SU (3) breaking. Different approaches in relating the long distance contributions in the B → K * /γ channels to the ones in the B → φ /γ channels have been used in the literature [69,70,98], which vary in the degree of correlation between the two. While Ref. [70] uses uncorrelated hadronic uncertainties, Refs. [69,98] have left the two contributions highly correlated noting that the spectator contribution is expected to be numerically small. We take an approach similar to the latter considering the insensitivity of the current data to such effects and use the same value of power corrections in B → K * and B s → φ amplitudes, even though this choice pertains to a quite oversimplifying optimistic attitude. We leave a more detailed analysis of this assumption by relaxing the correlation between the hadronic contributions in the two modes to a future work [88].

Experimental information considered
In this section we discuss the experimental measurements we use in our fit. Note that for the exclusive modes we make use of measurements in the large recoil region only. Our choice harbours on the fact that the QCD long distance effects in the low recoil region are substantially different from the large recoil regime [99][100][101][102] and would require a dedicated analysis. For the fit in this study we consider the following experimental information: For the B → K * μμ channel we use the LHCb measurements of CP-averaged angular observables extracted by means of the unbinned maximum likelihood fit, along with the provided correlation matrix [18]. Moreover, we employ the recent results for CP-averaged angular observables from ATLAS [29] and the ones measured by CMS [30,103]. 1 Finally, we use the CP-averaged optimised angular observables recently measured by Belle [75] 2 . Regarding the differential branching fractions, we use the recently updated measurements from LHCb [104] and the ones from CMS [103]. For the B → K * ee channel we consider the LHCb results from [105] and the Belle results from [75]. R K * observable is considered according to the recently presented measurements by LHCb [76] in both the low-q 2 and the central-q 2 bins; see also Eq. (3). 1 For all CMS data we use the 7, 8 TeV combined results, which can be found in https://twiki.cern.ch/twiki/bin/view/CMSPublic/ PhysicsResultsBPH13010 . 2 Belle measures the B 0 → K * 0 μμ and B + → K * + μμ channels together, without providing the mixture ratio. On the theoretical side, we can therefore use these measurements under the approximation that QCD power corrections differentiating the amplitudes of the two channels are small. We have numerically checked that the impact of known QCD power corrections [91] is indeed at the percent level in the observables of interest.
Our theoretical predictions are computed in the helicity basis, whose relevant expressions can be found in [21]; the same framework is employed to study B → K * γ , B s → φμμ, B s → φγ and B → K channels. For the latter, we use the full set of form factors extrapolated from the lattice results, along with the provided correlation matrix [106]; for the remaining channels, we use the full set of form factors estimated combining LCSR and lattice results, along with the correlation matrices [107].
For the factorisable and non-factorisable QCD power corrections, we refer to Sect. 2.2.
We include in our analysis the HFAG average for the branching fractions from [2].
We consider the LHCb CP-averaged angular observables and differential branching fractions measurements, along with the provided correlation matrix [33].
We use the LHCb measurement of the branching fraction from [108].
We employ the LHCb measurement of B → K ee differential branching fraction and R K from [35].
We use the HFAG average from [2]. We perform our theoretical computation at NNLO in α s and NLO in α em , following Ref. [

Results of the global fit
In this section we present the main results of our work. We perform this study using HEPfit [113] relying on its Markov Chain Monte Carlo-based Bayesian analysis framework implemented with BAT [114]. We fit to the data using 16 real free parameters that characterise the non-factorisable power corrections, as was done in [24], along with the necessary set of NP WCs. We assign to the hadronic parameters and the NP WCs flatly distributed priors in the relevant ranges mentioned in Sect. 2. The remaining parameters used in the fit are listed in Table 1. To better compare different scenarios, we use the Information Criterion [115,116], defined as where log L is the average of the log-likelihood and σ 2 log L is its variance. The second term in Eq. (9) takes into account the effective number of parameters in the model, allowing for a meaningful comparison of models with different number of parameters. Preferred models are expected to give smaller I C values.
The results for NP WCs for the several cases that we study can be found in Figs. 1, 2, 3, 4, 5 and 6, where the I C value for each model is also reported, and in Tables 2 and 3 in Appendix A. In Tables 4 and 5, we report the results of the fit for observables of interest. We observe that all cases have comparable I C values except cases (IV) and (V), which are disfavoured in the PMD approach while they remain viable in the PDD one. The main difference between the two approaches is that angular observables, in particular P 5 , call for NP in C NP 9,μ in the PMD approach, while they can be accommodated within the SM in the PDD one. Fig. 1 The two NP parameter fit using C NP 9,μ and C NP 9,e . Here and in the following, the left green panel shows the results for the PMD approach and the right red panel shows that for the PDD one. In the 1D distributions we show the 16th, 50th and 84th percentile marked with the dashed lines. In the correlation plots we show the 1, 2 and 3σ contours in decreasing degrees of transparency. The blue square and lines identify the values of the NP WCs in the SM limit. The numbers at the bottom left corner of the 2D plots refer to the correlation. We also report I C values for the two approaches (see Eq. (9)). Preferred models are expected to give smaller I C values Fig. 2 The two NP parameter fit using C NP 9,μ and C NP 10,μ . See caption of Fig. 1 for the colour coding and further details Let us discuss the various cases in more detail. It is important to stress that the evidence of NP in our fit for cases (I)-(V) is always larger than 3σ for one of the semileptonic NP WCs used in the analysis, given the need of a source of LFUV primarily from R K ,K * measurements. In particular, we remark that in the PMD scenarios of cases (I) and (II) we get evidence for NP at more than 5σ . However, looking at the corresponding PDD scenarios, the NP evidence gets significantly reduced, roughly between 3σ and 4σ . The reduction in the Fig. 3 The three NP parameter fit using C NP 7 , C NP 9,μ and C NP 9,e . See caption of Fig. 1 for the colour coding and further details significance comes from the larger hadronic uncertainties in the PDD approach which weaken the constraining power of the angular observables on the NP WCs.
Concerning case (III), we observe very similar findings to the ones obtained for case (I), since the effective coupling for the radiative dipole operator is well constrained, especially from the inclusive B → X s γ branching fraction.
Regarding case (IV), in which we vary the three NP parameters C NP 7 , C NP 10,μ and C NP 10,e , the model comparison between Fig. 4 The three NP parameter fit using C NP 7 , C NP 10,μ and C NP 10,e . See caption of Fig. 1 for the colour coding and further details the PDD and PMD realisation of this NP benchmark is quite informative: NP effects in the dipole operator and in the axial semileptonic currents cannot address at the same time R K ,K * ratios and the P 5 anomaly in a satisfactory way when we stick to small non-factorisable QCD power corrections; however, this is no longer true when we allow for a more conservative estimate of the hadronic uncertainties. In particular, the tension in the fit coming from the angular analysis of B → K * μμ can now be addressed by large QCD effects as Fig. 5 The three NP parameter fit using C NP 7 , C NP 9,μ , C NP 9,e and C NP 10,μ,e = −C NP 9,μ,e . See caption of Fig. 1 for the colour coding and further details those given in Eq. (8), while a C NP 10,e = 0 at about 3σ can successfully describe all the observational hints of LFUV showed by current measurements. This interesting possibility of axial lepton-flavour violating NP is not found in other global analyses [69][70][71][72], as it proceeds from the conservative treatment of hadronic uncertainties we proposed in Ref. [24].
Concerning Tables 4 and 5 of Appendix A, we would like to point out the pattern displayed by the transverse ratios R T K * and R T φ : cases (I)-(III) predict these values to be ∼ 1 Fig. 6 The five NP parameter fit using C NP 7 , C NP 9,μ , C NP 9,e , C NP 10,μ and C NP 10,e . See caption of Fig. 1 for the colour coding and further details with a small error, while the remaining cases give different predictions with the central value ranging between ∼ 0.7 and ∼ 0.8. Therefore, obtaining experimental information on transverse ratios may help in discerning between the different NP scenarios. We then show results for case (V), in which we vary C NP 7 , C NP 9,μ , C NP 9,e and correlate the semileptonic vector and axial currents according to C NP 9,μ = −C NP 10,μ and C NP 9,e = −C NP 10,e . In analogy with case (IV), only within the PDD approach we find for this NP benchmark a fairly good description of data, with C NP 9,μ = −C NP 10,μ compatible with zero at ∼ 2σ . Again, we are presented with the case where deviations in angular observables are addressed by large QCD power corrections, while LFUV is driven by semielectronic operators. Looking back at Tables 4 and 5, we note that, for this case, as well as for case (IV) and (VI), both transverse and longitudinal muon over electron ratios in the central-q 2 bin, namely R T K ,K * ,φ and R L K ,K * ,φ , are characterised by similar central values. We close our presentation with an analysis of case (VI) in which we float simultaneously C NP 7 , C NP 9,μ , C NP 9,e , C NP 10,μ , and C NP 10,e . As can be seen from Fig. 6, current measurements are informative enough to constrain, at the same time, all the NP WCs both in the PMD and PDD approaches. In particular, within the latter case, a nontrivial interplay among NP effects encoded both in C NP 9,μ and C NP 10,e , together with the hadronic contributions reported in Table 3, produces the weakest hint in favour of NP provided by our global analysis -sitting between 2σ and 3σ level -while allowing for a very good description of the entire data set, similar to the other cases. The power corrections we found are larger than those obtained in Ref. [90], but smaller than those required by the SM fit of B → K * μμ [24]. As discussed in detail in Refs. [27,89], the size obtained for the power corrections is compatible with the naive power counting relative to the leading amplitude. We stress (once again) that a more optimistic attitude towards the estimate of QCD power corrections (PMD approach) leads to the a much stronger claim in favour of NP, at a statistical significance larger than 5σ .
In Tables 2 and 3 we report mean and standard deviation for the NP WCs and absolute values of h λ for all the cases considered in the analysis. It is also relevant to observe that, once we switch on NP effects through C NP 9,μ in order to attempt at simultaneously explaining observables such as R K ,K * and P 5 in the PDD approach we find values for |h (1,2) λ | compatible with zero at ∼ 1σ . Conversely, if we set C NP 9,μ = 0 then a nonvanishing |h (2) − | is needed to account for the angular observables, as found in Ref. [24], showing that one cannot disentangle hadronic uncertainties and NP in B → K * μμ at present.

Discussion
In this work, we critically examined several BSM scenarios in order to possibly explain the growing pattern of B anomalies, recently enriched by the R K * measurement performed by the LHCb collaboration [76]. We carried out our analysis in an effective field theory framework, describing the non-factorisable power corrections by means of 16 free parameters in our fit along the lines of Ref. [24].
We performed all our fits using two different hadronic models. The first approach, labelled PMD, relies completely on the phenomenological model from Ref. [90] and corresponds to the more widely used choice in the literature. The second one, named PDD, imposes the result of Ref. [90] only at q 2 ≤ 1, 3 allowing the data to drive the hadronic contributions in the higher invariant mass region.
Regarding the NP contributions, we analyse six different benchmark scenarios, differentiated by distinct choices of NP WCs employed in the fits. Case (I) allows for C NP 9,μ and C NP 9,e , while case (II) considers the scenario with C NP 9,μ and C NP 10,μ ; case (III) studies NP effects coming as C NP 7 , C NP 9,μ and C NP 9,e , and case (IV) is the same as the latter but with C NP The comparison of different scenarios using the I C shows that all the considered cases are on the same footing except for cases (IV) and (V). These cases are strongly disfavoured in the PMD approach, as there is no C NP 9,μ in case (IV) to account for the deviation in P 5 , while C NP 9,μ is constrained by its correlation with C NP 10,μ and the measured value of BR(B s → μμ) in case (V).
In fact, from our analysis of radiative and (semi)leptonic B decays we identify two classes of viable NP scenarios: -The widely studied C NP 9,μ = 0 scenario: from Figs. 1, 2 and 3, we find a remarkable 5σ evidence in favour of C NP 9,μ = 0 in the PMD approach. It is indeed nontrivial that a single NP WC can explain all the present anomalies in b → s transitions [4,[68][69][70][71][72]. However, in the more conservative PDD approach, the significance of a nonvanishing C NP 9,μ drops to about 3σ , mainly driven by LFUV.
-An alternative scenario with nonvanishing C NP 10,e , which emerges in the presence of large hadronic corrections to 3 This choice is motivated in Ref. [24]. the infinite mass limit, namely our PDD approach. To our knowledge, a NP electronic axial current has not been studied in the literature, since it does not provide a satisfactory description of the angular observables within the commonly used PMD approach. We think that the present theoretical status of power correction calculations is not robust enough to discard this interesting NP scenario.
Finally the most general fit we performed, namely case (VI), confirms in the PDD approach that both scenarios above are viable, although a slight preference for C NP 9,μ = 0 is found. More data are needed to assess what kind of NP scenario (if the anomalies persist) is realised in Nature.

Appendix A: Numerical results
In this appendix we present the tables with the most relevant numerical results obtained from our global analysis. Mean and standard deviation for the NP WCs and h λ absolute values are reported in Table 2 for the PMD approach and in Table 3 for the PDD one. 4 In Table 4 we list the results in the PMD approach obtained for the key observables in the six NP scenarios. Analogous results for the PDD approach can be found in Table 5.