NLO corrections to $h\to b\bar b$ decay in SMEFT

We calculate the full set of next-to-leading order (NLO) corrections to $h\to b\bar{b}$ decay in the dimension-6 Standard Model Effective Field Theory (SMEFT). Our calculation forms the basis for precision studies of this decay mode in effective field theory, providing analytic and numerical results for contributions of the 45 dimension-6 operators appearing at NLO. On the technical side, we discuss several complications in NLO SMEFT computations which have not yet been addressed in the literature. These include subtleties in Higgs-$Z$ mixing, electric charge renormalization, and especially the treatment of tadpoles in SMEFT. In particular, we highlight the role of decoupling relations in eliminating potentially large tadpole corrections to the decay rate in hybrid renormalization schemes which employ the $\overline{\hbox{MS}}$~scheme for some Standard Model parameters (such as the $b$-quark mass and electric charge) and the on-shell scheme for others.


Introduction
While the discovery of the Higgs boson has been a triumph for the Standard Model (SM) of particle physics [1][2][3], the consistency of its properties, as currently measured, with those predicted by the SM (see the experimental analyses in [4][5][6][7][8] for example) has left few hints of new physics. An important property of the Higgs boson is its decay rate into b-quarks. Despite being the largest branching fraction of the Higgs, the process h → bb has only recently been observed by the ATLAS and CMS collaborations [9,10]. Considering the relative infancy of the Higgs measurements so far in the LHC program, as well as the prospect of future e + e − colliders for such studies [11,12], the possibility of uncovering new physics in the Higgs sector remains open. As such, the need for accurate theoretical predictions in order to correctly identify and parametrize any new physics which could be observed is paramount.
In the absence of the direct discovery of a new particle, one possible avenue along which to search for new physics is through the use of the Standard Model Effective Field Theory (SMEFT). In this approach the SM Lagrangian is supplemented with operators of mass dimension greater than four, each with its own Wilson coefficient. Provided the new physics is associated with a scale Λ NP which is much greater than the electroweak symmetry breaking (EWSB) scale and decouples [13], then its effect on processes at low energy is captured through non-zero values of these Wilson coefficients. This allows for a model independent approach in attempts to identify new physics: one calculates cross sections and decay rates within SMEFT and then fits the Wilson coefficients to data in order to extract limits or signals of new physics.
The SMEFT operators which can be written down at a given mass dimension are constructed out of SM fields and respect the usual SM gauge and Lorentz symmetries. A minimal basis of operators (though not unique) can be constructed by using the SM equations of motion [14] and techniques to quantify the minimal number of operators and their field content which appear at each mass dimension have already been developed [15][16][17]. At dimension-5 there exists only a single, lepton number violating operator, whose Wilson coefficient is heavily suppressed. On the other hand, at dimension-6 there are 59 independent operators for one generation of fermions excluding baryon number violating operators [18,19], giving a wide space in which to explore possible consequences for phenomenology.
Recently the inclusion of dimension-6 operators in NLO perturbative calculations has emerged. Some general features of these calculations have been described in e.g. [20][21][22][23], and the full 59 × 59 anomalous dimension matrix for the Wilson coefficients needed to perform a leading-logarithmic calculation has been calculated in [24][25][26]. At the moment, however, there is no automated tool to produce general NLO SMEFT predictions so these calculations are performed on a process-by-process basis. Because of the increased complexity of these calculations, results are available only for a handful of processes, and often contain a limited number of operators or are restricted to a particular set of corrections. There are many NLO SMEFT calculations which involve a subset of operators [27][28][29][30][31][32][33][34], or are restricted to QCD corrections only [35][36][37][38][39][40][41][42][43][44][45]. A calculation of Higgs pair production at NNLO in QCD involving dimension-6 operators which contain the Higgs field has also been performed [46]. A small set of processes has been computed at NLO including all relevant operators in both the tree and loop level diagrams. These include lepton decay [47,48] and Higgs decay into vector bosons [49][50][51][52][53][54].
In this paper we obtain the full set of NLO corrections from dimension-6 operators to the decay rate h → bb within SMEFT. This builds upon our previous NLO SMEFT calculations of weak corrections in the large-m t limit or those related to four-fermion operators [55], and QCD corrections [56]. On the practical side, our calculation forms the basis for a precision analysis of Higgs decay into b-quarks within SMEFT. However, even apart from that, calculating the full set of NLO corrections reveals features of SMEFT beyond tree level which have not been fully addressed in the literature. For instance, one encounters technical subtleties in the renormalization procedure concerning electric charge renormalization and Higgs-Z and Higgs-neutral Goldstone mixing. Moreover, when combining electroweak and QCD corrections it is natural to introduce hybrid renormalization schemes where some parameters are defined in the MS scheme and some in the on-shell scheme. In that case one must pay careful attention to tadpole contributions, not only including them in the renormalization procedure in order to obtain gauge-independent results, but also finding a renormalization scheme where enhanced electroweak corrections related to them are absent. In this work we address tadpole renormalization using the "FJ tadpole scheme" [57], which is especially convenient when performing loop calculations with automated tools, and advocate the use of decoupling relations in building a renormalization scheme which allows us to combine QCD and electroweak corrections in an optimal way. The organization of this paper is as follows. After giving an outline of the NLO calculation as a whole in section 2, we describe in detail the renormalization procedure in section 3, including our treatment of tadpoles. We discuss sources of enhanced NLO contributions to the decay rate in section 4, and explain how a hybrid renormalization scheme based on decoupling relations for the MS definition of the b-quark mass and electric charge is useful when combining QCD and electroweak corrections. In section 5 we present numerical results and examine uncertainties related to scale choices, and then conclude in section 6. We provide some details on the rotation of the SMEFT Lagrangian to the mass basis relevant for our NLO calculation in appendix A, including a novel treatment of gauge fixing in SMEFT, and give selected analytic results for the decay rate in appendix B. While the full analytic results are too long to print, we give them in electronic form in the arXiv submission of this article.

Outline of the calculation
The dimension-6 SMEFT Lagrangian may be written as L = L (4) + L (6) ; where L (4) denotes the SM Lagrangian, and L (6) depends on the dimension-6 operators Q i . We adopt the "Warsaw basis" [19] for these operators, which are listed in table 3, and the naming convention of the Wilson coefficients C i follows that of the corresponding operators. We define the Wilson coefficients such that they inherently carry two inverse powers of the new physics scale, Λ NP .
In this paper we study the decay rate for h → bb to NLO in SMEFT. We can write the perturbative expansion of the decay rate up to NLO in the form where the superscripts (0) and (1) refer to the LO and NLO contribution in perturbation theory respectively. Each of these can be split up into SM (dimension-4) and dimension-6 contributions with the notation The double superscripts (i, j) refer to the dimension-i contribution at j-th order in perturbation theory. In this counting each term in Γ (6,j) contains exactly one Wilson coefficient of a dimension-6 operator. In other words, we allow at most one insertion of a dimension-6 operator in a given Feynman diagram and keep the interference term of the dimension-6 amplitude with the SM, but drop the square of dimension-6 amplitude, which is formally a dimension-8 effect at the level of the decay rate.
It is useful to divide the NLO correction from dimension-6 operators into three pieces according to and analogously for the SM result Γ (4,1) , which was calculated in [58]. The definition of the three pieces, and the extent to which the dimension-6 corrections have been calculated in the literature, is as follows. First, Γ g,γ contains all virtual and real emissions involving gluons and photons. The QCD portion of this object was calculated in [56]. Second, Γ t contains virtual weak corrections in the large-m t limit. These were calculated in the on-shell renormalization scheme in [55], where they scale as αm 2 t /M 2 W . Finally, the object Γ rem contains the remaining virtual electroweak corrections. The only results available for these remaining contributions are those from four-fermion operators obtained in [55].
The main goal of the present work is to obtain the full NLO correction in SMEFT. To do this, we must calculate the UV-renormalized virtual corrections to the LO decay rate, and add them together with real emission corrections containing a photon or gluon. We then evaluate to NLO the formula where dφ i is the i-body differential Lorentz invariant phase-space measure. The 2-and 3body terms involving emissions of gluons or photons contribute to Γ g,γ . These contain IR divergences, which we regularize by performing the loop integrations and phase-space integrals in d = 4 − 2 dimensions. Most of the corrections involving photons can be extracted from the QCD calculation [56]. The exception is real and virtual diagrams containing a hγZ vertex which has no analogue in QCD. Analytic results for Γ g,γ are given in appendix B. The most challenging part of the calculation is to obtain the UV-renormalized 2-body matrix element M (1) (h → bb), which is needed to determine Γ rem . We do this by evaluating the expression where the terms on the right-hand side are the bare one-loop and counterterm amplitudes, respectively. The exact form of the counterterm and bare amplitude depends on the set of independent parameters in terms of which the SMEFT Lagrangian in the mass basis is expressed, and also the scheme in which these parameters are renormalized, as discussed in more detail below. We choose the parameters to be where α = e 2 /(4π) and α s = g 2 s /(4π) are the electromagnetic fine-structure and strong coupling constants respectively, and m f are the fermion masses. We allow for non-vanishing thirdgeneration masses m b , m t , and m τ , but set first-and second-generation fermion masses to zero. We work with the numerical approximation of a diagonal CKM matrix V ij = diag(1, 1, 1), but do not necessarily impose Minimal Flavour Violation (MFV); further details on this point can be found in appendix A.4.
To perform the NLO calculation, we follow the procedure set out in [55]. We first express the SMEFT Langrangian in the mass basis, using the parameters in (2.7). There are a number of differences in this procedure compared to the SM, the most significant of which involve gauge fixing, which are described in appendix A. We then trade the bare input parameters for renormalized ones in order to construct an explicit expression for the counterterm amplitude in (2.6). Here again there are a number of subtleties compared to the SM, especially in the structure of tadpole contributions. The full details of the renormalization procedure are covered in section 3. Finally, we must identify and evaluate the large number of one-loop Feynman diagrams which contribute to the bare matrix elements and UV counterterms. We have automated the procedure by implementing the SMEFT Lagrangian in the mass basis, including ghosts, into FeynRules [59], and then using the resulting model file to generate the diagrams with FeynArts [60] and compute them with FormCalc [61]. We have also made use of Package-X [62] when extracting analytic expressions for loop integrals.
The NLO correction Γ (1) obtained in this way is quite lengthy. In fact, we obtain contributions from 45 different dimension-6 operators when full mass dependence of third-generation fermions is kept. We give the result in symbolic form in the computer files available with the electronic version of this submission. We have performed three main checks on these results. The first is that the UV poles in the bare and counterterm matrix elements cancel against each other, and the related fact that the decay rate is independent of the renormalization scale µ up to NLO. The second is that the IR poles appearing in the 2-and 3-body contributions to Γ (1) g,γ cancel against each other. Finally, we have verified the gauge independence of our results by performing all calculations in both unitary and Feynman gauge.

The renormalization procedure
In this section we lay out the renormalization procedure used in our calculation. We draw on the methods used in [55] to construct the one-loop counterterm in section 3.1, but must deal with technical complications not present in the partial NLO calculation in the on-shell scheme performed there. We point out subtleties with charge renormalization in section 3.1.1 and with Higgs-Z mixing in section 3.1.2, before moving on to discuss tadpole renormalization in section 3.2.

The one-loop counterterm
The form of the NLO counterterm follows directly from the LO decay amplitude. We write the LO decay amplitude as which we split up as where the superscripts (4, 0) and (6, 0) refer to the dimension-4 and dimension-6 contributions respectively. In order to express results in terms of our choice of input parameters (2.7), it is These hatted quantities are defined in terms of masses and couplings as in the SM. After rotation to the mass basis following the steps in appendix A one finds (3.5) Our notation is such that C f H is the coefficient which contributes to the hf f coupling after rotating to the mass basis. Its precise definition in terms of the coefficients multiplying the weak-basis operators in as has already been done in (3.5). In contrast, in the G F -scheme, which was used when calculating the partial NLO results of [55], one employs We have found the choice (3.6) to be particularly convenient for the full NLO calculation, since it involves only parameters which appear in the Lagrangian, and no tree-level dependence on four-fermion operators contributing to muon decay is introduced. Of course, it is a simple matter to convert the results obtained here to other renormalization schemes, provided all finite shifts between input parameters are known completely to NLO in SMEFT -for instance, calculations needed to trade v T for G F as in (3.7) have been obtained in [49]. The NLO counterterm is obtained by interpreting the external fields and parameters in (3.4) and (3.5) as bare ones, which are then replaced by renormalized ones before expanding the resulting expression to NLO in the couplings. The bare and renormalized fields are related through wavefunction renormalization factors according to where the second equality on each line is valid to NLO. For the masses, electric charge, and Wilson coefficients we write where M is a generic mass. The bare quantities in (3.8) and (3.9) are labeled with a superscript (0) while the renormalized ones are not, and the counterterm for an arbitrary quantity X is denoted by δX. These NLO counterterms are calculated in perturbation theory and receive both dimension-4 and dimension-6 contributions, which we denote by δX (4) and δX (6) , respectively. Inserting these expressions into (3.4) and (3.5) and keeping only the linear terms in δX gives an expression for the NLO counterterm for the decay amplitude. Writing this as the dimension-4 counterterm is while the dimension-6 counterterm is From the definitions ofĉ w andŝ w in (3.3) one finds that (3.14) The NLO counterterms are computed by specifying a renormalization scheme and evaluating one-loop Feynman diagrams as appropriate in that scheme. For the Wilson coefficients, we use the MS scheme, where the counterterms involve only UV poles in the dimensional regulator = (4 − d)/2. 1 In that case, we can read off the NLO counterterms from the anomalous dimension calculation performed in [24][25][26]. The counterterms take the form where we have introducedĊ with γ ij the anomalous dimension matrix. In general γ ij is not diagonal, so any Wilson coefficient counterterm is a linear combination of many other Wilson coefficients in the chosen basis. The wavefunction, mass, and electric charge counterterms are determined by calculating a set of one-loop integrals in the mass basis. The construction of these counterterms in SMEFT closely follows the procedure used in the Standard Model, as outlined, for instance, in [63]. Most of the details needed for h → bb decay in the on-shell scheme were given in [55]. However, while wavefunction renormalization factors are always evaluated on-shell, in the present work we aim to be flexible in the treatment of mass and electric charge renormalization, allowing for hybrid schemes which define some of these parameters in the on-shell scheme, and some in the MS scheme. In that case we must pay careful attention to tadpole contributions, as explained in section 3.2. There are also some subtleties in electric charge renormalization and Higgs-Z mixing once dimension-6 effects are included, which we cover in sections 3.1.1 and 3.1.2 below.
When necessary, we distinguish parameters in the on-shell scheme from those in the MS scheme through the notation where O.S. indicates the on-shell scheme and we have made the µ dependence in the MS parameter X(µ) explicit. The counterterms in the two schemes have the same UV divergences, but differ in the finite parts: the UV-finite part is set to zero in the MS scheme and determined through on-shell renormalization conditions in the on-shell scheme. We can therefore facilitate conversion between the MS and on-shell schemes by writing where the notation splits the counterterm into UV-divergent (X div. ) and UV-finite (δX fin. ) pieces. Results in the on-shell scheme are picked out by setting c X = 1, while c X = 0 picks out the MS scheme. This notation allows us to suppress the extra labels in (3.17) and refer instead to a generic quantity X, with the understanding that the renormalization scheme can be specified by adjusting the value of c X and the numerical value of X appropriately. We use this notation in section 4 and appendix B.

Electric charge renormalization
The one-loop counterterm (3.12) involves both SM and dimension-6 contributions from electric charge renormalization. The SM calculation simplifies due to electroweak Ward identities, which relate the f f γ vertex function to two-point functions through gauge invariance. Adapting the notation of [63] to our conventions, these allow one to write where as usual the superscript (4) refers to dimension-4 contributions. The object Σ AA T (Σ AZ T ) is the transverse component of the γγ (γZ) two-point function. The γZ two-point function is needed for charge renormalization in the SM because the photon can mix into a Z-boson through loop corrections before coupling to the fermion, and it is for the same reason that the axial-vector (a f ) and vector (v f ) couplings of the Z-boson to fermions enter the expression.
To renormalize the h → bb decay amplitude we also need the dimension-6 counterterm δe (6) . We have determined this expression by renormalizing the f f γ vertices directly, without using the SM Ward identities. We find by explicit calculation that Although the counterterm can be obtained from two-point functions alone, one can verify that the term multiplying Σ AZ T differs from the form (v f − a f )/Q f through terms involving the class-7 operators Q Hf . In fact, since v T /4ĉ wŝw for these operators, a naive generalization of the SM result (3.19) would lead to the contradictory result that electric charge renormalization depends on the fermion charge Q f . An important check on this expression is that the UV poles in the NLO decay amplitude cancel once it is used.

Higgs-Z mixing
In general, the SMEFT Wilson coefficients contain imaginary parts even after writing the Lagrangian in the mass basis. While these drop out of the NLO the decay rate, they appear in the NLO decay amplitude and introduce complications into the renormalization procedure which are irrelevant in the SM. One of these is mixing of the SM Higgs field h with the longitudinal component of the Z-boson and the neutral Goldstone boson φ 0 (in R ξ gauge) at the one-loop level. Since h and φ 0 are the real and imaginary parts of the neutral component of the Higgs doublet H after electroweak symmetry breaking respectively (see e.g. (A.3)), this mixing must involve a complex coupling. However, in the SM neutral-current couplings are real after transformation to the mass basis, so there is no such mixing at NLO. In SMEFT, however, diagrams of the type shown in Figure 1 contribute to the h → bb decay amplitude, where f is any massive fermion. The sum of diagrams yields a gauge-invariant result proportional to where the . . . refer to contributions from second-and third-generation fermions, which take on the same structure. The loop integrals multiplying η 5 contain UV divergences which are exactly canceled by the piece in the Wilson coefficient counterterm (3.15) involvingĊ bH , which was calculated with the SMEFT Langrangian in the unbroken phase of the theory (i.e. when the vacuum expectation value of the Higgs field vanishes) in [25]. While in the unbroken phase it is unambiguous that the η 5 term arises from mixing of real and imaginary parts of the complex Higgs doublet, in the broken phase the exact origin (but not the result itself) depends on the gauge: in unitary gauge it is due entirely to Higgs mixing with the longitudinal component of the Z-boson, while in R ξ gauge it is due to the sum of graphs containing Z and neutral Goldstone bosons

Tadpoles
In the on-shell renormalization scheme tadpole contributions cancel between different terms in the renormalized amplitude. For this reason, no tadpoles were included in the partial NLO calculation in the on-shell scheme in [55]. However, if some parameters are renormalized in the on-shell scheme and some in the MS scheme, then tadpole cancellations only happen at the level of UV-divergent parts of the amplitudes. Tadpoles remain in the finite parts, and must be taken into account to arrive at a gauge-invariant result. In fact, only upon the inclusion of tadpoles are the one-loop matrix elements (including wavefunction renormalization factors) and also mass and parameter counterterms individually gauge invariant [64].
There are various schemes for the treatment of tadpoles available in the literature. We have chosen to perform our calculations using the so-called "FJ tadpole scheme" [57], an excellent discussion of which is given in [65]. 3 As explained in that paper, a property of the FJ tadpole scheme is that it is equivalent to a scheme where tadpoles are not renormalized. In other words, tadpole renormalization can be taken into account simply by including tadpole topologies into any n-point amplitude entering a given calculation. This scheme applies not only to the Standard Model, but rather to generic theories, therefore it extends to SMEFT with no essential complications. We find this scheme to be particularly convenient, since it means that instead of adding explicit tadpole counterterms to the already lengthy expression (3.12), we need only include tadpole topologies into our diagrammatic calculations, which in any case have been automated.
In h → bb decay within the SM, tadpole contributions appear in the two-point functions used for mass and parameter renormalization through the diagrams shown in Figure 2(a)-(c). In h → bb decay within SMEFT, tadpoles appear not only in the two-point functions, but also in the bare decay amplitude through the diagram shown in Figure 2(d). We can write any of these diagrams as the product of the one-point tadpole function T shown in Figure 2(e) with a tree level graph, provided we include the appropriate Higgs coupling and propagator. We write the result for the tadpole function where (4) and (6) represent the SM and dimension-6 contributions respectively. In unitary gauge one has where f refers to quarks (q) or charged leptons (l) with N q c = 3, N l c = 1, and For the dimension-6 contribution in unitary gauge we find 26) and in Feynman gauge An interesting feature of SMEFT is that, in contrast to the SM, tadpole diagrams contribute to electric charge renormalization through the γγ two-point function. These contributions are proportional to the hγγ coupling in SMEFT, which is induced by class-4 operators and involves the combination of Wilson coefficients (3.28) Direct calculation in unitary gauge of the piece of the electric charge counterterm as described in section 3.1.1 yields the result where the extra superscript "cl.4" indicates restriction to class-4 operators in table 3. The term proportional to the SM tadpole function T un. arises through diagrams of the type shown in Figure 2(b) with IJ = γγ. In Feynman gauge the division into tadpole and the remaining contributions reads instead Feyn. , (3.30) but the end result is the same due to (3.24). This example illustrates the general feature that parameter counterterms are gauge invariant only after including tadpoles. The same is true of the sum of bare matrix elements and wavefunction renormalization factors, which is also a gauge-invariant object. The mechanism through which tadpoles ensure this gauge invariance is rather non-trivial. For instance, in contrast to the SM, tadpoles contribute directly to bare matrix elements through diagrams of the type shown in Figure 2(d). They also contribute to wavefunction renormalization of the b-quark field. Evaluating the tadpole contribution to the b-quark self-energy shown in Figure 2(a) and using it to extract the wavefunction renormalization factor using the convention of [55], one finds where T is the tadpole function in the chosen gauge. While this purely imaginary contribution drops out of the NLO decay rate, it is needed to ensure gauge invariance of the sum of the NLO matrix element and the wavefunction renormalization factors, and also plays a role in the cancellation of tadpoles in the on-shell scheme. These examples illustrate that while the treatment of tadpoles in SMEFT is conceptually the same as in the SM, the exact structure of tadpoles in the diagrammatic calculations is more involved. We have calculated all tadpole contributions to the bare matrix elements and counterterms appearing in the h → bb decay amplitude at NLO in unitary gauge and in Feynman gauge, and confirmed that the gauge dependence in the tadpole functions cancels against that in other diagrams, such that the counterterms for mass and electric charge renormalization, as well as the sum of the bare matrix element and the wavefunction renormalization factors, are separately gauge invariant.
We have also confirmed that tadpoles completely cancel when all parameters are renormalized in the on-shell scheme. However, QCD corrections to the b-quark mass and electric charge are sensitive to energy scales much smaller than the Higgs mass if the on-shell scheme is used, so one would prefer to renormalize such parameters in the MS scheme. In that case tadpole cancellation can no longer occur, and tadpoles enter the finite parts of the renormalized decay rate, carrying along with them corrections scaling as m 4 t /(v 2 T m 2 H ), which can lead to sizeable weak corrections. It is thus a non-trivial problem to find a renormalization scheme which is well suited for combining electroweak and QCD corrections in SMEFT. We deal with this issue in the next section.

Enhanced NLO corrections and decoupling relations
The size of perturbative corrections to the decay rate depends on the renormalization scheme, and it is an important question whether it is possible to find a scheme which reduces the size of higher-order corrections. In section 4.1 we identify sources of enhanced NLO corrections to the decay rate, and in section 4.2 we emphasise the importance of decoupling particles with masses at the electroweak scale from the MS definitions of the b-quark mass and electric charge when combining QCD and electroweak corrections in SMEFT.

Structure of the NLO decay rate
The full NLO result for the decay rate, including mass dependence of third generation fermions, is quite lengthy. However, it is possible to identify two sources of parametrically-enhanced corrections and their dependence on the renormalization scheme. The first is logarithms of the small ratio m b /m H , which appear in the QCD-QED type corrections contained in the piece Γ g,γ defined in (2.4). The result for these corrections in the m b → 0 limit is given in appendix B.2, using the notation in (3.18) in order to keep the dependence on the renormalization scheme for the b-quark mass explicit. Setting µ = m H and keeping only the logarithmic corrections in the result, one has where c m b = 1 (c m b = 0) yields the result in the on-shell scheme (MS scheme) for m b . It is simple to show that the decay rate in SMEFT depends only on the real parts of the Wilson coefficients, to the order which we are working. We have therefore used the notation that Re(C i ) ≡ C i in writing (4.1), and do this whenever we write an expression for the decay rate in what follows. Evaluating (4.1) numerically using the inputs in table 1 below yields We see that the QCD corrections are dominated by the double logarithmic term on the first line of (4.1). This term is of IR origin and cannot be removed through a choice of renormalization scheme. 4 It would need to be treated with QCD resummation techniques which we do not explore here. The single logarithmic term in the second and third line of (4.1) arises from the finite part of the counterterm for b-quark mass renormalization in the on-shell scheme.
Although not as large as the double logarithmic term, it is still a −50% correction to the LO result, which can be removed from the explicit NLO correction and resummed by using the MS scheme for the b-quark mass. We conclude that the QCD-QED corrections to the decay rate are best behaved in the MS scheme for the b-quark mass, which is indeed standard in SM computations.
The second source of potentially large corrections to the decay rate are weak corrections enhanced by powers of m 2 t /v 2 T , which appear in the object Γ t defined in (2.4). We give explicit results for the SM and dimension-6 corrections to Γ t in appendix B.3, as above using the notation in (3.18) in order to study the dependence on the renormalization scheme. The results show that in the MS scheme for the b-quark mass and electric charge the dominant contributions are due to tadpoles and scale as m 4 t /(v 2 T m 2 H ). The appearance of such corrections in NLO SMEFT calculations which make use of the MS scheme has been emphasized in h → γγ decay in [52], and in the partial NLO calculation of Z → bb in [30]. In the on-shell scheme tadpoles are absent and the leading corrections scale as m 2 t /v 2 T . We translate this into numerical results using the SM as an example. Keeping only the leading terms in the large-m t limit, one finds in the MS scheme while in the on-shell scheme where we have set µ = m t as appropriate in the large-m t limit and again used the inputs in table 1. The correction in the MS scheme for the b-quark mass is a −15% correction to the LO result and thus anomalously large for a weak correction, while that in the on-shell scheme takes on a much smaller value, in line with naive expectations. The numerical results for the dimension-6 contributions differ from operator to operator, but it is still the case that the corrections tend to be larger in the MS scheme than the on-shell one due to tadpole corrections scaling as m 4 t /(v 2 T m 2 H ). The upshot of this discussion is that while the QED and QCD corrections are best behaved in the MS scheme for m b , the electroweak corrections are better behaved in the on-shell scheme for m b and e, where tadpole contributions from heavy particles such as the top quark cancel. At least in the SM, an apparent compromise would be to use the MS scheme for all parameters appearing in the tree-level result, be it quark masses, the electric charge, M W or M Z . This is however an imperfect solution, for although in that case no explicit tadpoles appear in the NLO corrections, they reappear in the RG equations. Moreover, in SMEFT it is not possible to remove all explicit tadpole contributions in this manner, since in contrast to the Standard Model they can also appear in the matrix elements for h →bb, through contributions such as that shown in Figure 2(d).
The resolution to this dilemma is to renormalize the b-quark mass and electric charge such that the QCD-QED corrections are treated in the MS scheme, while weak corrections involving the top quark and heavy electroweak bosons are treated in the on-shell scheme. In that way contributions from potentially large tadpole corrections cancel, but logarithms of m b /m H can still be resummed in the MS scheme. At the technical level, the simplest way to implement such a scheme is to make use of so-called "decoupling relations".

Decoupling relations
Decoupling relations connect MS-renormalized parameters in SMEFT with those defined in a low-energy theory where the top quark and electroweak bosons are integrated out. A detailed discussion of this in the SM for the b-quark mass defined in the MS scheme can be found in [66]. We shall consider only the dimension-4 piece of this low-energy theory, which we refer to hereafter simply as QED×QCD. This amounts to neglecting terms which scale as e.g. m 2 b /M 2 W , which are numerically negligible compared to the dimension-4 terms. We can then write the decoupling relations as where the parameters on the left-hand side are defined in SMEFT, and those on right-hand side, with the superscript , are defined in QED×QCD. These parameters obey the RG equations In what follows we will make use of the LO anomalous dimensions γ i , which read where N g = 3 is the number of fermion generations, Q u = 2/3 for up-type quarks, and α ( ) (µ) ≡ [e ( ) (µ)] 2 /(4π). The parameter m where α(M Z ) ≈ 1/129 compared the on-shell value α ≈ 1/137 (see e.g. [67]). The ζ i in eq. (4.5) are decoupling constants. They are determined by using the relation between the MS and on-shell parameters in the two theories. These take the form where we have used that the on-shell parameters e and m b are defined through non-perturbative renormalization conditions and do not depend on the Lagrangian. The z i factors are finite and determine the perturbative shifts between the on-shell and MS parameters. They fix the decoupling constants through the relations where i = e, b.
We write the perturbative expansion of the decoupling constants in SMEFT as where the superscripts (4, 1) and (6, 1) follow the notation of (2.3). At NLO the decoupling constants are proportional to the finite parts of heavy-particle contributions to the NLO renormalization constants. The expression for ζ e is compact. The SM expression is 12) and the SMEFT result reads where Q t = 2/3 is the charge of the top quark and the term on the second line of (4.13) is the UV-finite part of the class-4 electric charge counterterm (3.30) with m b → 0. The full results for ζ b are somewhat lengthy, and are given in computer files in the arXiv version of this paper. They simplify considerably in the large-m t limit, where they read where δb t is the UV-finite part of δm b in this limit and is given in (B.17).
The h → bb decay rate written in terms of the QCD×QED parameters m ( ) b and e ( ) , which we denote by Γ , is simple to obtain from the decay rate in terms of the parameters m b and e in the full SMEFT, which we denote by Γ. The LO results are the same up to a renaming of the parameters, and the NLO results are given by Eq. (4.15) is obtained by inserting (4.5) into Γ (0) and expanding to NLO. The same result can be obtained by replacing δm b /m b → δm b /m b + ζ b and similarly for δe/e in the NLO counterterms (3.11) and (3.12), and for this reason evaluating the decay rate using (4.15) is equivalent to using a new renormalization scheme. After splitting up the decay rate in this scheme as ,g,γ + Γ ,t + Γ ,rem , (4.17) it is possible to list a simple and illustrative result for the QCD×QED and large-m t limit of the weak corrections. In terms of the quantities defined in appendix B, we have The interpretation is that the QCD×QED corrections are calculated in the MS scheme, while contributions from top-quark loops are calculated in the on-shell scheme, where tadpoles cancel. This pattern holds for heavy gauge-boson contributions to the decay rate. In fact, after decoupling, heavy-particle contributions are effectively calculated in the on-shell scheme, so that the only non-vanishing tadpole contributions are suppressed by powers of light fermion masses and are negligible numerically.

Numerical results
In this section we present results for the h → bb decay rate at NLO in SMEFT. We first give numerical results with the default choice µ = m H in section 5.1, and then perform a study of perturbative uncertainties due to scale variations in section 5.2. Throughout the analysis we use the renormalization scheme defined in (4.15). Since the decoupling relations used in that scheme are valid in the limit where all fermion masses except the top-quark mass m t vanish, we shall use this approximation in presenting the numerical results. The dominant corrections to this limit scale as m 2 b /M 2 W and typically change the NLO corrections at the 1% level and are thus irrelevant for our discussion. The input parameters needed in the analysis are listed in table 1.

Results at µ = m H
To quote results for the dimension-6 contributions, we make the dependence on Λ NP explicit by defining dimensionless Wilson coefficients according tõ Contributions to the decay rate from dimension-6 operators are then suppressed by an explicit power ofv ( ) (µ) 2 /Λ 2 NP , which for the input parameters in table 1 leads to a roughly 5% suppression factor for Λ NP = 1 TeV andC i ∼ 1.
We shall present numerical results normalized to the LO SM decay rate. We thus define In quoting this result, we have kept a factor ofv/m b ∼ 80 multiplying theC bH contribution symbolic. We do this to highlight the fact that theC bH contribution to the decay rate scales as m b rather than m 2 b as in the SM, which can be seen explicitly in (B.2). The same is true of six additional coefficients which enter the decay rate at NLO:C bG ,C bW ,C bB ,C Htb ,C (1) qtqb andC (8) qtqb . It is worth mentioning that if MFV is imposed then all of these coefficients scale as y b ∼ m b /v, so that their contributions to the decay rate scale as m 2 b . However, our results are not limited to MFV, so keeping factors ofv/m b symbolic when multiplying the coefficients mentioned above is simply a matter of convenience. For the same reason, when quoting results from operators such as Q bB or Q bG where gauge bosons couple through field strengths rather than covariant derivatives, we keep enhancement factors of 1/e or 1/g s compared to the SM contributions symbolic. With these conventions, the NLO result can be written as Hq + − 7.9C Ht + 5.8C (1) Hq 22 Hl 22

+ C (3)
Hl 11 By far the largest NLO correction is fromC HG , which is a QCD effect enhanced by a double logarithm in m b /m H as described in section 4.1. Order 10% corrections (in units ofv 2 /Λ 2 NP ) arise fromC (1) Hq ,C Hq andC Ht . In total there are 16 operators which contribute at greater than a percent level to the decay rate, 12 of which first appear at NLO.
Generally speaking, an operator gives a significant contribution only if it involves QCD or large-m t corrections. To illustrate the relative importance of these two effects, we show in table 2 the division of the NLO corrections to operators appearing at tree level into QCD-QED corrections, large-m t corrections, and remaining corrections (denoted by Γ ,g,γ , Γ ,t , and Γ ,rem ). For the dimension-6 operators, the numbers are defined as the contribution of the Wilson coefficientC i to Γ (1) divided by its contribution to Γ (0) . The results show that while the QCD corrections are dominant, the electroweak corrections are non-negligible and depend strongly on the Wilson coefficient. For instance, the electroweak corrections fromC HD are −11%, while those fromC bH are +3%. Therefore, approximating the NLO corrections in SMEFT by multiplying the tree level result with a universal K-factor derived from the SM QCD corrections would be a poor estimate to the full calculation performed here. We also note that the large-m t corrections indeed make up the bulk of the electroweak corrections, although deviations from that approximation are between 10 − 40%. We have observed that this pattern holds for the other coefficients appearing in the NLO result.

Scale uncertainties
So far we have given results only at µ = m H . In this section we address two obvious questions concerning scale uncertainties: first, can the size of NLO corrections be reliably estimated through scale variations of the LO result, and second, what is the residual uncertainty beyond NLO?
We shall study these questions as typical in a perturbative analysis, namely by varying unphysical renormalization scales up and down by factors of two and taking the change in the decay rate as a measure of the uncertainty due to uncalculated, higher-order corrections. A difference in SMEFT compared to the SM is that while all parameters in the SM Lagrangian have been determined to good accuracy numerically, the exact values of the Wilson coefficients in SMEFT are largely unknown. Therefore, when performing scale variations, we give results symbolically in terms of the Wilson coefficients at a fixed reference scale. In our case, the natural choice of this reference scale is µ = m H , therefore our task is to express the Wilson coefficients C i (µ) in terms of the C i (m H ). This is achieved by solving the RG equations for the Wilson coefficients.
For variations of µ by factors of two, µ ∼ m H parametrically, so we can use the fixedorder expansion of the RG equations rather than the exact, exponentiated solution. In fact, the same holds for the SM masses and couplings renormalized in the MS scheme. Given that the anomalous dimensions of the Wilson coefficients are known only to one-loop, we use this same level of accuracy for the SM parameters throughout this section. The solutions of the RG equations to NLO in fixed order read The number of light quarks is n l = 5 and C A = 3. Results for γ e and γ b were given in (4.7), andĊ i was defined in (3.16). We have written (5.5) in a fashion which emphasizes that it is possible to use different renormalization scales µ C and µ R for the Wilson coefficients and the SM parameters, respectively. Until this point we have set µ C = µ R = µ, but in our scale uncertainty analysis it will be useful to consider independent variations of these scales. These scales appear not only implicitly in the Wilson coefficients, b-quark mass, and the strong and electromagnetic coupling constants, but also in explicit logarithms in the NLO decay rate. The explicit logarithmic dependence on the two scales in the NLO dimension-6 results can be reconstructed from the result at µ R = µ C = µ by using the RG equations along with the requirement that the decay rate is independent of the renormalization scales up to terms of order NNLO and higher. The results can be written as , α s (µ)} are the MS-renormalized parameters appearing in the calculation. By definition Γ (6,i) (µ, µ) = Γ (6,i) (µ).
With these pieces at hand, we obtain scale uncertainties using the following procedure. For the SM results, we vary the scale µ R up and down around its default value m H . For the dimension-6 results, we can vary both µ R and µ C using (5.7). The default setting is µ R = µ C = m H . We then assign an uncertainty to each scale individually by varying it up and down by a factor of two while leaving the other scale fixed, and add the resulting uncertainties from the independent µ R and µ C variations in quadrature to obtain a total uncertainty. The numerical values of the scale-dependent parameters at the different scales are determined in terms of their values at m H using (5.5). This results in (1) qtqb ± 0.03C tW e ( ) ± 0.03(C HW +C tH ) + . . . , (5.8) where the ellipses indicate dimension-6 terms which contribute less than 3% in units ofv 2 /Λ 2 NP . At NLO, we find where now the ellipses indicate terms with uncertainties smaller than 3%, other than those which appear already in eq. (5.8).
We see that the NLO calculation generally leads to a considerable reduction in the scale uncertainties compared to LO. For the operators already appearing at tree level, the NLO corrections are on the upper limits of what one would estimate through scale variations of the LO result. For operators which first appear at NLO, varying the scale in the LO results generally estimates the size of the NLO contribution quite well. A major exception is theC HG coefficient. In that case the size of the NLO correction is dramatically underestimated by scale variations in the LO result, and in fact the NLO result has a larger perturbative uncertainty associated with it than the leading one. This is not surprising, given that the large correction fromC HG is completely unrelated to RG running, as explained in section 4.1. A consequence of this is that a new coefficientC tG , which arises predominantly through the running ofC HG , is a significant source of uncertainty in the NLO calculation.
Needless to say, the uncertainties assigned to the decay rate through the above procedure are just estimates, and other methods for varying the scales are possible. The simplest one is to set µ R = µ C = µ and obtain uncertainties by varying the single scale µ up and down by a factor of two. Analytic results for the uncertainties in the LO result, which we denote by δΓ (i,0) , obtained in this way are quite simple: dropping terms of order NNLO and higher, one has δΓ (4,0) = ±2 ln(2)Γ (4,0) (γ b + γ e ) , where all scale-dependent quantities are to be evaluated at µ = m H . Compared to the results (5.8) using the quadrature method, only contributions from the dimension-6 coefficients appearing in the LO matrix elements are changed. Numerically evaluating (5.10) leads to the following result for those coefficients in units ofv 2 /Λ 2 NP : The result forC bH is almost identical to that obtained with the quadrature method, but the uncertainties assigned to the other coefficients are significantly smaller. Especially those for C H2 (3%) andC HD (2%) are artificially small uncertainties to assign to an LO calculation, and for this reason we have chosen the quadrature method by default. Even more conservative methods could be used, for instance a scan over µ C and µ R which takes into account simultaneous but uncorrelated variations to include choices such as µ R = m H /2, µ C = 2m H where neither scale is at its default value, but we do not explore such options here. Our main message is that it is important to assign uncertainties to the LO result, and these uncertainties are significantly reduced through the NLO calculation.

Conclusions
We have calculated the full set of NLO corrections to h → bb decay in SMEFT, obtaining contributions from the 45 dimension-6 Wilson coefficients which enter the decay rate at this order. These results form the basis for any future precision analysis of this decay in effective field theory. While the renormalization of the electroweak sector of SMEFT is conceptually similar to the SM, in section 3 we highlighted some technical differences regarding charge renormalization and also Higgs mixing with the Z and neutral Goldstone bosons. Moreover, the structure of tadpole cancellation in the h → bb decay amplitude in the on-shell renormalization scheme is rather intricate in SMEFT, since contrary to the SM, tadpole contributions to the matrix elements, b-quark wavefunction renormalization and electric charge renormalization must be taken into account.
Our calculation includes both electroweak and QCD corrections, which has led us to explore hybrid renormalization schemes where heavy particle masses are renormalized onshell while the b-quark mass and electric charge are renormalized in the MS scheme. In such schemes tadpoles do not cancel from the decay amplitude, need to be included in order to obtain gauge invariant decay results, and can lead to enhanced electroweak corrections. In section 4 we showed how these enhanced electroweak corrections can be removed from the decay rate by decoupling contributions from electroweak-scale masses from the running of MS renormalized parameters, which are then defined in a low-energy version of QED×QCD. We obtained the decoupling constants for the electric charge and b-quark mass to NLO in SMEFT, and used them to calculate the decay rates in a hybrid renormalization scheme which simultaneously avoids enhanced tadpoles corrections from the electroweak sector and resums UV logarithms in m b /m H in the QCD one.
In section 5 we gave numerical results in the aforementioned renormalization scheme with the scale choice µ = m H for all MS-renormalized parameters, namely the Wilson coefficients as well as the b-quark mass and electric charge. We also studied the perturbative uncertainties in the LO and NLO results as estimated through scale variations. We found that while in general the NLO corrections stabilize the scale dependence of the decay rate, genuine NLO effects inaccessible to an RG analysis based on scale variations can be significant. That said, we advocated introducing two renormalization scales, one for the Wilson coefficients and one for the MS renormalized b-quark mass and electric charge, and varying them independently in order to generate more reliable uncertainty estimates than those obtained from varying a common scale µ .
The analytic results for the NLO decay rate in SMEFT are rather lengthy and included in computer files with the arXiv submission of this article, both with the full m b dependence, which will be useful for future validations of our results, and in the m b → 0 limit, which is sufficient for phenomenology. We believe that the renormalization procedure and uncertainty analysis performed here can serve as a template for future NLO SMEFT calculations which aim to include electroweak and QCD corrections in a single framework. quantity λ in the Higgs potential can be eliminated in terms of the input parameters (2.7) according to

A.2 Gauge fields
In the following section we review the rotation to the mass basis of the gauge fields in SMEFT, closely following the procedure in [26]. We denote the covariant derivative in the electroweak sector of the SM by , and the generators are denoted (gτ ) a = (g 2 τ 1 , g 2 τ 2 , g 2 τ 3 , g 1 Y ), where τ I = σ I /2 with σ I the Pauli matrices and Y the hypercharge.
When including dimension-6 operators we must first redefine the gauge fields as to ensure correct gauge field normalization. Additionally, we modify the couplings as such that the combinations g 1 B µ =ḡ 1 B µ and g 2 W I µ =ḡ 2 W I µ remain unchanged. It can be shown thatḡ 1 andḡ 2 can be written in terms of the physical input parameters listed in (2.7) asḡ The class-4 operator Q HW B introduces a kinetic mixing term between the W 3 µ and B µ gauge fields not seen in the SM, which is of the form ∼ − 1 2 v 2 T W 3 µ B µ . This term can be removed by a linear shift in these fields, which proceeds as such that the new 'primed' gauge fields have diagonal kinetic terms. These are rotated to the mass basis according to whereÃ a µ comprises the physical gauge fields asÃ µ = (W + µ , W − µ , Z µ , A µ ), and R is given by (A.14) With this notation, the relation between the weak-basis fields A a µ and the mass basis fields A a µ is In terms of the input parameters in (2.7), the explicit definitions of the photon and Z-boson fields in terms of the weak-basis fields is Furthermore, the dimension-6 SMEFT covariant derivative in the mass basis is given by where Q = τ 3 + Y and τ ± = (τ 1 ± iτ 2 )/ √ 2.

A.3 Gauge fixing in R ξ gauges
Gauge fixing in SMEFT has been discussed in [68][69][70]. In this section we explain our own implementation, which we have used when verifying the gauge independence of the decay rate and counterterms with explicit one-loop computations. Throughout this section we follow closely the notation used for gauge fixing in the SM as presented in [71]. We parametrise the Higgs doublet in terms of real scalar fields as and use the real representation of the generators, T a = −iτ a , where the τ a were defined below (A.6). We expand each φ i about its vacuum expectation value, denoted φ i = φ 0 i as where χ i =4 are the Goldstone bosons, χ 4 is related to the physical Higgs boson, h, and In R ξ gauges one aims to remove the Goldstone-gauge boson mixing terms, which in the SM take the form where (gT ) a = (g 2 T 1 , g 2 T 2 , g 2 T 3 , g 1 T 4 ). The i = 4 component in (A.20) gives no contribution to the Lagrangian. We now include dimension-6 effects in SMEFT. We begin by defining the canonicallynormalized fields of the Higgs doublet in (A.3) in terms of those in (A.19) via the transformation such that the χ i are related to the fields in (A.3) by Moreover, we replace the gauge fields and couplings as in (A.7), (A.8) and (A.10) such that all the Goldstone-gauge mixing terms of the SMEFT Lagrangian may be written where the second term on the first line of (A.23) is the contribution arising from the explicit presence of the C HD Q HD term in the dimension-6 SMEFT Lagrangian. Here we have introduced the object (gT ) a , which is defined similarly to (gT ) a in (A.20), but with all instances of the gauge couplings replaced as g i → g i , and further defined 'primed' generators where M ab is given in (A.11), and also the object 25) where in the final line we have used that X has only diagonal elements, X 11 = X 22 = (X −1 ) 11 = (X −1 ) 22 = 1, (1 +v 2 T 2 C HD )X 33 = (X −1 ) 33 and that the X 44 component gives no contribution. In order to calculate the matrix (gF) a i we use, for example, that (gT ) 1 φ 0 equals g 2 v T /2 times a unit vector in the φ 1 direction. One finds (A.26) We follow the Faddeev-Popov gauge-fixing procedure such that the SMEFT gauge-fixed generating functional Z takes the form where G a is the gauge-fixing function and the object (α /g) b is defined below. We choose the gauge-fixing function in (A.27) as which defines the R ξ gauges in SMEFT. 5 We see that the form of the gauge-fixing function in (A.28) resembles that of the R ξ gauges in the SM with the gauge fields replaced by their primed counterparts and F replaced with F . The Goldstone-gauge boson mixing terms in (A.23) are then removed by the − 1 2 (G) 2 term in (A.27). Interactions of SM particles with ghost fields arise through the functional determinant in (A.27), for which we must determine the variation of G a under arbitrary gauge transformations. The gauge transformation of the scalar fields may be written where the second relation defines the object (α/g) a and the third relation defines the object We may use (A.21) and (A.29) to find the gauge transformation of χ i : where we have defined the object (gT ) a ij ≡ (X −1 ) ik (gT ) a kl X lj . Explicitly (gT ) a ij acts on χ i as (for brevity and as no other terms enter our calculation, we give only the Higgs contributions to this term) We may similarly write the transformation of the unprimed gauge fields as The object f abc = abc if a, b, c ∈ 1, 2, 3 and vanishes otherwise, which we have used to replace α b → g 2 (α/g) b in the above equation. The form of δA a µ in terms of the object (α /g) a is then found using (A.10), (A.30) and (A.33) We can now calculate the functional derivatives needed to evaluate (A.27) using the results in (A.31) and (A.34). First, one has (note that the gauge fields here are the unprimed gauge fields), where the explicit result is From (A.35) and (A.31), the variation of the gauge-fixing function, G a in (A.28) is Following the usual procedure the ghost Lagrangian is The ghost fields in (A.38) are given by c a = (c W 1 , c W 2 , c W 3 , c B ), and similarly for the fields in c a . The form of the ghost mass matrix in (A.38) is which is diagonalized by the matrix R in (A.13) such that The ghosts in the mass basis, denoted u a and u a , are thus related to those in the weak basis by where u a = (u W + , u W − , u Z , u A ), and similarly for u a . With the gauge fields A µ written in terms of the mass basis as described in (A.15), the ghost Lagrangian in the mass basis is therefore Although our derivation is rather different, we find that the Feynman rules produced by the Lagrangian in (A.42) exactly match those found in [68].

A.4 Yukawa sector
The fermion masses in SMEFT involve the Wilson coefficients of class-5 operators as well as the SM Yukawa matrices. The relevant part of the Lagrangian (following the convention used in [24][25][26]) is given by where the subscripts j and r i are SU (2) and generation indices respectively. In what follows we perform rotation to the mass basis using the down-type quarks as an example and suppress the explicit addition of the hermitian conjugate (+h.c.). After spontaneous symmetry breaking in unitary gauge and keeping only dimension-6 terms one finds where L yuk is defined as the term proportional to the hψψ operator. Additionally we have included the subscripts L and R on the quark fields to denote their handedness. As usual, we perform rotations on the quark fields to go to the mass basis where [m d ] = diag(m d , m s , m b ). After the field rotation, the hψψ term becomes Thus, in the mass basis the Wilson coefficients contributing to hψψ couplings are a linear combination of those in the weak eigenstate basis. Similar results can be derived for any Wilson coefficient C m i multiplying a mass-basis operator containing fermions. Note that in contrast to the SM, SMEFT contains flavour-violating Higgs couplings even in the mass basis. However, in our calculation we approximate the CKM matrix by the unit matrix, in which case these flavour-violating couplings do not contribute to the NLO h → bb decay rate at dimension-6. In fact, within this approximation there are no transitions between fermion generations, which allows us to introduce a compact notation for Wilson coefficients such as (A.49) which multiply mass-basis operators. First, for operators involving right-handed fields we can always indicate the generation by the explicit flavour. Examples of this are C bH ≡ C m and similarly for any fermion f . Some Wilson coefficients for operators containing left-handed fields use the subscripts q r and r , so it is not possible to indicate the doublet generation r through the flavours it contains. However, the third generation plays a prominent role in our calculation, so our convention is to suppress any dependence on r = 3 but display explicitly the flavour indices only on operators involving first-and second-generation fermions, which appear through electroweak boson self-energies and tadpoles. Examples of operators in this notation are where the first coefficient multiplies a mass-basis operator with field contentttbb and the second coefficient multiplies a mass-basis operators with fermion content cc and ss.
An important feature of SMEFT in the mass basis is that couplings between left and right-handed fields are not always associated with powers of the fermion mass, as in the SM.
For instance, the C bH operator contains a hbb coupling which is not proportional to the bquark Yukawa, which is y b ≈ √ 2m b /v T in the mass basis. For this reason h → bb offers an important probe on the flavour structure of SMEFT. However, in this work we are interested in the structure of NLO contributions in SMEFT rather than questions of flavour, so in our numerical analysis it is convenient to display results in such a way that all contributions to the decay rate multiply a symbolic factor of m 2 b /v 2 T as in the SM. We emphasize that this is not a restriction of our calculation but rather a matter of convenience. However, if the Wilson coefficients are generated by a new physics scenario which respects Minimal Flavour Violation (MFV) [72] it is something which occurs naturally. See refs [26,73] for further discussion on this in the context of SMEFT.

B Analytic results
In this section we give analytic results for the LO decay rate and the NLO QCD-QED corrections Γ g,γ in the small-m b limit used in our numerical analysis, as well as the large-m t corrections Γ t . We give results which can be easily converted between the on-shell and MS schemes for X ∈ {m b , e} using the notation in (3.18). We will also need to split the finite part of the counterterms in the on-shell scheme into QCD-QED, large-m t , and remaining pieces. To do so we define where the superscript i = 4, 6 labels the NLO contribution from dimension-i operators. We use this notation throughout the section.

B.1 LO decay rate
The LO contributions to the decay rate as defined in (2.2) are given by

B.2 QCD-QED corrections
The NLO result for the QCD-QED corrections in the SM can be written as Γ (4,1) g,γ = Γ (4,1) g,γ + 2c m b Γ (4,0) δb (4) g,γ , (B.4) while that in SMEFT takes the form Γ (6,1) g,γ = Γ (6,1) g,γ + 2c m b Γ (4,0) δb (4) g,γ  where δb (4) g,γ = − and we have used that δb (6) g,γ = 0 in the small-m b limit. The Γ g,γ are the QCD-QED corrections to the decay rates in the MS scheme for m b . The QCD corrections were obtained in [56]. Most of the QED corrections can be derived from those results by making appropriate replacements. The exception is the contribution proportional to the hγZ vertex in SMEFT, which arises from the real and virtual emission diagrams in Figure 3 and has no analogue in QCD. We have obtained the contributions from these diagrams to the decay rate by evaluating and adding together the virtual and real corrections as in (2.5). This new result together with the other QCD-QED corrections in the small-m b limit can be written as δê (6) t,tad = 8N c C HBĉ