Boosted top production: factorization and resummation for single-particle inclusive distributions

We study single-particle inclusive (1PI) distributions in top-quark pair production at hadron colliders, working in the highly boosted regime where the top-quark p_T is much larger than its mass. In particular, we derive a novel factorization formula valid in the small-mass and soft limits of the differential partonic cross section. This provides a framework for the simultaneous resummation of soft gluon corrections and small-mass logarithms, and also an efficient means of obtaining higher-order corrections to the differential cross section in this limit. The result involves five distinct one-scale functions, three of which arise through the subfactorization of soft real radiation in the small-mass limit. We list the NNLO corrections to each of these functions, building on results in the literature by performing a new calculation of a soft function involving four light-like Wilson lines to this order. We thus obtain a nearly complete description of the small-mass limit of the differential partonic cross section at NNLO near threshold, missing only terms involving closed top-quark loops in the virtual corrections.


Introduction
Nowadays, top-quark production is of great interest in elementary particle phenomenology at hadron colliders. This is due to the fact that top-quark physics is closely connected to the study of the recently discovered Higgs boson [1,2] and to the search for new particles. Millions of top-quark pair events have already been produced at the Large Hadron Collider (LHC). For this reason, the ATLAS and CMS collaborations were able to measure the topquark pair production cross section with remarkable precision, e.g. [3][4][5][6][7][8]. On the theoretical side, precise measurements require calculations of the measured observables which include corrections beyond the next-to-leading-order (NLO) in QCD. As an example, the total cross section, which can be measured with a relative error of approximately 5%, was recently evaluated at next-to-next-to-leading order (NNLO) in perturbation theory [9].
Differential distributions, such as the pair invariant mass distribution, the top-quark rapidity distribution, and the distributions with respect to the transverse momentum (both of the individual top quark or of the tt system) are also of great interest, especially in the search for new physics. The ATLAS and CMS collaborations already measured several differential distributions [10,11]. To date, the full set of NNLO QCD corrections to these observables is not known. However, studies of the soft gluon emission corrections to the toppair invariant mass distribution up to next-to-next-to-leading-logarithmic (NNLL) accuracy were presented in [12,13]. In those works, the resummation of the soft corrections was carried out in momentum space by employing methods developed in [14][15][16]. In the same -1 -

JHEP01(2014)028
papers, approximate formulas including all of the terms proportional to logarithmic (plus distribution) corrections up to NNLO were obtained starting from the NNLL resummation formulas. A study of the top-quark transverse momentum and rapidity distributions within the same approach was carried out in [17]. 1 The NNLL resummation of the transverse momentum distribution of the tt system, which presents additional technical complications with respect to the two distributions mentioned above, was considered in [20,21].
A kinematic situation of special interest for new physics searches is the so-called boosted regime, where the top quarks are produced with energies much larger than their mass. Examples of boosted top production include the differential distribution at high values of pair invariant mass M , or the high-p T tail of the top-quark transverse momentum distribution. The presence of new heavy particles decaying into pairs of energetic top quarks could generate bumps or more subtle distortions of differential distributions in this kinematic region. The LHC at center-of-mass energies of 7 TeV and 8 TeV has started to explore boosted top production experimentally, and more data will become available with the future 14 TeV run. At the same time, highly-boosted production is characterized by energy scales much larger than the top-quark mass, and QCD calculations of the differential cross sections must take this into account.
A factorization formalism appropriate for describing QCD corrections to the pair invariant mass distribution in the limit M ≫ m t was put forth in [22], opening up the opportunity to resum simultaneously soft-gluon corrections and small-mass logarithms of the form ln(m t /M ). This same formalism can be used as a way of simplifying the calculation of higher-order corrections in the small-mass limit, a fact exploited in [23] to obtain an NNLO soft plus virtual approximation to the invariant mass distribution in this limit. The goal of the present work is to develop the framework necessary for describing the highlyboosted limit of single-particle inclusive (1PI) distributions, for instance the p T ≫ m t region of the top-quark transverse momentum distribution. To this end, we derive a factorization formalism appropriate for describing the double soft and small-mass limit of the differential partonic cross section.
Our results have some common ground with those for the pair invariant mass distribution in [22], a fact which is particularly clear when deriving results in the double soft and small-mass limit using those in the soft limit as a starting point. In the soft limit, the partonic cross section for top-quark pair production factorizes into a hard function, related to virtual corrections, and a soft function, related to real emission in the soft limit [13,[17][18][19][24][25][26][27][28][29][30]. The hard function is common to both cases, while the soft function depends on the observable. The small-mass factorization of the virtual corrections to 1PI observables can thus be taken directly from [22]. It involves the virtual corrections calculated with m t = 0, in the form of a "massless hard function", and a second function encoding all m t dependence and related to collinear divergences in the small-mass limit.
On the other hand, the small-mass factorization of soft real radiation for 1PI observables in top-quark pair production has not yet been discussed in the literature, 2 and turns 1 Approximate NNLO formulas for the same observables obtained by means of standard Mellin space resummation methods can be found in [18,19]. 2 Note, however, that the single-hadron inclusive cross section at large values of the transverse momentum -2 -

JHEP01(2014)028
out to be rather different than that for the pair invariant mass distribution. A main result of our paper is that such real radiation factorizes into three component functions, as shown in (2.15) below. The physical interpretation is that soft radiation collinear to the observed top quark, soft radiation collinear to the unobserved anti-top quark, and wide angle soft emission are decoherent and factorize. We make a technical distinction between these different kinds of soft radiation in two ways: diagrammatically, through the method of regions, and at the operator level, in terms of Wilson loops. Our final results associate i ) wide-angle soft emission with a Wilson loop built out of four light-like Wilson lines and involving a delta-function constraint particular to 1PI observables; ii ) soft radiation collinear to the top quark with the Wilson loop defining the soft part of the heavy-quark fragmentation function [32] (this is equivalent to the partonic shape-function from B meson decays); and iii ) soft radiation collinear with the anti-top quark with the Wilson loop defining the heavy-quark jet function introduced in [33]. While the "massless" soft function involving four light-like Wilson lines is a matrix in color space, the two types of soft-collinear objects are color diagonal.
With this factorization at hand, one can resum soft and small-mass logarithms at the level of the differential partonic cross section by deriving and solving renormalization-group (RG) equations for the five component functions, or else use it as a tool for simplifying the calculation of higher-order corrections in this limit. In fact, of the five component functions mentioned above, only the massless soft function has not yet been calculated to NNLO; we build on the literature by performing this computation here. We thus achieve a nearly complete NNLO "soft plus virtual" approximation to the differential partonic cross section in the small-mass limit. The final missing piece is the NNLO virtual corrections involving closed heavy-quark loops and proportional to powers of n h = 1 for the top quark. We leave an analysis of these corrections to future work, emphasizing their potential complications on the factorization formalism in the small-mass limit. Even in their absence, our results represent the most complete fixed-order calculation of the p T distribution for boosted production performed so far. They go beyond the approximate NNLO formulas derived in [17,18] by determining the non-logarithmic (delta-function) coefficient in addition to the logarithmic plus distribution terms. They are also consistent with them, and the fact that the NNLO logarithmic plus distribution contributions obtained with the two methods are identical in the small-mass limit is a strong check on our factorization formalism.
The paper is organized as follows. In section 2 we introduce some notation and provide the factorization formula for the partonic cross section in the double small-mass and soft limit. We then devote section 3 to details of factorizing soft real radiation in the smallmass limit. In section 3.1 we analyze NLO phase space integrals in this limit using the method of regions, identifying three distinct momentum configurations which appear at leading power. In section 3.2 we discuss the all-order factorization of the soft function into a convolution of component parts related to these three momentum regions. In section 4 we give expressions needed for fixed-order expansions and present RG equations needed to resum logarithmic corrections. We discuss subtleties related to closed top-quark loops in of the produced hadron was recently studied in [31]. section 5, and conclude in section 6. Some details of the NNLO calculation of the massless soft function for 1PI kinematics are given in appendix A, while explicit expressions for the anomalous dimensions and matching coefficients are collected in appendix B.

Kinematics and factorization
We consider the scattering process where N 1 and N 2 indicate the incoming protons (at the LHC) or proton and anti-proton (at the Tevatron), while X represents an inclusive hadronic final state. In the Born approximation and also to leading order in the soft limit we will deal with later on, two different production channels contribute to the partonic scattering process (2.1): the quark-antiquark annihilation and gluon fusion channels. The partonic processes which we will analyze in detail are thus whereX contains any number of emitted partons. The relations between the hadronic momenta (P i ) and the momenta of the incoming partons (p i ) are p 1 = x 1 P 1 and p 2 = x 2 P 2 . At the hadronic level, we define the Mandelstam variables as while the corresponding quantities at the partonic level are given bŷ It will also be useful to introduce the variable Momentum conservation implies that s 4 = 0 at Born level (k = 0). We will be interested in the double differential distribution with respect to the transverse momentum p T and rapidity y of the top quark in the laboratory frame. Such 1PI observables are obtained by integrating over the phase space of the unobserved anti-top quark, along with any extra real radiation. The p T and rapidity are related to the hadronic invariants (2.3) according to where m ⊥ = p 2 T + m 2 t . Using (2.4) allows one to express the partonic Mandelstam variables in terms of the p T , y, x 1 , x 2 . Then, assuming factorization in QCD 3 and ignoring 3 We should mention however potential subtleties in two-to-two processes pointed out in [34][35][36][37]. power corrections in Λ QCD /m t , one can write the double differential distribution as where the f i/N are universal non-perturbative PDFs for the parton i in the hadron N and the hard-scattering kernels C ij are related to the partonic cross section and can be calculated perturbatively as series in the strong coupling constant. In addition, the lower limits of integration are given by The hard-scattering kernel is a function of the kinematic invariants needed to describe the differential cross section. As long as these invariants are parametrically of the same order, an expansion of the C ij in fixed orders of the strong coupling constant is appropriate. An interesting situation arises when there is a large hierarchy among two or more of the kinematic invariants. In that case it is often possible to factorize the hard-scattering kernel into a product of simpler functions depending only on a single mass scale, up to corrections in the small ratio of disparate scales. This factorization is useful for two reasons. First, the component functions are typically easier to calculate than the full hard-scattering kernels. Second, the factorization formula can be used as a starting point for resumming large logarithmic corrections in the ratio of scales which appear in the higher-order perturbative corrections.
An example often considered in the literature is the soft gluon emission limit, where the partonic invariants satisfy the parametric relation s 4 ≪ m 2 t ,ŝ,t 1 ,û 1 . In this limit the C ij factorize into a matrix product of a hard function H m ij and a soft function S m ij as follows: 4 Such a factorization was first derived in [38]. The hard and soft functions are two-by-two matrices for the qq channel, and three-by-three for the gg channel; the matrix structure is related to the mixing of a basis of color-singlet amplitudes through soft gluon exchange. The hard function is related to virtual corrections, and the soft function is related to real emission in the soft limit. Real emission in the soft limit is considerably easier to calculate than in the generic case. The eikonal factors related to soft gluon emissions exponentiate into Wilson lines, and at the level of the squared matrix element form a gauge invariant Wilson loop operator. Much is known about the perturbative properties of such Wilson loops. In position, Laplace, or Mellin space they contain a series of double logarithmic corrections.

JHEP01(2014)028
In momentum space, these translate into logarithmic plus distribution and delta-function corrections. In particular, defining expansion coefficients of the hard-scattering kernels as the n-th order term contains a tower of logarithmic plus distributions, a delta function term, and regular terms in the s 4 → 0 limit. Consider for instance the NNLO coefficient, which is currently not known. It has the form where the plus distributions Soft-gluon resummation at NNLL can be used to determine the coefficients D i of the plus-distribution contributions [17,18]. The delta-function coefficient, which is formally of NNNLL order, is unknown, as is the term R, which is non-singular in the s 4 → 0 limit and is related to hard gluon emission. In this paper we will discuss the application of the factorization formula (2.9) in the soft limit to the high-p T region of the double differential cross section, where m t ≪ p T . Producing the top quark with high transverse momentum requires that the partonic centerof-mass energy be large, so in this regime the generic situation is that s 4 ≪ m 2 t ≪ŝ,t 1 ,û 1 . We will refer to such a hierarchy of scales as the double soft and small-mass limit. In this limit it is possible to factorize the hard and soft functions themselves. We explain the form of this factorization in the remainder of the section. To simplify the discussion we ignore for the moment contributions from closed top-quark loops appearing in virtual corrections.
The factorization of the hard function in the small-mass limit was derived in [22]. It reads where x t ≡ −t 1 /ŝ and we used momentum conservation −û 1 /ŝ = 1 − x t (valid at Born level and also in the soft limit) in order to simplify dependence on the Mandelstam variables. This factorization can be thought of as a division of the virtual corrections into two -6 -JHEP01(2014)028 momentum regions. The hard matrix appearing on the right-hand side is related to the virtual corrections evaluated with m t = 0 and receives contributions from loop momenta whose virtuality is at the scaleŝ, while the coefficient function C D contains all the collinear singularities appearing in the limit m t → 0 and receives contributions from loop momenta with virtuality at the scale m 2 t . The factorization thus separates physics from the widely separated scales m 2 t ≪ŝ. Two different ways of deriving (2.14) were discussed in [22]. The first relied on the factorization of the heavy-quark fragmentation function in the soft limit [39][40][41][42], and the second used the factorization formula [43] (see also [44]) relating massive amplitudes in the small-mass limit to their massless counterparts.
The factorization of the soft function in the small-mass limit is more subtle. Compared to the factorization of the hard function and even the analogous factorization for the soft function appearing in the top-pair invariant mass distribution [22], a complication here is that the soft function in the small-mass limit is characterized by three distinct momentum scales rather than two. In the next section, we derive the following result: At leading order, each of the above functions is a delta function in its first argument. At higher orders, the three functions on the right-hand side are characterized by logarithmic corrections at the scale shown in their second argument and following from the parametric relation ω i ∼ s 4 . Before moving on, we discuss the interpretation of each of these component functions. First, the massless soft function S ij is related to wide-angle soft real emission corrections to the partonic processes (qq, gg) → QQ, where q and Q are massless distinct quarks. Such emissions are associated with a characteristic mass scale µ s ∼ s 4 / √ŝ . This massless soft function is analogous to that entering the factorization formula for the invariant mass distribution in the m t → 0 limit and calculated to NNLO in [45]. In fact, we will be able to construct results for the S ij to NNLO using calculations from that paper.
Second, the function S D describes soft emissions which are simultaneously collinear to the observed top quark. The characteristic scale for such soft-collinear emissions is µ d ∼ m t s 4 /ŝ. This function has to do with the soft part of the perturbative heavyquark fragmentation function. Field theoretically, it is related to a Wilson-loop operator closing at infinity and containing a finite segment with light-like separation. It is equivalent to the partonic shape function familiar from inclusive B decays, and was calculated to NNLO in [46].
Finally, the function S B describes soft emissions which are simultaneously collinear to the unobserved anti-top quark. The characteristic scale for this function is µ b ∼ s 4 /m t (note that s 4 ≪ m 2 t is important in this context). Our analysis shows that this function is the so-called heavy-quark jet function introduced in [33] and calculated to NNLO in [47].

JHEP01(2014)028
This function is very similar to Wilson-loop operator used in defining S D , the difference being that it contains a finite segment with a time-like separation instead of a light-like one.
By combining (2.14) and (2.15) one arrives at the factorized form of the hard-scattering kernel valid in the double soft and small-mass limit. The only subtlety is the treatment of terms proportional to n h = 1 and related to top-quark loops. For the counting s 4 ≪ m 2 t , such contributions appear only in virtual corrections and modify (2.14). We discuss them in more detail in section 5.
3 Factorizing soft real radiation in the small-mass limit In this section we discuss the factorization of soft real radiation in the small-mass limit. Such a factorization is equivalent to that of the massive soft function (2.15). We make this clear in the preliminary discussion below, introducing some notations and definitions in the process. We then approach the small-mass factorization in two steps. In section 3.1, we perform a diagrammatic factorization at NLO using the method of regions. Then, in section 3.2, we explain how to encode the all-order contributions from these regions in terms of three distinct Wilson-loop operators.
The massive soft function can be defined as where d R = N c in the quark annihilation channel and d R = N 2 c − 1 in the gluon fusion channel, with N c = 3 colors in QCD. The final state X is built of soft gluons in the massive theory (i.e. that relevant for the kinematic limit s 4 ≪ m 2 t ,ŝ,t 1 ,û 1 ), 5 and is a Wilson loop operator built out of soft Wilson lines where v i is the velocity vector associated with parton i. For i = 1, 2 we have v 2 i = 0, while for i = 3, 4 we have v 2 i = 1. We have made use of the basis-independent color-space formalism of [48] in our definitions. This allows us to deal simultaneously with the two cases (q a 1q a 2 , g a 1 g a 2 ) → t a 3t a 4 , where a i is the color index of the parton with velocity v i . Details on how to convert products T i · T j ≡ a T a i T a j to the basis-dependent matrices used in (2.9) can be found, for instance, in [45], and we do not repeat them here. For now we just mention that the amplitude of the scattering process is represented by an abstract color-space vector |M , and the generators T a i act on these vectors according to rules specific to whether i is a quark or gluon, and in the initial or final state. For example, 5 Here and in the remainder of the paper we avoid notational clutter by dropping the hat on the partonic stateX.
where we have used the identification of −∞ and ∞ to bring the definitions of the normal Wilson lines to the convention in the literature. Other important properties of the color generators are that T i · T j = T j · T i for i = j, and that T i · T i = C i , with C i = C F for quarks and antiquarks and C i = C A for gluons. In addition, amplitudes satisfy color conservation, While this is often expressed as a relation between the generators, i T a i = 0, it is important to keep in mind that it holds only when acting on a color-singlet vector, as above.
The soft function takes into account real radiation in the soft limit. We illustrate this by considering the structure of NLO phase-space integrals for single-particle inclusive observables. The three-body phase space for a final state containing the top-quark pair and a gluon with momentum k is We wish to integrate over the unobserved momenta p 4 and k. To do so, we use a technique introduced in [49]. The idea is to shift integration variables to p 4k = p 4 + k and then split the phase space up into two Lorentz invariant pieces: that for the two-to-two production p 1 + p 2 → p 3 + p 4k , and that for a subsequent two-body decay p 4k → p 4 + k. We thus write After a trivial integration the piece on the second line can be written as an integral over the unobserved gluon momentum: The piece on the first line can be arranged into a form appropriate for describing the double differential 1PI observables and is unimportant for what follows.
To evaluate the NLO real emission corrections to the differential cross section one integrates the squared matrix element over the phase space (3.8). The structure of these phase-space integrals simplifies in the soft limit k → 0, in which case s 4 ≪ m 2 t ,ŝ,t 1 ,û 1 . In this limit one can replace the squared matrix element by eikonal factors for a gluon -9 -JHEP01(2014)028 emission from each leg, approximate p 4k ∼ p 4 in the delta-function constraint, and drop any k dependence in the matrix element arising from the shift p 4 → p 4k − k. One must then evaluate integrals of the form 6 We have introduced factors convenient for the MS renormalization scheme, and absorbed them into the integral measure [dk] defined on the second line. The quantity ǫ = (4 − d)/2 is the dimensional regulator. These integrals are exactly those appearing in the NLO corrections to the soft function (3.1), which shows explicitly its connection with real radiation. In fact, the NLO bare soft function is calculated by associating a color factor T i · T j with each integral and summing over possible attachments to the partons i, j. A first step to factorizing the soft function in the small-mass limit is thus to understand the structure of the integrals (3.9). We turn to this problem in the following subsection, using the method of regions as a tool for performing a diagrammatic factorization. We end this section with some comments concerning the arguments of the massive soft function (3.1). The Wilson lines entering its definition depend on the velocity vectors v i , so the object on the left-hand side depends on invariants formed from the velocities and p 4 = m t v 4 . In order to keep contact with our physical picture of the soft function as representing soft real radiation, we express these scalar products in terms of the Mandelstam variables. However, in studying the properties of the integrals it is sometimes useful to keep the structure of the scalar products explicit. For instance, by considering properties of the integrals (3.9) under simultaneous rescalings of the different vectors and s 4 (see, for instance, [50]) one finds their general functional form is (3.10)

NLO phase space integrals and momentum regions
The NLO integrals (3.9) were evaluated for arbitrary m t in [17]. Here we are interested in the asymptotic expansion of those integrals in the small-mass limit, where m 2 t ≪ŝ,t 1 ,û 1 . To leading order in m 2 t /ŝ the results are ǫ , 6 The integrals Iij in (3.9) are connected to the position space integrals I ′ ij in eq. (20) of ref. [17] through relation -10 -

JHEP01(2014)028
In the above equations, we have definedx t = 1 − x t . These explicit results make clear that some of the integrals are characterized by a single mass scale, while some of them depend on more than one mass scale and contain logarithms of m 2 t /ŝ. We will now show how to reproduce these results using the method of regions [51]. This allows us to factorize the multiscale integrals into a sum of simpler, one-scale integrals. While this method was originally developed to construct the asymptotic expansions of loop integrals and is usually discussed in that context, it applies equally well to the phase-space integrals considered here. At the technical level, the reason for this is that integrals such as (3.9) are equivalent to loop integrals, since one can rewrite the delta-function constraint as the discontinuity of propagators (see for example [52]). Rather than actually doing this, one can simply apply the normal procedure for expanding loop diagrams by regions to the phase-space integrals directly. This proceeds as follows. First, one defines a region by associating a specific scaling to the components of the undetermined momentum k in terms of the external expansion parameter (in our case m 2 t /ŝ) . One then expands the integrand as appropriate for the particular momentum region, and integrates over the whole phase space. After finding all of the possible momentum regions which contribute at a given power, one adds their contributions together to obtain the asymptotic expansion of the full integral.
The exact scalings of the regions which contribute to the integrals (3.9) in the smallmass limit are perhaps not obvious at first sight. However, physical intuition suggests three possibilities: wide-angle soft emission, soft emission collinear to the observed top quark, and soft emission collinear to the unobserved anti-top quark. The regions analysis below shows that this is indeed correct, and moreover fixes the momentum scale associated with each of these regions.
To discuss the momentum regions, let us first introduce four light-like vectors n 1 , n 2 , n 3 and n 4 , whose space components are aligned with the momenta p 1 , p 2 , p 3 and p 4k , respectively. For convenience we normalize the vectors to satisfy n 1 · n 2 = n 3 · n 4 = 2. The other scalar products are then fixed to n 1 · n 3 = n 2 · n 4 = 2x t and n 1 · n 4 = n 2 · n 3 = 2x t . Picking two reference vectors n i and n j , we define the light-cone decomposition of an -11 -

JHEP01(2014)028
arbitrary four-vector k as In the following we drop the ij labels when there is no danger of confusion. A judicious choice of the light-cone vectors for a given integral can significantly simplify calculations, as will become evident in the examples below. For the discussion of regions, it is particularly convenient to choose i = 3 and j = 4. The scaling of the momentum p 4k in the limit s 4 ≪ m 2 t ≪ŝ is then given by where λ = m t / √ŝ . The delta function in (3.9) constrains the components of k to satisfy (3.14) We can use the relations to analyze the leading behavior of the propagators 1/(v i · k) and 1/(v j · k). This powercounting exercise is sufficient to identify the three relevant momentum regions: However, not every region contributes to each integral. For example, it is clear from powercounting that the sc (i.e. soft, collinear to the top) region only contributes to integrals involving v 3 , while the sc ′ (i.e. soft, collinear to the anti-top) region only contributes to integrals involving v 4 . In the following, we structure our discussion by analyzing how the three regions contribute to the list of integrals in (3.11).
Wide-angle soft emission. We first discuss the wide-angle soft region. In this region, to leading power in λ, we can approximate 2p 4k · k ≈ √ŝ n 4 · k in the delta function and also for the propagators. The contribution to the integral I m ij from the soft region is then given by the integral

JHEP01(2014)028
Note that the factor of √ŝ in the definition of the wide-angle soft region is a necessary condition for the delta function constraint to be satisfied, and in fact explains why this particular scaling appears.
The above integral is straightforward to evaluate. 7 It is instructive to write the result in the following way: The above result has several important features. First, the mass scale on which it depends is characterized by µ s ∼ s 4 / √ŝ ∼ |k s |. The scaling of k s enforced under the integrand determines the mass scale in the integral. Second, the dependence on the light-cone vectors n i is of the form required by (3.10). Finally, the result is non-zero only if i = j, and if i, j = 4, because otherwise one of the scalar products in (3.19) vanishes and the prefactor is zero in dimensional regularization. We can see this also in intermediate results. An explicit example is the following integral: The equality follows because the integral is scaleless, as one can verify by choosing n 1 and n 4 as basis vectors for the light-cone coordinate decomposition (3.12) and then integrating over the transverse and k s · n 4 components. The reason this happens is that when a gluon connects partons i and j, the more precise definition of wide-angle soft is If, say, n j = n 4 , one must drop its contribution inside the delta-function constraint, in which case the square of the soft momentum vanishes and the integral is scaleless. The total contribution from the wide-angle soft region to the NLO phase-space integrals is obtained by associating a factor of T i · T j with each integral I s ij and summing over legs. The result is proportional to The contributions above are derived from the general integrals after the replacement v i → n i . The time-like vectors are expanded out into light-like ones, which corresponds to calculating real emission corrections with massless partons. We use this fact in the next section to define the "massless" soft function, and calculate it to NNLO in appendix A. The NLO result for the bare soft function is exactly that given in (3.22), showing the direct correspondence between the operator definition and regions calculation. 7 A step-by-step derivation is given in [50].

JHEP01(2014)028
Soft emission collinear to the top quark. We next consider soft emission which is simultaneously collinear to the top quark. We call this region soft-collinear or simply sc. The scaling of soft-collinear momenta is k µ sc ∼ p µ 3 s 4 /ŝ. In contrast to the wide-angle soft region, to expand the integrand in the soft-collinear region we must keep the m t -dependence in the parameterization of p 3 , i.e, p µ 3 = √ŝ , no further specifications are needed. We can now consider contributions from the soft-collinear region, starting with that to I m 13 . Using the scalings in (3.15) and (3.16b) to perform the expansion under the integrand, we find to leading order in λ: where the second equality follows after a straightforward integration. Note that the integral is characterized by the single mass scale µ 2 ∼ k 2 sc , and that the scaling of k + sc is such that the two terms in the delta-function are of the same order. Moreover, the integral contains no information about the velocity v 1 . It is therefore easy to show that I sc 13 = I sc 23 = I sc 34 . The only other contribution from the soft-collinear region is to I m 33 , which is in fact saturated by that region: It is obvious that the soft-collinear region contributes only to integrals involving v 3 , so this completes the analysis. The total contribution from the soft-collinear region to the massive soft function is obtained by associating a factor of T i · T j with each integral I sc ij and summing over legs. The result is proportional to In the second equality we used color conservation (3.5), after which one sees that the contribution is diagonal in color space. Furthermore, the expansion in the soft-collinear region is such that the delta-function constraint has the form δ(s 4 − √ŝ n 4 · k sc ), i.e. the constraint vector n 4 is light-like. We will see in the next section that both of these features are important when identifying the contributions of this region with the soft part of the heavy-quark fragmentation function S D defined in (3.49) below. In fact, one can check that the NLO bare contributions to that function are exactly reproduced by (3.25), which is especially obvious after writing down the NLO integrals using the Feynman rules for Wilson lines and noting the correspondence with the integrands expanded in the softcollinear region.

JHEP01(2014)028
Soft emission collinear to the anti-top quark. Finally, we consider soft emission which is simultaneously collinear to the unobserved anti-top quark. The scaling of such sc ′ momenta is k µ sc ′ ∼ s 4 p µ 4 /m 2 t . To perform an expansion in this region we parametrize p 4 = √ŝ n 4 /2 + n 3 m 2 t /2 √ŝ , and then v 2 The analysis of the contributions from the sc ′ region to the integrals is very similar to that of the sc region. Using the scalings (3.15) and (3.16c), the contribution to I m 14 from this region is: The final equality follows from direct integration using the standard techniques. One sees that I sc ′ 14 = I m 14 . The integrals show familiar features: the particular scaling of k sc ′ (3.16c) ensures that the two terms in the argument of the delta-function scale the same, and the characteristic scale is µ 2 ∼ k 2 sc ′ . Moreover, the integral depends only on quantities related to parton 4. One can show that I sc ′ 24 = I sc ′ 34 = I sc ′ 14 . Furthermore, I sc ′ 44 = I m 44 . The total contribution of this region to the NLO soft function is thus proportional to As was the case with the emissions collinear to the top, the total contribution is color diagonal. The two regions do not, however, give identical contributions. The reason for this is that while after expansion in the sc region the delta-function constraint involves a light-like vector n 4 , in the sc ′ region the delta-function constraint involves the timelike vector v 4 . Instead of being related to the heavy-quark fragmentation function, these contributions are related to a different object, the heavy-quark jet function S B defined in (3.55) below. Here again one can check that the NLO contributions to S B are exactly those arising from (3.27).
Comments. To summarize, we have found three distinct momentum regions: soft, associated with the scale µ s ∼ s 4 / √ŝ ; soft and collinear to p 3 , associated with the scale µ sc ∼ m t s 4 /ŝ; and soft and collinear to p 4 , associated with the scale µ sc ′ ∼ s 4 /m t . Although all of these scales vanish in the limit s 4 → 0, the method of regions provides a technical way of separating out their contributions to the soft function, and identifying the exact mass scale associated with them.
One could use the regions method to prove the factorization formula (2.15) to all orders diagrammatically. It is convenient instead to use effective field theory to reorganize contributions from the different regions into field-theoretical objects encoding their all-order structure, a problem we turn to next. In either case, one might wonder if the three regions identified here are sufficient also at higher orders. Our explicit checks on factorization described below have shown that this is the case at least to NNLO. We have no proof -15 -

JHEP01(2014)028
beyond that, yet also see no physical effect (other than complications from heavy-quark loops we deal with later) that would give rise to other regions. This is an assumption in the "all-order" analysis that follows, and is common to most "proofs" of factorization relying on the regions method, effective field-theory based or not.

All-order factorization in the small-mass limit
Having identified the momentum regions which contribute to the phase-space integrals in the double soft and small-mass limit, we are now in position to explore their all-order structure. There are two possible routes to doing so. The first is to construct an appropriate version of soft-collinear effective theory and apply it to double differential cross sections for 1PI observables using a multistep matching procedure. Many of the steps of such a construction can be taken over from [13], for the soft limit, and from [33], for the boosted limit. A second, more direct route is to start from the definition of the soft function (3.1) for arbitrary m t , factorize the QCD gluon field appearing in the single Wilson-loop operator into a sum of fields whose Fourier components are restricted to certain regions, and then see how the different component fields factorize into operators. We pursue this second method here, and then comment on the alternate derivation at the end of the section.
Our aim is to decompose the operator definition of the massive soft function (3.1) into component operators whose diagrammatic expansions encode the contributions of the three distinct momentum regions. These operators are functions of gluon fields whose Fourier components are restricted to the scalings appropriate for a particular region. We write the decomposition as A a → A a s + A a sc + A a sc ′ . (3.29) By using the following identity for path-ordered exponentials: , (3.30) it is easy to show that We now need to perform a consistent power expansion in λ. This involves the expansion of the gluon fields themselves as well as their momenta. In the soft-collinear effective theory, the gluon fields scale the same as their momenta [53,54]. It is then clear that only the '+' component of the A sc field and its momentum needs to be kept when it interacts with the A s or the A sc ′ field. The same is true for the A s field when it interacts with the A sc ′ field. This is often called "multipole expansion" in the literature [54]. After this expansion, the velocity vectors in the Wilson lines Y v i andS v i on the right side of (3.32) can be replaced by their components along the plus direction, which in our reference system of choice is n 4 , i.e. v i → n 4 in (3.32). We then redefine the fields as is a Wilson line along the direction of n 4 but in the color representation of parton i, and similarly forS n 4 ,i . To show that these field redefinitions are actually the same for each i, we use the identity with (T c adj ) ab = −if cab . It then follows that the field redefinition of, e.g. the soft gluon field is which does not involve the particular color generator at all. After the field redefinitions, the two functions in (3.32) are just usual Wilson lines in terms of the new fields. We work with the redefined fields in the following, and drop the tilde on the two functions. In soft-collinear effective theory, such field redefinitions also remove the interactions among the various fields in the Lagrangian and are therefore JHEP01(2014)028 referred to as the "decoupling transformations" [55]. It is clear from (3.37) that the field redefinitions on the gluon fields made above are equivalent to those in [55].
We have now achieved a decomposition of the original Wilson lines into three separate Wilson lines, each involving gluon fields with a particular scaling. Moreover, the gluon fields with different scalings no longer interact. To factorize the matrix element, we then use that the inclusive state X can be written as a product of states involving s, sc and sc ′ gluons. Moreover, the 2p 4 · p X factor in the delta function can be written as a sum of the contributions from these modes. We use these facts to write (3.38) The above Wilson line operators (defined in (3.40), (3.53), and (3.56) below) arise after a multipole expansion appropriate for the particular momentum region, according to the rules which we explain below. This expansion ensures that the Feynman rules for the Wilson-line attachments in the different regions are such that they produce the correct, homogeneous expansion appropriate for the momentum region the gluon fields are restricted to. The form of this expansion is very similar to what appeared in the regions analysis in the previous subsection. We thus structure our discussion in a similar way, performing a region-by-region analysis which leads to operator definitions of the objects in (3.38).
The wide-angle soft region. Let us first consider the wide-angle soft region. We parametrize the external momenta as described in the previous subsection. The expansion inside the delta-function and Wilson lines then reads Note that for i = 3, 4, we have v µ i ≈ ( √ŝ /2m t )n µ i . However, an important property of Wilson lines is their invariance under rescalings of the reference vector n i → λn i , for an arbitrary number λ, which can be verified immediately from the definition (3.3) after a change of variables. We used this fact to eliminate factors of √ŝ /2m t . The expansion above implies that in the wide-angle soft region we can treat all partons as massless, replacing the time-like vectors v 3 and v 4 with light-like ones n 3 and n 4 .
From the above discussion, we are led to define the massless soft function as where O s (x) = S n 1 S n 2 S n 3 S n 4 (x) .

(3.41)
We calculate the massless soft function to NNLO in appendix A. Our explicit calculations show that integrals involving parton 4 vanish to this order. It seems likely to us that this is also true at higher-orders, but do not pursue a formal proof here.
Soft emission collinear to the top quark. Consider now the soft-collinear region, where p Xsc ∼ s 4 p 3 /ŝ. To set up a power counting we decompose external momenta as in the regions calculation. We can then expand the delta-function constraint the same way: wheren 3 = n 4 . As for the scalar products with the gluon field, for i = 3 we have We again use the invariance of Wilson lines under the scaling n i → λn i . With an appropriate choice of λ, the scalar product becomes irrespective of whether i = 1, 2, 4. It follows that we can replace the product of Wilson lines as where we have used color conservation. To understand the appearance of anti-path ordering in the second line, we note that color conservation i =3 T a i = −T a 3 only applies when acting on the color singlet amplitude directly, as in (3.5), so one must replace, e.g., where we have used that T i and T j commute when i = j. On the other hand, when i = 3 we just have a standard soft-collinear Wilson line for particle 3, and no further expansion is possible. Therefore, we can identify The squared matrix element involving soft collinear structure is then
Since the operator O sc only involves the color generator of parton 3, the right-hand side of the above equation is diagonal in color space. The function S D is then proportional to the unit matrix, namely S D ≡ 1 × S D . To make a connection with the literature, we use (3.4) to write with the Wilson line Y n defined as in (3.4) but with sc fields only. We now use the Fourier representation of the delta function to write whereP is an operator acting on the external states to pick up their momenta. This operator produces a term that cancels the p Xsc dependence in the Fourier exponent, allowing us to perform the sum over states to find 53) which can be shown to coincide with the soft part of the heavy-quark fragmentation function in eq. (50) of ref. [47] after appropriate replacements. In the above formula, the timeordering T and anti-time-ordering T are imposed to guarantee the correct ordering of the fields. 8 Soft emission collinear to the anti-top quark. Finally, we consider the sc ′ region, where p X sc ′ ∼ p 4 s 4 /m 2 t . For the scalar products v i · A a sc ′ , i = 4, we can perform exactly the same arguments as for the sc region. The sc ′ region then involves the operator wheren 4 = n 3 and the Wilson line Y ′ n is defined as in (3.4) but with sc ′ fields only. As for the the delta-function constraint, there is no possible expansion. Therefore, the matrix element for the sc ′ region is 8 See, e.g., appendix C of [16].

JHEP01(2014)028
The difference between S B and S D is the time-like vector in the delta-function constraint, as opposed to a light-like one. We can now go through the steps discussed above for S D to arrive at the result 56) which is consistent with the definition of the heavy-quark jet function in eq. (46) of ref. [47], after making the adaptions necessary to describe a final-state antiquark.
Comments. After inserting the matrix elements (3.40), (3.53), and (3.56) into (3.38) we arrive at the factorization formula (2.15) for the massive soft function in the small-mass limit. We achieved this by studying the factorization of the Wilson-line definition of the massive soft function in this limit. Another option would have been to analyze the differential cross section in soft-collinear effective theory through a multistep matching procedure, similarly to the analysis in [33], where energetic top-pair production in e + e − collisions was studied. In that case, after integrating out virtualities of orderŝ and m t , one is left with two copies of boosted HQET, which interact only through "soft cross talk". In our analysis, the sc-and sc ′ -momenta play the role of the residual momenta for the two copies of boosted HQET, and the soft momenta the role of the soft cross talk. It is then evident that many steps of an effectivetheory analysis could be carried over from [33] and lead to the same final result. We refer the interested reader to that work for the set-up that could be used in such an effectivetheory analysis.

Fixed-order expansions and resummation
The factorization formalism derived in this work can be used in different ways. The first is to view it as a tool for reformulating the calculation of complicated, multiscale higher-order corrections to the coefficient functions C ij in terms of much simpler one-scale calculations, up to corrections to the soft and small-mass limit. In that case, we need only fixed-order expansions of the component functions appearing in the factorization formula. However, in the limit where the mass scales characterizing the component functions are widely separated, for any choice of a common factorization scale µ f the fixed-order expansion of C ij contains large logarithms of scale ratios which can be resummed by deriving and solving RG equations for the component functions. In this section we collect results for the fixedorder expansions of the component functions to NNLO, and then discuss the structure of their RG equations.
It is simplest to discuss higher-order corrections and RG equations in Laplace space, where the distribution-valued functions related to soft real emission become simple functions, and convolutions reduce to multiplication. We define Laplace transforms of the component functions as We can then write the Laplace-transformed hard-scattering kernels where ξ = s 4 /ŝ, 9 as We now discuss the NNLO corrections to the various component functions above. The channel-independent functionss B ,s D , and C D , all related to (soft) collinear emissions, are particularly simple. For these, we can define coefficients as and similarly fors B and C D . Compact results for all of these functions to NNLO can be extracted from the literature, and are gathered in appendix B. For the channel-dependent, matrix-valued massless hard and soft functions, we define perturbative expansion coefficients Here and in the remainder of the section we suppress the subscript indicating the channel dependence of the hard and soft functions, as well as their explicit arguments. While the NNLO hard functions H (2) are unknown, the quantities Tr H (2)s(0) were recently extracted in [23], using NNLO corrections from massless two-to-two scattering obtained in [56][57][58][59][60] along with a subtraction procedure. The rather lengthy expressions can be found in electronic form with the arXiv version of that paper. The massless soft functions are not available in the literature, but we construct results up to NNLO in appendix A.
Here again the results are lengthy, and are included in Mathematica files with the arXiv submission of the present paper.

JHEP01(2014)028
We can make use of these results to form approximations to the Laplacetransformations of the expansion coefficients of the hard-scattering kernels defined in (2.10). To leading order in the soft and small mass limits, we havẽ The above result for the soft and small-mass limit of the NNLO coefficientc (2) is particularly interesting because the exact coefficient in fixed-order perturbation theory is unknown.
Approximations to this coefficient based on soft -gluon resummation to NNLL order for arbitrary m t were derived in [17,18]. To explain how the results given here go beyond those works, we define an explicit expansion of the coefficient function as The NNLO approximations from [17,18] determine the coefficients c (2,n) with n = 1 . . . 4, as exact functions m t . From the viewpoint of a fixed-order expansion, the results for these coefficients in the small-mass limit, determined from (4.6), do not offer an improvement. However, it is a very non-trivial check on the factorization formalism that the coefficients derived above agree with the small-mass limit of those from [17,18], a fact which we have confirmed.
The NNLO approximations from [17,18] do not determine c (2,0) , as this coefficient is formally of NNNLL order. The expansion (4.6) determines it in the small-mass limit (up to corrections involving heavy-quark loops, which we return to in the next section). It will thus be interesting to study the numerical implications of our results for high-p T top production, where corrections to the small-mass limit are negligible and the extra terms calculated here offer a clear improvement on the NNLO approximations from [17,18]. In fact, our results form the basis for a full NNLO soft plus virtual approximation in the small mass limit, meaning that they determine also the delta-function coefficient in (2.11).
The convergence of the fixed-order expansion discussed above can be invalidated when the logarithms of scale ratios are large. In that case, one must resum the logarithms by deriving and solving RG equations. The RG equations for the channel-independent functions read Perturbative results for the anomalous dimensions to order α 2 s using the above definitions (and that of the γ φ in (4.14) below) are listed in appendix B. The RG equations for the matrix valued hard and soft functions have a similar form. We write these as The anomalous dimension for the hard matrix is where A = 2C F γ cusp in the qq channel and A = (N + C F )γ cusp in the gg channel. That for the soft function has the form The coefficients A s and γ s can be derived from the results above, along with the condition that the µ-dependence of the partonic cross section cancels against that in the PDFs to give a µ-independent hadronic cross section. We write the Altarelli-Parisi kernels in the soft limit as (4.14) We then find A s = C F γ cusp in the qq channel and A s = C A γ cusp in the gg channel, and that The term proportional to ln x t (1 − x t ) is needed to cancel the µ-dependence of the PDFs and follows from the derivation given in section 3.2 of [17]. These RG equations can be solved in the standard way, and in fact many of the ingredients can be recycled directly from [22] after appropriate replacements. The perturbative components gathered here form a starting point for an analysis of a simultaneous smallmass and soft-gluon resummation to NNLL. We plan to return to a phenomenological analysis of resummation effects, and a comparison with the fixed-order approximations defined above, in future work.

Heavy-quark loops
A technical subtlety we now address is the treatment of closed heavy-quark loops. The way in which they contribute to the factorization formula is determined by our parametric counting s 4 ≪ m 2 t . Since in that case the heavy-quark mass is much larger than any of the scales characteristic of real radiation, it is not possible for a soft gluon to split into on-shell top quarks. Therefore, heavy-quark loops do not contribute to any of the functions in (2.15) related to soft real emission. Heavy-quark loops decouple from these functions, and the correct prescription is to evaluate them in a theory with five massless flavors.
For the virtual corrections, on the other hand, there is no such decoupling of heavyquark loops, and one must include them explicitly through diagrammatic computations. This leads to a modification of (2.14). In general, we can write The notation is such that the C ij h contains any explicit dependence on n h = 1 and represents the effects of closed top-quark loops, while the second and third factor in the r.h.s. of (5.1) are the same as in (2.14) and are evaluated with n l = 5 light flavors.
A few words about the renormalization schemes are in order. The hard functions are the Wilson coefficients arising when matching from the amplitudes in full QCD to the ones in the effective theory. To obtain the finite hard functions, we need to perform several renormalizations. These include the usual renormalization of the strong coupling constant, the quark masses, the quark and gluon propagators in full QCD, as well as the renormalization of the operators in the effective theory. We renormalize the masses and the propagators in the on-shell scheme, and renormalize the effective operators in the MS scheme. Furthermore, we renormalize the running coupling constant in the MS scheme with five active flavors (in practice by first performing renormalization with six active flavors and then applying a decoupling transformation).
An interesting question is whether the matching coefficient C h in (5.1) can be factorized into one-scale functions, whether it is diagonal in color space, and whether it depends on the channel. In the absence of heavy-quark loops, the m t dependence is contained solely in C D , and is related to regions of loop momenta collinear to p 3 and p 4 . Top-quark loops can introduce an m t dependence even in diagrams not involving p 3 or p 4 , and can otherwise change the regions analysis in such a way that not all m t dependence is related to collinear regions. A full analysis of these higher-order corrections is beyond the scope of the paper. However, we note that the two-loop corrections depending on n h were calculated in the small-mass limit in [61,62], and do not appear to factorize in a simple way, as pointed out in [43]. The same is true of the analogous NNLO corrections calculated for Bhabha scattering in [44].
In any case, the closed heavy-quark loops can be included in fixed-order perturbation theory by calculating the contributions involving powers of n h to the massive hard function in the small-mass limit. The NLO results can be extracted from [13]. To obtain the NNLO results as a matrix would be rather involved. However, we hope to extract the contribution -25 -of the n h -dependent part of NNLO hard function to the coefficient function (2.9) using the the method followed in [23] in future work. This requires two main pieces related to n h -dependent corrections. The first is the UV renormalized NNLO virtual corrections in the small-mass limit. The two-loop contributions were given in [61,62], but the one-loop squared pieces are not readily available in the literature. In addition, one must determine certain color-decomposed one-loop amplitudes to order ǫ 2 . With these building blocks in place, one can calculate the finite remainder of the n h -dependent terms and add them to the results gathered in the previous section to achieve a full soft plus virtual approximation at NNLO in the double soft and small-mass limit.

Conclusions
We have derived a novel factorization formula appropriate for the double soft and smallmass limits of single-particle inclusive cross sections in top-quark pair production at hadron colliders. This formula applies to double differential distributions in the rapidity and p T of the heavy top (or anti-top) quark within this limit. In the absence of closed heavy-quark loops in virtual corrections, we found that the partonic cross section factorizes into five component functions, each depending on a single momentum scale.
Our method was to start from the factorization formula (2.9) for top-quark pair production in the soft limit and then subfactorize the component parts as appropriate for the small-mass limit. For virtual corrections contained in the so-called hard function the method for doing this was already available in the literature. On the other hand, our result (2.15) for the factorization of the soft function, related to real emission in the double soft and small-mass limits, is new. We motivated the result by first performing a diagrammatic factorization of real emission using the method of regions in section 3.1. Our analysis revealed that three types of soft radiation contribute: radiation simultaneously soft and collinear to the observed top-quark, radiation soft and collinear to the unobserved anti-top quark, and wide angle soft emission. In section 3.2 we showed how to factorize the general Wilson-loop operator definition of the soft function into three component Wilsonloop operators related to these regions. We demonstrated explicitly that the two types of collinear operators are diagonal in color space and connected the Wilson loop operators with the soft part of the heavy-quark fragmentation function and the heavy-quark jet function introduced in [33]. The wide-angle soft emission goes into a "massless soft function" defined in (3.40). It involves a non-trivial matrix structure characteristic of soft emissions in two-to-two scattering.
Most of the component functions entering our factorization formula could be extracted to NNLO from results in the literature. We added to this literature by computing the NNLO massless soft function. We showed that the anomalous dimensions appearing in the renormalization-group equation for the NNLO function is consistent with the factorization formalism to this order, providing a strong consistency check on the factorized formula as well as our higher-order perturbative computation. An equivalent check, described in section 4, is that the NNLO logarithmic corrections obtained by expanding out the factorization formula agree with the small-mass limit of those determined by NNLL soft--26 -gluon resummation in [17,18]. Our results provide nearly all elements for the construction of an NNLO soft plus virtual approximation to the differential cross section in the smallmass limit. The missing piece is the NNLO virtual corrections related to closed heavy-quark loops, which we hope to calculate in future work.
We expect that the results obtained here will provide useful insight into the structure of higher-order QCD corrections to 1PI observables in the boosted regime. On the one hand, the NNLO soft plus virtual approximation can be used to study numerically to what extent logarithmic soft gluon corrections dominate over non-logarithmic delta-function terms, thus assessing the reliability of the NNLO approximations [17,18] for boosted production. On the other hand, our factorization formalism provides the starting point for a simultaneous resummation of small-mass and soft logarithms in the partonic cross section to NNLL.

A The massless soft function to NNLO
The calculation of the massless soft function (3.40) proceeds similarly to [45]. The difference in the present case is that the constraint vector is light-like instead of time-like, but since the basic phase-space integrals were calculated for an arbitrary constraint vector, this offers no additional complication. In fact, the end results are slightly simpler, and integrals involving parton 4 vanish.
We first go through the NLO calculation as an example. We obtain the bare NLO soft function through the following sum over legs: Explicit results for the matrices w ij can be found in [45]. The integral I 1 is defined as 10 while the stripped integralĪ 1 can be obtained from (A.2) through the relation

JHEP01(2014)028
The general result for the stripped integral resulting from gluon emissions associated with partons ij is where a ij = 2n i · n j n 4 · n i n 4 · n j . (A.5) (Note that a ij is different from the one defined in [45].) With this definition, we have The stripped integral can easily be expanded in ǫ (and shown to be equivalent to I 12 from [17], which contrary to first appearance does not depend on m t when Laplace transformed with respect to s 4 / √ŝ instead of s 4 /m t ). We observe that the ǫ-expansion of the integral I 1 (a ij ) involves only logarithms of the argument, since the 2 F 1 function does not depend on a ij :Ī The NNLO contributions read (taking account that attachments to parton 4 vanish) (1) 12 Ī 2 (a 12 ) + C AĪ6 (a 12 ) + C AĪ7,1 (a 12 ) + C AĪ7,2 (a 12 ) (A.8)
The bare function has poles in ǫ which must be removed through a renormalization procedure. This is most easily done in Laplace space. It is straightforward to obtain the bare Laplace-transformed function by performing the integral in the definition (4.1).

JHEP01(2014)028
The form of the RG equation (4.11) implies that the renormalized function can be found through the relations where the renormalization factor Z s reads We have defined expansion coefficients and similarly for the other anomalous dimensions. Note that γ s depends on the γ h through (4.15). An explicit expression for γ h (x t , α s ) can be found in appendix B. By evaluating (A.9) using the above expressions, and by renormalizing the bare coupling constant which appears ins bare in the MS scheme by replacing one finds that the renormalized soft function on the left-hand side is indeed finite. This provides a strong check on the validity of the factorization formula used to derive the RG equation, and also on our calculation of the bare NNLO function. Specifically, by expanding everything in powers of α s /4π, at NLO and NNLO one finds The renormalized soft functions can be written in terms of logarithms of arguments x t and 1 − x t . We list here the results for the renormalized NLO soft function in the qq and gg production channels. The specific form of the soft matrix depends on the choice of the color basis. In this paper we employ the s-channel singlet-octet basis (see for example (31) in [23]). For the sake of brevity we set N c = 3 and take into account that the soft function is symmetric. In the following, we indicate the element ij of the matrixs The matrix elements of the NNLO soft matrices have longer analytic expressions. The interested reader can find them in the form of mathematica input files included in the arXiv submission of the present work.

B Matching coefficients and anomalous dimensions
Here we list the matching coefficients and anomalous dimensions appearing in section 4.
We first list results for the coefficients in (4.4) as well the corresponding coefficients in the expansions ofs B and C D . Using [46], we find (for simplicity, here and below we set the number of colors to N c = 3 in the NNLO coefficient) The NNLO coefficient was originally extracted in [42] using the result fors D along with the NNLO result for the heavy-quark fragmentation function calculated in [63], and yields the above equation with δ C = 0 in the last line. It was extracted directly using the relationship between small-mass and massless amplitudes in the appendix of [22], which yields the above result with δ C = 1. The discrepancy between the two methods of extracting the NNLO coefficient remains unresolved. Finally, the Laplace transformed heavy-quark jet function is easily derived from results given in [47]. Explicitly, For γ S (4.9) one finds [41,42,46]  Similarly, the coefficients of the expansion of γ B in (4.8) are [47] γ B 0 = 2C F , The PDF anomalous dimensions are , (B.10) and for the gluon and quark PDFs respectively. In the s-channel singlet-octet basis, the matrix γ h (x t , α s ) appearing in the definition of γ s (x t , α s ) (4.15) is in the quark annihilation channel, while in the gluon fusion channel one finds γ h gg (x t , α s ) = 2 (γ g (α s ) + γ q (α s )) 1 + N c γ cusp (α s ) (ln x t + iπ) The anomalous dimensions γ q and γ g , entering in (B.12), (B.13) and, consequently in γ s in (4.15), are [15,65] γ q 0 = − 3C F , and [16,65] γ g 0 = − where we need only Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.