Renormalization and Scale Evolution of the Soft-Quark Soft Function

Soft functions defined in terms of matrix elements of soft fields dressed by Wilson lines are central components of factorization theorems for cross sections and decay rates in collider and heavy-quark physics. While in many cases the relevant soft functions are defined in terms of gluon operators, at subleading order in power counting soft functions containing quark fields appear. We present a detailed discussion of the properties of the soft-quark soft function consisting of a quark propagator dressed by two finite-length Wilson lines connecting at one point. This function enters in the factorization theorem for the Higgs-boson decay amplitude of the $h\to\gamma\gamma$ process mediated by light-quark loops. We perform the renormalization of this soft function at one-loop order, derive its two-loop anomalous dimension and discuss solutions to its renormalization-group evolution equation in momentum space, in Laplace space and in the"diagonal space", where the evolution is strictly multiplicative.


Introduction
Soft-collinear effective theory (SCET) offers a convenient framework for analyzing the factorization properties of cross sections and scattering amplitudes sensitive to different, hierarchical scales [1][2][3][4]. The corresponding factorization theorems contain hard functions, jet functions and soft functions, which receive contributions from different momentum regions in Feynman diagrams. The hard functions correspond to Wilson coefficients obtained when the full theory is matched onto SCET, while the jet and soft functions are defined in terms of matrix elements in the low-energy effective theory. Soft functions -matrix elements of non-local products of soft fields dressed by Wilson lines -play a particularly important role in the factorization theorems, because they often capture the physics at the longest relevant distance scales in a given process. In some cases the soft functions are non-perturbative objects, whereas in others they can be calculated using perturbation theory.
Recently, there has been a growing interest in understanding factorization at subleading power in scale ratios. In this case a large number of hard, jet and soft functions appear. In particular, while at leading power soft emissions are eikonal and can be described by soft Wilson lines, at subleading power the emission of soft fermions and power-suppressed emissions of soft gauge bosons need to be taken into account. Particularly interesting is the case of soft quark emission, which is absent at leading power. At subleading order in the SCET expansion there is a unique interaction that couples a soft quark to collinear quarks and gauge fields [3].
As a concrete example, we have considered in [5] the case of the h → γγ decay amplitude induced by b-quark loops. The relevant soft function is derived from the vacuum matrix element of the soft quark propagator dressed by two finite-length soft Wilson lines, i.e.

−
(4π) 1− N c e γ E µ 2 0| T Tr Sn(0, r 1n ) q s (r 1n )q s (r 2 n) S n (r 2 n, 0) |0 , where the trace is over color (but not spinor) indices, and = (4 − d)/2 is the dimensional regulator. The prefactor is chosen for later convenience. We have introduced two light-like reference vectors n µ andn µ (with n ·n = 2) aligned with the directions of the two photons in the h → γγ process. The soft quark fields are displaced from the Higgs vertex at position z = 0 by light-like distances along the n andn light cones. The Wilson lines S n and Sn connect the soft quarks with the Higgs vertex and ensure that the matrix element is gauge invariant. A different soft-quark soft function, which is relevant for inclusive cross sections, was introduced in [6] and has been studied further in [7,8].
In the context of the SCET, soft quarks couple to collinear fields via subleading interactions in the SCET Lagrangian [3]. In our case the soft matrix element in (1) is sandwiched between projection operators P n = / n/ n 4 from the right, where it connects to n-collinear fields, and Pn = / n/ n 4 from the left, where it connects withn-collinear fields. In the first step, we define a soft function S( + − ) via the Fourier transform of the above matrix element with respect to the coordinates r 1 and r 2 , N c e γ E dr 1 e ir 1 − dr 2 e −ir 2 + × µ 2 0| T TrPn Sn(0, r 1n ) q s (r 1n )q s (r 2 n) S n (r 2 n, 0) P n |0 . (2) Reparameterization invariance [9] ensures that the soft function S only depends on the product w = + − of the Fourier variables ± . At lowest order in perturbation theory (but not beyond) these variables can be identified with the light-cone components n · andn · of the soft momentum µ flowing through the soft quark propagator. In the next step we define a new soft function S(w) in terms of the discontinuity of the function S( + − ), i.e.
Our definition of the soft function differs from the one in [5] by a factor (−N c α b /π), where α b = α/9 is the electromagnetic coupling of the b quark.
In Section 2 we collect the expression for the bare soft function obtained at one-loop order. A heuristic derivation of the non-local renormalization factor, which subtracts the divergences in the bare function, is presented in Section 3 along with the expression for the renormalized soft function. In Section 4 we discuss the renormalization-group (RG) evolution equation satisfied by the soft function and present a conjecture for its anomalous dimension at twoloop order. We also construct an exact solution to the RG equation at next-to-leading order (NLO) in RG-improved perturbation theory. The asymptotic behavior of the RG-improved soft function for large values w m 2 b is studied in Section 5, where we introduce the concept of dynamical scale setting. Section 6 contains a discussion of the properties of the Laplace transform of the soft function, which satisfies a novel form of RG evolution equation. In Section 7 we study the soft function in the "diagonal space", in which its RG evolution is strictly local in w. To this end we generalize the concept of the "dual space" introduced in [10,11] to higher orders of perturbation theory. Finally, in Section 8 we analyze in detail the double convolution integral T 3 arising in the factorization theorem for the h → γγ decay amplitude [5] and show that, after a suitable rapidity regularization in the diagonal space, this quantity is RG invariant and free of endpoint divergences. Section 9 contains a summary of our main results and some conclusions. Several technical details of our analysis are presented in three appendices.

One-loop expression for the bare soft function
At one-loop order in perturbation theory the bare soft function has been calculated to all orders in the dimensional regulator = (4 − d)/2 in terms of hypergeometric functions [5]. The relevant Feynman graphs are shown in Figure 1. Due to the multipole expansion applied to soft fields in interaction terms with collinear fields, the soft momentum component + enters at the left lower vertex, while − goes out at the right lower vertex, as indicated in the first graph. In more complicated diagrams, such as the second graph, the assignment of momenta becomes non-trivial due to the presence of the Wilson lines (see [5] for a detailed discussion). As a result, one finds that the soft function S(w) defined via the discontinuity of the diagrams shown in Figure 1 has support for all values w > 0, even though at leading order (first graph) the discontinuity arises only if w > m 2 b . When expanded about = 0 the result reads Figure 1: Representative Feynman diagrams contributing to the soft function S( + − ) up to O(α s ) (taken from [5]). We use a double-line notation to represent the finite-length soft Wilson lines. The three points mark the vertices connecting to the incoming Higgs boson (top) and the outgoing photons (bottom). The soft function S(w) in (3) is given by the discontinuity of S(w). where (5) Here α s,0 is the bare QCD coupling and m b,0 denotes the bare mass of the b-quark. We have pulled out a factor e γ E / (4π) from the bare coupling, as appropriate in the MS scheme. Moreover, we have defined the dimensionless ratioŵ 0 = w/m 2 b,0 . To express the result in terms of physical parameters we renormalize the coupling according to α s,0 = µ 2 Z α α s (µ), where Z α = 1 + O(α s ). From now on α s ≡ α s (µ) always denotes the renormalized coupling. The most convenient scheme for the renormalization of the b-quark mass is the pole scheme, in which the mass is defined by the position of the pole in the renormalized quark propagator [12]. We denote the pole mass by m b . At one-loop order one finds that Renormalizing the quark mass in the pole scheme removes entirely the loop correction to the lowest-order term shown in the first line of the expression for S a (w) in (5). After the bare parameters have been renormalized, we obtain where where nowŵ = w/m 2 b . The remaining 1/ n poles must be removed by renormalizing the soft function itself. A feature that is at first sight surprising is the appearance of the 1/ pole in the function S (0) b , which vanishes at lowest order. It is obvious that this pole can only be removed by means of a renormalization factor Z S (w, w ; µ) that is non-local in the space of w values. In the following section we will present a conjecture for this renormalization factor and show that it indeed removes all remaining singularities.

Renormalization of the soft function
The soft-quark soft function appears in the analysis of the factorization theorem for the amplitude for the radiative Higgs-boson decay h → γγ induced by loops containing light b-quarks [5]. While this is not the dominant contribution to the decay amplitude (which instead comes from loops containing top-quarks or W bosons), this particular contribution has an interesting structure, since it receives large double-logarithmic corrections of order αα n s L 2n+2 with n ∈ N 0 , where L = ln(M 2 h /m 2 b ) − iπ [13][14][15][16][17]. M h m b denotes the mass of the Higgs boson. In a recent paper [5], two of us have derived a factorization theorem in SCET, with which these large logarithms can eventually be resummed to all orders in perturbation theory. The factorization theorem contains three terms: a term T 1 involving a hard function H 1 multiplied with the h → γγ matrix element of a local SCET operator containing a Higgs field and two photon fields, a second term T 2 involving the convolution of a hard function H 2 with a nonlocal SCET operator containing a Higgs field, a photon field and two collinear quark fields, and a third term T 3 containing a hard function H 3 multiplied with a double convolution of the soft-quark soft function with two jet functions. In the first two terms the SCET matrix elements live at the scale m b , however the matrix element in third term depends in addition on We have presented arguments in [5] suggesting that the third term should by itself be RG invariant. More accurately, we found that when this term is evaluated as a convolution over bare functions, with cutoffs implemented to remove endpoint divergences, almost all 1/ n poles cancel out. A single remaining pole term proportional to ζ 3 / can be attributed to the fact that the implementation of cutoffs on the bare functions is not strictly compatible with RG invariance (see also the discussion in Section 8).
We will now rely on these arguments to derive the renormalization condition for the soft function. In terms of bare functions, the third term of the factorization theorem for the bquark induced h → γγ decay amplitude takes the form (for now we omit the cutoffs required to defined this expression properly) Note that one of the jet functions J (0) (p 2 ) is evaluated at time-like momentum (p 2 > 0) and the other one at space-like momentum (p 2 < 0). The renormalization of the hard function is well understood, with the corresponding anomalous dimension being known to three-loop order [18]. The renormalization of the jet function J(p 2 ) was derived at one-loop order a long time ago [19] and has recently been extended to two loops [20]. One finds that the time-like (space-like) jet functions at different p 2 values mix with each other. The renormalization conditions can be written in the form where at one-loop order Here and below −M 2 h ≡ −M 2 h − i0 and −p 2 ≡ −p 2 − i0. We have introduced the symmetric distribution [19] The plus prescription is defined such that, when Γ(ω, ω ) is integrated with a function f (ω ), one must replace f (ω ) → f (ω ) − f (ω) under the integral. Using the inverse of the relations (10) we can express the first three functions in the double convolution for T 3 in (9) in terms of their renormalized counterparts. Requiring that the combined expression for T 3 is scale independent we then obtain where the renormalized soft function is defined as At one-loop order the inverse renormalization factors are given by the same expressions as in (11), but with the signs in front of α s reversed. The above renormalization condition can be simplified by carefully evaluating the convolution of the three Z factors. After some algebra we find the simple form where Note that the same plus distribution Γ(w, w ) as in (12) appears, and that the logarithmic terms contained in Z 33 and Z −1 J have conspired to generate a logarithm of the ratio w/µ 2 . This fact has already been anticipated in [5].
It is not at all obvious that this definition of the renormalized soft function ensures that all 1/ n pole terms in (8) are removed. Our "derivation" of the renormalization factor Z S is only a conjecture, since the convolution of the four bare functions in (9) contains endpoint divergences for ± → ∞ and is thus ill-defined. Nevertheless, applying the renormalization condition (15) to the bare soft function in (7) we find that the renormalization factor Z S (w, w ; µ) indeed removes all 1/ n pole terms. Note that (contrary to the one-loop renormalization of the jet function [19]) the plus distributions have a non-trivial effect, because the lowest-order soft function is not constant, see (8). For the renormalized soft function at one-loop order we obtain with We have defined L w = ln(w/µ 2 ) and L m = ln(m 2 b /µ 2 ). Figure 2 shows the renormalized soft function S(w, µ) in units of m b as a function of the dimensionless ratioŵ = w/m 2 b . We use m b = 4.8 GeV for the b-quark pole mass and choose µ = m b for the renormalization scale. Note that both curves are discontinuous atŵ = 1.

Renormalization-group evolution
The dependence of the renormalized soft function on the scale µ can be controlled by means of the RG evolution equation where the anomalous dimension is obtained from the coefficient of the single 1/ pole in Z S via We have checked explicitly that our expression for the renormalized soft function in (17) satisfies the evolution equation (19) at O(α s ).

Two-loop anomalous dimension
We can derive additional information by using the known RG equations obeyed by the hard and jet functions. To all orders in perturbation theory, the evolution equation for the hard function takes the form [18,21] where Γ cusp is the light-like cusp anomalous dimension in the fundamental representation of SU(N c ) [22], and γ q is the anomalous dimension of the quark field in light-cone gauge [21]. The cusp anomalous dimension is known to four-loop order while γ q is known to three loops. The anomalous dimension for the jet function has only recently been calculated at two-loop order [20]. The renormalized jet functions in (13) obey the RG evolution equation Compared with (11) we have introduced a dimensionless integration variable x, such that this relation holds for both time-like and space-like values of p 2 . The anomalous dimension is given by where The local terms (with x = 1) can to all orders by expressed in terms of the cusp anomalous dimension and an anomalous dimension γ (α s ). Since the plus distribution contained in Γ(1, x) is linked with the logarithmic term, it is also multiplied by Γ cusp . However, starting at twoloop order additional non-local terms arise, whose explicit form was obtained in [20] by using the RG invariance of the B − → γ −ν decay rate along with the calculation of the two-loop anomalous dimension of the B-meson light-cone distribution amplitude (LCDA) performed in [23]. With the help of the above expressions we can write the anomalous dimension for the soft function defined in (19) in the more general form where γ s (α s ) = 2γ q (α s ) + 2γ (α s ) .
Explicit expressions for the various anomalous dimensions are given in Appendix A. While the one-loop expression for the soft anomalous dimension shown in (20) is derived from the one-loop renormalization factor Z S , whose explicit form is checked by the fact that it properly removes all 1/ n pole terms in the bare soft function at one-loop order, we stress again that the two-loop form of the anomalous dimension shown above should be considered as a conjecture.

Exact solution to the RG equation
The RG equation (19) and the associated anomalous dimension γ S in (25) are closely related to the corresponding equations for the leading-twist LCDA φ B + (ω) of the B meson [24][25][26]. This correspondence is discussed in more detail in Appendix B. Structurally, the anomalous dimension γ S differs from the anomalous dimension for the LCDA by the argument of the logarithm (w/µ 2 versus ω/µ, because w has mass dimension 2 while ω has mass dimension 1) and a factor 2 in front of the non-local terms. In solving the evolution equation for the soft function we will use and extend some of the strategies developed in the literature on the B-meson LCDA. We will see, however, that the additional factor 2 in front of the non-local terms in the anomalous dimension give rise to major complications.
A closed analytic solution to the evolution equation (19) can be obtained based on the observation that, on dimensional grounds, any pure power w a provides an eigenfunction of the evolution kernel. To see this, note that the integrals [25] where H(a) = ψ(1 + a) + γ E is the harmonic-number function, and introduced in [20]. It follows that the ansatz with α s (µ α ) ≡ α provides a solution to the RG equation (19) with the initial condition (w/µ 2 s ) η at some matching scale µ = µ s , at which the soft function is free of large logarithms. Here β(α s ) = dα s /d ln µ is the QCD β-function, and we have introduced the RG functions (the first of which should not be confused with the soft function) They are the solutions to the equations with the boundary conditions S(µ s , µ s ) = 0 and a Γ (µ s , µ s ) = 0. The function a γs (µ s , µ) is defined analogously to a Γ (µ s , µ). Note that both S(µ s , µ) and a Γ (µ s , µ) take negative values if µ > µ s , because the cusp anomalous dimension is a positive quantity. Explicit expressions for these objects obtained at NLO in RG-improved perturbation theory are given in Appendix A.
Changing variables in the integral in the exponent in the first line of (30) from α to = a Γ (µ s , µ α ) allows us to perform this integral in closed form. Moreover, in the exponent of the term in the second line we can use that β(α s ) = −β 0 α 2 s /(2π) + . . . at leading order. We can thus rewrite (30) in the more compact form where we have defined the evolution function which satisfies U S (w; µ s , µ s ) = 1.
In order to apply this solution, we need to recast the initial condition for the soft function at the matching scale µ s in such a way that it is written as a superposition of pure powers (w/µ 2 s ) η . To this end, we express S(w, µ s ) in terms of its Laplace transform with respect to ln(w/m 2 b ), which in general is defined as In fixed-order perturbation theory the integral converges for 0 < η < 1. We thus write with 0 < c < 1. Since the right-hand side of (36) exhibits a power-law dependence on w, we can use the solution (33) to express the soft function at a different scale µ > µ s as In this expression the contour must be chosen such that 0 < c < 1 + min 0, a Γ (µ s , µ) , which in turn requires a Γ (µ s , µ) > −1. 1 Relation (37) provides an exact solution of the evolution equation (19) in the approximation where one ignores the yet unknown three-loop non-local terms in the soft anomalous dimension in (25), which would enter in the indicated higher-order corrections in the exponent of the last factor. However, this solution requires the Laplace transformS(η, µ s ) of the soft function at the matching scale µ s . Clearly, it would be more convenient to have a solution of the evolution equation that relates S(w, µ) in a more direct way with the initial condition S(w, µ s ). To obtain it, we use the fact that the exponent in the second line of (37) is of higher order in α s and hence can be expanded out. Eliminating the Laplace transformS(η, µ s ) by means of (35) we then obtain where we have used the definition of the function H(a) in (29). The integrand of the integral over η has single and double poles located at η = n + a Γ (µ s , µ) and η = −n for n ∈ N. As long as a Γ (µ s , µ) > −1, these two series of poles lie on different sides of the integration contour if we choose 0 < c < 1 + a Γ (µ s , µ). The integral can then be expressed in terms of Meijer G-functions (see e.g. [27,28] for a detailed discussion). This yields where a = a Γ (µ s , µ). This form of the solution is valid as long as one considers an upward scale evolution (µ > µ s ), for which a < 0. The Meijer G-function appearing above becomes singular when the last argument approaches 1, This is an integrable singularity only as long as a < 0. It is possible to relate the Meijer Gfunction to an expression involving a hypergeometric function and its derivatives by performing the contour integral in (38) using the theorem of residues. This is discussed in more detail in Appendix C.
There are two further simplifications that can be made, because the second term inside the brackets in relation (39) is of O(α s ). For this term it is thus sufficient to use the leading-order expressions for the soft function and for the RG function where β 0 = 23 3 is the one-loop coefficient of the β-function evaluated with n f = 5 light quark flavors (see Appendix A). We then find the final solution with r = α s (µ)/α s (µ s ) and a = a Γ (µ s , µ). The Meijer G-function enjoys the properties which can be used to restrict the last argument to values in the interval [0, 1]. This is particularly useful for a numerical evaluation. Relation (42) provides the desired result for the exact evolution of the soft function from the matching scale µ s to a different scale µ at NLO in RG-improved perturbation theory. In evaluating this result consistently one should use the one-loop expression (17) for the soft function S(w, µ s ) at the matching scale, the two-loop expression for the anomalous dimension γ s and the three-loop expression for the cusp anomalous dimension (see Appendix A). In the literature on double logarithmic resummations this approximation is often referred to as the next-to-next-to-leading logarithmic (NNLL) approximation. It correctly resums logarithms of the form α n s L k with 2n − 3 ≤ k ≤ 2n [29]. The above solution simplifies considerably if one works at leading order in RG-improved perturbation theory, corresponding to the NLL approximation, which resums the logarithms with 2n − 1 ≤ k ≤ 2n. One then uses the two-loop approximation for the cusp anomalous dimension, the one-loop approximation for γ s and the tree-level matching condition for the soft function at the scale µ s . In this approximation the integral over the Meijer G-function in (42) can be performed analytically, and we find In Figure 3 we show the renormalized soft function S(w, µ) at different values of the renormalization scale µ and fixed matching scale µ s = m b = 4.8 GeV. These solutions are derived by evolving the fixed-order expression at the matching scale µ s = m b obtained from (17), and shown in Figure 2, to the scale µ using the leading-order and NLO solutions to the RG equation presented in (44) and (42), respectively. Note that the discontinuous behavior of the soft function at w = m 2 b present in fixed-order perturbation theory (see Figure 2) has been smoothed out after RG evolution to a higher scale. However, for relatively low values of µ a prominent feature near w = m 2 b remains.

Asymptotic behavior and dynamical scale setting
In the solutions discussed in the previous section it is important that the matching scale µ s is chosen such that the initial condition S(w, µ s ) for the soft function is free of large logarithms. At one-loop order and beyond, the explicit expression for this initial condition involves the logarithms L w = ln(w/µ 2 ) and L m = ln(m 2 b /µ 2 ), see (18). The RG-improved solutions for the soft function are thus well defined as long as µ s ∼ m b ∼ √ w are all of the same order. Indeed, in Figure 3 we have set µ s = m b and varied w up to a maximum value such that √ w < √ 3 m b . Clearly, it would be desirable to be able to treat the two soft scales m b and w as independent and extend the solution for the soft function up to values w m 2 b . We now show how this can be accomplished.
We have seen in (16) and (25) that the renormalization factor and the anomalous dimension for the soft function do not contain any reference to the b-quark mass. Apart from the prefactor of m b in the normalization of the soft function, nothing prevents us from studying this function and its RG evolution in the limit m b → 0 or, more properly, w/m 2 b → ∞. Denoting the function in this limit by S ∞ (w, µ), we find from (17) One subtle point is that in this limit one should use the running b-quark mass m b (µ) in the prefactor, as shown in the second line. The reason is that the pole mass contains a logarithm L m in its definition, and hence the limit m b → 0 would be singular. Once this is done, the higher-order corrections to the soft function in the limit of large w m 2 b only contain powers of the logarithm L w and constants. This fact greatly simplifies the solution of its RG evolution equation. In a first step, we define a new function S ∞ (L w , µ) via From (33) it follows that the general solution to the evolution equation for S ∞ (w, µ) can be written in the form [20] S At NLO in RG-improved perturbation theory one can set η = 0 in the argument of the function H, because the tree-level term in S ∞ (L w , µ s ) is a constant. In this expression, which is indeed much simpler than the exact solution (42), one needs the explicit form of the function H(a) given in (29). The NLO solutions (42)  only. Note that for a fixed value of the renormalization scale µ there is always a range of w values for which µ s > µ. It is important in this context that the asymptotic solution is not restricted to the case where µ > µ s , unlike the solution in (42). Figure 4 shows different results for the soft function S(w, µ) obtained at NLO in RGimproved perturbation theory. The renormalization scale is fixed to µ = 60 GeV. We still normalize our results to the b-quark pole mass, such that the plots can be compared with those in Figure 3 . To obtain it, we rewrite (17) by pulling out the running mass m b (µ) rather than the pole mass. This has the effect of replacing the term (3L m + 8) in the expression for the function S a (w, µ) in (18) by 12. The fact that there is no such compensating effect for the function S b (w, µ), which starts at O(α s ), explains the small difference between the solid and dashed lines in the region of small w. We then use this modified form for the matching condition S(w, µ s ). The solid curve provides the optimal RG-improved solution for the soft function valid for small and large values of w. Note, however, that in practice we cannot evaluate this solution for arbitrarily large w values, because this would violate the condition µ s < µ, which is a prerequisite for the analytic solution shown in (42) to hold. We observe that the solution with the fixed choice µ s = m b of the matching scale breaks down for w > 3m 2 b , because the "large logarithms" ln(w/µ 2 s ) arising in this region are not resummed in this case. The asymptotic solution (48) provides an accurate description for w > 10m

Soft function in Laplace space
A much simpler solution to the RG evolution of the soft function is found in Laplace space. We define the Laplace transformS(η, µ) of the soft function at the scale µ as shown in (35). It then follows from (37) that Contrary to (39) no integral over the soft function is required. Instead, under a scale transformation the argument of the Laplace transform is shifted by an amount a Γ (µ s , µ). Recall that the Laplace transformS(η, µ s ) of the soft function at the matching scale µ s exists for 0 < η < 1, as is evident from the explicit formula shown in (51) below. It then follows from the above result that the Laplace transform at a higher scale µ > µ s exists for −a Γ (µ s , µ) = |a Γ (µ s , µ)| < η < 1.
It is not difficult to check that (49) is the solution of the partial differential equation where the functions F(a) and H(a) have been defined in (27) and (29), respectively. This generalized RG evolution equation forS(η, µ) can be derived from relation (36), with µ s replaced by µ, and the RG equation for the soft function following from (19) and (25). The Laplace transformS(η, µ s ) at the matching scale µ s can be derived from the NLO expression for the soft function given in (17). At one-loop order we obtaiñ As before H(a) is the harmonic-number function and L m = ln(m 2 b /µ 2 s ). Given this result, it is straightforward to work out the NLO solution forS(η, µ) at a different scale from (49).
The solution for the Laplace transform of the soft function derived here is more than just an academic result. In fact, it is well known that many factorization theorems involving complicated convolutions of jet and soft functions take on a particularly simple form in Laplace space (see e.g. [30,31]). The reason is that after RG improvement the jet functions J(p 2 ) can often be written as a derivative operator acting on a factor p 2 /µ 2 j a with some exponent a. Thus the convolution of one or more such jet functions with a soft function evaluates naturally to the Laplace transform of the soft function.

Soft function in the diagonal space
It has been shown in [10,11] that one can bring RG evolution equations of the type (19)in the approximation where the non-local two-loop contributions in the anomalous dimension γ S shown in the second line of (25) are neglected -into a simpler form using a suitably chosen integral transformation. The so-called "dual" soft function s(w, µ) obtained after this transformation obeys an evolution equation that is local in the variable w and hence much easier to solve. In the context of our problem, this equation would take the form We discuss the construction of the dual soft function in the above-stated approximation in Appendix D. It requires a straightforward extension of the technology developed in [10]. 2 When the non-local two-loop contributions neglected above are put back in place, we find that the dual soft function obeys the RG equation This follows from the analogy with the case of the B-meson LCDA discussed in Appendix B, for which the two-loop anomalous dimension in the dual space has been derived in [23]. Equation (54) is simpler than (25), because the term proportional to the plus distribution Γ(w, w ) has been eliminated. However, starting at two-loop order the dual anomalous dimension contains non-local terms with w = w , thus upsetting the original motivation for the construction of the dual space. The complexity of the original evolution equation is therefore reduced only marginally. In this sense the approach of [10,11] leads to a hybrid representation, in which the one-loop anomalous-dimension kernel is diagonalized but higher-order corrections to it are not. We will now generalize the approach of [10,11] and construct what we refer to as the "diagonal space", in which the RG equation for the soft function retains the simple local form (52) to all orders of perturbation theory.

Construction of the diagonal space
The key observation is based on the structure of the solution to the RG equation for the Laplace transform of the soft function given in (49). Shifting the Laplace variable η by −a Γ (µ s , µ) we can rearrange this solution in the form where we have used that a Γ (µ s , µ) + a Γ (µ, µ α ) = a Γ (µ s , µ α ). The auxiliary scale parameter ρ is arbitrary and only serves to split up the integral into two terms. It follows that the function has a particularly simple behavior under RG evolution, namely g w, η − a Γ (µ s , µ), µ, ρ = U S (w; µ, µ s ) g(w, η, µ s , ρ) .
Under scale evolution from µ s to µ the second argument of g gets shifted and the entire function is rescaled by the η-independent factor U S . If we integrate both sides of this equation over a suitably chosen contour in the complex η-plane, which runs parallel to the imaginary axis with Re(η) = c in the interval 0 < c < 1 + min 0, a Γ (µ s , ρ) , then the shift becomes unobservable. Indeed, defining the soft function in the diagonal space via the inverse Laplace transform where −a Γ (µ s , µ) < c < 1 + min 0, a Γ (µ, ρ) , we obtain the simple result s(w, µ, ρ) = U S (w; µ, µ s ) s(w, µ s , ρ) .
This function obeys the local RG equation (52) to all orders in perturbation theory. 3 In our discussion we have only included the non-local two-loop terms in the anomalous dimension (25), but the method described here can readily be extended to higher orders. When terms beyond the two-loop order are taken into account, one just needs to add the corresponding higher-order corrections in the exponential in (56). The soft function in the diagonal space depends on the auxiliary factorization scale ρ in addition to the renormalization scale µ. While its evolution in µ is local, as shown in (52), the evolution in ρ is more complicated. From (56) it follows that In all physical quantities the dependence on the scale ρ drops out. Note that logarithms of the scale ratio ρ/µ are already resummed in this expression, so there is no need to choose these two scales to be of the same order.

Transformation between momentum space and diagonal space
Without loss of generality we write the linear relations connecting the original soft function S(w, µ) with the soft function s(w, µ, ρ) in the diagonal space in the general form where the two transfer functions F S and F inv S explicitly depend on the renormalization scale µ and on the auxiliary scale ρ. Combining (58), (56) and (35) we find that (62) Likewise, from (36) and (56) we obtain where in an intermediate step we have defined the Laplace transform of the function s(w, µ, ρ) in analogy with (35). Changing the sign of the integration variable η brings the factor in front of the exponential to the same form as in (62).
We now use the fact that the arguments of the exponentials in both expressions can be expanded in a perturbative series in α s (µ), with coefficients that depend on the ratio α s (ρ)/α s (µ). Expanding the transfer functions at NLO in the form and similarly for F inv S (x, µ, ρ), we obtain To arrive at this form we have employed the leading-order approximation for a Γ (µ, µ α ) shown in (41) and the definition (29) relating H(a) to the function h(x) in (24). The solutions involve another example of a Meijer G-function. An alternative expression involving a hypergeometric function and its derivative can be derived by performing the integral in (62) using the theorem of residues. This is discussed further in Appendix C. While the transfer functions for the soft function are given by rather complicated expressions, they nevertheless enjoy some nice properties. Combining the two relations in (61) one can derive the orthonormality condition (with a, b > 0) To prove this identity directly one uses the useful relation which holds for purely imaginary (η − η) = it. We also find that powers of x are very simply transformed into the diagonal space, because Expanding this relation in powers of b yields the transformation rules for positive integer powers of ln x.

Convolution T 3 in the diagonal space
Remarkably, the factorization formula (13) retains its form in the diagonal space, i.e.
Here j(p 2 , µ, ρ) denotes the jet function in the diagonal space. To derive this relation we need the inverse transfer function for the jet function in the diagonal space, which we define as In analogy with (62) we find (drawing on expressions derived in [20]) (71) Expanding this function in a perturbative series as in (64) we obtain To establish that relation (69) holds we need to show that Using the representation for F inv J (x) given in (71) along with the identity (67) we find that the right-hand side of (73) indeed leads to the expression for the soft transfer function shown in (62). In Figure 5 we show the behavior of the leading-order transfer functions F S,0 (x) for the soft function and F J,0 (x) = J 1 (2 √ x) for the jet function. Both functions vanish in the limit x → 0 and oscillate for values x > 1. The transfer function for the jet function oscillates more rapidly, but apart from this the two functions exhibit a rather similar behavior. It is important that the jet functions and the soft function in (69) are evaluated at the same auxiliary scale ρ. Only then the dependence on ρ cancels out. In analogy with (60), we find that the jet function in the diagonal space satisfies the different equation It is then straightforward to check that the double convolution on the right-hand side of (69) is independent of ρ.

Leading-order soft function in the diagonal space
Given the transformations in (61) we can calculate the soft function s(w, µ, ρ) in fixed-order perturbation theory from the expression for the original soft function shown in (17). In this way we obtain s(w, µ s , ρ) in the form of a numerical integral, where µ s denotes the matching scale, at which a fixed-order expansion is well behaved. We can then evolve the obtained result to a different renormalization scale µ = µ s using the explicit solution (59) of the local RG equation in the diagonal space. At leading order in RG-improved perturbation theory it is possible to derive a simple formula for s(w, µ, ρ), which in this approximation is independent of the auxiliary scale ρ. Using that at the matching scale one can express the RG-evolved soft function in the diagonal space in terms of another Meijer G-function. We find Note that the last form is structurally similar to the momentum-space result shown in (44). In Figure 6 we show the leading-order solution for the soft function s LO (w, µ) at different values of the renormalization scale µ and fixed matching scale µ s = m b = 4.8 GeV. The behavior of the results in the region w > m 2 b is qualitatively similar to the behavior seen in momentum space, see Figure 3. However, for w m 2 b the soft function in diagonal space undergoes rapid oscillations, as shown in the right panel. These oscillations do not appear to have a physical interpretation, but they arise as a consequence of the peculiar integral transformation (61) that links functions in the diagonal space with those in momentum space. Beyond the leading order it is difficult to obtain results for the soft function in the diagonal space, because the oscillatory behavior of the transfer function makes it challenging to evaluate the slowly converging integral on the right-hand side of the first relation in (61) numerically.

RG-invariance of the convolution integral T 3
Several of the derivations performed in this paper rely on the assumption that one can pretend that the convolution integral is RG invariant, even though this quantity is ill defined. We will now present arguments indicating that this assumption can indeed be justified. However, we will also see that in order to properly define the quantity T 3 one needs to introduce a rapidity regulator, and that a generic rapidity regularization scheme is in conflict with RG invariance. To see what the problem is, let us evaluate the right-hand side of (77) using the expression for the jet function obtained at leading order in RG-improved perturbation theory [20], which reads where µ j is an appropriate matching scale. This leads to which is undefined. Note that RG resummation can cure the divergence of the first integral w → ∞, but it does not help to regularize the divergent integral over − . We could have chosen to use two different matching scales µ j− and µ j+ for the two jet functions, which would generate an extra factor a Γ (µ j− ,µ j+ ) − but still leaves a divergent integral. The divergence results from the fact that we are integrating over an infinite range of rapidities in the contribution T 3 and it requires regularization.

RG invariance and rapidity regularization
Two different rapidity regularization schemes were considered in [5], and in both cases it was observed that the convolution T 3 by itself is not RG invariant. In that paper the regulators were applied to expression (9), which is written in terms of bare functions. Here we will follow an analogous approach and apply the regulators to expression (77), which is formulated in terms of renormalized functions. In the first scheme, one introduces an analytic rapidity regulator δ and an associated scale ν [32,33]. It has been shown in [5] that a consistent way to do this, which preserves the analytic properties of the h → γγ decay amplitude, is to replace In the second scheme one uses hard cutoffs on the integrals over ± in such a way that The parameter σ should be thought of as being a positive quantity. Only after evaluating the integrals must one take the limit σ → (−1 − i0) by analytic continuation. When deriving the RG equation for the soft function based on the hypothesis of the RG invariance of T 3 it was important that the integrals in (77) run from 0 to ∞ and that the integration measures d ± / ± are invariant under rescalings of these variables. Imposing rapidity regulators as shown above is, in general, not compatible with this derivation. Indeed, using the RG equations for the hard, jet and soft functions given in (21), (22) and (19), along with the relations (25) and (26), we find that after the regulators are imposed the convolution integral T 3 is no longer RG invariant. In the analytic regularization scheme we obtain dT analytic where we have dropped the "−i0" prescription for brevity. In the cutoff scheme we obtain instead dT cutoff In both cases the kernel K is given by Note that "local" terms in the RG equations (those with x = 1) do not contribute in (82) and (83). To show that RG invariance is lost it suffices to work with the lowest-order approximations for all quantities involved, namely 4 We then obtain In the limit δ → 0 the first expression gives the same result as the second one. Can one find a better rapidity regularization scheme that preserves RG invariance? A third possibility, in which a cutoff |y| < y cut is placed on the rapidity y = 1 2 ln + − of the soft momentum, is explored in Appendix E. It offers a small improvement on the schemes discussed above in the sense that the breaking of RG invariance occurs first at two-loop order. Yet the fundamental question posed above remains. The key to answering it is the observation that the breaking of RG invariance is linked to the non-locality of the kernel function K(x, µ), i.e. the fact that it does not vanish for x = 1. In Section 7 we have constructed the diagonal space, in which the evolution is strictly local in w and hence K(x, µ) vanishes for x = 1. In this space the integrand of the convolution (69) is RG invariant at each point in the ( + , − ) plane, and thus RG invariance is preserved for all of our rapidity regularization schemes.
The fact that one can find RG-invariant regularization schemes implies that the logic followed in this paper is self consistent. Assuming that T 3 is RG invariant is then a sensible and well-defined thing to do, and from this hypothesis we have derived the two-loop evolution equations studied in this work. Good things never come without a price, though. While imposing rapidity regulators in the diagonal space does not destroy RG invariance, it does interfere with the cancellation of the dependence on the auxiliary scale ρ introduced in the construction of the diagonal space. The differential equations (60) and (74) are then no longer sufficient to guarantee that T 3 as defined in (69) is independent of ρ. This is not a serious problem, however, since T 3 by itself is not a physical quantity. It is only one out of three terms in a factorization theorem.

Regularized convolution T 3 in the diagonal space
We have just shown that it is possible to define an RG-invariant convolution T 3 by applying the rapidity regularization in the diagonal space. What we have not yet demonstrated is that this convolution is well defined. At fixed order in perturbation theory the integral over w diverges at infinity. 5 We will now show that this endpoint divergence is removed after RG improvement. It suffices to show this at leading order in RG-improved perturbation theory, because higher-order corrections do not change the functional form of the solutions in a significant way. For concreteness, we employ the analytic regularization scheme (80), where the regularized convolution T 3 in the diagonal space takes the form We now use the following solutions for the various component functions, valid at leading order in RG-improved perturbation theory: for the RG functions as well as relation (26) to bring the answer to the explicitly µ-independent form T analytic (90) 5 There is no divergence for w → 0, because the soft function vanishes at the origin.
The integral in the last line can be performed by changing variables from w to z = xw/m 2 b and using the property (68) of the transfer function. This gives Taking the limit δ → 0, we obtain the final result T analytic If the calculation is done in the cutoff regularization scheme (81) one finds the same result without the 1/δ pole and with ν 2 replaced by (−M 2 h ). Note that the integral over z in (91) converges for a + δ > 0 only, which for δ → 0 would require an unphysical choice µ j < µ s . Nevertheless, the above expression can be used to consistently define the result for any choice of scales by analytic continuation in a. The 1/δ pole and the associated logarithm involving the scale ν result from the rapidity divergence and need to cancel against corresponding terms in other contributions to the factorization theorem for the h → γγ decay amplitude.
To summarize this discussion, we find that the regularized convolution T 3 in the diagonal space is RG invariant to all orders of perturbation theory and free of endpoint divergences. This observation explains a posteriori why our heuristic arguments used in Sections 3 and 4 led to consistent results for the renormalization and scale evolution of the soft function. In the resummed expression for the convolution T 3 given above all large logarithms are contained in the RG functions S(µ 1 , µ 2 ) and a i (µ 1 , µ 2 ) except for a single power of the rapidity logarithm. For natural choices of the matching scales µ s ≈ m b and |µ h | ≈ M h the terms in the second line are free of large logarithms. As a final comment, we emphasize that expression (92) has no well-defined fixed-order expansion, because the limit a → 0 does not exist. The reason is that the endpoint divergence of the original integral in (13) has been removed by resummation. When one attempts to "undo" the resummation by expanding the result in powers of α s one encounters terms of order 1/α 2 s and 1/α s . These terms must cancel against similar terms contained in the other two contributions (T 1 and T 2 ) to the factorization theorem for the h → γγ amplitude.

Conclusions
In this work we have presented a detailed analysis of the renormalization and the scale evolution of the soft-quark soft function S(w, µ), defined in terms of the discontinuity of the soft quark propagator dressed by two finite-length Wilson lines connecting at one point. This function appears in the factorization formula for the h → γγ decay amplitude induced by loops of light quarks, recently worked out by two of us [5]. The results we have obtained and the techniques we have developed have a more general importance in the context of understanding SCET factorization theorems at subleading order in power counting, where soft functions containing soft quarks dressed by Wilson lines become a generic feature. The relevance of our results thus extends beyond the practical purpose of studying the h → γγ process.
A central argument of our work has been the hypothesis that the third term T 3 in the h → γγ factorization theorem, which involves a hard function and a double convolution of two radiative jet functions with the soft-quark soft function, should be RG invariant. This does not necessary have to be realized, but it is suggested by certain observations made in [5]. Based on this assumption we have derived the non-local renormalization factor Z S of the soft function at one-loop order in QCD, and we have shown that it successfully removes all the 1/ n poles of the bare soft function. The cancellation is highly non-trivial due to the non-local structure of the counterterms and the fact that the leading-order bare soft function is not simply a constant. This observation puts our hypothesis on firmer grounds. From this result we have derived the one-loop anomalous dimension of the soft function, which closely resembles the structure of the anomalous dimension of the leading-twist B-meson light-cone distribution amplitude. Pushing further, we have used existing results for the two-loop anomalous dimensions of the hard and jet functions to present a conjecture for the two-loop anomalous dimension of the soft function. It would be highly desirable to test this result by a direct two-loop calculation of the soft function.
Using our expression for the two-loop anomalous dimension we have presented an analytic closed-form solution to the non-local RG evolution equation satisfied by the soft function in momentum space. The result obtained at NLO in RG-improved perturbation theory, shown in (42), involves integrals over Meijer G-functions. We find that after RG evolution the discontinuous behavior of the soft function at w = m 2 b seen in fixed-order perturbation theory is smoothed out. We have also studied the asymptotic behavior of the soft function for large values w m 2 b and emphasized the need for a dynamical choice of the soft matching scale in this context. We have then studied the RG evolution of the soft function in Laplace space and in the so-called diagonal space, where its RG evolution is local in the w variable. The construction of the diagonal space in higher orders of perturbation theory requires a non-trivial extension of the dual-space formalism developed in [10,11]. We have explicitly constructed the transfer functions connecting the diagonal space with momentum space. Beyond the leading order in perturbation theory, functions in the diagonal space depend on an auxiliary scale ρ, which cancels in predictions for physical quantities. We have derived the differential equations governing the dependence on ρ.
Finally, we have studied the structure of the double convolution integral T 3 in more detail, finding that it requires rapidity regulators in order to be well defined. Using three different rapidity regularization schemes we have shown that introducing the regulators in momentum space breaks the RG invariance of T 3 . However, RG invariance can be preserved by imposing the rapidity regulators in the diagonal space. We have presented an explicit expression for T 3 at leading order in RG-improved perturbation theory, in which large logarithms of the scale ratio (−M 2 h /m 2 b ) are resummed to all orders of perturbation theory. The resummation improves the behavior of the convolution at large momenta and tames an endpoint divergence for w → ∞.
The results obtained in this work constitute an important step toward understanding the many subtleties and complexities of SCET factorization beyond the leading order in power counting. Not only will they help to resum the large logarithms on the h → γγ decay amplitude beyond the leading double-logarithmic approximation; more generally, we are confident that the techniques of solving the difficult non-local RG equations of soft functions developed in this work will find applications in many other factorization theorems of interest.
where r = α s (µ)/α s (µ s ). The function a γs (µ s , µ) is given by an analogous expression. Whereas the two-loop anomalous dimensions and β-function are required for a Γ and a γs , the expression for the Sudakov exponent S also involves the three-loop coefficients Γ 2 and β 2 .
We now list relevant coefficients of the anomalous dimensions and the QCD β-function, quoting all results in the MS renormalization scheme. For the convenience of the reader, we also give numerical results for N c = 3 and n f = 5. The two-loop cusp anomalous dimension Γ cusp was obtained long ago [22], while the three-loop coefficient was derived in [34]. The results are T F n f ≈ 36.8436 , The anomalous dimension γ s follows from (26) and can be determined up to two-loop order using the explicit expressions for γ and γ q given in [5] and [21,35], respectively. We find ≈ −151.280 . (A.4) Finally, the expansion coefficients for the QCD β-function to three-loop order are We choose to work with n f = 5 active quark flavors, because we are mainly interested in RG evolution to a scale µ between the mass scales of the bottom and top quarks.
In the analysis in Section 5 we need the ratio of the running b-quark mass m b (µ) defined in the MS scheme and the pole mass m b . At NLO in RG-improved perturbation theory this ratio is given by (A.6) where the one-and two-loop coefficients in the anomalous dimension of the quark mass are given by [12] γ m,0 = −6C F , γ m, B Light-cone distribution amplitude of the B-meson The RG evolution equation (19) and the associated anomalous dimension γ S in (25) are closely related to the corresponding equations for the leading-twist light-cone distribution amplitude φ B + (ω) of the B meson defined in heavy-quark effective theory [24]. This quantity obeys the differential equation where [20,23,25] γ + (ω, ω ; µ) = − Γ cusp (α s ) ln (B. 2) The anomalous dimension γ S in (25) differs from γ + by the argument of the logarithm (w/µ 2 versus ω/µ) and the factor 2 in front of the non-local terms. As a result of this, one finds that the solution to the evolution equation (B.1) can be written in terms of hypergeometric functions [26] rather than the more complicated Meijer G-functions in (53).

C Solution in terms of hypergeometric functions
In the derivations in Sections 4.2 and 7.2 we have written some results in terms of Meijer Gfunctions. Here we relate these functions to the more familiar hypergeometric functions. We begin with the solution to the RG equation of the soft function. For a Γ (µ s , µ) < 0 it is possible to perform the contour integral in (38) using the theorem of residues. The integrand contains double and single poles in the complex η-plane located at η = −n and η = n + a Γ (µ s , µ), where n ∈ N. For w > w (w < w ), we can close the contour in the upper (lower) half plane and pick up the residues of the single poles. We assume that a Γ (µ s , µ) > −1, such that all of the poles in the second series lie in the lower half plane. The expression for the residues of the single poles involves Γ-functions and their derivatives. Comparing the answer with (39) and denoting z = min(w, w )/max(w, w ) < 1, we find that Here H(a) = ψ(1 + a) + γ E is the harmonic-number function, and the functions G i are given by G 1 (z, a) = 4 F 3 (1+a, 1+a, 1+a, 2+a; 2, 2, 2; z) , G 2 (z, a) = 4 F 3 (1+a, 1+a, 2+a, 2+a; 1, 2, 2; z) , G 3 (z, a) = (∂ p 1 + ∂ q 1 ) 4 F 3 (1+a, 1+a, 2+a, 2+a; 1, 2, 2; z) .

(C.3)
The derivatives in the last expression act on the indices of the hypergeometric function 4 F 3 (p 1 , p 2 , p 3 , p 4 ; q 1 , q 2 , q 3 ; z). The functions G i (z, a) are real-valued for z < 1 and are singular in the limit z → 1 − . Using the properties of the hypergeometric functions we have derived the precise form of this singularity shown in (40). We have also encountered the Meijer G-function in the calculation of the transfer function F S (x, µ) needed for the calculation of the soft function in the diagonal space. At leading order one can ignore the exponential containing the integral over the function H in (61). The integrand of the η-integral then has double and single poles at values η = −n with n ∈ N. Evaluating the integral using the theorem of residues, we obtain G 2,0 0,4 1, 1, 0, 0 x = x 2 0 F 3 (2, 2, 2; x) − ln x + 4γ E + 4∂ q 1 0 F 3 (1, 2, 2; x) . (C.4)

D Soft function in the dual space
The dual soft function constructed following the approach of [10,11] is related to the original soft function via where the transfer function F S,0 (x) has been given in (65). As a consequence of the fact that no higher-order corrections to the transfer function are included in this approach, the function s dual (w, µ) obeys the non-local RG evolution equation ( 2) with the simplification that the Γ-functions appearing in (38) are now absent. The integral over η thus evaluates to a δ-function rather than a Meijer G-function, yielding s dual (w, µ) = U S (w; µ, µ s ) 3) It can readily be checked that this solution obeys the RG equation (53) at two-loop order. This result is simpler than the corresponding relation (42) obtained in momentum space, but not as simple as the solution (59) found in the diagonal space.

E Rapidity cutoff scheme
In Section 8.1 we have evaluated the convolution integral T 3 using two different rapidity regularization schemes. As an interesting and rather natural alternative we now consider a third scheme, in which a cut |y| < y cut is imposed directly on the rapidity y = 1 2 ln + − of the soft momentum. When this is done, the regularized convolution T 3 takes the form in this regularization scheme scheme.