Charming New $B$-Physics

We give a comprehensive account of the flavour physics of Beyond-Standard-Model (BSM) effects in $b \to c \bar{c} s$ transitions, considering the full set of 20 four-quark operators. We discuss the leading-order structure of their RG mixing with each other as well as the QCD-penguin, dipole, and FCNC semileptonic operators they necessarily mix with, providing compact expressions. We also provide the first complete results for BSM effects in the lifetime observables $\Delta \Gamma_s$ and $\tau(B_s)/\tau(B_d)$, as well as for the semileptonic CP-asymmetry $a_{sl}^s$. From a global analysis, we obtain stringent constraints on 16 of the 20 BSM operators, including the 10 operators $Q^{c\prime }_{1 \dots 10}$ involving a right-handed strange quark. Focussing on CP-conserving new physics, the constraints correspond to NP scales of order 10 TeV in most cases, always dominated by exclusive and/or radiative $B$-decays via RGE mixing. For the remaining four operators, including the two Standard-Model (SM) ones, larger effects are experimentally allowed, as previously noted in arXiv:1701.09183. We extend that paper's scope to the CP-violating case, paying attention to the impact on the decay rate and time-dependent CP-violation in $B_d \to J/\psi K_S$. We show that quantifiable constraints arise for new physics in either of the two SM operators, with the uncertain non-perturbative matrix element of the colour-suppressed operator determined from the data. For new physics in the SM suppressed coefficient $C^c_1$ we find two perfectly viable, narrow bands of complex Wilson coefficients. Somewhat curiously, one of them contains a region where the fitted matrix element for the colour-suppressed operator is in agreement with naive factorization, contrarily to a widely held belief that large non-factorizable contributions to $B_d \to J/\psi K_S$ are implied by experimental data.

Abstract: We give a comprehensive account of the flavour physics of Beyond-Standard-Model (BSM) effects in b → ccs transitions, considering the full set of 20 four-quark operators. We discuss the leading-order structure of their RG mixing with each other as well as the QCD-penguin, dipole, and FCNC semileptonic operators they necessarily mix with, providing compact expressions. We also provide the first complete results for BSM effects in the lifetime observables ∆Γ s and τ (B s )/τ (B d ), as well as for the semileptonic CP-asymmetry a s sl . From a global analysis, we obtain stringent constraints on 16 of the 20 BSM operators, including the 10 operators Q c 1...10 involving a right-handed strange quark. Focussing on CP-conserving new physics, the constraints correspond to NP scales of order 10 TeV in most cases, always dominated by exclusive and/or radiative B-decays via RGE mixing. For the remaining four operators, including the two Standard-Model (SM) ones, larger effects are experimentally allowed, as previously noted in [1]. We extend that paper's scope to the CP-violating case, paying attention to the impact on the decay rate and time-dependent CP-violation in B d → J/ψK S . Contrary to common lore, we show that quantifiable constraints arise for new physics in either of the two SM operators, with the uncertain non-perturbative matrix element of the colour-suppressed (or equivalently, colouroctet) operator determined from the data. For new physics in the coefficient C c 1 , suppressed in the SM, we find (in addition to CP-conserving new physics) two perfectly viable, narrow bands of complex Wilson coefficients. Somewhat curiously, one of them contains a region where the fitted matrix element for the colour-suppressed operator is in agreement with naive factorization, contrarily to a widely held belief that large non-factorizable contributions to B d → J/ψK S are implied by experimental data.

Introduction
A wide range of B meson decays are affected by b → ccs transitions, providing the opportunity to use a rich set of complementary observables to detect possible new physics (NP), or place constraints on such dynamics and its mass scale. These partonic transitions are generated at the tree level in the SM, and contribute again at tree level to lifetime observables such as ∆Γ s and τ (B s )/τ (B d ) [2], which stand out among others through their good theoretical control. Moreover, as we have shown in [1], there is the possibility that operators with this flavour structure can be involved in the rare B-decay anomalies [3][4][5][6][7][8][9][10][11][12][13][14] 1 . In this case, the contribution to the observables is through a charm loop, which radiatively generates the semileptonic Wilson coefficient C 9V , conjectured to lie behind the so-called P 5 anomaly. 2 That BSM b → ccs transitions can give important one-loop contributions to rare radiative and semileptonic decays should not really be surprising, given the same happens in the SM, where charm loops provide on the order of one-half of the b → sγ decay amplitude, as well as of C 9V . The size of these effects is amplified by strong renormalization-group (RG) running in the SM, which can be even stronger for certain BSM operators [1]. Neither does the NP scale need to be particularly low: because of strong RGE running, accounting for the P 5 anomaly merely requires a BSM contribution |C BSM 1 | ∼ 0.1, giving a naive NP scale For a weakly coupled, sufficiently leptophobic tree-level mediator, this may be allowed by high-p T LHC searches, and will not cause problems with electroweak precision observables. For strong coupling, the NP scale can be as high as 30 TeV and out of the reach of LHC direct searches altogether [1]. While the flavour anomalies provided some extra motivation for our earlier work [1], it is clear that our BSM operators can only provide a lepton-flavouruniversal contribution, and as such can only be a part of a more complete solution (see [16] and [17][18][19] for discussion of flavour-universal and flavour-non-universal combinations). Therefore, the scale suggested above is likely to be a lower bound if our operators contribute in this fashion. In the rest of this paper however, we will be focussing our attention on the BSM effects arising from our full basis of operators in the observables already discussed. A BSM model that generates new physics only in b → ccs is not realistic, as is already evident from the leading-order RG mixing we have mentioned. Thus the model-independent study presented here has to be considered as a building block in constructing or constraining UV models.
In the present paper, we extend the analysis of [1] in the following ways: • We study a full basis of 20 b → ccs operators, where [1] focused on those 4 that can generate the P 5 anomaly, and we obtain stringent constraints on the 16 others.
• We extend the case of C c 1−4 to allow for CP-violating NP and show how to include B d → J/ψK S data with minimal non-perturbative inputs.
To address the first item, we build the anomalous dimension matrix governing the RG evolution of the 34 relevant operators, and review the all-order block structure. We explain that, while some blocks first arise at one-loop order and others at two-loop order, an unambiguous leading-order evolution operator arises, and obtain a compact numerical expression for it. This can be useful in understanding how different four-quark operator coefficients affect the rare and semileptonic decays, but may also be useful to a reader wishing to consider a particular UV model, particularly when multi-operator correlations are important that go beyond the one-and two-parameter cases we consider in our modelindependent phenomenology.
We also compute the complete BSM contribution to the ∆B = 0 and ∆B = 2 lifetime observables, including the impact on the semileptonic CP asymmetry a s sl . Regarding the second item, a key point is that once CP violation is switched on, contributions to time-dependent CP violation in B d → J/ψK S are generated, the sine coefficient of which is precisely measured and usually taken to provide a clean determination of the CKM angle β; this no longer holds in the present context. However, we find that as long as NP only affects one of the two Wilson coefficients present in the SM, the global data set is sufficient to determine the (complex) ratio of the two relevant non-perturbative matrix elements jointly with the complex Wilson coefficient, such that theory input is needed only for the matrix element of the operator Q c 1 , which on grounds of large-N arguments alone should be close to its naive-factorization value (as is also borne out by QCD factorization [20][21][22], even though QCDF is not expected to provide reliable quantitative results -see our main discussion in Section 3.5).
As phenomenological results, for the 16 operators not previously considered we provide constraints for the CP-conserving case, almost all of which turn out to be of order 10 TeV. In addition, we discuss in detail the constraints on the two one-parameter CP-violating scenarios where we have a complex BSM contribution to one of the two SM coefficients C c 1,2 . As a perhaps curious result, we find that among the allowed regions is one where CP-violating BSM affects C 1 in such a way that the matrix element of Q c 2 is also close to its naive-factorization value.
We note that BSM effects in hadronic tree-level decays have previously been systematically studied in [23,24]. In these works all hadronic decay channels of the b quark were considered and the new contributions to the Wilson coefficients were also allowed to be complex, but only the SM operators Q c 1 and Q c 2 were investigated. The layout of our paper is as follows: In Section 2 we describe the setup of our model, specifying our operator basis and giving the RG evolution. Section 3 is devoted to our theoretical results for the lifetime observables, radiative and rare decay contributions, as well as our approach to B d → J/ψK S , and constitutes one of two main parts of our paper. As the second main part, we describe our phenomenological constraints in Section 4 (details of the inputs we use are is given in Section 4.1). Our conclusions can be found in Section 5.

Operator basis
For our study of BSM b → ccs effects we use the following weak effective Hamiltonian The first part H cc eff contains a complete set of 20 b → ccs four-fermion operators, reading where i, j are SU (3) colour indices and ψ L,R denotes the projections ψ L,R = (1 ∓ γ 5 )/2 · ψ.
The remaining 10 operators Q c i are obtained from those displayed by letting L/R → R/L. In the SM the operators Q c 1 and Q c 2 arise due to a tree-level W exchange. All other 18 operators are genuine BSM effects. For the overall normalisation of these operators we have chosen the CKM normalisation of the SM operators, which contains a tiny imaginary part.
For radiative and semileptonic decays b → sγ and b → s + − , further operators contribute. These comprise the radiative and semileptonic operators and their parity conjugates, again denoted by primes. m b denotes the MS b-quark mass, e is the electromagnetic coupling, and α = e 2 4π . These operators enter through the second part H rsl eff of (2.1) (2.5) We note for the reader that since we do not have any effects in the Q ( ) 10 operators, the decay B s → µµ is purely SM-like in our scenario.
Moreover, they receive contributions from QCD penguin operators, contained in the third term H QCD eff of (2.1), (2.6) Our QCD penguin operators are defined in [25] and include the chromomagnetic dipole operator, reading explicitly: where q runs over all active quark flavours in the effective theory, i.e. q = u, d, s, c, b, and a prime again denotes a chirality conjugate. Our complete operator basis comprises 34 operators and closes under QCD renormalization. We renormalize our operators as in [25][26][27][28]. We will neglect the tiny contribution |λ u | ≈ 0.00084 3 [29] (for similar results see [30]) and we thus get from unitarity of the CKM matrix λ c = −λ t . To isolate the BSM contribution, we split the Wilson coefficients into SM parts and BSM parts, We will neglect the small mass ratios m q /m W and m s /m b , which implies the vanishing of all primed Wilson coefficients in the SM, C ,SM i (µ) = 0 and in addition the vanishing of the four-quark coefficients C c,SM i (µ) = 0 for i = 2. For our phenomenology we will also assume that at an input scale µ 0 ∼ M W , the ∆C i vanishes for all but the four-quark operators. This corresponds to what we called 'Charming Beyond the Standard Model' (CBSM) scenario in [1]. As described below, RG evolution then generates BSM contributions to the penguin and dipole operators, which play a crucial role in the phenomenology of the CBSM scenario. The CBSM scenario should be viewed as a partial effective description of a more complete UV scenario, which will in general also involve nonzero initial values for the other ∆C i .

Renormalization-group evolution
As emphasized in our previous work [1], operator mixing can have a dramatic impact on the radiative and semileptonic Wilson coefficients, and as a result on the contributions of the b → ccs operators to the radiative and rare semileptonic decays. It is therefore crucial to include its effects, which is conveniently done through renormalization-group evolution from the BSM scale µ 0 where the Wilson coefficients are initially obtained to a scale µ ∼ m b appropriate to evaluating B-physics observables.
Collecting the 17 unprimed Wilson coefficients into a vector C, the coupled system of renormalization-group equations (RGE) governing their dependence on the renormalization scale µ can be written in matrix notation as where γ is the anomalous-dimension matrix. Note that the unprimed and primed Wilson coefficients do not mix at any order in QCD; moreover, the 17 primed Wilson coefficients, collectively denoted C , fulfil the same set of equations (2.9), with identical anomalousdimension matrix. Both statements follow directly from the parity invariance of QCD. The solution of (2.9) is a linear relation in terms of the so-called evolution operator U (µ, µ 0 ), which itself satisfies (2.9). Importantly, the SM and BSM parts C SM and ∆ C separately satisfy (2.10).
If the leading-order anomalous-dimension matrix has the form where γ (0) is a constant matrix, then the explicit form of the evolution operator to leading logarithmic order is where β 0 is the leading coefficient in the QCD β-function, In the CBSM setup, two sources of complication arise. The first is that the leading-order anomalous-dimension matrix does not have the form (2.11). The second is that part of the anomalous-dimension matrix first arises at 2-loop order. To discuss both issues further, let us decompose the Wilson coefficient vector C into subcomponents, (2.14) The anomalous-dimension matrix then has the block form which holds to all orders of QCD. 4 The element γ 99 vanishes for our normalization of Q 9V , which makes it proportional to a conserved quark current, but we will consider a different normalization shortly. γ 77 , γ 88 , γ cp , andγ cc first arise at one-loop order, giving rise to the form (2.11) at leading order. γ p7 and γ p8 arises first at two-loop order. This makes them scheme-dependent already at leading order, though the explicit coupling factors in the definitions of Q 7γ and Q 8g imply that the form (2.11) is maintained. On the other hand, some of the elements of γ c7 , γ c8 , γ c9 and γ p9 are non-zero at one loop, but with no accompanying factor of α s . To make this more explicit, let us decompose the set of (cb)(sc) Wilson coefficients further, C = ( C 1...6 , C 7...10 ) ≡ ( C A , C B ), and further decompose the anomalous dimension accordingly, (2. 16) In the second expression we have indicated in what blocks factors of α s /(4π) do or do not arise at leading order. The loop-and coupling-order of the leading contributions to each block are also given in Table 1 (the full algebraic expressions for each block are given in Appendix A). In (2.16) we have also made explicit which blocks vanish to all orders in QCD. Besides those blocks whose vanishing is already expressed through (2.15), these are the blocks concerning the would-be mixing of the four-quark operators in class B, i.e. Q c 7...10 , with the QCD penguin operators and Q 9V . 5 As a result of the structure of (2.16), C A and C B do not mix; their renormalization, including the mixing of each into the QCD penguins, dipoles, and C 9V , can therefore be considered independently. The evolution in the situation where only one of C A or C B is present is described by an anomalous dimension of the form (2.15), withγ cc → γ XX ,γ cp →γ Xp , γ c8 → γ X8 , γ c7 → γ X7 , and γ c9 → γ X9 , where X = A or B as applicable; let us refer to these limiting cases as 'case A' and 'case B', respectively.

Block
Loop Order Coupling Order γ AA ,γ BB ,γ Ap ,γ pp 1 α Table 1 Let us first consider case A, which includes the subset of operators discussed in [1] and 5 The absence of this mixing follows from the fact that the chiralities of the (s, b) pair of fields differ between these groups of operators. Chirality conservation of massless QCD prevents mixing of operators of the same dimension but distinct quark chiralities, even for nonvanishing quark masses. For the same reason, Q c 5,6 do not mix into Q9V .
in particular includes the SM case. Here, a rescaling is sufficient to bring the anomalous-dimension matrix into the form (2.11), such that the solution (2.12) applies (note that with the rescaling, γ 99 → −2β 0 ). This also shows that C 9V formally starts at order 1/α s when C A is nonzero (which includes the SM case). The blocks γ A7 , γ A8 , γ p7 , and γ p8 , due to two-loop diagrams, are scheme-dependent and induce scheme dependence of C 7γ and C 8g already at the leading order. The scheme-dependence ultimately cancels out in observables. It is convenient and customary to define effective dipole coefficients [31,32] C eff 7γ (µ) = C 7γ (µ) + y · C Ap , (2.18) in such a fashion that the leading-order expression for B(B → X s γ) is proportional to |C eff 7γ | 2 (or |C eff 7γ | 2 + |C eff 7γ | 2 , in the presence of primed operators). This ensures that C eff 7γ and C eff 8γ are scheme-independent to leading order. In the CBSM scenario, with our choice of operator basis and renormalization, including in particular anticommuting γ 5 , the vectors y, z and C Ap read extending the SM case [25,[33][34][35].
As in the SM case, one can choose to go to a special scheme where the vectors y and z vanish by means of a finite renormalization of the four-quark operators i is any of the four-quark operators whose Wilson coefficients appear in C Ap . In this special scheme, C 7γ and C 8g are (to leading order) equal to C eff 7γ and C eff 8g , respectively, by construction. The mixing among the four-quark operators is unaffected by the change of scheme, but the anomalous dimension elements involving the dipole operators change. The resulting leading-order anomalous dimension matrix in the special scheme is sometimes referred to as γ (0),eff .
Let us now turn to case B. In this case, the mixing of the CBSM operators into the dipoles through γ B7 and γ B8 arises at one-loop order, but there is no mixing into the semileptonic coefficient C 9V . To restore the form (2.11) we should now rescale the electromagnetic and chromomagnetic dipoles as This ensures that γ B7,8 begin at O(α s ); at the same time it makes the blocks γ p7 and γ p8 vanish at O(α s ). As was the situation for the semileptonic coefficients in case A, now the dipole coefficients are formally of order O(1/α s ). Note that this makes the Standard Model contribution formally subleading. This is reflected in the fact that we will find below that B(B → X s γ) poses extremely stringent constraints on C B , corresponding to BSM scales of tens of TeV. The leading-order dipole coefficients are now scheme-independent and are not complemented by finite leading-order contributions from the four-quark operators at leading order, i.e. there is no need to define effective dipole coefficients.
Since C A and C B do not mix under renormalization, and the mixing into the penguin and dipole coefficients is simply additive, there is no problem patching together the solutions for both cases into an unambiguous and scheme-independent evolution operator. The part of the evolution matrix relevant for our phenomenology corresponds to the subset of coefficients (2.25)

Remarks on computed and uncomputed ADM elements
As part of our earlier work [1], we made the first calculation of the mixing of the BSM operators Q c 3,4 into the operators P 3−6 , Q 7γ , and Q 9V . The mixing into the photon penguin Q 7γ is the most technically challenging, arising from two-loop diagrams. We did not compute the two-loop mixing of Q c 3,4 into the gluon penguin Q 8g , and neglect this mixing in our numerical results in both the previous work and this article -the corresponding result for the mixing of the SM operators Q c 1,2 is known to produce only a very small effect on the contribution of C 2 (M W ) to C eff 7γ (µ b ), and we expect a similarly small effect from the results we neglect. Further details can be found in the Appendix of [1].

Observables
In this section we collect a set of observables that are very sensitive to (cb)(sc) operators and allow us thus to constrain the possible size of new b → ccs contributions: the dominant weak annihilation contribution to the B s lifetime, τ (B s ), is given by a (cb)(sc) transition, which is also the leading term to the mixing induced decay rate difference of neutral B s mesons, ∆Γ s . Taking a (cb)(sc) operator and closing the charm quarks to a loop we get large penguin contributions, that sizeably affect b → sγ and b → s decays. Finally the gold-plated mode for CP-violation, B d → J/ψK S is triggered by a tree-level b → ccs decay. Below we determine the dependence of all these observables on the new four quark operators.
The total decay width of the B s meson can be expressed as with the transition operator According to the Heavy Quark Expansion (HQE) (see [36] for a review and early references) the transition operator can be expanded in inverse powers of the heavy b quark mass -each term in the expansion contains perturbative Wilson coefficients and non-perturbative matrix elements of ∆B = 0 operators. We will investigate the precisely measured lifetime ratio The SM value of τ (B s )/τ (B d ) will be taken from [37], which uses perturbative input from [38,39]. The leading contribution to this lifetime ratio is given by the third order in the HQE (see e.g. [36]) and the dominant contribution at this order to the B s decay rate is given by ∆B The ratio of charm and bottom quark mass is denoted by To avoid a double counting of the SM contribution due to Q c 1 and Q c 2 , we subtract explicitly the SM contributions proportional to C c,SM i (C c,SM j ) * . The terms Γ(i, j) are symmetric (Γ(i, j) = Γ(j, i)) and they can be further split up into contributions of eight different ∆B = 0 fourquark operators.
with XY = LL, LR. The matrix elements of these operators are parameterised as with the decay constant f Bs . The bag parameters B1,2 have been determined in [37], while for the remaining bag parameters we will use vacuum insertion approximation: For the individual contributions we get (where we define for the sake of brevity) (3.10) (3.12) All remaining terms can be extracted from those given via the following rules: (3.14) The interested reader can download a Mathematica program containing the full algebraic expressions from the arXiv version of this article . where Γ s 12 denotes the absorptive part of the mixing diagrams and M s 12 the dispersive part. The relative phase of these contributions is defined as Similar to the case of the total decay rate the off diagonal matrix element of the absorptive part of the B s −B s mixing matrix, Γ s 12 can be expressed in terms of double insertion of the effective Hamiltonian This matrix element can be split up into a SM contribution and a BSM contribution due to the new b → ccs transitions The numerical value of the SM part is based on [41][42][43][44], the BSM part is further decomposed as The arising four quark operators can be parametrised as which matches the definitions in [40] for Q andQ s and [45][46][47] for Q 4 and Q 5 ; these bag parameters have been recently determined in [45][46][47]. The coefficient functions read and The interested reader can download a Mathematica program containing the full algebraic expressions from the arXiv version of this article. At this point we would like to mention that neglecting the CKM structure λ u (as advocated in Section 2) might result in misleading conclusions for the semileptonic CP asymmetry -depending on the experimental precision of the corresponding measurement. The mixing observables can also be determined via the following relations (3.27) Using the notation in [40] and CKM unitarity (λ u + λ c + λ t = 0) we can further simplify (3.28) The real coefficients show a clear hierarchy c a b and the small ratio of CKM elements λ u /λ t has an imaginary part. Thus in the SM the value of ∆Γ s /∆M s is well approximated by the term proportional to c, while a s sl is well approximated by the term proportional to a. We see from (3.28) that the approximation λ u = 0 yields a vanishing semileptonic CP asymmetry. Since the current experimental uncertainty of a s sl is about a factor of 130 larger than the SM central value of this quantity, neglecting λ u gives reasonable results. Moreover complex Wilson coefficients can have large effects in the semileptonic CP asymmetry by creating an imaginary part in the coefficient c, that does not suffer from the CKM suppression due to the factor λ u /λ t .

The radiative decay B → X s γ
It is well known that weak radiative B meson decays are sensitive to BSM physics. The Standard-Model prediction of the branching ratio for B(B → X s γ) SM = (3.36±0.23)×10 −4 [48] is in good agreement with the current experimental average of B(B → X s γ) exp = (3.32 ± 0.15) × 10 −4 [49]. In accordance with Heavy Quark Effective Theory (HQET) we may express the inclusive decay rate for a B meson into a charmless hadron and a photon as Here the non-perturbative term δ np , for E γ > E 0 with the lower cut off of the photon energy E 0 = 1.6 GeV, is estimated to be at the (3 ± 5)% level [50,51]. Following the approach of [52,53] the branching ratio B(B → X s γ) can be expressed as where P (E 0 ) and N (E 0 ) denote, respectively, the leading-power perturbative contribution and non-perturbative corrections, and C is defined as We neglect BSM corrections to the non-perturbative part and split the perturbative term P (E 0 ) into an SM part and a BSM part, We similarly split the branching ratio, To zeroth order in α s and neglecting the strange quark mass, we have which follows from substituting C eff 7γ → C eff,SM 7γ +C eff 7γ in the SM expression. The barred coefficients are defined as where ∆C c x,y (µ) = 3∆C c x (µ) + ∆C c y (µ) and with a(z) = |z − 1| arctan 1 √ z−1 and z = 4m 2 c /q 2 . That is to say, in addition to the BSM corrections to the Wilson coefficients C ( ) eff 7γ , which arise from large logarithms in the charm loop and are included through leading-order RG evolution as described in Section 2.2, we also include the remainder of the charm loop (see left diagram in Figure 1). For the SM contribution and BSM contributions to C scheme-dependent. 6 However, the q 2 dependence itself is a leading effect, and in any case is a qualitatively new feature, which is why we present it here. For the coefficients C

Rare b → s decays
The rare decays B → K ( * ) + − are also calculable to leading order in the heavy-quark expansion [54]. Through zeroth order in α s , they receive CBSM contributions both through C ( ) eff 7γ (q 2 , µ) defined in the previous subsection and through further coefficients C ( )eff 9 (q 2 , µ) (see Figure 1 (right) for the contributing Feynman diagram). The CBSM contributions to the latter have the form This q 2 -dependent contribution to the rare b → s decays was a novel feature in our previous work [1], where we showed that, if the new physics scale is low enough, it can be potentially observable. This contrasts with statements elsewhere in literature where such a q 2 dependence was claimed to be an unambiguous criterion that the anomalies 6 This scheme dependence will cancel against a corresponding scheme-dependence in the (uncomputed) NLO correction to C eff ( ) 7 , which involves two-loop mixing as well as the UV contribution to C ( ) 7γ in the B → K * + − angular distribution are of a hadronic origin. In the present work, we stick to new physics scales above the weak scale (as in the "high-scale scenario" of [1]). In this situation, C 9V is formally O(1/α s ) already in the SM, although about half the numerical value originates from the O(1), formally subleading, W t loop. Our q 2dependent contributions are formally at the same order but turn out to be a numerically small correction. The q 2 -independent contribution to C ( ) 9V , however, can be dramatic. In particular, when ∆C c 1 or ∆C c 3 are nonvanishing, it is easy to generate an O(1) shift to C 9V , while satisfying other constraints, as already stressed in [1]. Similarly, ∆C c 1,3 strongly mix into C 9V , which will allow us to set stringent constraints on these coefficients in Section 4 below.
The scheme dependence of these effective coefficients enters through the h and y functions -in general the constant terms will vary depending on the choice of scheme. In principle there is also a model-dependent contribution arising from the matching of the UV theory, which in this paper we choose to ignore to keep our results model-independent.
When we consider the phenomenology in Section 4, we directly use C eff,BSM 9 (q 2 , µ) andC eff 7γ (q 2 , µ) as "observables" and show the constraints they impose upon the coefficients ∆C c i . Again, for the numerics in this work we use m c = m c,pole . We will take q 2 = 0 for C ( ) 7γ , but will consider C eff,BSM 9 at q 2 = 5 GeV 2 . This is because the former is primarily constrained from polarisation observables in radiative decay and the low-q 2 end of the B → K * + − dilepton mass distribution, while C eff,BSM 9 , if present, will become important away from the endpoint. We also note that C eff 9 and C eff 7 are negligible in the SM.

The hadronic decay B d → J/ψK S
The b → ccs operators listed in Section 2 trigger the neutral B meson decay B d → J/ψK S . As a peculiar feature, both B d andB d can decay into the J/ψK S final state, giving rise to a time-dependent CP-asymmetry via mixing-decay interference. The amplitudes of these colour-suppressed tree-level decays read

41)
where the hadronic matrix elements J/ψK S |Q i |B d contain the (CP-even) strong rescattering phases and the minus sign arises from η CP (J/ψK S ) = −1. Their determination is a non-perturbative problem which at present cannot be solved from first principles. We will develop a strategy to extract partial information on these matrix elements, jointly with information on the Wilson coefficients, from experiment in the following. Let us express the branching ratio and time-dependent CP-asymmetry in terms of the Wilson coefficients and hadronic matrix elements. Defining the branching ratio, obtained from the calculated B 0 d → J/ψK 0 d decay rate along with using the measured B 0 d lifetime (for a similar approach see for example [55]) reads where p c = 1.683 GeV [56] is the momentum of the final-state particles in the rest frame of the B d meson.
In writing (3.44) we have omitted penguin contributions, which in what follows will be negligible compared to the uncertainties stemming from the hadronic matrix elements, due to their small Wilson coefficients.
For the time-dependent CP-asymmetry, we neglect the tiny decay rate difference ∆Γ d , such that it takes the form where ∆M d is the precisely measured mass difference in the B d system. Defining also the mixing-induced asymmetry S J/ψK S and the direct CP asymmetry C J/ψK S are expressed as In writing (3.45), we have neglected CP violation in the Kaon system and, for the final expression, taken q B /p B = V * tb V td V tb V * td = e −2iβ+O(λ 4 ) , which amounts to omitting CP-violation in mixing (λ ≈ 0.22 is the Wolfenstein parameter). For reconstruction of (3.44), (3.45) and (3.46) see for example [40], noting that in their notation S JψK S = −A mix CP and C JψK S = A dir CP .
If the Wilson coefficients are real (recall their phase convention is fixed by (2.2)), the unknown hadronic matrix elements cancel out in (3.45) and, taking into account the tiny measured value of C J/ψK S , one obtains S J/ψK S ≈ Im(λ J/ψK S ) ≈ sin(2β) to percent-level accuracy, as in the SM. Once the Wilson coefficients become complex, this is no longer true and the S-and C-parameters become theoretically very uncertain. However, together with the branching fraction, they still comprise three observables, allowing to jointly determine up to three parameters. Consider now a scenario where new physics only affects C c 1 or C c 2 . In this case, the only nonperturbative input to the three observables is | Q c 1 | and r 21 , comprising three real parameters in total. Hence it is sufficient to have theoretical control over one of these parameters in order to obtain a constraint in the complex C c 1 (or C c 2 ) plane, generally a band.

To proceed, let us use a Fierz-transformed basis
The hadronic matrix element O c 1 = Q c 1 factorizes naively in the limit of a large number of colours. More precisely, where in the last expression we have taken 1/N 2 c = 1/9 and combined errors in quadrature. This provides the required theoretical constraint. The three observables can now be used to jointly constrain the complex Wilson coefficient and the complex matrix element ratio r 21 . In a global fit, this amounts to determining r 21 from data. This is the strategy which will be followed in our phenomenology below.
The fact that r 21 can (in the limited scenarios described) be determined from data motivates us to review expectations for the other operator matrix elements. Let us write them as the sum of a naive-factorization result and a correction term, The constant r 0 amounts to the naive-factorization term, r 0 = 1 for O c 1 and O c 3 and r 0 = 1/3 for O c 2 and O c 4 . QCD factorization in the heavy quark limit [21] implies that the coefficients r n , n ≥ 1, are perturbatively calculable but that there is also a power correction "suppressed" by only Λ QCD /(α s m c ). (See also [57,58].) As this power correction remains, at present, incalculable, one does not expect predictivity, beyond the large-N c argument already presented. From a large-N c perspective, the corrections to naive factorization are 1/N 2 c -suppressed for both O 1 and O 3 , but are unsuppressed for O 2 , O 4 . The importance of contributions that do not naively factorize (often called 'non-factorizable') for these two operators becomes clearest by rewriting The matrix elements of the colour-octet operators T c 1,3 vanish altogether in naive factorization, i.e. r 0 = 0 for them, and they appear with large coefficients in (3.54).
Therefore the naive-factorization values should be expected to receive large corrections, whereas r 31 = 1 + O(1/N 2 C ). The situation is further aggravated in the decay amplitude, which is proportional to the combination C c 1 (µ) + r 21 C c 2 (µ), the so-called colour-suppressed tree amplitude which entails a severe cancellation in the SM case.
Quantitative estimates of r 21 beyond naive factorization have been obtained using QCD factorisation and/or light-cone sum rules, [57][58][59]. In these approaches, the amplitude is usually parameterized as where N is a normalization and the scheme-and scale-independent combination a 2 is related to r 12 as ) corrections, which as we have said provide an (additional) uncertainty. By comparison, the LCSR computation in [59] finds values closer to NF. These numbers can be contrasted to the 'experimental' value r 21 (m b ) ∼ 0.46. It is important, however, to note that this assumes the absence of NP in the Wilson coefficients C c 1 , C c 2 (and in addition vanishing relative phase between the hadronic matrix elements.) As we will see later, this assumption is not implied by current experimental data when NP in b → ccs transitions is allowed.

Numerical Inputs
In this section, we describe all the numerical inputs that are used in our work, along with the experimental results and their corresponding uncertainties. We break these down into a set of fundamental inputs that are common to all our different observables, and then some specific input parameters for the individual observables. Recall that quark masses are in the MS scheme unless indicated otherwise.

Common inputs
We show in Table 2 input parameters that are common to all our theoretical predictions. These inputs are taken from the Particle Data Group (PDG) [56] and the CKMfitter group [29] (similar results for the CKM elements are also available from the UTfit collaboration [30]).

Lifetime ratio
The lifetimes of the B d and the B s meson are measured quite precisely. The Heavy Flavor Averaging Group (HFLAV) quotes [49]: The non-perturbative matrix elements of the ∆B = 0 operators were determined recently in [37] with Heavy Quark Effective Field Theory (HQET) sum rules, using the following notation: with the coefficients We use the following numerical values for the decay constant from the Flavour Lattice Averaging Group (FLAG) [62] and bag parameters from [37]: The tilded bag factors defined in Section 3.1 are expressed in terms of these as follows: For the SM value of the lifetime ratio we use the prediction from [37] τ (B s ) τ (B d ) = 0.9994 ± 0.0025 .

B → X s γ
For the inclusive radiative B s meson decay we use the experimental average obtained by HFLAV [49] and the theoretical prediction from Misiak et al. [48]: 14) The semileptonic branching ratio and the ratio C as defined in Section 3.3 are taken from [52]: The SM contribution to B → X s γ which interferes with our BSM contribution is given by [63] C eff,SM

Rare decays from BSM operators
In light of the recent anomalous measurements of b → s decays by LHCb, in particular the R K ( * ) results [5,64,65] there has been considerable work on fitting the semileptonic and radiative Wilson coefficients to data. While R K ( * ) are indicative of a lepton-flavournon-universal NP effect, a UV completion of the EFT will, in general, also include a leptonflavour-universal effect. Such a combined scenario has been shown to be consistent with (and even mildly preferred by) the data in [16] (see also [17][18][19]). In [1], we have studied the possible C 9V effects generated by C c 1 . . . C c 4 in detail. In the presence of the operators C c i involving right-handed strange quarks, the CBSM scenario produces an effect in C 9V as well, associated with the right-handed semileptonic operator Q 9V . We will treat C eff,BSM 9 (q 2 , µ) in (3.39) as a pseudo-observable and use the below value taken from the fit in [19] to constrain our model C eff,BSM,exp where the theoretical uncertainties associated with exclusive rare B decay are included in this value by the authors [19] and are taken from [66][67][68] compatible with [69,70]. Similarly, rare and radiative b decays can receive CBSM contributions via the coefficient C 7γ of the right-handed dipole operator Q 7γ . We treatC eff 7γ (q 2 , µ) in (3.36) as a second pseudo-observable, with our experimental value taken from the fit in [71], (4.20) (Some more recent fits [18,72] give very similar results for theC eff,exp 7γ coefficient that would not change our numerics significantly.) Recall that we take q 2 = 5GeV 2 in C eff,BSM 9 and q 2 = 0 inC eff 7γ . We have identified C eff 9 = C eff,BSM 9 and C eff 7γ =C eff 7γ due to the negligible SM contributions. from [49], and the branching ratio from [56]. As part of our theoretical calculation of S J/ψK S and C J/ψK S , we use the most recent CKMfitter [29] value sin 2β = 0.738 +0.027 −0.030 , (4.23) where the experimental measurement is not included in their fit. We note there is a very slight tension between the HFLAV average and the CKMfitter result, at the level of ∼ 1.1 σ.
For the non-perturbative decay constant of the J/ψ resonance, we take the value from the phenomenological study in [73] f J/ψ = (407 ± 6) MeV . (4.24) This value agrees well with the lattice determinations in [73][74][75]. The form factor can be determined via LCSR [76] or extrapolated from lattice simulations [77,78]. Both approaches have similar uncertainties, the LCSR values are slightly larger than the lattice results. Since the LCSR value is a direct evaluation at q 2 = M 2 J/ψ , we have chosen that as our input for the evaluation of the O c 1 matrix element. However we note for the reader that using an average instead would give very similar results, albeit with a slightly smaller central value and error, and therefore our choice is conservative. To determine confidence intervals for individual Wilson coefficients, we switch on BSM contributions in one Wilson coefficient at a time and set all others to their SM values. We construct a χ 2 test statistic from the experimental measurements of our chosen observables and our theoretical predictions, combining the experimental and theoretical errors in quadrature 7 :

Constraints on ∆C
where the index i runs over B(B → X s γ), ∆Γ s and τ (Bs) τ (B d ) and in addition our "pseudoobservables"C eff 7γ (q 2 , µ) and C eff,BSM 9 (q 2 , µ). The 1σ intervals implied by individual observables are displayed in Tables 3 and 4. To obtain combined constraints (Table 5), we sum up the individual χ 2 . In all cases we normalise to the best fit value by subtracting the relevant χ 2 value at the minimum, ∆χ 2 = χ 2 − χ 2 min . Considering ∆C 5 − ∆C 10 in the first column of Table 3, there are best fit ranges which correspond to those that pass through the SM point and those that do not. This can be understood by considering the functional form with whichC eff 7γ enters (3.34) and the impact that larger contributions from coefficients in C eff, BSM 7γ have upon reducing the parameter space allowed by radiative decay. The second column of Table 3 does not for any coefficient include the SM point, and this is simply due to the current disagreement between measurement and theory for the lifetime ratio. Column 3 containing ranges accommodated by ∆Γ s always includes the SM point. In Table 4 we show 1σ ranges for the primed coefficients, accommodated by pseudo observables C eff,BSM 9 andC eff 7γ . In the second column of Table 5 we show the 1σ allowed ranges for ∆C 5 − ∆C 10 which follow from combining constraints of all observables. The 1σ ranges for the primed coefficients ∆C 1 −∆C 10 however  . In the last two columns of Table 5 we re-express the combined bounds in terms of the scale of new physics Λ NP , defined through The lower bound on the NP scale corresponding to the negative boundary of the 1σ interval for ∆C i is denoted as Λ − , that corresponding to the positive boundary as Λ + . For ∆C 10 and ∆C 9 the 1σ interval only contains positive values and only Λ + is given, corresponding to the upper boundary. For ∆C 10 , the 1σ region is composed of two intervals, and Λ − and Λ + correspond to the smallest and largest Wilson coefficient values (−0.08 and 0.05, respectively). We see that with our definition of Λ NP , which is agnostic to the details of BSM  physics generating the ∆C c i , our observables can provide sensitivity to scales higher than 20 TeV in scenarios involving the tensor Wilson coefficient ∆C ( ) 9 . Other Wilson coefficients probe somewhat lower scales, but always at least 3 TeV.
For scenarios in which we consider NP in pairs of Wilson coefficients, setting all others to their SM values; the pattern of constraints may be divided into three categories: (i) Coefficients ∆C 5 − ∆C 10 : As explained in Section 2.2, the mixing of operators Q c 5 , ..., Q c 10 with Q 7γ first at one-loop gives rise to strong RG effects which enter C eff 7γ and result in the dominant constraint for such scenarios coming from radiative decay. As the coefficientC eff 7γ enters (3.34) both quadratically as well as linearly, it receives a contribution from C eff,SM 7γ . This combination results in a much narrower 1 sigma region, as is shown in all panels of Figure 2 as the blue shaded area. In the first and third panel, the presence of another purple band corresponds to contours wherē . In terms of the lifetime ratio, shown in the green shaded area, the contours slightly miss the SM point due to the current 1.4σ discrepancy between theory and experiment. For the B s −B s width difference the scenarios consisting of coefficients of operators with left and right handed vector currents and coefficients of tensor operators are the most restrictive. In all cases scenarios between even numbered coefficients are favoured owing to the 1 Nc suppression which always accompanies colour singlet operators in the calculations. Here the central value used is the best fit points for C 9V and C 7γ acquired from global fits to angular observables in [19], [71]. These are to be compared with the high scale scenarios considered already in [1] and discussed in Section 4.3. For combinations of these coefficients it is found that the strongest constraint comes from the experimental fits of angular observables to C 9V . This is due to a strong dependence of C BSM 9V (m b ) upon ∆C 1 (M W ) − ∆C 4 (M W ) which results in closely spaced contours (red dashes). The inclusive radiative decay branching ratio constraint (blue shading) is shown for extra information, in addition to that given by the contours ofC eff 7γ at fitted value ofC eff,exp 7γ (which includes inclusive branching ratio data). For these primed coefficients the radiative decay constraint is very much weaker than in the previous case due to the mixing of operators Q c 1 , ..., Q c 4 with Q 7γ occurring first at two-loop, and in addition, to the primed coefficientC eff 7γ only entering (3.34) Figure  4 and are the right handed counter parts in one to one correspondence with those of Figure 2. As in (i) the one-loop mixing under renormalization of Q c 5 , ..., Q c 10 with Q 7γ results in a stronger dependence ofC eff 7γ upon ∆C 5 − ∆C 10 and this results in strict constraints upon combinations of these coefficients. All plots show very similar constraints from the lifetime ratio and the width difference, although differences are more pronounced in ∆Γ s due to primed and unprimed coefficients not mixing in the theoretical prediction for Γ cc 12 and hence there is no linear contribution from SM  For the case of the C c 1−4 , the phenomenology for a scenario with purely real coefficients was covered in our previous paper [1] -we briefly recap our conclusions from that work, before expanding to a scenario with complex coefficients and the constraints that arise from B d → J/ψK S decays. Complex NP in C c 1 and C c 2 is studied in the recent work [79] where they compute constraints arising from τ (B s )/τ (B d ) and B(B → X s γ) as part of a sophisticated global fit. We performed a mutual check with the authors of that work, and found full agreement on the constraints arising from the previously mentioned observables (i.e. the blue and green bands in Figure 5 below).
In [1], we studied NP confined to the first four operators of the full basis (2.3), as these operators give a large contribution to a (lepton flavour universal) shift in the C 9V coefficient (since the RG running coefficients are O(1/α s ) in the logarithmic counting), while only being constrained by the radiative decays through two-loop RG mixing. In our study we found that while the SM is consistent with the lifetime ratio, width difference and radiative decay observables, there is also room for a shift in our (cb)(sc) Wilson coefficients without disagreement with data -see Fig. 3 in [1]. For instance, a shift to the C c 3 coefficient alone of order 0.2 could produce a shift of C 9V of order −1, and such an NP contribution is in fact slightly favoured with respect to the Standard Model, as it lessens the small tension present in τ (B s )/τ (B d ). NP in two coefficients simultaneously can also be accommodated, such as in the pair (∆C 2 , ∆C 4 ) = (−0.1, 0.3) which generates an O(1) contribution to C 9V (albeit with no change relative to the SM in the lifetime ratio). In light of the fact that current data supported a possible NP contribution in several different Wilson coefficients, the natural question was how to distinguish between these scenarios. We showed that an improvement in the future precision of both the lifetime ratio and the width difference ∆Γ s could identify the particular realisation of charming BSM physics in nature.
As introduced in Section 3.5, NP in (cb)(sc) operators can alter the B d → J/ψK S decay rate, as well as the related CP asymmetries S J/ψK S and C J/ψK S in the case of complex Wilson coefficients. In order to predict these three observables, hadronic operator matrix elements J/ψK S |O c i |B must be evaluated. In the SM, the NF approximation for them does not give good agreement data, and it is widely assumed that deviations from NF can bring the prediction in line with experiment; as reviewed in Section 3.5, large deviations from naive factorization are theoretically to be expected for O c 2 but not O c 1 . However, in this section we jointly consider the three observables, to constrain either ∆C 1 or ∆C 2 together with the uncertain matrix element ratio r 21 defined in Section 3.5, with the matrix element O 1 , for which violations of NF are expected to be small, taken in the range described there. Constructing a χ 2 test statistic out of the three observables and the O 1 range and profiling it over O 1 and r 21 results in a constraint in the complex ∆C 1,2 planes.
For the complex ∆C 1 plane, the resulting constraint is shown in Figure 5 (left) as red bands, implying a remarkably powerful constraint from this hadronic decay: the imaginary part must be close to either zero or ±0.2. In effect, we have determined the uncertain matrix element ratio r 21 from data. We show in the figure the 1, 2, and 3 σ regions for the B d → J/ψK S constraints. The Im C 1 = 0 band is very narrow, and has no 1σ region because of the slight discrepancy for sin2β mentioned in Section 4.1.6. In the same figure we overlay the constraints from a s sl (yellow), ∆Γ s (orange), and τ (B s )/τ (B d ) (green). We see that the lifetime ratio is the most constraining among the three, while the other two rule out some regions with larger real and/or imaginary NP contributions (B → X s γ provides no visible constraints at this scale). Combined with the constraint from B → J/ψK S , only a few small regions in the ∆C 1 plane are allowed.
In fact, the χ 2 -profiling allows, at each point in the plane, to simultaneously determine the value of r 21 that gives the best fit. For reasons of computational simplicity, we carry this out only along a circle in the middle of the (green) lifetime ratio ring, where the experimental central value of τ (B s )/τ (B d ) is obtained. Along this ring, the (complex) value of r 21 varies substantially. Along the two short black segments it is in excellent agreement with naive factorization. Surprisingly, these two segments happen to lie in two of the small regions in the plane that are allowed by all constraints. We have no explanation for this curious fact. But it certainly demonstrates that experimental data in B → J/ψK S does not imply large violations of naive factorization, in contrast to a widely held belief. (Similar results are found when varying the radius of the circle within the lifetime band, corresponding to different fixed values of τ (B s )/τ (B d ).) The results of repeating the study for complex ∆C 2 are displayed in Figure 5 (right). In this case, B d → J/ψK S data allows a band centred on the real, as well as a second band with negative real and imaginary parts that however is ruled out by the other constraints; among those, the most stringent constraints now come from the lifetime difference ∆Γ s and the radiative decay B → X s γ. Extending beyond the SM operators Q c 1,2 to NP in Q c 3,4 , the theoretical prediction for B d → J/ψK S requires more non-perturbative inputs in the form of r 31 , r 41 , alongside r 21 (as there is still a SM contribution to C c 2 ). Attempting to fit these from the B d → J/ψK S data does not bring any further insight, as there is sufficient freedom to always explain the measurements.

Conclusions
In this paper we have made a thorough study of the possible effects of new physics arising in tree-level b → ccs decays. This partonic decay mode contributes to a wide variety of different observables. Among them, the branching ratio for radiative B meson decay B(B → X s γ), the B meson lifetime ratio τ (B s )/τ (B d ), and the B s mixing observables ∆Γ s and a s sl stand out: They are inclusive decay modes that are well measured experimentally, and are theoretically controlled through the HQE. In addition, we have shown that, in the CP-violating case, new constraints which are only midlly affected by theoretical uncertainties arise from exclusive B d → J/ψK S observables. Taken together, effects in this set of observables are correlated in our "Charming BSM" scenario, and the observables provide very complementary constraints on it.
The space of NP contributing to b → ccs decays is spanned by 20 operators, defined in (2.3). We have calculated the contribution from the full basis to all of our observables; the most complex result being that obtained for mixing and the lifetime ratio for which we used tools from the Mathematica package FeynCalc [80,81]. Our full results are given in (3.20)-(3.24) and (3.4)-(3.14) (these results are also available as ancillary Mathematica files with arXiv version of this article). We have further calculated the renormalization group evolution for our basis. In the SM case we used known results for ADM entries available in the literature [25,35] and for the mixing of our operators as described in Section 2.2, elements from [33], our previous work [1], as well as those elements extracted from our b → s results either directly or through substitution of relevant colour factors. Our results are summarised in the full evolution matrix given in (2.25). An explicit recipe for making use of our results to place constraints on an arbitrary NP model is as follows: 1. Match the chosen BSM model onto our effective Hamiltonian (given in (2.2)) at the scale M W .
2. Use the RG evolution matrix (given in (2.25)) to run the effective coefficients down to the scale m b . We have extended our earlier results [1] in two important ways. Firstly, we complement our previous analysis of real new physics in Q c 1−4 by studying the possibility of CP-violating NP, and introduced observables from B d → J/ψK S as constraints. Secondly, we studied constraints from our main set of observables on the Wilson coefficients C c( ) 5−10 and in analogy with our study of C c 1−4 we have constrained from global fits to the "wrong chirality" coefficients C 9V and C 7γ , possible BSM effects in C c( ) 1−4 . When considering the introduction of new weak CP violating phases to SM coefficients, we have used the sine and cosine coefficients S J/ψK S and C J/ψK S of the time dependent CP asymmetry, alongside the branching ratio B(B d → J/ψK S ) to constrain the parameter space. By treating the uncertain hadronic matrix element ratio r 21 = O c 2 / O c 1 as a free parameter and profiling over it, we effectively determine it from data, reducing the required theoretical input to the colour-allowed matrix element O c 1 . We assume this to be close to its value in naive factorization, as is expected from large-N C counting. In this way we obtained constraints relying only mildly on theory, even though we are dealing with nonleptonic exclusive decays.
For NP in ∆C 1 , our result (shown in Figure 5) turns out to be very interesting. Firstly, the combination of inclusive and B d → J/ψK S constraints only leaves a few small regions of the complex C 1 plane where agreement is obtained between the full compliment of observables and their respective experimental averages, including some where C 1 has an imaginary part of close to ±0.2. Secondly, whereas naive factorisation is not expected to well describe class II colour suppressed decays, we find that in one of the allowed regions with complex C 1 , r 21 happens to be close to its NF prediction. This requires a small negative imaginary BSM shift ∼ −0.2i to C 1 . When considering ∆C 2 , we again observe a strong complementarity of the constraints. A broad band centred on real shifts is compatible with B d → J/ψK S data, as well as a diagonal region with negative real and imaginary parts. Unfortunately the other constraints we consider have no clear region of overlap where all the predictions can be brought into agreement with data, due to a mild tension between the lifetime ratio on the one hand and radiative decay and the B s width difference on the other hand. When considering ∆C 3 and ∆C 4 , there are too many non-perturbative parameters in play to obtain constraints from B d → J/ψK S . Turning now to our results for ∆C ( ) 5−10 and ∆C 1−4 , we group them into three categories exhibiting similar behaviour. Consider ∆C 1−4 , it was found that the strongest constraint upon these coefficients comes indirectly from angular observables through our displayed contours of constant C eff,BSM 9 at the fitted values of C eff,BSM,exp 9 ± 1σ, and to a lesser degree, contours of constantC eff 7γ at the fitted values ofC eff,exp 7γ ± 1σ. In all panels we would expect viable values of pairs of BSM coefficients to lie within the region bounded by these contours. We obtain no strong constraint from radiative decay for these coefficients as the small 2 loop mixing of C eff 7γ with this set of coefficients in addition to the purely quadratic dependence of the branching fraction uponC eff 7γ , leads to relaxed bounds. Of all the scenarios considered, we find that pairs involving ∆C 1 − ∆C 3 and ∆C 1 − ∆C 4 stand out as scenarios where agreement with all data can be found. Considering now ∆C 5−10 , in contrast to the above case, the mixing of C eff 7γ with ∆C 5−10 occurs at 1-loop and results in these coefficients being very highly constrained by radiative decay, and this indicates that models involving combinations of these coefficients are disfavoured by our study. Finally, we consider the coefficients ∆C 5−10 . These are constrained in a similar fashion to their unprimed counterparts by each of our observables except the radiative decay branching ratio, which we replace by contours ofC eff 7γ =C eff,exp 7γ ± 1σ. Graphical representations of the allowed parameter space are shown in Figures 2-4.
As a step towards converting our constraints into statements on more definite NP models, we considered what the equivalent NP scale Λ NP we are probing when we place limits on our Wilson coefficients, and our results were shown in Table 5. The tensor operators Q ( )c 9 are sensitive to the highest scales, with the best fit to those coefficients corresponding to scales in excess of 10 TeV. All our operators probe scales above 2 TeV, showing how our choice of observables can complement direct LHC searches for NP effects.

A Explicit expressions for anomalous dimensions matrices
Note that as mentioned in Section 2.2.1, the two-loop mixing of Q c 3,4 into the gluon penguin Q 8g has not been calculated -this corresponds to the zeros in the third and fourth elements of γ