Three-loop soft function for heavy-to-light quark decays

We compute the 1-jettiness soft function for the decay of a heavy quark into a light quark jet plus colorless particles at three-loop order in soft-collinear effective theory. The 1-jettiness measurement fixes the total small light-cone momentum component of the soft radiation with respect to the jet direction. This soft function is a universal ingredient to the factorization of heavy-to-light quark decays in the limit of small 1-jettiness. Our three-loop result is required for resummation at the N$^3$LL$^\prime$ level, e.g. near the endpoint in the photon energy spectrum of the $B \to X_s \gamma$ decay. It is also a necessary ingredient for future calculations of fully-differential heavy-to-light quark decay rates at N$^3$LO using the $N$-jettiness subtraction method, e.g. for semileptonic top decays. Using our result for the soft anomalous dimension we confirm predictions on the universal infrared structure of QCD scattering amplitudes with a massive external quark at three loops.


Introduction
Heavy quark (bottom, top) decays are phenomenologically important Standard Model (SM) processes. A prominent example is the decay b → sγ. In the SM it is loop-suppressed and represents a promising window to new physics. In particular, beyond SM interactions due to flavor changing neutral currents could add a measurable effect on the decay rate of B → X s γ on top of the small SM background. In the phenomenologically relevant region of large photon energies the decay rate Γ(B → X s γ) factorizes as [1] where ∆ = m b − 2E γ . 1 Within soft-collinear effective theory (SCET) [2][3][4][5][6][7] this factorization theorem was proven in ref. [5]. The hard function H encodes the short distance (electroweak) interaction and its virtual quantum corrections at and beyond the hard scale m b ∼ E γ . Explicit expressions up to two loops can be found in refs. [8,9]. The jet function J describes the collinear radiation in the final state jet initiated by the (massless) s quark and is governed by the virtuality scale √ m b ∆. In ref. [10] we computed the massless quark jet function to three-loop order (see also ref. [11]). Finally, S denotes the B-meson shape function [12,13] which describes the physics at scales smaller or similar to ∆. For ∆ Λ QCD nonperturbative effects are sizable and one can further factorize S into a purely perturbative 'heavy-to-light soft function' S hl and a nonperturbative shape function F [8]: S(ω, µ) = dω S hl (ω − ω , µ) F (ω ) , (1.2) with dω F (ω) = 1. The soft function S hl can be expressed as a partonic b-quark matrix element, see eq. (2.4), and was computed to two-loop order in ref. [14]. The functions J, S hl , and F vanish when their first argument is negative. This entails finite integration ranges in eqs. (1.1) and (1.2). The perturbative factorization functions H, J, and S hl depend individually on the common renormalization scale µ, while the total decay rate in eq. (1.1) is µ independent to the perturbative order one is working at. The renormalization group (RG) evolution of the hard, jet, and, soft functions is therefore not independent, but subject to a consistency relation (which will be relevant later). The combined RG running of the different functions to the common scale µ eventually resums large logarithms of the ratios between the hard, jet, and soft (matching) scales µ H ∼ m b , µ J ∼ √ m b ∆, and µ S ∼ ∆. A factorization theorem analogous to eq. (1.1) also holds for the decay B → X u ν. In particular the involved jet and shape functions are the same. The nonperturbative function F can thus be obtained from a fit to experimental data for the differential spectrum of one or both decays, see e.g. ref. [15], and then be used for theoretical predictions. For details we refer to ref. [8]. The current state of the art for such predictions includes resummation at the primed next-to-next-to-leading logarithmic (NNLL ) level [8,16], where the NNLL expression [17] is augmented with the next-to-next-to-leading order (NNLO) corrections to the factorization functions at their matching scales, i.e. to H(µ H ), J(µ J ), S hl (µ S ). 2 Still, the uncertainties from missing higher-order perturbative corrections represent a major contribution to the total error budget [15].
In order to reach N 3 LL accuracy of the decay rate in eq. (1.1) the three-loop corrections to the hard, jet, and soft, functions along with their anomalous dimensions are required. Our three-loop calculation of the jet function in ref. [10] represents a first step toward this goal. In the present paper we calculate the soft function S hl at three loops, while the three-loop hard function is left for future work. We also give explicit expressions for all three-loop (noncusp) anomalous dimensions necessary for N 3 LL ( ) resummation in eq. (1.1).
Another possible application of our result is within the context of N -jettiness subtractions [19,20]. In its simplest version the latter is an infrared (IR) slicing method which uses the observable N -jettiness T N [21] as an auxiliary resolution variable for soft and collinear real emissions. It was employed amongst others to compute the fully-differential decay rate of the semileptonic top decay t → W + (l + ν)b at NNLO in QCD [22]. In this case the resolution variable is T 1 (1-jettiness). For T 1 < T cut the decay rate is given by a factorization formula analogous to eq. (1.1), but with S = S hl , provided that T cut is small enough to neglect O(T cut /m t ) power corrections at the desired precision. For T 1 > T cut there is at least one additional hard parton in the final state. On the other hand quantum corrections to this part of the decay rate are only needed at one order lower in the perturbative expansion. In case of the NNLO t → W + (l + ν)b decay the T 1 > T cut piece can therefore be computed with standard (numerical) NLO technology. The soft function S hl we calculate in the present paper equals the soft function in the 1-jettiness factorization theorem for any heavy-to-light quark decay.
The new three-loop contribution is thus a necessary ingredient for future N 3 LO calculations of differential decay rates based on the N -jettiness method, not only for semileptonic top, but also any other heavy-to-light quark decay.
The outline of this paper is as follows. In sec. 2 we give four (slightly) different definitions of the heavy-to-light soft function S hl and show their equivalence. We give details on our threeloop calculation based on one of these definitions in sec. 3. In sec. 4 we present our results for the renormalized soft function and its anomalous dimension. We also use the latter to check the universal infrared structure of QCD scattering amplitudes that have a massive quark leg. We briefly summarize our findings in sec. 5.

Definitions
The 1-jettiness soft function for heavy-to-light decays is defined by the vacuum matrix element with the soft momentum operatorp µ and the Wilson lines T a is the (ultra)soft SCET (I) gluon field, v µ is the heavy quark velocity (v 2 = 1), n µ is the light-like jet direction (n 2 = 0) and P (P) denotes (anti-)path ordering of the A µ including their SU (N c ) color generators T a . The trace in eq. (2.1) is over color indices, and T[. . .] and T[. . .] represent time-and anti-time-ordered products of the field operators A a µ (x), respectively. The argument ω of S hl can be regarded as the (appropriately normalized) soft contribution to the 1-jettiness observable T 1 , cf. ref. [21].
The soft function in eq. (2.1) equals the perturbative contribution to the shape function in eq. (1.2): 3 where averaging over color and spin of the external HQET b-quark states is understood, the latter are normalized such that h v is the HQET heavy quark field with velocity v, and D µ = ∂ µ + igA µ (0). This b-quark matrix element was calculated to O(α s ) in refs. [23,24] and O(α 2 s ) in ref. [14]. The equivalence to eq. (2.1) can be seen as follows: In eq. (2.6) we could introduce the T and T symbols, because the field operators in Y − are already anti-time-ordered by default as a consequence of the anti-path-ordering, and Hermitian conjugation reverses the order. Similarly the T and T symbols in eqs. (2.7) and (2.1) are in fact redundant, but kept for clarity. In eq. (2.7) we performed the HQET field redefinition where the new (sterile) field h (0) v does not interact with soft gluons anymore, see e.g. ref. [5]. Note that given the (anti)-time-ordering the external b-quark states should be interpreted as 'in' states, The LSZ reduction formula relates the asymptotic 'in'/'out' b-quark states of the S-matrix element to weighted integrals over interpolating h v (x)/h v (x) field operators acting on the vacuum at macroscopically large negative/positive times. Performing the field redefinition in eq. (2.9) one also has to take into account factors of X + from these interpolating fields. For the 'in' states this factor is trivial and we can effectively replace [25] b where i, j are color indices in the fundamental representation. Here and in the following we assume the b quark to be at rest, i.e. v µ = (1, 0) for simplicity. The spatial position x of the endpoint of the Wilson line X + in eq. (2.10) is fixed by the position of the field operator, here h v (0), acting on the asymptotic state, because . Finally, the sterile HQET quark field operators in eq. (2.7) annihilate the sterile external quarks and the color averaging implicit in the b-quark matrix elements translates to 1/N c times the color trace in eq. (2.1).
For the actual calculation of the soft function it is convenient to express it as the imaginary part (discontinuity) of a (1 → 1) 'forward scattering' matrix element. Starting from eq. (2.6) and inserting a complete set of states we have (2.14) Here we use the usual light-cone (Sudakov) decomposition of four vectors: a µ = a − n µ /2 + a +nµ /2 + a µ ⊥ with n 2 =n 2 = 0 andn · n = 2. Note that the momentum operatorp + acting to the left on the external HQET state in eq. (2.12) vanishes, because the external heavy quarks are onshell and therefore have zero residual (soft) four-momentum. While for spacelike distance x 2 < 0 the field operators commute, the theta function θ(x − ) in eq. (2.13) implies t > 0 for time-like distances x 2 > 0. After combining the two Wilson lines using their unitarity property the remaining field operators are therefore automatically time-ordered. To make this manifest we explicitly inserted the T symbol in eq. (2.14).
We can now again perform the field redefinition in eq. (2.9). This time, however, the time ordering in eq. (2.14) implies that |b v = |b v , in and b v | = b v , out|. In contrast to eq. (2.10) the field redefinition now induces a non-trivial factor from the interpolating fields generating the 'out' state [25]: Using r = x − /2 as integration variable we thus obtain  Wilson line propagator with soft light-cone momentum ω: To conclude this section we comment on the relation of S hl to the analogous 1-jettiness soft functions where one or both Wilson lines in eq. (2.1) are changed from incoming to outgoing or vice versa. In the underlying full QCD processes the external heavy and light quark lines are correspondingly crossed from initial to final state or vice versa. Some of these soft functions are e.g. relevant for s-and t-channel single top production as well as charm production in deep-inelastic neutrino scattering ('light-to-heavy DIS'). For state-of-the-art fully-differential NNLO predictions we refer to ref. [26], refs. [27,28], and ref. [29], respectively. The soft function for the light-to-heavy DIS process is for instance simply given by interchanging X ↔ Y (i.e. v µ ↔ n µ in the Wilson lines) in eq. (2.1). Up to two loops the soft functions for the crossed processes can be shown to equal S hl as defined in eq. (2.1) in analogy to the massless case [30]. Unfortunately there is, to the best of our knowledge, no simple argument why this equality should hold at three loops and beyond, not even between heavy-to-light decay and light-toheavy DIS soft functions. 4 A dedicated three-loop analysis along the lines of ref. [30] would require to derive the analytic structure for the relevant two-loop single-emission and one-loop double-emission heavy-light soft currents, which is beyond the scope of this work.

Calculation
Our three-loop calculation of the soft function S hl is based on the definition in eq. (2.16) and performed very much along the lines of our jet function calculation in ref. [10], to which we refer for more details. We work in general covariant gauge with gauge parameter ξ, where ξ = 0 corresponds to Feynman gauge. Ultraviolet (UV) and (intermediate) IR divergences are regulated with dimensional regularization (d = 4 − 2 ).
We use qgraf [32] to generate all relevant three-loop (propagator-type) Feynman graphs with one internal light-like and two external time-like Wilson lines (one incoming, one outgoing), like the ones in fig. 2. 5 The diagrams are further processed by an in-house Mathematica code which assigns the corresponding Feynman rules and performs the necessary Dirac, Lorentz and color algebra. After that the diagrams are given by linear combinations of scalar Feynman integrals. These integrals can then be mapped onto 16 integral topologies with twelve linearly independent linear and quadratic propagators. The associated 16 integral families contain integrals with integer propagator powers ranging from minus three to plus five. The mapping of Feynman integrals onto the different topologies requires numerous multivariate partial fraction operations on products of linear Wilson line propagators followed by suitable shifts of the loop momenta. In order to automatize the extensive partial fractioning we implemented the algorithm outlined in ref. [33] in our code.
Next, we perform the integration-by-parts (IBP) reduction [34] of the integrals in each of the 16 families to a set of master integrals (MIs) using the public program FIRE5 [35]. 6 We then identify pairs of equal MIs of different families that are related by shifts of their loop momenta. The resulting total set of MIs across the 16 families still turns out to be redundant. We find 14 additional relations involving at least three MIs of different families due to partial fraction identities among their linear propagators. Finally, the three-loop contribution to the matrix element in eq. (2.16) can be expressed as a linear combination of 45 MIs belonging to nine different integral families. 7 At this point we already notice that the gauge parameter ξ manifestly cancels out in the sum of all diagrams indicating the correctness of our setup. The 45 contributing MIs can be cast into the form with the following (linearly-dependent) propagator denominators where the usual (causal) '−i0' prescription, i.e. D i → D i − i0, is understood. The nine integral families containing the 45 MIs are defined by their maximal topologies with twelve linearly independent D i . These topologies are determined by restricting the propagator powers in 5 Although the relation of S hl to the time-ordered product of Wilson lines in eq. (2.16) was not made explicit, also the NNLO computation of S hl in ref. [14] was performed in terms of the same type of loop diagrams. 6 The plain IBP reduction with FIRE5 yields an overcomplete set of MIs. To obtain a minimal MI basis for each family we employ the algorithm of ref. [33] to identify equal Feynman integrals. This algorithm is implemented in the FindRules command of FIRE5, which we apply to a large list of test integrals in each family. The output are identities among these integrals, which must also hold after IBP reduction. Demanding this yields another eight independent relations between MIs belonging to the same family, see also ref. [10]. 7 The total number of linearly independent MIs across all families is 64, but only 45 contribute to S hl . eq. (3.1), for instance by From the scaling properties of the integrand in eq. (3.1) for general time-like vector v µ and light-like vector n µ we conclude To solve the remaining 33 MIs we proceed in the same way as for our jet function calculation in ref. [10]. The method was inspired by refs. [37][38][39]. The key idea is to express the 33 MIs as a linear combination of quasi-finite integrals and known MIs. Quasi-finite integrals are free of (endpoint) divergences from the integrations in the Feynman parameter representation (at the Euclidean point ω = −1) for some (even) integer dimension. Starting from a given MI in d = 4 − 2 one can construct a corresponding quasi-finite integral by raising the spacetime dimension by an even number and/or increasing appropriate propagator powers by integer amounts. The former decreases (increases) the degree of IR (UV) divergence, whereas the latter decreases (but not necessarily increases) the degree of UV (IR) divergence. To systematically identify suitable quasi-finite integrals we employ the dedicated algorithm implemented in the public program Reduze2 [40]. For our purposes we find 18 integrals that are quasi-finite in 4 − 2 and 15 integrals that are quasi-finite in 6 − 2 dimensions. To compute them in the respective dimension we first expand their nonsingular integrands in the Feynman parameter representation to high enough order in . We then perform the integrations with the help of HyperInt [41], a powerful computer algebra package for the analytical evaluation of convergent linearly reducible (Feynman) integrals in terms of multiple polylogarithms. The quasi-finite integrals (in their respective dimension) are related to the original MIs (in d = 4 − 2 ) by dimensional recurrence [42][43][44] and IBP reduction. To determine the relevant dimensional recurrence relations between integrals in d and d + 2 dimensions we use the public code LiteRed [45,46]. Our choice of the 33 quasi-finite integrals is such that their results together with the 12 already computed MIs uniquely determine the remaining 33 MIs. We successfully verified all analytic expressions for the MIs obtained in this way numerically using the sector decomposition program FIESTA4 [47]. Finally we insert the results for the 45 MIs in the IBP reduced expression for each three-loop Feynman diagram contributing to S hl and expand to the required order in , see below. We also repeated the calculation for the relevant lower-order graphs using the same setup.

Results
After computing the relevant Feynman diagrams as described in the previous section we take their imaginary part according to eq. (2.16) using (4.1) Adding the contributions of all diagrams (including the lower-order ones) we obtain the bare soft function 8 in terms of the bare coupling α bare s = Z α µ 2 α s with n f being the number of light (massless) quark flavours. The color constants of the SU (N c ) gauge group are C A = N c , C F = (N 2 c − 1)/(2N c ), and T F = 1/2. The coefficients K X of each color structure are given in app. B. For illustration we show in fig. 2 for each of the six three-loop K X coefficients one sample Feynman diagram (arranged in the corresponding order) that contributes to it. Throughout this work we employ the MS renormalization scheme. The relevant terms of the strong coupling renormalization factor Z α are For the expansion of eq. (4.2) we employ the distributional identity with the usual plus distributions defined as The bare and renormalized soft functions are related by where the ⊗ symbol denotes a convolution of the type Convolutions among the plus distributions L n take the form A generic expression for V mn k can be found in ref. [8].

Anomalous Dimension
The RGE of our 1-jettiness soft function reads with the anomalous dimension For the loop expansion of the anomalous dimensions we adopt the notation (4.13) With the soft renormalization factor Z S determined from our bare results in eq. (4.11) we obtain γ S 0 = 4C F , (4.14) in addition to the known terms of the cusp anomalous dimension given in app. A. The oneand two-loop results in eqs. (4.14) and (4.15) agree with those in ref. [14] (after adapting to their conventions). In the following we relate the soft anomalous dimension to corresponding collinear and hard anomalous dimensions in SCET factorization in order to verify our results.
As we will see γ S 2 can thus also be determined indirectly, i.e. without a dedicated three-loop calculation, using known results. However, to the best of our knowledge, the explicit expression in eq. (4. 16) has not been given in the literature so far.
The anomalous dimension associated with the virtual IR singularities due to strong interactions among onshell partons in a squared QCD scattering amplitude can be understood as the anomalous dimension of a corresponding hard function in SCET. It is therefore intrinsically tied to the UV divergences of soft and collinear operator matrix elements in SCET by RG consistency. The generic all-order structure of the anomalous dimension for QCD amplitudes involving massive quarks was derived in ref. [48]. For the heavy-to-light decay with one massless and one massive external quark the it is given by 9 where p is the (outgoing) four-momentum of the massless quark (p 2 = 0, v · p = m b /2) and Γ q cusp (α s ) is the light-like cusp anomalous dimension in the fundamental representation of SU (N c ). The noncusp anomalous dimensions γ q and γ Q are associated with each massless and massive external quark, respectively. They accordingly contribute to the anomalous dimension of any QCD scattering amplitude with multiple quark legs [48] and are in that sense universal. The renormalized SCET hard function H hl corresponds to the finite part of the respective QCD amplitude squared, i.e. where all IR and UV divergences have been subtracted. We thus have For the QCD amplitude with two external heavy quarks (one outgoing, one incoming, where p 2 i = m 2 i , p 1 · p 2 > 0) the (hard) anomalous dimension reads Here the angle-dependent cusp anomalous dimension Γ Q cusp (β, α s ) with (Minkowskian) cusp angle β = arccosh( p 1 ·p 2 m 1 m 2 ) is defined such that in the large angle expansion, 10 there is no O(β 0 ) term. As the large angle limit corresponds to the limit where the mass of one or both of the quarks vanishes it is not surprising that the coefficient of the leading term in eq. (4.20) equals the light-like cusp anomalous dimension [49,51]. For completeness and comparison we also recall the corresponding anomalous dimension for a QCD amplitude with two massless quarks (p 2 i = 0, p 1 · p 2 > 0): (4.21) 9 Here and in the following we suppress a +i0 accompanying the scalar product in the argument of the logarithm, which is necessary for the analytic continuation to other kinematical situations, where e.g. both quarks are outgoing/incoming. 10 Note that in the literature traditionally often the full Γ hh is referred to as the angle-dependent cusp anomalous dimension, see e.g. refs. [49,50].
We stress that Γ q cusp (α s ) and γ q are the same as in eq. (4.17). Renormalization group invariance of the decay rate in eq. (1.1) requires For the noncusp anomalous dimensions this implies 11 The three-loop contribution γ q 2 was obtained from the calculation of the three-loop massless quark form factor [52] via eq. (4.21). In ref. [10] we directly computed the massless quark jet function anomalous dimension γ Jq 2 . It was initially derived indirectly from the RG invariance of the factorized DIS cross section in the threshold region [53] using the three-loop results of refs. [52,54]. The heavy quark noncusp anomalous dimension γ Q 2 can be extracted from the three-loop result of Γ hh in ref. [50] using eqs. (4.19) and (4.20). In fact it can be read off directly from the (nonlogarithmic) constant in the large-angle expansion of (−Γ hh ) explicitly performed in appendix B of ref. [55]. We have We give the explicit expressions for γ q 2 and γ Jq 2 in app. A. We can now solve eq. (4.23) for γ S 2 and find exact agreement with eq. (4.16). This serves as a valuable cross check of our three-loop calculation of S hl . At the same time it confirms the prediction [48] regarding the two-parton correlation part of the IR singularity structure of QCD scattering amplitudes with massive external quarks according to eqs. (4.17) and (4.19).

Renormalized results
Upon MS renormalization the coefficients in the loop expansion of the 1-jettiness soft function for heavy-to-light quark decays take the form By iteratively solving the RGE in eq. (4.10) as an expansion in α s the terms depending on the renormalization scale µ, i.e. the coefficients S (m) n with n ≥ 0, are completely determined by the lower-order constants S (l<m) −1 and anomalous dimension coefficients. To three-loop order 11 In our convention the jet function RGE is analogous to eq. (4.10).
we have Our explicit calculation of S hl (ω, µ) perfectly reproduces eqs. (4.27) -(4.38), which serves as a cross check. In addition it yields the delta function coefficients The expression for S −1 is new and represents together with the three-loop soft anomalous dimension in eq. (4.16) the main result of this work.

Summary
In this paper we calculated the 1-jettiness (T 1 ) soft function for heavy-to-light quark decays at N 3 LO. The renormalized result is given in sec. 4.2. The three-loop delta-function coefficient in eq. (4.40) and the three-loop contribution to the soft noncusp anomalous dimension in eq. (4.16) represent the genuinely new information at this order. In app. A we also collect all other noncusp anomalous dimensions required for N 3 LL resummed heavy-to-light decay rates that are differential in either T 1 or closely related observables like the photon energy in B → X s γ or the jet invariant mass in B → X u ν. We explicitly checked the relation between hard, soft, and jet anomalous dimensions required by RG consistency. This also confirms the predicted universal structure [48] of the IR singularities of QCD amplitudes due to two-parton interactions involving massive external quarks at three loops. That is because we used this prediction to derive the three-loop hard anomalous dimension for heavy-to-light decays from the known three-loop IR singularities of the massive (heavy-heavy) and massless (light-light) quark form factors.
For N 3 LL accuracy also the three-loop contributions to the hard, jet, and soft functions in the corresponding T 1 -type factorization theorems for the decay rates are needed. Our new soft function result represents together with the three-loop contribution to the jet function, which we computed in ref. [10], the two universal (i.e. process-independent) ingredients at this order. As such they also play a crucial role in the calculation of differential N 3 LO heavy-to-light quark decay rates using the N -jettiness IR subtraction (slicing) method, e.g. for t → W + (l + ν)b. The three-loop calculations of the corresponding process-dependent heavy-to-light hard functions are presumably feasible using state-of-the-art multi-loop technology and may be performed in the not too far future.

A Hard, collinear, and cusp anomalous dimensions
For completeness we collect here the explicit expressions for all anomalous dimensions other than γ S in eqs. (4.14) -(4.16) relevant for the heavy-to-light quark decay up to three-loop order. The convention for the loop expansion of the listed anomalous dimensions is analogous to eq. (4.13).
The hard noncusp anomalous dimension associated with massless external partons appearing in eq. (4.17) is given up to three loops by [53,63]  The hard noncusp anomalous dimension associated with the massive external quarks in eqs. (4.17) and (4.19) has the coefficients The one-and two-loop terms in eqs. (A.7) and (A.8) can be found in ref. [48]. The three-loop contribution is copied for completeness from eq. (4.24). The known terms of the noncusp quark jet function anomalous dimension are [17,53]

B Bare data
Here we present our expressions for the coefficients of the different color structures in the bare soft function, eq. (4.2). We show the results as an expansion in = (4 − d)/2 to the order required for the calculation of the renormalized three-loop soft function using eqs. (4.3), (4.5), and (4.7):