{'fraction_non_alphanumeric': 0.03586711273632615, 'fraction_numerical': 0.011525936134976825, 'mean_word_length': 4.195288139508026, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Factorization and Resummation of Higgs Boson Differential Distributions in Soft-Collinear Effective Theory

We derive a factorization theorem for the Higgs boson transverse momentum (p_T) and rapidity (Y) distributions at hadron colliders, using the Soft Collinear Effective Theory (SCET), for m_h>>p_T>>\Lambda_{QCD} where m_h denotes the Higgs mass. In addition to the factorization of the various scales involved, the perturbative physics at the p_T scale is further factorized into two collinear impact-parameter Beam Functions (iBFs) and an inverse Soft Function (iSF). These newly defined functions are of a universal nature for the study of differential distributions at hadron colliders. The additional factorization of the p_T-scale physics simplifies the implementation of higher order radiative corrections in \alpha_s(p_T). We derive formulas for factorization in both momentum and impact parameter space and discuss the relationship between them. Large logarithms of the relevant scales in the problem are summed using the renormalization group equations of the effective theories. Power corrections to the factorization theorem in p_T/m_h and \Lambda_{QCD}/p_T can be systematically derived. We perform multiple consistency checks on our factorization theorem including a comparison with known fixed order QCD results. We compare the SCET factorization theorem with the Collins-Soper-Sterman approach to low-p_T resummation.

The Higgs boson is the last undiscovered particle of the Standard Model (SM), and its discovery is a major goal of both the Tevatron and the Large Hadron Collider (LHC) physics programs. If a scalar particle is discovered at either collider, the measurement of its properties will be crucial to determine whether the particle found is the Standard Model Higgs boson, or whether it hints at physics beyond the SM. The theoretical community has devoted significant effort to understanding precisely the production cross section and decay widths of the SM Higgs particle in order to facilitate such studies, as reviewed recently in Ref. [1]. The dominant production mode at both the Tevatron and the LHC is the partonic mechanism gg → H proceeding through a top-quark loop [2][3][4][5][6]. The perturbative QCD corrections are known through next-to-leading order in full QCD [7,8], while the corrections in the m t → ∞ limit are known through next-to-leading order [9] and next-to-next-to-leading order [10][11][12]. Resummation of logarithmically-enhanced threshold corrections to the cross section has been studied [13][14][15][16]. The inclusion of such theoretical calculations is crucial for experimental searches for the Higgs boson, as they increase the predicted cross section in the SM by a factor of two at the LHC and by more than a factor of three at the Tevatron.
The study of differential distributions of the Higgs boson is also needed in experimental analyses. For example, for a SM Higgs in the mass range 130 GeV ≤ m h ≤ 160 GeV, one of the most promising discovery modes is through the partonic process gg → h → W + W − → + ν −ν . Since the final state contains two neutrinos, reconstruction of the Higgs mass peak is not possible. An understanding of the kinematic distributions for both signal and backgrounds is needed in this search channel. The NNLO differential distributions in the m t → ∞ effective theory were obtained in Refs. [17][18][19][20], and detailed studies of the effects of experimental cuts on Higgs boson cross sections have been performed [21][22][23]. However, a large reducible background comes from pp → tt → bW +b W − → + ν −v + jets. Such backgrounds are brought under control with a series of cuts which include a jet veto so that any process involving a jet with high transverse momentum, taken to be roughly p T > 20 GeV [24][25][26], is rejected. Such a cut selects Higgs boson with primarily low transverse momentum, and therefore a proper implementation of such jet vetoes requires a good understanding of the Higgs differential distributions at low p T where resummation of large p T /m h logarithms is necessary. The study of the low-transverse momentum region in hadronic collisions has been under investigation since the early days of QCD [27][28][29]. It has been studied for the Higgs boson following the seminal analysis of Collins, Soper, and Sterman (CSS) [30,31] in several works [32][33][34][35][36].
The purpose of this paper is to derive a factorization theorem for the Higgs transverse momentum p T and rapidity Y distribution, in the regionQ ∼ m h p T Λ QCD , using the Soft Collinear Effective Theory(SCET) [37][38][39]. HereQ and m h denote the partonic center of mass energy and the Higgs mass respectively. Although we focus on Higgs production, our methods and results can be immediately generalized to the differential distributions of any one or more color neutral particles. The factorization theorem we derive has the form (1) where Q is the hadronic center of mass energy, H is the hard Wilson coefficient arising matching QCD onto SCET, and f i/P is the standard parton distribution function (PDF) for taking a parton of species i from the proton. G ij is a perturbative coefficient at the p T scale that has the form βω × I αβ n;g,i ( x 1 x 1 , t + n , b ⊥ , µ T ) I σω n;g,j ( which is a convolution over the collinear functions I αβ n,n;g,i and the Inverse Soft Function (iSF) S −1 . Logarithms of m h /p T are summed by the Renormalization Group (RG) equations in SCET and are encoded in H(x 1 , x 2 , µ Q ; µ T ) which is the hard coefficient evolved from the renormalization scale µ Q ∼ m h down to µ T ∼ p T . The logarithms of Λ QCD /p T are summed via the standard DGLAP evolution of the PDFs and are encoded in the PDFs evaluated at µ T ∼ p T . The factorization formula in Eq.(1) is derived by matching QCD onto a sequence of effective field theories EFT: which is shown graphically in Fig. 1. The first step QCD(n f = 6) → QCD(n f = 5) denotes the usual procedure of integrating out the top quark to get an effective coupling of the Higgs to gluons. The Higgs production mechanism then proceeds via this effective coupling. The hard scaleQ ∼ m h is then integrated out by matching onto SCET p T , which describes the dynamics of soft and collinear modes with transverse momenta of order p T . The factorization theorem in SCET p T takes the form where the collinear functionsB αβ n,n are the Impact-parameter Beam Functions(iBFs). The iBFsB αβ n,n are extensions of the beam functions that appear in [40,41] and reduce to them for b ⊥ = 0 after contraction of the transverse indices α and β. The beam functions of Ref. [41] were shown to have wide applicability to the analysis of observables at the LHC. The iBFs are proton matrix elements evaluated at the scale µ T ∼ p T . The iBFs are matched onto the standard QCD PDFs by performing an OPE in Λ QCD /p T and the logarithms of Λ QCD /p T are summed via the standard DGLAP equations used to evaluate the PDFs at the scale µ T ∼ p T . This is shown schematically in Fig. 1 and gives the final form of the factorization theorem shown in Eqs. (1) and (2) where the collinear functions I αβ n,n;g,i are just the iBF to PDF matching coefficients.
While the factorization and resummation of transverse-momentum distributions has been studied extensively in the QCD literature, and SCET analyses [42,43] have been performed in the past, our analysis contains several interesting differences that we believe are worth further investigating. A summary of the main points of this paper is given below: 1. We derive a factorization theorem for the Higgs transverse momentum and rapidity distributions using effective field theory methods. A clear separation of the dynamics associated with the scalesQ ∼ m h p T Λ QCD into perturbative Wilson coefficients and standard QCD PDFs is achieved. Large logarithms of ratios of the relevant scales are summed using RG equations in the effective theories. Power corrections in p T /m h and Λ QCD /p T can be systematically derived by going to higher orders in the power counting of the effective theories.
2. In addition to the factorization of the scales m h p T Λ QCD , the perturbative physics of the p T scale is further factorized into an iSF S −1 and two distinct collinear functions I αβ n;gi and I αβ n;gi . This additional factorization simplifies the structure of higher order perturbative corrections at the p T scale. They can now be obtained through higher order computations of the simpler perturbative functions S −1 , I αβ n;gi and I αβ n;gi . The QCD (n f = 5) theory obtained after integrating out the top is matched onto SCET p T at the µ ∼ m h scale followed by RG running down to the µ ∼ p T scale summing logarithms in m h /p T in the process. Using the soft-collinear decoupling property of the leading order SCET p T Lagrangian, the cross-section is factorized into a n-collinear iBF(Impact-parameter Beam Function), an-collinear iBF, and an iSF(Inverse Soft Function). The iBFs are then matched onto the standard QCD PDFs at the µ ∼ p T scale and the logarithms of Λ QCD /p T are summed via the DGLAP equations which determine the PDFs at the µ ∼ p T scale.
subtraction terms required to avoid double counting soft emissions that are already contained in the soft function. However, the equivalence of soft zero-bin and soft subtractions allow us to rewrite the factorization in terms of the standard QCD PDFs. 4. We give expressions for the factorization formula in both impact-parameter space (bspace) and momentum space and discuss how they are related. The factorization can be formulated entirely in momentum space; see Eqs.
(1), (53), and (52). However, the matching coefficients I αβ n;gi and I αβ n;gi are obtained by matching the iBFs in b-space onto the standard QCD PDFs. 5. The iBFsB αβ n,n (x, t, b ⊥ , µ) that appear in the SCET p T factorization theorem in Eq.(4) are extensions of the beam functions introduced recently in [40,41]. The additional functional dependence on b ⊥ found here is required to facilitate resummation of the low p T region, and our functions reduce to those studied previously for b ⊥ = 0. 6. The iBFs are similar to transverse-momentum dependent PDFs studied previously in the literature [44][45][46][47], except for an additional dependence on a residual light-cone momentum component. The iBFs are defined in covariant or non-singular gauges and are invariant under non-singular gauge transformations. In singular gauges, like the light cone gauge where the gauge potential is non-vanishing at infinity, the iBFs must be modified [44][45][46][47] with a transverse gauge link to be gauge invariant under the singular gauge transformations. However, the convolutions over perturbative quan-titiesB αβ n ⊗Bn αβ ⊗ S −1 is fully gauge invariant since it is obtained by applying a soft-collinear decoupling transformation to a manifestly gauge invariant operator in SCET p T . Correspondingly, after matching the iBFs onto gauge invariant PDFs, the quantity G ij = I αβ n;g,i ⊗ In ;g,jαβ ⊗ S −1 is also gauge invariant. As a result, even though the iBFs are defined in covariant gauges, the factorization theorem is completely gauge invariant as expected. The additional functional dependence on the residual light-cone momentum avoids singularities that occur when this variable is inclusively integrated over, as it is when considering the usual transverse-momentum dependent PDFs.
7. The iBFs come with perpendicular indices α, β which are contracted between the ncollinear and then-collinear iBFs. This structure has not been discussed in previous SCET analyses. This non-trivial index structure encodes dot products between the momenta of final state gluons emitted in different directions, and is needed to separately calculate contributions from n andn gluons.
8. The Landau poles associated with impact-parameter integrations in the CSS formulation can be avoided in the effective field theory approach employed here. As already mentioned, the factorization theorem can be formulated entirely in momentum space with no reference to any impact-parameter or impact-parameter integrations. The perturbative function G ij and the PDFs in Eq.(1) are naturally evaluated at the renormalization scale µ ∼ p T so that one never encounters the Landau pole associated with the limit µ → 0.
9. The iBFs are not specific to Higgs production, but have universal applicability to the study of differential distributions in gluon-initiated processes at colliders. Analogous functions exist for quark-initiated processes, and can be derived following the approach discussed in this work.
10. In the non-perturbative region p T ∼ Λ QCD , where the effective field theory power counting in Λ QCD /p T breaks down, a phenomenological model can be introduced for the quantity G ij analogous to what is done in the standard approach.
11. We perform several consistency checks on the factorization theorem including a comparison with known fixed order QCD results. We ensure that the renormalization scale dependence properly cancels among the various factored objects in the factorization formula, which is a non-trivial check on the structure of the factorization theorem.
Our paper is organized is follows. In Section II we review the calculation of Higgs production through the gluon-fusion mechanism in perturbative QCD. We derive our factorization formula in Section III. We discuss the factorization formula and compare to other results in the literature in Section IV. In Section V we calculate the various quantities in fixed-order perturbation theory, while Section VI is devoted to resummation of logarithms. We perform a series of consistency checks on our factorization formula including a comparison with fixed-order perturbative QCD in Section VII. Finally, we conclude in Section VIII.

II. HIGGS PRODUCTION IN QCD
We begin by reviewing the gluon-initiated production of a Higgs boson in QCD. The coupling of the Higgs boson to gluons arises primarily from a top-quark loop. For m h < 2m t , we can integrate out the top quark to derive an effective coupling of the Higgs boson to gluons [2,5,7,9,58]. The effective Lagrangian is given by where C GGh , the Wilson coefficient in the MS scheme, is known through O(α 4 s ) [13,[59][60][61][62]. Calculations of the total cross section at higher orders in QCD perturbation theory using this effective Lagrangian have been shown to reproduce the result of the full theory to percentlevel accuracy when m h < 2m t if the effective-theory cross section is normalized by the full top-quark mass-dependent leading-order cross section [8,63]. Additional corrections to the transverse momentum spectrum of O(p 2 T /m 2 t ) are also present in the full theory. These affect the low-p T region at the percent level, and have recently been studied in Refs. [64,65]. An effective field theory approach for studying such corrections at higher orders has recently been developed [66,67]. The scale µ at which the coupling constant is evaluated should be chosen as µ ≈ m t to minimize logarithms that appear in the O(α 3 s ) expression for C GGh . For notational ease we define and express later results using c. The Higgs boson must recoil against at least one parton in order to have non-zero transverse momentum. At leading-order in perturbative QCD, three partonic processes contribute to Higgs production at non-zero p T : gg → hg, q(q)g → hq(q), and qq → hg. We focus here on the dominant process gg → hg. The other partonic production channels are suppressed by smaller PDF luminosities, and can be straightforwardly included in the analysis if desired. The diagrams contributing to this process are shown in Fig. 2. We denote the incoming proton momenta as p 1 , p 2 , the outgoing parton momentum as p 3 , and the outgoing Higgs momentum as p h . We define the light-cone vectors n andn through FIG. 2: The leading order diagrams in QCD that contribute to Higgs production with non-zero p T .
The momenta of the initial hadrons and the Higgs boson are given by where Y = 1/2 ln(n · p h /n · p h ) is the rapidity of the Higgs boson. The partons entering the hard-scattering process have momenta and virtualities given bŷ where x 1,2 are the usual Bjorken momentum fractions. Several distinct kinematic regions contribute to production of a Higgs boson at low transverse momentum. In the first, the Higgs recoils against a n-collinear gluon with the following momentum in light-cone coordinates: This corresponds to a Higgs boson produced at high rapidity and low transverse momentum. A similar region exists with the Higgs recoiling against ann-collinear gluon with scaling Finally, the Higgs boson may recoil against a gluon with an "soft" momentum p 3 ∼ (p T , p T , p T ). This corresponds to the production of a Higgs at central rapidity and low transverse momentum. Since the final-state gluons in the low-momentum region are restricted to be either soft or collinear to the initial partons, large logarithms of the form ln (m h /p T ) appear in the perturbative expansion and must be resummed to all orders. The n-collinear,n-collinear and soft gluon modes will be used to construct the effective theory that facilitates this resummation.
For future use we reproduce here the differential and total cross sections for the gluongluon scattering process. The hadronic cross section can be expressed as a convolution of parton distribution functions and partonic cross sections: whereŝ,t, andû are the usual partonic Mandelstam variables. For production of the Higgs with non-zero p T , the differential partonic cross section is given by [68] The total partonic cross section for gg → h through next-to-leading order in QCD perturbation theory is [7,9] where z = m 2 h /ŝ. This result assumes the scale choice µ 2 =ŝ. The dependence of the partonic cross section on the renormalization and factorization scales can be restored by using the known renormalization group running of the cross section. The result is presented in Ref. [9].

III. EFT FRAMEWORK
We derive a factorization theorem via a sequence of effective theories which factorize the physics associated with the different scales Q ∼ m h p T Λ QCD into calculable perturbative functions and standard QCD PDFs. As we are assuming that the mass of the Higgs is sufficiently small (m h < 2 m t ), we can integrate out the top quark in the matching step QCD(n f = 6) → QCD(n f = 5) to obtain an effective coupling of the Higgs boson to gluons. The cross sections obtained using this effective theory were described in Sec. II. To derive a renormalization group equation allowing resummation of large logarithms ln (m h /p T ) that appear at low transverse momenta, the matching to SCET p T is required. The soft-collinear decoupling property of the leading order SCET p T Lagrangian also leads to a factorization of the soft and collinear sectors, which simplifies calculations of the cross section in the low p T region. Finally, the matching to SCET Λ QCD expresses the cross section in terms of the standard parton distribution functions. We describe in this section the details of each stage in the matching in QCD(n f = 5) → SCET p T → SCET Λ QCD .
A. QCD to SCET p T As already mentioned, the perturbative expansion in QCD for the transverse momentum spectrum of the Higgs contains logarithms of m h /p T . In the low transverse momentum region Λ QCD p T m h , these logarithms become large and must be resummed to all orders in perturbation theory. In the effective theory formulation, this is done by matching QCD onto the effective theory SCET p T , which describes the dynamics of the degrees of freedom recoiling against the Higgs, and solving the RG equations of the effective theory operators. The effective theory SCET p T is formulated in terms of collinear and soft modes with momentum scalings where p n , pn, and p s denote typical momenta for the n-collinear,n-collinear and soft modes respectively. The effective theory has a well defined power counting in the parameter η and has distinct quark and gluon fields for each of these modes. We note that ultrasoft gluons with scaling p us ∼ m h (η 2 , η 2 , η 2 ) would not produce a Higgs boson with sufficient transverse momentum, and therefore do not need to be included as separate modes in the effective theory. The gluon fields A µ n,pn (x), A μ n,pn (x), and A µ s,q (x) destroy n-collinear,n-collinear, and soft gluons respectively. The presence of distinct collinear and soft gluons requires the theory to be invariant under collinear and soft gauge transformations [38,69]. The momenta of the effective theory modes are separated into labelp and residual k parts Derivative operators are similarly separated into label and residual operators so that, for example, a derivative acting on the n-collinear field takes the form such that the label operators act on the label momentum subscripts and the residual derivative operator acts on the residual co-ordinate dependence x µ . We note that such a field with label momenta can be written explicitly as a Fourier transform of a standard quantum field. As an example, a field with no dependence on residual coordinates can be expressed as As already discussed, after integrating out the top quark, the gg → h process is mediated by the effective QCD operator where L mt is given in Eq. (5). In SCET p T , the leading order operator that mediates this process is [15,42] O(ω 1 , where we have suppressed the perpendicular labels on the soft and collinear fields, the sum of which are constrained to be the negative of the Higgs transverse momentum, and we use the standard notation and the big brackets in Eq. (21) denote the fact that any label operators appearing inside do not act outside the brackets. The B field is defined as [70] and the soft Wilson line S n in position space is given by with an analogous expression for Sn. The covariant derivatives are dressed with momentum space Wilson lines so that and D µ n are the usual covariant derivatives in · D n =P n + gn · A n,p , iD ⊥µ n = P µ n⊥ + gA ⊥µ n,p , in · D n = in · ∂ + gn · A n,p . (26) We match the QCD operator onto the the effective SCET p T operator via and determine the Wilson coefficient by computing the gg → h process on the LHS in QCD and on the right side in SCET p T . At tree level the Wilson coefficient is The process gg → gh where the transverse momentum of the Higgs is of order m h η ∼ p T is mediated through an extra soft or collinear emission from the Wilson lines in O(ω 1 , ω 2 ). We note that one can instead implement a two-step matching procedure by first matching onto an intermediate SCET with scaling parameter ξ ∼ p T /m h . The relevant modes in this theory are collinear and ultrasoft gluons with momenta p n ∼ m h (ξ 2 , 1, ξ), pn ∼ m h (1, ξ 2 , ξ), and p us ∼ m h (ξ 2 , ξ 2 ξ 2 ). The leading operator in this theory takes the schematic form where Y n,n are ultrasoft Wilson lines built up from eikonal interactions with the collinear gluons. The ultrasoft gluons of the intermediate theory just become the soft modes of SCET p T , and we replace Y n,n → S n,n to obtain Eq. (21). The Wilson coefficient from matching these operators is unity to all orders as also seen in [42]. A more detailed explanation of this point can also be found in [71].
For convenience, in deriving the factorization theorem, we will work with soft fields that are in position space and collinear fields that are in position space conjugate to the transverse label momentum so that no label momenta of the transverse momentum appear in the soft and collinear fields. This amounts to not separating from momentum components of order m h η a residual part of order m h η 2 . The separation into label and residual components is employed only for the hard light cone collinear momentum components of order m h .

B. Factorization in SCET p T : iBFs and the soft function
It is easier to work in terms of the hadronic Mandelstam variables u and t instead of p T and Y , which correspond to the Higgs transverse momentum and rapidity respectively. These two sets of variables are related by The transformation between these sets of variables has a rather simple Jacobian given by Thus, a restriction on the u and t Mandelstam variables is equivalent to a restriction on the p T and Y of the Higgs. The double differential cross-section in the Mandelstam variables can be written in SCET as where O and C denote the SCET p T operator and the matching coefficient respectively. The ⊗ symbol denotes a convolution in the label momenta ω 1,2 as in Eq. (27). Note that the constraint delta functions δ(u − (p 2 − p h ) 2 ) and δ(t − (p 1 − p h ) 2 ) restrict the final states to those that satisfy u = (p 2 − p h ) 2 and t = (p 1 − p h ) 2 , or equivalently pick out the states with the corresponding values of p T and Y . The states X n , Xn, X s correspond to final state particles with the n-collinear,n-collinear and soft momentum scaling respectively. It is only the states with such momentum scalings that will have a non-zero overlap with the SCET p T operator O(ω 1 , ω 2 ). The overall factor of 1/4 in square brackets in Eq. (31) is from the average over the initial proton spins. Using the fact that the soft and collinear modes are decoupled in the leading order SCET p T Lagrangian we arrive at the factorization formula where the jet and soft functions are defined as T is the time-ordering symbol, andT denotes anti-time ordering. Details of the derivation of this formula are given in appendix A. The above factorization theorem can be brought into a more concise form involving a simpler convolution structure so that We recast this factorization theorem in terms of jet and soft functions that have momentum space light cone coordinates as where we have defined the hybrid-fourier space jet and soft functions as The hybrid jet functions B αβ n,n (ω, k ± , b ⊥ , µ) are similar to the functions that appeared in [40] and more recently in [41], but differ because of their dependence on the impact parameter b ⊥ and because their perpendicular indices α, β are not contracted with each other. We will refer to these functions as impact-parameter Beam Functions (iBFs) in analogy to the Beam Functions of [41] which have b ⊥ = 0. These iBFs are implicitly defined with a zero-bin [72] subtraction in order to avoid double counting the soft region already present in the soft function S(ω 1 ,ω 2 , b ⊥ , µ). For clarity we will refer to the zero-bin subtracted iBF as a purely collinear iBF. We will refer to the iBF defined without a zero-bin subtraction as the naive iBF or simply the iBF 1 when the context is clear. These purely collinear iBFs that appear in the factorization theorem are in general gauge dependent quantities. This is seen from their 1 Note that the iBF is still always implicitly defined with an "ultra-soft" zero-bin subtraction [41] to avoid double-counting between the collinear and the ultra-soft modes with scaling p us ∼ M (η 2 , η 2 , η 2 ). In the paper, the phrase "zero-bin subtraction" refers to the "soft" zero-bin subtraction, to avoid doublecounting between the collinear and soft modes, unless otherwise specified. The ultra-soft and soft zerobin subtractions are distinct since the soft modes are constrained by p T while the ultra-soft modes are unaffected. The soft zero-bin subtraction in the beam functions can be recast in terms of the two iBFs and an inverse soft function (iSF) as shown in the next section. Now the iBF is defined without the soft zero-bin subtraction which is instead taken into account by the presence of the iSF. However, there is still a double-counting left over between the collinear and the ultra-soft modes (which are unaffected by p T constraints) which is avoided by the ultra-soft zero-bin subtraction in the iBF. This ultra-soft zero-bin subtraction in the iBF should be understood to be implicitly part of its definition. Thus, at b ⊥ = 0, the iBF reduces exactly to the beam functions in [41]. In pure dimensional regularization, the ultrasoft zero-bin subtraction is scaleless and vanishes. We thank I.W. Stewart and F. Tackmann for detailed dependence on the impact parameter b ⊥ = 0 which leads to a spatial separation between the fields in the matrix element not connected by a Wilson line. However, this additional gauge link can be placed at infinity along the light-cone, and it does not contribute in covariant gauges where the gauge potential vanishes at infinity. The iBF is thus well-defined in covariant gauges. This is similar to what occurs for transverse-momentum dependent PDFs in QCD [44][45][46][47]. In light-cone gauge, this additional gauge link at infinity is required due to the asymptotic behavior of the gauge potential. It is possible that Glauber modes are responsible for building up this extra contribution in SCET [47]. We note that the total convolution over the hard Wilson coefficient |C(ω 1 , ω 2 , µ)| 2 , the purely collinear iBFs, and the soft function in Eq. (35) is just equal to the total perturbative cross-section for gluon-initiated Higgs + multi-parton production and thus gauge independent as required.

C. Equivalence of zero-bin and soft subtractions
The purely collinear iBFs B αβ n,n (ω, k ± , b ⊥ , µ) defined with a zero-bin subtraction can be written as whereB αβ n,n (ω, k ± , b ⊥ , µ) is the naive iBF or simply the iBF defined without a soft zerobin subtraction but still with an implicit ultrasoft zero-bin subtraction. The functions B αβ {n0,n0} (ω, k ± , b ⊥ , µ) denote the soft zero-bin limit of the iBFs. It has been demonstrated in several processes [73][74][75] that the soft zero-bin subtraction is equivalent to a subtraction of the soft function. The same holds true in this case allowing us to recast the factorization in terms of the iBFs as opposed to the purely collinear iBFs. In appendix B we show that the convolution over the purely collinear iBFs and the soft function can be replaced with a convolution over the naive iBFs with an Inverse Soft Function (iSF) so that It is useful to recast the factorization theorem in the latter form in terms of iBFs and an iSF and eventually the iBFs will be matched onto the standard QCD PDFs via an operator discussions which led to an understanding of this point while at the Aspen Center for Physics summer workshop on Forefront QCD and LHC Discoveries, 2010. We also thank A. Jain, M. Procura, and W. Waalewijn for an inquiry which prompted us to explain this point here explicitly.
product expansion. The factorization theorem in terms of the iBFs and an inverse soft function takes the form For later convenience we switch to the variables so that the factorization theorem takes the form We have definedB and In the above equation, the iBFs are written in terms of t ± n,n instead of k ± n,n to make it manifest that the iBFs actually depend on t ± n,n , as demanded by reparameterization invariance. The choice of the scale µ T ∼ p T will become manifest once we perform the Higgs phase space integrals and rewrite the u and t Mandelstam variables in terms of p T and Y. We will do this in the next section. The RG evolved H(x 1 x 2 Q 2 , µ Q ; µ T ) hard function sums up logarithms of m h /p T . The iBFs are proton matrix elements and will give rise to logarithms of Λ QCD /p T in the perturbative cross-section that must be resummed. For this reason, as discussed in the next section, the iBFs will be matched onto PDFs and the logarithms of Λ QCD /p T will be resummed via the standard DGLAP evolution equations.

D. iBFs to PDFs
The matching of the iBF onto the PDF is given bỹ where I αβ g,i ( z z , t + n , b ⊥ , µ) is the matching coefficient and the gluon pdf is defined as so that the leading order perturbative expression is normalized as A matching equation analogous to Eq. (43) holds for then-collinear iBFB αβ n . Note that in the iBF to PDF matching in Eq. (43), the iBF can match onto quark PDFs beyond tree level. In the initial analysis presented in this paper, we ignore this effect. It is simple and straightforward to include the effects of the quark PDFs if desired. By noting the the PDF is scaleless and that the infrared structure of the iBF and PDF match, as discussed further in Section V, one obtains the following all orders expression for the Wilson coefficient We explicitly check this expression at one order beyond tree level in Section V. An analogous expression holds for then-collinear Wilson coefficient I βᾱ n;g,i . Using Eq. (43) in Eq. (41) we arrive at the factorization theorem Next we perform the Higgs phase space integrals to get where we have defined the u and t dependent function In terms of the p T and Y variables, related to the Mandelstam u and t variables as in Eq. (29), we have where we have defined the p T and Y dependent perturbative function

E. Momentum space vs impact-parameter space
As seen in the last section, the factorization formula for the p T , Y and u, t dis- (49), respectively. These functions are defined in terms of impact-parameter space Wilson coefficients I βα n,n;g,i (z, t, b ⊥ ) and the impactparameter space iSF S −1 (k − , k + , b ⊥ ). We can now introduce momentum space Wilson coefficients and a momentum-space iSF via so that the function G ij (x 1 , x 1 , x 2 , x 2 , p T , Y, µ T ) can be written as with no reference to b ⊥ . The expression for G ij (x 1 , x 1 , x 2 , x 2 , u, t, µ T ) can be trivially obtained from Eq. (53) by rewriting the p T and Y variables in terms of the u and t variables by using Eq. (29). Eq. (50) together with Eq. (53) gives the factorization theorem entirely in momentum space with no reference to the impact parameter. This form of the factorization makes manifest the p T dependence and how it is related to the momenta of the final state particles involved. However, the impact-parameter space formulation is necessary for the matching of the iBF onto the PDF. This matching must be done in impact-parameter space as in Eq. (43). It is only in impact-parameter space that the infrared singular structure of the iBFs and the PDFs becomes manifest. Thus, in order to obtain the momentum space factorization, one must first go to the impact-parameter space formulation and obtain the iBF to PDF matching coefficients I αβ n,n;g,i (z, t, b ⊥ ) and then fourier transform to momentum space.

A. General discussion
In this section we summarize the main points and give a more detailed discussion about the structure of the factorization theorem. We also compare our results with both previous SCET studies of low-p T resummation and with the standard QCD approach. In the remaining sections we present fixed order results for the various perturbative functions appearing in the factorization theorem, anomalous dimension computations and RG equations, various internal consistency checks of the factorization theorem, and checks of the factorization by comparing with known fixed order QCD results.
To recap, the factorization theorem for the Higgs transverse momentum p T and rapidity Y has the form where H and G ij are perturbative coefficients and f i/P are the standard QCD proton PDFs. The indices i, j run over the gluon, quark, and anti-quark species. In our initial analysis we have ignored contributions from quark and anti-quark initiated Higgs production since they are numerically small, but it is straightforward to generalize the analysis and include these effects if desired.
The factorization theorem achieves a clear separation of the physics associated with the scalesQ ∼ m h p T Λ QCD which are now encoded in the hard function (27) and (42), followed by RG running down to the scale µ T ∼ p T . The RG running sums up the logarithms of m h /p T . The perturbative function G ij (x 1 , x 1 , x 2 , x 2 , p T , Y, µ T ) is a convolution of two perturbative collinear functions and an iSF(inverse Soft Function). The result for G ij involves logarithms of p T /µ T and it is thus natural to choose µ T ∼ p T to minimize these logarithms. The QCD PDFs f i/P (x, µ T ) are also evaluated at this scale µ T ∼ p T and the logarithms of Λ QCD /p T are summed via the standard DGLAP evolution of the PDFs. This picture is roughly summarized in Fig. 1. We note that in principle, the iSF could depend on an additional scale µ S . In this analysis we set µ S = µ T , since we see no evidence of large logarithms associated with the additional scale Q that must be resummed via an additional renormalization-group evolution in µ S . We show this when calculating the iSF in Section V. If a higher-order analysis indicates such a need for a separate µ S , it is straightforward to include.
The perturbative function G ij defined in Eq. (51) has an intricate structure on which we now elaborate. As seen in Eq. (51), G ij is a convolution, in light-cone momenta and an impact parameter, over a n-collinear function I αβ n;g,i ,n-collinear function I αβ n;g,i and the iSF S −1 . The I αβ n,n;g,i functions are obtained from matching the iBFs (impact-parameter Beam Functions)B αβ n,n onto the standard QCD PDFs as defined in Eq. (43) and are related to the iBFs as in Eq. (46). The iBFs are defined without a soft zero-bin subtraction but with an ultrasoft zero-bin subtraction. This is a consequence of the equivalence of soft zero-bin and soft subtractions [73][74][75] which we have discussed in some detail in Section III C. If one works with iBFs defined with a soft zero-bin subtraction, the purely collinear iBFs defined earlier in Section III C, then these objects would naturally match onto purely collinear PDFs which differ from the standard QCD PDFs by soft zero-bin subtraction terms. For more details on the relationship between standard PDFs and purely collinear PDFs see [74,75].
The n-collinear andn-collinear iBFs describe the physics of collinear emissions in the n andn directions respectively but in addition, as already mentioned, they each encode the same soft region due to the absence of a soft zero-bin subtraction. As a result the soft region, which encodes the physics of soft emissions and soft virtual effects, ends up being double-counted between the two iBFs. The iSF plays the role of precisely subtracting a soft region so that in the convolution over the two iBFs and the iSF there is only one soft region as required. After the iBFs are matched onto PDFs as in Eq. (43), we end up with a convolution over the two collinear functions I αβ n,n;g,i and the iSF which defines G ij as in Eq. (51).
The iBFsB n,n are defined in covariant or non-singular gauges where the gauge potential vanishes at infinity. Under non-singular gauge transformations, the iBFs are invariant. However, the iBFs are not invariant, as defined, under singular gauge transformations with boundary conditions where the gauge potential is non-vanishing at infinity. In singular gauges like the light-cone gauge, the definition of the iBF must be modified [44][45][46][47] with a transverse gauge link in order to be invariant under singular gauge transformations. However, the restriction of the iBFs to covariant gauges suffices since the quantity G ij which is a convolution over the two iBFs and the iSF is fully gauge invariant. This can be seen in two ways. First, since the iBFs and the iSF were obtained after applying the soft-collinear decoupling property of the SCET p T Lagrangian on the manifestly gauge invariant operator O defined in Eq. (21), the total convolution over the iBFs and the iSF must again be gauge invariant. Second, the product of the hard coefficient H and the perturbative function G ij , appearing in the factorization theorem in Eq. (54), just gives the perturbative cross-section for Higgs + multi-parton production as is manifest before applying the soft-collinear decoupling transformation. Together with the gauge invariance of H, the quantity G ij must again be gauge invariant. The additional functional dependence of the iBF on the variable t n also removes the singularity that occurs when this variable is fully integrated over, as it is in the transverse-momentum dependent PDF. In our factorization theorem this variable is instead constrained by the external kinematics of the process.
The factorization of G ij into the convolution I αβ n;g,i ⊗ In ;g,jαβ ⊗ S −1 can be quite useful in the implementation of higher order radiative corrections in α s (p T ). This is because the computation of the individual functions I αβ n,n;g,i , defined in covariant gauges, and S −1 is much simpler than the full computation of G ij without factorization.
Another property of the iBFsB αβ n,n and the matching coefficients I αβ n,n;g,i is their index structure. The indices α, β run over the transverse directions. Note that these indices remain uncontracted in each iBF due to their non-trivial dependence of the iBFs on the impact parameter b α ⊥ . Instead, these indices are contracted between the n-collinear and thē n-collinear functions as seen in Eqs. (41) and (51). This contraction of the indices between the n andn collinear functions is necessary to reproduce dot products between the transverse momenta of final state particles in the n andn directions which will appear at higher orders in perturbation theory.
We next discuss the formulation of factorization in impact-parameter space vs momentum space. So, far we have discussed the factorization theorem in terms of the impact-parameter space iBFs and iSF as seen in Eq. (51). In particular, the matching of the iBFs onto the standard QCD PDFs, as defined in Eq. (43), must be done in impact-parameter space. It is only in impact-parameter space that the infrared singularities of the iBFs become manifest as poles in dimensional regularization, and can therefore be seen to be the same as those of the standard PDFs. However, once the impact-parameter space matching coefficients are obtained one can Fourier transform to momentum space as in Eq. (52) and obtain the perturbative function G ij entirely in terms of the momentum space matching coefficients and iSF as seen in Eq. (53). In this momentum space factorization there is no reference to any impact-parameter. This form also makes manifest that the Landau poles which appear in the impact-parameter space CSS formulation [31] are avoided here. We discuss this point further below in our comparison with the standard approach.

B. Comparison with CSS approach
The classic QCD analysis of resummation of low transverse momentum logarithms expresses the cross section in the low p T region as [31] The sum is over parton species labeled by a, b, while x A,B denote the equivalent parton fractions x A,B = e ±Y m h /Q respectively. The functions A, B, and C have perturbative expansions in α s , while b 0 is an arbitrary constant chosen for computational convenience. One significant difference between this result and our approach outlined in the previous section is the appearance of the Landau pole of the strong coupling constant when µ 2 = 0 in the exponent. To deal with this singularity, several modifications of this formula are employed, including a deformation of the b ⊥ integration contour [36,76,77] and the introduction of a phenomenological model to cut off the b ⊥ → ∞ region [78]. In our approach the most natural choice for the scale which controls the lower limit of the RG evolution is µ L = p T . This can also be understood by noting that the perturbative function G ij is independent of the impact parameter, in both the impact-parameter and momentum-space formulations of the factorization theorem, and depends on p T and µ and no other dimensionful scales. Furthermore from the structure of the factorization theorem, we see that the logarithms of m h /p T are summed by the RG evolution of the hard coefficient H(x 1 , x 2 Q 2 , µ) which multiplies the function G ij and also has no reference to an impact parameter. In the effective theory, non-perturbative effects such as those indicated by the appearance of the Landau pole are encoded in operators suppressed by Λ QCD /p T . When p T ∼ Λ QCD , the expansion in this parameter breaks down, and a model of G ij fit to data can be used analogous to the standard approach. However, no reference to a non-pertubative function is needed above Λ QCD . Previous comparisons of b-space and momentums-space resumma-tion formalisms have indicated numerical agreement between the obtained results down to p T ∼ few GeV [79]. At this stage, power-suppressed operators presumably give important contributions. The use of SCET allows such effects to be studied in a systematic way. The avoidance of the Landau singularity also simplifies the matching of the resummed result to the fixed-order expression. In the usual approach, a large cancellation between the resummed component and the fixed-order QCD contribution occurs, leading to potential instabilities in the matched distribution. This cancellation typically occurs numerically because of the introduction of a non-perturbative model for the large b ⊥ region. Since it can be arranged analytically if the b ⊥ integrals can be done, avoidance of the Landau pole is useful for this purpose also (we note that the matching to fixed-order QCD results can be made smoother in the CSS approach by an appropriate change in the logarithm resummed [35]). We note that because of the effective theory expansion in Λ QCD /p T , it is not possible to calculate exactly the p T = 0 result for very largeQ as in the b-space approach [28]. However, the model of the non-perturbative region can be arranged to reproduce this result [79].

C. Comparison with previous SCET analyses
Resummation of the of the Higgs production transverse momentum distribution has been studied previously in [43] by a perturbative matching calculation onto SCET followed RG running. An explicit all orders factorization formula at the p T scale, corresponding to our Eq.(4), was written in [42]. The structure of our factorization theorem differs from that in [42] in several key ways. Our factorization theorem involves a convolution over the residual light-cone momenta that connects the collinear and soft sectors necessary by momentum conservation which is absent in [42]. This additional dependence is required in order to extend the factorization and resummation to more exclusive quantities where final-state hadronic quantities are also restricted, as is clear from the dependence on the light-cone momenta in the beam functions of Ref. [41]. The iBF in our factorization theorem has transverse indices which are contracted between the n-collinear andn-collinear sectors necessary to generate dot products between the transverse momenta of partons in the n and n directions in the final state. This index structure is not present in [42]. We disagree with this result, as it is not the structure that arises from the SCET analysis. In our formulation, the matching onto the standard QCD PDFs was done separately in the n andn sectors making full use of the SCET factorization. The work of [42] does the matching to PDFs by matching the full product of soft and collinear functions onto the QCD result, as required if the index structure of the collinear sectors is dropped. Finally, our work incorporates the formalism of zero-bin subtractions necessary to avoid double counting the soft region which was not addressed in Ref. [42].

V. FIXED ORDER ANALYSIS
In this section we give details of the perturbative computation of the hard Wilson coefficient, the n-collinear andn-collinear iBFs, and the iSF which appear in the factorization theorem. We emphasize that the iBF and iSF derived here are universal, and can be applied to other gluon-initiated processes at hadron colliders.

A. Calculation of the QCD → SCET p T Wilson coefficient
We begin by discussing the matching of QCD onto SCET p T in order to extract the Wilson coefficient C(ω 1 , ω 2 , µ). The Wilson coefficient can be extracted from the relation presented in Eq. (27) by computing radiative corrections to the matrix elements of both the QCD and SCET operators and encode their difference in C(ω 1 , ω 2 ). For the tree level and one loop matching one can compute the matrix elements h|O QCD |p 1 ,p 2 and h|O|p 1 ,p 2 in QCD and SCET p T respectively wherep 1 µ =n ·p 1 n µ 2 andp µ 2 = n ·p 2n µ 2 denote the momenta of the initial state n-collinear andn-collinear gluons. The diagrams contributing at next-to-leading order in α s in SCET p T are shown in Fig. 3. Labeling these graphs from left to right they take the form  so that the SCET p T operator is multiplicatively renormalized. With on-shell external gluons and using Feynman gauge, the quantities V a,b,b0 take the form We note that the collinear graphs in Figs. 3b and 3c require a zero-bin [72] subtraction given by the V b0 term in order to avoid over-counting the soft region. These integrals are all scaleless and vanish in dimensional regularization. The Wilson coefficient can therefore be extracted directly from the finite part of the full QCD result in dimensional regularization and is well known in the literature [15,80]. Through next-to-leading order, it is given by where the first term is just the result of tree level matching quoted earlier in Eq. (28) and we have used the property C(ω 1 , ω 2 , µ) = C(ω 1 ω 2 , µ) to write the LHS above. We note that the 11/2 in the next-to-leading order expression arises from integrating out the top quark loop to produce an effective Higgs-glue vertex. The hard Wilson coefficient H(x 1 x 2 Q 2 , µ) = |C(x 1 x 2 Q 2 , µ)| 2 appearing in the factorization theorem can be obtained from the above equation at next-to-leading order.

B. Calculation of the iBF
In this section we present results for the calculation of the iBFs as defined in Eqs. (36), (33), and (42). We compute the n-collinear iBF by inserting a complete set of states as and then computing the product of matrix elements. Recall thatB αβ n denotes the iBF without a soft zero-bin subtraction as opposed to B αβ n which is defined with a soft zero-bin subtraction. An analogous expression exists for then-collinear iBF. In this section we focus on the n-collinear iBF, since then-collinear iBF can be calculated in an analogous fashion. The lowest order result for the iBF is obtained by choosing |X n = |0 and computing the tree level matrix elements to get Higher order radiative corrections to the term in the iBF with |X n = |0 correspond to pure virtual corrections with no real emissions in the final state. The one loop virtual correction corresponds to the first diagram in Fig. 4 and its conjugate and is given bỹ where V b is given in Eq. (57) and w.f. denotes the external wave-function graphs. Thus, this virtual correction is also scaleless and vanishes in pure dimensional regularization. Next we compute the single real gluon emission contribution to the iBF corresponding to a single gluon with momentum p g in the final state: |X n = |p g . The single gluon emission contribution to the iBF corresponds to the last three diagrams in Fig. 4. Only the second diagram in Fig. 4 contributes if we use the physical gluon polarization sum pols.
Direct calculation of this graph gives where we work with the convention d = 4 − 2 . After performing the integrations over the light cone momenta p + g and p − g the single gluon emission contribution to the iBF can be brought into the form Three-dimensional quantities have been denoted by a vector symbol. The result is the same in any R ξ gauge due to the structure of the Wilson-line coupling. We note that if t + n was integrated over, as it is when considering transverse-momentum dependent PDFs, singularities would occur at the integration boundaries. This is avoided in the iBF. The variable t + n is instead set by the external kinematics through the convolution structure of the factorization theorem, as we will see in a later section when comparing to fixed-order QCD.
In order to perform the integral appearing in Eq. (64), we require the following two integrals: and we have used the notation | b ⊥ | = b ⊥ . To utilize these results, it is most useful to expand the integral over the second piece of Eq. (64) in terms of form factors. We can schematically write this as We can express these in terms of the two integrals calculated in Eq. (65): .
Furthermore it is convenient to parameterize the iBF in terms of two form factors F 1 and F 2 as Note that if the indices α, β are contracted the term proportional to F 2 vanishes. Using Eqs. (64), (65), (66), and (67) we find the following results for the form factors F 1 and F 2 The form factor F 1 has ultraviolet and infrared divergences that are regulated by in our pure dimensional regularization calculation. On the other hand, the form factor F 2 is finite. It is convenient to remove the fractional powers of the Bessel-function arguments by writing them in terms of hypergeometric functions in the following way: Since 0 F 1 (ν, 0) = 1, and therefore J 2 (z) ∼ z 2 , it is clear that F 2 is finite and that we can safely set = 0 to get For the form factor F 1 , we first rewrite its expression in Eq. (69) as where we have have switched to an M S definition of µ by replacing µ 2 → µ 2 (4π) − e γ . We have also used the partonic center-of-mass energy squaredQ 2 = x 1 x 2 Q 2 to form dimensionless ratios where required; this is convenient for the comparison with fixed-order QCD calculations, which is discussed in a later section. We note that the ratio of scales which appears in this expression is t + n /µ 2 , indicating that we should choose µ ∼ p T to minimize logarithms since in the partonic calculation t + n ∼ p 2 T . We write the expansion of the form factor F 1 in as and derive the following expressions The result for the single gluon emission contribution to the n-collinear iBF in pure dimensional regularization is given by Eqs. (68), (71), (73), and (74). The expression for thē n-collinear iBF is obtained by repacing t + n → t − n and x 1 → x 2 in the expression for the n-collinear iBF.

C. Soft function: real emission
In this section we give results for the computation of the iSF which was defined earlier as where the position space soft function that appears on the RHS above is defined in Eq. (33). The iSF in the factorization theorem, before doing the Higgs phase space integrals, has the argumentsω as seen in Eq. (47). For convenience we introduce the notation which we will often use in this section. We compute the iSF by inserting a complete set of soft states in the position space soft function as Higher order corrections to the term with |X s = |0 corresponds to virtual graphs with no real emissions in the final state. At one loop, the virtual corrections correspond to the first diagram and its permutations in Fig. 5 which gives the result where I s is the scaleless integral and vanishes in pure dimensional regularization. Next we compute the contribution to the iSF from the real emission of an soft gluon corresponding to choosing |X s = |k for a gluon of momentum k, as shown in the second diagram of Fig. 5. Explicit computation gives Switching to an M S definition of µ and performing integrals as before, we can derive the following expression: 4 .
(83) The expansion in proceeds identically to that for the iBF. Defining the expansion we arrive at the result We have used the scale Q to define dimensionless ratios in logarithms and plus distributions, as in the calculation of the iBF. We note that when this result is plugged into the full expression for G ij in Eq. (49), the logarithm that appears is ln tû /(µ 2Q2 ) , wheret,û are the partonic Mandelstam variables. This indicates that we should choose µ ∼ p T in the soft function, and is the motivation for our decision to set µ S = µ T in Section IV.

D. iBFs to PDFs
In order to be able to sensibly match the iBFs onto the PDFs with a finite Wilson coefficient, as indicated in Eq. (43), the infrared divergences in the iBF and PDF must match. In this section, we show that this can be manifestly seen from the integrand level expressions for the iBF and PDF at the order we are working. We first note that the tree level expression for the iBF to PDF matching is given by which can be straightforwardly verified from Eq. (43). Next we consider the integrandlevel expression for the real gluon emission contribution to the iBF in Eq. (64) where the integrals over the light cone momentum components of the final state gluon have already been performed. The infrared divergences now occur in the limit where the perpendicular momentum of the final state gluon vanishes p g⊥ → 0. In this limit, the iBF real-emission contribution in Eq. (64) becomes (87) A calculation of the real-emission contribution to the PDF defined in Eq. (44) can be done inserting a single-gluon state between the field operators in the same way as was done for the iBF. The result at the integrand level is from which we see that Thus, the infrared divergences which appear upon integration over p g⊥ are identical for both the iBF and the PDF, up to an overall multiplicative factor. Similarly, we note from Eqs. (61) and (60) that the contribution of one loop virtual corrections to the iBF take the form so that once again the infrared structure is the same for the iBF and PDF virtual contributions up to the same multiplicative factor. Thus, the infrared divergences in the sum of the real and virtual graphs for the iBF can be matched to those in the PDF as in Eq. (43) with the finite Wilson coefficient given in Eq. (86). Additional contributions to the Wilson coefficient will come from the finite terms in the iBF. In fact, since the PDFs are scaleless and we have shown that the infrared divergences of the iBF and the PDF match up, the Wilson coefficient for the iBF to PDF matching is just given by the finite part of the iBF in pure dimensional regularization as indicated in Eq. (46). For a more detailed explanation of this point see [81]. Using Eqs. (46) and (68) we can now write the iBF to PDF matching coefficient at the first order beyond tree level as where the form factors F 1;0 and F 2 are given in Eq. (74) and (71) respectively. We perform an explicit check of this expression in section VII by using this expression in the factorization theorem and comparing the fixed order result in full QCD.

VI. RUNNING
In the factorization formula of Eq. (54), the logarithms of m h /p T are summed via the RG evolution of the hard coefficient H(x 1 x 2 Q 2 , µ) between the scales µ Q ∼ m h and µ T ∼ p T . The factor H(x 1 x 2 Q 2 , µ Q ; µ T ) in the factorization formula denotes this RG evolved hard coefficient. The logarithms of Λ QCD /p T are summed via the DGLAP evolution of the PDFs up to the µ T ∼ p T scale. The DGLAP evolved PDFs are denoted by f i/P (x, µ T ). In this section we give details of this RG running above and below the p T scale.

A. Running above p T
The anomalous dimension of the hard coefficient H can be obtained from the anomalous dimension of the operator which was first defined in Eq. (21). The RG equation for O has the form where the anomalous dimension γ O is given by and Z O is the renormalization constant that relates the bare operator O b to the renormalized operator O by We use the notation The anomalous dimension of O(n ·p 1 , n ·p 2 ) is typically written in the general form where the cusp (Γ O ) and non-cusp (γ O ) pieces have expansions in α s as Equivalently, one can think of the running of the Wilson coefficient C of the operator O which has the RG equation The running of absolute value squared of the Wilson coefficient H(Q 2 , µ) = |C(Q 2 , µ)| 2 , which is the quantity that appears in the factorization theorem, is given by so that where Γ H and γ H have the α s expansions The result for the anomalous dimension γ O is well-known [16,43,80,82,83]. Since the SCET p T graphs that renormalize the operator O are scaleless, as discussed in Section V, the anomalous dimension can be read off [81] from the infrared pole structure of loop corrections [80,82] to the QCD operator O QCD defined in Eq. (20). In this section we reproduce this known result at one loop by explicit calculation of the SCET p T graphs using off-shell external gluons to regulate the infrared divergences and dimensional regularization to regulate the ultraviolet divergences. We do this calculation in both Feynman gauge and using the background field method as an additional check on our calculation.

Anomalous dimension computation in Feynman gauge
At one loop, there are three graphs in SCET p T , in addition to the wave function graphs, which determine the anomalous dimension γ O = γ c , as shown in Fig. 3. The first graph involves a soft gluon loop emanating from the soft Wilson lines and the remaining two graphs involve a (n,n)-collinear gluon originating from the (n,n)-collinear Wilson line and attaching to the (n,n)-collinear gluon. These diagrams are given by where the quantities V a (p 1 ,p 2 ), V b (p 1 ), and V b0 (p 1 ) with off-shell momentap 1,2 are given by and note that in the limit p 2 1 , p 2 2 → 0 these integrals reduce to the expected results of Eq. (57). Note that the (n,n)-collinear graphs require the zero-bin [72] subtraction terms V b0 (p 1 ) and V b0 (p 2 ) respectively. These zero-bin subtractions are necessary to remove the soft region in the collinear graphs and thus avoid double counting the soft region. The results for these graphs and the zero-bin subtraction term are where we have used an off-shell infrared regulator ∆ 2 =p 2 1 =p 2 2 . We note that the offshellness of the external gluons does not completely regulate the infrared divergences individually in V b and V b0 . Furthermore, V b contains a UV-IR mixed double pole. However, in the zero-bin subtracted combination V b − V b0 the infrared divergences are completely regulated by off-shell external gluons and all pole terms are pure UV. This is critical for being able to properly extract the anomalous dimension from the pole structure. The sum of these graphs is Fig. 3a + Fig. 3b + Fig. 3c Note that the mixed UV-IR terms proportional to 1 UV ln ∆ 2 cancel between the soft and collinear graphs so that the UV and IR divergences are clearly separated in the sum of all graphs. Along with the gluon wave-function and strong coupling renormalization constants, we can now extract the anomalous dimension of O. The bare operatorÔ b obtained by expanding the bare B b ⊥ fields in the strong coupling and keeping the leading term, can be related to the renormalized operatorÔ aŝ so that the renormalized operator can be written aŝ where Z 3 and Z g are the gluon wave function renormalization and coupling renormalization constants respectively and we have defined At one loop, in the MS scheme, the gluon wave function and the strong coupling renormalization constants are known [84]: We see from the above results that in the MS scheme δZ O at one loop is given by where the subscript "div" above indicates the UV divergent part of the sum of graphs. Putting these results together we get the one loop result Using N c = C A = 3 this result simplifies to from which we get the well known result for the one loop anomalous dimension γ In the notation of Eqs. (96) and (97) this result for the one loop anomalous dimension corresponds to Thus, the non-cusp piece γ O 0 vanishes. Putting this together with the notation in Eqs. (99),(100), and (100) the anomalous dimension of the hard coefficient H at one loop is by and also contains only a cusp contribution.

Anomalous dimension computation in background field method
One can also compute the anomalous dimension using the background field method with Feynman gauge for internal lines in the loop graphs, which we show as an additional check on our SCET computational techniques. In this method, the combination gA µ n⊥ is not renormalized. In this case, the result for the soft graph V bfg a is unaffected The collinear graph V b is modified to and the zero-bin subtraction term is unaffected One can now extract the anomalous dimension following the same procedure as in the Feynman gauge computation of the last section. The bare operator obtained at leading order by expanding the bare B b ⊥ field can again be written in terms of the renormalized operator as in Eq. (106), so that the renormalized operator can be written aŝ (120) Note that in the background field method, compared to Eq. (107), the wave function (Z 3 ) and coupling constant (Z g ) renormalization constants cancel out and do not appear. In other words, the combination gA µ n,n⊥ does not undergo renormalization in the background field gauge. Thus, we get which gives the same result as in Eq. (113) and thus leads to the same anomalous dimension as derived from the Feynman gauge calculation.

B. Running below p T
Logarithms of Λ QCD /p T are summed via the standard DGLAP equations for the PDFs. For completeness we write down the RG equation for the gluon PDF at leading order where P g←g (z) is the standard gluon splitting function and we have ignored contributions from the quark and antiquark PDFs. This DGLAP equation is used to evolve the gluon PDF up to the µ T ∼ p T scale. These evolved PDFs, denoted by f g/P (x, µ T ), are the ones that appear in the factorization theorem as in Eq. (54).

VII. CONSISTENCY CHECKS
A. Cross-section comparisons with full QCD In this section we show that the iBFs and iSF computed in fixed-order perturbation theory in Section V reproduce the QCD cross section up to terms suppressed by the ratio p T /m h .

Leading-order cross section
We begin by reproducing the leading-order cross section of perturbative QCD for gg → h, with no final state gluon emissions, using our factorization formula in Eq. (50). We utilize the leading-order expressions for the Wilson coefficient H, the iBFs and the iSF: Inserting these results into Eq. (50), the integrals over x 1,2 , t n , tn, and b ⊥ can be easily performed to give which measures the amount of energy released into final-state particles besides the Higgs (for z = 1, all energy in the partonic scattering process goes to the Higgs). It is convenient for this expression and for later results to switch to partonic Mandelstam variablesû,t defined byû − m 2 , and introduce a partonic differential cross section via The partonic cross section at leading order takes the form Upon integration over the variablesû andt, this reproduces the leading order result of Eq. (13).

Next-to-leading order cross section
The next-to-leading order prediction of QCD has contributions from virtual corrections to gg → h and from the emission of one final state gluon gg → gh. To compare with this next-to-leading order prediction of QCD, we expand out all factors appearing in Eq. (50) through next-to-leading order in perturbation theory. After expansion, we can write the cross section in the schematic form (127) The next-to-leading order expressions for the terms appearing in this cross section can be found in Eqs. (58), (85), and (91). Before explicitly computing the various contributions, we note several simplifying features that appear when the iSF and iBF terms are combined.
• When the combinations I n S −1(1) is performed, several terms that appear in the individual expressions cancel. In particular, it is simple to check that the following combination that appears in the the cross section vanishes upon integration over t n and tn: After making these simplifications, we have four contributions to the cross section at nextto-leading order. We now discuss each contribution in detail. To simplify the notaton, we will set µ 2 =Q 2 = x 1 x 2 Q 2 , where x 1,2 are the arguments of the PDFs.
After integration over the Mandelstam variables, the contribution of this term to the total partonic cross section is clearlŷ 2.
n S −1(1) : After implementing the simplifications described above, the only terms remaining from the iSF in Eq. (85) that contribute to the cross section are given by The first term gives a contribution proportional to the leading order cross section. When computing the contribution of the second term to Eq. (50), the following identity is needed: This follows immediately from writing the hypergeometric function in terms of a Bessel function using Eq. (70), and then applying the orthogonality relation between J 0 functions of different arguments. After performing the b ⊥ integral using the identity above, and expressing t max n,n in terms of Mandelstam invariants using their definitions in Section V, we find the following contribution to the cross section: To find the total cross section resulting from this term, we note that the range of thê t integration is −Q 2 (1 − z) ≤t ≤ 0. The definition of the plus distribution 1/[x] + as the O( 0 ) term in the expansion of x −1− allows us to derive the integral Using this result, it is straightforward to derivê 3. H (0) I (1) n I (0) n S −1(0) : The terms from the next-to-leading order Wilson coefficient I n in Eq. (91) that contribute to the cross section in Eq. (50) are The calculation of the first term in this expression proceeds identically to the previously discussed case. Upon substitution of the second term into the cross section in Eq. (50), the delta-functions constraints set x 1 /x 1 = z. Using the relation between hadronic and partonic Mandelstam invariants, it can also be shown that they setû = 0. Since the Higgs boson transverse momentum can be written in terms of Mandelstam invariants as p 2 T =tû/Q 2 , the second terms contributes only for zero transverse momentum. The final contribution of the I (1) n term to the differential cross section is Upon integration over the Mandelstam invariants, we find the following contribution to the total cross section: An identical calculation as that outlined above yields the following results for the differential and total partonic cross sections contributions from In: where the differential cross-section term above differs from the corresponding contributions from I n only by the interchangeû ↔t.
We now compare the differential cross section for p T > 0 and the total cross section with the full QCD expressions given in Eqs. (12) and (13). Summing the four contributions to the differential cross section, and keeping only terms which contribute for p T > 0, we obtain the following result: We have dropped the plus prescription on the Mandelstam invariants since we are only considering the region p T > 0. This expression agrees with the QCD result in Eq. (12) up to terms which are finite asû,t → 0. Such terms can arise from suppressed operators in the EFT, and not from emission of collinear or soft gluons. Our result therefore correctly reproduces that of QCD in the limit of soft or collinear gluon emission, as expected.
The total cross section obtained by combining the four contributions above iŝ The terms with δ(1 − z) and ln (1 − z) in this result agree with the QCD cross section in Eq. (13). Terms non-singular in the limit z → 1 can arise also from emission of hard jets with the Higgs, and therefore agreement in the non-singular polynomial in 1 − z is not expected. The description of these pieces requires power-suppressed operators in the EFT. Our total cross section therefore properly reproduces the expected result from fixed order QCD.
B. Top down vs bottom up running: impact-parameter space Section VI A was devoted to the RG running of the hard Wilson coefficient H(x 1 x 2 Q 2 , µ). First H(x 1 x 2 Q 2 , µ Q ) is obtained from a matching calculation at the hard scale µ Q ∼ m h . Next, after computing the anomalous dimension of H(x 1 x 2 Q 2 , µ), it is RG evolved down to the µ T ∼ p T scale to obtain H(x 1 x 2 Q 2 , µ Q ; µ T ) which has the logarithms of m h /p T resummed. We refer to this procedure as top-down running. Equivalently, one can compute the anomalous dimensions of the iBFs and the inverse soft function and perform an RG running from the µ T ∼ p T scale, where the iBFs and inverse soft function live, up to the hard scale µ Q . We refer to this procedure as bottom-up running. Scale invariance of the cross-section implies that the total anomalous dimension of the convolution of the iBFs and inverse soft function in the t ± n,n variables must exactly cancel the anomalous dimension of H(x 1 x 2 Q 2 , µ). In this section, we use this consistency condition [85,86] to obtain the anomalous dimension of the iBF which we will compare with what was obtained from the fixed order calculation. We will derive the anomalous dimension of the iBF in a third way in section VII D as yet another check.
We define the quantity D as where X(µ) is the convolution of the iBFs and the inverse soft function The factorization formula for the differential cross-section at the p T scale in Eq. (41) contains the quantity D. Furthermore, all the renormalization scale dependent quantities are contained in D. Scale invariance of the cross-section then requires that D has zero anomalous dimension or equivalently that the µ-dependence cancels between H(x 1 x 2 Q 2 , µ) and X(µ).
We begin by defining the renormalization constants which relate the bare and renormal-ized quantities as where the sunscript b on the objects on the RHS denotes bare quantities and the objects on the LHS are the renormalized quantities. Note that the renormalization of the hard coefficient of H is multiplicative. On the other hand, the renormalization of the the iBFs and the inverse soft function involve a convolution. This implies that the hard coefficient has local running while the iBFs and the inverse soft function undergo convolution running. The consistency of top-down and bottom-up running implies that the total convolution of the iBFs and the inverse soft function undergoes local running that precisely cancels the running of H(x 1 x 2 Q 2 , µ).
The scale invariance of D leads to a consistency condition that relates the various renormalization constants as To derive this consistency equation, we first write D entirely in terms of bare quantities and use the first equation in Eq. (145) to replace H b . Next we write D entirely in terms of renormalized quantities and then use the last three equations in Eq. (145) to replace the renormalized iBFs and the iSF in terms of the corresponding bare quantities. Finally by equating these two ways of writing D we arrive at Eq. (146). We can make the content of Eq. (146) more transparent by noting that we can write the renormalization constants as a sum of a leading order piece and perturbative corrections as Zn(t) = δ(t) + δZ n (t), where We have used the fact that the dependence on t n /Q and tn/Q of the soft renormalization constant is separable [86] as was seen in the explicit calculation in section V. The consistency condition in Eq. (146) can now be written as The leading order results for the soft and hard coefficient renormalization constants from sections V and VI A respectively are from which we can solve for the one loop renormalization constant of the iBF as δZ (1) n,n (t) = which is in agreement with the ultraviolet pole terms for the iBF derived in section V. Thus, we have an explicit demonstration of the equivalence of the top-down and bottom-up running at leading order. Eq. (151) implies that the leading order anomalous dimension of the iBF is given by In section VII D, we provide another independent check by computing the anomalous dimension by taking the derivative with respect to µ of the RHS of Eq. (43) which expresses the iBF in terms of a matching coefficient and the standard QCD PDF.

C. Top-down vs bottom-up running: momentum space
The consistency of the top-down and bottom-up running can be checked in yet another way by looking at the momentum space factorization. We wrote down the momentum space factorization earlier in Eqs. (50) and (53). However, in order to demonstrate the equivalence of the top-down and bottom-up running above the p T scale, we need the momentum space analog of Eq.(41) which is given by The momentum space functions are given in terms of the impact-parameter space functions by the Fourier transforms B αβ n,n (x, t, k ⊥ , µ) = We have made use of Eq. (38) to recast the factorization theorem in terms of purely collinear iBFs, defined with a zero-bin subtraction as in Eq. (37), and the soft function. In this momentum space version of the factorization formula at the p T scale, the running of the purely collinear iBFs and the soft function is determined entirely by scaleless virtual graphs as shown in Figs. 4 and 5. The real gluon emission graphs for the momentum space purely collinear iBFs and the soft function are finite and do not contribute to the running. This is in contrast to the impact parameter space factorization, examined in the previous section, in which the real emission graphs have UV divergences and do contribute to the running. The consistency of top-down and bottom-up running should work out in both impact-parameter and momentum space. In this section we show how the virtual graphs for the purely collinear iBFs and the soft function cancel the running of H(x 1 x 2 Q 2 , µ) to maintain the µ-independence of the differential cross-section. For one loop running, we show that where F is defined to be the product of the purely collinear iBFs and the soft function The bare and renormalized F are related as and writing the left-hand side in terms of renormalized fields and couplings we get where the quantities B αβ n , Bn αβ , S are all written in terms of renormalized fields and couplings. We suppress the arguments in Eq. (156) for notational simplicity.
The zeroth order expression for F is where and the contribution of the virtual graphs in Figs. 4 and 5 to F is The one loop virtual graphs contributing to the purely collinear jet functions and the soft function are given by where expressions for V a (p 1 ,p 2 ), V b (p 1 ), and V b0 (p 1 ) were given earlier in Eq. (57) and in Eq. (103) with off-shell momentap 1,2 to regulate infrared divergences. One can rearrange the terms in Eq. (162), at the level of the integrand, using Eq. (163) to get The first line in curly brackets is precisely the sum of graphs in Fig. 3 which determine the anomalous dimension at one loop for O(n ·p 1 , n ·p 2 ) and the second line in curly brackets is just the sum of conjugate diagrams. From Eqs. (159), (164), (105), and Eq. (110) we get so that as expected for the scale invariance of Eq. (153).

D. Running after matching
For completeness, we will also verify that the running after matching our results to the standard PDFs cancels the running of the hard Wilson coefficient H(Q 2 , µ 2 ). To do so, we begin with the relevant quantity appearing in the factorization formula of Eq. (47): We have rearranged the x i integrals appearing in Eq. (47) into a form convenient for this calculation. The quantity X must obey the renormalization group equation in order for the cross section to be scale invariant and for the operator O defined in Eq. (92) to have the correct anomalous dimension. It is simplest to derive the running of X by differentiating directly the beam function defined in Eq. (43). This accounts for the running of the product I ×f g/p in Eq. (167). From Eqs. (43), (46), (86), and (91), we can check that the terms appearing in the n-collinear beam function that contribute to the µ dependence arẽ The running of the coupling constant and PDF are given by the standard expressions We now differentiate the beam function in Eq. (169) with respect to µ, and arrive at the result (171) We note that this result is in agreement with the anomalous dimension found in Eq. (151). An identical calculation for the othern-collinear beam function yields We next differentiate the finite part of the iSF defined in Eq. (85), and find We note that it is possible to manipulate the result of Eq. (85) such thatQ rather than Q appears in the logarithms and distributions, and we have done so in writing this expression.
We are now ready to plug all of this into X and determine its running. The leading order expression for X can be easily determined forom the results of Section V to be Upon substitution of the derivatives computed above in Eqs. (171), (172), and (173) into the derivative of X, we find the following result through next-to-leading order: The computation is simplified by noting that all plus distribution terms appearing in this expression cancel due to the following identity: such that only the terms in Eqs. (171), (172), and (173) with an explicit ln(Q 2 /µ 2 ) contribute to the running. This verifies that the running of X matches the required form in Eq. (168).

VIII. CONCLUSION
We have derived a factorization theorem for the transverse momentum and rapidity distributions of the Higgs boson at hadronic colliders, in the region m h p T Λ QCD , employing the methods of effective field theory. The factorization theorem naturally allows for a resummation of the low p T region and can be straightforwardly generalized to other processes that involve the production of one or more color neutral particles, such as the Drell-Yan process. The problem of resummation in the low p T region has been extensively studied in the QCD literature [27-29, 31, 35] and more recently with effective field theory methods [42,43]. However, the form of our factorization theorem and the method of our analysis contains several new and interesting elements that are worth pursuing further.
The derived factorization theorem takes the form shown in Eq. (1) which is a convolution between the perturbative coefficients H(x 1 x 2 Q 2 , µ Q ; µ T ) and G ij (x 1 , x 1 , x 2 , x 2 , p T , Y, µ T ) and standard QCD PDFs. This form was obtained by matching onto a sequence of effective field theories. First the top quark is integrated out to obtain the standard effective coupling of the Higgs to gluons. The hard coefficient H(x 1 x 2 Q 2 , µ Q ; µ T ) is obtained by matching the effective Higgs operator onto an effective SCET p T operator, integrating out the hard scale m h in the process, followed by RG evolution between the scales µ Q ∼ m h and µ T ∼ p T . The RG evolution of the hard coefficient, determined in the effective theory SCET p T , sums up logarithms of m h /p T . The perturbative coefficient G ij , evaluated at the µ T ∼ p T scale, takes the form in Eq. (2) which involves a convolution over the n andn collinear functions I αβ n,n;g,i and the inverse Soft Function S −1 . This form of G ij into two collinear and one soft sectors is the result of the soft-collinear decoupling property of the leading order SCET p T Lagrangian. The collinear functions I αβ n,n are obtained from an OPE in Λ QCD /p T of the impact-parameter Beam Functions (iBFs)B αβ n,n which appear in the factorization theorem at the p T scale as in Eq. (41). At leading order in Λ QCD /p T the iBFs match onto the standard QCD PDFs with Wilson coefficients I αβ n,n;g,i as in Eq. (43). The form of the factorization theorem is pictorially summarized in Fig. 1.
The iBFs are universal objects in that they will appear in the p T distributions for any gluon initiated process of color neutral particles. The iBFs defined here will have analogs for quark-initiated scattering as well, and the generalization to such processes is straightforward. Similarly, the iSF is also universal and again has a straightforward analog for quark-initiated processes. The iBFs are defined without a soft zero-bin subtraction, so that in addition to describing the emissions and virtual effects of collinear gluons they also do the same for soft gluons. The presence of two iBFs in the factorization theorem leads to this soft region being double-counted. The iSF plays the role of subtracting a soft region so that there is only one soft region in the full factorization theorem as required. Each iBF is still defined implicitly with an ultrasoft zero-bin subtraction to avoid double counting with ultrasoft modes which are unaffected by p T constraints. Each iBF and correspondingly the collinear functions I αβ n,n;g,i are defined in the class of covariant or non-singular gauges. The iBF is not invariant under singular gauge transformations where the gauge potential is non-vanishing at infinity and must be modified. However, the quantity G ij in the factorization theorem which is a convolution over the two iBFs and an iSF is fully gauge invariant. The factorization of G ij into two iBFs and an iSF can be quite useful for the computation of higher order radiative corrections at the p T scale since the computation of each of the factored objects is much simpler.
The factorization theorem can also be formulated entirely in momentum space with no reference to an impact-parameter. This form of the factorization theorem is given in Eq. (1) with the quantity G ij given in Eq. (53) entirely in terms of momentum-space objects. In this effective field theory formulation one avoids Landau poles associated with impact-parameter integrations in the usual CSS formulation. The use of SCET allows power-suppressed corrections of the form Λ QCD /p T and p T /m h to be systematically included. When p T ∼ Λ QCD , a model of the non-perturbative region can be included as in the usual approach.
The structure of our factorization theorem differs also from previous SCET analyses. We are led to gluon iBFs that have an index structure, which at higher orders is needed for the separate calculation of effects from each collinear sector.
We believe our study provides a new framework for the investigation of the low-p T region in hadron collider processes. Several interesting differences from previous approaches arise naturally in our derivation within SCET that require further investigation. We look forward to future research investigating the implications of the factorization theorem derived here.
Appendix A: Details of deriving factorization in SCET p T After matching n f = 5 QCD onto SCET p T as in Eq. (27), the differential cross-section in the Mandelstam variables u and t is given by where the delta functions in u and t pick out the appropriate final states. The overall factor of 1/4 in square brackets comes from the average over initial proton polarizations. The final state hadrons have been broken up into n-collinear,n-collinear, and soft states as |X = |X n , Xn, X s , since these are the states with non-zero overlap with the SCET p T operator O(ω 1 , ω 2 ) O(ω 1 , ω 2 ) = g µν h Tr T {S n (gB (0)µ n⊥ ) ωn S † n Sn(gB where S n (x µ ) = P exp ig s 0 ∞ ds n · A a s (x µ + sn µ ) (A4) with an analogous expression for S n . Since the collinear and soft fields are decoupled from each other, we can factorize the matrix element of O(ω 1 , ω 2 ) into collinear and soft sectors as hX n XnX s |O(ω 1 , ω 2 )|pp = hX n XnX s |g µν h Tr T {S n (gB µ n⊥ ) ωn S † n Sn(gB ν n⊥ ) ωn S † n } |p 1 p 2 = X s |Tr T {S n T A S † n SnT B S † n } |0 X n | B αA n⊥ ω 1 |p 1 Xn| B B n⊥α ω 2 |p 2 .
The absolute value squared of the matrix element now takes the form | hX n XnX s |C ⊗ O|pp | 2 where we have simplified the contraction of color indices in the last two equalities. Next we insert the identity operator 1 = d 4 p n d 4 pn d 4 p s δ (4) (p n − P Xn )δ (4) (pn − P Xn )δ (4) (p s − K Xs ), and we decompose the momentum components of order m h for the final state collinear particles into label and residual parts so that whereP − Xn ,P + Xn ∼ m h and K − Xn , K + Xn m h and we write the remaining momentum components as P +,⊥ Xn = K +,⊥ Xn , P −,⊥ Xn = K −,⊥ Xn , P µ Xs = K µ Xs .
Similarly, we write p − n =p − n + k − n , p + n =p + n + k + n , p +,⊥ n = k +,⊥ n , p −,⊥ n = k −,⊥ n , p µ s = k µ s , where againp − n ,p + n ∼ m h andk − n ,k + n m h . The delta functions in Eq. (A7) for the large collinear momentum components break up into Kronecker delta functions over the hard collinear label momenta and residual delta functions which we can write using the integral representation to get 1 = (A11) The exponential factors above involving the momenta of the final state collinear and soft states allow us to perform a shift in coordinates and perform a sum over the states X n,n,s states. Furthermore we break up the Higgs momentum as n · p h = n ·p h + n · k h , n · p h =n ·p h +n · k h , where n ·p h ,n ·p h ∼ m h and n ·k h ,n ·k h , k ⊥ h m h . The Higgs phase space integrals can now be written as a sum over label momenta and integrals over the residual momenta as We also break up the momentum conserving delta function into Kronecker deltas over label momenta and delta functions over the residual residual momenta as After all of the above manipulations, the SCET p T double differential cross-section can be brought into the form initial pols. Xn,Xn,Xs e ikn·x e ikn·y e iks·z × dω 1 dω 2 dω 1 dω 2 C * (ω 1 , ω 2 )C(ω 1 , ω 2 ) Next we perform the integral over x + to get the fourier transform of the n-collinear field B A 1n⊥β (x) with respect to x + as with an analogous equation for then-collinear field B B 2n⊥β (y). Next we combine the sum over labels and integral over residual as usual so that and absorb the residual momenta k − n and k + n into ω ( ) 1 and ω ( ) 2 respectively. We arrive at 2n⊥α ω 2 (0)|p 2 × δ (2) (k s + k n⊥ + kn ⊥ + k h⊥ ) Next we note that the n-collinear matrix element is proportional to δ(ω 1 − ω 1 ) which can be seen as follows: where we have inserted a complete set of states and used momentum conservation and then performed the sum overs states again. Analogously, then-collinear matrix element is proportional to δ(ω 2 − ω 2 ) so that one can perform the integral over ω 1,2 to get The iBFs, the zero-bin subtraction terms, and the soft function at first order in perturbation theory receive contributions from virtual (V) and real (R) emission graphs so that relation Z H (µ)δ(s n − t n )δ(sn − tn) = 1 Q 2 dt n dt n Z −1 n (t n − s n , µ)Z −1 n (t n − sn, µ) × Z −1 S −1 ( −t n + tn Q , −t n + t n Q , µ).