B-decay discrepancies after Moriond 2019

Following the updated measurement of the lepton flavour universality (LFU) ratio R_K in B ->Kll decays by LHCb, as well as a number of further measurements, e.g. R_K* by Belle and B_s ->mu mu by ATLAS, we analyse the global status of new physics in b ->s transitions in the weak effective theory at the b-quark scale, in the Standard Model effective theory at the electroweak scale, and in simplified models of new physics. We find that the data continues to strongly prefer a solution with new physics in semi-leptonic Wilson coefficients. A purely muonic contribution to the combination C_9 = -C_10, well suited to UV-complete interpretations, is now favoured with respect to a muonic contribution to C_9 only. An even better fit is obtained by allowing an additional LFU shift in C_9. Such a shift can be renormalization-group induced from four-fermion operators above the electroweak scale, in particular from semi-tauonic operators, able to account for the potential discrepancies in b ->c transitions. This scenario is naturally realized in the simplified U_1 leptoquark model. We also analyse simplified models where a LFU effect in b ->sll is induced radiatively from four-quark operators and show that such a setup is on the brink of exclusion by LHC di-jet resonance searches.


Introduction
In recent years, several deviations from standard model (SM) expectations have been building up in B-decay measurements. While each of them could be a first sign of physics a e-mail: jason.aebischer@tum.de b e-mail: waltmann@ucsc.edu c e-mail: diego.guadagnoli@lapth.cnrs.fr (corresponding author) d e-mail: meril.reboud@lapth.cnrs.fr e e-mail: peter.stangl@lapth.cnrs.fr f e-mail: david.straub@tum.de beyond the SM, statistical fluctuations or underestimated experimental or theoretical systematic uncertainties cannot be excluded at present. These deviations -or "anomalies"can be grouped into four categories that have very different experimental and theoretical challenges: (i) Apparent suppression of various branching ratios of exclusive decays based on the b → sμμ flavourchanging neutral current (FCNC) transition [1,2]. The uncertainties are dominated by the limited knowledge of the B to light meson hadronic form factors [3][4][5]. (ii) Deviations from SM expectations in B → K * μ + μ − angular observables [6][7][8][9] (also based on the b → sμμ transition), where form factor uncertainties are much smaller than for the branching ratios, but hadronic uncertainties are nevertheless significant [10,11]. (iii) Apparent deviations from μ-e universality in b → s transitions in the processes B → K and B → K * (via the μ/e ratios R K [12] and R K * [13], respectively).
(Note that e-μ universality in b → c ν transitions is tested to hold at the percent level [26][27][28].) While the deviations in (i) and (ii) could be alleviated by more conservative assumptions on the hadronic uncertainties, it is tantalizing that a simple description in terms of a single nonstandard Wilson coefficient of a semi-muonic operator like (sγ ρ P L b)(μγ ρ μ) or (sγ ρ P L b)(μγ ρ P L μ) leads to a consistent description of (i), (ii), and (iii), with a best-fit point that improves the fit to the data by more than five standard deviations compared to the SM (for a single degree of freedom) [29][30][31][32][33][34]. Moreover, it was shown that simplified models with a single tree-level mediator can not only explain (i), (ii), and (iii), but even all four categories of deviations simultaneously without violating any other existing constraints [35][36][37][38][39][40]. Taken together, these observations explain the buzz of activity around these deviations and the anticipation of improved measurements of the theoretically clean ratios R K and R K * . The purpose of this article is to examine the status of the tensions after inclusion of a number of updated or newly available measurements, in particular: • The new measurement of R K by the LHCb collaboration combining Run-1 data with 2 fb −1 of Run-2 data (corresponding to about one third of the full Run-2 data set).
The updated measurement finds [ where the first uncertainty is statistical and the second systematic, and q 2 is the dilepton invariant mass squared. The SM predicts lepton flavour universality, i.e. R SM K is unity with uncertainties that are well below the current experimental sensitivities. While the updated experimental value is closer to the SM prediction than the Run-1 result [12], the reduced experimental uncertainties imply a tension between theory and experiment at the level of 2.5σ , which is comparable to the situation before the update.
• The new, preliminary measurement of R K * by Belle [42].
(2) 1 In our numerical analysis, we use the full one-dimensional numerical likelihood provided in [41], which is markedly non-Gaussian, rather than symmetrizing the uncertainties in (1). 2 Belle provides also results for the q 2 bins [0.045, 1.1] GeV 2 and [1.1, 6] GeV 2 . The data from these bins is contained in the q 2 bin [0. 1,8] GeV 2 used in our analysis, except for the region between 0.045 GeV 2 and 0.1 GeV 2 , which, however, is affected by relatively large theoretical uncertainties due to its proximity to the di-muon kinematical threshold [14]. We also include the q 2 bin [15,19] GeV 2 , which is integrated over a sufficiently large q 2 range to make theoretical predictions in this kinematic region valid. For a detailed discussion of the operator product expansion and violation of quark hadron duality at high q 2 , see e.g. [43][44][45].
Given their sizable uncertainties, these values are compatible with both the SM predictions (R SM K * approximately unity) and previous results on R K * from LHCb [13] R K * = BR(B → K * μμ) BR(B → K * ee) = 0.66 +0.11 −0.07 ± 0.03, for 0.045 GeV 2 < q 2 < 1.1 GeV 2 , 0.69 +0. 11 −0.07 ± 0.05, for 1.1 GeV 2 < q 2 < 6 GeV 2 , that are in tension with the SM predictions by ∼ 2.5σ in both q 2 bins. • One further, important piece of information included in our study is the 2018 measurement of B s → μμ by the ATLAS collaboration [46], that we combine with the existing measurements by CMS and LHCb [47][48][49]. • Finally, we discuss possible connections to the anomalies in the b → c ν decays. We take into account the latest updates of R D and R D * from Belle [50,51] which, compared to previous results, have central values closer to the SM predictions. Taking into account these updates, the latest world averages from HFLAV are [52] R D = 0.340 ± 0.027 ± 0.013, corresponding to a combined deviation from the SM by 3.1σ .
In this paper we will explore the implications of all these, as well as other data, to be described in fuller detail in the next section, in the context of global fits to effective-theory new-physics scenarios, whether at the weak scale or at higher scales, identify those scenarios that lead to a good description of the data, and discuss possible realizations in terms of simplified new-physics models. Our numerical analysis is entirely based on open-source software, notably the global likelihood in Wilson coefficient space provided by the smelli package [53], built on flavio [54] and wilson [55]. As such, our analysis is easily reproducible and modifiable.
The rest of this work is organized as follows.
• In Sect. 2, we describe our statistical approach and the experimental measurements we employ in our numerical analysis. • In Sect. 3, we perform a global analysis of b → s transitions, first in the weak effective theory (WET) below the electroweak (EW) scale, then in the SM effective field theory (SMEFT) above the EW scale, which allows us to extend the discussion to the charged-current deviations and to incorporate constraints from electroweak precision tests and other precision measurements. • In Sect. 4, we discuss a number of specific simplified new-physics (NP) models that are favoured by the current data, assuming the deviations to be due to NP. • Sect. 5 contains our conclusions.

Setup
Our numerical analysis is based on a global likelihood function in the space of the Wilson coefficients of the WET valid below the EW scale, or the SMEFT valid above it. Theoretical uncertainties (for observables where they cannot be neglected) are treated by computing a covariance matrix of theoretical uncertainties within the SM and combining it with the experimental uncertainties (approximated as symmetrized Gaussians, using the same procedure as in previous studies by some of us [29,30,[56][57][58][59]). The main assumption in this approach is that the sizes of theory uncertainties are weakly dependent on NP, which we checked for the observables included. This approach was first applied to b → s transitions in [59]. The theoretical uncertainties in exclusive B-decay observables stem mainly from hadronic form factors, which we take from [3] for B to light vector meson transitions and from [5] for B → K , as well as unknown nonfactorizable effects that are parametrized as in [3,54,59] (and are compatible with more sophisticated approaches [10,11]). Additional parametric uncertainties (e.g. from CKM matrix elements) are based on flavio v1.3 with default settings [54]. For more details on the statistical approach and the list of observables and measurements included, we refer the reader to [53].
Here we highlight the changes in observables sensitive to b → s transitions included with respect to the recent global analyses [29,30] by some of us.
• We include the LHCb update of R K [41] (cf. footnote 1) and the new, preliminary measurement of R K * by Belle [42]. The Belle results are available for various q 2 -bin choices, separately for B ± and B 0 decays and in an isospin averaged form. In our numerical analysis we use the 0.1 GeV 2 < q 2 < 8 GeV 2 and 15 GeV 2 < q 2 < 19 GeV 2 bins, separately for B ± and B 0 decays.
, that we combine with the CMS and LHCb measurements [47][48][49]. This combination is discussed in detail in Appendix A. Our combination is in slight tension with the SM prediction of BR(B s → μ + μ − ) by approximately 2σ . • We include the updated LHCb measurement of forwardbackward asymmetries in b → + − [60] as well as its branching ratio [61]. For the theory predictions of the baryonic decay we follow [62,63]. • Here we are working with the global likelihood described in [53] (i.e. including as many observables sensitive to the Wilson coefficients as possible), while in [29] we focused on observables sensitive to the b → s transition only. This means e.g. that we also include all the observables sensitive to the b → sγ, g dipole transitions studied in [64]. In addition, the global likelihood also includes observables that do not directly depend on the Wilson coefficients of interest but whose theory uncertainties are strongly correlated with those of the directly dependent observables. This is in particular relevant for the b → sμμ observables. In our figures, we indicate the set of observables consisting of b → sμμ, b → sγ, g, and other correlated observables as "b → sμμ & corr. obs.". Like in the previous analysis [30] by some of us, we again include the LFU differences of angular observables 3 D P 4 and D P 5 . To this end, we have added them to the global likelihood in version 1.3.0 of the smelli package.

Effective-theory analysis
Having at hand the global likelihood in the space of NP Wilson coefficients, L( C), we perform a numerical analysis within effective field theories, the WET and the SMEFT, by discussing in detail certain, well-motivated, one-and twocoefficient scenarios (see e.g. [66]), and collecting a fuller set of them in Appendix E. This analysis proceeds in two steps: 1. We first investigate the Wilson coefficients of the WET at the b-quark mass scale. This analysis can be seen as an update of earlier analyses (see e.g. [29][30][31][32][33][34]) where our considered scenarios have already been discussed. Such analysis is completely general within our assumed oneor two-coefficient assumption, and barring new particles lighter than the b quark (see e.g. [67][68][69][70]). 2. Next, we embed these results into the SMEFT at a scale above the electroweak scale. This is based on the additional assumptions that there are no new particles beneath and that EW symmetry breaking is approximately linear (see e.g. [71]). This allows us to correlate NP effects in b → s , within a general, common formalism, with other sectors like EW precision tests or b → c transitions (cf. [53,[72][73][74]).

b → s observables in the WET
We start by investigating the constraints on NP contributions to the | B| = | S| = 1 Wilson coefficients of the WET at the b-quark scale μ b ≈ m b that we take to be 4.8 GeV. We work with the effective Hamiltonian with the normalization factor The dipole operators are given by 4 where σ μν = i 2 [γ μ , γ ν ], and the semi-leptonic operators We have omitted from H bs eff, NP semi-leptonic tensor operators, which are not generated at dimension 6 in theories that have SMEFT as EW-scale limit, as well as chromomagnetic and four-quark operators. Even though the latter can contribute via one-loop matrix elements to b → s processes, their dominant effects typically stem from renormalization group (RG) evolution above the scale μ b , and we will discuss these effects in the SMEFT framework in the next section. For the same reason, we have constrained the sum over lepton flavours to e and μ: semi-tauonic WET operators can contribute via QED RG mixing, but their direct matrix elements are subleading [75]. 4 The sign of the dipole coefficients C ( ) 7 are fixed by our convention for the covariant derivative D μ ψ = ∂ μ + ieQ ψ A μ + ig s T A G A μ .

Scenarios with a single Wilson coefficient
We now consider the global likelihood in the space of the above Wilson coefficients. We start with scenarios where only a single NP Wilson coefficient (or a single linear combination motivated by UV scenarios) is nonzero. The best-fit values, 1 and 2σ ranges, and pulls for several such scenarios are listed in Table 1. For the 1D scenarios, the pull in σ is defined as We make the following observations.
• Like in previous analyses, two scenarios stand out, namely a shift to C bsμμ 9 by approximately −25% of its SM value (C SM 9 (μ b ) 4.1), or a shift to the combination C bsμμ 9 = −C bsμμ 10 by approximately −15% of its SM value. However, at variance with previous analyses, it is the second scenario, rather than the first one, to have the largest pull. Given our assumptions about hadronic uncertainties, the pull exceeds 6σ . The main reason why the combination C bsμμ 9 = −C bsμμ 10 performs better is the ∼ 2σ tension in BR(B s → μμ), which remains unresolved in the C bsμμ 9 scenario. We will comment further on this finding in Appendix C.
• New physics in C bsμμ 10 alone also improves the agreement between theory and data considerably. However, tensions in B → K * μμ angular observables remain in this scenario.
• Muonic scenarios with right-handed currents on the quark side, C bsμμ 9 and C bsμμ 10 , or the lepton side, C bsμμ 9 = C bsμμ 10 , do not lead to a good description of the data.
• Scenarios with NP in bsee Wilson coefficients only, while able to accommodate the discrepancies in R K ( * ) , do not help for the rest of the data. Given the pulls in Table 1, such scenarios are, all in all, less convincing.
The scalar Wilson coefficients C bsμμ S/P and C bsμμ S/P are strongly constrained by the B s → μμ decay. In theories that have SMEFT as their EW-scale limit, they satisfy the constraint [77] In this case, they cannot lead to significant modifications in semi-leptonic b → sμμ transitions [76]. However, the preference of the combination discussed in Appendix A for a suppressed B s → μ + μ − branching ratio means that a destructive interference of these Wilson coefficients with the SM contribution to the leptonic decay can lead to a moderate improvement of the likelihood.

Scenarios with a pair of Wilson coefficients
Next, we consider the likelihood in the space of pairs of Wilson coefficients. The results in Table 1 suggest that NP in both C bsμμ 9 and C bsμμ 10 ought to give an excellent fit to the data. The left plot of Fig. 1 shows the best fit regions in the C bsμμ 9 -C bsμμ 10 plane. The orange regions correspond to the 1σ constraints from b → sμμ observables (including B s → μ + μ − ) and observables whose uncertainties are correlated with those of the b → sμμ observables (cf. last point in Sect. 2). In blue we show regions corresponding to the 1σ (right plot) and 2σ (left plot) constraints from the neutral-current LFU (NCLFU) observables R K , R K * , D P 4 , and D P 5 . In the right plot, the 1σ constraints from only R K (purple) and only R K * (pink) are shown. The combined 1 and 2σ region is shown in red. The dotted contours indicate the situation without the Moriond-2019 results for R K and R K * . The best fit point C 0.40 has a χ 2 = 6.6, which, corrected for the two degrees of freedom, corresponds to a pull of 6.3σ . In this scenario a slight tension between R K and R K * remains, as it predicts R K R K * while the data seems to indicate R K > R K * . In addition, there is also a slight tension between the fit to NCLFU observables and the fit to b → sμμ ones, especially in the C bsμμ 9 direction.
Overall, we find a similarly good fit of the data in a scenario with NP in C bsμμ 9 and C bsμμ 9 . The scenario is shown in the right plot of Fig. 1 0.47. The χ 2 = 6.4 corresponds to a pull of 6.0σ . Interestingly, in this scenario a non-zero C bsμμ 9 is preferred at the 2σ level. The right-handed quark current allows one to accommodate the current experimental results for the LFU ratios, R K > R K * . This scenario cannot address the ten- Other two-coefficient scenarios (including dipole coefficients, scalar coefficients, and electron specific semileptonic coefficients) are discussed in Appendix E.

Universal vs. non-universal Wilson coefficients
In view of the updated R K ( * ) measurements, which are closer to the SM prediction than the Run-1 results, our fit in C bsμμ 9 and C bsμμ 10 shows a tension between the fit to NCLFU observables and the fit to b → sμμ ones, especially in the C bsμμ 9 direction. Therefore, it is interesting to investigate whether lepton flavour universal new physics that mostly affects b → sμμ observables but none of the NCLFU observables is preferred by the global analysis. In Fig. 2 we show the likelihood in the space of a LFU contribution to C 9 vs. a purely muonic contribution to the linear combination C 9 = −C 10 , i.e. we consider a two-parameter scenario where the total NP (right). Solid (dashed) contours include (exclude) the Moriond-2019 results for R K and R K * . As R K only constrains a single combination of Wilson coefficients in the right plot, its 1σ contour corresponds to χ 2 = 1. For the other fits, 1 and 2σ contours correspond to χ 2 ≈ 2.3 and 6.2, respectively Wilson coefficients are given by 6 The best fit values in this scenario are C univ. = −0.44 with a χ 2 = 6.8 that corresponds to a pull of 6.5σ . The updated values of R K ( * ) favour a nonzero lepton flavour universal contribution to C 9 in this scenario.
One qualification is in order at this point. It is conceivable that a new effect in C 9 , and all the more the C univ. 9 contribution discussed above, is mimicked by a hadronic SM effect that couples to the lepton current via a virtual photon, for example charm-loop effects at low q 2 and resonance effects at high q 2 , see e.g. [81][82][83]. In our analysis, this possibility is taken into account in the uncertainty attached to the relevant observables that contribute to the (yellow) b → sμμ region in Fig. 2. Specifically, non-factorizable effects are parameterized as in [59] which, at 1σ , envelops the hadronic 6 Such decomposition was adopted for the first time in [80], to which we refer the reader for additional scenarios beyond the one we consider. We note that a shift in C univ. 10 would not produce a good overall fit. This may be appreciated from Fig. 1 (left). A C univ. 10 shift would only move the (yellow) b → sμμ region vertically, hence it would not help reach better agreement with the (blue) NCLFU region. We therefore set nonmuonic C 10 contributions to zero for simplicity. , ∀ , and a muonspecific contribution to the linear combination C 9 = −C 10 (see text for details). Solid (dashed) contours include (exclude) the Moriond-2019 results for R K and R K * effects identified in [10,84]. With such 'standard' procedure (adopted e.g. also in [85][86][87]), the global fit in Fig. 2 requires a non-SM C univ.
More conservative assumptions about the aforementioned hadronic uncertainties would not impact any of the observables that have been updated since Moriond 2017, when some of us performed a detailed study [29] of enlarged hadronic uncertainties, building on a similar study in Ref. [59]. By assuming non-form-factor hadronic uncertainties that are doubled with respect to the assumption denoted above as 'standard' and used in the present study, the global fit would be compatible with a vanishing NP contribution to C univ. 9 at the 1σ level.
With the above qualification in mind, in the following Sections we will address the question whether the above scenario, with a muon-specific contribution, C bsμμ 9 = −C bsμμ 10 , plus a universal contribution, C univ. 9 , can be justified from the UV point of view, assuming that both contributions are due to new physics. We start from a general discussion of how such contributions can arise in the SMEFT (Sect. 3.2) and discuss in detail how a lepton-flavour universal NP effect, C univ. 9 , can arise from RG effects in Sect. 3.3.

The global picture in the SMEFT
In this Section we summarize the interpretation of the above results within the SMEFT. In contrast to the discussion in WET at the b-quark scale, more Wilson coefficients become relevant in SMEFT due to RG mixing above [88][89][90] and below [91,92] the EW scale. Due to the pattern of Wilson coefficients preferred by the global fit, we focus on SMEFT Wilson coefficients that either contribute to the semimuonic Wilson coefficients in the form C bsμμ 9 = −C bsμμ 10 or induce a LFU effect in C bs 9 . Our notation in the following will be such that refers to leptons below the EW scale and l to the lepton doublets above the EW scale. Furthermore, we will work in a basis where generation indices for RH quarks are taken to coincide with the mass basis [93], which can be done without loss of generality.
The direct matching contributions to C 9,10 at the EW scale are well known [94,95], where the normalization factor N is defined in (7), the Z penguin coefficient c Z is and ζ = 1−4s 2 w ≈ 0.08 is the accidentally suppressed vector coupling of the Z to charged leptons. Using the notation of [96], the corresponding operators are given by [O (1) where q i , and l i are the left-handed SU (2) L doublet quarks and leptons and e i are the right-handed lepton singlets, ϕ is the SM Higgs doublet, and τ I are the Pauli matrices. Equations (19) and (20) highlight the well-known fact that a LFU contribution to C 9,10 induced by the SMEFT coefficients [C lq ] 2223 contribution, there can also be contributions from other Wilson coefficients. However, such additional contributions need, in general, be small (see e.g. [97]).
Apart from the direct matching contributions, additional SMEFT Wilson coefficients can induce a contribution to C 9 at the b mass scale through RG evolution above or below the EW scale, as pictured in Fig. 3. In view of the size of the effect preferred by the data, we can identify three qualitatively different effects that can play a role: from electroweak running above the EW scale. However, this solution is seriously challenged by EW precision tests [73]. (22) and (23), that induce a LFU contribution to C bs 9 from gauge-induced running both above and below the EW scale [75,98].
• Four-quark operators (defined below in Sect. 3.3.2) that also induce a LFU contribution to C bs 9 analogously to the semitauonic ones [99].
The case of semitauonic operators is particularly interesting as it potentially allows for a simultaneous explanation of the anomalies in neutral-current b → s transitions and in q 3 Diagrams inducing a contribution to C 9 through RG running above (left panel) and below (right panel) the EW scale. A sizeable contribution to C 9 is obtained when f = u 1,2 , d 1,2,3 or l 3 , see text for details lq ] 2223 at 2 TeV. All other Wilson coefficients are assumed to vanish at 2 TeV. Solid (dashed) contours include (exclude) the Moriond-2019 results for R K , R K * , R D , and R D * . For sets of data that effectively constrain only a single Wilson coefficient (namely R D ( * ) and NCLFU observables), 1σ contours correspond to χ 2 = 1. For the other data (b → sμμ and the global likelihood), 1 and 2σ contours correspond to χ 2 ≈ 2.3 and 6.2, respectively charged-current b → c transitions involving taus [53,98]. We now discuss these two possibilities in turn, from an effectivetheory point of view. Specific realizations in terms of simplified models will be discussed in Sect. 4.

Semi-tauonic operators
Intriguingly, a large value for [C (3) lq ] 3323 , that can explain the hints for LFU violation in charged-current b → c transitions (R D and R D * ), also induces a LFU effect in C 9 that goes in the right direction to solve the b → sμμ anomalies in branching ratios and angular observables. An additional lq ] 2223 (a = 1 or 3) of similar size can accommodate the deviations in R K and R K * . Since the linear combination [C (1) lq ] 3323 generates a sizable contribution to B → K ( * ) νν decays [95] that is constrained by B-factory searches for these modes, such models are only viable if the semitauonic singlet and triplet Wilson coefficients are approximately equal. 7 For reference, the leadinglog contribution from [C (1,3) lq ] 3323 to C 9 due to gauge mixing is given by 8 Figure 4 shows the likelihood contributions from R D ( * ) , NCLFU observables, b → sμμ observables, and the global likelihood in the space of the two Wilson coefficients [C (1) lq ] 2223 at the renormalization scale μ = 2 TeV. It is interesting to note that before the Moriond 2019 updates (indicated by the dashed contours), for a purely muonic solution with [C (1,3) lq ] 3323 = 0 (corresponding to the vertical axis), the best-fit values for NCLFU and b → sμμ data were in perfect agreement with each other (even though R D ( * ) cannot be explained in this case). Including the R K ( * ) updates, the best-fit point of the NCLFU and b → sμμ data instead lies in the region with non-zero semitauonic Wilson coefficients, just as required to explain the R D ( * ) anomalies. In fact, the agreement between the 1σ regions for R K ( * ) & D P 4,5 , R D ( * ) , and b → sμμ improves compared to the case without the R K ( * ) updates. We note that a further improvement of the fit is achieved by taking into account the Moriond 2019 update of R D ( * ) by Belle [50,51], which moves the 1σ region for R D ( * ) slightly closer to the SM value, exactly to the region where the contours of NCLFU and b → sμμ observables overlap. The best fit values in this scenario are [C (1,3) lq ] 3323 = −5.0 × 10 −2 TeV −2 1 that corresponds to a pull of 7.8σ . The pull is considerably larger in the present scenario than in those discussed in Sect. 3.1 since it can also explain discrepancies in b → c transitions. Plugging the above values into Eq. (25) one finds a good agreement with the universal-C 9 shift found in Sect. 3.1.3.
It is interesting to use the global fit in this scenario as the basis for predictions of several observables that are sensitive to the Wilson coefficients [C (1,3) lq ] 3323 and [C (1,3) lq ] 2223 and are supposed to be measured with higher precision in the near future. We collect predictions for LFU ratios, angular observables, and branching ratios in B and B s decays in Table 2.
In the previous discussion we have focused on the global consistency between NCLFU and b → sμμ data, which after Moriond 2019 would call for a universal C 9 contribution, which in turn is quantitatively consistent with the NP contribution demanded by R D ( * ) . An interesting question 9 is the comparison of this SMEFT scenario with the SM when taking into account LFU data alone, i.e. omitting b → sμμ data. This question is of interest e.g. because LFU data are deemed the theoretically cleanest.
To address this question we performed a fit taking into account only LFU observables in b → s and b → c ν. The two sets of observables are effectively orthogonal and therefore lead to fully consistent best-fit regions in the plane of the Wilson coefficients [C (1) lq ] 2223 , for both pre-and post-Moriond 2019 data. With the Moriond 2019 updates we find best fit values of the Wilson coefficients [C (1,3) lq ] 3323 = −5.0 × 10 −2 TeV −2 and [C (1,3) lq ] 2223 = 3.4 × 10 −4 TeV −2 corresponding to χ 2 = 5.3 and a pull of 4.9σ . 9 We thank the Referee for the suggestion.

Four-quark operators
Four-quark SMEFT operators can induce a LFU contribution to C 9 through gauge-induced RG running above and below the EW scale from diagrams like in Fig. 3. Since the flavourchanging quark current in O 9 is left-handed, these operators must contain at least two q fields and the other current must be a SU (3) c singlet. Neglecting CKM mixing (i.e. to zeroth order in the Cabibbo angle), the following operators could thus play a role: In the discussion to follow, we will consider each one of the above operators, and discuss whether they can produce at all a LFU C 9 contribution of the right size, in the light of all known low-energy constraints.
qq ] 2333 and [O (1) qu ] 2333 induce a LFU contribution to C 10 that is much bigger than the one in C 9 by a factor (1 − 4s 2 W ) −1 ∼ 13 through asbZ coupling generated by a top-quark loop, so they cannot explain the data. qq ] 23ii operators lead to enormous NP contributions to CP violation in D 0 -D 0 mixing that are excluded by observations [52]. Note that this happens even for real SMEFT Wilson coefficients since the CKM rotation between the mass basis for left-handed down-type quarks (relevant for b L → s L transitions) and uptype quarks (relevant for u L ↔ c L ) is itself CP violating. For example, the operator [O (1) qq ] 2311 contributes to D 0 -D 0 mixing through the operator qq ] 2311 where λ is the Cabibbo angle. For a typical value of [C (1) qq ] 2311 ∼ (1 TeV) −2 we find C D ∼ (10 3 TeV) −2 , which exceeds existing bounds [100] by an order of magnitude. Working instead in the basis where the up-type quark mass matrix is diagonal 10 and only using operators [Ô (a) qq ] ii j j that do not contribute to D 0 -D 0 mixing, it turns out that all these operators either generate excessive contributions to CP violation in K 0 -K 0 mixing 11 or do not generate an appreciable contribution to C 9 .
• [O (1) qu ] 2311 and [O (1) qd ] 2311 can induce a NP contribution to / [101,102] (i.e. direct CP violation in K 0 → ππ) through RG running above the EW scale, but for a low enough scale they can lead to a visible effect in C 9 without violating this bound.
• [O (1) qd ] 2322 can induce a NP contribution to M s , the mass difference in the B s -B s system, through RG running above the EW scale, but for a low enough scale it can lead to a visible effect in C 9 without violating this bound.
• An effect in C 9 generated by [O (1) qu ] 2322 and [O (1) qd ] 2333 is not strongly constrained from the point of view of lowenergy data, and is thus allowed from these data alone. We will discuss possible contributions to these operators from the point of view of UV completions in Sect. 4.2, see Eqs. (43)- (44).
To summarize, a LFU contribution to C 9 could, barring cancellations, be generated by the SMEFT Wilson coefficients The generic size of these Wilson coefficients required for a visible effect in C 9 is in the ballpark of 1/TeV 2 . Consequently, realistic model implementations of such an effect have to rely on tree-level mediators with sizeable couplings to quarks and masses potentially in the reach of direct production at LHC. We will discuss such simplified models in Sect. 4.2. 10 Operators and couplings in such up-aligned basis are here and henceforth denoted with a hat. 11 As well known, see [100] for the most recent study, this constraint is violated for new-physics scales below O(10 4 ) TeV, which is way above the scale required by B discrepancies.

Simplified new-physics models
In this section we consider simplified models with a single tree-level mediator multiplet giving rise to the Wilson coefficient patterns that are in agreement with the above findings in the EFT.
In Sect. 4.2, we discuss realizations of LFU contributions to C 9 via RG effects from four-quark operators as discussed in Sect. 3.3.2. We will show that there is a single viable mediator, a scalar transforming as (8, 2) 1/2 under the SM gauge group, and that it is strongly constrained by LHC dijet searches.
4.1 Explaining the data by a single mediator: the U 1 leptoquark solution As anticipated above, the only single mediator that can yield non-zero values for lq ] 2223 is the U 1 vector leptoquark, which transforms as (3, 1) 2/3 under the SM gauge group. We define its couplings to the left-handed SM fermion doublets q and l as From the tree-level matching at the scale = M U , one finds Consequently, for a given leptoquark mass, a τ -channel contribution to R D ( * ) depends only on g 32 lq and g 33 lq , while a μchannel contribution to R K ( * ) depends only on g 22 lq and g 23 lq . We perform an analysis at leading-logarithm accuracy, i.e. we take into account the one-loop RG running and mixing and perform the matching at tree level. As has been shown in [98,116], the one-loop matching of a simplified 12 U 1 leptoquark model would lead to • a 13% enhancement [116] of the tree-level contribution to semileptonic operators, • contributions to leptonic operators inducing τ → μνν, • contributions to semi-leptonic operators inducing B → K ( * )ν ν, • contributions to electric and chromomagnetic dipole operators in the WET.
Except for the last point, all of these contributions either only correct an already included tree-level result or are small compared to the current experimental sensitivity or compared to RG effects that contribute to the same operators (cf. [98]). We thus neglect them in the following. On the other hand, the one-loop matching contributions to the electric and chromomagnetic dipole operators can lead to relevant shifts in the Wilson coefficient C 7 at the b-quark scale, which are constrained by measurements of b → sγ observables (cf. [64]). In order to be sensitive to this possibly important effect, we will take the one-loop matching contributions to the SMEFT quark-dipole operators into account. ' as [96] [ The matching result depends on the couplings of the U 1 vector leptoquark to the SM gauge bosons, which can be written as where These couplings are determined by SM gauge invariance except for the two parameters k s and k Y . In the following, we make the choice k s = k Y = 1, which leads to a cancellation of divergent tree-level diagrams in U 1 -gluon and U 1 -B-boson scattering [117] and further avoids logarithmically divergent contributions to the dipole operators [37], making them finite and gauge independent. We note that k s = k Y = 1 is automatically satisfied in any model in which the U 1 leptoquark stems from the spontaneous breaking of a gauge symmetry but can also be realized for a composite U 1 [118]. We perform the one-loop matching to quark-dipole operators at the scale = M U by computing the diagrams shown in Fig. 5. Working in the basis in which the downtype Yukawa matrix is diagonal, and using the conventions mentioned above, we find the Wilson coefficients of the EW dipole operators where Y b and Y s denote the Yukawa couplings of the b and s quark respectively and a summation over the lepton index is implied. The Wilson coefficients of the chromomagnetic dipole operators at the matching scale read Using the tree-level matching conditions from SMEFT onto WET [119,120], we have checked that these results are consistent with the findings in [98].

R D ( * ) and indirect constraints
Within the defined framework we now search for a region of the leptoquark parameter space that explains all the B anomalies while at the same time avoids indirect low-energy constraints. We perform a fit with fixed M U = 2 TeV in the space of tauonic couplings g 32 lq and g 33 lq , which we take to be real for simplicity. This allows us to determine the region in which R D ( * ) can be explained by the semi-tauonic operators discussed in Sect. 3.3.1. The results of the fit are shown in Fig. 6 and our findings are as follows: • The strongest constraints are due to -leptonic tau decays τ → νν, which receive a contribution due to RG running, 13 13 Our analysis includes RG-induced logarithms. Note that the interactions in Eq. (29) and (33) provide a simplified model and not a complete UV theory of the U 1 -leptoquark. In such a UV theory, it could in principle be possible that the RG-induced logarithms are (partially) canceled by finite terms, which are not present in the simplified model. See e.g. discussion in [118]. Barring cancellations, and in view of the renormalization-scale independence of the full result -logarithms plus analytic terms -the contributions from the RG-induced logarithms usually provide a realistic estimate of the size of the effects. (assuming vanishing right-handed couplings), which is compatible with findings in [113]. This puts some tension on models based on a U (2) q flavour symmetry [37,106,108,113], where the natural expectation for the size of Based on the above results, we select a benchmark point from the best-fit region in the fit to tauonic couplings, which is also shown in Fig. 6. We then perform two fits in the space of muonic couplings g 22 lq and g 23 lq shown in Fig. 7: one for vanishing tauonic couplings (left panel) and one at the benchmark point g 32 lq = 0.6, g 33 lq = 0.7 (right panel). Our findings are as follows: • For vanishing tauonic couplings (left panel of Fig. 7), the data available before Moriond 2019 leads to a very good agreement between the fits to b → sμμ (orange contour) and NCLFU observables (dashed blue contour), while the R D ( * ) measurements cannot be explained in this scenario. Taking into account the updated and new measurements of R K ( * ) presented at Moriond 2019, one finds a slight tension between the fits to b → sμμ (orange contour) and NCLFU observables (solid blue contour).
14 Such contributions are, however, model-dependent. For example, they will be quite different in models with additional vector-like fermions running in the loops [105,106], as shown explicitly in Ref. [121]. Likelihood contours from different observables in the space of the tauonic U 1 leptoquark couplings g 32 lq and g 33 lq at 2 TeV. The grey areas are excluded at the 2σ level. R D ( * ) data and leptonic τ decays select a well-defined region in the g 32 lq versus g 33 lq plane. For R D ( * ) , which only constrain one degree of freedom, 1σ contours correspond to χ 2 = 1, while for others (the global likelihood, leptonic τ decays, BR(B → X s γ )), 1 and 2σ contours correspond to χ 2 ≈ 2.3 and 6.2, respectively. The dashed contour refers to pre-Moriond data of the corresponding solid contour This is analogous to the tension mentioned in Sects. 3.1.2 and 3.1.3.
• The tension disappears if one considers non-zero tauonic couplings that can also explain R D ( * ) , which is exemplified in the right panel of Fig. 7 for the benchmark point g 32 lq = 0.6, g 33 lq = 0.7. As discussed in Sect. 3.2, the semitauonic operators obtained from the tree-level matching (cf. Eq. 30) induce a lepton-flavour universal contribution to C 9 , which affects the predictions of b → sμμ observables and makes the fits to b → sμμ and NCLFU observables again compatible with each other at the 1σ level. Consequently, the deviations in neutral current and charged current B-decays can all be explained at once. This very well agrees with our findings in the SMEFT scenario in Sect. 3.3.1. • Given the presence of non-zero values for the tauonic couplings at the benchmark point, the strongest constraint on the muonic couplings g 22 lq and g 23 lq is due to LFV observables, in particular τ → φμ and B → K τ μ. The region in the space of muonic couplings that is excluded at the 2σ level by these observables is shown in gray in the right panel of Fig. 7.
Having identified a viable benchmark point, we conclude that the U 1 vector leptoquark can still provide an excellent description of the B anomalies while satisfying all indirect constraints.

Comparison between indirect and direct constraints
In addition to indirect constraints, highp T signatures of models containing a U 1 leptoquark have been discussed in detail considering current and future LHC searches [114,[124][125][126][127][128]. In this section, we compare direct constraints found in the latest study, [128], to the strong indirect constraints discussed in Sect. 4.1.1. To this end, we adopt the notation of [128] and use the parameters β i j L and g U , which are related to our notation by We perform a fit with fixed g U = 3, β 33 L = 1 (i.e. g 33 lq ≈ 2) in the space of M U and β 23 L . These values are chosen to allow for a direct comparison with the constraint from pp → τ τ shown in Fig. 1 of [128] and pp → τ ν shown in Fig. 6 of [128]. We include both of these direct constraints in our Fig. 8 as hatched areas. In addition, we show the results from our fit, namely the constraint from leptonic τ decays and the region preferred by the R D ( * ) measurements. Our findings are as follows: • The indirect constraint from leptonic τ decays is stronger than the direct constraints in nearly all of the parameter space shown in Fig. 8, except for a small region at large β 23 L 0.75, where the constraint from pp → τ ν is the strongest one.
• In the region where R D ( * ) can be explained, the indirect constraint from leptonic τ decays is considerably stronger than the direct ones.
• Small values for as naturally expected in models based on a U (2) q flavour symmetry [37,106,108,113] require a relatively small mass M U to explain R D ( * ) . Thus, as also pointed out in [113,124], there is already some tension between this natural expectation and the direct searches.
We note that the direct constraints shown in Fig. 8 depend on the coupling strength g U . While the assumptions g U = 3, β 33 L = 1 lead to a lower bound on the leptoquark mass M U 2.7 TeV, this bound does not apply to the scenario discussed in  39)). The green region is the preferred region from R D ( * ) , while the gray shaded area is excluded by leptonic tau decays at the 2σ level. The hatched areas are excluded by LHC searches recasted in [128] Sect. 4.1.1, which features considerably smaller couplings. 15 Latest direct constraints from U 1 pair production that are independent of the coupling strength g U only exclude masses M U 1.5 TeV [114,128]. Therefore, the scenario discussed in Sect. 4.1.1 is currently not constrained by direct searches.

Lepton flavour universal C 9 from leptophobic mediators
As discussed in Sect. 3.3.2 in the SMEFT context, a lepton flavour universal contribution to C 9 can also be induced from a four-quark operator via RG effects. However, the fourquark operator would realistically have to be generated by the tree-level exchange of a resonance with mass not exceeding a few TeV and O(1) couplings. Such resonance could then be in reach of direct LHC searches, apart from other indirect constraints.
Since it was shown in Sect. 3.3.2 that the only viable operators are of type O (1) qu and O (1) qd , the conceivable tree-level mediators, excluding fields that admit baryon number violating diquark couplings, 16 are (see e.g. [129]) • (1, 1) 0 with spin 1 (Z ), 15 The partonic cross section relevant for the direct constraints in Fig. 8 scales as σ ∼ (g U /M U ) 4 [128]. 16 Note that baryon number violating diquark couplings do not necessarily lead to tree-level proton decay.
Such leptophobic mediators are discussed here mostly because they are a logical possibility -in particular, they are not necessarily motivated from a UV perspective. It is immediately clear that the spin-1 mediators Z and G are not viable: to generate the Wilson coefficients [C (1) qu ] 23ii or [C (1) qu ] 23ii with sufficient size, they would require a sizeable flavour-violating coupling to left-handed down-type quarks that would induce excessive contributions to B s -B s mixing. Thus the only potentially viable models are the scalar mediators.
The new scalar bosons have the following Lagrangian couplings to quarks, The SMEFT Wilson coefficients that are generated by a treelevel scalar exchange and can contribute to C 9 via RG effects read where R = u, d and (c H , c ) = (−1/6, −2/9). Clearly, to generate one of the down-type Wilson coefficients [C (1) qd ] 23ii in (28), at least one flavour-violating coupling in the down-aligned basis has to be present, leading to excessive contributions to B 0 -B 0 or B s -B s mixing. Thus we assume vanishing down-type couplings y i j Xqd = 0 in the following.
For the up-type couplings y Xqu , it is instead convenient to work in the up-aligned basis (cf. footnote 10 for notation), as settingŷ 12 Xqu =ŷ 21 Xqu = 0 allows avoiding dangerous contributions to D 0 -D 0 mixing. We then obtain the following non-vanishing matching conditions relevant for RG-induced contributions to C 9 ,  Numerically, it turns out that a visible lepton flavour universal effect in C 9 requires O(1) couplings for a scalar mass of the order of 2 TeV even for terms without CKM suppression. Thus, LHC searches for di-jet resonances are sensitive to the new scalars. If the NP effect in C 9 is generated through the terms in (43), the scalar has sizeable couplings to up quarks, leading to an excessive production cross section at the LHC.
To avoid this, we will further assumeŷ 1i Xqu =ŷ i1 Xqu = 0 and use the terms in (44). Production in pp collisions is still possible via charm quarks. The leading-order cross sections of the charged and neutral components read where (c H ,c ) = (1, 4/3), √ s being the center of mass energy and L i j are the parton luminosities as defined in [130] and we have neglected contributions suppressed by CKM factors.
We confront these cross sections with exclusion limits from ATLAS [131] and CMS [132]. Our procedure to obtain constraints on the scalar model parameters is detailed in Appendix D. For definiteness, we choose X = in the following. In fact, according to the results of [99], the H case is considerably more constrained by contributions to B → X s γ introduced radiatively at the two-loop level.
The left plot in Fig. 9 shows contours of C univ. 9 in the plane ofŷ 22 qu vs. the color-octet scalar mass for a scenario in whichŷ 32 qu = −1. We also show the 95% C.L. exclusion from dijet resonance searches at LHC. Obviously, a visible negative contribution to C univ. 9 of at most −0.3 can only be generated in a thin sliver of parameter space for masses above 2 TeV and withŷ 22 qu ∼ 1. This scenario is on the brink of exclusion. 17 The bending of the C univ. 9 contours at low masses in the left plot of Fig. 9 is due to the fact that there is a CKM-suppressed contribution even forŷ 22 qu = 0 from the third term in (44). To investigate whether such scenario, where production is only possible through a b quark PDF, is allowed, in the right plot of Fig. 9 we show the C univ. 9 contours and the LHC exclusion forŷ 32 qu vs. the color octet scalar mass settingŷ 22 qu = 0. Clearly, generating an appreciable contribution to C univ. 9 is excluded by di-jet searches in this scenario.

Conclusions
Motivated by the updated measurements of the theoretically clean lepton flavour universality tests R K and R K * by the LHCb and Belle experiments, as well as by additional measurements, notably of B s → μμ by the ATLAS collaboration, we have updated the global EFT analysis of new physics in b → s transitions. A new-physics effect in the semimuonic Wilson coefficient C bsμμ 9 continues to give a much improved fit to the data compared to the SM. However, compared to previous global analyses, we find that there is now also a preference for a non-zero value of the semi-muonic Wilson coefficient C bsμμ 10 , mostly driven by the global combination of the B s → μ + μ − branching ratio including the ATLAS measurement. The single-coefficient scenario giving the best fit to the data is the one where C bsμμ 9 = −C bsμμ 10 , which is known to be well suited to UV-complete interpretations, and indeed is predicted in several new-physics models with tree-level mediators coupling dominantly to left-handed fermions.
We have also studied the possibility of a simultaneous interpretation of the b → s data and the discrepancies in b → cτ ν transitions in the framework of a global likelihood in SMEFT Wilson coefficient space. We find one especially compelling scenario, characterised by new physics in all-lefthanded semitauonic four-fermion operators. These operators can explain directly the discrepancies in b → cτ ν transitions, and, at the same time, radiatively induce a lepton flavour universal contribution to the b → s Wilson coefficients. An additional nonzero semimuonic Wilson coefficient then allows accommodating the R K ( * ) discrepancies. Such picture can be quantitatively realized in the context of the U 1 leptoquark simplified model, and we find that indeed an excellent description of the data can be obtained, including the deviations in b → cτ ν transitions.
Another logical possibility to generate a lepton flavour universal NP effect in C 9 is via RG effects from a four-quark operator. We have investigated this possibility in the SMEFT and in simplified tree-level models. We find that the only setup potentially viable in the light of indirect constraints is a colour-octet scalar. However, due to its TeV-scale mass and large coupling to quarks, it is strongly constrained by di-jet resonance searches at the LHC, which put it on the brink of exclusion.
Our study illustrates how the theoretical picture has evolved as a consequence of crucial measurement updates, and how this picture stays coherent in spite of the numerous constraints. We collect in Table 2 a number of predictions directly related to the discussed SMEFT scenario. The situation will only get more exciting due to the host of new analyses using the full Run-2 data set, as well as the Belle-II data set, to which we look forward.
Acknowledgements D.S. warmly thanks the organisers of the Rencontres de Moriond 2019 on "Electroweak Interactions and Unified Theories" for the opportunity to present the results in this paper prior to their appearance on the arXiv. We also acknowledge useful remarks from Ben Allanach and Sébastien Descotes-Genon. The work of D. S. and J. A. is supported by the DFG cluster of excellence "Origin and Structure of the Universe". The research of W. A. is supported by the National Science Foundation under Grant No. NSF 1912719. The numerical analysis has been carried out on the computing facilities of the Munich Computational Center for Particle and Astrophysics (C2PAP).

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study, with no experimental data generated.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .
In all three cases, the measurements are correlated, since the B 0 and B s have a similar mass, such that the signal regions in dimuon invariant mass squared overlap. ATLAS and CMS provide two-dimensional likelihood contours, from which we interpolate the two-dimensional likelihoods, while LHCb directly provides the two-dimensional likelihood numerically. The three resulting likelihoods are shown as thin lines in Fig. 10 and compared to the SM central values.
Next, we assume the three experiments to be uncorrelated (which we assume to be a good approximation given the dominance of statistical uncertainties) and combine the two-  Fig. 10. For our global likelihood in Wilson coefficient space, we need to make an additional approximation, namely that the experimental likelihood is approximately Gaussian (see [53] for a discussion). Thus we fit a two-dimensional Gaussian to the product likelihood. The resulting contours are shown as thick dashed lines in Fig. 10.
Since throughout our numerical analysis, we never consider NP effects in b → dμμ transitions, it is also interesting to compare the combined confidence regions for the B s → μ + μ − branching ratio, fixing B 0 → μ + μ − to its SM central value or profiling over it. We find We stress that the similarity of these two numbers is not trivial, as for individual measurements, especially the CMS and ATLAS ones, the best-fit value for BR(B 0 → μ + μ − ) deviates strongly from the SM prediction. Conversely, for The values for the Gaussian approximation only differ from these numbers in a negligible way. Compared to the SM predictions 18 : we then find the following one-dimensional pulls: 19 • if both branching ratios are SM-like, 2.3σ , 20 18 For the SM values, we have used flavio v1.3 with default settings. The B s → μ + μ − branching ratio refers to the time-integrated one, see Refs. [133,134] for state-of-the-art discussions on the SM uncertainty. The decay constants are taken from the 2019 FLAG average with 2 + 1 + 1 flavours [135], V cb = 0.04221(78) from inclusive decays, and V ub = 0.00373 (14) from B → π ν. 19 Here, the "one-dimensional pull" is −2 times the logarithm of the likelihood ratio at the SM vs. the experimental point, after the experimental uncertainties have been convolved with the covariance of the SM uncertainties. 20 Converting the likelihood ratio to a pull with two degrees of freedom, one obtains 1.8σ ; this is why the star in Fig. 10 is inside the 2σ contour.

B. Global likelihood of the standard model
In view of the sizable pulls in various NP scenarios considered in this work, an interesting question is how good the overall agreement of the SM with the data is. To this end, we consider the value of the global likelihood at the SM point, L( 0), or χ 2 SM ≡ −2 ln L( 0). Here, the normalization of the likelihood is such that χ 2 = 0 corresponds to the case where all measurements are in exact agreement with the theoretical predictions. By means of Wilks' theorem, this χ 2 value can be translated to a p-value, quantifying the goodness of fit.
However, the intepretation of this global χ 2 value is not straightforward as it is subject to several ambuguities.
• Since the global likelihood contains many observables not sensitive to the Wilson coefficients that we consider in our analysis, which focuses on discrepancies in B physics, this p-value depends on the number of observables included in the test. • The likelihood contains measurements that are actually averages (e.g. by PDG or HFLAV) of several measurements, such that the number of degrees of freedom is reduced and a constant contribution to the χ 2 coming from a tension between the averaged measurements is concealed. • Due to the statistical approach, where theory uncertainties are combined with the experimental ones and no explicit nuisance parameters are present, the contribution to the absolute χ 2 from different sectors can depend sensitively on the central values chosen for the parameters. For instance, a lower V cb value would lead to better agreement of b → s branching ratios (via lower |V tb V * ts |) but worse agreement with b → c ν transitions. While this issue does not affect the χ 2 used in our numerical analysis, it makes the interpretation of the absolute χ 2 for individual sectors difficult.
With these caveats in mind, we provide in Table 3 the absolute χ 2 values obtained with smelli for various subsets of observables. The number N obs in the first column refers to the number of observations, i.e. independent measurements of an observable, which is greater than or equal to the number of observables. In our case, N obs is to be interpreted as the χ 2 's number of degrees of freedom. Through the chi2_dict method introduced in v1.3, it is possible to explore the χ 2 also for different observable sets or parameter inputs. C. C 9 vs. C 9 = −C 10 As already discussed in Sect. 3.1.1 and summarized in Table  1, we find a preference for the C bsμμ 9 = −C bsμμ 10 scenario over the C bsμμ 9 -only scenario. Since the opposite result was found in previous analyses, like e.g. [29], some further comments are in order. The change in preference is not related to the updated measurements of R K and R K * but can be traced back to the inclusion of several observables into the likelihood that were not considered in [29]. 21 Among the newly included observables mentioned in Sect. 2, the following are in particular relevant here.
Like the theoretical predictions for BR(B s,d → μ + μ − ), also the predictions for these observables depend on the B s,d -meson decay constants F B s,d and the CKM parameters. As detailed in [53], these nuisance parameters are "integrated out" and enter the likelihood in terms of a covariance matrix that contains all experimental and theoretical uncertainties together with their correlations. Due to these correlations, experimental information on F = 2 observables can reduce the theoretical uncertainties of BR(B s,d → μ + μ − ). Such a reduction of theoretical uncertainties has been discussed explicitly for M s,d and BR(B s,d → μ + μ − ) in models with minimal flavour violation in [140]. Via the correlations, F = 2 observables can have an impact on a fit to the Wilson 21 Also other post-Moriond 2019 global fits that appeared slightly before and after the preprint of this paper use different sets of observables and found a preference for the C bsμμ 9 -only scenario [97,[136][137][138][139]. 22 For the full list of observables, see [53].   and C bsμμ 10 even if they do not directly depend on these coefficients themselves.
To illustrate the effect of including the above listed observables into the likelihood, we show 1σ contours for several fits to subsets of observables in Fig. 11.
• The pink contour corresponds to a fit including only b → sμμ observables excluding BR(B s,d → μ + μ − ). This fit is equivalent to the one shown in [29] and clearly prefers the C  -only scenario and the C bsμμ 9 = −C bsμμ 10 scenario perform similarly well and neither of them is clearly preferred.
• The red contour corresponds to a fit that additionally includes F = 2 observables. As described above, taking them into account can have an effect on the uncertainties of B s → μμ and we find that this leads to an increased preference for a non-zero NP contribution to C bsμμ 10 . In this case, C Other observables that are correlated with the b → sμμ observables (e.g. b → sγ ) only have a marginal impact on the fit. Therefore, the orange contour in Fig. 11 essentially coincides with the orange contour in Fig. 1 left.
Note that the main effect of the correlations between the C bsμμ 10 -dependent observables and the F = 2 observables vanishes if the theory predictions of the F = 2 observables exactly equal the corresponding experimental values. This   All other Wilson coefficients are assumed to vanish. Solid (dashed) contours include (exclude) the Moriond-2019 results for R K and R K * assume the theory predictions of F = 2 observables to be SM-like. 23 23 Many NP models that explain the deviations in rare B decays also predict a NP contribution to the Wilson coefficient C bsbs V L L of the operator (s L γ μ b L ) (s L γ μ b L ), which modifies the theory prediction of M s . A large class of models predicts M s > M SM s [141][142][143]. We have explicitly checked that a contribution to C bsbs V L L that increases M s up to its 2σ experimental bound can only marginally decrease

D. Di-jet bounds on leptophobic mediators
As discussed in Sect. 4.2, generating a LFU contribution to C 9 radiatively from four-quark operators requires a fairly light scalar mediator with strong coupling to quarks, that can lead to a signal in di-jet resonance searches at the LHC. In our numerical analysis, we employ two recent di-jet resonance searches.
• A search at √ s = 13 TeV by CMS based on 36 fb −1 and covering the mass range from 600 to 8100 GeV [132]. We apply the 95% C.L. constraint on the production cross section times branching ratio times acceptance, assuming −0.14 −0.12 −0. 10  100% decay into di-jets and a constant acceptance of 57%.
• A search at √ s = 13 TeV by ATLAS based on 80 fb −1 employing initial state radiation to cover the low-mass range not covered by CMS [131]. We apply the 95% C.L. constraint on the production cross section times branching ratio times acceptance times efficiency, assuming 100% decay into di-jets and using the mass-dependent efficiency and acceptance for the Z model given in the publication.
The 95% C.L. bounds under these assumptions are shown in Fig. 12.

E. WET and SMEFT scenarios with Wilson coefficient pairs
The aim of this appendix is to present additional projections of the likelihood onto pairs of Wilson coefficients both in WET and SMEFT.
• The plots of Fig. 13 complement the ones in Fig. 1 and show Wilson coefficient spaces involving muonic axialvector currents C bsμμ 10 and C bsμμ 10 . • Fig. 14 shows the space of electronic Wilson coefficients C bsee 9 and C bsee 10 . R K and R K * can be explained simultaneously if both C bsee 9 and C bsee 10 are present.   (Fig. 17). The fits prefer C bs 7 to be SM-like but show a slight (∼ 1σ ) preference for a positive shift in C bs 7 . • In Fig. 18 we consider the effects of non-zero scalar and pseudo-scalar Wilson coefficients that obey the SMEFT  Fig. 19 show scenarios with various combinatios of SMEFT Wilson coefficients. The left plot in Fig. 19 contains only tauonic Wilson coefficients, while the right plot shows a scenario with lepton flavour universal coefficients. The latter case is strongly constrained by electro-weak precision observables. • Finally, Fig. 20 shows that the b → sμμ anomalies can be explained by SMEFT 4-quark operators that enter the rare semi-leptonic decays at the loop level.