Collinear expansion for color singlet cross sections

We demonstrate how to efficiently expand cross sections for color-singlet production at hadron colliders around the kinematic limit of all final state radiation being collinear to one of the incoming hadrons. This expansion is systematically improvable and applicable to a large class of physical observables. We demonstrate the viability of this technique by obtaining the first two terms in the collinear expansion of the rapidity distribution of the gluon fusion Higgs boson production cross section at next-to-next-to leading order (NNLO) in QCD perturbation theory. Furthermore, we illustrate how this technique is used to extract universal building blocks of scattering cross section like the N-jettiness and transverse momentum beam function at NNLO.


Introduction
Our knowledge of the structure of quantum field theory (QFT) is rapidly advancing. On the one hand this steady progress allows us to answer fundamental questions about the interactions of nature by deriving precise predictions for the outcome of scattering experiments that can be compared with experimental observation. On the other hand we learn about the mathematical structures that underly this description. Progress in QCD perturbation theory has allowed us to venture to predictions at nextto-next-to-next-to leading order (N 3 LO) in the strong coupling constant for select inclusive and differential cross sections at the Large Hadron Collider (LHC) [1][2][3][4][5][6][7][8][9]. Resummation of kinematic limits of cross sections has reached the similarly astounding precision for a multitude of observables [10][11][12][13][14][15][16][17][18]. Nevertheless, the difficulty of describing the scattering of fundamental particles is ever rising with increasing demand for precision and for more complex observables. To overcome seemingly insurmountable complexity, parametric or kinematic expansions have proven highly effective. For example, expanding the gluon fusion Higgs boson production cross section around the production threshold of the Higgs boson allowed for the computation of the first hadron collider cross section at N 3 LO in QCD perturbation theory [1].
In this article we detail a technique for the efficient expansion of differential partonic cross sections for the production of a color singlet final state h in hadron-hadron collisions in the kinematic limit that all radiation produced alongside h is collinear to one of the collision axis of our scattering process. The method outlined here is based on the work mentioned before and extends existing technology. It also shares many similarities with the method developed in refs. [85][86][87][88][89] to expand cross sections around the limit of all radiation being soft. Our expansion is carried out at the integrand level, i.e., before loop or phase space integrals are carried out. The resulting expressions can be interpreted diagrammatically. This in turn greatly simplifies the analytic computation of matrix elements by employing powerful loop integration techniques like the reverse unitarity framework [90][91][92][93][94] or integration-by-part (IBP) identities [95,96]. Our expansion is systematically improvable as we can compute to arbitrarily high power in our expansion parameter. The mathematical functions that appear in each term of the expansion are determined by the first few expansion coefficients.
The collinear expansion of cross sections can find many applications in the computation of higher order corrections to scattering processes. Cross sections for the production of hard probes h can be approximated by performing a systematic collinear expansion. Recently, an all-order factorization theorem was derived for the first order in this collinear expansion [26]. While the usefulness of our expansion technique depends on the specific observable in question, it is obvious that key observables like the rapidity or transverse momentum of a hard probe are amenable to such an expansion. We demonstrate the applicability of our collinear expansion to the rapidity distribution of the Higgs boson produced via gluon fusion. By calculating the collinear expansion to its second order, we demonstrate the excellent convergence of our series towards the full result at NNLO in perturbation theory.
Kinematic limits of cross sections can also be used to identify universal structures of quantum field theories. Our expansion technique allows to gain access to splitting functions or integrated counter terms that may find application in subtraction algorithms used for the computation of fully differential cross sections. Universal building blocks that find their application in the resummation of perturbative cross sections can be accessed efficiently using this expansion technique.
One example of such universal building blocks are so-called beam functions [72,97] which arise in SCET and play a crucial role in factorization theorems of hadronic observables. We demonstrate how to relate beam functions to the kinematic limit of our perturbative cross sections and how they can be extracted efficiently. Specifically, we investigate the transverse momentum (q T ) dependent beam functions and N -jettiness (T N ) beam function. We illustrate our method by computing these quantities through NNLO, up to the second order in the dimensional regularization parameter , confirming recent results in the literature [98][99][100]. These results are necessary input for the calculation of aforementioned beam functions at N 3 LO in QCD, where much progress has been already made for the quark T N beam function [101][102][103], and which has already been achieved for the quark q T beam function and TMDPDF [104]. In our companion papers [105,106], we complete this task by computing the q T and T 0 beam functions in all channels at N 3 LO based on the methods outlined in this article.
In recent years the universal structure of cross sections beyond leading power in kinematic expansions within SCET have been explored . As this avenue of research is still growing rapidly, our expansion techniques may provide analytic information towards the structure of cross sections at higher power. In fact, the method developed in this paper is inspired by the calculation of power corrections in fixed order SCET for T 0 [117,118,129] and q T [130]. It will be interesting to extend these studies to higher order in α s and the power expansion. We hope that our techniques will provide readily accessible tools for the computation of yet unknown universal building blocks.
This article is structured as follows: In section 2 we setup a parameterization for differential cross sections for color singlet production at hadron colliders. This will mainly serve to develop a notation and to identify the objects that we aim to expand. In sec-tion 3 we introduce the general strategy of expanding differential hadronic cross sections around the collinear limit, identifying the relevant kinematic regions and formally defining what we intend by collinear expansion. We then continue the discussion about collinear expansions in section 4 by showing in practice how to perform the collinear expansion for squared matrix elements. We will show explicit examples of the expansion of two loop cut diagrams at leading and beyond leading power, both for real radiation as well as for loop corrections. In section 5 we explain how our collinear expansion of cross section is related to the effective field theory framework of SCET and in particular to the factorization of hadronic differential cross sections. In section 6 we review the role of SCET beam functions in the factorization of hadronic differential cross sections and we show that they are naturally connected to the leading term of our collinear expansion of cross sections. We discuss in detail how to obtain beam functions both in the case of q T and T N . In section 7 we apply our formalism to compute the rapidity spectrum of the Higgs in gluon fusion at NNLO in QCD via the collinear expansion of the partonic cross section. We conclude in section 8. Our two-loop results for the q T and T N bare beam functions are provided as ancillary files.

Setup for differential cross sections
In this section, we develop the notation for differential cross sections at hadron colliders. In section 2.1, we introduce our generic notation for the production of a colorless hard probe h in a proton-proton collision. In section 2.2 we provide a detailed derivation of the required differential phase space.

General setup and notation
We consider the production of a colorless hard probe h and an additional hadronic state X in a proton-proton collision. Examples of such processes are the gluon fusion production cross section of a Higgs boson or the hadronic production of a Z boson or virtual photon (Drell-Yan). (2.1) Here, P 1,2 are the momenta of the incoming protons, which in the hadronic center-of-mass frame are given by where S = (P 1 + P 2 ) 2 is the hadronic center-of-mass energy and the protons are aligned along the directions n µ = (1, 0, 0, 1) ,n µ = (1, 0, 0, −1) .
In eq. (2.1), p h is the momentum of the hard probe h, and k is the total momentum of the hadronic state X, and as indicated both momenta are taken to be incoming.
The hadronic process in eq. (2.1) receives contributions from the partonic processes where i and j are the flavors of the incoming partons, and their momenta are given by such that the partonic center of mass energy is given by In eq. (2.4), X n is a hadronic final state consisting of n ≥ 0 partons with momenta {p 3 , · · · , p n+2 } and total momentum k µ ≡ i>2 p µ i . We are interested in describing processes that are differential in the four momentum p µ h , which we parameterize in terms of its rapidity Y and virtuality Q, and by momentum conservation its transverse momentum p µ h⊥ is fixed to be p µ h⊥ = −k µ ⊥ . The momentum k µ is parameterized in terms of the variables We refer to the hadronic cross section differential in the above variables as the general differential cross section, Here, the sum runs over all possible initial state configurations i, j, the f i (x) denote the parton distribution functions, and dη ij /(dQ 2 dw 1 dw 2 dx) is the general partonic coefficient function. Eq. (2.9) is normalized by σ 0 , which contains all constant factors appearing in the Born level cross section. The Bjorken momentum fractions x 1,2 can be expressed in terms of the variables introduced above.
where the momentum fractions appearing at Born level are given by where τ = Q 2 /S and we use the functions At Born level, k µ = 0, such that the momentum fractions x 1,2 reduce to x B 1,2 , while in the presence of real radiation the kinematic constraint k µ < 0 dictates that x 1,2 ≥ x B 1,2 . The general partonic coefficient function in eq. (2.9) is given by Here, the sum runs over all hadronic final states X n consisting of n partons, and dΦ h+n is the phase space measure of the h + X n final state which will be discussed in more detail in section 2.2. |M ij→h+Xn | 2 is the associated squared matrix element summed over final and initial state colors and helicities. We have also pulled out the overall normalization factor N ij , related to the spins and polarizations of the incoming partons. Depending on the initial state, it is given by (2.14) Here, g, q (q) and q (q ) indicate a gluon, (anti-)quark, and (anti-)quark of different flavor than q, respectively. We expand the general partonic coefficient function in α s as Here, η V ij contains the Born cross section and purely virtual corrections, and can itself be expanded in α s /π with the first term η V (0) ij = δī j for flavour diagonal processes like Drell-Yan or Higgs production. The η ( ,m,n) ij are separately holomorphic in the vicinity of w 1 = 0 or w 2 = 0.
The differential cross section for a specific observable T that only depends on p µ h and k µ is obtained from our general differential cross section given in eq. (2.9) as Here, the convolution integral is defined as (2.17) The corresponding partonic coefficient function differential in T is given by where T (Q, Y, w 1 , w 2 , x) picks out the value of the observable at a given phase space point. Note that in the above equation the variables z i are still functions of w 1 , w 2 and x as specified in eq. (2.12). The partonic coefficient function in (2.16) contains ultraviolet (UV) and infrared (IR) divergences. We regulate such divergences using conventional dimensional regularisation by extending the space time dimension by an infinitesimal amount to be d = 4 − 2 . UV divergences are removed by renormalization in the MS scheme. IR singularities are removed by the standard mass factorization redefinition of the PDFs. Specifically, the unsubtracted PDF f i (x) can is given in terms of the finite PDF in the MS scheme f R i (x) as where the sum runs over all parton flavors j, Γ ij is the PDF counterterm that is known through three loops [140,141], and we suppress the associated factorization scale µ. This allows us the write the hadronic differential cross section of eq. (2.16) in terms of finite quantities, dσ

Differential phase space
To derive the phase space differential in the variables defined in eq. (2.8), we start from the generic expression for the phasespace of the h + X n system, where and k µ = n+2 i=3 p µ i is the total momentum of X n . Next, we separate the integration over p h and k by inserting the unity This splits the h+X n phase space measure into an integral over the phase space Φ m 2 for two massive particles and the phase space Φ 0 n for n massless partons of total invariant mass µ 2 , The two phase space measures are defined as The on-shell constraint for p h is used together with the definition of the rapidity of eq. (2.7) to define the born momentum fractions x B 1,2 . Transforming from k µ to the variables introduced in eq. (2.8), we obtain the desired result for the differential phase space, In the special case of having zero or one final state parton, eq. (2.27) becomes (2.28)

Collinear expansion of color-singlet cross sections
In this section, we introduce the general strategy of expanding cross sections around the collinear limit. We begin by identifying the key kinematic regions in which we want to expand cross sections in section 3.1. Next, we define the collinear expansion of hadronic cross sections in section 3.2. Finally, we comment on the use of different coordinates in performing a collinear expansion in section 3.3. We will provide explicit examples on how to implement this in practice for matrix elements in section 4. In this section, it will be very convenient to work with light-cone coordinates 1 . We decompose a momentum p µ as where the p ± components are explicitly given by and p ⊥ is the remaining transverse component. 1 Note that another popular conventions in the literature defines light-cone coordinates through the

Power counting and modes
Hadronic color singlet cross sections for infrared and collinear safe observables are finite quantities. The perturbative description of such observables becomes inadequate when the value of the observable forces hadronic radiation produced on top of the colorless final state to be in the infrared or collinear regime. For example, it is well known that when imposing an infrared-sensitive measurement T , the cross section receives contributions of up to two logarithms L = ln(Q/T ) per order of the coupling constant, i.e. σ ∼ α n s L 2n . In the limit of T vanishing the presence of such logarithms hence signals the sensitivity to infrared physics and truncated perturbation theory does not yield an accurate description of the observable, as the large logarithms L can overcome the suppression in α s .
In order to understand cross sections in their infrared and collinear kinematic regimes, it is necessary to identify and characterize these regimes. Simply speaking, one can classify two particles with momenta p i and p j as collinear to each other when p i · p j → 0, while a particle with momentum p i is considered soft when p µ i → 0. More precisely, particles should be classified as soft and collinear relative to the hard scale of the process.
In the case of sufficiently inclusive color-singlet processes, there are only two designated directions, namely the lightlike directions n µ andn µ of the incoming protons. Employing the lightcone notation introduced in eq. (3.1), we can thus classify a momentum p µ i = (p + , p − , p ⊥ ) in the different dominant kinematic regions as hard : p µ i ∼ Q (1, 1, 1) , n-collinear : p µ i ∼ Q (λ 2 , 1, λ) , n-collinear : Here, λ 1 is an auxiliary power counting parameter indicating the suppression of the different modes relative to the hard scale Q. Let us discuss eq. (3.4) in more detail: • The hard region describes momenta directly associated with the production of the hard probe h. Since h has invariant mass p 2 h = Q 2 , parametrically hard momenta also have virtuality p 2 i ∼ Q 2 . For example, virtual corrections to the partonic process, i.e. the form factor, are sensitive to this scaling.
• The n-collinear region describes a momentum where n · p i n · p i , and hence p i is aligned with the n-direction. The scaling of the transverse component follows by noting that for on-shell particles • The soft region describes low-energetic, but isotropic radiation, as is manifest from the homogeneous scaling in eq. (3.4). The choice of m in eq. (3.4) depends on whether the observable T under consideration is sensitive to the lightcone momenta only (m = 2) or also to transverse momenta (m = 1). In the SCET literature, these two cases are referred to as ultrasoft and soft, respectively, but in this work we simply refer to both cases as soft.
Note that in more general cases, such as also measuring final-state jets or complicated observables, more modes may arise, and this has given rise to a plethora of "scaling hierarchies" in the literature, see for example refs. [82,[142][143][144][145][146][147][148][149][150][151][152][153][154][155][156][157][158][159]. For sufficiently inclusive observables as considered in this paper it suffices to only consider eq. (3.4). Above, we have only given heuristic arguments for the observation that the scalings in eq. (3.4) are the only relevant ones. In fact, it is precisely the modes stated in eq. (3.4) that arise in proofs of QCD factorization. These proofs are based on the insight that there is a one-to-one correspondence between the momentum regions giving rise to the large Q behavior and mass divergences in massless perturbation theory [160,161]. These mass divergences arise at pinch-singular surfaces, i.e. momentum regions when loop momenta can not be deformed away from singularities in the appearing propagators, and thus correspond to classically allowed scatterings. The position of these singular surfaces can be determined using the Landau criteria [162], and one can then derive power counting rules for these singular surfaces to approximate amplitudes in the vicinity of the surface, and these rules precisely lead to the momentum regions shown in eq. (3.4). For a more comprehensive discussion of these proofs, we refer to refs. [24,163,164]. Similarly, these singular surfaces also arise when analysing Feynman diagrams using the method of regions [59]. They are also the basis for the formulation of SCET [67][68][69][70][71], an effective field theory to describe QCD in the infrared limit that separates quark and gluon fields into modes corresponding to eq. (3.4). We will discuss this connection in more detail in section 5.

Expanding cross sections around the collinear limit
Having characterized the relevant kinematic regions for infrared-sensitive observables, we now discuss the expansion of hadronic cross sections of eq. (2.16) around the particular limit where all final state radiation becomes collinear to one of the incoming proton momenta.
Let us define our collinear expansion: we want to expand around the limit where all real momenta are treated as n-collinear, and thus the total momentum k µ of the hadronic final-state is n-collinear as well, i.e. it it scales as We then want to expand the hadronic differential cross section in eq. (2.16) to obtain a power series in λ, Here, the leading-power cross section σ (0) scales as λ −2 , 2 the next-to-leading power (NLP) cross section 3 σ (1) as λ −1 , and so forth. Depending on the observable T , this series may start at higher orders in λ, but for the infrared-sensitive observables discussed in this paper we always encounter a leading O(λ −2 ) term.
It is desirable that Born quantities like Q 2 and Y are unaffected by the expansion we want to carry out, as they set the hard scales of the process. The importance of this for expansions at subleading power in λ was already stressed in refs. [129,130]. As a consequence, the Bjorken momentum fractions given in eq. (2.10) need to be expanded. Expressing them in terms of hard quantities and the momentum k we find Since the momentum fractions enter as arguments of the PDFs, a pure hadronic expansion to higher orders in λ will automatically involve derivatives of PDFs, as firstly noted for T 0 in ref. [117]. Furthermore, the variables w 1 and w 2 we introduced in eq. (2.8) must also be expanded, where x B 1,2 are the momentum fractions at Born level, see eq. (2.11). As a consequence we find that the n-collinear limit of eq. (2.18) becomes lim n−coll.
where w 1,2 are evaluated according to eq. (3.8). The definition of our observable T itself may not be invariant under rescaling according to our power counting. In order to achieve a pure expansion of the hadronic cross section we may either expand the observable constraint or solve the constraint using one of the remaining integration variables and expand subsequently. We address how the general partonic coefficient function Constructing a collinear expansion can be done with different objectives in mind. One objective can be to obtain a pure series expansion of the hadronic cross section as discussed above. Another objective can be to simplify the computation of the partonic coefficient function which does not require a pure expansion of the hadronic cross section. In the latter scenario one would only expand the partonic coefficient function η ij on the right-hand side of eq. (2.18), but not expand the w 1,2 and the momentum fractions x 1,2 as presented above. This approach can also serve as a suitable proxy to a collinear expansion, where parts of the cross section are kept exact.

Expansions using different coordinates
So far, we defined our power counting such that the invariant mass Q 2 and rapidity Y of the produced hard probe h scale homogeneously as O(λ 0 ) and the lightcone components of the total momentum k µ of the hadronic final state have a homogeneous power counting in λ. This is reasonable, since one can only measure directly the final-state particles in the hadronic collision, which are then used to define the power counting. In particular, Q 2 and Y are the only hard scales in the considered hadronic cross section dσ/(dQ 2 dY dT ).
This setup immediately implies that the momenta p 1 and p 2 of the incoming partons do not have a homogeneous power counting, as it is evident from their explicit expressions in terms of Q 2 , Y and k µ , see eq. (2.10). Thus, p 1 and p 2 give rise to an infinite tower of power corrections in λ, which in turn requires an expansion of w 1 and w 2 used to define the general partonic coefficient function, as shown in eq. (3.8).
Since one has access to all incoming and outgoing momenta in the calculation of the partonic coefficient function, the collinear expansion can be also be defined by assigning a homogeneous power counting to p 1 , p 2 and k. Since this assignment is only meaningful for the partonic process, we refer to it as partonic collinear expansion. In this approach, the rescaling appropriate for the collinear limit is given by The key advantage of this assignment is that w 1 and w 2 , which are the natural variables to express the partonic partonic coefficient function, now have homogeneous power counting and do not need to be re-expanded themselves. Thus, the collinear expansion has been reduced to an expansion in w 2 . The drawback of the partonic collinear expansion is that the rapidity Y of hard probe h no longer uniformly scales as O(λ 0 ), which is evident from the expression which now is a quantity derived from p − 1 and p + 2 , rather than fixing these as in eq. (3.10). Comparing the rescalings in eqs. (3.8) and (3.11), we see that both approaches agree at leading power, but differ at subleading power. Since there is a well-defined relation between the two approaches, one can easily obtain one expansion from the other, but care has to be taken to consistently apply the power expansion.
In practice, each choice of defining the expansion has its advantages and disadvantages. We can discuss these by classifying the expansions according to the choice of independent variables used to express the partonic coefficient function, which by Lorentz invariance only requires four independent invariables. It is useful to summarize the above observations for the following possibilities: This parameterization has the advantage that is entirely expressed in terms of information about the final state momenta, including the Born measurements Q and Y of the hard probe h. As properties of the final state, Q and Y need not be expanded, and the collinear expansion is a strict expansion in k µ . Since partonic matrix elements are typically more concise when expressed in terms of the incoming momenta and the final-state radiation, the main drawback of this expansion is that it leads to lengthier expressions for the expanded matrix element. Furthermore, measuring the rapidity Y fixes a reference frame for all momenta, such that boost invariance is not manifest anymore.
• (Q 2 , w 1 , w 2 , x): Since w 1 and w 2 are defined as ratios of Lorentz scalars, boost invariance is manifest in this parameterization. Its disadvantage is that w 1 and w 2 do not have manifest power counting in terms of the observables Q 2 and Y , and instead must be expanded in k µ according to their definitions in eq. (2.8). Alternatively, one can assign homogeneous power counting to w 1 and w 2 using eq. (3.11), which then requires to expand the rapidity Y in λ.
Here, we trade Q 2 and Y for the lightcone momenta (p − 1 , p + 2 ) of the incoming partons. This parameterization has the advantage of expressing everything in terms of the momenta of massless particles, i.e. the incoming momenta and the hadronic radiation. A disadvantage of this parameterization is that p − 1 and p + 2 do not have manifest power counting in terms of hadronic variables Q 2 and Y , and thus must be expanded in λ.
These parameterizations are of course equivalent, and in practice the preferred parameterization depends on the intended application. While the general illustration of the power expansion is made most manifest using (Q 2 , Y, k + , k − , x), expanding the partonic cross sections is simplified using (Q 2 , w 1 , w 2 , x). Finally, we give the explicit relation between the different parameterizations. We can change variables from (w 1 , w 2 ) to (k + , k − ) using where the required variable transformations are given by and for brevity we keep implicit that k µ is parameterized in terms of (k + , k − , x). Note that since k + and k − are defined in the hadronic center-of-mass frame, manifest boostinvariance is lost, and thus eq. (3.13) becomes explicitly Y -dependent. Eq. (3.13) makes it clear that defining the power counting in terms of Q 2 and k requires a expansion of w 1 and w 2 on the right hand side.
We can further change variables from (Q 2 , Y ) to (p − 1 , p + 2 ), dη ij dp − 1 dp , (3.15) where the required variable transformations are given by Here, fixing the power counting of p − 1 , p + 2 and k requires to expand Q 2 and Y accordingly.

Collinear expansion of matrix elements
In this section we show how the technique of collinear expansions developed in the previous section is applied in practice. To setup our conventions for this section, we first discuss the phase space volume for producing the hard probe h with additional emissions in section 4.1, before illustrating the collinear expansion of matrix elements explicitly for both real radiation in section 4.2 and for loop integrals in section 4.3.
Throughout this section, we will consider the scenario where k µ is collinear to the incoming parton with momentum p µ 1 = p − 1 n µ /2. According to eq. (3.4), this implies that we assign the following scaling to k: In order to obtain a strict power series expansion of the hadronic cross section it is necessary to expand the partonic momentum components p − 1 and p + 2 around the collinear limit. For the purpose of this section we instead perform a partonic collinear expansion (see section 3.3), treating p 1,2 as external variables and thus as O(λ 0 ) quantities. All final results are functions of k − /p − 1 and k + /p + 2 and one can straightforwardly recover a pure expansion in terms of hadronic observables following section 3.3.

Collinear phase space
The phase space volume for producing the hard probe h with two emissions, as defined in eq. (2.25), is given by It follows immediately that the phase space volume in eq. (4.2) scales as Φ h+2 ∼ λ 2−4 .
As usual, we take all momenta to be incoming, and denote the total momentum of all outgoing partons by k = p 3 + p 4 . In the above diagram and those below, the dashed line indicates the on-shell constraints of the final state particles, with the solid lines representing massless partons and the double line representing the heavy color-singlet state h.
The scaling of Φ h+2 ∼ λ 2−4 can be easily deduced without calculating the actual phase space integrals. Since k is treated collinear to p 1 , both final-state momenta p 3 and p 4 must be collinear to p 1 as well. The associated integration measures and δ functions entering eqs. (2.25) and (2.26) transform as As a consequence, the double-real phase space measure scales as which is precisely the scaling observed in eq. (4.2). Similarly, it follows that the more general case of the h + n real emission phase space has the scaling (4.5)

Collinear limit of real radiation
We consider the following example of a more complicated, purely real Feynman integral, (4.6) Let us first consider the case where both p 3 and p 4 are collinear to p 2 . In this scenario, since both propagators in eq. (4.6) only involve collinear momenta, and thus scale homogeneously as λ 2 under then-collinear rescaling of eq. (3.4) and no expansion of eq. (4.6) in λ is needed.
In contrast, if we consider p 3 and p 4 to be collinear to p 1 , then the second propagator is not homogeneous in λ anymore, as it contains both n-collinear andn-collinear momenta. To expand the propagator in this limit, we apply the n-collinear rescaling of eq. (3.4) to p 3 and p 4 , With these rescalings, it is now straightforward to expand the second propagator in eq. (4.6) in λ, For n = 0, this propagator is linear in the real momenta p 3 and p 4 , and thus corresponds to an eikonal propagator. Higher orders in λ only involve pure powers of the eikonal propagator, thus yielding a relatively simple structure of the expansion. Together with eq. (4.4), the integral in eq. (4.6) can thus be expanded as This expansion can be represented diagrammatically as Here, the dotted line indicates the expanded (eikonal) propagator and the dot on the line represents higher powers of this propagator. The label denotes the additional kinematic factor arising from eq. (4.8).
The expansion in eq. (4.9) results in several advantages. First, we observe that each term in the expansion is homogeneous under the n-collinear rescaling transformation in eq. (4.7). As a consequence, we may directly determine the functional dependence of each term in the expansion on k + similarly to the case of the phase space volume. In other words, the resulting functions will be simpler since they only depend on k + via a multiplicative pre-factor. Second, the structure of expanded Feynman integrals is amenable to IBP reduction techniques via the framework of reverse unitarity [90][91][92][93][94]. The benefit is that the appearing integrals can be related to so called master integrals. In our example we find the IBP relations (4.11) Clearly, it is very advantageous that any higher order term in our expansion is related to the same master integrals as the first, which in our example is just the phase space volume. The unexpanded integral of our example in eq. (4.6) is itself related to the phase space volume by an IBP identity, From this we can easily see that the coefficients obtained in eq. (4.11) and eq. (4.12) correspond exactly to the coefficients of the expansion of the exact result. In summary, we outlined a procedure that allows us to perform an expansion of real radiation integrals around the limit of all final state partons becoming collinear to an initial state momentum. This expansion is carried out by simply performing the appropriate collinear rescaling transformation of eq. (3.4) on all final state parton momenta and subsequently expanding the integrand of our real radiation integral in the artificial parameter λ, prior to actually evaluating the integral. Each term in the expansion in λ then corresponds to exactly one term in the expansion of the integral in k + . The computation of the terms in the expansion is greatly facilitated by applying techniques like IBP identities via the reverse unitarity framework. In particular, any term appearing at higher orders in the expansion will be expressible in terms of master integrals that appear already in the first few terms.

Expansion of loop integrals
In contrast to the phase space integral over real momenta considered in section 4.2, where the requirement of k being collinear to p 1 restricted p 3,4 to be collinear to p 1 as well, such a restriction does not appear for loop momenta. Despite this, it is still useful to expand loop integrals in a similar fashion around the hard, collinear and soft regions. As discussed in section 3.1, for factorization proofs this is crucial to separate these different regions into distinct matrix element, while in the method-of-regions approach of ref. [59] it used to simplify loop integrals by expanding the integrand in all relevant limits and combining their individual results.
Here, we will show for a simple example how one can easily approximate and expand loop integrals in the discussed regimes, and that the sum of all regions indeed reproduces the full result. This will be illustrated using the following real-virtual diagram, (4.14) Here, the total final-state hadronic momentum is k = p 3 , and thus the p 3 integral is actually fixed. In eq. (4.14), (a) n = Γ(a + n)/Γ(a) is the (rising) Pochhammer symbol, and we abbreviate common loop factors by As before, we consider the limit where k is collinear to p 1 , such that k + ∼ O(λ 2 ) and k − ∼ O(λ 0 ). This immediately implies that eq. (4.14) scales as (4.16) The second line has homogeneous scaling in λ −2−4 , and is the dominant contribution in the limit λ → 0. We will see below that this result is entirely from the region where the loop momentum is collinear to p 1 . In other words, the leading-power limit of I RV arises from the region where both loop and real momenta are collinear to p 1 . The first line in eq. (4.16) does not scale homogeneously in λ, but is suppressed at least as O(λ 2 ) compared to the leading-power limit. We will see that this line entirely arises from the region where the loop momentum is hard. In particular, the two contributions have different fractional scalings in λ −2 and λ −4 , respectively. These scalings arise entirely from the loop integral measures, and thus can be easily distinguished between the different contributions.

Collinear limit
We first consider the loop momentum p 4 to be collinear to the incoming parton with momentum p 1 . According to eq. (3.4), we hence transform The first three propagators in eq. (4.14) scale homogeneously as O(λ −2 ) under this rescaling, while the last propagator in eq. (4.14) is not homogeneous in λ and must be expanded, Together with eq. (4.5), this allows us to expand the integrand in eq. (4.14) as The overall scaling in λ −4 arises from d d p 4 ∼ λ 4−2 and dΦ h+1 ∼ λ −2 , and thus is independent of the structure of the diagram itself. Each order of the expanded integrand is now homogeneous in λ. The expansion in eq. (4.19) can be illustrated graphically as Here, the dotted lines are linear (eikonal) propagators, and the dot on the line denotes that the propagator is raised to one power. Note that in the second diagram, the explicit 1/p 2 4 propagator is canceled, indicated by the contracted vertex. Being able to represent collinear expanded diagrams again in a diagrammatic fashion is extremely useful. In particular, the structure observed in the collinear loop expansion makes it possible to use IBP techniques for the computation of loop and phase space integrals.
The leading-power integral in eq. (4.19) can be evaluated as and thus correctly reproduces the last line of eq. (4.16). Note that the higher-order terms in λ, such as the second integral in eq. (4.19), can be shown to vanish identically in dimensional regularization.

Hard limit
The hard region is characterised by treating the loop momentum as uniformly larger than our expansion parameter λ, while the final state momentum p 3 is still treated as collinear to p 1 . Only one propagator in eq. (4.14) involves p 3 , and can be expanded in λ as .

(4.22)
All other propagators in eq. (4.14) scale as O(λ 0 ) and are not expanded. Together with the rescaling of the phase space measure according to eq. (4.5), the leading-power hard limit of eq. (4.14) becomes (4.23) The overall scaling in λ −2 arises entirely from the phase space measure, as the hard loop measure scales as d d p 4 ∼ λ 0 . This shows that hard and collinear loops never have the same dependence on , and thus can be easily distinguished by their overall scalings. Despite the modified propagator in this integral, it can still be subjected to the usual loop integration techniques like IBPs and differential equations. The same holds true for all higher order terms in the expansion of the full loop integral. The explicit example in eq. (4.23) is also easily performed using Feynman parameters. We obtain This result exactly agrees with the infinite sum over n in eq. (4.16) evaluated for m = 0, i.e. the O(λ −2 ) to eq. (4.16). Furthermore, every higher-order term in the expansion in k + of the second to last line of eq. (4.14) corresponds to exactly one term in the integrand expansion of I RV in λ. Terms proportional to odd powers of λ drop out identically. Since higher order terms in the expansion essential just modify the powers of the propagators at the integrand level according to eq. (4.22) it is particularly convenient to use IBP techniques in such a computation.

Anticollinear limit
For completeness, we also consider the limit where p 4 is collinear to the incoming parton with momentum p 2 , in contrast to k which is chosen collinear to p 1 . According to eq. (3.4), we hence transform With this rescaling, we need to expand two propagators of the integrand in eq. (4.14), For brevity, we only show the two leading terms each. Together with eq. (4.5), this allows us to expand the integrand in eq. (4.14) as which is scaleless and thus vanishes in pure dimensional regularization. Note that the integral from expanding the propagators through O(λ n ) scales as λ n−4 . Since the only term with this dependence in eq. (4.16) is fully given by the n-collinear limit of eq. (4.21), the p 2 -collinear limit must in fact vanish to all orders in λ.

Soft limit
We want to compare the result of the collinear expansion to a soft expansion of Feynman diagrams. To obtain the purely soft region, we rescale the loop momentum p 4 in eq. (4.14) as To obtain the soft-collinear overlap, we first rescale p 4 as collinear, followed by a subsequent soft rescaling, Let us explicitly discuss the transformation of two of the propagators in eq. (4.14) under eq. (4.29), . (4.30) In the first case, rescaling the collinear limit as soft only amounts to an overall rescaling by λ −2 , but does not change the relative scaling of the two terms in the propagator. In the second case, we observe both that only one term in the denominator gets rescaled in the soft limit, and thus one will encounter a different kinematic structure when expanding this propagator in λ than in the collinear limit. However, in both cases shown in eq. (4.30), it is easy to see that the soft-collinear limit is identical to taking the soft limit directly. The same holds for the two other propagators in eq. (4.14) that are not explicitly shown here.
In conclusion, we find that at the diagram level, the soft-collinear overlap is equal to soft limit itself. Finally, we note that the leading-power soft limit of eq. (4.14) is given by

(4.31)
This integral is scaleless and vanishes in dimensional regularization.

Discussion
To summarize the key results of this section, we have shown how Feynman diagrams can be systematically expanded in their collinear limit by assigning the appropriate scalings to all loop and real momenta, which allows one to expand the integrand in λ. In particular, the expanded integrands allow for a diagrammatic representation and are amenable to standard integral techniques such as IBP [95,96] relations or the method of differential equations [165][166][167][168][169]. This significantly simplifies evaluating the expanded integrals compared to the exact integral, and thus provides a convenient strategy to approximate Feynman diagrams in the collinear limit. The illustrated methods are conceptually very simple, and thus easily extend to more complicated diagrams with additional external partons or multi-loop integrals.
In the case of real radiation, the requirement that the total real momentum k µ is collinear implies that all real momenta are collinear individually. This does not apply for loop momenta, which are not confined to be in the collinear region. As a consequence we need to follow the method of regions [59] and compute the regions where the loop momenta are hard and where they are collinear. The sum of both regions yields the correct expansion of our Feynman integrals. The results of the different regions give rise to different scalings as λ −4 and λ −2 , respectively. This difference is entirely due to the loop measure, and thus hard and collinear contributions can be easily identified by their scaling exponent. In other words, since the expansion of the loop integrand itself is a simple Laurent series in λ, the loop measure fully determines the non-integer powers of λ −n .
We also discussed the soft limit of matrix elements. We found in an explicit example that the soft region of a loop integral can be obtained by first computing the collinear region of this integral and subsequently taking the soft limit. As a matter of fact this property holds more generally. The soft-collinear overlap of a partonic coefficient function can either be computed by first performing the collinear expansion and then the soft expansion, or vice versa. More precisely, expanding the first n terms in the collinear expansion around the production threshold up to n terms will correctly reproduce the n th power in the threshold expansion. This provides a stringent test of the collinear expansion by comparing to existing analytic results for which a threshold expansion was performed. It also provides a considerable simplification for the calculation of collinear master integrals. For example, if the method of differential equations is utilized to compute master integrals for the collinear expansion, then the boundary conditions for these differential equations can be chosen to be the threshold-expanded integrals. For the computation of threshold expanded integrals see for example refs. [85,86,170].

Kinematic expansions and SCET
Differential cross sections may be kinematically enhanced in all different momentum regions shown in eq. (3.4). Above we only discussed the expansion cross sections around one particular limit, namely the collinear limit. However, in order to perform a physically sensible and consistent expansion of a hadronic cross section we need to expand in observable quantities. A collinear expansion of a hadronic cross section alone typically does not satisfy this requirement. In order to obtain a physical expansion in an observable all momentum regions where the observable is kinematically enhanced must be considered. Depending on the observable of interest, the necessary ingredients to achieve this goal may vary.
Soft-Collinear Effective Theory (SCET) [67][68][69][70][71] provides an excellent tool to organize the expansion in such kinematic limits, and we discuss in section 5.1 how the tools developed in the previous section connect to factorization theorems derived in SCET. In such factorization theorems, it is crucial to account for the overlap of regions when combining multiple kinematic expansions, which we address in section 5.2.

Kinematic expansions and factorisation theorems
The momentum regions shown in eq. (3.4) are precisely the basis for the formulation of SCET, which is an effective field theory describing QCD in its collinear and soft limits, i.e. the leading infrared region. Schematically, the SCET Lagrangian is expanded as Here, the superscript (0) indicates the leading-power (LP) terms in the expansion in λ 1, where as before λ is an auxiliary power counting parameter. The L (k) indicate subleading power Lagrangians [107-111, 120-123, 126, 127] that are suppressed by λ k w.r.t. to the leading power. The leading-power SCET Lagrangian can be organized as Here, L h contains the hard scattering operators mediating the underlying hard interaction, and L (0) n,n,s are the SCET Lagrangians for n-collinear,n-collinear and soft fields as defined by eq. (3.4), respectively. 4 More generally, in the presence of multiple collinear directions as required e.g. for multijet processes, eq. (5.2) contains a sum over all relevant collinear directions {n i }. SCET also allows for a treatment of Glauber modes, which appear as non-local potentials in L (0) G , the leading power Glauber Lagrangian [172]. In SCET, factorization is achieved by a field redefinition of soft and collinear fields which decouples the soft and collinear Lagrangians from each other [70]. These modes can still interact with each other through the Glauber Lagrangian L (0) G , which thus can break factorization. In this work we will consider observables where the Glauber contributions from L (0) G either cancel identically [19-21, 23, 24, 164, 173] or start contributing to higher perturbative orders that the one we consider in this work [174,175].
In SCET, the leading kinematic regions are made manifest and decoupled from each other at the Lagrangian level, and thus greatly simplifies the derivation of factorization formulas. For suitable factorizable infrared-sensitive observables T , which we take to vanish as T → 0 in the Born limit, the hadronic cross section eq. (2.20) can be factorized as [20,72] dσ As usual, Q and Y are the invariant mass and rapidity of the colorless final state. The sum runs over all flavor combinations (i, j) contributing at Born level, ij → h, and σ 0 is the corresponding Born partonic cross section. 5 The hard function H ij encodes virtual corrections to the Born process ij → h, i.e. it is given as the corresponding renormalized form factor. The beam functions B i (x, T ) encode the probability to extract a parton of type i with momentum fraction x from the proton, together with the contribution from collinear radiation to the observable T , while the soft function S(T ) encodes the effect of soft exchange between the protons. Since S(T ) only differs between quark-and gluoninduced processes, we suppress an explicit flavor label. Finally, ⊗ denotes a convolution in T , whose precise structure depends on the chosen observable T , and often can be made multiplicative in a suitable conjugate space. Note that in eq. (5.3) we suppress explicit renormalization scales that are present in all functions. The factorisation of degrees of freedom at the Lagrangian level makes the ingredients for the various functions in eq. (5.3) evident. The hard function, n-collinear andn−collinear beam functions and the soft function are each defined in terms of only hard, n-collinear, n-collinear and soft degrees of freedom, respectively. This implies that the expansion techniques developed in this article are perfectly suited to determine beam functions from a perturbative computation using a pure collinear expansion of cross sections. Here, it is important that both real and loop momenta are expanded as collinear. We will provide explicit examples by obtaining the NNLO beam functions for T = q T and T = T N (Njettiness) in section 6. We also note that in a similar fashion, one can also obtain the soft function by considering a purely soft expansion.
We also stress that since SCET is an effective field theory, it can be systematically extended by including the power-suppressed Lagrangians L (k>0) in eq. (5.1). This is the EFT analog of expanding cross sections to subleading order in λ about the soft and collinear limits. However, at subleading powers, collinear and soft interactions do not simply factorize similar to eq. (5.3) anymore, and factorization theorems and the resummation of large logarithms become much more involved [112][113][114][115]. Since our expansion technique allows us to perform collinear expansions of partonic cross sections to arbitrary order in λ, we hope that it will also provide insights into the structure of factorisation theorems beyond the leading power, and that it can be used to determine universal quantities like generalizations of soft and beam functions at subleading power.

Soft-collinear overlap and zero-bin subtractions
In order to obtain a full description of a cross section in its infrared limit, we need to combine all collinear and soft regions. Schematically, we expand where the σ (n,n,s) correspond to the expansion of the cross sections where all emissions are treated as n-collinear,n-collinear and soft, respectively. The ellipses denote mixings of these cases, as well as power-suppressed corrections. Note that here in the following, we do not consider the hard region. While it is required to obtain an infrared-finite cross section, it corresponds to physics at the hard scale µ 2 ∼ Q 2 , and does not affect the soft-collinear overlap discussed in the following. In practice, eq. (5.4) is often too naive, as there is a nontrivial overlap between the collinear and soft regions. This arises because the soft limit of a squared matrix element is equal to the soft limit of the collinear limit of the same matrix element. As discussed and illustrated in more detail in section 4.3.4, this can be understood since the soft limit can be equivalently obtained by either directly rescaling 5) or by first rescaling into the collinear limit with a subsequent soft rescaling, Since the second rescaling only lowers the scaling of each component, no information is lost, and eqs. (5.5) and (5.6) produce the same expansion of a matrix element.
Consequently, when one integrates over a collinearly-rescaled momentum, the integral will always contain contributions from the soft region. Schematically, if we write for the collinear expansion f (n) of an arbitrary integrand f , then clearly the integration range extends into a region where the assumed collinear scaling is not justified. In particular, the p − integral extends to p − → 0, which corresponds to a soft region. This contribution to the soft region can be identified and extracted by further expanding f (n) as indicated in eq. (5.6).
In conclusion, the collinear limit of the cross section has an overlap with the soft limit, which can be extracted by an additional reexpansion in the soft limit, which has been demonstrated explicitly for a mixed real-virtual integral in section 4.3.4. We thus need to modify eq. (5.4) as where the soft limit of the collinear cross sections are denoted by σ (n→s) and σ (n→s) , respectively. The terms in brackets hence correspond to the true n-andn-collinear limits of the cross section. Note that in general, σ (n→s) = σ (s) , because the observable T itself has to be expanded in the collinear and soft limits. Let us connect these observations to the corresponding treatment in SCET. As a modal EFT, SCET is built to separately describe soft and collinear modes, and hence as a matter of principle collinear momenta are not allowed to overlap with the soft sector. In practice, it is not feasible to introduce a cutoff between soft and collinear modes. Instead, one follows the same strategy outlined above: after calculating a collinear integral, one subtracts its soft limit to obtain the pure collinear result. This procedure is referred to as zero-bin subtraction [176], and is crucial to a well-defined separation of modes in the EFT. In practice, the zero-bin subtractions are often absent in dimensional regularization or equal to the soft function itself, and thus can be easily taken into account.

Beam functions from the collinear limit
In this section we show how we the collinear expansions can be used to compute beam functions. We briefly introduce the notion of beam functions in section 6.1, and then show in section 6.2 how they are related to the collinear expansion of cross sections developed before. Our method is briefly contrasted with other methods of calculating beam functions in section 6.3. We show explicitly how to obtain the N -jettiness and the q T beam functions at NNLO using this method in section 6.4 and section 6.5, respectively.

Beam functions
Beam functions are defined as gauge-invariant hadronic matrix element that measure the large lightcone momentum entering the hard interaction, as well as the contribution to the observable T from collinear radiation. For example, the quark beam function B q is defined in SCET as [72] Here, the χ n = W † n q are collinear quark fields defined in SCET as quark fields dressed with collinear Wilson lines W n , p n (P ) is a proton moving along the n-direction with momentum P , andn · P is the SCET momentum operator that determines the lightcone momentum of all fields to its right. By boost invariance, the beam function only depends on the momentum fraction x = p − /P − . Similarly,T is the appropriate measurement operator determining the observable T in terms of all momenta of the fields to its right.
Comparing eqs. (6.1) and (6.2), it is evident that the beam function extends the PDF by measuring an additional observable T on top of the longitudinal momentum fraction carried by the struck parton. Both beam functions and PDFs are in general intrinsically nonperturbative matrix elements. For perturbative T Λ QCD , one can perform an operator product expansion of the beam function onto the PDF [72], Here, the only nonperturbative input is given in terms of the PDFs, while the matching kernel I ij are perturbatively calculable. For completeness, we remark that PDFs and beam functions can also be defined without invoking SCET by expressing the collinear quark fields χ n in terms of standard quark fields and collinear Wilson lines W n , which are defined as path-ordered exponentials of the gluon field projected onto the appropriate collinear direction. Beam functions are also often written as the Fourier transform of a position-space correlator, where the separation between the quark fields corresponds to the exchanged momentum and often avoids the need for the momentum operator P in eq. (6.1). PDFs and TMDPDFs were originally defined in this way [20,177], and the equivalence of both formulations was also discussed in the context of T N beam functions in refs. [72,97]. Note that the study of parton distributions from lattice QCD requires the definition in position space, see e.g. refs. [178][179][180][181][182][183][184][185][186] for recent progress towards calculating TMDPDFs on lattice, and refs. [187,188] for a more general overview of parton physics from lattice QCD. For perturbative calculations, both formulations are equivalent.

General strategy
In section 5.1, we discussed that the hard, beam and soft functions in the factorized cross section in eq. (5.3) are each defined only in terms of the hard, collinear and soft modes of eq. (3.4), respectively. Hence, in the limit where all loop and final-state momenta are treated as n-collinear, the hard function, then-collinear beam function, and the soft function only contribute at their respective tree level, where they are normalized to unity and flavor diagonal. Thus, the strict n-collinear limit of eq. (5.3) is given by where we remind the reader that the flavor sum runs over all flavors contributing at Born level, ij → h, and σ 0 is the associated Born partonic cross section. We remark that eq. (6.4) is to be understood at the bare level, as it for example does not encode scale independence. Indeed, as we will see, even after UV renormalization and IR subtraction one encounters leftover poles in , which in the full factorized cross section in eq. (5.3) cancel with the other ingredients.
In the following, we assume that the Born process is diagonal in flavor, i.e. only the gg channel (as in Higgs production in gluon fusion) or the qq,qq channels (as in Drell-Yan or bb initiated Higgs production) contribute, where q is an arbitrary quark flavor. With this assumption, we can fix j =ī in eq. (6.4), which allows us to easily read off the bare beam function by comparing with the n-collinear limit of the cross section given in eqs. (2.16) and (2.18), By fixing the flavor of then-collinear parton asī, we extract the correct beam function for the flavor i in a flavor-diagonal process. Eq. (6.5) has precisely the structure of eq. (6.3), and we can immediately read off the bare matching kernel as the collinear limit of the partonic coefficient function, We stress that the partonic coefficient function here is limited to strictly collinear modes only. In contrast, in the collinear expansion for cross sections discussed before, we also included non-collinear modes when computing loop integrals. However, we showed that the collinear and non-collinear modes can easily be separated by looking at their respective generalized scaling behaviour. Extracting the required parts is consequently easy. In the strictly collinear limit the general partonic coefficient function of eq. (2.15) becomes lim strict n−coll.
dη jī dQ 2 dw 1 dw 2 dx (6.7) The strict collinear limit for the partonic coefficient function for the observable T in eqs. (6.5) and (6.6) is then obtained in analogy to eq. (3.9). A special case of eq. (6.6) is the bare matching kernel differential in w 1 , w 2 and x itself, from which one can project out all other beam functions we interested in. In fact, this double-differential beam function can be related to the fully unintegrated parton distribution first formulated in refs. [189,190] and within SCET in refs. [191,192], where it is also known as double-differential beam function (dBF). Importantly, in general the projection of (w 1 , w 2 , x) onto the desire observable T only holds at the bare level, not after renormalization of the dBF [192]. The renormalization of the dBF is also significantly more complicated than that of the T N and q T beam functions we are interested in, see refs. [193,194] for explicit results at NNLO.
I bare ij still contains infared poles that cancel upon PDF renormalization in eq. (6.5). Even after α s renormalization, this still leaves divergences that cancel in the cross section when combining the n-collinear limit with then-collinear and soft limits, but are remnant in the bare matching kernel. In the EFT, these divergences are of ultraviolet origin and thus can be absorbed in the standard fashion through a counterterm. Subtracting both IR and UV poles in this manner, we obtain the renormalized matching kernel as where ⊗ T denotes the appropriate convolution in T . According to eqs. (6.6) and (6.8), we can obtain the beam function matching kernel as follows: 1. Obtain the bare kernel I bare ij from the strict collinear limit of the partonic cross section.
2. Apply α s renormalization throughẐ αs , which renormalizes the bare coupling constant α b s in the MS scheme.
3. Subtract the EFT UV divergences with the beam-function counterterm Z i B . This renormalization does not change the parton flavor i, and only differs between quark and gluons, but is independent of the quark flavor. In general, this counterterm enters through a convolution in T , which can be trivialized by going to suitable conjugate space.
4. Subtract IR divergences by convolving with the PDF counterterm Γ jj , which as usual mixes parton flavors.
Since the Γ jj and Z i B commute, one can freely rearrange their order in eq. (6.8). Since the beam function counter term Z i B gives rise to the renormalization group equation of the beam function, in practice one can either predict Z i B from the RGEs and check that this cancels all poles in , or determine Z i B by absorbing all poles remaining after QCD UV and IR subtraction and verify that it reproduces the RGE dictated by the EFT. For the T N and q T beam functions, this is discussed in more detail in our companion papers [105,106].
In eq. (6.6), we assumed that the partonic coefficient function is taken in the strict n-collinear limit. As discussed in section 5.2, for certain observables there can be overlap with the soft limit, which in the factorized cross section in eq. (5.3) is already accounted for by the soft function. In such instances, one has to subtract off the soft-collinear overlap, The second term in the above equation denotes that the collinear limit is further reexpanded in the soft limit.

Comparison to alternative methods
Before illustrating our method for the T N and q T beam functions in sections 6.4 and 6.5, we briefly contrast our approach to methods previously used in the literature. Here, we focus on how to calculate the bare beam function, since the renormalization and subtraction of UV and IR divergences always proceeds in the same fashion. Most calculations of beam functions explicitly calculate matching coefficients from matrix element of the beam function operator, see e.g. refs. [97-99, 104, 193-201]. Let us explain some features of this approach for the concrete example of a quark beam function as shown in eq. (6.1), whose bare matching kernel I qj is obtained by evaluating the matrix element in eq. (6.1) with an external on-shell parton of flavor j. In ref. [97], the analytic structure of these matrix elements was discussed in detail for the T N beam function, and it was shown that one can calculate it by taking the discontinuity of matrix elements of the time-ordered operator.
Firstly, this implies that the beam function can be calculated using SCET Feynman rules. Since a single collinear sector in SCET is equal to a boosted copy of QCD, one can equivalently employ QCD Feynman rules. In this case, eikonal vertices arise from the Feynman rules of the Wilson lines W n that are part of the collinear quark fields χ n = W † n q. These can be avoided in lightcone gauge, wheren · A = 0 such that W n = 1, but similar terms arise from the gluon propagator in ligthcone gauge. Since the beam function is defined as a gauge-invariant matrix element, both approaches yield equal results.
The discontinuity can be obtained by using the Cutkosky rules [202] (see also ref. [203]), which corresponds to taking particles exchanged between the quark fields in eq. (6.1) onshell. This is analogous to our approach, where we explicitly consider on-shell radiation into the final state. Alternatively, one can not apply an on-shell constraint and integrate over all particles, and explicitly take the discontinuity afterwards. Both approaches are discussed in more detail in ref. [196], where they are referred to as on-shell and dispersive method, respectively.
An alternative method that does not directly rely on the definition of the beam function in SCET was pointed out in ref. [204], where it was shown that one can equivalently calculate the beam function from phase-space integrals over QCD splitting functions. This approach was used in refs. [100][101][102][103], where the required splitting function at N 3 LO was obtained following the method of ref. [205]. This approach requires to use a physical gauge where gluons are explicitly transverse, for example the lightcone gaugen · A = 0.
Similar to ref. [204], our method does not rely on directly calculating SCET matrix elements. However, our approach is manifestly gauge invariant as it is based on a physical cross section, similar to the direct calculations. The connection of our approach to these previous methods can be understood as follows: Prior to integrating over real radiation, the collinear expansion reproduces precisely the collinear limit of QCD, which in the SCET approach is immediately encoded in the structure of the SCET matrix element, whereas in the approach of ref. [204] it is obtained from the QCD splitting function. In practice, one advantage of our method is that it can be easily integrated with standard methods of generating Feynman diagrams. One can then use standard methods to evaluate the integrals using IBPs [95,96] and the method of differential equations [165][166][167][168][169] in the reverse unitarity framework [90][91][92][93][94] over the real radiation phase space, keeping only the total momentum k fixed. This intermediate result, dη ij /(dQ 2 dw 1 dw 2 dx), is the bare fully differential beam function, from which one can then project out the desired beam functions.

Factorization
N -jettiness is an inclusive event shape that yields an N -jet resolution variable. It was first introduced in ref. [206], and its factorization was derived using SCET in refs. [72,206,207]. Since the same beam function appears for all T N , we focus only on the simplest case T 0 , also known as beam thrust, that is relevant to color-singlet processes. Beam thrust is defined as [206,207] Here, q 1,2 are the Born-projected momenta of the incoming partons, given by where as before Q and Y are the invariant mass and rapidity of the color-singlet final state h, respectively. The sum in eq. (6.10) runs over all final-state particles excluding h, and as usual all final-state momenta are taking as incoming. The Q a,b are measures that determine different definitions of 0-jettiness. The original definitions are [72,208] leptonic: 12) The precise choice does not affect the calculation of the beam function, but it becomes important for the calculation of power corrections [129]. We note in passing that at subleading power, the leptonic definition is clearly preferred as it gives rise to smaller power corrections that the hadronic definition [117,129]. At small T 0 Q, the cross section can be factorized as [72] dσ As indicated, this factorization holds up to power corrections suppressed by T 0 /Q that were studied in refs. [117,118,129,209,210] and the relevant SCET operators have been derived in refs. [120][121][122]. In the case of fiducial cuts applied to the decay products of h, these corrections can be enhanced as O( T 0 /Q) [211]. Furthermore, starting at N 4 LO it also receives contributions from perturbative Glauber-gluon exchanges that are not captured by eq. (6.13) [174,175]. The beam function B i (t, x, µ), sometimes also referred to as the virtuality-dependent beam function, appears in the factorization of all T N [206], deep-inelastic scattering [212], and in the factorization of color-singlet processes in the generalized threshold limit [26]. It is known at NNLO [97,195,196,208], and we compute it at N 3 LO for all partonic channels in our companion paper [106]. Previous progress towards the calculation of the quark beam function at N 3 LO was made in refs. [101][102][103].
In eq. (6.13), the beam functions are defined to measure the Q a,b -independent combinations t a = −q − 1 k + and t b = −q + 2 k − , while the measurement-dependent normalization factors Q a,b only arise in the convolution in eq. (6.13). This definition naturally arises because T 0 simplifies in the n-collinear limit to lim n−coll.
and similarly in then-collinear limit. The soft function in eq. (6.13) only differs between quark annihilation and gluon fusion, but is independent of quark flavors, and we suppress the explicit color index in eq. (6.13). S(T , µ) is a hemisphere soft function for two incoming lightlike Wilson lines. Through NNLO, it is equal to the hemisphere soft function for e + e − → dijets [72,213], which itself is known at NNLO [72,97,[214][215][216][217][218].

Calculation of T N -dependent beam functions
Since the collinear limit of T 0 given in eq. (6.14) only depends on the total momentum k µ of all real emissions, the T N beam function can be calculated using the method outlined in section 6.2. In contrast, the soft limit of eq. (6.10) requires knowledge of all individual momenta {k i }, and thus can not be calculated in this fashion.
Using eqs. (6.6) and (6.14), we can calculate the bare beam function kernel as The zero-bin for the T N beam function is known to be scaleless and thus vanishes in pure dimensional regularization, and hence need not be included explicitly [206]. Note that w 2 > 0 implies t > 0, which we keep implicit. The bare kernel contains UV divergences from the limit w 1 = 1 − z → 0 and w 2 = t/Q 2 → 0, which are both regulated using dimensional regularization. The divergences from small t can be made manifest through the standard identity 1 µ 2 where [ln n x/x] + is the standard plus distributions. Following section 6.2, we obtain the renormalized matching as where the structure of convolution in t [97,208] is made explicit. In practice, it is more useful to perform the renormalization in Fourier or Laplace space, where the convolution in t turns into a simple product. In particular, the structure of Z i B can be easily predicted from the beam function RGE in Fourier space. For details on this, we refer to ref. [106].
We have implemented the described procedure at one loop through O( 4 ) and at two loops through O( 2 ), as required for the calculation of the three-loop beam function. We use the collinear limit of the cross sections for Higgs and Drell-Yan production to extract the gluon and quark beam functions, respectively. As intermediate checks, we verified that the UV and IR counterterms correctly cancel all appearing divergences. The final renormalized results agrees with the NNLO results reported in refs. [195,196], and the higher-order terms in agree with ref. [100]. Our bare results are provided as ancillary files.
6.5 q T beam functions

Factorization
The factorization of the transverse-momentum ( q T ) distribution of a colorless probe h in the limit q T Q was first derived by Collins, Soper, and Sterman (CSS) in refs. [20,219,220] and elaborated on in refs. [33,34,37,164]. The factorization was also discussed using SCET in refs. [73][74][75]221]. The factorized cross section is commonly formulated in Fourier (impact parameter) space, with b T being Fourier-conjugate to q T , as this significantly simplifies the resummation of large logarithms [222]. We write the factorized q T spectrum as It receives power corrections suppressed by q 2 T /Q 2 , which were studied at fixed order in perturbation theory in ref. [130]. The study of their all-order structure has been initiated using the SCET operator formalism in refs. [116,[120][121][122]223], and and their nonperturbative structure has been explored in refs. [224,225]. These corrections are enhanced as O(q T /Q) when applying fiducial cuts to h [211], but for Drell-Yan and Higgs production can be uniquely included in the factorization theorem [226], and are also linear when one includes radiation from massive final states [227].
TMD factorization is complicated by the fact that the bare beam and soft functions not only contain IR and UV divergences, but also so-called rapidity divergences. These must be regularized using a dedicated rapidity regulator, and after removing the regulator this gives rise to the rapidity renormalization scale ν. Several such regulators are known in the literature [73-75, 130, 164, 172, 219, 221, 228, 229], leading to several equivalent schemes for defining TMD beam and soft functions. It is also common to combine beam and soft functions into a ν-independent TMDPDF as where ζ i ∝ ω 2 i is known as the Collins-Soper scale [219,220]. The TMD beam and soft functions appearing in eq. (6.18) are known at NNLO in various regulators [98,99,197,198,200,201,[230][231][232][233]. The quark beam function and the soft function are also known at N 3 LO [104,234] using the exponential regulator of ref. [221].
An important remark is in order concerning differences between quark-and gluoninduced processes. In the quark case, eq. (6.18) exactly applies, while the gluon beam functions can also depend on the gluon helicity due the vectorial nature of b T . As first pointed out in ref. [37], the gluon beam function can be decomposed into a polarizationindependent piece B 1 and a polarization-dependent piece B 2 as Here, b µ ⊥ is a Minkowski four vector with b 2 ⊥ = −b 2 T , and g ρλ ⊥ is the transverse component of the metric tensor. In this case, the hard function in eq. (6.18) also depends on the helicities of the colliding gluons.
We will only focus on the production of scalar particles such as a Higgs boson, where the hard function has the trivial helicity structure Thus, the only combination that enters the factorized cross section in this case is where we suppress the arguments of all functions for brevity. SinceB 2 describes a spin flip of the incoming gluon, it vanishes at tree level, and thus theB 2B2 term first contributes at O(α 2 s ). Thus, for a scalar process, theB 2B2 term does not show up in the strict n-collinear limit, which hence can be used to calculateB 1 in the same fashion as for the quark case. Nevertheless,B 2 could be calculated with the same technique for a different process that induces a cross termB 1B2 , for example the production of a pseudoscalar probe h. We also note that sinceB 2 is already known at NNLO [99,201], theB 2B2 term in eq. (6.22) is already known at N 3 LO.

Calculation of q T -dependent beam functions
In our setup for the differential hadronic cross section, eq. (2.9), we measured the transverse momentum of h indirectly through as by momentum conservation k T = q T . Both are defined as the magnitude of a d − 2dimensional vector, with the associated solid angle already integrated over in the phase space measure. The q T measurement can also be defined in different schemes to account for extending the transverse vector into d − 2 dimensions, but the scheme dependence must cancel in the renormalized beam functions. For a more detailed discussion, see e.g. ref. [233]. Using eqs. (6.6) and (6.23) together with the leading-power relation Q 2 = zs, we obtain the matching kernel of the beam function as Here, the superscript naive indicates that this is not yet the final result for the bare matching kernel, as it requires further manipulation. First, we note that eq. (6.24) contains divergences as x → 1 or z → 1 that are not regulated by dimensional regularization, and are a manifestation of the aforementioned rapidity divergences. In our setup, we must regulate these with a regulator that acts only on the total radiation momentum k µ , but not on individual emissions. The only such regulator known in the literature is the exponential regulator of ref. [221], where one inserts a factor exp[2τ e −γ E k 0 i ] into the phase of each real emission k i . Inserting this regulator into eq. (6.24) and solving the δ functions, we obtain where we defined the so-called label momentum of the beam function as ω = Qe Y . In eq. (6.25), all divergences as x → 1 and z → 1 are manifestly regulated by the exponential, and any leftover divergences are regulated by dimensional regularization. As indicated, the limit τ → 0 should be taken before the limit → 0.
To proceed, we Fourier transform to the conjugate b T space, which trades convolutions in q T for simple products in Fourier space. In d − 2 dimensions, the Fourier transform reads We can then apply the zero-bin subtraction to subtract overlap with the soft function, see section 5.2, which for the exponential regulator is equivalent to dividing by the soft function in Fourier space [98]. This in turn completes the manipulations that forced us the introduce the label naive before their execution. Next, we can apply the usual UV and IR counterterms to obtain the renormalized matching kernel as where following ref. [234] we identify the rapidity renormalization scale as ν = 1/τ . The all-order structure of the beam function counter termZ i B can be predicted from the beam function RGE, which we show in detail in ref. [105].
We have implemented the described procedure at NLO through O( 4 ) and at NNLO through O( 2 ), as required for the calculation of the three-loop beam function. We use the collinear limit of the cross sections for Higgs and Drell-Yan production to extract the gluon and quark beam functions, respectively. Since the bare soft function required in eq. (6.27) has not been published beyond NLO, we have similarly calculated it from the soft limit of the cross section. Our bare results agree with those of refs. [98,99], and the renormalized beam functions also agree with ref. [233]. We provide these beam functions, as well as the bare two-loop soft function, in ancillary files.

Collinear expansion of rapidity distributions
Computing analytic coefficient functions at high orders is a complicated task, and finding suitable approximations can be vital. Here we demonstrate that our expansion techniques have the potential to approximate the rapidity spectrum of color neutral hard probes. We perform a computation of the first two terms in the collinear expansion of the rapidity distribution of the Higgs boson produced via gluon fusion at NNLO. This application also demonstrates that our technique allow one to relatively easily obtain predictions beyond leading power of the kinematic expansion.
The required partonic matrix elements were calculated exactly in ref. [53], and the differential distribution was obtained for example in ref. [93]. Currently, this observable is known at N 3 LO computed via a threshold expansion [5] and via an approximate differential computation [4]. The exact computation of the partonic coefficient function is still elusive due to its extreme difficulty, and a collinear expansion of the same could provide a useful ingredient in future phenomenological studies.
We integrate out the degrees of freedom of the top quark and work in an effective theory that couples the Higgs boson directly to gluons [235][236][237][238][239][240][241][242]. We generate all required Feynman diagrams with QGRAF [243] and perform their collinear expansion up to second term as illustrated in section 4. We then employ IBP identities [95,96] to reduce the expanded diagrams to master integrals, which we then compute using the framework of reverse unitarity [90][91][92][93][94] and the method of differential equations [165][166][167][168][169]. With this we obtain the bare partonic coefficient function dη ij dQ 2 dw 1 dw 2 dx w 2 ∼λ 2 (7.1) expanded up to the second term in ω 2 . Next, we perform a variable transformation from (ω 1 , ω 2 ) → (z 1 , z 2 ) via eq. (2.10) and (2.18), and replace the variable x by ξ via x = ξ(z 1 + z 2 ) 2 ξz 1 (1 − z 2 ) + z 2 (1 + z 1 ) ξz 2 (1 − z 1 ) + z 1 (1 + z 2 ) . (7. 2) The expansion in ω 2 is comparable to an expansion inz 2 = 1 − z 2 , as can be seen by applying the rescaling transformation of eq. (3.4). We find Introducing the variable ξ has the advantage that its integration domain is independent from z 1 and z 2 and ranges from 0 to 1. Next, we expand the partonic coefficient function after this change of variables up the second power inz 2 and integrate over ξ. The result is an approximation for the partonic coefficient function of eq. (2.18) with the observable integrated out. We perform UV renormalisation and combine our partonic matrix elements with collinear counter terms in order to obtain a finite partonic coefficient function through NNLO.
Obtaining the equivalent expansion inz 1 can easily be done by simply relabelling the variables,z 1 ↔z 2 . With this we obtain the following approximation for the full renormalized partonic coefficient function, dη R, approx. The last term in the above equation removes the overlap in the two expansions. The hadronic cross section expanded to this order is then obtained by inserting eq. (7.4) into eq. (2.20), Note that the leading-power limit of eqs. (7.4) and (7.5) precisely correspond to the leadingpower generalized threshold factorization theorem of ref. [26], cf. their eqs. (17) and (18). We have implemented the approximate partonic coefficient function in eq. (2.20) in a private C++ code. Note that we only expanded the NNLO correction to the partonic coefficient coefficient, but keep the lower orders exact. To illustrate our results numerically, we evaluate eq. (7.5) for the LHC with a center-of-mass energy of 13 TeV using the MMHT14 parton distribution functions [244]. Figure 1 shows the rapidity distribution obtained with this collinear expansion normalized to the exact results obtained from ref. [245]. The green line shows our result using only the first term in the collinear expansion, while the red line shows the result including also the second term in the collinear expansion. The blue band in the figure represents the variation of the cross section under a variation of the factorization and renormalisation scale by factor of two around their central values µ F = µ R = m H /2. We observe that the collinear expansion approximates the shape of the rapidity spectrum quite well, in particular towards large values of |Y |. This is kinematically expected, as large rapidities enforce all final-state radiation to be collinear to the corresponding incoming parton, such that the collinear expansion is in fact the correct kinematic limit, see also ref. [26]. In addition, including the second-order term in the expansion clearly improves the results, illustrating that the collinear expansion indeed can be used to produce systematically improvable approximations of key collider physics observables.
The computation of the expanded partonic coefficient functions was greatly simplified compared to the computation of the exact result obtained e.g. in ref. [245]. Explicitly, the complexity of the analytic formulae is greatly reduced, and the function space required to express the coefficient function is much simpler. We expect that a similarly drastic simplification will also occur when applying our method at N 3 LO, which is a natural application of this research.

Conclusions
We have developed a method to efficiently expand differential cross sections for the production of colorless final states in hadron collisions around the particular kinematic limit that all hadronic final-state radiation becomes collinear to one of the colliding hadrons. This yields a generalized power expansion in a power counting parameter λ characterizing this limit.
A key feature of our method is that the expansion is systematically improvable, as it allows to compute to arbitrary order in the power counting parameter λ. Furthermore, λ is treated as a purely symbolic power counting parameter agnostic of the actual observable. This greatly simplifies the expansion, as it can be carried out at the integrand level, i.e. before any phase space or loop integrations are carried out. Subsequently, carrying out phase space and loop integrals is greatly facilitated as integrands become simpler as a result of the expansion. Moreover, the expanded integrands have again a diagrammatic nature very much like the original Feynman integrands they were derived from. This observation makes it manifest that widely used and powerful loop integration techniques like IBP relations and the method of differential equations are applicable to the coefficients of the collinear expansion. We also stress that the basic functions (the so-called master integrals) required in the computation of higher orders in the expansion are already obtained in the lowest few nontrivial orders of the expansion.
Our method also sheds light on the connection between the collinear limit of hadronic cross sections and factorization theorems derived in SCET. The latter include so-called beam functions, universal quantities defined as hadronic matrix elements of collinear fields in SCET, which can be related to standard light-cone PDFs through convolutions with perturbative matching kernels. We have shown that these kernels are precisely given by the first term in a strict collinear expansion of hadronic cross sections. As a first application of this, we reproduced the matching kernels for the N -jettiness and q T beam functions at NNLO from a collinear expansion of the NNLO cross sections for the Drell-Yan process and for Higgs boson production in gluon fusion. The analytic results of this computation are provided as ancillary files together with the arXiv submission of this article.
As another application of the collinear expansion, we have demonstrated its usefulness to efficiently calculate approximate hadron collider cross sections. By combining the collinear expansion with the limit where one partonic momentum fraction becomes equal to its Born value, x i → x B i , we obtained the first two terms in the collinear expansion of the rapidity distribution of a Higgs boson produced in gluon fusion through NNLO in QCD perturbation theory. This example illustrates not only that key collider observables can be approximated with high accuracy using our technique, but also that results beyond the leading power can be easily obtained.
In summary, the method of collinear expansions is a great tool to study the infrared limit of QCD. At leading power in the collinear expansion, it provides access to the universal beam functions governing the collinear limit, which we employ to calculate the T N and q T beam functions at N 3 LO in two companion papers [105,106]. We also believe that the collinear expansions will similarly shed light on the universal structure of hadron collision processes beyond the leading power. Finally, it provides a powerful tool to achieve cuttingedge phenomenological predictions at very high orders in perturbation theory.