Factorization at subleading power and endpoint divergences in h → γγ decay. Part II. Renormalization and scale evolution

Building on the recent derivation of a bare factorization theorem for the b-quark induced contribution to the h → γγ decay amplitude based on soft-collinear effective theory, we derive the first renormalized factorization theorem for a process described at subleading power in scale ratios, where λ = mb/Mh « 1 in our case. We prove two refactorization conditions for a matching coefficient and an operator matrix element in the endpoint region, where they exhibit singularities giving rise to divergent convolution integrals. The refactorization conditions ensure that the dependence of the decay amplitude on the rapidity regulator, which regularizes the endpoint singularities, cancels out to all orders of perturbation theory. We establish the renormalized form of the factorization formula, proving that extra contributions arising from the fact that “endpoint regularization” does not commute with renormalization can be absorbed, to all orders, by a redefinition of one of the matching coefficients. We derive the renormalization-group evolution equation satisfied by all quantities in the factorization formula and use them to predict the large logarithms of order ααs2Lk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha \alpha}_s^2{L}^k $$\end{document} in the three-loop decay amplitude, where L=ln−Mh2/mb2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ L=\ln \left(-{M}_h^2/{m}_b^2\right) $$\end{document} and k = 6, 5, 4, 3. We find perfect agreement with existing numerical results for the amplitude and analytical results for the three-loop contributions involving a massless quark loop. On the other hand, we disagree with the results of previous attempts to predict the series of subleading logarithms ∼ααsnL2n+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sim {\alpha \alpha}_s^n{L}^{2n+1} $$\end{document}.


Introduction
Soft-collinear effective theory (SCET) [1][2][3] provides a convenient framework for addressing the problems of scale separation and factorization in high-energy physics using the powerful tools of effective field theory. Much recent work has focused on exploring the structure of factorization at subleading order in power counting -a problem that turns out to be unexpectedly subtle and full of complexities. Specific applications discussed in the literature include the study of power corrections to event shapes [4] and transverse-momentum distributions [5,6], the threshold factorization for the Drell-Yan process [7,8], and the factorization of power-suppressed contributions to Higgs-boson decays [9,10]. One finds that such factorization theorems contain a sum over convolutions of Wilson coefficients with operator matrix elements, where the relevant SCET operators mix under renormalization. Several new complications arise, which do not occur at leading power. The most puzzling one is the appearance of endpoint-divergent convolution integrals over products of component functions each depending on a single scale [6,[8][9][10][11][12][13][14]. In some sense, such endpoint divergences indicate a failure of dimensional regularization and the MS subtraction scheme, because some of the 1/ n pole terms are not removed by renormalizing the individual component functions, and hence naive scale separation is violated. Standard tools of renormalization theory are then insufficient to obtain well-defined, renormalized factorization theorems involving convergent convolutions over renormalized functions. Indeed, for none of the abovementioned examples is it currently known how to formulate a theoretically consistent renormalized factorization theorem. This would, however, be needed in order to fully establish SCET as a versatile tool and apply it to several observables of phenomenological interest.
In a recent paper [9], two of us have started a detailed discussion of factorization at subleading power within the framework of SCET. As a concrete example, we have considered the decay amplitude for the radiative Higgs-boson decay h → γγ mediated by the Higgs coupling to bottom quarks. In the limit m b M h one finds for the corresponding contribution to the decay amplitude 3 ) denote the mass, Yukawa coupling and electromagnetic coupling of the b quark, and ε µ (k i ) are the polarization vectors of the two photons. At leading doublelogarithmic order, the large logarithms in this amplitude were resummed a long time ago using conventional tools of perturbative quantum field theory [15,16] (see [17,18] for more recent related work), and it was found that However, we will see later that the sum of this series provides a poor numerical approximation to the decay amplitude. In order to go beyond the leading double-logarithmic approximation it is necessary to start from a consistent all-order factorization theorem,

JHEP01(2021)077
which properly separates the relevant scales in this process. In [9] we have factorized the decay amplitude and expressed it in terms of convolutions of bare matching coefficients with unrenormalized operator matrix elements. The convolution integrals contain endpoint divergences that require both dimensional and rapidity regulators. We have shown that, to all orders of perturbation theory, the endpoint divergences cancel out in the sum of all contributions to the factorization theorem, and they can be removed by suitably rearranging the factorization formula.
Here we continue our study of the b-quark induced h → γγ decay amplitude and derive a renormalized factorization theorem for this process. The establishment of such a factorization formula -the first of its kind -is the main accomplishment of the present work. In section 2 we recall the factorization theorem for the b-quark induced h → γγ decay amplitude as derived in [9]. A crucial step in the derivation of this formula made use of two D-dimensional refactorization conditions, which were derived from the requirement that in the sum of all terms the dependence on the rapidity regulator must cancel to all orders of perturbation theory. In section 3 we prove these refactorization conditions using techniques from SCET. The renormalized form of the factorization formula is derived in section 4. We first discuss the renormalization of the various operators and matching coefficients and derive their renormalization-group (RG) evolution equations using standard tools of quantum field theory. There is one important subtlety related to the fact that the presence of hard cutoffs on the convolution integrals, which is a consequence of the regularization of endpoint divergences, does not commute with renormalization. In section 4.5 we show that, to all orders of perturbation theory, the additional terms encountered when the cutoffs are moved from the bare to the renormalized convolutions amount to an extra contribution to one of the matching coefficients. The evolution equations satisfied by the renormalized matching coefficients and matrix elements are derived in section 5. Using these equations, we predict in section 6 the large logarithms of order α b α 2 s L k with k = 6, 5, 4, 3 in the three-loop decay amplitude, finding complete agreement with existing multi-loop results in the literature [19,20]. However, we do not confirm previous predictions for the series of subleading logarithms of order α b α n s L 2n+1 [16,21], which were based on conventional resummation techniques. Section 7 contains our conclusions. Several technical details are relegated to a series of appendices. A short letter summarizing our main results has recently appeared in [22]. There we have discussed the resummation of large logarithms at next-toleading logarithmic order. A complete resummation of large logarithms in RG-improved perturbation theory is left for future work.

Factorization formula in terms of bare objects
As shown in (1.1), the b-quark induced contribution to the h → γγ decay amplitude receives large logarithms of the form α b α n−1 s L k with k ≤ 2n. In order to resum these logarithms in a systematic way it is necessary to factorize the amplitude into objects depending only on one of the three relevant scales, set by the Higgs-boson mass M h , the mass m b of the bottom quark, and the intermediate scale In [9], two of us have derived a "bare" factorization theorem, which accomplishes this. It contains three terms consisting JHEP01(2021)077 contains a Higgs field coupled to two collinear gauge fields describing photons moving along opposite light-like directions n 1 and n 2 ≡n 1 . Here and below, fields without an argument are located at the spacetime point x = 0. The canonical choice of the reference vectors is n µ 1 = (1, 0, 0, 1) and n µ 2 = (1, 0, 0, −1). This operator descents from full-theory graphs in which all internal momenta are hard, of order M h . Next, the operator contains a Higgs field, an n 2 -collinear photon field and two n 1 -collinear b-quark fields, which share the momentum k 1 of the other photon. This operator is generated by fulltheory graphs in which a loop momentum is collinear with the photon direction n 1 and JHEP01(2021)077 carries virtuality of order m b . The factor 2 in front of this contribution in (2.1) arises because there is an analogous contribution with n 1 and n 2 interchanged. The symbols A µ n 1 and X n 1 in the above definitions denote effective photon and b-quark fields defined in SCET (the so-called "gauge-invariant building blocks" [2,23]), which differ from the ordinary quantum fields A µ and ψ in that they contain collinear Wilson lines in their definition and that they obey the constraintsn 1 · A n 1 = 0 and / n 1 X n 1 = 0. Note that the Feynman rule for the vector field A µ contains a factor of e b , which is the reason why we have divided by e 2 b in the definition of O 1 . The symbol ⊥ on 4-vectors indicates the components orthogonal to the light-cone basis vectors n 1 and n 2 . Finally, the operator contains the time-ordered product of the scalar Higgs current with two subleading-power terms in the SCET Lagrangian [3], in which hard-collinear fields are coupled to a soft quark field. It arises from full-theory graphs containing a soft quark propagator between the two photons, with all momentum components of order m b . In terms of gauge-invariant building blocks, the relevant subleading-power terms in the SCET Lagrangian read where G µ n 1 is the building block for the hard-collinear gluon field. In interactions of hardcollinear fields with soft fields the soft field operators must be multipole expanded for consistency [3,24], and we denote x µ − = (n 1 · x) n µ 1 2 and y µ + = (n 2 · y) n µ 2 2 . The h → γγ matrix element of O 3 in (2.1) can be factorized further into a convolution of two radiative jet functions with a soft function [9], i.e. (for simplicity we use the default choices of the reference vectors n 1 and n 2 , such thatn 1 Contrary to the 4-vectors x − and y + introduced above, the integration variables + and − correspond to the light-cone components n 1 · and n 2 · of a soft 4-momentum µ . Here and below we omit the photon polarization vectors when presenting expressions for operator matrix elements. The radiative jet function J(p 2 ) has been studied first in the calculation of the decay amplitude for the rare exclusive decay B − → γ l −ν l in the context of QCD factorization [25]. The soft function S(w) is defined in terms of the discontinuity of a soft quark propagator dressed with soft Wilson lines oriented along the light-like directions n 1 and n 2 . For a more detailed discussion about the derivation of the factorization theorem and the precise definition of the various SCET fields, operators, jet and soft functions the reader is referred to [9], where we have shown that the three operators O i form a basis of O(λ 3 ) SCET operators contributing to the decay h → γγ, and that the sum of the JHEP01(2021)077 three terms in (2.1) correctly reproduces the decay amplitude at two-loop order. Here λ ∼ m b /M h is the expansion parameter of the effective theory.
Major complications arise from endpoint-divergent convolution integrals in the second and third term in (2.1), which need to be properly identified and regularized. The integral over z in the second term contains singularities at z = 0 and z = 1. Likewise, the integrals over + and − in (2.6) contain singularities at ± = ∞. Some of these endpoint divergences are regularized by the dimensional regulator D = 4 − 2 , but others require an additional rapidity regulator [26][27][28]. In [9] we have regularized the rapidity divergences by means of an analytic regulator imposed on the convolution variables z and ± . The singular contributions in the rapidity regulator cancel in the sum of the second and third term of the factorization formula, but this requires that for z → 0 (or 1) these two terms must have closely related structures to all orders of perturbation theory. Indeed, we have shown that this condition gives rise to the D-dimensional "refactorization conditions" which must hold to all orders of perturbation theory. The symbol [[f (z)]] means that one retains only the leading terms of a function f (z) in the limit z → 0 and neglects higher power corrections. We have rewritten the original function H , (2.8) where the new functionH Compared with [9] we have rewritten the second term in a different but equivalent way. The subtraction terms involving the [[. . . ]] symbol remove the singularities at the endpoints z = 0 and z = 1, such that the limit δ → 0 is smooth. Note that both the matching coefficient H 2 (z) contain terms that are singular for z = 0, 1 and the two subtraction terms properly remove the singularities of the product of these two quantities. This generalizes a simple "plus-type" subtraction prescription for the bare operator proposed in [12,29], which works only for cases where the relevant matching coefficient approaches a constant plus power-suppressed terms as z → 0.  We stress that the form of (2.9) is independent of the particular rapidity regularization scheme used to regularize the divergent convolution integrals in (2.1) and (2.6), as long as the same regularization scheme is applied consistently to all terms in the factorization theorem. The refactorization conditions ensure that the integrands of the second and third term in the factorization formula are identical in the singular regions (up to powersuppressed terms), and hence the endpoint divergences can be removed by a rearrangement of these terms, which leads to (2.9).
Removing the endpoint divergences in the way described above comes at the price of introducing hard upper limits on the integrals over + and − in the last term of the factorization formula (2.9), which originally have power counting ± = O(m b ). 1 This gives rise to additional large rapidity logarithms in the matrix element of the operator O 3 . They are a consequence of the so-called collinear anomaly, which results from the fact that a classical symmetry of the effective theory SCET II under rescalings of the light-cone vectors n 1 and n 2 is broken by quantum effects [26]. The presence of the upper limits also leads to some power-suppressed contributions of O(m 2 b /M 2 h ) to the third term, which should be dropped for consistency. Moreover, to obtain the correct result for the decay amplitude the upper limit on the (positive) integration variable + must be analytically continued from M h to −M h −i0 after the integral has been evaluated, as indicated by the limit σ → −1. In figure 2 we show a graphical illustration of the rearrangement of the factorization formula that eliminates the endpoint divergences. The subtractions performed on the second term remove the shaded gray regions from the integrals in the third term. Note that in this process the hard region in which | ± | ≥ M h is subtracted twice. This over-subtraction needs to be corrected by adding back the "infinity bin" in the form of a contribution It is an open question whether it is possible to formulate an alternative "endpoint regularization scheme" avoiding the hard cutoffs. If it exists, such a scheme would need to commute with the operation of renormalization in dimensional regularization, and hence in particular it would need to respect gauge invariance.

JHEP01(2021)077
to the matching coefficient of the operator O (0) ∞ (w) denotes the soft function in the limit where w/m 2 b,0 → ∞, as given in relation (4.46) below. In [9] we have presented explicit expressions for all quantities appearing in the factorization formula (2.9) at next-to-leading order (NLO) in α s , corresponding to two-loop order for the decay amplitude. For completeness, the corresponding expressions are collected in appendix A. The main goal of the present work is to turn the bare factorization theorem into a formula involving renormalized matching coefficients and operator matrix elements. As we shall see this is a highly non-trivial task. The resulting formula provides the basis for a systematic resummation of the large logarithms ln(M 2 h /m 2 b )−iπ to all orders of perturbation theory.

SCET derivation of the refactorization conditions
In this section we derive the refactorization conditions (2.7) using methods from SCET. This discussion is more technical than that in the remaining sections and draws heavily on SCET jargon as well as notations introduced and results derived in [9]. Our arguments should be relatively easy to follow for SCET practitioners. Readers not interested in a technical proof of the relations (2.7) can skip this discussion and proceed directly with section 4. The bare matching coefficient H (0) 2 (z) has been calculated in [9] by computing the on-shell h → b(z k 1 )b(z k 1 ) γ(k 2 ) amplitude for the decay of a Higgs boson into a pair of n 1 -collinear bottom quarks and an n 2 -collinear photon. Here and below we sometime use the symbol z ≡ 1 − z for brevity. To simplify the matching calculation one sets the b-quark mass to zero and assigns momentazk 1 and zk 1 to the outgoing quark and anti-quark, respectively, where k 2 1 = 0. Then the only relevant momentum invariant is hard, 2k 1 · k 2 = M 2 h . In the absence of any non-zero low-energy scale the matrix element of the bare operator O (0) 2 is given by its tree-level expression. Hence, one finds that

Refactorization condition for [[H
where in the last step we have used thatn 1 · k 1 = M h . By calculating the same amplitude in the full theory and comparing the answer with this expression we have derived the result for the bare matching coefficient H 2 (z) given in (A.9) of appendix A. In SCET, the momenta of n 1 -collinear particles in the basis (n 1 · p, n 2 · p, p ⊥ ) have the generic scaling (λ 2 , 1, λ) in units of the hard scale M h . In the case at hand the large O(1) components of the b-quark momenta are equal to z andz. Now consider the limit where z ∼ λ 1. Then the scaling of the outgoing anti-quark momentum becomes soft, (λ 2 , λ, λ). The process is now characterized by two different scales: the hard scale  order contribution to the decay amplitude in this limit can be written in the form The first operator in the time-ordered product describes the decay of the Higgs boson into an n 1 -hard-collinear quark and an n 2 -hard-collinear anti-quark. Hard matching corrections to this vertex are accounted for by the coefficient H (0) 3 . The second operator is an insertion of the subleading-power SCET Lagrangian given in (2.5), which couples a soft quark to an n 2 -hard-collinear quark. There are no hard matching corrections to this Lagrangian [3].
We now decouple soft gluons from the hard-collinear fields by performing the usual field redefinitions [1], but we do not change the names of the fields for simplicity. This leads to In this matrix element the different types of fields (n 1 -hard-collinear, n 2 -hard-collinear and soft) no longer interact with one another. We now match this result onto the low-energy effective theory called SCET II by integrating out the hard-collinear fields. In this step we use the definition of the bare jet function J (0) (p 2 ) given in appendix B to obtain (using that n 2 =n 1 and n 1 =n 2 ) where k 1− =n 1 · k 1 n µ 1 2 , and hence (k 2 + zk 1− ) 2 = zn 1 · k 1n2 · k 2 = zM 2 h . A graphical illustration of this result is shown in figure 3. Note the important fact that, since we perform the calculation on-shell and for massless quarks, there is no soft or n 1 -collinear scale in the problem, and hence the soft and n 1 -collinear matrix elements are equal to their tree-level expressions. In particular, the soft Wilson lines do not give rise to any non-trivial contributions, and the soft matrix element simply provides a factor e izk 1 ·y + v(zk 1 ). Matching this result with (3.1) we obtain where we have used (2.8) in the first step. This establishes the first relation in (2.7).

Refactorization condition for [[ O
We now replace the n 1 -collinear field X n 1 by a soft quark field q s and perform the soft decoupling transformation. This leads to (using thatn 1 = n 2 ) The structure of the soft Wilson lines follows from the fact that the operator on the righthand side derives from the amplitude in (3.3) after integrating out the n 2 -hard-collinear fields. We now need to evaluate the on-shell h → γγ matrix element of this operator. To this end, we need an insertion of the subleading-power SCET Lagrangian, which turns the soft quark field back into an n 1 -hard-collinear quark field. We thus obtain of the jet function and of the soft-quark soft function collected in appendix B. Taking into account a minus sign from an odd number of interchanges of fermion fields, we find where + = n 1 · and − = n 2 · , and the function S (0) ( + − ) is defined as (3.10) Figure 4 shows a graphical representation of the result (3.9). We now switch back to momentum space and define We then find where we have relabeled the integration variable + → − + for later convenience. The singularities of the integrand in the complex + plane are shown in figure 5. There is a pole at + = i0 from the denominator, a cut infinitesimally above the real axis for negative values of + from the jet function, and a cut infinitesimally below the real axis for positive values of + from the soft function. The integral is thus non-zero only if + ≥ 0 and we can deform the contour such that it wraps around the cut in the lower half-plane. This leads to

Renormalized factorization formula
The main goal of this work is to establish the renormalized factorization formula which is structurally equivalent to the bare formula (2.9). We have omitted the external states in the matrix elements for brevity. It is not at all evident that such a formula exists, because, as we will show, the presence of cutoffs on some of the integrals does not commute with the operation of renormalization. We will show that this complication gives rise to some additional contributions, which to all orders of perturbation theory can be absorbed into the definition of the matching coefficient H 1 (µ).

Parameter renormalization
In a first step, we relate the bare parameters entering the decay amplitude M b to the corresponding renormalized parameters. These are the mass of the b quark (which enters in the matrix elements of the operators O i ), its Yukawa coupling (which enters in the expressions for the matching coefficients H i ), as well as the gauge couplings of QCD and QED. We write the relevant renormalization conditions in the MS subtraction scheme as The factor of µ from the renormalization of the Yukawa coupling multiplies the entire decay amplitude. It can be ignored, since after parameter renormalization the amplitude is finite and the limit → 0 is smooth. In our analysis we consider QCD radiative corrections only. To first order in α s ≡ α s (µ), we then have and Z α = 1. Here β 0 = 11 3 C A − 4 3 T F n f is the first coefficient of the QCD β-function, with n f = 5 being the number of active quark flavors.

JHEP01(2021)077
W.F.R. + + + Figure 6. One-loop diagrams contributing to the calculation of the diagonal renormalization factor Z 22 . In the second and third graph the gluon is emitted from the Wilson lines contained in the definition of the collinear quark fields. These graphs must be supplemented by wave-function renormalization.
In our analysis we will sometimes use the b-quark pole mass m b instead of the running mass m b (µ). At NNLO the relation between the two quantities is given by [30] , and for the purposes of this work we do not need the scaleindependent two-loop contribution denoted by the dots. This relation, as well as the relations in (4.3), are known to very high orders of perturbation theory.

Operator renormalization
The matrix elements of the bare operators O (0) 1,2 as well as the bare jet and soft functions J (0) and S (0) contain ultraviolet (UV) divergences not eliminated by the renormalization of the bare parameters. These divergences must be removed by renormalizing the operators themselves, allowing for the possibility of operator mixing. In recent work we have studied the renormalization properties of the jet function [31] and the soft function [32]. We now discuss the renormalization of the remaining operators at first order in α s . We find that the operators O 1 and O 2 mix under renormalization, in such a way that From the definition of the operator O 1 in (2.2) it follows that The relation Z 11 = Z −1 m holds to all orders of QCD perturbation theory, since the quantum fields in the definition of O 1 do not carry color. This factor is known to very high orders of perturbation theory.

JHEP01(2021)077
The diagonal renormalization factor Z 22 can be derived by studying the UV divergences of the matrix element of the operator O 2 (z) defined in (2.3) between an initial-state Higgs boson and a final state consisting of a photon with momentum k 2 and a pair of n 1 -collinear bottom quarks sharing the total momentum k 1 . The relevant one-loop diagrams are shown in figure 6. From a straightforward calculation, we obtain at NLO in QCD where z, z ∈ [0, 1], and the plus distribution is defined in the usual way. The above result is the well-known Brodsky-Lepage kernel [33,34]. This is not surprising, because the colored fields in the operator O 2 in (2.3) have the same structure as the quark fields entering the definition of the leading-twist light-cone distribution amplitude of a transversely polarized vector meson. Before Fourier transformation, they form the non-local structure [9] where [0, tn 1 ] denotes a finite-length Wilson line along the light-like directionn 1 , and ψ(x) is the conventional quark field in QCD. The two-loop expression for Z 22 can in principle be derived from results available in the literature [35][36][37][38][39], but it is not needed for our purposes in this work. It will be useful to rewrite the above result in the equivalent form (4.9) The operator O 2 is not only renormalized multiplicatively, but it mixes with O 1 under renormalization. The factor Z 21 (z) can be derived from the condition that the h → γγ matrix element of the renormalized operator O 2 (µ) in (4.5) be free of UV divergences. Starting from the expression for the matrix element of the bare operator O (0) 2 (z) presented in [9] and shown explicitly in (A.2), applying the factor Z 22 to this expression and renormalizing the bare parameters α b,0 , α s,0 and m b,0 , we find that some 1/ poles remain, which must be removed by the counterterm Z 21 O (0) 1 . In this way we obtain can be obtained in two equivalent ways: either by taking the limit z → 0 in the expression for (4.11) -14 -

JHEP01(2021)077
The renormalization factor [[Z 22 (z, z )]] can be obtained from (4.9) by taking the limit z, z → 0 (with z ∼ z ), leading to Notice the important fact that, contrary to (4.5), the convolution integral now runs over the interval z ∈ [0, ∞). The reason is that z and z are treated as variables which take infinitesimally small values. According to (2.4) the operator O 3 is defined in terms of a time-ordered product of the scalar current J S = hX n 1 X n 2 with two subleading-power Lagrangian insertions. In this product only the current gets renormalized. The corresponding renormalization factor is known to three-loop order and given in eq. (A.2) of [41]. At first order in α s it reads where Here and below we use the implicit definition −M 2 h ≡ −M 2 h − i0 to fix the sign of the imaginary part of the logarithm. When the h → γγ matrix element of O 3 is expressed in terms of a convolution over jet and soft functions, as shown in (2.6), relation (4.14) implies a relation between Z 33 and the renormalization factors Z J and Z S of the jet function and the soft function, which has been given in [32]. It can be written as which despite appearance is independent of the choice of + . This result implies a simple relation between the anomalous dimensions of the matching coefficient H 3 and of the jet and soft functions, which was conjectured in [32] and will be given in relation (5.14) in section 5.1. In general, it is well known that time-ordered products such as in (2.4) can mix under renormalization with local operators of the same order in power counting (see e.g. [40] for a discussion in the context of SCET). In our case, however, the presence of the UV cutoffs on the integrals over + and − in (2.9) and (4.1) prevents such a mixing. The renormalized matrix element of O 3 , which depends on several different physical scales, is expressed in (4.1) as a double convolution over renormalized jet and soft functions. The calculation of the renormalized jet function J(p 2 ) at two-loop order and the study of its RG evolution equation have been discussed in [31], while the renormalization of the soft function S(w) at one-loop order and the derivation of its two-loop RG equation have been studied in [32]. When the JHEP01(2021)077 renormalized expressions are used, we find that (omitting the limit σ → −1 for brevity) It is tempting to interpret this extra contribution as a mixing of the operator O 3 with O 1 , but in reality its origin lies in the fact that imposing the upper cutoffs on the convolution integrals over + and − does not commute with renormalization. It is thus more appropriate to treat the extra term as a contribution to the bare matching coefficient H 1 . We will come back to this issue in section 4.5.

Renormalized matrix elements
With the renormalization factors fixed, we now proceed to derive the h → γγ matrix elements of the renormalized operators in the factorization formula (4.1). For the case of O 1 we trivially obtain For the matrix element of O 2 we find where (4.20) Note that we use the running b-quark mass m b (µ) in the prefactor but the pole mass m b in the argument of the logarithm L m = ln(m 2 b /µ 2 ). Besides being convenient this is not unnatural, because the linear factor of m b (µ) in each matrix element plays the role of a running coupling, whereas the quantity m b appearing in the arguments of the logarithm L m is due to phase-space effects. If desired, one can always switch back from one choice to JHEP01(2021)077 the other using relation (4.4). In the limit z → 0 the above expression simplifies to Note that the same result is obtained using relation in (4.11) along with the renormalization factors given in (4.12) and (4.13).
To obtain the renormalized matrix element of O 3 , we start from the expressions for the renormalized jet and soft functions. They are [25,31] and [32] Hereŵ = w/m 2 b and L w = ln(w/µ 2 ). The result for the function S a (w, µ) takes this relatively simple form only if one uses the pole mass in the argument of the θ(w − m 2 b ) distribution in (4.23). When these expressions are used in the double convolution integral shown in the third term of (4.1), one obtains In other examples where the collinear anomaly appears, the rapidity logarithms take on a simpler form and (typically) exponentiate [26]. In the present case their structure is more complicated, because the rapidity logarithms arise from a double integral over a rather complicated integrand. In order to resum these logarithms it is necessary to factorize the matrix element into a convolution over jet and soft functions, each of which depends on a different scale, and then solve the RG evolution equations of these various functions.

Renormalized matching coefficients
Ignoring the cutoffs on the convolutions integrals in (4.1) for a moment, one would conclude that the UV divergences of the bare matching coefficients H (0) i are removed by applying the inverse matrix of renormalization factors from the right to the row vector of bare matching coefficients (H 3 ), where Specifically, one would then derive as well as At O(α s ) most elements of the inverse matrix Z −1 can be obtained from the corresponding renormalization factors Z ij by simply changing the sign in front of α s . The only exception are the entries Z −1 21 and [[Z −1 21 ]], for which we find The three relations in (4.28) indeed provide the correct renormalization conditions for the corresponding matching coefficients. Using the expressions for the bare matching JHEP01(2021)077 coefficients derived in [9] and collected in (A.9) and (A.10), we find In order to obtain a well-behaved expression we need to restrict the integration to the interval z ∈ [0, 1], like in the first term. We thus define where the sum of the three terms in the second line is now well defined and free of endpoint singularities, such that the limit δ → 0 is smooth. The quantity δ H (0) 1 accounts for the mismatch of integration limits, which one encounters when equating the factorization formula (4.1) expressed in terms of renormalized quantities with formula (2.1) expressed in terms of bare quantities (recall that both correctly reproduce the decay amplitude). After a straightforward calculation we find that δ H (4.33) Starting from the expressions for the bare quantities derived in [9] and given in appendix A, and using our results for the various renormalization factors, we find that δ H . It is then straightforward to obtain from (4.32) (4.34) JHEP01(2021)077

Higher-order analysis of cutoff effects
In our discussion so far we have glanced over an important subtlety related to the cutoffs on the various terms in the bare and renormalized factorization theorems (2.9) and (4.1). In (4.32) the renormalized matching coefficient H 1 (µ) is expressed in terms of bare quantities and renormalization factors. However, it is far from obvious that the sum of the terms on the right-hand side is indeed only a function of the hard scale −M 2 h and independent of the soft scale set by the b-quark mass. Indeed, the definition of δ H  .23) and (A.6). However, we will now show that the sum which enters in (4.32), is independent of the b-quark mass to all orders of perturbation theory. This combined quantity is thus truly a hard subtraction term. From the definition (4.16) it follows that where the factor m b,0 on the left-hand side stems from the matrix element O 1 . Throughout this section we do not write out the limit σ → −1 explicitly, but it is understood in all expressions where σ occurs. Rewriting the renormalized jet and soft functions in the first line in terms of the corresponding bare functions, using relations derived in [31,32] and expressing the renormalized matching coefficient H 3 (µ) in terms of the bare one using the last relation in (4.28), the right-hand side of this equation can be put in the form (4.37) The quantity Z J denotes the renormalization factor of the jet function defined as At one-loop order one finds [25,31] is the symmetric Lange-Neubert kernel [42]. Note that the quantity Z J satisfies the symmetry relation As shown in (4.15), the renormalization factor of the soft function can be expressed in terms of the same object [32]. If it was not for the cutoffs on the integrals, the quantity in (4.37) would evaluate to zero, because the integrals over the products of Z J factors in the second and third lines would yield δ(ρ ± − ω ± ).
This allows us to express the product of renormalization factors in the second line of (4.33) in terms of Z J and Z −1 J . We find The two terms in this result are structurally similar to those appearing in (4.37). Indeed, it is possible to rearrange these terms in such a way that For the first integral, the variables ± are in the hard region, and this forces the variables ρ ± and ω ± to be in the hard region as well (see appendix C for more details). In fact, we can recast relation (4.44) in the alternative form where in the first term the dimensional regulator = (4 − D)/2 must be kept in place after renormalization in order to regularize the divergences for ρ ± → ∞, as indicated by the superscript "( )". This form makes it explicit that the quantity δH (4.46) with C 1 ( ) given in (A.8). In order to calculate the quantity δH (0), tot 1 at leading order in α s we use the lowest-order expressions for the bare jet and soft functions. We also note that the integral in the last line of (4.44) vanishes at first order in α s . After a straightforward calculation we find 47) where H( ) = ψ(1 + ) + γ E is the harmonic-number function. This generalizes relation (4.17) to higher orders in .
In terms of this quantity, the correct all-order definition of the renormalized matching coefficient H 1 (µ) is obtained as

Contributions to the decay amplitude
As a cross check we evaluate the three terms T i shown in the three lines of the renormalized factorization theorem (4.1) using the results for the matrix elements and matching JHEP01(2021)077 coefficients obtained above. This yields where we have defined Adding up the three contributions we correctly reproduce the QCD amplitude up to higher-order corrections in α s and m 2 b /M 2 h .

RG evolution equations
In the previous section we have accomplished the main goal of this work: the establishment of the renormalized factorization formula (4.1). We have derived explicit expressions for the renormalized matrix elements and matching coefficients to first order in α s (corresponding to the two-loop order for the decay amplitude), and we have shown that to all orders of perturbation theory the effects of the upper cutoffs on the convolution integrals over + and − (and on the z integral over the functions inside the braces [[. . . ]]) is not in conflict with factorization. We will now derive the RG evolution equations for the various objects in (4.1), which set the basis for the systematic resummation of the large logarithms in the decay amplitude.
In terms of the various renormalization factors Z ij derived in section 4.2 we obtain the corresponding anomalous dimensions γ ij in the usual way, i.e.

JHEP01(2021)077
where Z (1) ij denotes the coefficient of the single 1/ pole in Z ij . In this way, we obtain the diagonal elements as well as the off-diagonal elements Note that both Z ij and γ ij are scale-dependent quantities, but we suppress this dependence for the sake of simplicity of the notation.
The diagonal elements of the anomalous-dimension matrix are also known to higher orders in α s . Relation (4.6) implies that γ 11 = −γ m is determined in terms of the anomalous dimension of the quark mass, defined as The quantity γ m is known to five-loop order [43]. The two-loop expression for γ 22 (and with it [[γ 22 ]]) can in principle be derived from [35][36][37][38][39]. The anomalous dimension γ 33 , which according to (4.14) is equal to the anomalous dimension of two-jet current operators in SCET, can to all orders be written in the form [41] where Γ cusp is the light-like cusp anomalous dimension in the fundamental representation of SU(N c ) [44], and γ q is the anomalous dimension of the quark field in light-cone gauge. The cusp anomalous dimension has recently been calculated to four-loop order [45], while γ q is known to three loops. It can be determined from the three-loop expression for the divergent part of the on-shell quark form factor in QCD [46,47].

RG equations for the operator matrix elements
From the renormalization conditions (4.5) and (4.11) it follows that the matrix elements of the renormalized operators satisfy the RG evolution equations We have checked that these equations are satisfied to O(α s ).
As mentioned earlier, in order to resum all large logarithms contained in the matrix element of the operator O 3 one must factorize the matrix element in the form leading power (5.7) and solve the RG equations for the jet and soft functions separately. We have derived the corresponding evolution equations at two-loop order in two recent papers. For the jet function one finds [25,31] which in this form holds for both space-like and time-like values of p 2 . The anomalous dimension is given by where The local terms (with x = 1) can to all orders be expressed in terms of the cusp anomalous dimension and an anomalous dimension γ (α s ), which was recently obtained at two-loop order [31]. Since the plus distribution contained in Γ(1, x) is linked with the logarithmic term, it is also multiplied by Γ cusp . However, starting at two-loop order additional nonlocal terms arise, whose explicit form was obtained in [31] by using the RG invariance of the B − → γ l −ν decay rate along with the calculation of the two-loop anomalous dimension of the B-meson light-cone distribution amplitude (LCDA) performed in [48]. The RG equation for the soft function is tightly linked to that of the jet function [32].
with γ s (α s ) = 2γ q (α s ) + 2γ (α s ) . (5.13) Via this relation the quantity γ s is known to two-loop order. As defined above, the anomalous dimensions of the matching coefficient H 3 and of the jet and soft functions obey the simple relation which is a consequence of relation (4.15). In this form it is easy to see that the variable + drops out from the result.

RG equations for the matching coefficients
The renormalized matching coefficients obey the evolution equations where the quantity in the first equation results from non-trivial effects of the cutoffs on the scale evolution. Once again, we have checked that all of these equations are satisfied to O(α s ). The evolution equation for the matching coefficient H 1 (µ) calls for a more careful discussion. We have seen in section 4.5 that the definition of this coefficient in higher orders is quite subtle and requires a careful treatment of the effects of the cutoffs on the various convolution integrals. In this section we will derive the evolution equation for H 1 beyond one-loop order using the RG invariance of the decay amplitude on the left-hand side of (4.1). Explicitly, we evaluate

JHEP01(2021)077
where T 2 and T 3 denote the second and third lines in the factorization formula (4.1). The scale dependence of the third term has been studied in our recent work [32], where we have shown that where the kernel K(x) is given by This quantity appears in the non-local terms of the anomalous dimensions γ J and γ S in (5.9) and (5.12). The two terms on the right-hand side of relation (5.18) are, in fact, identical, as can be seen by redefining the integration variables according to ± → σ ∓ . Note the important fact that the result is not given by a hard function. Indeed, at leading order in perturbation theory the last integral evaluates to 20) which is sensitive to the low scale m b . Using the evolution equations (5.6) and (5.15) it is straightforward to show that the scale dependence of the second term is given by

JHEP01(2021)077
we obtain for the left-over terms (5.24) The refactorization theorem for the bare operator O (5.25) suggests that this result closely resembles the structure of (5.18). We would thus like to establish a similar relation that holds after renormalization. However, the integral on the right-hand side contains endpoint divergences, which are regularized by the dimensional regulator . In other words, the matrix element of the bare operator [[ O (0) 2 (z) ]] contains some 1/ poles not contained in the bare jet and soft functions. Their presence makes the derivation of the renormalized relation non-trivial. After a careful analysis we find that All integration variables are in the hard region, and hence the bare soft function can be replaced by its asymptotic form S (0) ∞ , see (4.46). The two terms on the right-hand side of this relation are UV divergent, but by construction their sum is finite when expressed in terms of renormalized parameters. After a somewhat cumbersome calculation, we find Like the renormalized matrix element on the left-hand side of (5.26), the quantity ∆ 21 contains single-logarithmic terms of order α b α n s L n+1 h in higher orders of perturbation theory. Combining the results (5.18) and (5.24), and using relation (5.26), all terms involving the soft function cancel out, and we obtain the exact formula
In the second step we have used the renormalized refactorization condition (5.23). It is now explicit that this quantity depends only on the hard scale L h . Performing the integrals over x and z, and using the explicit expressions for the quantities K(x) and ∆ 21 given above, we can compute the first two expansion coefficients of D cut (µ) in powers of α s . We find Note also that instead of calculating D cut directly one can recast the second relation in (5.30) in the form where The two-loop coefficient is straightforward to calculate and contains polylogarithms of fourth order. An advantage of the form (5.34) is that it brings the evolution equation for H 1 in (5.15) to a more canonical form. However, one finds that the new anomalous JHEP01(2021)077 dimension γ cut α b (α s L h ) n contains higher powers of the logarithm L h in higher orders of perturbation theory, and it is therefore not of the Sudakov type. This logarithmic behavior in higher orders has an important implication for the solution of the RG evolution equations. As long as it is not known how to resum the logarithmic terms in γ cut or D cut , it is impossible to systematically integrate the evolution equation for H 1 (µ) from a high matching scale down to a scale µ M h . We are thus forced to choose a value for the factorization scale µ that is of order the Higgs-boson mass. The challenge is then to solve the evolution equations for the operator matrix elements in (5.6), (5.8) and (5.11) in order to evolve these matrix elements up to the same scale µ ∼ M h .
Let us mention an interesting and indeed fortunate coincidence in this context. In [32] it has been found that the solution of the RG evolution equation for the soft function S(w, µ) requires that the factorization scale µ must be larger than the matching scale µ s , at which the initial condition for the soft function can be calculated in fixed-order perturbation theory. Because of the fact that the argument w = + − of the soft function is integrated from the soft region (w ∼ m 2 b ) up into the hard region (w ∼ M 2 h ), it is necessary to perform a dynamical scale setting µ 2 s ∼ + − under the integral when solving the RG equation [22]. Hence, it is necessary, also from this point of view, that the factorization scale µ is chosen of order the hard scale M h , so as to ensure that µ > µ s for all values of the integrand. As a final comment, we stress that for the method of dynamical scale setting to be consistent, it is crucially important that the renormalized soft function S(w, µ) in (4.23) does not develop any large logarithms ln n (w/m 2 b ) in the limit where w m 2 b . Otherwise, it would be necessary the refactorize the soft function in the limit where the variable w approaches the hard region.

Large logarithms in the three-loop decay amplitude
The RG evolution equations established in the previous section, along with the explicit expressions for the relevant anomalous dimensions, provide the basis for a systematic resummation of the large logarithms L = ln(−M 2 h /m 2 b ) in the b-quark induced h → γγ decay amplitude. In the factorization theorem (4.1) contribution from the last term is enhanced by at least two powers of logarithms with respect to the other terms, because the integrals over + and − provide a logarithmic enhancement. A careful analysis reveals that, for generic values of µ, the first two terms in the factorization formula, T 1 and T 2 , yield terms of O(α b α n s L n+1 ) to the decay amplitude, while the third term, T 3 , yields terms of O(α b α n s L 2n+2 ). In particular, starting at first order in α s the series of leading and subleading logarithms receives contributions from this last term only. For this reason, it is often stated in the literature that the large double-logarithmic corrections arise from the region in which the quark propagator connecting the two photons carries a soft momentum ("soft quark contribution" [15][16][17]).
With the results derived in this paper it is possible to perform a resummation of large logarithms at (almost) NLO in RG-improved perturbation theory, in which all terms enhanced by large logarithms are exponentiated, while contributions not in the exponent are suppressed by powers of α s . This requires the O(α s ) expressions for all matching JHEP01(2021)077 coefficients and matrix elements as well as for the jet and soft functions, each evaluated at its characteristic scale. The various anomalous dimensions in the evolution equations are needed at two-loop order in QCD, while the three-loop expressions are needed for the cusp anomalous dimension and the β-function. With the exception of γ 21 , all these ingredients are known or, in the case of γ 22 , can be derived from existing results.
While conceptually straightforward, performing the resummation at NLO in RGimproved perturbation theory is technically challenging because of the complicated structure of the RG evolution equations for the soft and jet functions, the need to perform a dynamical scale setting for which the matching scales float with the integration variables + and − [32], and the analytic continuation σ → −1 that needs to be performed in T 3 and requires extending the running coupling α s (µ) into the complex plane. We leave a detailed discussion of these technicalities for future work, mentioning however that for the case of T 3 the resummation has been studied at LO in RG-improved perturbation theory in [22]. Instead, here we will use the RG equations derived in section 5 to predict the terms in the three-loop h → γγ decay amplitude that are enhanced by at least three powers of the large logarithm L. Solving the RG evolution equations (5.6), (5.8), (5.11) and (5.15) iteratively in perturbation theory, it is straightforward to derive the necessary higher-order logarithmic contributions to the various operator matrix elements and matching coefficients.

Higher-order logarithms in the matrix elements
The renormalized matrix elements of the operators O i , the jet function and the soft function have been derived in section 4.3 at first non-trivial order in α s . The result for the matrix element of O 1 in (4.18) is exact without any higher-order corrections. For the matrix element of O 2 , we extend relation (4.19) in the form There is no need to derive the higher-order logarithmic terms for the jet function, because the complete two-loop expression for J(p 2 , µ) has been obtained in [31]. We thus turn directly to the case of the soft function. The iterative solution of the RG equation (5.11) involves some rather complicated integrals, which need to be simplified using various identities for polylogarithms. We present the results as higher-order contributions to the coefficient functions S a and S b defined in (4.23). Recall that in this equation we use the running b-quark mass in the prefactor but the pole mass everywhere else. We JHEP01(2021)077 parameterize our results in the form where the functions s ia (ŵ) vanish in the limitŵ → ∞, while the functions s ib (ŵ) vanish forŵ → 0. For the coefficients r i we obtain T F n f , while the functions s ia (ŵ) are given by In the derivation of these results one needs the two-loop coefficients of the anomalous dimensions Γ cusp , γ s and γ m , which are collected in appendix D. Similarly, for the functions JHEP01(2021)077 − 12 lnŵ ln(1 −ŵ) + 24 lnŵ ln 2 (1 −ŵ) + 8 ln 2ŵ ln(1 −ŵ) To compute the remaining O(α 2 s L 0 w,m ) terms in the soft function would require a complete two-loop calculation.

Higher-order logarithms in the matching coefficients
The renormalized matching coefficients have been discussed in section 4.4. Higher-order corrections to the coefficients H 2 and H 3 can be derived straightforwardly by perturbatively solving the corresponding evolution equations in (5.15). We extend the first relation in (4.31) in the form where f 2 (z) = 2C F ln 2 z − β 0 ln z . (6.8) To derive the terms of O(α 2 s L h ) would require knowledge of the two-loop coefficient of the anomalous dimension γ 22 in (5.6). The matching coefficient [[H 2 (z, µ)]] can be readily derived by taking the limit z → 0 in the above expression.
The second relation in (4.31) can be extended as

9) with
T F n f , T F n f .

JHEP01(2021)077
In deriving these coefficients we have used the two-loop expression for the anomalous dimension γ 33 in (5.5), which is given in appendix D. To determine the O(α 2 s L 0 h ) contribution would require a complete two-loop calculation of H 3 .
Given the higher-order logarithmic corrections to H 2 shown in (6.7), it is straightforward to integrate the first evolution equation in (5.15) perturbatively. In this way we obtain where The derivation of the O(α b α 2 s L 2 h ) term would require knowledge of the O(α 2 s L h ) contribution to H 2 (z, µ), which in turn needs the two-loop coefficient of the anomalous dimension γ 22 .

Higher-order logarithmic contributions to the decay amplitude
Given the above higher-order results for the matrix elements and matching coefficients, it is straightforward to derive the higher-order logarithmic corrections to the b-quark induced h → γγ decay amplitude at O(α 2 s L k ) with k ≥ 3. We find where the dots refer to three-loop terms containing less than three powers of logarithms (L or L m ). The higher-order expansion coefficients are T F n f . (6.14)

JHEP01(2021)077
The amplitude is scale-independent to the order we are working, meaning that the µ dependence of the running couplings y b (µ), m b (µ) contained in M 0 (µ) and of α s (µ) is compensated by the terms containing L m = ln(m 2 b /µ 2 ). As a cross check of our results, we now compare expression (6.13) for the decay amplitude with the results of previous calculations. To this end we need to perform transformations to different renormalization schemes. First, we express the running parameters m b (µ) and y b (µ) = √ 2 m b (µ)/v in the prefactor M 0 (µ) in terms of the b-quark pole mass, using relation (4.4). We then eliminate the remaining scale dependence by making the choice in the running coupling α s (µ). In this "on-shell scheme" (OS), we find that the amplitude takes the form T F n f , and the dots refer to terms containing less than three powers of logarithms. The contributions to the decay amplitude of O(α b α 2 s n f ) have been calculated in closed analytic form in [19], and we find full agreement with the results obtained by these authors. Moreover, recently the entire three-loop gg → h amplitude has been calculated in numerical form [20]. The authors of this paper have kindly repeated their calculation for the case of h → γγ and found the following numerical result for the three-loop coefficient inside the rectangular bracket in (6.15) (to much higher accuracy than indicated here): is given by the β-function for n l = 4 massless quarks, i.e. α (n l ) s (µ). Since in our case the massive b-quark is much lighter than the Higgs boson, it is more appropriate to work in a scheme in which one uses the running coupling defined with n f = n l + 1 active quark flavors, as we have done in (6.15). The relevant conversion relation is [49] Performing this scheme transformation we find that relation (6.17) is replaced by where n f = n l + 1 = 5. Our analytic result in (6.15) is in perfect agreement with this expression. We emphasize that the coefficient of the L 3 term is sensitive to the two-loop anomalous dimensions of the jet and soft functions. The observed agreement thus presents a highly non-trivial cross check of our conjecture for the two-loop anomalous dimension of the soft function made in [32]. As a side remark, we stress that the above results indicate that the leading double logarithms generally do not provide the dominant contributions to the decay amplitude. Using m b = 4.8 GeV for the b-quark pole mass and M h = 125.1 GeV for the mass of the Higgs boson, and indicating powers L n with the help of a subscript n, we find numerically is itself a complex number. It follows that in order to obtain reliable results it is essential to perform the resummation of logarithmic terms beyond the leading double-logarithmic approximation using RG-improved perturbation theory, such that all large logarithms are exponentiated.

JHEP01(2021)077
The above comparison provides a highly non-trivial cross check of our factorization theorem (4.1) derived in SCET. Previous authors have analyzed the leading and subleading logarithms in the h → γγ decay amplitude using more traditional resummation techniques [15,16]. They worked in a different renormalization scheme, in which the factor m b (µ) in the prefactor M 0 (µ) is eliminated in favor of the b-quark pole mass, whereas the Yukawa coupling y b (µ) and the running coupling α s (µ) are evaluated at µ = M h . In this scheme OS , our result for the decay amplitude assumes the form The leading double logarithms of O(α b α n s L 2n+2 ) have been correctly obtained in [15,16]. In [16] also the next-to-leading logarithms (NLL) of O(α b α n s L 2n+1 ) were analyzed and an all-order formula for them was proposed. The result reads .
When expanded to O(α 2 s ) this formula yields Interestingly, the coefficient of the C 2 F α 2 s L 5 term (marked in bold face) does not agree with our result shown in (6.21). It is difficult to trace the origin of this discrepancy, given that JHEP01(2021)077 in [16] the derivation of the subleading logarithmic contributions is only sketched. The authors start from an analysis of the off-shell Sudakov form factor [50], i.e. the quark form factor of a vector current in the limit where Q 2 = |q 2 | |p 2 i |, and then take the limit where the two external legs go on-shell (p 2 i → 0). They also need to account for the kinematic differences between quark scattering off a vector current and photon scattering off a scalar (Higgs) current. In general, a consistent framework based on effective field theory, such as the SCET approach developed here, is certainly helpful to derive consistent results for corrections appearing beyond the leading order in both logarithmic counting and power counting in λ = m b /M h . In our approach, the coefficients of the leading and subleading logarithms in (6.21) are determined in terms of one-loop coefficients of anomalous dimensions. We find that 25) where Γ 0 = 4C F , γ q,0 = −3C F , γ 0 = 0 and γ m,0 = −6C F are the one-loop coefficients of the cusp anomalous dimensions and the anomalous dimensions of the collinear quark field, the jet function and the running quark mass, respectively, and we have used relation (5.13) to eliminate the one-loop anomalous dimension γ s,0 of the soft function. In [22] we have extended the prediction of the NLL terms to higher orders of perturbation theory, finding that (6.23) must be replaced by (n + 1) 2 (2n + 3)(2n + 5) , (6.26) where ρ = C F αs(M h ) 2π L 2 , and in the prefactor m b denotes the pole mass. The second term inside the brackets in the second line is not in agreement with (6.23).
In the recent paper [21] the resummation approach of [16] was extended to predict the leading and subleading logarithms in the gg → h amplitude. The authors showed that in an appropriate "abelian limit" their result reproduces the formula for the subleading logarithms shown in (6.23). We thus believe that also in this work the subleading logarithms at three-loop order and beyond are not correctly accounted for. Matching their results with the calculation of [20], the authors have concluded that the coefficient of the L 5 term in equation (C.1) of this paper should take the form (C A − C F ) 11 9 C A − 3 2 C F − 2T F n b /640 ≈ 0.0017361111. Adjusting the color factors in our result (6.21) to those relevant for the gg → h amplitude, we find instead the result

JHEP01(2021)077 7 Conclusions
We have derived the first renormalized factorization theorem for an observable described at subleading order in SCET power counting. We have focused on the contribution to the radiative Higgs-boson decay amplitude induced by the Higgs coupling to light bottom quarks; however, the methods we have developed are more general and can be applied to other subleading-power factorization theorems. Endpoint-divergent convolution integrals arising when the factorized decay amplitude is expressed in terms of bare matching coefficients and bare operator matrix elements have been tamed by introducing rapidity regulators on the convolution integrals. We have proved two D-dimensional refactorization conditions for the bare matching coefficient H 2 (z) in the endpoint region z → 0, which ensure that the dependence on the rapidity regulator cancels out to all orders of perturbation theory. With the help of these relations the factorization formula can be recast into a form where the endpoint divergences are removed by means of suitably chosen subtraction terms (for T 2 ) and cutoffs on the convolution integrals (for T 3 ).
The main accomplishment of the present work has been to show that such an endpointregularized factorization theorem can be consistently formulated in terms of renormalized matching coefficients and operator matrix elements. This is a highly non-trivial point, because endpoint regularization and renormalization do not commute. We have derived the RG evolution equations satisfied by the renormalized matching coefficients and operator matrix elements and derived most of the anomalous dimensions at two-loop order (and the remaining ones at one-loop order). We have then used our results to predict in analytic form the logarithmically enhanced three-loop contributions to the b-quark induced h → γγ decay amplitude of O(α b α 2 s L k ) with k = 6, 5, 4, 3, finding perfect agreement with a numerical computation of these terms performed by the authors of [20]. On the other hand, our findings for the structure of the coefficient of the subleading term (with k = 5) disagrees with the predictions of previous authors [16,21], who had attempted to study the structure of the subleading logarithmic contribution using conventional tools. This demonstrates the usefulness of having a fully systematic approach based on effective field theory to study factorization beyond the leading power in scale ratios.
We are confident that the results presented in this work are a major step forward in the quest for a consistent formulation of SCET factorization theorems at subleading power. For the particular example considered -the b-quark induced h → γγ decay amplitudethey form the theoretical basis for a systematic resummation of large double and single logarithms beyond leading order in RG-improved perturbation theory. The technical details of this resummation will be discussed in future work. It will also be important to generalize our analysis to the non-abelian case of the Higgs-boson production in the gluon-gluon fusion channel gg → h, extending the approach of [17,18] to higher logarithmic accuracy.

JHEP01(2021)077
One of us (M.N.) thanks Gino Isidori, the particle theory group at Zurich University and the Pauli Center for hospitality during a sabbatical stay. This research has been supported by the Cluster of Excellence PRISMA + funded by the German Research Foundation (DFG) within the German Excellence Strategy (Project ID 39083149). The research of Z.L.L. is supported by the U.S. Department of Energy under Contract No. DE-AC52-06NA25396, the LANL/LDRD program and within the framework of the TMD Topical Collaboration.

A Bare matching coefficients and matrix elements
For completeness, we list here the expressions for the h → γγ matrix elements of the bare operators O 2 (z) |h ]] one needs to take the limit z → 0 in the above expressions. In [9] this limit has been obtained in closed form in the dimensional regulator . One finds 3 , as given in the third line of (2.9), one needs the expressions for the bare jet and soft functions at NLO in α s . For the bare jet function one obtains This function has a cut along the positive p 2 axis starting at p 2 = 0 and extending to infinity. The bare soft function, which is defined in terms of the discontinuity of a soft quark propagator dressed with Wilson lines, can be written in the form At first order in α s one finds

JHEP01(2021)077
The first term restricts both ρ − and ω − to be in the hard region. For the second term this is not obvious but still true, because the factor ρ − in front of the plus distribution Γ(ρ − , ω − ) removes the factor 1/ρ − in the measure of the integral. We thus get (focussing only on the integrals over "minus" momenta) In these expressions we can drop the plus prescription, because the integrand of the integral over ω − contains a term θ(M h − ω − ) if we write this integral as ∞ 0 dω − · · · . The plus prescription then gives a subtraction term involving θ(M h − ρ − ), which vanishes since the integral over ρ − runs from M h to infinity. The first integral on the right-hand side of (C.9) is clearly in the hard region. For the second integral ω − must be treated as a hard variable of O(M h ), because the jet function does not contain any reference to the mass of the b quark and ρ − is in the hard region. In other words, the region where ω − = O(m b ) gives rise to a power-suppressed contribution and hence must be dropped. (Recall that we must only keep the leading-power terms in the result.) Finally, for the third integral ω − is in the hard region, and the region where ρ − = O(m b ) or smaller gives rise to a power-suppressed contribution. For this to be true, it is important that the measure is dρ − and not dρ − /ρ − . An analogous argument holds for the second integral over Z J factors in (C.7).
With all integration variables restricted to the hard region, we can replace the bare soft function S (0) (w) by its asymptotic form S (0) ∞ (w) defined in (4.46). When the dust settles, we obtain from (C.7) where at this order we can use the lowest-order expressions for S (0) ∞ (w) given in (4.46). Performing the integrals then leads to (4.47).

JHEP01(2021)077 D Two-loop anomalous dimensions
We define the expansion coefficients of the cusp anomalous dimensions as Γ cusp (α s ) = Γ 0 α s 4π + Γ 1 α s 4π 2 + . . . , (D. 1) and similarly for all other anomalous dimensions. Below we list relevant expansion coefficients needed in section 6 in the MS renormalization scheme. The expansion coefficients of the cusp anomalous dimension Γ cusp are given by [44] Γ 0 = 4C F , For the coefficients of the anomalous dimension of the running quark mass, which is related to the anomalous dimension of the operator O 1 (µ) by γ 11 = −γ m , one finds [30] γ m,0 = −6C F , 3) The anomalous dimension γ q of the collinear quark field entering in (5.5) has coefficients [41,46] γ q,0 = −3C F , (D.4) while the coefficients of the anomalous dimension γ entering in (5.9) read [9] γ 0 = 0 , The anomalous dimension γ s then follows from (5.13). One finds the expansion coefficients [32] γ s,0 = −6C F , (D.6) Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.