Factorization at subleading power and endpoint-divergent convolutions in h → γγ decay

It is by now well known that, at subleading power in scale ratios, factorization theorems for high-energy cross sections and decay amplitudes contain endpoint-divergent convolution integrals. The presence of endpoint divergences hints at a violation of simple scale separation. At the technical level, they indicate an unexpected failure of dimensional regularization and the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} subtraction scheme. In this paper we start a detailed discussion of factorization at subleading power within the framework of soft-collinear effective theory. As a concrete example, we factorize the decay amplitude for the radiative Higgs-boson decay h → γγ mediated by a b-quark loop, for which endpoint-divergent convolution integrals require both dimensional and rapidity regulators. We derive a factorization theorem for the decay amplitude in terms of bare Wilson coefficients and operator matrix elements. We show that endpoint divergences caused by rapidity divergences cancel to all orders of perturbation theory, while endpoint divergences that are regularized dimensionally can be removed by rearranging the terms in the factorization theorem. We use our result to resum the leading double-logarithmic corrections of order αsnln2n+2−Mh2/mb2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha}_s^n{\ln}^{2n+2}\left(-{M}_h^2/{m}_b^2\right) $$\end{document} to the decay amplitude to all orders of perturbation theory.


Introduction
Factorization theorems play an important role in understanding the structure of observables sensitive to several different mass or distance scales. They provide the means to separate long-distance from short-distance phenomena and resum large logarithmic corrections to all orders of perturbation theory. Whereas at leading power in scale ratios a typical factorization theorem consists of a product (often in the convolution sense) of functions each associated with a single scale, at subleading power several complications arise. Yet, in times when high-precision calculations of scattering and decay processes are state-of-the-art, going beyond the leading order in scale ratios has become a necessity in many cases.
For example, obtaining more precise fixed-order calculations of collider cross sections has been an important goal for many years. A major difficulty in these calculations is the proper handling of the infrared singularities that arise in both virtual and real contributions.

JHEP04(2020)033
A method based on N -jettiness (T N ) slicing [1,2] allows one to obtain the next-to-nextto-leading order (NNLO) QCD result from a much easier next-to-leading order (NLO) calculation, combined with information about the singular dependence of the cross section on the T N resolution variable [3]. Calculations of the dominant power corrections in T N /Q have helped to improved the numerical stability for several processes [4][5][6]. In other cases, the quantity of interest may itself be a power-suppressed process. An example is the decay of a heavy scalar particle into a pair of light fermions [7,8].
Soft-collinear effective theory (SCET) [9][10][11][12][13] provides a convenient framework for addressing the problem of factorization and scale separation using the powerful tools of effective field theory. Much recent work has focused on exploring the structure of factorization theorems at subleading order in power counting. In a first step, this involves the construction of a basis of subleading SCET operators (see e.g. [14][15][16][17][18][19]) and the calculation of their anomalous dimensions [7,16,17,20,21]. Specific applications that have been discussed include the threshold resummation at subleading power for the Drell-Yan and Higgs production processes [22][23][24][25][26][27], the study of power corrections to event shapes [28] and transverse-momentum distributions [29,30], and an analysis of power corrections to the exclusive B-meson decay B → γlν l [31]. One finds that such factorization theorems contain a sum over several products of Wilson coefficients times operator matrix elements, where the relevant SCET operators mix under renormalization. At the same time 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. Such endpoint divergences were first observed in the analysis of exclusive, non-leptonic and radiative weak decays of B mesons in the context of QCD factorization [32][33][34][35]. While the corresponding decay amplitudes were shown to factorize at leading order in Λ QCD /m b , some of the power-suppressed corrections were found to contain endpoint-divergent convolution integrals of hard-scattering kernels with meson light-cone distribution amplitudes. These divergences have been interpreted as a breakdown of factorization in the non-perturbative domain. To the best of our knowledge, the first examples of endpoint-divergent convolution integrals in a perturbative setting have been found in the study of the exclusive decay B → χ cJ K [36] and of non-local power corrections to the inclusive radiative decay B → X s γ [37].
It is by now well known that endpoint-divergent convolution integrals are a generic feature of SCET factorization theorems beyond leading power. In simple cases, where the relevant momentum modes are hard, collinear and ultrasoft and the relevant effective theory is the so-called SCET-1, the endpoint divergences are regularized by the regulator = (4 − D)/2 of dimensional regularization. Nevertheless, these divergences indicate an unexpected failure of dimensional regularization and the MS subtraction scheme, because the 1/ n pole terms arising in integrals over products of functions depending on different mass scales cannot be removed by means of the renormalization factors of these individual functions. Hence, a simple scale separation is spoiled by these divergences. Standard tools of renormalization theory are then no longer sufficient to obtain a well-defined, renormalized factorization theorem involving convergent convolutions over renormalized functions. In more complicated cases, in which the relevant momentum modes are hard, collinear and JHEP04(2020)033 soft and the relevant effective theory is the so-called SCET-2, one encounters the additional complication that there exist endpoint divergences which are not regularized by the dimensional regulator but require an additional rapidity regulator.
In this work, we consider the b-quark induced amplitude for the decay of a Higgs boson into two photons, h → (bb) * → γγ, as a concrete example of a quantity obeying a subleading-power factorization theorem. The power suppression arises from the need for an insertion of the b-quark mass, which implies a suppression by one power in the expansion parameter λ = m b /M h 1. (The corresponding amplitude for a heavy vector boson decaying into two light particles would not feature such a suppression.) This power suppression introduces sensitivity to a kinematic region in which a soft quark is exchanged between the two photons. Contrary to soft gluons, soft quarks are power suppressed in SCET, and their emission cannot be described in terms of soft Wilson lines [38,39]. The factorized decay amplitude in SCET is subject to both types of endpoint divergences: those regularized with the dimensional regulator and those regularized by a rapidity regulator. This example therefore illustrates the most complicated situation in which endpoint divergences can arise.
The b-quark induced h → γγ amplitude has been analyzed long ago using conventional tools of perturbative quantum field theory [40,41], and it was shown how to resum the leading double-logarithmic QCD corrections in the ratio M 2 h /m 2 b to all orders of perturbation theory. More recently, the scope of this approach has been extended significantly by validating its validity up to three-loop order for several relevant processes [42][43][44]. Here we present the first SCET analysis of the b-quark induced h → γγ process. Our goal is to establish an all-order factorization theorem for the decay amplitude and develop a framework that will in principle allow for the resummation of large logarithms at arbitrary order in perturbation theory. This paper takes the first steps toward this goal. In section 2 we derive a factorization theorem in terms of products (in the convolution sense) of bare Wilson coefficients with matrix elements of bare SCET operators. In sections 3 and 4 we explicitly calculate the matching coefficients and operator matrix elements at NLO in the QCD coupling α s , corresponding to the two-loop approximation for the decay amplitude. While formulating the proper procedure of MS-like renormalization in the presence of endpoint divergences is still an open challenge (see e.g. the recent works [7,27,39]), we take important steps in the direction of solving this problem. We introduce a suitable analytic regulator for the rapidity divergences in section 5, where we also derive the exact conditions required for the cancellation of rapidity divergences to all orders of perturbation theory. This cancellation is highly non-trivial in our case, since it arises between two terms in the factorization formula depending on different physical scales (hard and collinear versus hard, hard-collinear and soft). In section 6 we succeed to rearrange the factorization formula in such a way that all endpoint divergences -those requiring analytic regulators and those regularized dimensionally -are removed. Our main result, a factorization formula free of endpoint divergences, is presented in (6.1). In section 7 we identify the leading double logarithms contributing to the decay amplitude and resum them to all orders of perturbation theory, carefully taking into account the time-like kinematics of the h → γγ process. We thus resum logarithms of M 2 h /m 2 b along with iπ terms [45,46], and we also include JHEP04(2020)033 the scale evolution of α s (µ), thereby extending the results obtained in [40][41][42][43][44] in several important ways. Finally, in section 8 we summarize our results and discuss some aspects that will be relevant for extending the resummation beyond the leading double-logarithmic approximation.
As this paper deals with a rather sophisticated application of scale factorization, we shall assume that the reader is familiar with the basic formalism of SCET as described, e.g., in the founding papers [9][10][11][12][13].

Bare factorization theorem
In this paper, the h → γγ amplitude is always to be understood as the contribution mediated by virtual b-quarks. We neglect the (numerically dominant) contribution from the top-quark loop, which does not give rise to the large logarithms investigated here. At NLO in QCD perturbation theory, the amplitude for the process h → (bb) * → γγ is given by [47] − iπ is the large logarithm in the problem, and Here α b ≡ e 2 b /(4π) (with e b = −e/3) is the electromagnetic coupling of b-quarks. The running b-quark mass m b , the Yukawa coupling y b and the QCD coupling α s are defined in the MS scheme. The question how to set the renormalization scale in these quantities has been a subject of some discussion (see e.g. [47]), and it is most pressing for the definition of m b , which enters as an argument in the large logarithm L and as a prefactor. The difference between the b-quark pole mass m pole b and the running mass m b (M h ) is almost a factor of 2, and hence scale ambiguities can have a large impact on the decay rate. This issue was particularly pressing in times when low-energy supersymmetry was still a viable option, because in such models the Higgs coupling to b-quarks could be enhanced by a large factor tan β 1, making it numerically comparable with the top-quark induced contribution, which is the leading contribution in the Standard Model (SM). But even in the context of the SM an accurate prediction for the b-quark induced contribution is desirable. The uncertainty associated with scale setting in the b-quark mass and Yukawa coupling is not unimportant in the context of otherwise very precise, higher-order calculations of the h → γγ decay amplitude and the related gg → h production amplitude (see [48] for a review).

JHEP04(2020)033
In order to address this question rigorously, one should factorize the two hierarchical scales M h m b and resum the large logarithms L to all orders of perturbation theory. This resummation has first been accomplished in [40,41] using conventional tools of perturbative QCD, but only at the level of the leading double logarithms α n sL 2n+2 , wherê where 2 F 2 (. . . ) is a generalized hypergeometric function. Note that in this approximation the question of the scale choice is not resolved, because the leading double-logarithmic (LDL) approximation is not a consistent approximation in the sense that it ignores terms that are parametrically larger than O(1).
Analyzing the QCD diagrams giving rise to the decay amplitude at one-and two-loop order using the method of regions [49][50][51], we find that the amplitude receives leading contributions from the following momentum regions of loop momenta: We use a power counting where λ = m b /M h is the small parameter in the construction of the effective field theory. When soft exchanges are present, one also needs regions with hard-collinear scaling: Here n µ 1 and n µ 2 are light-like momenta aligned with the directions k µ 1 and k µ 2 of the external photons, respectively, and satisfying n 1 · n 2 = 2. In the rest frame of the Higgs boson, we can choose n µ 1 = (1, 0, 0, 1) and n µ 2 = (1, 0, 0, −1). Above we have indicated the scalings of the components (n 1 · , n 2 · , ⊥ ) of the loop momentum in the decomposition It will be useful below to introduce also the conjugate vectorsn µ 1 ≡ n µ 2 andn µ 2 ≡ n µ 1 . To illustrate the region analysis, let us consider a one-loop graph and a sample twoloop diagram contributing to the h → γγ decay amplitude in the full theory, see figure 1. We find that the first graph receives leading-power contributions from exactly four regions, in which the loop momentum flowing between the two photon vertices is either hard (h), n 1 -collinear (c), n 2 -collinear (c) or soft (s). Note that in the case where is soft the two remaining propagators carry hard-collinear momenta. The calculation of the individual contributions from the four regions requires rapidity regularization. Using the regularization scheme introduced in section 5, we find The scale µ is associated with the dimensional regulator , while ν is associated with the rapidity regulator η. Adding up the four contributions, we recover the correct result for the diagram, , and (s, s). In section 5.1 we will calculate all contributions from the various regions at two-loop order and show that their sum correctly reproduces the QCD amplitude. This serves as a cross check that we have identified the leading momentum regions correctly.

Matching onto SCET-1
One can describe the transition from QCD to SCET as a two-step process, as illustrated in figure 2. At a hard matching scale of order µ h ∼ M h , we match the SM (after electroweak symmetry breaking and with top-quarks integrated out) onto a version of SCET called SCET-1, which contains hard-collinear fields interacting with soft fields. The relevant momentum regions in this first effective theory are shown, at one-loop order, in figure 3. From this it is straightforward to identify the associated SCET operators. The decay amplitude can be written in the form

9)
JHEP04(2020)033 where some of the products of hard matching coefficients and operator matrix elements are meant in the convolution sense. In the first diagram the loop momentum is hard, and so in the effective theory all propagators are shrunk to a point. The relevant SCET-1 operator O 1 contains the Higgs field coupled to two hard-collinear photon fields. In the second diagram the loop momentum connecting the two photons is n 1 -hard-collinear, and hence the propagator between the second photon (with momentum k 2 ) and the Higgs boson is hard. In the effective theory this propagator is shrunk to a point. The relevant SCET operator O 2,n 1 contains the Higgs field, an n 2 -hard-collinear photon and two n 1hard-collinear quark fields. The latter two fields share the large momentum component n 1 · k 1 of the photon, and hence the operator depends on a variable z ∈ [0, 1] denoting the momentum fraction carried by one of the two hard-collinear fields. Finally, in the third contribution the propagator between the two photons is soft and the two remaining propagators are hard-collinear. While soft gluon interactions with hard-collinear fields are contained in the leading-order SCET Lagrangian [9,10], interactions of soft quarks with hard-collinear quarks and gluons first arise at subleading power in SCET-1, more precisely at O(λ 1 2 ) [12]. The relevant operator O 3 thus contains the time-ordered product of a scalar current made up of two hard-collinear quarks (coupled to the Higgs field) with two insertions of the subleading SCET Lagrangian coupling these quarks to soft quarks.
The relevant SCET-1 operators needed to describe these contributions are  breaking) is denoted by h, whereas fields denoted by calligraphic letters are the so-called "hard-collinear building blocks" of SCET [11,52]. They are composite objects invariant under "n i -hard-collinear gauge transformations", which preserve the scaling of the particle momenta shown in (2.5). These composite fields consist of hard-collinear fields dressed by hard-collinear Wilson lines for QCD and QED interactions, respectively. The n i -hard-collinear photon field is defined as [11,52] where iD µ n i is a covariant derivative containing an n i -hard-collinear gauge field. Note that the building block A µ n i includes a factor of the b-quark electric charge e b in its definition. The n i -hard-collinear b-quark field is defined as n i . Since we work to lowest order in QED interactions and the external photons have transverse polarization, we can however ignore the electromagnetic Wilson line in this definition. The advantage of using hard-collinear building blocks JHEP04(2020)033 is that gauge invariance is explicit despite the fact that SCET is intrinsically non-local (through the appearance of the Wilson lines). Note the important fact that SCET fields living in different hard-collinear sectors (labeled by n 1 and n 2 ) can only interact with each other via the exchange of soft particles.
The operator O 2,n 1 (t) in (2.10) contains two hard-collinear quark fields in the same sector (n 1 ). When this appears, such fields do not need to live at the same spacetime point but can be delocalized along then 1 lightcone [12]. The reason is that the derivative (in 1 · ∂) = O(λ 0 ) is of leading order in SCET power counting and thus can be inserted an arbitrary number of times in such an operator. It is convenient to move from position to momentum space and define (2.14) where The δ-function enforces that the longitudinal momentum componentn 1 · p of the outgoing hard-collinear anti-quark is equal to ω. In SCET these large momentum components of hard-collinear particles are always positive and add up to the total hard-collinear momentumn 1 ·P n 1 in that sector, which in our case is the momentumn 1 ·k 1 of the external photon. We can thus parameterize , provided we also rescale ω → ξ ω. It follows that the h → γγ matrix element of O 2,n 1 (ω) can only be a function of the invariant ratio z 1 = ω/n 1 · k 1 . Hence, the convolution can be expressed in terms of functions H 2 (z 1 ) and γγ| O 2 (z 1 ) |h , whose only dependence on the reference vectors n 1 andn 1 is through their argument z 1 = ω/n 1 ·k 1 . The same functions appear in the convolution of H 2,n 2 (ω) and γγ| O 2,n 2 (ω) |h in the n 2 -hard-collinear sector, but in this case they are evaluated at z 2 = ω/n 2 · k 2 . As a side remark, we note that the photon matrix element of the n 1 -hard-collinear fields contained in O 2,n 1 (t) is related to the leading-twist light-cone distribution amplitude φ γ (u) of the photon [54,55], Here χ m denotes the magnetic susceptibility and b b is the b-quark vacuum condensate. It follows that γγ| O 2 (z) |h ∝ φ γ (z). It can be shown rigorously that in the limit µ → ∞ the renormalized distribution amplitude approaches the asymptotic form We will come back to the significance of this observation in section 8.

JHEP04(2020)033
The power-suppressed Lagrangian insertions in the definition of O 3 describe the coupling of a soft quark to hard-collinear fields at O(λ 1 2 ) in the SCET expansion parameter. In the notation of [12], one has where the soft quark fields q s must be multipole expanded for consistency [12,13], and we denote Here ξ n 1 is a collinear quark spinor subject to the constraint / n i ξ n i = 0, but not yet dressed by a Wilson line. When transformed into hard-collinear building blocks, the corresponding expressions are where A µ n 1 and G µ n 1 are the building blocks for the hard-collinear photon and gluon fields, respectively. G µ n i is defined in analogy with the first relation for A µ n i in (2.12). Note that in the subleading-power Lagrangians the different hard-collinear fields appear at the same spacetime point, unlike the more general case in (2.10). Also, it can be shown that the subleading-power Lagrangians are not renormalized [12].
Let us now analyze the power counting of the three operators in (2.10). With our power counting, the hard-collinear quark, (transverse) photon and gluon fields all scale like λ 1 2 , whereas the soft quark field scales like λ 3 2 . Using that m b ∼ λ, it follows that O 1 ∼ λ 2 and O 3 ∼ λ 2 , whereas it would appear that O 2,n i ∼ λ 3 2 . However, the h → γγ matrix element of O 2,n i requires a mass insertion of m b in order to be non-zero, and an O(λ) quark mass term in the hard-collinear Lagrangian corresponds to a power-suppressed interaction of O(λ 1 2 ). We refrain from writing this interaction in the form of a time-ordered product, since after matching to SCET-2 the quark mass term will be promoted to a leading-order interaction. As a consequence, it follows that O 2,n i ∼ λ 2 has the same scaling as the other operators.
Let us note, for comparison, that the coupling of a heavy vector boson Z to a pair of hard-collinear quarks would be described by the SCET-1 operator Z µX n 1 X n 2 ∼ λ. Hence, it follows that the three operators in (2.10) are of subleading order in SCET power counting. It is only for this reason that we obtain a factorization theorem involving a sum of three types of operators, as indicated in (2.9).
At leading order in the SCET expansion, the hard-collinear quark and gluon fields in the operators in (2.10) can interact with soft gluons (and photons) through eikonal couplings involving the components n i · G s (and n i · A s ) of the soft gauge fields. These couplings can be removed from the Lagrangian by means of the field redefinitions [10]

JHEP04(2020)033
and similarly for the n 2 -hard-collinear fields with x − replaced by x + . Here S n i are soft Wilson lines defined as where we have written out the color (i, j, k, l) and spinor (α, β, γ) indices explicitly and taken into account a factor (−1) arising from the anti-commutators of fermionic fields.

Matching onto SCET-2
Eventually, we are interested in evaluating the h → γγ matrix element in (2.9) using onshell photon states. The various operators and Wilson coefficients in (2.9) can be evolved to scales below the hard matching scale µ h ∼ M h using renormalization-group (RG) evolution equations in SCET-1. At a hard-collinear matching scale of order µ hc ∼ √ M h m b , we match the theory constructed in the previous subsection onto another version of SCET called SCET-2. It contains collinear fields and soft fields, which do not interact with each other. For the operators O 1 and O 2,n i this matching is trivial in the sense that we simply need to replace the various hard-collinear fields with collinear fields, which are defined in a completely analogous way. Note that, as mentioned above, the quark mass now appears in the leading-order SCET Lagrangian. The collinear quark and transverse gauge fields all scale like λ in SCET-2, and hence the operators in (2.9) now scale like λ 3 .
The matching relation for the operator O 3 is non-trivial, because the coupling of a soft quark to a collinear photon produces hard-collinear quanta, which must be integrated out at the hard-collinear scale. This is illustrated in figure 4. In this process two jet functions arise, one in each hard-collinear sector. There is a crucial difference here with regard to the case of electroweak Sudakov resummation at leading power, as studied in [56,57]. The exchange of soft gluons between two hard-collinear quarks is entirely described in terms of emissions from soft Wilson lines, without the appearance of non-trivial matching coefficients. In the present case, however, non-trivial matching coefficients appear (see also [38,39]). Integrating out the hard-collinear fields gives rise to a generalized jet func- tion J, which we define as such that J = 1+O(α s ). Here p denotes the momentum leaving the Higgs vertex at position r, and (p − k 1 ) is the momentum flowing into the vertex at point x. The function J can depend on the two invariant scalar products p 2 and (p − k 1 ) 2 , where k 1 is the momentum of the external photon, and thus k 2 1 = 0. Note that the above definition refers to the massless theory, because the quark mass must be set to zero for the hard-collinear fields.
The corresponding matrix element of the hard-collinear fields in the second line of (2.23) can be obtained from the expression above by hermitian conjugation, exchange of n 1 ↔ n 2 and k 1 ↔ k 2 , and using crossing symmetry (i.e. moving the photon from the initial to the final state). This yields (2.25) We now insert the expressions (2.24) and (2.25) into the matrix element (2.23) and perform the relevant contractions of color and spinor indices. We are then led to define the soft matrix element in the form where the trace in the first line is over color indices, and we have introduced the finite-length JHEP04(2020)033 soft Wilson lines 27) and similarly for the electromagnetic Wilson lines. Here the path ordering is such that the fields appear in the same order as the arguments in the functions. For example, in the expression for S n 1 (x − , 0) the fields are ordered according to t-values, such that fields at t =n 1 · x/2 appear to the left and fields at t = 0 appear to the right. Reparameterization invariance [53] implies that the structure functions S i ( ) can only depend on the invariants n 1 · n 2 · and 2 ⊥ . The prefactor e 2 b /π is included for later convenience. The multipole expansion ensures that only those components of the soft momentum which are commensurate with the relevant hard-collinear momenta enter at the vertices x and y [12,13]. Concretely, the momentum (n 1 · )n µ 1 2 leaves the vertex at x, where the soft quark connects to n 1 -hard-collinear fields, while (n 2 · )n µ 2 2 flows into the vertex at y, where the soft quark connects to n 2 -hard-collinear fields.
When the soft and jet functions are combined according to (2.23), one finds that all structures except for S 1 vanish. Hence, after matching onto SCET-2 the h → γγ matrix element of O 3 assumes the factorized form (omitting the photon polarization vectors for simplicity) The matrix element of the hermitian conjugate operator in (2.23) is obtained by interchanging the subscripts 1 and 2 everywhere. The two matrix elements are related to each other by changing the sign of the loop momentum , which leaves the soft function S 1 ( ) unchanged. Hence, the term (1 ↔ 2) simply leads to a factor 2. Note that the second argument of the function J is zero in both cases, since the relevant momentum is equal to the multipole expanded soft momentum, e.g. (n 1 · )n µ 1 2 for the case of the first jet function, which is light-like. We will from now on only display the first argument for brevity. Finally, because of the multipole expansion, the transverse components of the loop momentum only enter in the soft function. We can thus rewrite the answer in the form (with + ≡ n 1 · and − ≡ n 2 · ) This result clearly shows that the matrix element depends on both the soft and the hardcollinear scale, in accordance with figure 4. In figure 5 we show the analytic structure of the integrand in (2.29) in the complex − plane. The jet function J(−n 2 · k 2 − )/(− − + i0) has a pole at − = i0 and a cut along (and slightly above) the negative real axis starting at the origin. The integrated soft function has a cut starting at + − = −i0 and extending to infinity. For + < 0, this cut is also located above the real axis. The contour can then be closed in the lower half-plane (the dimensional regulator ensures that the contribution from the half-circle at infinity vanishes) and hence the integral yields zero. For + > 0, however, the discontinuities of the soft function lie in the lower half-plane and thus the integral is non-zero. It follows that we can express the matrix element in the form

JHEP04(2020)033
is the discontinuity of the soft function S( + − ). An analogous discussion can be made in the complex + plane. In this case one obtains (2.31) withn 1 · k 1 ↔n 2 · k 2 and + ↔ − interchanged, i.e. with different signs in the arguments of the jet functions. It follows that both expressions are equivalent. However, we will see later that the matrix element of O 3 is ill-defined without an additional rapidity regulator, and this regulator would break the exchange symmetry if we would use either (2.31) or the corresponding other result. The symmetry can be restored by using the symmetric form At this stage, the factorization theorem for the b-quark induced contribution to the h → γγ decay amplitude can be written in the form (2.34) with γγ| O 3 |h given in (2.33). Subtleties arise from endpoint-divergent convolution integrals in all but the first term, which need to be properly identified and regularized. To see how the problem arises, we now present explicit expressions for the Wilson coefficients and operator matrix elements at NLO in α s .

Hard matching coefficients
We begin by describing the calculation of the Wilson coefficients H i , which contain the hard matching contributions obtained when matching QCD to SCET. The relevant Feynman diagrams up to one-loop order are shown in figure 6. Evaluating these diagrams in dimensional regularization with D = 4 − 2 spacetime dimensions and in the MS scheme, we obtain , JHEP04(2020)033 . (3.1) Here y b,0 is the bare b-quark Yukawa coupling, and α s,0 and α b,0 denote the bare QCD and electromagnetic couplings, respectively. These expressions are exact to all orders in . Note that the factor α b,0 included in H 1 arises because we have defined O 1 without including any electric charges in the operator. When the bare couplings are renormalized according to this ensures that the three Wilson coefficients depend on the dimensionless ratio We will see in the next section that the matrix element γγ| O 2 (z) |h is symmetric under the exchange z ↔ (1 − z), as is the case for H 2 (z). We can therefore rewrite We will also see that the matrix element is z-independent at leading order. It is then obvious that the integrals over z i in (2.34) diverge at the endpoints z i = 0 and z i = 1.
With the help of (3.3), the divergence is concentrated at z i = 0. Moreover, we will find later that the integrals over ± in (2.33) diverge when + → ∞ or − → ∞ at fixed + − . These divergences are not regularized by the dimensional regulator and require additional (rapidity) regulators. This will be discussed in more detail in section 5.

Operator matrix elements
In the next step, we need the matrix elements of the SCET operators in (2.10). We have calculated these matrix elements at NLO in α s in both SCET-1 and SCET-2. In SCET-1, the hard-collinear scaling can be enforced by taking the external photons off-shell, with virtualities (−k For our purposes, it is more important to have the matrix elements in SCET-2 obtained with physical (on-shell) photon states. To all orders in QCD perturbation theory, the matrix element of O 1 is given by (omitting the polarization vectors for simplicity) where m b,0 is the bare b-quark mass. The reason is simply that O 1 does not contain any fields with color charges, and hence there are no QCD corrections to the matrix element. The Feynman rule for the collinear photon fields A ⊥ n i ,µ gives e b ε ⊥ * µ (k i ). The relevant Feynman diagrams for the matrix element γγ| O 2 (z) |h are depicted in figure 7. This matrix element is symmetric under the exchange of z and (1 − z). We have not been able to calculate it in closed form for arbitrary values of . Expanding the NLO correction in powers of , we obtain (with 0 ≤ z ≤ 1) It is obvious from this result that at NLO in α s the matrix element contains terms that are singular at z → 0 or z → 1. The former terms are contained in K(z), while the latter ones are contained in K(1 − z). In order to compute the convolution of H 2 with O 2 , it is important to calculate these leading singularities exactly, without performing an expansion in . In the form (3.3), we only need to work out the behavior near z = 0. We obtain Note that, after renormalization of the bare couplings according to (3.2), the matrix element in (4.2) depends on the dimensionless ratio m 2 b,0 /µ 2 . Concerning the matrix element of the operator O 3 in (2.33), we now present our explicit expressions for the jet and soft functions at NLO in α s . The diagrams contributing to the jet function are shown in figure 8. We obtain This function has a cut along the positive p 2 axis starting at p 2 = 0 and extending to infinity. To the best of our knowledge, this expression (expanded in powers of ) was first obtained in the discussion of the exclusive B-meson decay process B − → γl −ν l in [58]. The calculation of the soft function is more complicated. The relevant Feynman graphs are shown in figure 9. Due to the multipole expansion applied to soft fields in interaction terms with hard-collinear fields, the soft momentum component + enters at the left lower vertex, while − exists at the right lower vertex, as indicated in the first graph. In a more complicated diagram such as the second graph, the assignment of momenta becomes non-trivial. For example, if the gluon emitted from the Wilson-line segment belongs to S  figure 9 has support for all values + − > 0, even though at leading order (first graph) the discontinuity arises only if where Here it suffices to quote the answer in term of series expansions in . In both expressions the dots refer to terms of O( ) and higher, which vanish for r → 0 orŵ → 0, respectively, and which are therefore not needed to reproduce the amplitude at NLO in α s . Note that the factor α b,0 in (4.6) ensures that after coupling renormalization the soft function depends on the dimensionless ratios m 2 b,0 /µ 2 , r andŵ. We conclude this section with an important remark. The term proportional to C 2 ( ) in (4.7) exhibits a stronger singularity at w = m 2 b,0 than the remaining terms. However, this term is absorbed entirely when the bare mass parameter m b,0 inside the leading-order term in S a (w) is renormalized and replaced by the b-quark pole mass.

Rapidity divergences and analytic regulators
The unregularized endpoint divergences of the convolution integrals involving the operators O 2 and O 3 can be tamed by introducing a rapidity regulator. Since the h → γγ process involves time-like kinematics, however, none of the regulators suggested in the literature can be used consistently in our case. Regularizing gluon emissions from Wilson lines [59] is not sufficient in our case, since the divergences arise already at zeroth order in α s . Using a non-analytic regulator such as |2 z | η , as proposed in [60], would destroy the analytic properties of the integral involving the soft function in (2.29), and hence the steps that led to (2.33) would no longer hold. Also, with this regulator one would loose certain minus signs in the arguments of logarithms, in such a way that the analytic properties of the decay amplitude are not correctly reproduced. In our case the rapidity logarithms in the sum of all terms in the factorization theorem (2.34) must conspire to give ln(−M 2 h /m 2 b −i0) with a non-zero imaginary part.

Introducing rapidity regulators
In our work, we will employ an analytic regulator on the convolution variables z 1,2 and ± in the factorization theorems (2.33) and (2.34). Figure 10 illustrates the physical origin of the endpoint divergences. In the n 1 -collinear region the variable z 1 = − /(n 1 · k 1 ) is formally of O(1), but it is integrated over the interval z 1 ∈ [0, 1], and the region near z 1 = 0 gives an unsuppressed contribution. The same applies in the n 2 -collinear region, where z 2 = + /(n 2 · k 2 ). In the soft region, on the other hand, the variables ± are formally of O(m b ), but they are integrated up to infinity. We thus see that the collinear and soft regions overlap, and this overlap must be avoided by introducing appropriate regulators. Rather than introducing the regulators in the SCET Feynman rules, in our problem it is simplest to impose them on the convolution variables themselves. This is similar in spirit to the approach of [61], where the regulator was imposed on the phase-space integrals for real emissions. For the soft contribution, we use the regulator under the integral in (2.29), which gives rise to another branch cut above the real axis in the complex − plane and thus does not invalidate the arguments that led to (2.33).
Here η is an infinitesimal analytic regulator and ν denotes the associated mass scale. For

JHEP04(2020)033
the collinear contributions, either + or − is of O(M h ), whereas the other variable is of O(m 2 b /M h ). We must therefore multipole expand the regulator term and use under the integral in (3.3). Using that M 2 h =n 1 · k 1n2 · k 2 , we see that the regulators indeed match up at the boundaries of the soft, n 1 -collinear and n 2 -collinear regions.
Introducing the regulators as just described, the bare factorization formula (2.34) takes the form We can now split up the double integral in the last contribution into regions where − > + and + > − . In the first region ( − > + ) rapidity divergences arise for − → ∞, and they cancel against rapidity divergences of the n 1 -collinear contribution in the second line, which arise for z 1 = − /(n 1 · k 1 ) → 0. Likewise, in the second region ( + > − ) rapidity divergences arise for + → ∞, and they cancel against rapidity divergences of the n 2collinear contribution in the second line. The two combinations of terms give the same result, so that we can write the final form of the bare factorization theorem as . It is important to perform the expansion in the analytic regulator η before the expansion in the dimensional regulator . Renormalizing the bare b-quark Yukawa coupling and quark mass in the MS scheme using 1 we write the answer for the three terms in the factorization formula (5.4) in the form with M 0 defined in (2.2). We find Explicit expressions for the NLO coefficients k i are listed in appendix C. Adding up the various contributions we recover the leading terms in the decay amplitude given in (2.1). Note that the 1/η poles cancel in the sum of all terms, in accordance with the physical picture discussed above. We now discuss the nature of this cancellation in more detail and derive the precise conditions under which it holds.

Cancellation of rapidity divergences
The cancellation of the 1/η poles in the sum of the contributions T 2 and T 3 in (5.7) is highly non-trivial, because the first term is a product of functions depending on the hard and collinear scales −M 2 h and m 2 b , whereas the latter term also depends on the hard-collinear scale M h m b . How is it conceivable that such a cancellation could persist in higher orders of perturbation theory? We will show now that this is possible because the hard matching coefficientH 2 (z) and the matrix element γγ| O 2 (z) |h in the second term in (5.4) can be factorized in the limit z → 0 into products which closely resemble the structure of the third term. 2 We introduce the notation [[f (z)]] to denote all those terms in a function f (z) that are of leading power (modulo logarithms) in the limit z → 0. For the bare functions we will consider, this means that we keep all terms of order z n with n ∈ Z. Using this notation, we find that the required refactorization conditions read

JHEP04(2020)033
and (5.9) where we have identified z = − /(n 1 · k 1 ) and + = w/ − in the last step. Our explicit expression given in the previous sections confirm that these relations indeed hold at NLO in α s . Concretely, relation (5.8) follows directly from (3.4) and (4.5), whereas relation (5.9) can be derived from (3.1), (4.5) and (4.6). But indeed these relations must hold to all orders of perturbation theory, and it would be desirable to derive them rigorously from SCET.
In order to show that the conditions (5.8) and (5.8) ensure that the rapidity divergences cancel to all orders of perturbation theory, we rewrite the factorization formula (5.4) in the short-hand form (5.10) We also adopt the default choice for the reference vectors n 1 and n 2 , so thatn i · k i = M h for i = 1, 2. Introducing the new variable w = + − = O(m 2 b ) and rearranging the order of integrations, we obtain for the third term (5.11) In the second term rapidity divergences arise from the region where z → 0. To evaluate them, we can use the refactorization conditions (5.8) and (5.9). This leads to (5.12) To proceed, we expand the bare jet function (4.5) in a power series in the bare coupling, with c 0 ( ) = 1. Using this expansion in (5.11), we obtain a double sum over terms proportional to c m ( ) c n ( ). We find that the terms with m = n cancel out in the sum of the two contributions, while those with m = n give (5.14) The terms inside the bracket contain the rapidity divergence and rapidity logarithms, whereas the genuine hard-collinear logarithms ln(−wM 2 h ) in the second line appear first JHEP04(2020)033 at O(α 2 s ). The corresponding expression derived from (5.12) reads (5.15) In the sum of the two contributions the 1/η poles cancel to all orders of perturbation theory, and we obtain At first sight, the fact that the rapidity divergences cancel between the second and third term in the factorization theorem (5.3) appears like a miracle, because these terms depend on different scales. The second term is a convolution of a hard function (depending on the scale −M 2 h ) and a collinear function (depending on the scale m 2 b ). The third term, instead, is a product of a hard function (depending on −M 2 h ), two hard-collinear functions (depending on the scale M h m b ) and a soft function (depending on the scale m 2 b ). What happens is that in (5.11) the two jet functions combine in such a way that their convolutions depend only on the products of scales (±M h w/ − )(∓M h − ) = −wM 2 h . This is parametrically equal to the product of the hard and collinear scales. This non-trivial fact makes the all-order cancellation possible.
At this point it is instructive to compare our results (5.4) and (5.16) with the corresponding expressions obtained in the simpler case of Sudakov resummation at leading power, as discussed in detail in [56,57]. An example would be the decay Z → qq, where Z is a color-neutral heavy vector boson. In this case there is no analogue of the operator O 1 , because the particles produced at the decay vertex of the Z boson are the same as the particles in the external state. More importantly, however, the are no non-trivial jet functions, because soft gluon exchanges between the (hard-)collinear quarks are described by soft Wilson lines, to all orders of perturbation theory. This leads to drastic simplifications in (5.16), as we can see by setting all expansion coefficients c n ( ) of the jet function to zero, except for c 0 = 1. This leads to (5.17) The two integrals over the soft function define two soft matrix elements living at the scale m 2 b,0 . The structure that multiplies the hard matching coefficient H 3 then contains two JHEP04(2020)033 terms: one associated with a single large rapidity logarithm and one without such an enhancement. In other words, the relevant soft matrix element is a linear function of the large rapidity logarithm, in agreement with the findings of [56,57].

Subtraction of endpoint divergences
We would now like to find a way to reproduce the answer (5.16) from an expression that is structurally similar to the original expression (5.3) but free of endpoint divergences. This is accomplished by the rearrangement leading power of the factorization formula. Note that, by construction, the second term is free of endpoint divergences. The regulator δ is introduced only in an intermediate step, such that both integrals are well defined. Also, as mentioned earlier, the region where + − → 0 in the last term does not give rise to endpoint divergences. The limit σ → −1 in this term is meant in the sense of analytic continuation. The integral over + must be evaluated assuming that the upper limit σM h > 0, and the result must then be analytically continued to σM h = −M h − i0. Note that the term in the third line of (6.1) is structurally similar to the resummation approach of [40][41][42][43][44], where the cutoffs on the integrals over ± were implemented by hand. However, we have provided precise definitions of the hard, jet and soft functions in this term, which are valid to all orders of perturbation theory. Also, we have derived the remaining terms in the factorization formula, which have not been discussed previously in the literature.
The contribution ∆H 1 is a "hard subtraction" defined by 2) The function S ∞ (w) is obtained by setting m b,0 → 0 in the expression for the soft function in (4.7), which leads to As a result, ∆H 1 depends on the hard scale (−M 2 h ) only. After a straightforward calcula-JHEP04(2020)033 tion, we obtain at NLO in α s .

(6.4)
Note that the third contribution in (6.1) contains some power-suppressed terms of order m 2 b /M 2 h arising from the upper limits of the integrations over ± . These must be dropped for consistency. In the definition of ∆H 1 , we have omitted such power corrections by using the limiting function S ∞ (w) instead of S(w) in (6.2). The projection onto the leading-power terms can be accomplished in a similar way be writing The last expression is free of power-suppressed terms.
We now obtain for the three terms in the rearranged factorization theorem (6.1), written as in (5.6), where in the last term L = L h − L m . Adding up the three contributions we correctly reproduce the QCD amplitude in (2.1). Note that the sum of the first two terms (T 1 + T 2 ) as well as the third term T 3 by itself are "almost" RG invariant. We believe that the remaining ζ 3 / terms are removed when the bare component functions in the factorization theorem are expressed in terms of renormalized functions. Some comments are in order concerning the structure of the various terms in (6.1). The subtraction of the leading singular terms in the second line, which removes the endpoint divergences of the z integral, generalizes a subtraction prescription proposed in [7,63], which works only for cases where the hard matching coefficient corresponding toH 2 (z) JHEP04(2020)033 approaches a constant plus power-suppressed terms as z → 0. In the present case the situation is more complicated, becauseH 2 (z) contains powers of ln z and hence is more singular. The form of the third term in (6.1) appears to violate strict factorization, because the soft momentum components + and − scale like m b in the soft region, but the integrals over these components extend into the hard region. In other word, this third term contains large (rapidity) logarithms beyond those arising from the RG evolution of the various component functions. This is a result of the collinear anomaly [64], which provides an extra source of large logarithms (so-called rapidity logarithms [59]) not related to conventional scale evolution. In simpler cases such as the q T resummation in Drell-Yan [65] or Higgs production [66], the collinear anomaly introduces a dependence of the product of two collinear beam functions on the hard scale, thereby violating strict factorization. Using simple differential equations, one can show that to all orders of perturbation theory this additional source of large logarithms takes the form of a linear function of ln(Q 2 /q 2 T ) in the logarithm of the cross section. In the present case, the anomalous dependence on the hard scale enters via the upper limits on the integrals in the third term in (6.1) and is thus of a more complicated nature.

Resummation of the leading double logarithms
In (6.1) we have obtained a factorization theorem for the b-quark induced h → γγ decay amplitude which is free of endpoint divergences. The collinear anomaly provides a source of large rapidity logarithms, which are contained in the third term. They arise from the fact that the soft momentum components ± are integrated up to values of O(M h ). Nevertheless, these logarithms would be under control if one knew how to obtain reliable results for the jet and soft functions J(p 2 ) and S(w) for arbitrary values of their arguments.
In its present form, the factorization theorem is written in terms of bare Wilson coefficients and operator matrix elements. The next step would be to replace these objects in terms of renormalized component functions. Note that except for the hard function H 3 all quantities must be renormalized "non-locally", by integrating them with Z-factors depending on two variables z and z or ± and ± . There is also a non-trivial operator mixing, such that the operators O 2 and O 3 require counterterms involving O 1 . While we have been able to calculate the relevant renormalization factors for the first two terms in the factorization theorem (including their mixing), the renormalization factors for the jet and soft functions are currently still unknown. We thus leave a detailed discussion of renormalization for future work.
The situation simplifies drastically if we are only interested in the leading doublelogarithmic corrections to the decay amplitude. They arise from the α s / 2 pole terms in the expressions of the bare functions contributing to the third term in the factorization formula (6.1), because only this term contains an extra enhancement by two rapidity logarithms. These leading poles can readily be identified using our results for the bare functions H 3 , J and S given in (3.1), (4.5), (4.6) and (4.7). Note that the leading poles (α s / 2 ) n in higher orders are all governed by the same coefficients. After renormalization of the bare JHEP04(2020)033 quark mass and Yukawa coupling, we obtain The leading double poles and the associated single poles multiplying the so-called cusp logarithms of the corresponding hard, hard-collinear and soft scales in (7.1) can be subtracted by means of local counterterms. The corresponding anomalous dimensions governing the scale evolution of the renormalized functions are of the form where the dots refer to terms without logarithmic enhancement, and is the cusp anomalous dimension [67]. Here n f = 5 is the number of light quark flavors.
To resum the leading double logarithms to all orders of perturbation theory, we solve the above evolution equations by ignoring all but the leading logarithmic terms. At the natural scales for each function, we use the lowest-order matching conditions, i.e.
In the approximation where the running of the strong coupling is ignored, and where one uses the one-loop expression for the cusp anomalous dimension, this leads to corresponding to a simple exponentiation of the correction given in (7.2). Note that the resummation shuts itself off when the soft momentum component − approaches M h or + approaches −M h . Thus, the exponential is well behaved even outside the region of soft momenta ± . Using the approximation (7.5) in the third term in the factorization formula (6.1), we can compute the infinite tower of leading double logarithms by performing the integral It is a straightforward task to include the running of the coupling in the solution of the evolution equations (7.3). Defining the Sudakov exponent S(µ 2 0 , µ 2 ) (not to be confused with the soft function) as [68] S(µ 2 0 , µ 2 ) = − 1 4 we obtain the solution

JHEP04(2020)033
Despite appearance, this expression is independent of the reference scale µ 2 . Once again, the resummation shuts itself off when the soft momenta ± approach their upper values ∓M h . Note that some terms in (7.8) involve the running coupling α s (µ 2 ) evaluated at timelike argument µ 2 = −|µ 2 | − i0. This quantity is well defined and can be computed using standard expressions for α s (µ 2 ). Indeed, using Sudakov exponents with time-like arguments is a well-known method to consistently resum large logarithms along with iπ terms to all orders of perturbation theory [45,46]. The resummation of the leading logarithms including the running of α s can be obtained by using expression (7.8) instead of (7.2) under the integral in (7.6). We defer a detailed numerical discussion of resummation effects to future work, where we will also resum large logarithms beyond the leading double-logarithmic approximation.

Conclusions and outlook
In this work we have started a detailed discussion of factorization at subleading power in SCET. At this order, factorization theorems for high-energy cross sections and decay amplitudes contain endpoint-divergent convolution integrals. Their presence indicates a violation of simple scale separation. These endpoint divergences cannot be removed using standard techniques of dimensional regularization and MS subtractions. They thus hint at an unexpected failure of conventional techniques used in renormalization theory. With the example of the b-quark induced radiative decay h → γγ of the Higgs boson we have worked out a rather complicated example, in which endpoint-divergent convolution integrals give rise to additional singularities that require both the dimensional regulator and a rapidity regulator η. We have derived a factorization theorem for the decay amplitude in terms of bare Wilson coefficients and operator matrix elements, and we have introduced a new rapidity regulator which preserves the analytic properties of the decay amplitude, taking into account the time-like kinematics of the h → γγ process. We have derived the conditions (5.8) and (5.9) under which the endpoint divergences caused by rapidity divergences cancel to all orders of perturbation theory. Moreover, we have shown that endpoint divergences that are regularized dimensionally can be avoided by rearranging the terms in the factorization theorem. While we have not yet derived the full set of renormalization factors and anomalous dimensions for the various objects entering the factorization theorem, we have succeeded in resumming the leading double-logarithmic corrections of order α n s ln 2n+2 (−M 2 h /m 2 b ) to all orders of perturbation theory. Our findings have important consequences for the general theory of electroweak Sudakov resummation in SCET. In the standard approach to this problem, as formulated in [56,57], one matches the full theory (the SM or an extension thereof) to SCET at a high scale Q v, where v ≈ 246 GeV (the equivalent of m b in our approach) denotes the scale of electroweak symmetry breaking, and Q (the equivalent of M h in our case) can for instance be the large center-of-mass energy of a collider process or the mass of a new heavy particle in a decay process. Evolving the effective theory from the high scale Q to the low scale v, one can resum large (double) logarithms of the scale ratio Q/v 1. In this evolution the particles of the SM are treated as massless. At a low scale of O(v), one then switches to an JHEP04(2020)033 effective theory with massive SM particles (and broken electroweak symmetry) by performing a second matching step. It was shown in [56,57] that the SCET matrix elements in this low-energy effective theory contain one power of the large rapidity logarithm ln(Q/v) due to the collinear anomaly, but that this does not invalidate the resummation procedure just described. Thus, the problem can be treated as a relatively simple two-scale process. We indeed find that this procedure is correct for processes described at leading power in the effective theory; see relation (5.17) and the discussion surrounding it. However, it needs to be generalized for processes such as h → γγ, which arise at subleading order in the expansion in scale ratios. These processes are characterized by three hierarchical mass scales: the hard scale Q, the low scale v, and the intermediate hard-collinear scale √ Q v. As can be seen from our result (5.16), the hard-collinear scale enters in an essential way, and the consistent resummation of all large logarithms requires that one takes this intermediate scale properly into account.
We close this paper with a speculation. Our final form of the factorization formula in (6.1) is explicit and well defined, however it is structurally rather different from the original form of the bare factorization theorem as given in (2.33) and (2.34). On the other hand, we have noted earlier in section 2.1 that the h → γγ matrix element of the operator O 2 (z) is proportional to the light-cone distribution amplitude of the photon, and that this quantity approaches the asymptotic form (2.18) when evolved to a high renormalization scale. As long as this matrix element vanishes at the endpoints z = 0 and z = 1, the integral in the second term in (2.34) is convergent and well defined. If an analogous statement holds for the product of the renormalized jet and soft functions in the third term, such that they tame the behavior of the integrand as ± → ∞, then also this integral would be well behaved. These observations motivate us to wonder whether it might be possible to derive a renormalized version of the factorization theorem, which is structurally equivalent to (2.34), i.e.  Choosing a high value of the renormalization scale, µ 2 = O(M 2 h ), would then cure the endpoint divergences by virtue of resummation. Also, with such a scale choice the first term on the right-hand side would not contain any large logarithms and could be calculated in fixed-order perturbation theory. While it is tantalizing to speculate about this possibility, we have so far not succeeded in deriving a formula such as (8.1) from first principles. The bottleneck is that in the regularized theory (before resummation) the various convolution integrals give rise to endpoint divergences, which are regularized either dimensionally (giving rise to 1/ poles) or by means of the rapidity regulator (giving rise to 1/η poles). In an MS-like subtraction scheme, one would remove these singularities by applying renormalization factors to the operators and hard functions. In the presence of endpoint singularities, this is however not sufficient to remove all poles. One would thus need to generalize the JHEP04(2020)033 MS scheme and define Z factors for the convolution integrals themselves. The systematics of such an approach will still have to be developed. JHEP04(2020)033 When these expressions are combined with the hard matching coefficients given in (3.1), one finds that all endpoint divergences are regularized dimensionally. After renormalizing the bare quark mass and Yukawa coupling using (5.5), we obtain the finite, off-shell h → γγ amplitude

B Exact analytic expression for the soft function
In (4.7) we have presented some terms in the expression for the soft function using an expansion in powers of to the order required for our analysis. We have, however, also obtained closed analytical expressions for the component functions S a (w) and S b (w) in terms of generalized hypergeometric functions. We find

C NLO coefficients in the bare decay amplitude
Here we list our explicit findings for the coefficients k i accounting for the NLO corrections in (5.7). We obtain JHEP04(2020)033 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.