An Exponential Regulator for Rapidity Divergences

Finding an efficient and compelling regularization of soft and collinear degrees of freedom at the same invariant mass scale, but separated in rapidity is a persistent problem in high-energy factorization. In the course of a calculation, one encounters divergences unregulated by dimensional regularization, often called rapidity divergences. Once regulated, a general framework exists for their renormalization, the rapidity renormalization group (RRG), leading to fully resummed calculations of transverse momentum (to the jet axis) sensitive quantities. We examine how this regularization can be implemented via a multi-differential factorization of the soft-collinear phase-space, leading to an (in principle) alternative non-perturbative regularization of rapidity divergences. As an example, we examine the fully-differential factorization of a color singlet's momentum spectrum in a hadron-hadron collision at threshold. We show how this factorization acts as a mother theory to both traditional threshold and transverse momentum resummation, recovering the classical results for both resummations. Examining the refactorization of the transverse momentum beam functions in the threshold region, we show that one can directly calculate the rapidity renormalized function, while shedding light on the structure of joint resummation. Finally, we show how using modern bootstrap techniques, the transverse momentum spectrum is determined by an expansion about the threshold factorization, leading to a viable higher loop scheme for calculating the relevant anomalous dimensions for the transverse momentum spectrum.


I. INTRODUCTION
Many phenomonologically important observables of Quantum Chromodynamics (QCD) are transverse momentum sensitive. That is, they are defined as a measurement that directly puts a constraint on the momentum flowing perpendicular to some fiducial jet axis, without a corresponging cut on the rapidity. Examples include the transverse momentum ( Q T ) distribution of generic high-mass color-neutral systems (Drell-Yan, Higgs, vector boson pair,. . . ) in hadron-hadron collisions, semi-inclusive fragmentation of hadrons, the scalar sum of transverse momentum magnitudes as found in jet or beam broadening [1][2][3][4][5][6][7], and vetoes on the transverse momenta of clustered jets [8][9][10] 1 . An universal feature of all such transverse momentum sensitive factorizations is the presence of rapidity divergences. In a naive soft or collinear sector, one encounters integrals over the light-cone components of the participating partons, which dimensional regularization fails to regulate. This is due to the fact that dimensional regularization breaks any possible dilatation invariance of a theory, but not the Poincare invariance. Thus classes of momenta with differing invariant masses can be distingushed via their relative scaling with respect to the dimensional regularization mass scale µ. However, classes of momenta differing only by their boost are not distinguished by the boost invariant dimensional regularization. These rapidity divergences are then a necessity when soft and collinear modes exist at the same invariant mass scale, as found in a transverse momentum sensitive observable. From a more practical point of view, these rapidity divergences (an artefact of using a factorized form of the cross-section) are directly tied to a large logarithm of the fixed order QCD calculation (QCD itself is rapidity divergence free). Much literature has been devoted to the issue of a convenient scheme to regulate the light-cone integrals. Once regulated, the divergences are isolated, cancelled at the level of the physical cross-section, and the residual logarithms left from the isolation are exponentiated either by hand or by evolution-equations, thereby controlling the large logarithm found in fixed-order perturbation theory.
The transverse momentum spectrum for color-singlets in particular has long been a critical quantity to understand the factorization properties of QCD. One simply wishes to know the differential cross-section for the relative momentum of the color singlet object with respect to the beam axis, while being inclusive over all other radiation in the event. Taking the transverse momentum to be small relative to the hard scale (Q) involved in the production of the observed particles (the invariant mass scale of the Drell-Yan pair or the Higgs boson, for example) implicitly constrains how the recoiling QCD radiation moves with respect to the beam axis. In this limit, the cross-section is dominated by either soft radiation, or emissions collinear to the beam. Such soft and collinear emissions constitute the infra-red structure of perturbative QCD, which presents an underlying universality due to the fact that the emissions in these different kinematic regimes (hard, collinear, or soft) factorize from each other, and do not quantum mechanically interfere 2 . Due to this infra-red sensitivity, the fixed order expansion for the cross-section becomes dominated by large logarithms of the hard production scale Q to the infra-red recoil scale Q T . The factorization allows one to resum these large logarithmic contributions to all orders. When this resummation is combined with the fixed order distribution that is not singularly enhanced, one achieves a remarkably precise description of the QCD spectrum, giving a benchmark for theory versus 3 experimental predictions. If we denote the potentially large logarithm as L Q T = ln Q T Q , and assign the scaling α s ∼ L −1 Q T when the logarithm is large, then the current state of the art for these resumations is the N 2 LL+NLO accuracy, where one in the fixed order result includes contributions for up to two final state recoiling partons, and has resummed all logarithms up to contributions that scale as α 2 s L Q T [21][22][23][24][25][26][27] 3 . In this paper, we introduce a new methodology for calculating the control quantity for rapiditiy divergences, the rapidity anomalous dimension, see Refs. [1,2] 4 . We exploit the fact that cross-sections often have multiple singular regions with distinct scalings of the lowscale modes, leading to distinct factorization formula, see Ref. [35]. These factorizations, however, while resumming different logarithms, must be consistent with each other at any fixed order in perturbation theory, since they describe the same cross-section. This allows us to calculate in the threshold region the transverse momentum spectrum, and through consistency with the more standard transverse momentum dependent parton distribution function (TMD-PDF's) factorization, extract the rapidity anomalous dimension. We can then combine the technology developed for threshold calculations (see Refs. [36][37][38][39][40][41][42]) with modern bootstrap techniques from amplitudes to push the calculation to three-loop order. In a companion paper, two of us will present the full three-loop rapidity anomalous dimension phenomologically relevant for collider experiments. Though we deal mainly with the transverse momentum spectrum in hadron-hadron collisions, however, we believe that our approach is widely adaptable to many transverse momentum sensitive observables, at least where one can understand the analytic structure of the fixed order calculation. We refer to this methodology as the "exponential regulator," since it implements an exponential cutoff in the total energy of the final state. Alternatively, one may think of it as a threshold regulator, where one imposes in addition to the transverse momentum observable a constraint on the total energy of the soft radiation crossing the cut. It effectively acts as a gauge invariant cut-off on the rapidity integrals.
The factorization approach to resummation we will adopt is that of Soft Collinear Effective Field Theory, [43][44][45], which gives a precise set of rules for determining the all-orders form of the factorized formulae. In general one seeks to write a cross-section sensitive to a singular region of phase-space as a product of functions of the form: The H denotes calculation of the hard process, the B a beam function for the initial state radiation off of the colliding hadrons 5 , J i functions for any possible final state jets, and the S the contribution from soft wide angle radiation. Each function has a field theoretic operator definition, and the ⊗ denotes a convolution over the contribution from each sector to the relevant observable and any possible momentum recoil. The factorized functions summarize the contribution from "on-shell" modes of QCD with a specified scaling, and can be calculated independently. This convolution structure is to be expanded according to the scaling of the modes to produce a formula that is homogenous in the power counting, a procedure known as the multipole expansion [46,47]. Since the multipole expansion enforces a homogeneous power counting in each convolution, one is prevented from developing a large 4 logarithm in the effective theory matrix element. Instead, one is often rewarded with potentially multiple-divergences of the naive function. That is, one is trading the large logarithm of the perturbative expansion for an explict divergence in the calculation of effective theory matrix element. Several variations of the SCET formalism has been applied to transverse momentum distributions before, see Refs. [2,[48][49][50]. In the operator based factorization literature, three approaches have appeared to accomplish this task, the collinear anomaly [48], the Collins-Soper equation [19,33,34,51], and finally the framework of the rapidity renormalization group (RRG) [1,2]. The outline of the paper is as follows: first we review the topic of transverse momentum resummation in the SCET formalism. For a general regularization scheme, we show that as long as the regulator is implemented symmetrically with appropriate subtractions in the different sectors, the rapidity resummation's scheme dependence is fixed by the hard function's scheme dependence. The subtractions themselves are regulator dependent. Since the hard function is free of rapidity divergences, this necessarily implies a universality to the rapidity anomalous dimension regardless of regulator+subtractions. Having established the factorization framework, we introduce the exponential regulator with the example of the oneloop calculation, which is defined by taking the limit of the fully differential soft functions. We then discuss the relation of transverse momentum and threshold factorizations from their connection with the fully differential functions, showing that the exponential regulator necessarily calculates a rapidity renormalized transverse momentum soft function. Lastly we show that the utility of the regulator at higher loops lies in reducing the calculational problem of the integrals to that of the threshold soft function. Using the extensive work done on this subject, we make an ansatz of the fully differential soft function in terms of harmonic polylogarithms (HPL's), and demonstrate how to reproduce existing results at one and two loops by bootstrapping. We also present partial results at three loops using this technique, and the full result is deferred to a companion paper. Finally, we conclude with thoughts on future directions. Some technical details are collected in appendicies.

II. REVIEW OF TRANSVERSE MOMENTUM FACTORIZATION
The factorization theorem for the Drell-Yan transverse momentum distribution takes the form 6 : dσ Q 2 is the invariant mass of the Drell-Yan pair, and σ 0 is the leading-order (LO) cross section. The light-cone coordinates are defined with respect to the beam axis, and satisfy: n 2 =n 2 = 0, n ·n = 2 , 6 Several equivalent definitions can be found in the literature. They differ ultimately with regards to how the rapidity divergences are regulated, and the necessary subtractions that must take place given the form of the regulation.
When switching to the x A,B momentum fractions of the partons in the hard collision, we intrepret the light-cone momentum components of the Drell-Yan pair as: is the hadronic center of mass of the two colliding nucleons N A and N B , the momenta of which can be written as P A = P + A n/2 and P B = P − Bn /2 respectively. The H qq (which we will shorten to H) is the hard function that contains all the virtual correction to the LO contribution, and is obtained by matching from QCD to SCET. The B n,q/N (which we will shorten to just B n ) are beam functions encoding the energetic emissions along the beam axis, and S is a soft factor encoding the contribution from soft states. In the fully differential form, these functions have the following operator definitions: where b = (b + , b − , b ⊥ ) and d a = N c = 3 for the Drell-Yan process. The χ n field is a gauge invariant quark field operator dressed with a collinear wilson, and together with the soft and collinear wilson lines have the respective definitions: In the factorization of the transverse momentum distribution, these functions compute the contribution to the observable from modes with momentum scaling: As a result, some light-cone coordinates need to be set to zero for proper power counting of the multipole expansion, and the relevant beam 7 and soft functions become, Often the momentum modes are called "on-shell," since their dispersion relation satisfies p 2 n = p 2 n = p 2 s = Q 2 λ 2 homogenously, and as λ → 0, they scale to exactly on-shell emissions. These modes have the important property that they are all at the same invariant mass-scale, as depicted in Fig. 1, so that the appropriate effective field theory is SCET II , and are distinguished only with size of their relative energy scale or typical rapidity. Since dimensional regularization is invariant under boosts, one cannot distinguish these modes from each other in an integral with dimensional regularization alone. In so called SCET I theories, modes are distinguished by their invariant masses, and since dimensional regularization breaks dilatation invariance, it suffices to regulate the theory completely and seperate the modes. A further regulator is needed when integrating over the light-cone variables in a typical diagram, and several have been proposed in the literature, each their various strengths and weaknesses. They may be classed into analytic style regulators [2,52], deformations of the wilson line directions [33,53], or finally a mass added to the eikonal propagator (the "δ" regulator) [31,54]. Beyond the obvious requirement of ease of calculational use, one would also demand the regulator preserve gauge invariance, non-Abelian exponentiation [55,56], and a democratic treatment of sectors (at least up to terms that vanish as the regulator is taken to its singular limit). For all regulators that have an explicit mass scale associated, like deformations of the wilson line direction or the δ-regulator, the zero-bin subtraction will not be zero [54,57].
It is important to note the origin of the light-cone singularities. In the factorization theorem of Eq. (2), the TMD-beam functions are localized at either b + = 0 or b − = 0, while the soft function is localized at both. This prevents momentum sharing in these small momentum components, since if we were to perform the fourier transform in Eq. (2), no recoil convolutions would appear in either the n orn directions. The rapidity and mass of the Drell-Yan pair sets these momentum scales once and for all, to leading power. This is a direct consequence of the multipole expansion and the scaling of Eq. (11), and such an expansion is necessary to guarantee no large logarithms appear in the EFT.
Introducing a regulation scheme with the appropriate subtractions, one will have a generalized renormalized definition of the TMD-beam function and soft function (which can be combined together to form a TMDPDF). One removes systematically the light-cone and the ultra-violet divergences from each function. Removing these divergences will introduce a scale at which the divergences are subtracted, which we will generically call ν for the rapidity divergences, and µ for the ultra-violet. Since the physical cross-section is finite, these divergences will cancel between the various functions in the factorization formula, and the variation under the scale where the divergences are subtracted are controlled by generalized renormalization group equations: and for the rapidity renormalization: where b 0 = 2e −γ E . Since divergences cancel in the physical cross-section, the ultra-violet anomalous dimensions satisfy the constraint, Similar constraint is manifestly written for the rapidity renormalization group. We have used the fact that the anomalous dimensions of the two TMD-PDFs should be the same up to relabeling n andn, given that the regularization procedure treats the two beam sectors identically. The arguments of these anomalous dimensions are dictated by the factorization structure in (2). The hard production scale Q 2 ∼ x A x B P + A P − B gets factorized into the large momentum components of the beams sectors, the x A P + A or x B P − B of the n orn collinear sectors respectively. This hard production scale appears in the hard function H, including its anomalous dimension, and so to cancel it, it must reappear in the low scale EFT matrix elements. However, no propagator in the low scale matrix elements has virtuality at this hard scale by construction, so that in the beam sectors, the scale Q can only appear associated with the large light-cone momentum component. Yet it is precisely integrals over these light-cone components that give rise to the rapidity divergence, so that the beam functions must depend "anomalously" on the ratio ν The soft function's ν-dependence is 8 then constrained by the fact the anomalous dimensions sum to zero. Importantly, Lorentz invariance dictates that in the anomalous dimension, the logarithm of x A P + A must combine with the logarithm of x B P − B to form the scale Q 2 , so that at most one logarithm of the ν scale can appear in the logarithm of the renormalized functions, see Refs. [58,59]. Thus the rapidity scale ν does not appear in the anomalous dimension of Eq. (15), and the ultra-violet anomalous dimensions have the form: By consistency, the µ dependence of the rapidity anomalous dimension 8 is controlled by the cusp anomalous dimension for wilson lines: This leads to an all-orders form for the rapidity anomalous dimension: Lastly, we comment on the scheme dependence of how the rapidity divergences are isolated and removed from the bare functions. Schematically, the cross-section in coordinate space has the form: Then we may consider the derivative: where we have used the fact that the total derivate with respect to Q 2 acts as partial derivative on the momentum components on the Drell-Yan pair: and x A P + A and x B P − B only appear in B n,⊥ and Bn ,⊥ with ν in combinations of ν according to the renormalization group equations. The left-hand side of Eq. (23) is independent of how the low-scale matrix elements are regulated. Indeed, the hard-function 9 is the same in a wide variety of infra-red observables both with and without rapidity divergences, and thus can be calculated with or without intermediate rapidity regularization. Then the cancellation of rapidity divergences between the zero-bin subtracted collinear and soft functions allows us to conclude: Thus all scheme dependence of the rapidity anomalous dimension is directly controlled by the scheme dependence of the hard function, and the anomalous dimension is independent of the regularization procedure 9 . There is a further scheme dependence involved with the decomposition into the cusp and non-cusp contributions to the rapidity anomalous dimension, however, this depedence is completely controlled as an initial condition to the solution of the differential equation (20). For all functions and parameters appeared here, we default to a fixed order expansion around α S /4/π as: Note that for soft functions, we have assumed non-Abelian exponentiation when performing the expansion [55,56]. We gather in App. A each of anomalous dimensions to the highest known perturbative order.

III. THE EXPONENTIAL REGULARIZATION PROCEDURE
We now explain how one can calculate the rapidity anomalous dimension of the soft function of Eq. (13) through the exponential regulated soft function. We first note that the origin of the divergences lay in the multipole expansion between the beam and soft sector's light-cone components in Eq. (3). We are free to consider then not a strict expansion, but the limiting behavior of the functions as the light-cone components are localized. Specifically, we consider the soft function in coordinate space: Since no information is lost by taking b + = b − = ib 0 τ /2 -as the fully-differential soft function is always a function of the product b + b − by the RPI III transformations of the effective theory (see Ref. [60]) -we use the same notation, i.e. S with no subscript for the soft function here. A picture for the coordinate space soft function is depicted in Fig. 2 with τ = 1/ν. We will show later that the ν = 1/τ is indeed the artificial scale appearing in the rapidity regularization once we take the limit of τ → 0. By imposing the energy constraint on the momentum crossing the cut in the diagrammatic expansion, we regulate the integrals over the light-cone components of momenta. This can be seen more clearing in the momentum space, where the measurement of τ forms an exponential damping factor for the rapidity divergence. It is in this sense we call the function the exponential regulated soft function. This deformation of the transverse momentum dependent (TMD) soft function is particularly revealing, since it is directly relatable to the threshold soft function of Ref. [61], which we define as: The limit b ⊥ = 0 can be taken smoothly, both before and after renormalization, since no onshell singularities are probed in this limit. One can perform the deformation to any SCET II soft functions, forming a general regularization procedure for these theories. The relation to the standard threshold soft function, and the fact that the limit is smooth implies several important features about the exponential regulated soft function. First, its UV anomalous dimension is the same as the threshold soft function: To qualify as a valid regularization scheme for the TMD soft function, it also has to satisfy the following condition in the τ → 0 limit (for a derivation of this result, we refer to Sec. IV): from which we can derive constraints on the function form of S( b ⊥ , τ ; µ). To make our statement explicit, let's first use Eq. (35) to write the exponential regulated soft function as 10 : where the second term on the RHS is the µ-independent part of the soft function, and has a well-behaved series expansion about b ⊥ = 0. By demanding Eq. (36) holds, we obtain the following equation: This is a non-trivial constraint, since at each order in pertubation theory, the double logarithmic contribution to the τ → 0 behavior of the µ-independent part must be fixed by the cusp anomalous dimension, and higher order logarithms are determined from the betafunction, both of which form important checks on any calculation of the function. The same regularization can be easily adapted to regulate the rapidity divergence in the TMD-PDFs through the use of fully differential beam function, which we will defer to future work. This regularization of the TMD soft function has several features to commend it. Firstly, since it is defined via a measurement constraint on the final state radiation, it is manifestly gauge invariant. Non-abelian exponentiation also follows trivially, which we have used in writing down Eq. (37), since the measurement factorizes in its Laplace form to act on each final state parton. Lastly, as seen from Eq. (34), we can actually realize the exponential regulated soft function from its Taylor series expansion about the threshold limit, where all integrals will be reducible to known master integrals. As explained in Sec. V, this means by matching the Taylor series to an ansatz of special functions, we can deduce the full transverse-space dependence of the function from a finite number of terms. Being able to deduce the full transverse-space dependence is critical to being able to construct the rapidity anomalous dimension. In the all-orders form of the exponential regulated soft function in Eq. (37), transverse-space dependence is entirely controlled by its µ-independent part, which depends on its arguments solely through the scaleless ratio of x = − b 2 ⊥ /b 2 0 /τ 2 (neglecting the scale dependence in α S ). It is the Taylor series about x = 0 that is probed by the threshold limit, but it is the x → ∞ that controls the rapidity anomalous dimension in Eq. (36).

12
Technically, an infinite number of terms would be necessary, assuming an infinite radius of convergence. However, the space of functions appearing in perturbative calculations is tightly constrained, allowing the full dependence to be deduced from only a finite number of terms even when the taylor series has a finite radius of convergence. It is fascinating that there is a mother function relating both threshold resummation to the transverse momentum resummation: both can be obtained by taking appropriate limits of a single function.
To illustrate how the regulator actually works, we take the one-loop calculation of the soft function as an example. The relevant diagrams are depicted in Fig. 3. For light-like Wilson lines, Fig. 3(b) vanishes and we only need to consider Fig. 3(a) and its conjugate. The bare exponential regulated soft function is given by the integral where d = 4 − 2 . We work in the MS scheme by a redefinition of the bare scale µ 2 0 = µ 2 e γ E /(4π). Note that b ⊥ is in two dimension, while k ⊥ is in 2 − 2 dimension. Due to rotation invariance in the ⊥ plane for Drell-Yan production, we let b ⊥ = | b ⊥ |(0, 1). Without loss of generality, we can parameterize k ⊥ as It is also convenient to use light-cone coordinate for the integral measure, where is the area of unit sphere in n dimension. Integrating out θ, making the following change of variables 13 with the Jacobian 1/(2v), and using the on-shell delta function, we arrive at where J n (z) is the Bessel function of the first kind. The variable v is related to the rapidity of soft gluon by v = exp(−2Y ). It is clear from Eq. (44) that without the threshold regulator factor, the v integral diverges at the end points of infinite rapidity. This is the so-called light-cone/rapidity singularity. The exponential regulator provides an exponential damping factor at infinite rapidity. The resulting v and k ⊥ integrals can be done in closed form, giving It is straightforward to expand the above expression using, e.g., HypExp [62] to arrive at The renormalized fully differential soft function at one-loop is then obtained by removing the poles, The exponential regulated Q T soft function is obtained by taking the τ → 0 limit and keeping only the non-vanishing terms, Once we identify ν = τ −1 , we can make smooth connection with the rapidity RG formalism, and check that Eq. (48) satisfy the µ and ν RG equation.

IV. TRANSVERSE MOMENTUM AND THRESHOLD FACTORIZATION
In contrast to the factorization of Eq. (2), a distinct formula was proposed in Ref. [63], which did not perform the multipole expansion: This factorization utilized fully differential beam and soft functions, that are sensitive to the total momentum flow crossing the cut in the diagrammatic expansion of the functions. Since the multipole expansion was not performed, large logarithms may still remain in its perturbative expansion, even after renormalization group evolution. One can find consistent factorization theorems that utilize these fully differential functions in a multi-differential measurement of beam thrust and Drell-Yan transverse momentum, see Refs. [35,64,65] 11 . The relative values of b + , b − are unimportant, since they always appear in the product of b + b − as explained earlier. Thus to examine the limit to the TMD soft function, we set b + = b − = t and write: Now the exponential regulated soft function is connected to the fully differential by the analytic continuation from τ → −2it/b 0 . One can equally well use this Fourier transformed function as a definition of the exponential regulated soft function, instead of the Laplace. However, given that much work on soft threshold integrals has been done in Laplace space, as well as to avoid a proliferation of imaginary numbers, we found it convenient to adopt the Laplace space definition, i.e. using τ as the new argument.
To understand the origin of our central result, Eq. (36), we simply approach the factorization for the differential spectrum of the Drell-Yan pair from two different limits. For a Drell-Yan pair, the allowed phase space at zero rapidity in terms of the transverse momentum and the residual partonic energy scale is plotted in Fig. 4. First we consider the factorization starting with the standard inclusive differential Drell-Yan cross-section, at large to moderate Q T , moving along the upper line in Fig. 4. To avoid convolutions, it is simplest to work in position space, and the standard inclusive Drell-Yan cross-section admits a factorization into collinear PDFs as: whereσ is the inclusive hard coefficient, i.e. partonic cross section, and the standard collinear PDF is related to the fully differential beam functions by taking the transverse and small lightcone component to zero. This cross-section admits a further factorization in the threshold region, where the hard inclusive coefficient splits into the form factor derived hard function and a soft factor that is fully differential as in Ref. [63], and the PDFs are taken to their threshold expressions, that is, taking the Bjorken scale x → 1 in the momentum space:σ Substituting these functions into Eq. (51), we achieve the threshold factorization for the differential spectrum of Drell-Yan: Thus the fully differential or exponential regulated soft function does appear in a factorization theorem with homogeneous power counting, when the modes are organized as in Fig. 5. Alternatively, we may approach the threshold regime already assuming small transverse momentum. Let's rewrite Eq. (2) as, The TMD-beam functions can then be further factorized into an additional soft factor and the threshold PDF, reminicent of Ref. [69]: Both the PDF's and coefficients for TMD matching to PDF's get expanded.
both sides of these equations have the same rapidity divergences, which on the right hand side are carried by the soft factor alone. This is the same soft factor appearing in the SCET + factorization of the multi-differential beam thrust and transverse momentum phase space, see Ref. [65]. By substituting Eq. (57) in to Eq. (56), we again achieve another threshold factorization for the Drell-Yan process, where now all functions have been refactorized in the threshold power counting as in Fig. 6. Demanding consistency between these two factorizations in their overlapping domain of validity, we conclude: This equality holds at the level of renormalized functions. The left-hand side is free from rapidity divergences, but in the limit b + , b − → 0 (the small τ limit) has a large logarithm at each order in perturbation theory (the limit to zero light-cone position is not smooth). This corresponds to the fact that each factor on the right is naively rapidity divergent. With appropriate regularization and subtractions, these divergences will cancel, making way for the RRG. Following the arguments of the logarithms of the intermediate rapidity renormalization, we are then lead to Eq. (36) similarly to how we concluded Eq. (25). That is, since the rapidity divergences cancel between the three soft functions, we can interprete the fully differential soft function as a direct calculation of the rapidity renormalized soft function. Note that the expansion is very important. When we factorize the threshold region in the inclusive hard coefficient, we perform no expansion between the energy of the hadronic final state and its transverse momentum. In contrast, when further factorizing the small transverse momentum factorization in the threshold limit, an expansion between the energy of the final state and its transverse momentum has already been performed to arrive at Eq. (2). The expansion in (58) is the common region of validity between these two approaches to the threshold region 12 . The smoothness of the limit | b ⊥ | → 0 is also seen from the threshold factorization of Eq. (51) using Eq. (52). If we fourier transform Eq. (51) with respect to Q T , and take the limit | b ⊥ | → 0, we recover the traditional factorization of the threshold Drell-Yan spectrum, see for instance Ref. [70]. Since this factorization has no singularities associated with its localization at zero impact parameter, we conclude the limit | b ⊥ | → 0 is smooth to all orders, which is born out by explicit calculations up to and including three loops. Again, this is not surprising since the resummation structure driven by the renormalization group for the threshold factorization is resumming large logarithms associated with the light-cone variables n · b andn · b, not the transverse momentum.
Similar functions appear in joint resummation (see Refs. [71][72][73]) that seeks to combine threshold and transverse momentum resummation. In particular, a similar refactorization to Eqs. (57) and (58) was considered. There the authors sought to combine into a single formula the resummation for both the threshold logarithms and the transverse momentum spectrum. Our aim has been distinct, which was to provide a new method for calculating all quantities needed for resummation from a single fully differential function. However, the family of factorizations we have derived would also allow us to examine the structure for genuine joint resummation. We find that there are three distinct factorization theorems, each of which is seperately consistent under ultraviolet and rapidity renormalization, Eqs.
(2), (55) and Eq. (55) with the substitution of Eq. (58). One can consider a merging scheme as derived by [65] that would also attempt to combine both threshold and transverse momentum resummation, such that the scheme is accurate to N 3 LL in all limits. One could also include small-x resummation following [74,75].

V. BOOTSTRAPPING THE FULLY DIFFERENTIAL SOFT FUNCTION
At first sight, the one-loop calculation using exponential regulator in section III doesn't seem to simplify the calculation. Even worse, exponential regulator introduce an extra nontrivial scale τ = 1/ν into the problem, which leads to the appearance of non-trivial analytic function in the one-loop calculation. However, such seeming weaknesses will be shown to be strengths, once we examine the two-loop calculation for fully differential soft function already performed in Ref. [76], where the results are given in terms of polylogarithms up to weight four with rational coefficients. In this section, we shall show that the simple structure of the results in Ref. [76] allows us to calculate the fully differential soft function without actually calculating the corresponding Feynman integrals.
As defined in Sec. II, we can expand the renormalized fully differential soft function in the following exponential, thanks to the on-Abelian exponentiation theorem: The results in Ref. [76] then can be rewritten in terms of Harmonic Polylogarithms (HPLs) 12 Though the appropriate threshold factorization is of the differential spectrum (51).
of Remiddi and Vermaseren [77], taking into account the exponentiation in Eq. (59): where we have only kept the scale independent part by setting µ = τ −1 . c s i are scale independent constant in threshold resummation, whose explicit formula are collected in the appendix. We use C a to denote the Casimir of the initial parton. C a = C F for Drell-Yan production, and C a = C A for Higgs production. H w ≡ H w (x) are HPLs with weight vector . We have used the shorthand notation for the weight vector of HPLs [77] 13 . The exceedingly simplicity of Eq. (60) makes one wonder whether there is simpler way to obtain them, instead of the brute-force calculation done in Ref. [76]. Indeed, we found that the results in Eq. (60) can be obtained using bootstrap method, which we shall explain below.
The bootstrap program is extremely successful in calculating scattering amplitudes in planar N = 4 SYM, in particular for six-point amplitudes. Briefly speaking, for a L-loop planar amplitudes with the Bern-Dixon-Smirnov [78] factored out, one can make an ansatz consists of rational linear combination of transcendental function of transcendental weight 2L. In general the ansatz contains a large number of unknown coefficients. Remarkably, in the case of six-point planar amplitudes, they can be uniquely fixed by expanding the ansatz in the boundaries of phase space, where prediction exist thanks to knowledge of resummation and integrability. This approach is so powerful that even planar five-loop NMHV amplitudes in N = 4 SYM can be obtained in tis way [79].
On the other hand, examples of application of bootstrap method in QCD calculation are less common. The reason is that, in QCD the ansatz is usually much more complicated than in N = 4 at given loop order. For example, the transcendental functions in the ansatz can be multiplied by non-trivial ratio function of kinematics variables, and the transcendental weight can ranged from 1 to 2L in an L-loop amplitude. Furthermore, integrability is lost in QCD, therefore the number of boundary data for fixing the ansatz are much smaller than in N = 4 SYM.
Nevertheless, the simplicity of Eq. (60) is hard to ignored: • At one and two loops, the results are given solely in terms of HPLs with rational coefficients. Furthermore, the indices of weight vector are drawn from the set {0, 1} only.
• The last entry of the weight vectors, or the first entry of the symbols [80], are always 1. This is ensured by that the threshold limit b 2 ⊥ → 0 is a smooth limit. Fully differential soft function admits a simple Taylor series expansion over b 2 ⊥ in that limit. We will explain more about this in the following.
• The first entry of the weight vectors are always 0, at least through to two loops. 13 In this notation, weight vector with n trailing zeros to the left of a 1 is written as n + 1. For example,

19
These observations lead us to make the following ansatz for the scale independent part of the fully differential soft function at L loop: where c s L is the L-loop scale independent constant of threshold soft function. r i ∈ Q are rational numbers. F i (x) are transcendental function with transcendental weight 2 ≤ [F i (x)] ≤ 2L. These can include single HPL H 0, w n−2 ,1 , where w n−2 is a weight vector of length n − 2 with entries drawn from {0, 1}. We also allow F i (x) to be product (multiple) zeta value of weight 2 ≤ m ≤ 2L − 2 and a HPL H 0, w n−m−2 ,1 . The summation is over all possible F i (x). With the ansatz at hand, what remains is to fixed the rational coefficients r i using all possible constraints. We identify two such constraints.
The first constraint comes the fact that rapidity divergence is only a single logarithmic divergence at each order on the exponential, Eq. (59), but the scale independent term in Eq. (61) can contains higher order logarithmic divergence in τ → 0. Therefore, there must be non-trivial cancellation for the higher order logarithmic divergence between the scale dependent part and the scale independent part. As an concrete example, consider the oneloop ansatz: where the scale dependent part is uniquely fixed by RG equation, and we have used c s 1 = 2C a ζ 2 . The linearity of ln τ divergence demands that the ln 2 τ divergence in Eq. (62) should cancel. Using we find that which is in agreement with the result in Eq. (47). In general, at L-loop, the logarithmic divergent terms ln k τ with 1 < k ≤ 2L must cancel between the scale independent part and the scale dependent part. Beyond one-loop, the linearity in the logarithmic rapidity divergence is not enough to completely fix the unknown coefficients. The remaining degrees of freedom have to be fixed using the second constraint, which comes from the expansion in the threshold limit b 2 ⊥ → 0. Following Eq. (33), we notice that an expansion in b 2 ⊥ is possible at the integrand level by expanding the exponential, whereŜ(k, µ) is the fully differential soft function in the momentum space, with k being the total momentum of radiation crossing the cut. We recognize that the first term in Eq. (65) is simply the threshold soft function, Starting from n = 1, we encounter terms of the form Note that S(k, µ) is a function of k + k − and k 2 . By Lorentz invariance, k µ becomes n µ or n µ after the d d k integral. Using b ⊥ · n = b ⊥ ·n = 0, we obtain that only even n survives in Eq. (65). The first non-vanishing term start at n = 2, where= means that the equality only holds after integrating over k. For arbitrary positive integer m, we have where the function f (2m) is given by with d ⊥ = d − 2, and we also define f (0) = 1. We therefore obtain that the fully differential soft function is given by the expansion that is, by insertion of numerator (k + k − −k 2 ) m into the integrand for threshold soft function. By expressing k as the sum of soft momentum crossing the cut k µ = i∈cut k µ i , the numerator insertion can be reduced to scalar master integrals relevant for the threshold soft function via Integration-By-Parts identities (IBP) [81]. At high power of m, the reduction can be very computational-demanding. Fortunately, some trick can be used to ease the effort, which is explained in the App. B. We note that these integrals have been computed to three loops [36-39, 41, 42], in an effort to obtain the soft-virtual corrections for Higgs production at N 3 LO [40,82]. Eq. (71) is one of our main result, which in principle provides any number of boundary conditions for fixing the unknown coefficients r i 14 . As a concrete example, we present below the expansion in x = − b 2 ⊥ /b 2 0 τ 2 through to x 5 at one and two loops for the scale independent part: where we have given the results for SU(N c ) gauge theory with n f number of Dirac Fermion, n s number of complex scalar, both in fundamental representation. δ R is a parameter choosing the specific regularization scheme. δ R = 0 for Four-Dimensional-Helicity (FDH) scheme [83], and δ R = 1 for 't Hooft-Veltman scheme. β gen 0 is the corresponding beta function: Using a generic two-loop ansatz from Eq. (61), it is straightforward to determine all the coefficients by combining the linearity constraint for the ln τ divergence, and the data from threshold expansion. For example, the one and two-loop result in FDH scheme for N = 4 SYM can be bootstrapped from these constraints by setting n f → 2, n s → 3, C A , T f = T s → N c , and δ R = 0: where the constants c s,N =4 i are given in the App. A. By setting n s = 0 and δ R = 1, we also reproduce the QCD result in Eq. (60). We note that the result in N = 4 SYM agrees the leading transcendental part of the QCD result 15 , while there is no lower transcendental weight part in N = 4 SYM. Such relation between N = 4 SYM and QCD was first observed for anomalous dimension of twist two operator [84]. In the context of soft function in SCET, this is true for threshold soft function through to three loops [82].
As outlined above, the calculation of fully differential soft function doesn't require calculating a single Feynman diagram. At least through to two loops, it only relies on the constraints imposed by linearity of ln τ divergence, and more importantly boundary data from threshold soft function. At higher orders, it is limited by the availability of threshold data, the amount computation resource required for IBP reduction for Eq. (71), and potentially also the completeness of the ansatz in Eq. (61) 16 . To illustrate the power of this approach, we consider the calculation of fully differential soft function in QCD at three loops for two relatively simple color structure, C a n 2 f and C a C F n f . These contributions are simple because they correspond to self-energy insertion of one-loop diagram. A representative diagram which contribute to C a C F n f color structure is depicted in Fig. 7. We find that the ansatz in Eq. (61) is sufficient for bootstrapping these color structure at three loops. The results reads The full QCD results require significant more work and we defer them to a second publication.n nn n FIG. 7: Representative diagram for color structure C a C F n f at three loop in QCD.
As explain in section III, the fully differential soft function acts as a mother function to both the TMD soft function and the threshold soft function. It is an easy exercise to obtain TMD soft function by taking the limit τ → 0, and keeping the non-vanishing terms only. Then one can identify ν = τ −1 . We find that as expected, the TMD soft function satisfy the µ and ν evolution equation. The corresponding rapidity anomalous dimension is extracted to be for QCD, and for N = 4 SYM. The QCD result in Eq. (76) agrees with those obtained using different methods [30][31][32]. The corresponding result for the scale independent part are given by + C a n f 10ζ 2 3 + 28ζ 3 9 − 328 81 (78) for QCD, and for N = 4 SYM.

VI. CONCLUSIONS
We have introduced a new method to calculate naively rapidity divergent soft functions, by deforming the soft function's measurement into one that is calculable in dimensional regularization, and can be reconnected to the naive rapidity divergent soft function. The most practical benefit of this method will be the three-loop non-cusp anomalous dimension needed for transverse momentum resummation for Drell-Yan and Higgs production. This is the largest source of certainty on the N 3 LL analysis of these spectra, given that the exact value of four loop cusp anomalous dimension has been found to have negligible impact. This anomalous dimension for QCD processes and the N 3 LL + N N LO transverse momentum spectrum for Higgs production will be presented in companion papers.
Ultimately underpinning this procedure is the factorization of multi-differential crosssection, where the same cross-section in different singular regions experiences different factorizations. These factorizations must be consistent with each other, even after resumming logarithms, allowing us to derive results on the other factorizations from ones more amiable to calculation. Other observables sensitive to rapidity divergences can also be treated this way, perhaps simplifying certain two or even three-loop calculations. Beyond perturbation theory, what is particularly appealing about this type of multi-differential factorization is that it gives an in-principle non-perturbative regularization of the rapidity divergences, something that hitherto was only feasible by tilting wilson-lines off the light-cone. For our specific processes, it suggests that the non-perturbative corrections to both transverse momentum resummation and threshold resummation coming from soft radiation are intimately connected, and work (see for instance [85][86][87]) on such corrections in transverse momentum resummation should be revisited in this light. Indeed, given the rapid development of integrability technology ( [88]) in planar N = 4 SYM for null polygonal wilson loops, it is reasonable to suppose that the fully differential soft function can be calculated exactly in that theory. In this example, the typical model for non-perturbative corrections -being inspired from renormalon related loop corrections probed by the running coupling [87,[89][90][91] -should fail, since the beta-function vanishes in that theory. The theory being conformal has no new structures arising in the deep infra-red.
Another benefit of this multiple region factorization is that it clarifies the structure of the transverse momentum spectrum in the threshold region, and in particular, how one can perform joint resummation using the techniques of [35,65]. Again we have two boundary theories, the threshold region and the TMD-PDF region, that are connected by an intermediate theory where functions refactorize. Evolving this intermediate theory to its natural scales would give a natural joint resummation formula that reduces to transverse momentum resummation or threshold resummation in the correct phase-space regions. Going beyond [71], we can also naturally include the full low-scale matrix elements consistently, giving a N 3 LL formula uniformly valid across phase-space. It would be interesting to see how such a resummation changes the transverse momentum spectrum, particularly at colliders with lower center of mass energies. which becomes nearly impossible when m approaches 15 or larger. Fortunately, the reduction of k + and k − to arbitrary powers can be done without resorting to IBP identities. To see it, we first insert some delta functions to separate the above integration into two steps, For the first integration, let's consider the following integral which is schematical representation of the actual integrand, where P 1,2,... , Q 1,2,... , K 1,2,... are linear combinations of loop and emission momenta, and i 1,2,... , j 1,2,... , k 1,2,... are the indices of propagators. The second line follows from dimensional analysis and rescaling dependence on n andn. F(d) is a function of the space-time dimension d that can only be obtained through direct calculation; l is given by the number of loops n l and soft emissions n e , l = (n e (d − 2) + n l d)/2 17 . It is straightforward to write down the result of the second integration, where B is the standard Beta function, B(x, y) = Γ(x)Γ(y)/Γ(x + y). So it becomes clear that we can write any integral with non-zero a or b in term of the corresponding integral with zero a and b by multiplying it the proper factor, . . = (2k 0 ) a+b B(a + l − (i 1 + i 2 + . . . ) − (k 1 + k 2 + . . . ), b + l − (j 1 + j 2 + . . . ) − (k 1 + k 2 + . . . )) B(l − (i 1 + i 2 + . . . ) − (k 1 + k 2 + . . . ), l − (j 1 + j 2 + . . . ) − (k 1 + k 2 + . . . )) which essentially realizes the reduction of insertions of k + and k − of arbitrary powers. As for the reduction of k 2 to arbitrary powers, we find no alternative but relying on the full 17 In practice, inverse unitarity allows us to treat δ(k 2 1,2,... ) in the phase space integration as propagators 1/k 2 1,2,... , and the definition of m becomes l = (n e + n l )d/2 accordingly.