Two-Loop Radiative Jet Function for Exclusive $B$-Meson and Higgs Decays

The rare radiative $B$-meson decay $B^-\to\gamma\ell^-\bar\nu$ and the radiative Higgs-boson decay $h\to\gamma\gamma$ mediated by light-quark loops both receive large logarithmic corrections in QCD, which can be resummed using factorization theorems derived in soft-collinear effective theory. In these factorization theorems the same radiative jet function appears, which is a central object in the study of factorization beyond the leading order in scale ratios. We calculate this function at two-loop order both in momentum space and in a dual space, where its renormalization-group evolution takes on a simpler form. We also derive the two-loop anomalous dimension of the jet function and present the exact solution to its evolution equation at two-loop order. Another important outcome of our analysis is the explicit form of the two-loop anomalous dimension of the $B$-meson light-cone distribution amplitude in momentum space.


Introduction
Soft-collinear effective theory (SCET) is a convenient tool to study the factorization properties of cross sections and scattering or decay amplitudes sensitive to several hierarchical energy scales [1][2][3][4][5]. Typically, 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, whereas the jet or soft functions are defined in terms of matrix elements of collinear or soft fields in the low-energy effective theory.
Jet functions -matrix elements of non-local products of collinear fields -play an important role in these factorization theorems. They often live at an intermediate scale, which lies between the hard-scattering scale of the process and lowest energy scale it is sensitive to. The most familiar jet function is defined as the spectral function (i.e., the discontinuity) of the quark propagator dressed by a light-like Wilson line connecting the two quark fields [6]. This jet function appears in a large variety of phenomenological applications. For example, it enters in the factorization theorems for inclusive B-meson decays into light final-state particles, such asB → X s γ andB → X u −ν [7,8]. It also appears in the resummation of threshold logarithms in deep inelastic scattering [9,10].
Recently, there has been a growing interest in understanding factorization at subleading power in scale ratios. Beyond the leading power a large variety of hard, collinear 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. In the notation of [4], it reads where ξ n is a collinear quark spinor subject to the constraint / n ξ n = 0, W n is a collinear Wilson line, D µ n is a covariant collinear derivative (containing collinear gauge fields) acting on collinear fields, and q s describes a soft quark. The collinear particles carry large momentum flow along a light-like direction n µ . The soft quark field must be multipole expanded for consistency, and we denote x µ − = (n · x) n µ 2 , wheren µ is a conjugate light-like vector satisfying n ·n = 2 (see [11] for a pedagogical introduction to SCET).
New jet functions can be defined in terms of the matrix elements of collinear fields in the presence of one insertion of this subleading Lagrangian [12][13][14]. These functions are called "radiative jet functions" [15][16][17], since they involve the emission of a soft particle from inside a jet. The radiative jet functions are of general interest not only because they appear in the description of power corrections to established factorization theorems. There exist interesting physical processes that are sensitive to soft quark exchange already at leading order in the expansion of the relevant decay or scattering amplitude. One of the first appearances of a radiative jet function J(p 2 ) related to collinear interactions with a soft quark appeared in the theoretical description of the exclusive, radiative B-meson decay B − → γ −ν [18,19]. In this process the soft spectator quark of the B meson couples to a collinear photon and an off-shell collinear quark, which then connects to the weak-interaction vertex, where it annihilates the b-quark and turns into a virtual W − boson. Interestingly, the same radiative jet function has recently been encountered in a completely different context: in the theoretical description of the contribution to the radiative Higgs-boson decay h → γγ that is induced by loops of light b-quarks [20]. (This is not the dominant contribution to the decay amplitude, but it is a particularly interesting one with regard to its factorization properties.) Following our recent work [20], we define the radiative jet function J(p 2 ) in terms of the matrix element where p s is the momentum carried away by the soft quark in (1), and p ≡ k + p s+ . We denote by e q the electric charge of the quark, which emits the final-state photon. Note that p s · x − = p s+ · x with p µ s+ = (n · p s )n µ 2 . The jet function J(p 2 ) depends on the only non-trivial kinematic invariant (note that k 2 = 0 and p 2 s = 0) In Section 2 we present a detailed discussion of some general properties of the radiative jet function J(p 2 ), with a special focus on its behavior under renormalization-group (RG) evolution. We derive the two-loop anomalous dimension of the jet function and present exact solutions to its RG evolution equations both in momentum space and in the so-called dual space. The technique we employ for obtaining these solutions is general and can be applied to other radiative jet functions, too. Section 3 contains a description of the calculation of the jet function at two-loop order in QCD. The renormalization of the jet function is discussed in Section 4. In Section 5 we briefly comment on phenomenological implications of our results in the context of B − → γ −ν decay. We then present our conclusions.

General properties the radiative jet function
We begin by reviewing and deriving some general properties the jet function, some of which are based on insights that were uncovered a long time ago, while several others are new.

One-loop expressions
At one-loop order, we find that the bare jet function in d = 4 − 2 spacetime dimensions reads Renormalizing the bare coupling in the MS scheme, where β 0 = 11 3 C A − 4 3 T F n f is the first coefficient of the QCD β-function, one obtains Here and below α s ≡ α s (µ) always denotes the renormalized coupling. For simplicity, we will from now on drop the "−i0" prescription, which defines the sign of the imaginary part of the jet function in the time-like region, where p 2 > 0. While at one-loop order one could renormalize the jet function by means of a local counterterm, the correct renormalization factor has a more complicated non-local form. 1 The proper renormalization condition has been derived from the consistency of the factorization formula for the B − → γ −ν decay amplitude, requiring that the amplitude be independent of the renormalization scale [19]. In this process, the known RG equations for the B-meson light-cone distribution amplitude (LCDA) [21] and some other quantities have been used. The jet function depends on a single argument p 2 , which can be either time-like or space-like. The time-like (space-like) jet functions belonging to different p 2 > 0 (p 2 < 0) values mix under renormalization, but there is no mixing between the time-like and space-like jet functions. Treating both cases at the same time, we write the renormalization condition in the form At one-loop order, one finds where the symmetric distribution arises in the so-called Lange-Neubert kernel for the B-meson LCDA [21] (see also [22]). The plus prescription is defined such that, when Γ(y, x) is integrated with a function f (x), one must replace f (x) → f (x) − f (y) under the integral. At one-loop order the plus distribution has no effect when the renormalized jet function is derived from (7). One finds This result was first obtained in [18,19]. One of the main goals of this paper is to calculate the two-loop corrections to this formula.

Renormalization-group evolution
The renormalized jet function obeys the RG evolution equation where the anomalous dimension is as usual obtained from the coefficient of the single 1/ pole in Z J via Without loss of generality, we express the result in the form [19] where we have identified the coefficient of the logarithmic term with the light-like cusp anomalous dimension Γ cusp (α s ) in the fundamental representation of SU(N c ), a central quantity in the theory of the renormalization of Wilson loops with cusps [23,24]. Since the plus distribution is linked with the logarithmic term, it is multiplied by the same quantity. The non-logarithmic term γ (α s ) of the local part of the anomalous dimension vanishes at one-loop order [19]. Its two-loop expression will be derived here for the first time.
The first three terms in the anomalous dimension γ J retain their form to all orders in perturbation theory. The most general form of the local terms is a linear function of ln(−p 2 /µ 2 ), and the cusp anomalous dimension is the coefficient of the logarithmic term. However, the form of the non-local terms (with x = 1) is presently only known at one-loop order. The O(α 2 s ) corrections indicated in (13) thus refer to higher-order non-local terms, which are presently unknown. For the more familiar jet function entering inclusive processes such as the rare inclusive decayB → X s γ, such higher-order corrections are known to be absent [6], i.e. the functional form of the one-loop anomalous dimension is preserved in higher orders, and the non-local terms are determined completely by the cusp anomalous dimension. We will see, however, that further non-local higher-order terms do exist in the case of the exclusive jet function in (2).
For the remainder of this section, we will ignore the unknown higher-order non-local terms in (13), but we will keep the remaining terms at arbitrary order in perturbation theory. It is then possible to derive exact solutions to the RG evolution equation (11) using a technique developed in [19,21]. It is based on the observation that on dimensional grounds the integral (the variables x and y can carry any mass dimension) evaluates to a dimensionless function of the exponent a. Here ψ(z) = Γ (z)/Γ(z) is the digamma function. It can then be checked that the ansatz with α s (µ α ) ≡ α provides a solution to the RG equation with the initial condition (−p 2 /µ 2 j ) η at some matching scale µ = µ j , at which J(p 2 , µ j ) is assumed to be free of large logarithms. Here β(α s ) = dα s (µ)/d ln µ is the QCD β-function, and we have defined the RG functions which are the solutions to the equations The function a γ (µ j , µ) is defined analogously to a Γ (µ j , µ). Note that both S(µ j , µ) and a Γ (µ j , µ) take negative (positive) values if µ > µ j (µ < µ j ), since the cusp anomalous dimension is positive. Explicit expressions for these objects obtained at next-to-next-to-leading order (NNLO) in perturbation theory can be found in the appendix of [10]. At any fixed order in perturbation theory, the renormalized jet function at the matching scale depends on p 2 only via powers of the logarithm L p = ln(−p 2 /µ 2 j ). We can generate these logarithms by taking derivatives with respect to η. Hence, with the definition the exact solution of the evolution equation can be written in the closed form [19] We now change integration variables in the exponent of the last term from α to = a Γ (µ j , µ α ). This yields where a ≡ a Γ (µ j , µ). This leads to the final result [19] An alternative solution of the RG equation (11), which works for more general initial conditions (even though this is not needed for the case at hand), can be obtained by taking a Fourier transform of the jet function with respect to ln(−p 2 /µ 2 ), i.e.
Since the right-hand side of this equation exhibits a power-like dependence on (−p 2 /µ 2 ), one can use the technique described above (with η replaced by it) to show that for general initial condition the evolution equation (11) is solved by This solution forms the basis of the construction of the jet function in the so-called dual space (see below). For the case of downward scale evolution, for which µ < µ j , it is possible to evaluate the integral over t in closed form. This is discussed in Appendix A.

Jet function in the dual space
For the case of the B-meson LCDA, it has been shown in [25] that one can bring an RG equation with an anomalous dimension of the type shown in (13) -in the approximation where unknown, non-local contributions to the anomalous dimension arising at two-loop order and higher are neglected -to a much simpler form using a suitably chosen integral transform. Adapted to our case, the key observation based on (23) is that the function has a particulary simple behavior under RG evolution. Shifting the integration variable in (23) from t to t = t − ia Γ (µ j , µ), one finds that Defining a dual jet function j(p 2 , µ) via the Fourier transform we then obtain This dual function obeys the local RG equation We stress again that this equation only holds in the approximation where the unknown nonlocal contributions to the anomalous dimension (13) are neglected. Equation (41) below shows the generalization required when these terms are included at two-loop order.
The relation between the original function and the dual function can be derived by combining (26) and (A.1). This gives The integrand of the t-integral has poles at values t = in with n ∈ N. Evaluating the integral using the theorem of residues, one obtains where J 1 (x) is a Bessel function. One thus finds the integral transforms [25] where the second relation follows from the orthonormality condition At one-loop order, we find that In the beautiful papers [26,27] it was shown that the Lange-Neubert kernel for the Bmeson LCDA can be written in a remarkably compact form as a logarithm of the generator of special conformal transformations along the light-cone. Using tools from conformal field theory, the above-mentioned transformation of the evolution equation (11) to the local form (28) was rederived. In subsequent work by the same authors [28,29] the evolution equation was extended to two-loop order. It was found that, starting at O(α 2 s ), non-local terms appear in the dual space as well. We will use these results in our two-loop analysis below.

Two-loop evolution of the jet function
The two-loop RG equation obeyed by the jet function in the dual space can be derived from the QCD factorization theorem for the decay B − → γ −ν valid in the region of large photon energy (E γ m b /2). At leading power in Λ QCD /m b , the corresponding decay amplitude can be written in the factorized form [18,19] where F B is related to the B-meson decay constant in the heavy-quark limit (F B ≈ f B √ m B modulo radiative corrections), H is a hard-scattering function depending on the short-distance scales m b and 2E γ , and φ B + denotes the leading-twist LCDA of the B meson depending on a variable ω = O(Λ QCD ) [22]. J is the jet function discussed above, with p 2 < 0 in the space-like region. This function depends on an intermediate scale of order 2E γ ω ∼ m b Λ QCD . In the dual space, the right-hand side of (34) takes an identical form [25], i.e.
where the dual function j is related to the original jet function by the first relation in (31), and ρ + is defined in an analogous way as In [29], the RG evolution equation for the function η + (s, µ), which is related to ρ + (ω, µ) via s η + (s, µ) = ρ + (1/s, µ), was derived at two-loop order. Interestingly, it was observed that at this order the evolution of the dual function is no longer local. Rather, it was found that (we present the equation for ρ + rather than η + ) where arises from conformal symmetry breaking. We can now use the known RG equations for the B-meson decay constant in heavy-quark effective theory [30] and of the hard-scattering function [31][32][33][34] to derive the two-loop evolution equation for the jet function in the dual space. We obtain where The two-loop expressions for the anomalous dimensions on the right-hand side of this relation are listed in Appendix C. Using these results, we obtain This quantity is genuinely non-abelian; it starts at two-loop order and has no C 2 F term. Interestingly, we find that γ (α s ) = −γ W (α s ) coincides, up to a sign, with the anomalous dimension γ W of the Drell-Yan soft function derived in [35]. It would be interesting to explore the nature of this connection in more detail.
Starting from (41), we can apply the second transformation rule in (31) to derive the explicit form of the RG evolution equation (11) for the jet function in momentum space at two-loop order. The relevant anomalous dimension γ J is given by where γ dual J denotes the anomalous dimension in the dual space, which according to (41) is given by Using the orthonormality condition (32), we obtain from (44) The integral in the second line diverges for x → 1 and must be evaluated in the sense of distributions by studying its action on a smooth test function f (x). We find Using this result, we obtain This is the desired extension of relation (13) to two-loop order.
As an important tangential outcome of our analysis, we now derive the explicit form of the two-loop RG evolution equation of the B-meson LCDA in momentum space. Because of the structural similarity of the evolution equations (37) and (41) in the dual space, we can apply the same method as above to obtain This is the desired extension of the Lange-Neubert kernel to two-loop order.

Solutions to the two-loop evolution equations
The evolution equation (41) in the dual space can be solved using the same technique we have adopted in Section 2.2. The key observation is that any power (−p 2 ) a of the momentum squared is an eigenfunction of the evolution kernel. To see this, note that defines a dimensionless function of the exponent a, where It follows that the function provides a solution to (41) with the initial condition (−p 2 e −2γ E /µ 2 j ) η at the matching scale µ = µ j . We now use the fact that the initial condition j(p 2 , µ j ) depends on p 2 only through powers of the logarithmL p = ln(−p 2 e −2γ E /µ 2 ). Writing we find in analogy with (21) that the general solution to the evolution equation (41) in the dual space is obtained as where we have used that β(α s ) = −β 0 α 2 s /(2π) + . . . at leading order. It is not difficult to transform this solution back to momentum space. Using the second relation in (31), we find that In the final step, we employ the relation between the functionsĴ(∂ η , µ j ) defined in (53) and J(∂ η , µ j ) defined in (18), which follows by setting µ = µ j in the above solution. We conclude that the general solution to the momentumspace RG evolution equation (11) is given by The explicit expressions for the two-loop anomalous dimensions (45) and (48), as well as the explicit solutions of the corresponding evolution equations (54) and (57), are among the most important new results obtained in this paper.

Two-loop calculation of the bare jet function
The radiative jet function can be directly evaluated from its definition in (2). In doing so, we work to leading order in the electromagnetic coupling e but include higher-order QCD corrections. It is useful to recast the original definition in terms of so-called "gauge-invariant collinear building blocks" defined as [3,36] Note that the effective photon field A µ n and the gluon field G µ n , which contain the gauge couplings in their definition, are separately gauge invariant. The definition of the jet function in (2) therefore defines two gauge-invariant objects J A (p 2 ) and J G (p 2 ) via A third contribution involving the derivative / ∂ ⊥ acting on X n (x) vanishes, because this field carries vanishing transverse momentum.
Since we work to leading order in electromagnetic interactions, the function J A is equivalent to the collinear quark propagator dressed by light-like Wilson lines. This object was studied in [6], where its discontinuity at two-loop order was presented. The function J A can be determined from this result in a straightforward way. The two-loop calculation of the new jet function J G can be performed using similar methods. As shown in [6], the propogator in (2) can be rewritten in terms of standard QCD fields, because the collinear SCET Lagrangian (without couplings to soft fields) is equivalent to the original QCD Lagrangian [4]. We therefore use QCD Feynman rules for convenience.
In principle, the matrix element for J G can be calculated in a general covariant gauge. However, the Feynman rules for vertices derived from the collinear gluon field G µ (x) are rather complicated due to the Wilson lines contained in its definition. In light-cone gauge n · A(x) = 0, on the other hand, the field G µ (x) = g s A µ (x) takes on a very simple form, since the Wilson lines become trivial (W n = 1). The smaller number of Feynman diagrams and the absence of ghost contributions result in a more efficient computation of J G in this gauge. The free gluon propagator with momentum l µ in light-cone gauge is given by where we do not adopt the Mandelstam-Leibbrandt prescription to regularize the singularity atn · l = 0 (see [37] for more details). Figure 1 illustrates the non-vanishing two-loop Feynman diagrams contributing to J G in light-cone gauge. After performing simplifications of the Dirac, Lorentz and color algebras, these diagrams can be transformed into linear combinations of scalar Feynman integrals belonging to one of five different integral topologies, each of which contains up tp seven linearly independent squared propagators and up to two linear propagators. Mapping the Feynman integrals to specific integral topologies requires partial-fraction decompositions on linear propagators followed by shifts of the loop momenta. The five integral topologies can be cast into the form where A 1 = 9 i=1 a i and A 2 = 12 i=10 a i . The denominators D n are defined as (omitting the "−i0" terms for brevity) The five integral topologies are given by restricting the propagator index in (61) as follows: topology 1: a 5 = a 8 = a 12 = 0 , topology 2: a 5 = a 8 = a 11 = 0 , topology 3: a 5 = a 8 = a 10 = 0 , topology 4: a 3 = a 8 = a 12 = 0 , topology 5: a 3 = a 6 = a 12 = 0 . (63) We use the public program FIRE5 [38] to perform the integration-by-parts (IBP) reduction of the integrals in each of the five topologies. Besides the ten IBP identities in each topology, obtained by inserting differential operators in front of the integrands, additional linear algebraic identities can be derived by means of a kinematic constraint. The fact that the soft momentum p s+ = (p − k) in (2) is light-like leads to the relation By contracting both sides of this identity with the loop momenta l µ 1 and l µ 2 , two classes of linear identities can be derived. They are a − 1 − a − 7 + p 2 n · p a − 10 = 0 for topologies 1, 2, 4 and 5 , and −a − 1 + a − 3 + a − 7 − a − 9 + p 2 n · p a − 11 = 0 for topologies 1 and 3 , Here the operator a − n lowers the index a n on the n th denominator in (61) by one unit. After implementing the above linear identities in FIRE5, the scalar Feynman integrals are reduced to ten master integrals (MIs), all of which can be remapped onto topology 4 by shifting the loop momenta.
To compute the MIs analytically, we use a method inspired by [39][40][41], which has also been employed in the recent three-loop calculations of the soft and jet functions arising in the factorization theorems for the inclusive decaysB → X s γ andB → X u −ν [42,43]. The basic strategy is to first map each MI in d = 4 − 2 spacetime dimensions to an analytically calculable quasi-finite integral in higher dimension d = 6 − 2 or d = 8 − 2 , and then determine the linear relations between the MIs and the corresponding quasi-finite integrals by dimensional recurrence relations [44][45][46] and IBP reduction. The quasi-finite integrals are free of divergences from the Feynman parameter integrations. They can be found by observing that raising the dimension by an even number decreases the degree of infrared divergences, while increasing appropriate propagator indices by integer amounts decreases the degree of ultraviolet divergences. The expansions of the quasi-finite integrals are linearly reducible and can be analytically evaluated by the powerful Maple package HyperInt [47]. We use the public code LiteRed [48,49] to determine the dimensional recurrence relations. After further IBP reduction, we construct the linear relations between the MIs in d + 2 and d dimension in the form I(d + 2) = A(d) · I(d). This allows us to build linear relations between the MIs in d = 4 − 2 and the corresponding quasi-finite integrals in higher dimension. Then the analytical results of the MIs are obtained by solving the linear equations.
With the integrals at hand, the jet function J G is obtained by evaluating the Feynman diagrams shown in Figure 1. Adding to this result the contribution from J A , we obtain for the bare jet function at two-loop order where the two-loop coefficients are given by (68) In Appendix B we present our results for the two jet functions J A and J G separately, including terms up to and including O( 2 ).
Given the above result for the bare jet function in momentum space, it is straightforward to derive the corresponding expression for the bare jet function in the dual space. The first relation in (31) implies that this function can be obtained from (67) by means of the replacement dimension, one can construct the renormalization factor Z dual J in the dual space, defined in analogy with (7), using a general relation valid for Sudakov problems, in which the anomalous dimension contains an explicit dependence on ln µ [50]. Applied to our case, it reads .
This exact solution must be expanded in powers of 1/ (this is formally an expansion about = ∞) to generate the renormalization factor in the MS scheme. When one exponentiates this expression to obtain Z dual J itself, one encounters higher powers of the anomalous dimension. These involve additional integrals, e.g.
At two-loop order, we obtain whereL p = ln(−p 2 e −2γ E /µ 2 ). The relevant expansion coefficients Γ n and γ n of the anomalous dimensions are listed in Appendix C. Applying this result to the bare jet function in the dual space, we find that indeed all 1/ poles cancel out. This provides a non-trivial consistency check on the two-loop results for γ η (α s ) in (42) and h(x) in (38), which were obtained in [29]. The result for the renormalized dual jet function reads Given the above result, we can obtain the jet function in momentum space by applying the second integral transformation in (31). It follows that we must perform the replacementŝ where now L p = ln(−p 2 /µ 2 ). This gives the final result with The two-loop expressions (73) and (76) for the radiative jet functions in the dual space and in momentum space are important new results obtained in this paper.
We have derived the two-loop result (76) for J(p 2 ) in a second, totally independent way by performing a two-loop calculation of the hard matching coefficients H 2 (z) and H 3 in the factorization formula for the h → γγ decay amplitude derived in [20]. Further details about this rather difficult calculation will be presented elsewhere. In the limit where z → 0, these coefficients were shown to obey the refactorization formula The two-loop result we have obtained for the bare jet function from this relation is in complete agreement with the expression given above. This not only confirms our result for the jet function but also provides a non-trivial test of the refactorization formula (78).

Phenomenological impact of the two-loop corrections
The two-loop jet function we have calculated in this paper is the last missing ingredient for a calculation of the leading-power contributions to the B − → γ −ν decay amplitude at NNLO in RG-improved perturbation theory. This is important, since this process provides the most direct information about the properties of the B-meson LCDA [51]. The potential impact of power corrections in Λ QCD /m b has been studied in [52][53][54].
In order to estimate the potential impact of the two-loop corrections to the jet function, we recall that, as shown in (34), the decay amplitude is proportional to the convolution where E γ is the photon energy in the B-meson rest frame. For the purposes of illustration, we fix the renormalization scale at a value µ = µ j ≈ 1.5 GeV, corresponding to a typical matching scale for the jet function. The LCDA naturally lives at a lower scale µ 0 of order 1 GeV, but since in practice there is no large scale hierarchy, we will not resum logarithms of the ratio µ j /µ 0 . Following [51], we define the hadronic matrix elements where µ 0 = 1 GeV is a fixed reference scale, which is part of the definition of the logarithmic moments. For simplicity, we choose the matching scale µ j such that ln(2E γ µ 0 /µ 2 j ) = 0. We then obtain at two-loop order (with n f = 4 light quark flavors) where all hadronic parameters are defined at the scale µ j . We have chosen α s (µ j ) = 0.35 as a reference value for the strong coupling, which corresponds to µ j ≈ 1.5 GeV. On general grounds, one expects the logarithmic moments σ n to be of O(1), even though in concrete models for φ B + (ω) such as the simplest exponential model [22] the higher moments defined with µ 0 = 1 GeV tend to take on much larger values. We observe that the two-loop corrections are not significantly smaller than the corrections arising at one-loop order and hence should be included in future analyses. For example, using the central values of the default choices σ 1 = 1.5 ± 1 and σ 2 = 3 ± 2 adopted in [51], one finds We leave a complete analysis of the B − → γ −ν decay rate including higher-order perturbative corrections for future work.

Conclusions
We have presented a detailed study of the radiative jet function J(p 2 ) defined in (2), which plays a central role in the theory of factorization at subleading power in scale ratios. This object appears in factorization theorems for important exclusive processes such as the rare B-meson decay B − → γ −ν and the contributions to the radiative Higgs-boson decay h → γγ mediated by light-quark loops. In the first case, in particular, a precise knowledge of the jet function is a prerequisite for extracting accurate information about the B-meson LCDA from experimental data in the region of high photon energy. The B-meson LCDA itself is a central quantity in the theory of QCD factorization applied to exclusive decays of B mesons [55,56]. We have calculated the radiative jet function at two-loop order both in momentum space and in a dual space, in which its RG evolution equation takes on a particularly simple form. We have further derived the anomalous dimensions for the jet functions in momentum and the dual space, including for the first time the so-far unknown two-loop contributions not controlled by the cusp anomalous dimension. Our derivation of these terms has relied on a recent calculation of the corresponding contributions to the anomalous dimension of the B-meson LCDA in the dual space [29]. We find that the quantity γ (α s ) in (13) obeys an unexpected relation with the anomalous dimension of the Drell-Yan soft function [35], which deserves further study. Finally, we have obtained analytic solutions to the two-loop RG evolution equations of the radiative jet function in both spaces. The technique we used for obtaining these solutions is general and can be applied to other radiative jet functions as well.
The results presented in this paper will play an important role in the renormalization of the factorization theorem for the light-quark induced h → γγ decay amplitude derived in [20]. Using the evolution equations of the radiative jet function derived here and of a hard function well known in SCET, we will be able to derive the two-loop evolution equation satisfied by the soft-quark soft function [57], another central object in the theory of factorization beyond the leading power.

A Evolution with a general boundary condition
For the case of downward scale evolution, for which µ < µ j and hence a Γ (µ j , µ) > 0, it is possible to evaluate the integral over t in (23) in closed form using the theorem of residues, adopting a technique developed in [58] for the case of the B-meson LCDA. The integrand contains single poles in the complex t-plane located at t = in and t = −i[n − a Γ (µ j , µ)], where n ∈ N is a positive integer. Expressing the Fourier imageJ(t, µ j ) in terms of the original jet function usingJ we see that for x < 1 (x > 1) the contour can be closed in the upper (lower) half plane, and one then needs to perform the infinite sum over the residues of the poles. Let us assume that a Γ (µ j , µ) < 1, such that all poles in the second series lie in the lower half plane. This condition is always satisfied in practical calculations. At leading order, it is equivalent to the statement that α s (µ)/α s (µ j ) > exp 2β 0 where a ≡ a Γ (µ j , µ). This elegant formula relates the jet function at the scale µ < µ j to the jet function defined at the matching scale µ j , irrespective of what the initial condition is. We stress that this formula cannot be applied to the case of upward evolution (µ > µ j ). In the limit x → 1 and for a < 1 2 , the hypergeometric function behaves like 2 F 1 (1 − a, 2 − a; 2; x) ∼ (1 − x) −1+2a , which generates a non-integrable singularity if a < 0.

B Two-loop results for the functions J A and J G
Here we present our results for the jet functions J A and J G defined in (59), keeping terms of higher order in . We find

C Anomalous dimensions
Here we list expressions for the relevant anomalous dimensions up to two-loop order. We define the perturbative expansion coefficients via Γ cusp (α s ) = Γ 0 α s 4π + Γ 1 α s 4π 2 + . . . , (C. 1) and similarly for all other anomalous dimensions. The coefficients needed in (72) are Γ 0 = 4C F , T F n f , γ 0 = 0 , (C. 2) The coefficient of the anomalous dimensions in (43) are