Subleading Power Rapidity Divergences and Power Corrections for $q_T$

A number of important observables exhibit logarithms in their perturbative description that are induced by emissions at widely separated rapidities. These include transverse-momentum ($q_T$) logarithms, logarithms involving heavy-quark or electroweak gauge boson masses, and small-$x$ logarithms. In this paper, we initiate the study of rapidity logarithms, and the associated rapidity divergences, at subleading order in the power expansion. This is accomplished using the soft collinear effective theory (SCET). We discuss the structure of subleading-power rapidity divergences and how to consistently regulate them. We introduce a new pure rapidity regulator and a corresponding $\overline{\rm MS}$-like scheme, which handles rapidity divergences while maintaining the homogeneity of the power expansion. We find that power-law rapidity divergences appear at subleading power, which give rise to derivatives of parton distribution functions. As a concrete example, we consider the $q_T$ spectrum for color-singlet production, for which we compute the complete $q_T^2/Q^2$ suppressed power corrections at $\mathcal{O}(\alpha_s)$, including both logarithmic and nonlogarithmic terms. Our results also represent an important first step towards carrying out a resummation of subleading-power rapidity logarithms.

There has been significant interest and progress in studying power corrections [52][53][54][55][56][57][58][59] both in the context of B-physics (see e.g. refs. [60][61][62][63][64][65][66][67][68][69]) and for collider-physics cross sections (see e.g. refs. ). Recently, progress has been made also in understanding the behaviour of matrix elements in the subleading soft and collinear limit [92] in the presence of multiple collinear directions using spinor-helicity formalism. In ref. [93] the first all-order resummation at subleading power for collider observables was achieved for a class of powersuppressed kinematic logarithms in thrust including both soft and collinear radiation. More recently in ref. [94] subleading power logarithms for a class of corrections in the threshold limit have also been resummed. In both cases the subleading power logarithms arise from widely separated virtuality scales, and their resummation make use of effective field theory techniques. Given the importance of observables involving nontrivial rapidity scales, it is essential to extend these recent subleading-power results to such observables, and more generally, to understand the structure of rapidity logarithms and their evolution equations at subleading power.
In this paper, we initiate the study of rapidity logarithms at subleading power, focusing on their structure in fixed-order perturbation theory. We show how to consistently regularize subleading-power rapidity divergences, and highlight several interesting features regarding their structure. In particular, power-law divergences appear at subleading power, which give nontrivial contributions and must be handled properly. We introduce a new "pure rapidity" regulator and an associated "pure rapidity" MS-like renormalization scheme. This procedure is homogeneous in the power expansion, meaning that it does not mix different orders in the power expansion, which significantly simplifies the analysis of subleading power corrections. We envision that it will benefit many applications.
As an application of our formalism, we compute the complete O(α s ) power-suppressed contributions for q T for color-singlet production, which provides a strong check on our regularization procedure. We find the interesting feature that the appearing power-law rapidity divergences yield derivatives of PDFs in the final cross section. Our results provide an important ingredient for improving the understanding of q T distributions at next-to-leading power (NLP). They also have immediate practical applications for understanding and improving the performance of fixed-order subtraction schemes based on the q T observable [95].
To systematically organize the power expansion, we use the soft collinear effective theory (SCET) [96][97][98][99], which provides operator and Lagrangian based techniques for studying the power expansion in the soft and collinear limits. The appropriate effective field theory for observables with rapidity divergences is SCET II [100]. In this theory, rapidity logarithms can be systematically resummed using the rapidity renormalization group (RRG) [8,9] in a similar manner to virtuality logarithms. The results derived here extend the rapidity renormalization procedure to subleading power, and we anticipate that they will enable the resummation of rapidity logarithms at subleading power.
The outline of this paper is as follows. In sec. 2, we give a general discussion of the structure and regularization of rapidity divergences at subleading power. We highlight the issues appearing for rapidity regulators that are not homogeneous in the power-counting parameter, focusing on the η regulator as an explicit example. We then introduce and discuss the pure rapidity regulator, which is homogeneous. In sec. 3, we derive a master formula for the power corrections to the color-singlet q T spectrum at O(α s ), highlighting several interesting features of the calculation. We also give explicit results for Higgs and Drell-Yan production, and perform a numerical cross check to validate our results. We conclude in sec. 4.

Rapidity Divergences and Regularization at Subleading Power
Rapidity divergences naturally arise in the calculation of observables sensitive to the transverse momentum of soft emissions. In a situation where we have a hard interaction scale Q and the relevant transverse momentum k T of the fields is small compared to that scale, λ ∼ k T /Q 1, the appropriate effective field theory (EFT) is SCET II [100], which contains modes with the following momentum scalings n−collinear : k n ∼ Q (λ 2 , 1, λ) =⇒ k − /Q ∼ 1 , (2.1) Here we have used lightcone coordinates (n · k,n · k, k ⊥ ) ≡ (k + , k − , k ⊥ ), defined with respect to two lightlike reference vectors n µ andn µ . For concreteness, we take them to be n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1). Unlike SCET I where the modes are separated in virtuality, in SCET II the modes in the EFT have the same virtuality, but are distinguished by their longitudinal momentum (k + or k − ), or equivalently, their rapidity e 2y k = k − /k + . This separation into modes at hierarchical rapidities introduces divergences, which arise when k + /k − → ∞ or k + /k − → 0 [9,[101][102][103][104]. These so-called rapidity divergences are not regu-lated by dimensional regularization, which is boost invariant and therefore cannot distinguish modes that are only separated in rapidity. Rapidity divergences can be regulated by introducing a rapidity regulator that breaks boost invariance, allowing the modes to be distinguished, and logarithms associated with the different rapidity scales to be resummed. The rapidity divergences cancel between the different sectors of the effective theory, since they are not present in the full theory. They should not be thought of as UV, or IR, but as arising from the factorization in the EFT. By demanding invariance with respect to the regulator, one can derive renormalization group evolution equations (RGEs) in rapidity. In SCET, a generic approach to rapidity evolution was introduced in refs. [8,9]. These rapidity RGEs allow for the resummation of large logarithms associated with hierarchical rapidity scales.
At leading power in the EFT expansion, the structure of rapidity divergences and the associated rapidity renormalization group are well understood by now, and they have been studied to high perturbative orders (see e.g. ref. [105] at three-loop order). Indeed, in certain specific physical situations involving two lightlike directions, rapidity divergences can be conformally mapped to UV divergences [104,[106][107][108][109], giving a relation between rapidity anomalous dimensions and standard UV anomalous dimensions. However, little is known about the structure of rapidity divergences or their renormalization beyond the leading power. 1 In this section, we discuss several interesting features of rapidity divergences at subleading power, focusing on the perturbative behavior at next-to-leading order (NLO). At subleading power there are no purely virtual corrections at NLO, and so we will focus on the case of the rapidity regularization of a single real emission, which allow us to identify and resolve a number of subtleties. After a brief review of the structure of rapidity-divergent integrals at leading power in sec. 2.1, we discuss additional issues that arise at subleading power in sec. 2.2. We discuss in detail the behavior of the η regulator at subleading power, highlighting effects that are caused by the fact that it is not homogeneous in the power expansion. In sec. 2.3, we introduce the pure rapidity regularization, which regulates rapidity instead of longitudinal momentum and which we find to significantly simplify the calculation at subleading power. Finally, in sec. 2.4, we discuss the distributional treatment of power-law divergences, which arise at subleading power.

Review of Rapidity Divergences at Leading Power
We begin by reviewing the structure of rapidity divergent integrals at leading power. As mentioned above, we restrict ourselves to the case of a single on-shell real emission, which suffices at NLO. Defining δ + (k 2 ) = θ(k 0 )δ(k 2 ), its contribution to a cross section sensitive to the transverse momentum k T of the emission is schematically given by 1 For some interesting recent progress for the particular case of the subleading power Regge behavior for massive scattering amplitudes in N = 4 super Yang-Mills theory, see ref. [110].
Here, we have extracted the overall 1/k 2 T behaviour, and g(k) is an observable and process dependent function, containing the remaining phase-space factors and amplitudes. The precise form of g(k) is unimportant, except for the fact that it includes kinematic constraints on the integration range of k ± , For our discussion we take k T > 0 such that we can work in d = 4 dimensions. In the full theory, eq. (2.2) is finite, with the apparent singularities for k ± → 0 or k ± → ∞ being cut off by the kinematic constraints in eq. (2.3). In the effective theory, one expands eq. (2.2) in the soft and collinear limits specified in eq. (2.1). This expansion also removes the kinematic constraints, k ± min → 0 soft and collinear limits , such that individual soft and collinear contributions acquire explicit divergences as k ± → 0 or k ± → ∞. This is actually advantageous, since the associated logarithms can now be tracked by these divergences. To regulate them, we introduce a regulator R(k, η), where η is a parameter such that lim η→0 R(k, η) = 1. By construction, inserting R(k, η) under the integral in eq. (2.2) does not affect the value of dσ( k T ) when taking η → 0 in the full calculation. To describe the limit k T Q, we expand eq. (2.2) in the soft and collinear limits described by the modes in eq. (2.1). To be specific, the soft limit of eq. (2.2) is obtained by evaluating the integrand together with the regulator R(k, η) using the soft scaling k s of eq. (2.1), and expanding in λ, Since the leading-power result must scale like 1/k 2 T , the LP soft limit g s (k µ = 0) must be a pure constant, which implies that the kinematic constraints in eq. (2.3) are removed. This introduces the aforementioned divergences as k − → 0 or k − → ∞, which are now regulated by R(k, η).
The analogous expansion in the collinear sectors is obtained by inserting the k n or kn scalings of eq. (2.1) into eq. (2.2), and expanding in λ, In this case, only the lower bound on k ± is removed by the power expansion, while the upper limit is given by the relevant hard scale Q. The expansion of g(k n ) in the collinear limit can still depend on the momentum k − /Q ∼ O(λ 0 ), as indicated by the functional form of g n (k − /Q), and likewise for then-collinear limit. Without the rapidity regulator, the integrals in eqs. (2.5) and (2.6) exhibit a logarithmic divergence as k ± → 0 or k ± → ∞, which is not regulated by dimensional regularization or any other invariant-mass regulator. Since k + k − = k 2 T is fixed by the measurement, this corresponds to a divergence as the rapidity y k = (1/2) ln(k − /k + ) → ±∞. The rapidity regulator R(k, η) regulates these divergence by distinguishing the soft and collinear modes. To ensure a cancellation of rapidity divergences in the effective theory, it should be defined as a function valid on a full-theory momentum k, which can then be expanded in the soft or collinear limits. Since there are no divergences in the full theory, this guarantees the cancellation of divergences in the EFT expansion.
At leading power a variety of regulators have been proposed. Since the divergences are only logarithmic, and the focus has not been on higher orders in the power expansion, there are not many constraints from maintaining the power counting of the EFT. Therefore, a variety of regulators have been used, including hard cutoffs [47,48,50,102], tilting Wilson lines off the lightcone [111], the delta regulator [112], the η regulator [8,9], the analytic regulator [41,113,114], and the exponential regulator [115].
At subleading power, we will discuss in more detail the application of the η regulator, which can be formulated at the operator level by modifying the Wilson lines appearing in the SCET fields as [8,9] where S n and W n are soft and collinear Wilson lines. The operator P picks out the large (label) momentum flowing into the Wilson line, ν is a rapidity regularization scale, η a parameter exposing the rapidity divergences as 1/η poles, and w a bookkeeping parameter obeying Note that at leading power, one can replace |2P z | → |n · P| in eq. (2.8), as employed in refs. [8,9], while at subleading power we will show that this distinction is actually important.
The η regulator was extended in ref. [116] to also regulate Glauber exchanges in forward scattering, where regulating Wilson lines alone does not suffice.

Rapidity Regularization at Subleading Power
We now extend our discussion to subleading power, where we will find several new features. First, while at leading power, rapidity divergences arise only from gluons, at subleading power rapidity divergences can arise also from soft quarks. Soft quarks have also been rapidityregulated to derive the quark Regge trajectory [117]. Here, since we consider only the case of a single real emission crossing the cut, this simply means that we must regulate both quarks and gluons. More generally, one would have to apply a rapidity regulator to all operators in the EFT, as has been done for the case of forward scattering in ref. [116]. It would be interesting to understand if these subleading rapidity divergences can also be conformally mapped to UV divergences of matrix elements, as was done for the rapidity divergences in the leading power q T soft function in refs. [104,109]. Second, the structure of rapidity divergences becomes much richer at subleading power, placing additional constraints on the form of the rapidity regulator to maintain a simple power expansion. This more interesting divergence structure follows directly from power counting. For example, the subleading corrections to the soft limit can be obtained by expanding the integrand in eq. (2.5) to higher orders in λ. The power counting for soft modes in eq. (2.1) implies that the first O(λ) power suppression can only be given by additional factors of k − /Q or k + /Q in eq. (2.5). At the next order, O(λ 2 ), one can encounter additional factors (k + /Q) 2 , (k − /Q) 2 . The possible structure of rapidity-divergent integrals in the soft limit up to O(λ 2 ) is thus given by 2 where it is understood that k + = k 2 T /k − . We can see that the O(λ 0 ) limit only produces logarithmic divergences, while the power-suppressed corrections give rise to power-law divergences. The prototypical rapidity-divergent integral encountered in the soft limit is thus given by where α counts the additional powers of k − . 2 We can also have integrals with an additional factor of kT /Q or k 2 T /Q 2 , which however do not change the structure of the integrand and can thus be treated with the same techniques as at leading power.
A similar situation occurs in the collinear sectors. In the n-collinear limit, k ∼ Q(λ 2 , 1, λ), the large momentum k − is not suppressed with respect to Q, such that the power suppression can only arise from explicit factors of k 2 T . (Of course, k + ∼ O(λ 2 ) can also give a suppression, but it can always be reduced back to k + = k 2 T /k − .) Similarly, in then-collinear limit k + is unsuppressed, and power suppressions only arise from k 2 T . However, the structure of the collinear expansion of g(k) is richer than in the soft case, because there is always a nontrivial dependence on the respective unsuppressed ratio k ∓ /Q. To understand this intuitively, consider the splitting of a n-collinear particle into two on-shell n-collinear particles with momenta (2.12) The associated Lorentz-invariant kinematic variable is given by Expanding any function of s 12 in k T thus gives rise to additional factors of the large momentum k − . Thus, in general, expanding g(k n ) in the collinear limit can give rise to both positive and negative powers of k − that accompany the power-suppression in k 2 T . These factors are of course not completely independent, as the sum of all soft and collinear contributions must be rapidity finite, i.e., any rapidity divergences induced by these additional powers of k − must in the end cancel against corresponding divergences in the soft and/or other collinear contributions. In summary, the generic form of integrals in the collinear expansion is given by 14) Here, g n (x) and gn(x) are regular functions as x → 0. At LP, only α = 0 contributes, which gives rise to logarithmic divergences, while at subleading power for α = 0 we again encounter power-law divergences. As we will see in sec. 2.4, these power-law divergences have a nontrivial effect, namely they lead to derivatives of PDFs in the perturbative expansion for hadron collider processes. The presence of power-law divergences at subleading power also implies that more care must be taken to ensure that the regulator does not unnecessarily complicate the power counting of the EFT. For example, with the exponential regulator [115], or with a hard cutoff, power-law divergences lead to the appearance of powers of the regulator scale, and hence break the homogeneity of the power expansion of the theory.
Furthermore, at leading power one also has the freedom to introduce and then drop subleading terms to simplify any stage of the calculation. While this may seem a general feature and not appear very related to the regularization of rapidity divergences, we will see in a moment that this freedom, explicitly or not, is actually used in most of the rapidity regulators in the literature.
In summary, having a convenient-to-use regulator at subleading power imposes stronger constraints than at leading power. In particular, we find that the regulator • must be able to regulate not only Wilson lines, but all operators, including those generating soft quark emissions, • must be able to deal not only with logarithmic divergences, but also with power-law divergences without violating the power counting of the EFT by inducing power-law mixing, • and should be homogeneous in the power-counting parameter λ to minimize mixing between different powers.
The first requirement means one cannot use regulators acting only on Wilson lines, such as taking Wilson lines off the light-cone as in ref. [111], the δ regulator as used in refs. [13,112], and the η regulator as used in refs. [8,9], while the η regulator as modified and employed in refs. [116,117] and the analytic regulator of ref. [114] can be used. The second requirement is satisfied by all dimensional regularization type regulators, such as the η regulator or analytic regulator, but not by those that are more like a hard cutoff, including the exponential regulator [115]. To highlight the last point, in the following we discuss in more detail the properties of the η regulator at subleading power.

The η Regulator at Subleading Power
In the η regulator, one regulates the k z momentum of emissions through the regulator function (see eq. (2.7)) For a single massless emission this corresponds to regulating its phase-space integral as (2.17) In the soft limit k + ∼ k − ∼ λQ, the regulator is homogeneous in λ and therefore does not need to be expanded. The prototypical soft integral in eq. (2.11) evaluates to Symmetry under α ↔ −α implies that This reflects the symmetry under exchanging k − ↔ k + , which is not broken by the η regulator. One can easily deduce the behavior as η → 0 from eq. (2.18). Since sin(η) ∼ η, a pole in η can only arise if both Γ functions have poles, which requires α = 0. A finite result is obtained if exactly one Γ function yields a pole, which requires α to be even. For odd α, the expression vanishes at η = 0. Hence, the exact behavior for η → 0 is given by In particular, since the η regulator behaves like dimensional regularization, it is well-behaved for power-law divergences and the soft integrals only give rise to poles from the logarithmic divergences.
In the collinear sector, the behavior is more complicated at subleading power, because the regulator factor 2k z = k − − k + is not homogeneous in λ. At leading power [8,9,116], one takes advantage of the fact that 2k z → k − in the n-collinear limit and 2k z → k + in thencollinear limit, so that the expanded result correctly regulates the collinear cases, and makes it symmetric under the exchange n ↔n. A fact that will be important for our analysis is that this power expansion induces higher order terms. These terms have never been considered in the literature since they are not important at leading power. However, at subleading power one can no longer neglect the subleading component of the regulator. Implementing the η regulator at subleading power in the collinear limits thus requires to expand the regulator eq. (2.16) itself, Applying this to the general LP integral in the n-collinear sector, eq. (2.14) with α = 0, we obtain and analogously for I Here, the first line is the standard LP integral, while the second line arises from expanding the regulator and is suppressed by While it is also proportional to η, the remaining integral can produce a 1/η rapidity divergence to yield an overall finite contribution.
In sec. 3, we will see explicitly that these terms from expanding the regulator are crucial to obtain the correct final result at subleading power. However, in practice they are cumbersome to track in the calculation and yield complicated structures. To establish an all-orders factorization theorem, the mixing of different orders in the power expansion due to the regulator becomes a serious complication. Hence, it is desirable to employ a rapidity regulator that is homogeneous in λ. We will present such a regulator in the following sec. 2.3.

Pure Rapidity Regularization
We wish to establish a rapidity regulator that is homogeneous at leading power such that it does not mix LP and NLP integrals, as observed in sec. 2.2.1 for the η regulator. This can be achieved by implementing the regulator similar to the η regulator of refs. [8,9,116], but instead of regulating the momentum k z with factors of w|2k z /ν| −η/2 , one regulates the rapidity y k of the momentum k µ , where To implement a regulator involving rapidity we use 3 factors of Here we have defined a rapidity scale υ (\upsilon) which is the analog of the scale ν (\nu) in the η regulator. Although υ is dimensionless, in contrast to the dimensionful ν, it still shares the same properties as pure dimensional regularization. In particular, it will give rise to poles in η that can be absorbed in MS-like rapidity counterterms. To ensure υ independence of eq. (2.26), we introduced a bookkeeping parameter w = w(υ) in analogy to the bookkeeping parameter w(ν) in the η regulator, see eq. (2.9) and ref. [9].
3 Note that we can implement the pure rapidity regulator in terms of label and residual momentum operators for example as where the label momentum operator P picks out the large O(λ 0 ) momentum component of the operator it acts on, while ∂ picks out the O(λ) or O(λ 2 ) components. In this case, the operator picks out the rapidity of the operator it acts on.
We call eq. (2.26) the pure rapidity regulator, and pure rapidity regularization the procedure of regulating rapidity divergences using eq. (2.26). When only the 1/η poles are subtracted we then refer to the renormalized result as being in the pure rapidity renormalization scheme.
If we want to make the rapidity scale υ into a true rapidity scale Υ, then we can change variables as (2.27) With this definition eq. (2.26) becomes 28) and the factor regulating divergences depends on a rapidity difference between the scale parameter Υ and y k . It is interesting to consider the behavior of amplitudes regulated with eq. (2.28) under a reparameterization transformation known as RPI-III [52], which takes n µ → e −β n µ and n µ → e βnµ for some, not necessarily infinitesimal, constant β. For a single collinear sector, this can be interpreted as a boost transformation. Since RPI transformations can be applied independently for each set of collinear basis vectors {n i ,n i } they in general constitute a broader class of symmetry transformations in SCET. Prior to including a regulator for rapidity divergences all complete SCET amplitudes are invariant under such transformations. All previous rapidity regulators violate this symmetry. For the pure rapidity regulator in eq. (2.28) we have y k → y k + β, so the transformation is quite simple. 4 It can be compensated by defining the rapidity scale to transform like a rapidity, Υ → Υ+β. Therefore, the υ η factor in the regulator does for RPI-III what the usual µ factor does for the mass-dimensionality in dimensional regularization.
As an example of the application of this new regulator, we consider again a real emission with momentum k µ . The regulator function R(k, η) that follows from eq. (2.26) is given by The real-emission phase space is then regulated as A peculiar feature of the pure rapidity regulator is that it renders the prototypical soft integrals scaleless such that they vanish. That is, using eq. (2.29) in eq. (2.11), we obtain The final integrals are scaleless and vanish for all integer values of α, just like scaleless integrals vanish in dimensional regularization. 5 Considering the collinear sectors, the prototypical collinear integrals in eq. (2.14) with R Y (k, η) become Although the regulator does not act symmetrically in the n-collinear andn-collinear sectors, the asymmetry is easy to track by taking η ↔ −η and υ ↔ 1/υ when swapping n ↔n and is homogeneous in λ, it does not generate any subleading power terms, in contrast to eq. (2.22) for the η regulator. In particular, the LP integral becomes where we used the standard distributional identity 1/x 1+η = −δ(x)/η+L 0 (x)+O(η) to extract the 1/η divergence. (See sec. 2.4 below for a more general discussion.) Taking η → −η, the analogous 1/η pole in then-collinear sector has the opposite sign, such that the 1/η poles cancel when adding the n-collinear andn-collinear contributions. This is a general feature in all cases where the soft contribution vanishes as in eq. (2.31). Some comments about the features of the pure rapidity regulator are in order: • It involves the rapidity and therefore breaks boost invariance as required to regulate rapidity divergences. The boost invariance is restored by the dimensionless υ rapidity scale, analogous to how the dimensionful mass scale µ in dimensional regularization restores the dimensionality.
• Rapidity divergences appear as 1/η poles, allowing the definition of the pure rapidity renormalization scheme as a dimensional regularization-like scheme.
• At each order in perturbation theory, the poles in η and the υ-dependent pieces cancel when combining the results for the n-collinear,n-collinear, and soft sectors.
• The pure rapidity regulator is homogeneous 6 in the SCET power counting parameter λ. Therefore it does not need to be power expanded, and hence does not mix contributions at different orders in the power expansion.
• For the case of a single real emission considered here: -Soft integrals and zero-bin [102] integrals are scaleless and vanish.
-It follows that the η poles and the υ dependent pieces cancel between the n-collinear andn-collinear sectors.
-The results for the n-collinear andn-collinear sectors are not identical but are trivially related by taking η ↔ −η and υ ↔ 1/υ when swapping n ↔n.
The introduction of this new pure rapidity regulator allows us to regulate rapidity divergences at any order in the EFT power expansion, while maintaining the power counting of the EFT independently at each order.
Although in this paper we will only use pure rapidity regularization for a single real emission at fixed order, we note that one can derive a rapidity renormalization group for the pure rapidity regulator by imposing that the cross section must be independent of υ. Similar to the η regulator, this regulator is not analytical and can also be used to properly regulate virtual and massive loops. This will be discussed in detail elsewhere.
To conclude this section we note that the pure rapidity regulator can be seen as a particular case of a broader class of homogeneous rapidity regulators given by where c = 1 is an arbitrary parameter governing the antisymmetry between the n-collinear andn-collinear sectors. As for the pure rapidity regulator, this regulator is homogeneous in λ and renders the same class of soft integrals scaleless. However, it requires an explicit dimensionful scale ν to have the correct mass dimension. Note that for c = 1, eq. (2.35) only depends on the boost invariant product k + k − and therefore does not regulate rapidity divergences. For c = −1, it recovers the pure rapidity regulator and the dependence on ν cancels. Lastly, for c = 0 and massless real emissions, eq. (2.35) essentially reduces to the regulator of ref. [114].

Distributional Treatment of Power Law Divergences
To complete our treatment of rapidity divergences at subleading power, we show how their distributional structure can be consistently treated when expanded against a general test function. In particular, we will see that the power-law rapidity divergences lead to derivatives of PDFs. In the collinear limit at NLP, we obtain divergent integrals of the form which appear for both the η regulator (with a = 1 − α = 1, 2, 3) and the pure rapidity regulator (with a = 1 − α = 1, 2). The function g n (k − /Q) is defined to be regular for k − /Q → 0. If it is known analytically, we can in principle evaluate the integral in eq. (2.36) analytically and expand the result for η → 0 to obtain the regularized expression. However, g n (k − /Q) is typically not given in analytic form. In particular, for pp collisions it contains the parton distribution functions (PDFs) f (x). Therefore, to extract the rapidity divergence, we need to expand 1/(k − ) a+η in η in a distributional sense. To do so, we first change the integration variable from k − to the dimensionless variable z defined through (2.37) In eq. (2.37), the rapidity divergence arises as z → 1. For a = 1, it can be extracted using the standard distributional identity where L 0 (y) = [θ(y)/y] + is the standard plus distribution and we remind the reader that its convolution against a test functiong(z) is given by For a > 1, these distributions need to be generalized to higher-order plus distributions subtracting higher derivatives as well. For example, for a = 2 one obtains where the second-order plus function L ++ 0 (1−z) regulates the quadratic divergence 1/(1−z) 2 . Its action on a test functiong(z) is given by a double subtraction, In appendix B, we give more details on these distributions, generalizing to arbitrary a ≥ 1.
Note that the second-order plus function has also appeared for example in ref. [118]. Eq. (2.40) implies the appearance of derivatives of delta functions, δ (1 − z), which will induce derivatives of the PDFs that are contained ing(z). The appearance of such derivatives in subleading power calculations was first shown in ref. [78] in the context of SCET I -like observables. However, in such cases they arose simply from a Taylor expansion of the momentum being extracted from the PDF. Here, they also arise from power-law divergences, a new mechanism to induce derivatives of PDFs. Recently, power-law divergences inducing derivatives of PDFs have appeared also in the study of SCET I -like observables involving multiple collinear directions at subleading power [92]. We believe they are a general feature of calculations beyond leading power.
In practice, the higher-order distributions can be cumbersome to work with. Instead, we find it more convenient to use integration-by-parts relations to reduce the divergence in eq. (2.37) to the linear divergence 1/(1 − z), which yields explicit derivatives of the test function. For the cases a = 2 and a = 3 we encounter in sec. 3, this gives Equations (2.42) and (2.43) can be used to write the kernels fully in terms of a standard L 0 , but they must be applied within the integral to directly yield derivatives of the test functioñ g(z). In our application in sec. 3,g(z) will always involve the PDF f (x/z) and vanish at z = x. We can thus also write eqs. (2.42) and (2.43) as operator equations, Note that the second relation is quite peculiar, as we have to add the boundary term proportional to g (x), and thus cannot be interpreted as a distributional relation. In our calculation in sec. 3, this term will not contribute due to an overall suppression by η, such that only the divergent term in eq. (2.45) needs to be kept.

Power Corrections for Color-Singlet q T Spectra
In this section we use our understanding of rapidity regularization at subleading power to compute the perturbative power corrections to the transverse momentum q T in color-singlet production at invariant mass Q, which is one of the most well studied observables in QCD. Schematically, the cross section differential in q T can be expanded as where σ (0) is the leading-power cross section and σ (2n) the N n LP cross section. These terms scale like and hence only the LP cross section is singular as q T → 0. In particular, σ (0) contains Sudakov double logarithms log 2 (Q/q T ). The factorization of σ (0) in terms of transverse-momentum dependent PDFs (TMDPDFs) was first shown by Collins, Soper, and Sterman in refs. [2][3][4] and later elaborated on by Collins in ref. [111]. Its structure was also studied in refs. [119][120][121]. The factorization was also studied in the framework of SCET by various groups, see e.g. refs. [9,12,13]. Using the notation of ref. [9], the factorized LP cross section for the production of a color-singlet final state L with invariant mass Q and total rapidity Y in a proton-proton collision can be written as 7 where x a,b = Qe ±Y /E cm are the momentum fractions carried by the incoming partons. In eq. (3.3), H ij is the hard function describing virtual corrections to the underlying hard process ij → L, theB i are TMD beam functions in Fourier space andS is the TMD soft function in Fourier space. While H ij only depends on the MS renormalization scale µ, the beam and soft functions also depend on the rapidity renormalization scale ν.
Recently, there has been some progress towards a nonperturbative factorization of the NLP cross section dσ (2) /dq 2 T , which involves higher twist PDFs [81,84]. Here, we are interested in studying the perturbative power corrections to the NLP terms, where one can perform an OPE to match onto standard PDFs. At subleading power, the perturbative kernels also involve (higher) derivatives of distributions, which can always be reduced to standard distributions acting on derivatives of PDFs. The NLP cross section at O(α n s ) thus takes the form whereσ LO is the LO partonic cross section which serves as an overall normalization. The C (2,n) ab are perturbative coefficients, expressed in terms of distributions, and we suppress the explicit Q and Y dependence in the kernels C (2,n) ab . Their general logarithmic structure is More explicitly, at NLO they have the form i.e. they only contain a single logarithm ln(Q 2 /q 2 T ) and a q T -independent piece. (Note that due to the dependence on z a,b , it will yield a Q 2 and Y dependence.) We emphasize that in the form given here, all logarithms have been extracted, and the q T distribution is directly expressed in terms of PDFs and their derivatives.
In the following, we will derive a master formula to obtain the NLO NLP kernels C (2,1) ab for arbitrary color-singlet processes, as well as the explicit results for Higgs and Drell-Yan production. The study of higher perturbative orders, and the derivation of a factorization and resummation is left to future work. However, we do wish to comment on one complication which occurs for q T at higher orders, that we have not addressed. Unlike for beam thrust, at NNLO and beyond, one can have power-suppressed contributions at small q T from two hard partons in the final state that are nearly back-to-back such that their transverse momenta balance to give a small total q T . At NNLO, this is at most a constant power correction, since it is not logarithmically enhanced. but at higher orders it can have a logarithmic contribution. These power corrections are of a different nature than those discussed here, and are not captured as an expansion in the soft and collinear limits about the Born process. The remainder of this section is organized as follows. In sec. 3.1, we derive the master formula for the NLP corrections using the η regulator, showing in particular that the terms from expanding the regulator contribute. In sec. 3.2, we rederive this master formula in pure rapidity regularization, which will be simpler due to the fact that one does not have additional terms from the expansion of the regulator, and due to the fact that the soft sector is scaleless. In sec. 3.3, we then apply the master formula to derive explicit results for Drell-Yan and gluon-fusion Higgs production. In sec. 3.4, we discuss our results and compare them with the known NLP results for beam thrust. Finally in sec. 3.5, we provide a numerical validation of our results.

Master Formula for Power Corrections to Next-to-Leading Power
We consider the production of a color-singlet final state L at fixed invariant mass Q and rapidity Y , measuring the magnitude of its transverse momentum q 2 T = | q T | 2 . The underlying partonic process is where a, b are the incoming partons and X denotes additional QCD radiation. Following the notation of ref. [90], we express the cross section as dσ dQ 2 dY dq 2 Here, the incoming momenta are given by

10)
k = i k i is the total outgoing hadronic momentum, and q is the total leptonic momentum.
In particular, k T = i k i,T is the vectorial sum of the transverse momenta of all emissions. Since the measurements are not affected by the details of the leptonic final state, the leptonic phase-space integral has been absorbed into the matrix element, The matrix element M also contains the renormalization scale µ 2 , as usual associated with the renormalized coupling α s (µ), and may also contain virtual corrections. There is an important subtlety when measuring the transverse momentum q T using dimensional regularization, as the individual transverse momenta k i,T are continued to 2 − 2 dimensions. The measurement function δ(q 2 T − | k T | 2 ) in eq. (3.9) can thus be interpreted either as measuring the magnitude in 2 − 2 dimensions or the projection onto 2 dimensions. This scheme dependence cancels in the final result, but can lead to different intermediate results. At the order we are working, both choices give identical results, so for simplicity of the following manipulations we specify to measuring the magnitude in 2 − 2 dimension. For detailed discussions, see e.g. refs. [127,128].
The δ functions measuring the invariant mass Q and rapidity Y fix the incoming momenta to be Equation (3.9) can now be simplified to dσ dQ 2 dY dq 2 where we introduced the abbreviation (3.14) This emphasizes that the squared matrix element depends only on the Born measurements Q and Y , which fix the incoming momenta through eqs. (3.10) and (3.12), and the emission momenta k i . The restriction that ζ a,b ∈ [0, 1] is kept implicit in eq. (3.13) through the support of the proton PDFs.

General Setup at NLO
For reference, we start with the LO cross section following from eq. (3.13), and A LO is the squared matrix element in the Born kinematics, see eq. (3.14). For future reference, we also define the LO partonic cross section,σ LO (Q, Y ), by At NLO, the virtual correction only contributes at leading power and is proportional to δ(q 2 T ). At subleading power, it suffices to consider the real correction, given from eq. (3.13) by In the following, we will mostly keep the symbol k + often leaving the use of the relation k + = k 2 T /k − to the end, since this makes the symmetry under k + ↔ k − manifest. The integral in eq. (3.18) is finite as the physical support of the PDFs, 0 ≤ ζ a,b ≤ 1, cuts off the integral in k − . As discussed in sec. 2.1, these constraints will be expanded for small q T Q, after which the integral becomes rapidity divergent. To regulate the integral, we use the η regulator where one inserts a factor of w 2 |2k z /ν| −η into the integral, We now wish to expand eq. (3.19) in the limit of small λ ∼ q T /Q 1. Using the knowledge from the EFT, this can be systematically achieved by employing the scaling of eq. (2.1), for the momentum k. By inserting each of these scalings into eq. (3.19) and expanding the resulting expression to first order in λ, one precisely obtains the soft and beam functions as defined in the η regulator. This illustrative exercise is shown explicitly in appendix A. Here, we are interested in the first nonvanishing power correction, which occurs at O(λ 2 ) ∼ O(q 2 T /Q 2 ). We will explicitly show that the O(λ) linear power correction vanishes. To compute the O(λ 2 ) result, we will consider the soft and collinear cases separately, deriving master formulas for all scalings applicable to any color-singlet production.

Soft Master Formula for q T
We first consider the case of a soft emission k ∼ Q(λ, λ, λ). In this limit, the incoming momenta from eq. (3.12) are expanded as where as usual k + = k 2 T /k − , x a,b = Qe ±Y /E cm as in eq. (3.16), and the terms in square brackets correspond to O(λ 0 ), O(λ 1 ), and O(λ 2 ), respectively. It follows that the PDFs and flux factor are expanded as Here, (sym.) denotes simultaneously flipping a ↔ b and letting k − → k + , Y → −Y . For brevity, we introduced the abbreviation Φ (n) for the O(λ n ) pieces. Note that we expanded to the second order in λ, as the O(λ 1 ) piece will vanish and the first nonvanishing correction in fact arises at O(λ 2 ). The expansion of the matrix element is process dependent, and we define the expansion in the soft limit through The LP matrix element scales as A s ∼ λ 0 . The next two matrix elements are each suppressed by an additional order in λ relative to the one before.
Plugging the expansions eqs. (3.22) and (3.23) back into eq. (3.19) and collecting terms in λ, the soft limit through O(λ 2 ) is obtained as The first term in curly brackets is the leading-power result, the second term the O(λ) contribution, and the last line contains the O(λ 2 ) contribution. Since each of these terms has a homogeneous scaling in λ, they can only contribute integer powers of k − , yielding integrals of the form I

Leading Power [O(λ 0 )]
The leading soft limit of the squared amplitude A is universal and given by where µ MS is the renormalization scale in the MS scheme and C = C F , C A is the Casimir constant for the qq and gg channel, and the limit vanishes for any other channel. The cross section at LP thus becomes In sec. A.1, we use this to compute the known bare LP soft function at NLO as a cross check.

O(λ)
Here, we show that power corrections at O(λ) ∼ O(q T /Q) vanish at NLO. At this order, we can let → 0 to obtain the cross section from eq. (3.24) as From eq. (3.22), the expansion of the phase space is given by Hence, this contribution to eq. (3.27) is proportional to I s of the matrix element is suppressed by O(λ) relative to A LO , which from power counting can only be given by either k − or k + = k 2 T /k − . Hence, the Φ (0) A (1) term is also proportional to I (±1) s (R z ) = 0 and vanishes as well. More generally, power counting combined with the behavior of the integrals in eq. (2.20) shows that at NLO, the power expansion is in q 2 T /Q 2 . It would be interesting to extend this proof to higher perturbative orders. We also remark that the collinear limit will not have a O(λ) expansion at all, and thus the consistency condition that rapidity divergences cancel between soft and collinear sectors already implies that the soft NLP result cannot contribute to the leading logarithm.

Next-to-Leading Power [O(λ 2 )]
The first nonvanishing power correction thus arises at O(λ 2 ) ∼ O(q 2 T /Q 2 ). To derive a general master formula at this order, we decompose the expansion of the matrix element according to the possible dependence on k ± , which follows from power counting and mass dimension, The expansion is defined such that all A (i) have the same mass dimension. We now only need to plug eq. (3.29) back into eq. (3.24), collect the powers of k − (using that k + = k 2 T /k − ) and apply eq. (2.18). Only terms proportional to I (0) s (R z ) will yield a divergence in η, and thus constitute the LL correction at NLP, while all other terms contribute at NLL. We find dσ (2),LL s dQ 2 dY dq 2 and dσ (2),NLL s dQ 2 dY dq 2 exhibit an explicit rapidity dependence, which is surprising for the boost-invariant observable q T . In fact, we will see explicitly that the full soft expansion exactly cancels against rapiditydependent terms in the collinear expansions, yielding a rapidity-independent final result. This behavior is expected since the rapidity dependence arises from the rapidity-dependent regulator, and therefore we expect that they should cancel in the final regulator independent result.

Collinear Master Formula for q T
We next consider the case of a n-collinear emission k ∼ Q(λ 2 , 1, λ), from which one can easily obtain then-collinear case from symmetry. Here, it is important to consistently expand the rapidity regulator in eq. (3.19) in the n-collinear limit, Applying this to eq. (3.19) yields We now expand all pieces in λ. The incoming momenta from eq. (3.12) are expanded as where we grouped the terms of common scaling together and defined The superscript (2) denotes the suppression by λ 2 . Expanding the PDFs and flux factors in λ, we obtain The expansion of the matrix element is process dependent, and we define it by Note that in contrast to the soft limit, there is no O(λ) suppressed term here. Next, we switch the integration variable in eq. (3.33) via , (3.37) where the lower bound on the z a integral follows from the physical support of the PDF f a (x a /z a ). Inserting eqs. (3.35) -(3.37) into eq. (3.33) and collecting the O(λ 0 ) and O(λ 2 ) pieces, we obtain the leading n-collinear limit as The corresponding result in then-collinear case reads As discussed in sec. 2.4, a striking feature of eqs. (3.39) and (3.40) is the appearance of power divergences 1/(1 − z) 2+η and even 1/(1 − z) 3+η , which can be regulated using higherorder plus distributions, see also appendix B. Here, we find it more convenient to employ the integration-by-parts relations in eqs. (2.42) and (2.43) to write the kernels fully in terms of standard plus distributions, at the cost of inducing explicit derivatives of the PDFs. In order to apply these relations, we need to identify all divergences in 1/(1 − z) 2 and 1/(1 − z) 3 . To do so, first note that the LP matrix element scales as where P is the appropriate splitting function in d = 4 − 2 dimensions, which itself scales like P (z, ) ∼ 1/(1 − z). Due to the overall prefactor of k − ∼ (1 − z), the LP matrix element is finite as z → 1. Power counting implies that the subleading matrix element element can at most yield one additional pole 1/(1 − z). Motivated by these two observations, we write the expanded squared amplitude as 42) and likewise for An in then-collinear limit. The power suppression of A (2) n is made manifest by extracting the factor k 2 T /Q 2 . For brevity, we suppress any dependence of A Here, we used that the LL result is proportional to δ(1 − z a ) to cancel the z a integral in eq. (3.38), and the A Here, all terms with an explicit rapidity dependence arise from the expansion of the regulator itself, see eq. (3.32). In practice, they will exactly cancel against the soft NLL result eq. (3.31).

Derivation of the Master Formula in Pure Rapidity Regularization
In sec. 3.1, we used the η regulator of the form |2k z /ν| −η to derive the master formula. In this section, we repeat the derivation of the master formula using the pure rapidity regulator introduced in sec. 2.3. As discussed there, this regulator has the advantage that it is homogeneous in the power expansion, which reduces the number of terms at subleading power. Furthermore, it renders the soft sector scaleless. The result using the generalization of the pure rapidity regulator, eq. (2.35), is shown in appendix C for completeness. The derivation of the n-collinear expansion proceeds similar to the calculation shown in sec. 3.1.3. In eq. (3.39), one has to replace the regulator factor by and drop the terms in η/e 2Y , as they are fully induced by the expansion of the regulator. The NLP LL result is then easily obtained from eq. (3.43) by replacing ν → q T υ, In then-collinear limit, one has to replace the regulator factor and drop terms in η/e −2Y in eq. (3.40). The NLP LL result is then obtained from eq. (3.43) by replacing η → −η, ν → q T /υ and exchanging a ↔ b as dσ (2),LL n dQ 2 dY dq 2 Summing eqs. (3.46) and (3.48), the poles in η precisely cancel, and the dependence on e Y and υ cancel as well to yield a pure logarithm in ln(Q/q T ). This cancellation has to occur between the two collinear sectors, as there are no contributions from the soft sector. The NLP NLL result for the pure rapidity regulator is identical to that in eq. (3.44) upon dropping all rapidity-dependent pieces, which we have explicitly verified by repeating the derivation in sec. 3.1.3 using the pure rapidity regulator. This provides a highly nontrivial check of our regularization procedure, and our understanding of subleading-power rapidity divergences.

Next-to-leading Power Corrections at NLO
In this section, we give explicit results for the full NLP correction at NLO for gluon-fusion Higgs and Drell-Yan production in all partonic channels. Since both are s-channel processes, their power corrections are always proportional to their Born cross sections, and we express the NLP result as Here, we suppress the explicit Q and Y dependence in the kernels C (2,n) ab . The required H + j and Z + j amplitudes are conveniently expressed in terms of the Mandelstam variables which allows us to straightforwardly obtain the LP and NLP expansions in both the soft and collinear limits, as required by the collinear and soft master formulas. In the following, we only give the final results after combining soft, n-collinear, andn-collinear power corrections. The results were computed separately using both regulators, which provides a highly nontrivial check of our calculation.

Gluon-Fusion Higgs Production
We first consider on-shell Higgs production in gluon fusion in the m t → ∞ limit, for which the LO partonic cross section is given bŷ The LO matrix element in d = 4 − 2 dimensions is given by [129,130] At NLO, there are three distinct partonic channels, gg → Hg, qq → Hg, and gq → Hq, which we consider separately. Here, we calculate the full LL and NLL kernels for all channels. The LL results will be summarized in sec. 3.4.

gg → Hg
The spin-and color-averaged squared amplitude for g(p a ) + g(p b ) → H(q) + g(k) is given by [129] A gg→Hg (Q, The full result from combining the soft, n-collinear, andn-collinear contributions is given by Substituting these results into eq. (3.49) yields the NLP cross section for gg → Hg at NLO.

gq → Hq
The gq → Hq channel has power corrections at both LL and NLL. The spin-and coloraveraged squared amplitude for g(p a ) + q(p b ) → H(q) + q(k) is given by [129] A gq→Hq (Q, The full result from combining the soft, n-collinear, andn-collinear contributions is given by Substituting these results into eq. (3.49) yields the NLP cross section for gq → Hq at NLO.

qg → Hq
The result for qg → Hq can be obtained from eq. (3.56) by exchanging f q ↔ f g and a ↔ b, Substituting these results into eq. (3.49) yields the NLP cross section for qg → Hq at NLO.

qq → Hg
The qq → Hg channel has no leading logarithms and thus only contributes at NLL. The spinand color-averaged squared amplitude is given by [129] A qq→Hg (Q, The results for the kernels are given by Substituting these results into eq. (3.49) yields the NLP cross section for qq → Hg at NLO.

Drell-Yan Production
We next consider the Drell-Yan process pp → Z/γ * → + − , and for brevity denote it as pp → V . In contrast to on-shell Higgs production, it is important to be able to include off-shell effects. The LO partonic cross section is given bŷ where Q is the dilepton invariant mass, v ,q and a ,q are the standard vector and axial couplings of the leptons and quarks to the Z boson, and the + − phase space has already been integrated over. At NLO , there are two distinct partonic channels, qq → V g and qg → V q, which we consider separately. Here, we calculate the full LL and NLL kernels for all channels.
The LL results will be summarized in sec. 3.4.

qq → V g
We first consider the qq → V g channel, for which the spin-and color-averaged squared amplitude is given by [131] The full result from combining the soft, n-collinear, andn-collinear contributions is given by Substituting these results into eq. (3.49) yields the NLP cross section for qq → V g at NLO.
qg → V q The spin-and color-averaged squared amplitude for the qg → V q channel is given by [131] A qg→V q (Q, (3.63) The full result from combining the soft, n-collinear, andn-collinear contributions is given by Substituting these results into eq. (3.49) yields the NLP cross section for qg → V q at NLO.
The result for gq → V q can be obtained from eq. (3.63) by exchanging a ↔ b and f q ↔ f g , Substituting these results into eq. (3.49) yields the NLP cross section for gq → V q at NLO.

Discussion
Since the full calculation of the power corrections is rather involved, and contains a number of moving pieces, here we highlight several interesting features of the calculation, and compare them to the perturbative power corrections for beam thrust. For the purposes of this discussion, it is convenient to recall the form of the LL power corrections for the Born partonic configurations dσ (2),LL gg→Hg are identical up to switching of the labels on the PDFs. For the channels with a quark emission, we have dσ (2),LL gq→Hq are again identical up to the switching of the labels on the PDFs. First, we note that these results involve a more complicated structure of derivatives than the power corrections to the SCET I beam thrust observable, where at most a single derivative appeared in a given term [78,79,82]. Furthermore, for beam thrust, at LL there are no derivatives for the channels involving quark emission. Interestingly, the explanation for this arises from very different reasons in the soft and collinear sectors. In the soft sector, it is a simple consequence of the modified power counting of the soft modes, which implies that they must be expanded to two orders in the power counting. In the collinear sector, where the power counting is the same for q T and beam thrust, it arises from the presence of the power law singularities, which must be expanded against the PDFs. The cancellation of rapidity divergences between the soft and collinear sectors therefore exhibits a much more nontrivial relationship.
Another feature of the LL power corrections is the independence from explicit factors of the color-singlet rapidity Y , suggesting that the expansion parameter is indeed q 2 T /Q 2 , as is expected from the fact that q T is boost invariant. In fact, the rapidity dependence is induced purely by the PDFs and their derivatives. This is particularly interesting for the case of Drell-Yan, where the only terms that contribute arise from derivatives acting on the PDFs, which leads to a more nontrivial rapidity dependence, and in particular, a rapidity dependence that is different from that at leading power. This has potentially interesting implications for power corrections for q T subtractions, and we will show this rapidity dependence numerically in sec. 3.5. It is also interesting to discuss the universality of these results between Higgs and Drell-Yan production. For the case of beam thrust, the LL results are related by a Casimir scaling, C A ↔ C F . Here we see explicitly that this is not the case for q T . However, we see that all terms involving the derivatives of the PDFs are universal up to exchanges of the partonic indices. Here this is easily understood from the soft sector. In the soft sector, these terms arise solely from expanding the momenta entering the PDFs, which add power suppression to the LP matrix element. The LBK theorem guarantees that there is no term arising from an interference between two O(λ) suppressed terms. This then immediately implies the universality of this component of the result. The non-universality arises only in the coefficient of the f f term, which arises from corrections to the amplitude itself. It would be interesting to understand this in more detail, in particular how it extends to other processes, and to higher orders.

Numerical Results
In this section, we validate our results by numerically comparing the NLP spectrum to the full q T spectrum, which we obtain by numerically integrating eq. (3.18). For Drell-Yan production, we fix Q = m Z = 91.1876 GeV and use α s (m Z ) = 0.118. For Higgs production, we work in the on-shell limit with Q = m H = 125 GeV and α s (m H ) = 0.1126428 corresponding to a three-loop running from α s (m Z ). In both cases, we use E cm = 13 TeV and the NNPDF31 NNLO PDFs [132] with fixed factorization and renormalization scales µ f = µ r = Q. We also fix the rapidity to Y = 2 to have a nontrivial test of the rapidity dependence of our results and to break the degeneracy between the qg and gq channels.
We compare the nonsingular cross section at NLO 0 , 8 which is obtained by subtracting all singular terms which diverge as 1/q 2 T from the full q T spectrum, against our predictions for the NLP cross section. The dependence of the nonsingular cross section on q T is given by where c 1 is predicted by the LL term at NLP and c 0 is predicted by the NLL term at NLP. Note that c 0 is independent of q T , but has a nontrivial dependence on Q and Y . The O(q 2 T ) corrections arise at subsubleading power.
In fig. 1, we show the q T spectrum for all channels contributing to Higgs production. The corresponding results for Drell-Yan production are shown in fig. 2. In the left panel, we compare the nonsingular q T spectrum (solid red) against the NLP LL (green dashed) and full NLP (blue dashed) predictions. For all channels, the NLP NLL result is an excellent approximation of the nonsingular spectrum up to q T ∼ 10 GeV. The solid green line shows the nonsingular spectrum minus the NLP LL correction, which in all cases is almost perfectly constant up to q T ∼ 10 GeV, as expected from the structure of eq. (3.73). The solid blue line shows the nonsingular spectrum minus the full NLP correction, which vanishes as q 2 T for small q T as expected from eq. (3.73). This provides a strong numerical check of our analytic results of the NLP contributions. The right panels of figs. 1 and 2 compare the nonsingular spectrum q 2 T dσ/dq 2 T with the NLP LL and NLP NLL approximations. Again, we find excellent agreement up to q T ∼ 10 GeV.
In fig. 3, we show the rapidity dependence of the power corrections for the gg and qg channels for Higgs production and for the qq and qg channels for Drell-Yan production. We show the individual NLP terms as given in eq. (3.73), with the LL term proportional to c 1 shown in green and the NLL term proportional to c 0 shown in blue. Since their q T dependence is trivial, we fix q T = 1 GeV, which only affects the overall size of the LL term, and we normalize the results to the LO rapidity spectrum. Despite the fact that the kernels have no explicit rapidity dependence, we observe a nontrivial rapidity dependence due to the PDF derivatives, and in the case of the qg channels also because they involve different PDFs than the Born process. This is different than the case of beam thrust, which for certain definitions has an explicit rapidity dependence through factors of e ±Y in both the LL and the NLL kernels [78,82,90]. The rapidity dependence is particularly interesting for Drell-Yan production, where the term proportional to the PDFs themselves vanishes, see eq. (3.66), and so the power corrections are determined solely by the structure of the PDF derivatives.
At large values of |Y |, this leads to a relatively large dependence of the power corrections on the rapidity. For Higgs production this effect is more moderate due to the appearance of a term proportional to PDFs as present at LO, which dominates the rapidity dependence. This observation, which we believe is likely to persist at higher perturbative orders, could have important implications in the context of q T subtractions [95], where it is important to understand the rapidity dependence of the power corrections. Our results suggest that the rapidity dependence may be well behaved for the case of Higgs production but could be more problematic for Drell-Yan production. We leave the investigation of the structure at higher perturbative orders to future work. . Rapidity dependence of the LL (green) and NLL (blue) power corrections for Higgs and Drell-Yan production at NLO, relative to the LO rapidity dependence. The qq channel for Higgs production is not shown, as its LL power corrections vanish.

Conclusions
In this paper, we have studied in detail the structure and consistent regularization of rapidity divergences at subleading order in the power expansion. We have discussed several new features appearing at subleading power that put additional requirements on the rapidity regulator. As a result, most of the rapidity regulators that have been used in the literature at leading power become either unsuitable or inconvenient at subleading power. In particular, we have shown that the η regulator, which in principle can be applied at subleading power, is not homogeneous in the power expansion, which leads to undesirable complications at subleading power. We have introduced a new pure rapidity regulator, which is homogeneous in the power counting. It allows us to regulate rapidity divergences appearing in q T distributions at any order in the power expansion, while respecting the power counting of the EFT. This significantly simplified the analysis of rapidity divergences and the associated logarithms at subleading power. It would be interesting to study its application to other physical problems of interest and to further study its properties.
We have also found a rich structure of power-law divergences at subleading power, which can have a nontrivial effect on the final NLP result. Furthermore, at subleading power, rapidity divergences arise not only from gluons, but also from quarks. It would be interesting to further understand their formal properties.
As an explicit application of our formalism to a physical observable, we considered the q T spectrum for color-singlet production, for which we computed the complete NLP corrections, i.e., including both the logarithmic and nonlogarithmic contributions, at fixed O(α s ). This provides a highly nontrivial test of our regulator. In this case, the power-law rapidity divergences have the effect of inducing derivatives of the PDFs in the final NLP result for the q T spectrum. We also find that unlike for the case of beam thrust, where the LL power corrections for Higgs and Drell-Yan production are related by C A ↔ C F , this is not the case for the LL power corrections for q T , which have a different structure for these two processes.
Our results represent a first important step in systematically studying subleading power corrections for observables with rapidity divergences. It opens the door for addressing a number of interesting questions. It will be important to extend our results and to better understand the structure of subleading-power rapidity divergences at higher perturbative orders. As a particularly interesting application, the power corrections for the q T spectrum can be used to improve the numerical performance and to better understand the systematic uncertainties of q T subtractions, whose feasibility at next-to-next-to-next-to-leading order has recently been demonstrated in ref. [133] for Higgs production. We also hope that recent advances in the renormalization at subleading power, which has enabled the all-orders resummation of subleading-power logarithms, can also be extended to enable the resummation of subleading-power rapidity logarithms, with possible applications in a variety of contexts.

A NLO Results for q T at Leading Power
In this section we derive the LP beam and soft functions using the η regulator as a validation of our general setup.

A.1 Soft Function
The bare soft function at LP can be calculated using the known LP soft limit of a matrix element given in eq. (3.25), (A.1) We suppress that this limit only exists if either ab = gg or ab = qq. Inserting into eq. (3.26) and using eq. (2.18), we have Here, we also replaced the MS scale µ MS in terms of the MS scale µ using Choosing instead µ 2 = (4π) Γ(1− ) µ 2 MS would modify the O( 0 ) piece by π 2 /3. The divergence as q T → 0 is regulated using the distributional identity The terms in brackets yield the one-loop soft function integrated over the azimuthal angle of q T . The fully differential result can be read of as where the two-dimensional plus distributions are defined as in ref. [16], This result agrees exactly with the result in ref. [127]. + L 0 ( q T , µ) P gg (z) 2 + δ(1 − z) ln ω ν , (A.14) where ω = Qe Y . The finite part agrees with ref. [9], and thus after renormalization will give the same renormalized beam function kernel. Also note that the η poles cancel with the soft function eq. (A.6) after adding then-collinear beam function. The P gg (z)/ pole cancels with the UV divergence from the bare gluon PDF. The remaining pole and the 2 pole in the soft function, eq. (A.6), only cancel after taking virtual corrections into account.

B Higher-Order Plus Distributions
Subleading power corrections often involve divergences of the form In sec. 2.4 we encountered the two cases a = 2 and a = 3, which were treated using integration by parts to relate them to the case a = 1, where one can use the relation Here L n (x) = ln n x/x 1 + is defined in terms of standard plus distributions, which regulate functions g(x) with support x ≥ 0 diverging less than 1/x 2 as x → 0. The defining properties of such plus distributions are g(x) where g(x) has support x ≥ 0 and diverges less than 1/x 1+a as x → 0. For a = 1, this naturally reduces to eq. (B.3). For a = 2, one obtains the ++ distributions used e.g. in ref. [118].
The distributions defined in eq. (B.4) can be integrated against any test function f (x) that is at least a−1-times differentiable at x = 0. To be specific, consider the example integral where we assume x 0 > 0 and f (k) (0) is the k-th derivative of f (x) at x = 0. In eq. (B.5), we used that the term in square brackets in the first integral behaves as O(x a ) and thus cancels the divergent behavior of g(x) as x → 0, which allows us to drop the plus prescription in the first integral in the last line. In the second integral, we used eq. (B.4) to change the integration bounds from [0, x 0 ] to [x 0 , 1]. In the latter interval, g(x) is regular and the plus prescription can be dropped. The power-law divergence in eq. (B.1) can be regularized in terms of the higher-order plus distributions in eq. (B.4) as This result can be verified by integrating both sides against a test function (1 − z) m with m < a, and treating η as in dimensional regularization to render all integrals finite. In eq. (B.6), δ (k) (1 − z) is the k-th derivative on δ(1 − z), which thus induces a sign (−1) k in an integral over z and picks out the k-th derivative of any test function it acts on. Note that only the k = a − 1 term in eq. (B.6) diverges for η → 0, so irrespective of the power a, any power law divergence (1 − z) −a−η has exactly one single pole.
C Derivation of the Master Formula for Generic c In secs. 3.1 and 3.2, we derived master formulas for the NLP correction to the q T spectrum using the η regulator and the pure rapidity regulator, respectively. In sec. 2.3, we also introduced a class of homogeneous rapidity regulators spanned by a parameter c = 1. Here, we give the master formulas for this regulator for generic c = 1. In this regulator, the soft contribution is scaleless and vanishes, similar to the pure rapidity regulator. Thus, one only needs to consider the n-collinear andn-collinear limits. The derivation of the n-collinear expansion proceeds similar to the calculation shown in sec. 3.1.3. One can also obtain it from the result for the pure rapidity regulator, eq. (3.46), using the replacement This result is well-defined for all c = 1, whereas one encounters two explicit poles as c → 1. This behavior is expected because for c = 1 the regulator depends on the boost-invariant product k + k − = q 2 T and therefore does not regulate rapidity divergences, as explained at the end of sec. 2.3. For c = −1 we recover the result of pure rapidity regularization of eq. (3.46). In this case, the ν dependence in the regulator eq. (2.35) cancels, which is reflected by the vanishing of the coefficient of ln(ν/q T ) in eq. (C.2).
In then-collinear limit, the regulator for arbitrary c = 1 is obtained from the pure rapidity regulator through n (1) − 2A  Summing eqs. (C.2) and (C.4), the poles in η precisely cancel, and the dependence on c, υ and e Y cancels as well to yield a pure logarithm in ln(Q/q T ). As for the pure rapidity regulator, this cancellation has to occur between the two collinear sectors, since the soft sector does give a contribution. The NLP NLL result is identical to that in pure rapidity regularization, which is given by eq. (3.44) upon dropping all regulator-dependent pieces, as explained in sec. 3.2. This provides another check of our regularization procedure.