Drell-Yan $q_T$ Resummation of Fiducial Power Corrections at N$^3$LL

We consider Drell-Yan production $pp\to V^* X \to L X$ at small $q_T \ll Q$. Experimental measurements require fiducial cuts on the leptonic state $L$, which introduce enhanced, linear power corrections in $q_T/Q$. We show that they can be unambiguously predicted from factorization, and resummed to the same order as the leading-power contribution. We thus obtain predictions for the fiducial $q_T$ spectrum to N3LL and next-to-leading-power in $q_T/Q$. Matching to full NNLO ($\alpha_s^2$), we find that the linear power corrections are indeed the dominant ones, and the remaining fixed-order corrections become almost negligible below $q_T \lesssim 40$ GeV. We also discuss the implications for more complicated observables, and provide predictions for the fiducial $\phi^*$ spectrum at N3LL+NNLO. We find excellent agreement with ATLAS and CMS measurements of $q_T$ and $\phi^*$. We also consider the $p_T^\ell$ spectrum. We show that it develops leptonic power corrections in $q_T/(Q - 2p_T^\ell)$, which diverge near the Jacobian peak $p_T^\ell \sim Q/2$ and must be kept to all powers to obtain a meaningful result there. Doing so, we obtain for the first time an analytically resummed result for the $p_T^\ell$ spectrum around the Jacobian peak at N3LL+NNLO. Our method is based on performing a complete tensor decomposition for hadronic and leptonic tensors. In practice this is equivalent to often-used recoil prescriptions, for which our results now provide rigorous, formal justification. Our tensor decomposition yields nine Lorentz-scalar hadronic structure functions, which directly map onto the commonly used angular coefficients, but also holds for arbitrary leptonic final states. In particular, for suitably defined Born-projected leptons it still yields a LO-like angular decomposition even when including QED final-state radiation. We also discuss the application to $q_T$ subtractions.


Introduction
The neutral and charged Drell-Yan processes, pp → Z/γ * → and pp → W → ν, are important benchmark processes at the LHC. We are interested in the kinematic region where the vector boson is produced with small or moderate transverse momentum q T , which contains the bulk of the total cross section. In this region, differential distributions can be measured to sub-percent precision [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15], allowing for high-precision tests of the electroweak sector of the SM, including the precise measurement of the W boson mass [7] and the weak mixing angle [3,14].
For small transverse momentum q T Q, where Q is the invariant mass of the colorsinglet final state, the differential cross section admits an expansion in q T /Q dσ dq 2 (1.1) The dominant term scales as dσ (0) /dq 2 T ∼ 1/q 2 T and is referred to as the leading-power (LP) contribution. The additional terms dσ (n) are suppressed by (q T /Q) n relative to dσ (0) , and are referred to as power corrections or subleading-power contributions.
The power corrections dσ (n) in eq. (1.1) are classified by their relative (q T /Q) n suppression, and we refer to dσ (1) as the next-to-leading power (NLP) term, dσ (2) as NNLP etc. Due to their suppression, they are less relevant at small q T Q, and are included by matching to the full fixed-order calculations, which amounts to numerically extracting the complete set of power-suppressed terms at a given fixed order in α s . They are in principle known to O(α 3 s ) from the NNLO V + 1-parton calculations [93][94][95][96][97][98][99][100]. Nevertheless, the subleading-power terms also contain logarithms α n s ln m (q T /Q), and so in principle should be resummed as well to maintain their power suppression relative to the resummed LP term. Hence, given the high precision reached at LP, it is important to investigate the resummation of the subleading-power corrections to avoid them limiting the theoretical precision. First progress towards this direction has been made in ref. [101], where the power corrections were explicitly calculated at NLO, and in ref. [102], where the resummation at subleading power in a related, simpler context was studied. In ref. [101], it was explicitly shown that the linear NLP corrections for the inclusive q T spectrum are absent, i.e. dσ (1) = 0, consistent with earlier numerical observations, see e.g. ref. [103]. On the other hand, in ref. [104], it was shown explicitly that linear corrections do generically arise once fiducial cuts on the final-state leptons are applied. 1 In this work, we consider the generic Drell-Yan process pp → V X → LX, with the intermediate vector boson decaying to the "leptonic" (color-singlet) final state L. We study the origin and resummation of power corrections that arise from applying fiducial cuts or performing measurements on L, which we will refer to as fiducial power corrections. While our primary application will be to Z/γ * → and W → ν, most of our general analysis, which is carried out in section 2, will be for generic L. Our analysis and general results also immediately apply to the simpler case of an intermediate color-singlet scalar, such as Higgs production, though we will not consider this case explicitly here.
We encounter two classes of fiducial power corrections in our analysis: 1. Linear fiducial power corrections in q T /Q arise when azimuthal symmetry is preserved by the leptonic measurement at leading power, but is broken at O(q T /Q). For such measurements, the linear fiducial power corrections constitute the complete NLP corrections dσ (1) , and can be unambiguously predicted from factorization, and resummed to the same logarithmic order as the LP term dσ (0) .
The prototypical example is the q T spectrum in the presence of fiducial cuts on L, which generically break azimuthal symmetry and induce linear power corrections. It also applies to other more complicated q T -like observables, that resolve the recoil of the leptonic final state and vanish at Born level, e.g. the φ * observable or the scalar p T -imbalance p T 1 − p T 2 .
2. Leptonic fiducial power corrections in q T /p L arise when the leptonic measurement is sensitive to the edge of Born phase space, with p L corresponding to the distance to the Born edge. In the bulk of the leptonic phase space p L ∼ Q, and the discussion in point 1) applies. As p L gets smaller, the leptonic power corrections become enhanced, and for q T ∼ p L they become O(1) and must be retained exactly to all powers to obtain the actual LP result.
The prototypical example is the lepton p T spectrum close to the Jacobian peak p T = Q/2, with p L = Q − 2p T . Close to the Jacobian peak q T ∼ p L Q, fixed-order predictions are not reliable, which is a well-known effect. The resummation at strict LP is also not sufficient as it neglects the O(1) corrections in q T /p L . Hence, in this limit the resummation including all leptonic power corrections is required.
The inclusion of the fiducial power corrections in the q T factorization is derived in section 2. As we will see, the fiducial power corrections are a property of the leptonic decay and are independent of the underlying production of the decaying vector boson. This allows one to include them in the factorization theorem by treating the leptonic vector-boson decay exactly in q T and consequently makes it possible to resum them at the same level of precision as the singular cross section dσ (0) . In particular, this yields a resummation of the NLP terms dσ (1) to N 3 LL.
Our derivation in section 2 is general and independent of the specific method to perform the actual resummation, and of whether q T is treated as a perturbative scale or not. It is based on performing a Lorentz decomposition of the hadronic and leptonic tensors, which encode the production and decay of the intermediate vector boson. The basic idea of Lorentz-decomposing the hadronic tensor is of course not new and has been used before, typically to analyze the angular dependence for lepton pair production, see e.g. [99,[105][106][107][108]. Here, we use it for both hadronic and leptonic tensors to discuss the power counting at small q T . The tensor decomposition is discussed in section 2.3. It is constructed in a fully Lorentz-covariant way based on minimal requirements on symmetry and to make the small-q T limit maximally transparent, which leads to a direct equivalence with the Collins-Soper (CS) frame.
Our tensor decomposition holds for any leptonic final state L. In section 2.4, we show that for the specific cases of Z/γ * → and W → ν it directly maps onto the angular decomposition of the fully differential cross section in terms of CS angles. In section 2.4.4, we discuss that Born leptons have a well-defined theoretical interpretation as a Born projection of the full leptonic final state, and that in this case an analogous angular decomposition in terms of generalized angular coefficients also holds for generic L, in particular when including QED final-state radiation. This implies that the use of sodefined Born leptons is theoretically preferred over other lepton definitions in this context.
Our main power-counting analysis of both linear and leptonic fiducial corrections and their inclusion in the factorization is given in section 2.5. Some of the more technical details, such as the required power-counting of the hadronic tensor, are discussed in section 2.6 using soft-collinear effective theory (SCET) [109][110][111][112][113], which provides a systematic expansion of QCD in q T /Q ∼ λ. Our analysis does not rely on the precise formalism to factorize dσ (0) , and thus provides formal justification for existing approaches in the literature that include the exact lepton kinematics in the factorized cross section [65-67, 73, 74, 86, 87], as discussed in section 2.7.
In section 3, we summarize our specific q T resummation setup, implemented in the C++ library SCETlib [114], which we use to obtain numerical results for all factorized cross sections at fixed order and including resummation up to N 3 LL. Some additional details on the numerical inputs and computational setup can be found in section 4.1. In section 4, we discuss and illustrate in more detail the different sources of fiducial power corrections and the mechanism for their resummation. We consider three concrete examples, the q T spectrum with fiducial cuts (section 4.2), the p T distribution near the Jacobian peak (section 4.3), and the φ * distribution (section 4.4). In all cases, we validate numerically that the fiducial power corrections are indeed captured by the q T factorization, that their resummation significantly improves their perturbative stability, and that the size of remaining fixed-order power corrections is significantly reduced. In addition, we provide for the first time the resummed p T spectrum at N 3 LL+NNLO accuracy.
In section 5, we discuss the immediate implications of our findings for the q T subtraction method for fixed-order calculations [115], which is briefly reviewed in section 5.1. By including the fiducial power corrections predicted from q T factorization in the subtractions, their numerical performance improves tremendously. In the presence of fiducial cuts, it reduces the size of missing power corrections by an order of magnitude or more. Moreover, it makes the subtractions applicable also near the edges of Born phase space, where it otherwise breaks down due to uncontrolled leptonic power corrections. In section 5.2, we demonstrate this explicitly for the example of the p T spectrum in the vicinity of the Jacobian peak p T ∼ Q/2. In section 5.3, we discuss the example of Drell-Yan production with symmetric lepton cuts, for which large corrections due to a sensitivity to small q T have been observed before [116,117]. In fact, some of our numerical results in sections 4 and 6 rely on q T subtractions with fiducial power corrections to obtain stable results.
In section 6 we compare our resummed predictions at N 3 LL+NNLO for the fiducial q T and φ * distributions in pp → Z/γ * → + − with measurements by ATLAS [4] and CMS [15]. We compare the results both with the fiducial power corrections at fixed order as well as resummed, illustrating the improvement from resumming the fiducial power corrections and the fact that this significantly reduces the impact of the remaining fixedorder matching corrections.
To summarize, our general analysis is given in section 2, with the general setup and definitions in section 2.1, a review of the factorization for the inclusive q T spectrum in section 2.2, the hadronic tensor decomposition in section 2.3, the leptonic tensor and relation to angular coefficients in section 2.4, the main power-counting analysis in section 2.5, some of the more technical details in section 2.6, and the relation to other approaches in section 2.7. In section 3, we summarize our q T resummation setup. In section 4, we provide a detailed analysis of fiducial power corrections for the fiducial q T spectrum, the p T distribution, and the φ * observable. Section 5 discusses the applications to q T subtractions. In section 6, we compare our resummed N 3 LL+NNLO results for q T and φ * to ATLAS and CMS measurements. We conclude in section 7. In appendix A, we collect the hard functions and leptonic tensors required for Drell-Yan. Some additional data comparisons are provided in appendix B.

Factorizing production and decay
We consider the production of a single (generically off-shell) electroweak vector boson V in unpolarized proton-proton collisions and its subsequent decay into a set of colorless particles, which we refer to as the "leptonic" final state L. Throughout this paper we will work at leading order in the electroweak interactions. At this order, the matrix element for this process factorizes as where M µ V →L is the amplitude for V to propagate and decay into the leptonic final state L, X is any additional hadronic radiation in the final state, and J µ V is the electroweak qq current that couples to V , including electroweak charges and couplings. The polarizations of the hadronic current and of the propagating vector boson are encoded in the Lorentz index of J µ V and M µ V →L . The currents for V = γ and V = Z read 3) The differential cross section for pp → V X → LX in the lab frame, which we take to be the hadronic center-of-mass frame, is given by the square of eq. (2.1), integrated over phase space, and factorizes as (2.5) Here, q is the total momentum of the leptonic final state L (i.e., the momentum carried by V or V ). Parametrizing it in terms of the total leptonic invariant mass Q = q 2 and the rapidity Y and transverse momentum q T defined in the lab frame and with respect to the beam axis along the z direction, we have Importantly, in addition to q, the cross section depends on a set of differential observables O measured on L, as well as a set of fiducial cuts or angular projections Θ applied to L. The sum over V, V in eq. (2.5) runs over all intermediate vector bosons that contribute to the observed final state. In particular, it encodes the interference of V = γ with V = Z for neutral-current Drell-Yan. In the following, we suppress the dependence on V, V , as in the second line of eq. (2.5), unless it is of some relevance. The hadronic tensor W µν V V describes the QCD dynamics of the proton-proton collision. It is integrated over any additional hadronic radiation X and independent of the measurement and cuts performed on L, where the matrix elements of J µ V are implicitly averaged over proton spins. In addition to q, it depends on the incoming proton momenta P a,b . In the lab frame (and neglecting proton masses), P µ a = E cm 2 (1, 0, 0, 1) , P µ b = E cm 2 (1, 0, 0, −1) , (2. 8) where E 2 cm ≡ (P a + P b ) 2 is the hadronic center-of-mass energy. The leptonic tensor L µν V V describes the propagation and decay of the intermediate vector boson, (2. 9) In addition to q and the polarization of V encoded in the Lorentz indices, it depends on the measurement and cuts acting on the leptonic phase space point Φ L . The leptonic phase-space measure with total momentum q is defined as (2.10) It is straightforward to extend this setup to the collision of generic hadrons h a (P a ) and h b (P b ) with nonzero, possibly distinct masses m a and m b . This is relevant for treating proton or ion mass corrections in pp → XL, pA → XL, or AA → XL, where A and A are ions with these atomic numbers. In this case, the differential cross section in the lab frame becomes where the incoming hadron momenta in the lab frame are given by i.e., E a,b and v a,b are the lab-frame beam energies and beam velocities. Here, we also allow E a v a = E b v b , so the lab frame does not necessarily have to coincide with the hadronic center-of-mass frame. The leptonic tensor in eq. (2.11) is unchanged and given by eq. (2.9). The hadronic tensor is given by replacing the proton states |pp by hadron states |h a h b in eq. (2.7). For notational simplicity in the following we will always write the flux factor as 1/(2E 2 cm ) with the obvious replacement as in eq. (2.11) for the massive case understood. Lorentz invariance dictates that any Lorentz-scalar functions of q, P a , P b can only depend on Lorentz invariants formed out of q, P a , P b . There are six independent invariants, three of which contain nontrivial kinematic information, which we choose as (recall m T = (2. 13) In the second step we plugged in eqs. (2.6) and (2.12) to write s aq and s bq in terms of lab-frame quantities, and in the last step we took the limit of massless, center-of-mass collisions, which corresponds to taking v a,b → 1 and E a,b → E cm /2. It is clear that these are in one-to-one correspondence to the three kinematic variables Q, Y , and q 2 T = q 2 T . The three remaining invariants only encode the beam parameters, (2.14) The conservation of the vector current in QCD, ∂ µ J µ γ = 0, implies The same relation for V = Z does not automatically follow from gauge invariance, because the axial-vector current is not conserved in QCD due to finite quark masses and because of the Adler-Bell-Jackiw axial anomaly [118][119][120]. In the unbroken electroweak theory, the axial anomaly cancels in all gauge currents thanks to the anomaly cancellation in the SM. Since the anomaly coefficient is mass independent, it also does not contribute to the divergence of J µ Z after electroweak symmetry breaking, namely it still cancels between uptype and down-type quarks due to their opposite T u,d 3 = ±1/2. However, the nonzero quark masses now explicitly break the axial-vector current conservation. Therefore, we have the non-conservation relation 2 16) In practice, neglecting all but the top-quark masses, we thus have the chiral Ward identity At the partonic level, the leading contribution to this relation (without explicit top quarks in the final state) is the gluon-fusion top-quark triangle diagram. To isolate these nonconserved contributions, we write the hadronic matrix element as where the first term is "conserved" by construction, i.e., it vanishes when contracted with q µ , while the second term ∼ q µ contains the non-conserved pieces in eq. (2.17). Similarly, we can write the hadronic tensor as where the conserved part W µν cons arises from squaring the conserved parts of the currents. In practice, the non-conserved pieces rarely matter for various reasons: First, for a real, on-shell massive vector boson with physical polarization ε, they vanish due to q · ε = 0. As a result, for an off-shell vector boson near the resonance, they are suppressed by 1−q 2 /m 2 V . This is easy to see in unitary gauge, where all Goldstone bosons have been eaten up and the vector-boson propagator is proportional to g µν −q µ q ν /m 2 V . (In 't Hooft-Feynman gauge, the second term is generated by the exchange of Goldstone bosons.) Second, we can repeat the analogous discussion on the leptonic decay side, and split the leptonic tensor into conserved parts, q µ L µν cons = q ν L µν cons = 0, and non-conserved parts. The non-conserved parts of W µν are ∝ q µ,ν , and thus they only survive when contracted with the non-conserved parts of the leptonic tensor. However, considering leptonic decays (i.e., with the intermediate vector boson coupling to a leptonic current) the non-conserved leptonic parts are proportional to the lepton masses and can thus be neglected. 3 Therefore, for simplicity, we will ignore 2 As discussed, this relation is not anomalous. It also holds after suitable renormalization that preserves the non-renormalization of the axial anomaly [120], see e.g. refs. [121,122] for a detailed discussion.
3 A notable exception is associated Higgs production, which has a gg → Z * → ZH contribution. As a consequence of Yang's theorem, the ggZ vertex vanishes if all three bosons are real and on shell. Therefore, for real, on-shell gluons, the effective gg → Z contribution via a top-quark triangle is purely ∝ q µ and thus the gg → Z * → ZH process proceeds entirely via the non-conserved parts in eq. (2.17). Starting at O(α 2 s ), one or both gluons are off shell, and the ggZ vertex also contributes to the conserved parts, and therefore also to the Drell-Yan process Z → [17,123]. the non-conserved contributions for the most part, though we emphasize that they do not pose any additional conceptual problems and could be straightforwardly included in our analysis.

Factorization for the inclusive q T spectrum
If the measurement on L is inclusive, i.e., if we integrate over O and setΘ(q, Φ L ) = 1 in the leptonic tensor in eq. (2.9), it reduces to (2.20) In the second equality we used leptonic current conservation and the fact that after the integration L µν (q) can only depend on q. Explicit expressions for the scalar coefficients L(q 2 ) in case of Drell-Yan are given in appendix A.2. Inserting eq. (2.20) into eq. (2.5) yields where all QCD dynamics are encoded in the Lorentz-scalar inclusive hadronic structure function By Lorentz invariance, W incl can only depend on the three kinematic invariants q 2 , s aq , s bq , which are in one-to-one correspondence to the three kinematic variables Q, Y , q 2 T = q 2 T , see eq. (2.13). (In addition, W incl depends on the beam invariants in eq. (2.14), which we keep implicit.) In particular, since L(q 2 ) only depends on q 2 , the entire dependence on Y and q 2 T in eq. (2.21) is carried by W incl , and there is no dependence on the direction of q T . Hence, W incl encodes the inclusive (without fiducial cuts) q T distribution for fixed Q, Y .
We are interested in the region of small transverse momentum q T Q. In this limit, W incl satisfies the factorization theorem [56][57][58][59][60][61][62]64] (2.23) As indicated, W incl receives power corrections in (q T /Q) 2 and (Λ QCD /Q) 2 , but remains valid in the nonperturbative regime q T ∼ Λ QCD . In eq. (2.23), H ab (Q 2 , µ) denotes the hard function, which encodes the production of the vector boson in the underlying hard interaction ab → V . The MS scheme result for H ab can be obtained either by matching QCD onto SCET or as the IR-finite part of the corresponding form factor using dimensional regularization to regulate IR divergences. Explicit expressions for different vector bosons are given in appendix A.1. In practice, the leptonic prefactor L(q 2 ) is often included in the hard function in the inclusive case.
The second factor in eq. (2.23) encodes physics at the low scale µ ∼ q T , and can be written in several equivalent forms, In eq. (2.24a), the beam functions B i (x, k T , µ, ν/ω) describe the extraction of an unpolarized parton i with longitudinal momentum fraction x and transverse momentum k T from an unpolarized proton, the soft function S( k T , µ, ν) encodes soft radiation with total transverse momentum k T , and δ 2 ( q T − · · · ) encodes momentum conservation in the transverse plane. Eq. (2.24b) shows the equivalent result in Fourier space, whereB i andS are the Fourier conjugates of B i and S. Equivalently, one can write this as shown in eq. (2.24c), where the transverse-momentum dependent beam and soft functions have been combined into transverse-momentum dependent PDFs (TMDPDFs) A key feature of both transverse-momentum dependent beam functions and TMDPDFs is their explicit dependence on the energy of the colliding parton, encoded either in its lightcone component ω or in the Collins-Soper scale ζ, where where the lightcone components of the hadron momenta are given by and we also indicated the massless, center-of-mass limit. Accounting for the mass dependence of P − a and P + b implicit in the velocities v a,b ≤ 1 captures kinematic hadron-mass corrections to the factorization theorem in eq. (2.23). The factors of e ±Y in P ± a,b come from our lightcone conventions, see eq. (2.105), which imply that in the lab frame p ± = e ±Y (p 0 ∓p z ).
The ν/ω dependence of the beam functions or the ζ dependence of the TMDPDFs is a remnant of so-called rapidity divergences [56,[60][61][62][124][125][126]. Their regularization and renormalization induces an additional scale in the individual beam and soft functions in eq. (2.24a), here denoted as ν, analogously to the appearance of the MS scale µ from renormalizing UV divergences. Importantly, the ν dependence cancels between the beam and soft functions, such that eq. (2.24) is independent of ν. This fact allows one to combine beam and soft functions into ν-independent TMDPDFs as shown in eq. (2.25), where the Collins-Soper scale ζ is the remnant of the rapidity divergences.
In principle, the beam and soft functions (or TMDPDFs) are nonperturbative objects, and thus allow for a rigorous field-theoretic treatment of the q T spectrum in the nonperturbative regime q T ∼ Λ QCD . For perturbative | k T | Λ QCD , the beam functions (or TMDPDFs) can be matched perturbatively onto collinear PDFs [58,127], while the soft function is perturbatively calculable. The required perturbative results are known at N 3 LO [128][129][130][131][132][133][134][135][136][137]. In the perturbative regime, eqs. (2.23) and (2.24) allow one to resum large logarithms ln(q T /Q) arising to all orders in α s . In section 3 we review this procedure and describe the specific resummation setup used for the numerical results in this paper.
We note that there are various approaches in the literature on how to perform this resummation. While they all aim to describe the same inclusive hadronic tensor W incl and must ultimately all be based on the factorization theorem in eq. (2.23), they can differ in practice, e.g., due to differences in the rapidity regularization scheme, the different equivalent forms of eq. (2.24), different mathematical methods of performing the actual resummation, and different choices for the precise form of the logarithms that are being resummed. Crucially, all our results in this section 2 are general and hold independently of how precisely the resummation is carried out, and thus immediately apply to all formulations in the literature. 4 This is because they only rely on general arguments, such as Lorentz invariance and power counting, and the general structure of the hadronic and leptonic tensors.
The factorization theorem for the inclusive q T spectrum in eq. (2.23) receives corrections that are suppressed by powers of q T /Q relative to the leading term. As indicated in eq. (2.23), the leading corrections scale as (q T /Q) 2 , while linear power corrections are absent. This can be understood intuitively from the azimuthal symmetry of W incl , i.e., the fact that it only depends on the Lorentz invariants in eq. (2.13), which in turn only depend on q 2 T . The absence of linear power corrections in W incl has been verified explicitly by analytic O(α s ) calculations at next-to-leading power [101]. More formally, an argument for their absence to all orders in the inclusive case is presented in section 2.6. In the remainder of this section, we discuss how eq. (2.23) is extended to the case where the decay products are resolved and, notably, linear power corrections arise.

Hadronic tensor decomposition
We now return to the generic, fiducial cross section in eq. (2.11), and bring it into a form suitable for factorization at small q T . The manipulations of this section are exact in q T , i.e., we do not yet expand in q T Q. The key idea is to decompose the hadronic tensor W µν (q, P a , P b ) into Lorentz-scalar projections with respect to four orthogonal unit fourvectors that are constructed from the four-vectors P µ a,b and q µ and their invariants, and by imposing reasonable symmetry constraints.
For the decomposition to be complete, we should pick one timelike vector t µ and three spacelike vectors x µ , y µ , z µ , (2.28) Motivated by eq. (2.19), we take the timelike vector to be such that the conserved and non-conserved parts of W µν will get projected onto orthogonal components. The spacelike vectors must be given by linear combinations of P µ a,b and q µ . It will prove convenient to take z µ to lie in the plane spanned by P µ a and P µ b , where λ a and λ b are scalar functions of the kinematic invariants. Imposing t · z = 0 and z 2 = −1 then uniquely fixes z µ to up to a conventional overall sign. The s ij are all positive definite, as can be seen from their explicit expressions in eqs. (2.13) and (2.14), and s ab s aq s The choice for the remaining x µ and y µ is degenerate in principle. To reflect the fact that interchanging the initial-state hadrons is equivalent to a 180 • rotation about an axis in the transverse plane, we require x µ to be invariant under P a ↔ P b and y µ to only change sign. All together we then have We can write x µ as a linear combination of q µ and P µ a,b , where we chose the q µ coefficient to be positive to fix the overall sign of x µ . Imposing t · x = z · x = 0 and x 2 = −1, we find for the scalar coefficients and normalization factor (2.34) Finally, y µ is chosen to complete a righthanded coordinate system where we use the convention 0123 = +1. For completeness, the results for the unit vectors in the massless limit are (2.36)

Reference frame interpretation
The four-vectors t µ , x µ , y µ , z µ are orthogonal and normalized, and thus uniquely define a reference frame, namely the frame in which they have components t µ = (1, 0, 0, 0), x µ = (0, 1, 0, 0), y µ = (0, 0, 1, 0), and z µ = (0, 0, 0, 1). Since t µ = q µ / q 2 , this frame is automatically a frame where the vector boson is at rest, i.e., where q µ = ( q 2 , 0, 0, 0). A goal of this section is to show that this frame turns out to be the well-known Collins-Soper (CS) frame [138]. We will also find and discuss some subtleties in the massive case due to the fact that different CS-frame definitions that are equivalent in the massless case are no longer equivalent in the massive case. Let us first remind the reader that the vector-boson rest frame is not unique in itself because different rest frames can still differ by spatial rotations, i.e., by their orientation of the x, y, z-axes. There are many ways to perform a sequence of pure boosts to go from a given frame, say the lab frame, to the rest frame, and the difference between them precisely corresponds to an overall spatial rotation in the rest frame. Hence, a unique way to define a specific vector-boson rest frame is to specify the precise boost sequence to go from the lab frame to the rest frame. We will discuss how to rotate between different rest frames in section 2.6.3 below.
Intuitively, the CS frame is defined such that its z-axis points into the same direction as in the lab frame and its x-axis points in the direction of q T in the lab frame. In terms of boosts from the lab frame, the CS frame is defined by performing two boosts (see figure 3): 1. A longitudinal boost by Y in the beam direction (taken to be the z-axis) that takes us to the leptonic frame in whichỸ = 0 andq z = 0. Here and in the following, the tilde denotes the same physical quantity but evaluated in the leptonic frame.
2. A transverse boost in the direction of q T (taken to be the x-axis) with boost parameters which takes us from the leptonic frame to the rest frame.
Under these boosts a generic four-vector p µ transforms as where we explicitly indicated by a subscript in which frame the component-form is given, with p 0,x,y,z always denoting the lab-frame components andp 0,x,y,z always denoting the leptonic-frame components. To illustrate the boosts, applying them to q µ itself, we obtain Hence, we indeed arrive in the vector-boson rest frame, which is of course how eq. (2.37) was chosen in the first place.
We can now use this definition of the CS frame to make contact with our unit vectors t µ , x µ , y µ , z µ . To do so, we perform the same exercise for them, i.e., evaluate them in the lab frame and then boost them to the CS frame. For t µ = q µ /Q, this would just repeat eq. (2.39). For z µ , evaluating its general covariant expression in eq. (2.31) in the lab frame and applying the two boosts to the CS frame, we obtain Similarly, starting from the expression for x µ in eq. (2.33), we obtain This shows explicitly that the frame defined by t µ , x µ , y µ , z µ is equivalent to the CS frame (in its boost definition), and that this equivalence also holds in the general massive case. It is quite pleasing to see that the CS frame naturally appears in a covariant way only by imposing eq. (2.30) and the symmetry constraints in eq. (2.32). Another definition of the CS frame [138], which is also often used in practice, is to consider P a and P b in the vector-boson rest frame, and to define the z-axis to bisect the angle between P a and − P b , while the x-axis is chosen to lie in the plane defined by P a and P b . Denoting individual components in the CS frame (defined via the above boosts) by hats, we have P µ a = E a (1, 0, 0, +v a ) lab = (P 0 a ,P x a , 0,P z a ) CS , where explicit expressions for the components can be straightforwardly obtained from eq. (2.38). The angles γ a,b between P a,b and the z-axis (see figure 1 right) are given by (2.43) The bisector criterion amounts to requiring these two angles to be equal, i.e., This can only be satisfied for generic Y if and only if v a = v b = 1, i.e., both hadrons are massless. This means the bisector definition of the CS frame is only equivalent to the above boost definition for massless hadrons, for which both definitions where originally introduced in ref. [138], while for nonzero hadron masses the two definitions are no longer equivalent. 5 The key advantage of our construction of t µ , x µ , y µ , z µ and the corresponding boost definition of the CS frame is that they are symmetric under interchanging the beams (see eq. (2.32)) and furthermore are manifestly independent of the beam parameters, i.e., they only depend on q µ without reference to the beam momenta beyond the beam direction itself. In the rest of the paper, we will always use this definition, unless stated otherwise.

Helicity decomposition
Using t µ , x µ , y µ , z µ , we can define polarization vectors for the vector boson in a fully covariant way as which correspond to positive/negative helicity and longitudinal polarization with respect to z µ . Using these, we project the hadronic tensor onto the entries of a helicity density matrix [106], Since the ε µ ±,0 span the space orthogonal to t µ = q µ /Q, this decomposition fully captures the conserved part of the hadronic tensor, see eq. (2. 19). (To also account for the nonconserved parts, we would just have to include the fourth time-like polarization t µ .) From its definition in eq. (2.7), it is clear that W µν is hermitian, W * µν = W νµ , so its symmetric (antisymmetric) components are purely real (imaginary). Therefore, the nine helicity components W λλ are fully specified by a total of nine real-valued, Lorentz-scalar hadronic structure functions. We will use the following linear combinations: The reason for the somewhat odd numbering and normalization will become apparent shortly. In the second equality, we have given the projections in terms of x µ , y µ , z µ , corresponding to linear vector-boson polarizations. The inclusive structure function from eq. (2.22) is given by Since the projections of W µν that define the W i are orthogonal, we can easily invert them and write W µν in terms of the W i as where the P µν i are the same projections as in eq. (2.47) up to a trivial difference in normalization, for example, Contracting the leptonic tensor L µν with W µν decomposed as in eq. (2.49), we have with the corresponding leptonic structure functions defined as (2.52) The cross section in eq. (2.11) in terms of the scalar structure functions now becomes (2. 53) which generalizes the inclusive cross section in eq. (2.21) to arbitrary leptonic observables and fiducial cuts. As for W incl before, Lorentz invariance implies that the hadronic structure functions W i only depend on the three kinematic invariants q 2 , s aq , s bq , or equivalently the three kinematic variables Q 2 , Y , q 2 T , see eq. (2.13). In particular, they do not depend on the orientation of q T . Since the x µ , y µ , z µ reduce to the spatial coordinate axes in the CS frame, the structure functions correspond to the individual tensor components of the hadronic tensorŴ µν evaluated in the CS frame, e.g., W −1 =Ŵ xx +Ŵ yy , W 0 = 2W zz , etc. We will refer to eqs. (2.47) and (2.49) as the CS tensor decomposition.
We note that one may also decompose the hadronic tensor in terms of Lorentz structures directly formed out of g µν − q µ q ν /q 2 and its contractions with P µ a,b , see e.g. refs. [99,105,107,108]. This automatically ensures that the projectors are covariant combinations of q µ and P µ a,b and that the corresponding coefficients are Lorentz-scalar functions. This is usually not manifest when one considers the individual tensor components in the CS frame (or any other rest frame). However, as we have seen, the CS-frame components are reproduced by the CS tensor decomposition in a manifestly covariant manner as the Lorentz-scalar structure functions W i that only depend on Lorentz invariants. Hence, there is no formal preference for either decomposition and the two are related by a straightforward change of basis. We will see in the following sections that the physics at small q T Q becomes particularly transparent when using the CS tensor decomposition.

Leptonic decomposition and relation to angular coefficients
In this subsection, we discuss the leptonic decay in more detail. For the most part, we specifically consider the leading-order Drell-Yan decays neglecting lepton masses, m 1,2 = 0, and summing over lepton polarizations. These are the primary application we are eventually interested in. The kinematics of the process in the lab and CS frames are illustrated in the left and right panels of figure 1. In section 2.4.4, we discuss the extension to more complicated leptonic final states, e.g. including QED final-state radiation, which is important at the precision of current Drell-Yan measurements. In particular, there we show to what extent the LO discussion carries over for measurements that are performed in terms of suitably defined Born leptons.

Definition of CS angles
It is convenient to introduce spherical coordinates (cos θ, ϕ) in the CS frame, in terms of which we can parametrize p 1,2 , as illustrated in the right panel of figure 1, as p µ 1,2 = Q 2 t µ ± x µ sin θ cos ϕ ± y µ sin θ sin ϕ ± z µ cos θ . In the lab frame, the incoming hadron momenta are head-to-head (assuming the lab frame and hadronic center-ofmass frame coincide), while the vector boson has nonvanishing three-momentum q. The scattering p(P a ) p(P b ) → V (q) X(p X ) defines the hadron plane (green). In the CS frame (right), the vector boson is at rest. The leptons are produced back to back in the lepton plane (blue). The magnitudes of the hadron momenta in general differ for Y = 0, but their angles γ a,b with respect to the z axis (indicated by the double arcs) become equal in the limit of vanishing hadron masses. The Collins-Soper angles θ and ϕ are defined as indicated.
The angles θ, ϕ are known as Collins-Soper angles. 6 From eq. (2.55), one can easily derive their explicit expressions in terms of lab-frame quantities E 1,2 , p x,y,z 1,2 , Note that we have arbitrarily chosen the positive orientation of the z axis by having hadron a move in the z direction in the lab frame. As a result, the negatively charged lepton moves into the same rest-frame hemisphere as hadron a for cos θ > 0. In experimental measurements at the LHC, where the choice of a and b is arbitrary, hadron b is often taken to be the one closer to the vector boson in rapidity to ensure that angular distributions do not average out when integrating over rapidity, see e.g. refs. [3,5,11,14]. The resulting angles θ * and ϕ * , which are often also referred to as Collins-Soper angles, are then related to eq. (2.56) by On the other hand, eq. (2.56) does not depend on the chosen orientations of the x and y axes in the lab frame as long as they form a right-handed coordinate system. The advantage of eq. (2.55), or equivalently the boost definition to define the CS frame, is that it stays true regardless of whether hadron masses are included or neglected, and thus also any relations like eq. (2.56) that are derived from it are independent of any beam parameters. On the other hand, with the bisector construction including hadron masses, eq. (2.56) no longer holds, see also the discussion in section 2.3.1.

Leptonic decay parametrization by angles
The fully-differential leptonic tensor for the 1 → 2 Drell-Yan decays in eq. (2.54) at tree level has the form Only the contribution proportional to L + (L − ) survives the contraction with the symmetric (antisymmetric) P µν i corresponding to the parity-even (parity-odd) hadronic structure functions W −1,0,1,2,5,6 (W 3,4,7 ). The normalization is chosen such that L + (q 2 ) = L(q 2 ) agrees with the inclusive coefficient in eq. (2. 20), and such that L − (q 2 ) = L + (q 2 ) for W decays, where parity is maximally violated. Explicit expressions for the L ± (q 2 ) are given in appendix A.2.

Relation to angular coefficients
From eq. (2.62), we can write the fully-differential cross section in the CS angles as where the angular coefficients A i are given in terms of the helicity cross sections or the hadronic structure functions as . (2.66) We deliberately chose the numbering and normalization of the W i in eq. (2.47) to match the often used form of the cross section in eq. (2.65). The only exception is the inclusive cross section, which is split into orthogonal contributions from W 0 and W −1 . For the same reason, we refrained from normalizing the spherical harmonics in eq. (2.61). We remind the reader that both numerator and denominator in eq. (2.66) in general involve a sum over the intermediate vector bosons, , the parity-even leptonic prefactors L + (q 2 ) ≡ L +V V (q 2 ) do not in general cancel in the ratio in eq. (2.66). A priori, eq. (2.62) or eq. (2.65) simply provide a convenient way to parametrize the fully-differential Drell-Yan cross section for massless 2-body decays. For this purpose, it is irrelevant whether or not the CS angles can be reconstructed experimentally. Similarly, the choice of the CS tensor decomposition is a priori arbitrary, and we could have used another decomposition. Of course, the combination of using the CS tensor decomposition for the hadronic tensor together with using the CS angles to parametrize the leptonic tensor is what leads to the simple angular dependence in eq. (2.62). If we were to choose a different tensor decomposition W i , we would also choose polar coordinates cos θ , ϕ with respect to its corresponding rest frame, and arrive at eqs. (2.65) and (2.66) in terms of cos θ , ϕ , A i , and W i . On the other hand, when cos θ and ϕ are explicitly measured, or when eq. (2.65) is used as a template to measure the A i , it obviously does matter with respect to which frame they are defined. It is also straightforward to relate the W i or A i for different frames, see section 2.6 below.

Extension to more complicated leptonic final states
Up to now, our discussion in this subsection assumed the leading-order dilepton final states in eq. (2.54), and so in particular eq. (2.65) is derived in this limit. For a generic leptonic final state L, e.g. when including QED final-state radiation (FSR) or for more complicated electroweak decays like V * → V H or V * → V 1 V 2 , there is a priori no reason that the L i are proportional to spherical harmonics g i (θ, ϕ) any longer, in which case one cannot use eq. (2.65) to define the A i beyond this LO.
On the other hand, as we saw in eq. (2.66), the A i are in one-to-one correspondence with the underlying hadronic structure functions W i . The W i are by construction independent of L (apart from its total momentum q µ ) and thus well-defined for arbitrary L. The physical reason for the appearance of nine independent structures in both cases is exactly the same, namely the spin-1 nature of the intermediate vector boson (and the fact that we ignore the non-conserved parts). Hence, the cross section in the CS tensor decomposition in eq. (2.53) should be considered as the generalization of the LO angular decomposition in eq. (2.65) to arbitrary leptonic final states and measurements. One could also use eq. (2.66) as the all-order definition of the A i in terms of the W i and conventional LO weak couplings and propagators included in the L ± (q 2 ). One could then easily rewrite eq. (2.53) in terms of the so-defined A i multiplied by generic leptonic coefficients L i (q, O, Θ), which in the simplest case reduce to L ± (q 2 )g i (θ, ϕ) as in eq. (2.60), but in general can also be more complicated. Although at that point, it is easier and perhaps less confusing to directly work in terms of the W i and eq. (2.53) as it is.
Nevertheless, in the context of Drell-Yan measurements, the LO relation in eq. (2.65) is very useful in practice because the g i are orthogonal spherical harmonics. This allows one to directly measure the A i (or W i ) by performing a fit to the angular dependence of the (θ, ϕ) distribution or by projecting out different terms by taking suitably weighted angular integrals of it [5,11]. This procedure has received some criticism, since it seemingly relies on a LO QED interpretation of the angular dependence, while QED final-state radiation can be relevant at the level of precision reached by Drell-Yan measurements. In fact, even the definition of the CS angles (θ, ϕ) themselves becomes nonobvious, because with additional QED radiation in the final state, the lepton momenta generically no longer add to the full vector-boson momentum q µ . Instead, we now have where p µ 1,2 are the measured lepton momenta, which depend on the lepton definition, and k µ is the remaining momentum not included in the definitions of p 1,2 . We stress that here we are not concerned with the experimental methods to reconstruct and calibrate the leptons or to recover photon radiation. The "measured" lepton momenta p µ 1,2 refer to the truth-level lepton definition to which the raw reconstructed momenta are corrected or unfolded. This truth-level definition must be theoretically well-defined to have a meaningful measurement that can be compared to theoretical calculations, and one can consider the question whether certain truth-level definitions are theoretically preferred or not. 7 Obviously, the (θ 1 , ϕ 1 ) angles describing the orientation of p 1 now depend on the lepton definition and also on whether they are defined in the full vector-boson rest frame (where q µ or equivalently the full L is at rest) or the dilepton rest frame (where only p µ 1 + p µ 2 is at rest). Especially in the latter case, there is no guarantee (in fact it seems quite unlikely) 7 On the other hand, whether a specific truth-level definition receives more or less associated experimental uncertainties is a separate, experimental question, to which we have nothing to say here. While these two questions are not entirely unrelated, they should nevertheless be kept well separated. We thank Daniel Froidevaux for discussions on this issue. that the angular distribution in (θ 1 , ϕ 1 ) will still admit a decomposition in terms of the nine spherical harmonics g i (θ 1 , ϕ 1 ).
For "bare" leptons, p µ 1,2 are defined without including any FSR photons. This means infrared QED singularities are regulated by the lepton mass leading to potentially large logarithms of the lepton mass. This effect is reduced by defining "dressed" leptons, which include all photons radiated within a cone of some size around the leptons, and hence can be theoretically thought of as QED "lepton jets". With either definition, the remaining momentum k µ in eq. (2.67) is nonzero and so the dilepton and vector-boson rest frames are no longer equivalent.
Another option is to include all k µ into p µ 1,2 , i.e., the lepton momenta are (partially) defined by the condition q µ = p µ 1 + p µ 2 . This is basically what "Born" leptons are. Their full definition corresponds to defining an IR-safe projection of the full leptonic final state L onto a Born-like 2-body final state. In principle there are many ways to do so, but as long as the Born projection is well defined so are the Born leptons.
To illustrate this, let us consider an explicit example: We start by defining the leptonic thrust axis n L of the full leptonic final state L in its rest frame. The thrust axis n L is defined in the usual way as the axis n, with n 2 = 1, that minimizes where the sum runs over all particles in L, including in particular all QED FSR, and E i , p i are defined in the rest frame of L. The overall positive (negative) orientation of n L can be fixed by convention, e.g., to point into the hemisphere that contains the lepton (antilepton). 8 Imposing the condition p µ 1 + p µ 2 = q µ , requiring massless on-shell momenta, p 2 1,2 = 0, and using n L to define the direction of p 1 = − p 2 , then uniquely determine (recall In the second step we defined the unit vectors where t µ is the same as before, and n µ L describes the overall orientation of L. More generally, we can also carry out the construction in two steps, first constructing q µ = P µ 1 + P µ 2 with massive P 2 1,2 = 0 and then projecting them onto massless p 1,2 . Here, we first cluster all emissions with either the lepton or antilepton based on whose hemisphere they are in, which yields the massive hemisphere momenta P µ 1,2 , (2.71) Next, we project P µ 1,2 onto massless momenta p µ 1,2 by preserving the three-momentum direction, p 1 /| p 1 | = P 1 /| P 1 | = n L , and the total energy, p 0 1 + p 0 2 = P 0 1 + P 0 2 = Q, which yields eq. (2.69). The spherical coordinates (θ L , ϕ L ) of n L in the CS frame now provides a unique, all-order definition of the CS angles (θ, ϕ) ≡ (θ L , ϕ L ), i.e., This is the generalization of eq. (2.55), where the CS angles now describe the overall orientation of L in the CS frame.
It is easy to see that the above definitions are IR safe and reduce to the respective LO definitions. In principle, any other IR-safe way of clustering the emissions into P µ 1,2 is possible. Other ways to project them onto massless p µ 1,2 are also possible, as long as the projection is IR safe and preserves the total leptonic momentum q µ = P µ 1 + P µ 2 = p µ 1 + p µ 2 . In practice, defining the projection by keeping the orientation fixed is the most natural and also the easiest, as it avoids any confusion about which particular direction is used to define the CS angles.
The advantage of Born leptons is that they do admit an analogous LO-like angular decomposition as we will now show. More generally, it is sufficient to restrict to leptonic measurements that can be written in terms of P µ 1,2 , For such measurements we can write the general leptonic tensor in eq. (2.9) as where F µν (P 1 , P 2 ) is the projection of the full leptonic decay L µν (Φ L ) onto the massive 2-body (P 1 , P 2 ) phase space, with L µν (p 1 , p 2 ) given by the LO result in eq. (2.58).
The key point is that F µν (P 1 , P 2 ) is defined in a Lorentz-covariant way, and therefore obeys the following Lorentz decomposition (ignoring as before the non-conserved parts) where F ±,0 ≡ F ±,0 (q 2 , P 2 1 , P 2 2 ) are Lorentz-scalar functions that can only depend on three independent invariants formed out of P 1,2 , which we chose as q 2 = (P 1 + P 2 ) 2 and P 2 1,2 . The decomposition in eq. (2.77) is chosen so that at LO, comparing to eq. (2.58), we have The leptonic structure functions are now obtained as defined in eq. (2.52), by contracting eq. (2.77) with the projectors P µν i and performing the phase-space integrals in eq. (2.74), with the underlying leptonic structure functions given by L 1,2,5,6 (q 2 , θ, ϕ, P 2 1 , P 2 2 ) = 3 16π F + (q 2 , P 2 1 , P 2 2 ) g 1,2,5,6 (θ, ϕ) , Eqs. (2.79) and (2.80) are the generalization of eq. (2.60) to an arbitrary Born-projected leptonic final state L. The (θ, ϕ) dependence is still completely described by the same g i (θ, ϕ) in eq. (2.61). The L i for i ≥ 1 are still given by their own respective g i times a common leptonic form factor F + for i = 1, 2, 5, 6 and F − for i = 3, 4, 7. On the other hand, the angular dependence of L −1 and L 0 now gets mixed up by F 0 , which enters with a flat (θ, ϕ) dependence corresponding to g −1 + g 0 = 2. If the measurements are defined in terms of massless Born leptons, then they are also independent of P 2 1,2 , such that the P 2 1,2 integrals in eq. (2.79) can be performed to give L i (q 2 , θ, ϕ) that are given by the same expressions as in eq. (2.80) but in terms of corresponding integrated F ±,0 (q 2 ) = dP 2 1 dP 2 2 λ(q 2 , P 2 1 , P 2 2 ) F ±,0 (q 2 , P 2 1 , P 2 2 ) . (2.81) Removing all leptonic measurements, the inclusive q T spectrum is now given by i.e., in terms of the same inclusive hadronic structure function multiplied by a generalized inclusive leptonic function. We remind the reader that also here there is always an implicit sum over intermediate vector bosons, The cross section differential in the CS angles becomes where we suppressed the arguments of the structure functions for brevity, and the analogous expression also holds differential in P 2 1,2 . To make contact with eq. (2.65), in the second step we split the flat contribution from F 0 as (3/2)(1 + cos 2 θ) + (1/2)(1 − 3 cos 2 θ) = 2, and in the last step we factored out the inclusive cross section, denoting the resulting normalized coefficients of the angular dependence asÃ i , These are the generalization of the A i in eq. (2.66) for an arbitrary Born-projected final state. They implicitly depend on the specific Born projection used because the CS angles (θ, ϕ) implicitly depend on it. TheÃ i are the angular coefficients that are measured by decomposing or projecting the (θ, ϕ) dependence defined in terms of Born-projected leptons. It would be interesting to precisely identify the underlying Born projection that is effectively used in the measurements [5,11]. Generically, the QED corrections to F + , F − , and F 0 will differ and thus not cancel in eq. (2.84). In other words, even though Born-projected leptons admit a well-defined LO-like angular decomposition as shown in eq. (2.83), the resultingÃ i in eq. (2.84) still differ by QED FSR corrections from the LO A i in eq. (2.66). These corrections should be of generic O(α em ) size, i.e., neither enhanced by soft or collinear photon emissions nor suppressed near the Z pole. In the limit of an on-shell Z boson, they would produce the QED corrections to the Z decay rate to leptons. For i ≥ 1, the hadronic contributions to A i andÃ i are the same. As we will discuss below, W 0 is suppressed by O(q 2 T /Q 2 ) relative to W incl at small q T , such that at LO in QED A 0 vanishes like q 2 T for q T → 0. Interestingly,Ã 0 receives an additional contribution F 0 W incl , and therefore it no longer vanishes for q T → 0 but goes to a calculable O(α em ) constant.
Ref. [139] considered QED radiation off massive final-state leptons, and found linear power corrections even in the inclusive case. Since their massive leptons correspond to bare leptons, this is not entirely surprising. It would be interesting to identify the precise source of linear power corrections, i.e, whether the bare leptons induce linear corrections in the leptonic tensor itself, or populate additional leptonic structure functions that come with linearly suppressed hadronic structure functions, or both.
Finally, while most of the above discussion was phrased in terms of QED FSR corrections to Drell-Yan, it applies to an arbitrary Born-projected final state L. For example, keeping the P 2 1,2 dependence, it applies to Drell-Yan-like electroweak diboson production V * → V H or V * → V 1 V 2 if one remains inclusive over the decays of the final-state bosons.

Factorization for fiducial power corrections
We now investigate the structure of power corrections in the limit q T Q in the presence of measurements on the leptonic final state. To expand in q T Q, we introduce a formal power-counting parameter The leptonic measurementsÔ andΘ in eq. (2.52) are functions of the total four-momentum q of the final state, and admit an expansion in λ aŝ We refer to the corrections in λ in these expansions as fiducial power corrections. For observables that exist at Born level, e.g., cuts on the lepton momenta, the leading-power (LP) observablesÔ (0) andΘ (0) are simply obtained by taking the Born limit q T → 0. For q T -like resolution variables like φ * that scale like q T itself and vanish at Born level,Ô (0) or Θ (0) are given by the leading, nontrivial contribution in the q T → 0 limit.

Linear fiducial power corrections
We first assume that the leptonic measurement does not induce any additional nontrivial dynamic scale p L , such that the power expansion in eq. (2.86) is genuinely an expansion in q T /Q. We can then focus on the linear O(λ) fiducial power corrections. Let us consider leptonic measurements that are azimuthally symmetric at leading power, which we will indicate by L (0) (/ ϕ) and define more precisely in a moment. We will show that for such measurements the only linear O(λ) power corrections that arise are due to the linear fiducial power corrections in eq. (2.86). As a result, the O(λ) power corrections can be uniquely predicted and resummed in terms of leading-power hadronic structure functions.
For measurements that can be parameterized in terms of CS angles θ, ϕ, which includes our default Drell-Yan cases in eq. (2.54), azimuthal symmetry means that they do not depend on ϕ. Azimuthal symmetry at leading power then simply means that the LP measurementsÔ (0) (q, θ) andΘ (0) (q, θ) are ϕ independent, which implies that they average out against cos(nϕ) and sin(nϕ), In particular, the integration against all spherical harmonics g i (θ, ϕ) in eqs. (2.60) and (2.61) vanishes, except for i = −1, 0, 4, which do not depend on ϕ. More generally, we define a generic leptonic measurement as azimuthally symmetric if it only contributes to L −1,0,4 , such that azimuthal symmetry at leading power is defined by Note that this definition is also natural from the point of view of the CS tensor decomposition in eq. (2.47). Azimuthal symmetry corresponds to symmetry under rotations of the x and y axes around the z axis. The projections for i = −1, 0, 4 are precisely those that are invariant under azimuthal rotations (corresponding to the norm and cross product, or are independent of x and y), which is the physical reason why their corresponding g i (θ, ϕ) do not depend on ϕ.
A primary example is a fiducial cut on the lepton transverse momenta p T 1,2 ≥ p min T , for which we havê In words, the leptons are exactly back-to-back at leading power, and whether they pass the cut only depends on their rest-frame energy Q/2 and scattering angle θ. We discuss this case as well as a cut on the lepton rapidity in more detail in section 4.2. On the other hand, angular asymmetries that are designed to project out the cos ϕ or cos(2ϕ) dependence in the angular distribution (by construction) do not qualify under eq. (2.88). Power expanding the leptonic structure functions, which includes the power expansion of the measurement, we have (2. 90) where with the assumption in eq. (2.88) only L i . We also need to power-count the hadronic structure functions W i , (2.91) The λ scaling of the first nonzero contributions W , and is derived more carefully using SCET in section 2.6. From table 1, we see that the only nonvanishing LP structure functions W (0) i are for i = −1, 2, 4, 5. The physical reason is that at LP, angular momentum conservation works the same way as at tree level, i.e., as in the collision of two massless partons with p a +p b = q, p T a = p T b = q T = 0. In this limit, the CS frame coincides with the leptonic frame, and the longitudinal polarization vector is given by Table 1. Scaling of the hadronic structure functions W i in the CS tensor decomposition in eq. (2.47) in the limit λ ∼ q T /Q 1 relative to the leading W −1,4 ∼ 1/q 2 T . In this table we generically count Λ QCD ∼ q T . In several cases we only derive bounds on the scaling, where ∼ λ ≥m means the W i is suppressed by at least λ m . We group the structure functions by parity and whether they arise from the dispersive (i = −1 . . . 4) or absorptive parts (i = 5 . . . 7) of the production amplitude [140]. The second-to-last column shows the corresponding angular dependence g i (θ, ϕ) on the Collins-Soper angles for 2-body decays, and in the last column we indicate whether there is a nonvanishing LP leptonic tensor L (0) i for observables that are azimuthally symmetric at Born level.
It is easy to see that projections of the tree-level partonic matrix element onto ε µ 0 vanish, for any polarization s a (s b ) of the quark (antiquark), and similarly for the axial-vector current. It follows that structure functions W i that involve contractions with ε µ 0 vanish at tree level. We will see in section 2.6 that to all orders, each contraction with ε µ 0 is in fact penalized by at least one power of λ.
Suppressing the arguments of L i and W i , the strict LP cross section is given by (2.94) The i = 2, 5 contribution does not survive because L (2. 95) In the first term, only i = −1, 2, 4, 5 contribute to the sum due to table 1. For the second term, assuming LP azimuthal symmetry, only i = −1, 0, 4 contribute. From table 1, indicates that this contribution does not match onto leading-twist collinear PDFs and is suppressed for Λ QCD q T . Gray boxes indicate contributions to the LP and linear NLP cross sections σ (0) and σ (1) . The latter solely arise through the leptonic tensor as L vanish. Crosses indicate that a contribution vanishes. In the bottom right panel, "rest" refers to structure functions i = 1, 3, 6, 7.
W 0 ∼ O(λ 2 ), and as we will argue in section 2.6, all power corrections to W −1,4 are quadratic in λ, so we have For W −1 and W 0 , this statement is equivalent to the absence of linear power corrections in the inclusive cross section ∝ W −1 + 1 2 W 0 . For W 4 , it is equivalent to the absence of linear power corrections in the inclusive forward-backward asymmetry. Hence, the second term in eq. (2.95) vanishes, and we arrive at (2.97) We have thus shown that for leptonic measurements that are azimuthally symmetric at leading power, all linear power corrections uniquely arise from linear fiducial power cor-rections L (1) i multiplying the leading-power hadronic structure functions W (0) −1,2,4,5 . The power-counting logic leading to eq. (2.97) is illustrated in figure 2.

Leptonic fiducial power corrections
We now turn to leptonic fiducial power corrections that arise from the presence of an additional, physical scale p L induced by the leptonic measurement. In this case, the power expansion of the measurements a priori receives power corrections in both q T /Q and q T /p L , The case of linear fiducial power corrections discussed above corresponds to p L ∼ Q. We refer to the q T /p L corrections as leptonic fiducial power corrections. For q T p L Q, they become enhanced compared to the q T /Q corrections and for q T ∼ p L they become O(1) and cause the naive expansion in q T to break down.
Generically this happens when the leptonic measurement is close to an edge of Born phase space that is sensitive to additional radiation, such that a nonzero q T opens up new phase space beyond the Born edge, with p L ∼ q T parametrizing the distance from the Born edge. We will demonstrate this effect in detail in section 4.3 for the important example of the p T spectrum near the Jacobian To expand in such regions, it is necessary to count both which explicitly avoids expanding in q T /p L ∼ O(1) and thereby retains all leptonic power corrections exactly to all powers. Expanding the leptonic measurements in this limit where with a slight abuse of notation the superscript (0) now refers to the leading-power term in λ with the modified power counting in eq. (2.99), and the q T /p L argument is meant to remind us that we have not expanded in this ratio. The cross section at leading power in λ in this limit is given by where the L i (q T /p L ) arise from the modified LP leptonic measurements in eq. (2.100). The hadronic tensor does not know anything about p L , and so its power expansion is unaffected by eq. (2.99). However, as we are now keeping some terms in the leptonic measurement that we would otherwise drop in the strict q T → 0 limit, the azimuthal symmetry we might have in the strict q T → 0 limit is typically lost now, and so we do not require it. As a result, also the W 2,5 contribute at LP in eq. (2.101). Furthermore, there are now generically linear power corrections to eq. (2.101) from both L (1) i and W (1) i .

Generic fiducial power corrections
Since eq. (2.101) relies on counting p L /Q ∼ λ it is not valid for p L ∼ Q. Hence, to cover the full leptonic phase space, we have to satisfy two competing conditions from the different regions. For q T p L ∼ Q we must not expand in p L /Q, while for q T ∼ p L Q we must count q T ∼ p L to avoid uncontrolled power corrections in q T /p L . The natural way to satisfy both requirements is to expand the leptonic measurements neither in q T nor p L and thus keep the exact leptonic tensor, (2.102) Of course, we still need to expand the hadronic tensor in q T /Q ∼ λ, and all four LP hadronic structure functions in principle contribute. For q T p L ∼ Q, eq. (2.102) obviously captures the linear power corrections as in eq. (2.97), while for q T ∼ p L Q it captures as required all leptonic power corrections as in eq. (2.101). In the following, we always use the notation dσ (0+L) to denote the inclusion of the exact leptonic tensor as in eq. (2.102).
Eq. (2.102) is our final master formula. By treating the leptonic tensor exactly, it in fact incorporates all fiducial power corrections that multiply the leading-power hadronic structure functions. The leptonic tensor does not produce small-q T logarithms, which solely arise from the hadronic tensor. Therefore, eq. (2.102) automatically resums all logarithms in fiducial power corrections to the same order as the resummation is included for the hadronic tensor. All further power corrections to dσ (0+L) are obtained by working to subleading power in the hadronic structure functions, and arise purely from subleadingpower QCD dynamics.
One might argue that we could have immediately kept the leptonic tensor exact from the start, just because there is no reason or benefit to expanding it, and so one should not. On the other hand, one might argue that by doing so one keeps a seemingly arbitrary set of power corrections in the cross section, and there is a priori no guarantee that doing so would make things better and not worse, and so one should expand the leptonic tensor in order to have a consistent power expansion for the cross section at each order in the power expansion. Both arguments are found in the literature.
Our analysis provides several formal justifications for keeping the exact leptonic tensor. First, for the common case of leptonic measurements that are azimuthally symmetric at Born level (and generic p L ∼ Q), it uniquely predicts all linear O(λ) next-to-leading power corrections in the cross section. In other words, in this limit the retained power corrections are not arbitrary but provide an unambiguous, systematic improvement in the power expansion of the cross section, including their logarithmic resummation. Second, it retains all leptonic power corrections, which as argued is mandatory to correctly obtain the actual leading-power result for q T ∼ p L . Third, the actual regions and relevant scales p L depend on the leptonic measurement and identifying them can be quite involved. Keeping the exact leptonic tensor is by far the simplest (and perhaps only sensible) way to guarantee that all such possible regions are correctly treated. Finally, it ensures that any in-between regions q T p L Q are smoothly covered.
The LP hadronic structure functions entering in eq. (2.102) are given by (2.103) The cases i = −1, 4 are a straightforward generalization of the standard inclusive factorization theorem in eqs. (2.23) and (2.24), where B a,b and S are the same beam and soft functions, and only the hard functions H i ab depend on the projection i. They are collected in appendix A.1. The W 2,5 contribution, which corresponds to the cos(2ϕ) and sin(2ϕ) angular modulations of the cross section, are proportional to a (weighted) convolution of two Boer-Mulders functions h ⊥ 1 in the transverse plane [141][142][143], where h ⊥ 1a measures the net transverse polarization of flavor a, longitudinal momentum fraction x, and given transverse momentum k T within an unpolarized proton [144]. It does not match onto leading twist-2 collinear PDFs, i.e., for Λ QCD q T , each h ⊥ 1 is suppressed by at least one power of Λ QCD /q T [145] relative to the leading-power beam functions B a in W −1,4 , which do match onto leading twist-2 PDFs. The matching of h ⊥ 1 onto subleading twist-3 PDFs was carried out in ref. [146]. On the other hand, the first contribution to W 2,5 that does match onto leading twist-2 PDFs is suppressed by q 2 T /Q 2 relative to W −1,4 [107]. For these reasons, we will neglect the i = 2, 5 contributions in our numerical results. However, it should be stressed that for q T ∼ Λ QCD they do become formally leading-power contributions.
Perturbatively, eqs. (2.102) and (2.103) allow us to resum fiducial power corrections to the same order to which the LP hadronic structure functions are known. We stress that due to the different sum over flavors with different weights H i ab , the resummation effects do not in general cancel in the ratio W −1 . This is relevant when computing the angular coefficient A 4 , corresponding to the forward-backward asymmetry, at small q T .

Uniqueness of linear power corrections
There are several loose ends in the nontechnical discussion of the previous subsection that we now tie up to establish that eq. (2.97) uniquely and unambiguously captures all linear power corrections for LP-azimuthally symmetric observables.
1. In section 2.6.1, we derive the power counting of the hadronic structure functions in table 1.
3. In section 2.6.3, we show explicitly that the linear power corrections in eq. (2.97) are unique, i.e., switching to a different basis induces only quadratic power corrections.

Power counting hadronic structure functions
To derive the λ scaling of the W i in table 1, we use SCET [109][110][111][112][113], which is the effective theory of QCD in the limit λ 1 and provides a systematic expansion of QCD in powers of λ. In SCET, each collinear sector is parametrized by two lightlike reference vectors n µ i andn µ i , where n i ·n i = 2. In our case, the relevant collinear sectors describe the incoming hadrons, the colliding partons and collinear radiation off of them. We choose the reference vectors along the beam axis in the leptonic frame as in terms of which any four-vector p µ can be decomposed as With this choice, the space perpendicular to n a,b coincides with the transverse plane. It is useful to define the transverse metric and antisymmetric tensor, (2. 106) In addition, we need to distinguish a direction in the transverse plane, which we take to be where q µ ⊥ ≡ g µν ⊥ q ν , and we remind the reader that we aligned the x axis in the leptonic frame with the transverse component of q µ .
To discuss the power counting of the hadronic structure functions in SCET, we first write t µ , x µ , z µ in terms of n µ a,b and n µ ⊥ . From their explicit expressions in the leptonic frame in eqs. (2.40) and (2.41), we have where as before = q T /Q ∼ λ and γ = √ 1 + 2 = 1 + O(λ 2 ). It is straightforward to expand eq. (2.108) in λ, Note that the relations for y µ and z µ are exact and do not receive power corrections, which is a direct consequence of the symmetry we imposed on z µ . The simple form of eq. (2.109) motivates our choice of n a,b in the leptonic frame in eq. (2.104). If we had chosen n a,b as (1, 0, 0, ±1) lab in the lab frame instead, there would be additional factors of e ±Y in eq. (2.109). The power counting of the hadronic structure functions is determined by the order in λ at which contractions of eq. (2.109) with the hadronic current are populated when expanding the hadronic currents J µ V in eqs. (2.2) and (2.4) in an explicit power expansion in λ in terms of corresponding SCET current. Schematically, where x is the current's spacetime position, α, β are spinor indices, and the second sum runs over the flavor labels q, q carried by the SCET operator, which in general are distinct from the flavors coupling directly to the vector boson. The ω 1,2 and n 1,2 are (at this point arbitrary) large label momenta and directions. The leading-power hard-scattering operator O where χ q n,ω (x) is an n-collinear field with total label momentum ω, with color indices implicit. It is defined as where ξ q n is an n-collinear quark field with flavor q, and W n is a Wilson line constructed from n-collinear gluons, such that the product W † n ξ n is invariant under n-collinear gauge transformations. The δ function picks out the total large label momentum ω. The sign conventions for ω in eqs. (2.112) and (2.113) are chosen such ω 1,2 > 0 for incoming particles [147].
In principle, eq. (2.111) also receives a contribution from a corresponding leading-power O (0) gg operator, whose matching coefficient is proportional to q µ [147]. It precisely captures the non-conserved part of the current, see eq. (2.17) and the discussion below it. Since it does not contribute to Drell-Yan for massless leptons, we do not consider it here.
When evaluating proton matrix elements of eq. (2.111), momentum conservation requires n 2,1 = n a,b and ω 2,1 = ω a,b in the case where parton a is a quark. Making use of these identifications and the fact that ω a ω b = q 2 , the hard matching coefficients for V = γ, Z, W are given by [147,148] where the vector and axial-vector contributions have the same flavor-diagonal matching coefficient C q (q 2 ) but different singlet coefficients C vf (q 2 ) and C af (q 2 ). The latter arise from closed quark loops coupling to the vector boson, and thus involve an electroweak coupling different from the external quark flavors. Here, V qq is the CKM-matrix element for q ∈ {u, c, t} and q ∈ {d, s, b} (and we take it to vanish in all other cases). Importantly, the spin structure of the leading-power hard matching coefficient is proportional to γ µ ⊥ = g µν ⊥ γ ν , and therefore satisfies Using eq. (2.109), it is easy to see that contractions with the longitudinal polarization vector ε µ 0 = z µ vanish to all orders at the level of the amplitude, This is the all-order analogue of eq. (2.93) in the limit λ 1. It follows that projections onto ε µ 0 in eq. (2.47) are only populated by matrix elements of the subleading-power currents J (j)µ V with j ≥ 1 in eq. (2.110), and are penalized by at least one power of λ. This implies that only W −1,2,4,5 , which do not involve longitudinal polarizations, can scale as O(λ 0 ), while W 1,3,6,7 are suppressed at least by O(λ), and W 0 = 2W 00 is suppressed by at least O(λ 2 ). This completes the derivation of table 1.
For W −1,2,4,5 , our power-counting argument agrees with the well-known scaling of the leading contributions given by eq. (2.103), while for W 0 and W 1 it reproduces the known scaling at fixed O(α s ) [107]. For the remaining W 3,6,7 , our argument provides a lower bound on the degree of power suppression. To our knowledge, this is the first time that the scaling of W 3,6,7 at small q T has been explicitly considered for generic currents.
We also point out that starting from eq. (2.109), it is straightforward to identify the subleading-power SCET currents that populate a given W i . For example, W 1,3,6,7 can only receive their leading contributions from the interference of J (1)µ V with itself due to eq. (2.116). The hard-scattering operators to O(λ 2 ) relevant for colorsinglet production have been constructed in refs. [149][150][151] using the approach of helicity operators [148,152], and the list of operators contributing to J (1)µ V is fairly short. Due to the explicit power suppression from the current, it should be possible to derive factorization theorems for these W i in the q T Q limit using SCET. This would be relevant e.g. to understand the degree to which resummation effects are universal between W i and W −1 , and hence to what extent they cancel in predictions for the angular coefficients A i . A conjecture for the factorization of W 1 at small q T was given in ref. [153], and it would be interesting to analyze it using the systematic organization of subleading operators in SCET.

Vanishing
We next discuss the absence of linear power corrections in W −1,4 , cf. eq. (2.96). The projectors defining W −1,4 involve x µ , which in principle receives a linear power correction, see eq. (2.109). However, this O(λ) correction is proportional to n a + n b and thus orthogonal to the leading-power SCET current in eq. (2.111) due to eq. (2.115), similar to the longitudinal polarization vector discussed above. We therefore have up to quadratic power corrections, The question then reduces to why −g ⊥µν W µν and 2i ⊥µν W µν do not receive linear power corrections relative to the contribution from the squared LP current. It is well known that for e + e − → dijets event shapes such as thrust, the leading O(λ) corrections vanish [150,[154][155][156][157]. The explicit proof in refs. [150,157] relies on invariance under rotations about the axis defined by the lightlike directions that parametrize the collinear sectors for the outgoing jets. The analogous statement here is that −g ⊥µν W µν and 2i ⊥µν W µν are indeed invariant under rotations about n µ a −n µ b . To see that this implies the absence of linear power corrections, we discuss the possible sources of power corrections in turn: 1. Subleading hard-scattering operators were shown not to contribute to the thrust spectrum at O(λ) in ref. [150], using the rotational symmetry. While thrust is described by SCET I and q T is a SCET II observable, the operator basis involving only collinear fields is identical and has manifest crossing symmetry, so for these contributions the argument carries over to the case of Drell-Yan.
2. On the other hand, contributions involving soft fields differ between SCET I and SCET II , and occur both through subleading hard-scattering operators and subleading Lagrangian insertions. For SCET I -like event shapes the vanishing of such terms at O(λ) was demonstrated in refs. [150,157]. The analysis of O(λ) terms in SCET II is more difficult due to the non-locality of the theory and existence of O(λ 1/2 ) operators [158]. In ref. [159], subleading-power Lagrangians and hard-scattering operators involving soft fields in SCET II are constructed, and it is demonstrated that soft O(λ) contributions are absent for the inclusive Drell-Yan small q T spectrum, including the forward-backward asymmetry.
3. For the choice of n a,b in eq. (2.104), the measurement function for q T is the vectorial sum of the perpendicular momenta of all particles in the hadronic final state. Unlike the case of e + e − event shapes, the sum factorizes into n a -collinear, n b -collinear, and soft contributions without approximation, so power corrections from the q T measurement are absent. Since fundamentally the hadronic structure functions only depend on the Lorentz invariants in eq. (2.13), the measurement can be marginalized over the azimuthal angle of q T and thus preserves the rotational symmetry.
4. A source of power corrections absent in the e + e − case are the Born measurements on Q and Y that set the arguments of the PDFs. As has been discussed in detail in refs. [101,160], these give rise to new nonperturbative functions such as derivatives of the PDFs at subleading power in q T . It can easily be seen from the exact result that these corrections are quadratic, Recall that W incl = W −1 + W 0 /2 and we already showed that W 0 ∼ λ 2 , so the absence of linear power corrections for W incl and W −1 is equivalent. We can thus conclude that W −1 and W 4 do not receive linear power corrections.

Choice of tensor decomposition is O(λ 2 )
In section 2.3, we defined a set of reference vectors t µ , x µ , y µ , z µ to perform the tensor decomposition into hadronic structure functions, which turned out to be equivalent to the CS frame (using the boost definition). The x µ , y µ , z µ were uniquely determined by imposing eqs. (2.30) and (2.32), but these constraints are not technically required. In general, we can also pick a different set x µ , y µ , z µ of orthonormal, spacelike reference vectors. These in turn define a vector-boson rest frame related to the CS frame by a rotation r ∈ SO(3) that in general depends on q, (2.119) The corresponding hadronic structure functions W i are related to the W i by a corresponding orthogonal transformation (2.120) In terms of the W i , the fully differential cross section is given by CS n a i v e lab frame rest frame Figure 3. Sequence of boosts definining different vector-boson rest frames. In the general case (gray), we first boost by a rapidity y along the lab-frame beam axis, and then boost into the rest frame. Different choices of y lead to a relative Wigner rotation of the resulting rest frames. Relevant special cases are y = Y (Collins-Soper frame, red), y = 0 (the naive rest frame, blue) and y → ±∞, (Gottfried-Jackson frame, green).
Note that the parametrization of the lepton phase space used to evaluate the L i is irrelevant here. Of course, if corresponding angles θ , ϕ are considered, the corresponding spherical harmonics are related by the same rotation R ij (q), such that their angular coefficients are given by ratios of the W i .
First, let us point out that there is never any frame ambiguity to the order in the power expansion we are working in, because to the working order different frame choices simply amount to a specific choice of basis or coordinate system, which the final result cannot depend on. An ambiguity can only arise in the higher-order terms that are partially retained and partially neglected, which in general do depend on the frame choice.
To remove the trivial effect from a mere basis choice at LP, we start from a common LP rest frame at q T = 0, which is the leptonic frame. A convenient way to parametrize the different possible frames for nonzero q T is by the sequence of boosts starting from the lab frame as shown in figure 3. Specifically, we first boost by a rapidity y along the beam direction and then directly into the rest frame. Some cases of interest are the "naive" rest frame, obtained by performing a single direct boost from the lab frame (y = 0) into the rest frame, the CS frame (y = Y ), the Gottfried-Jackson (GJ) frame defined by the limit y → −∞, or the GJ frame obtained by taking y → +∞. Another way to conceptualize these frames is through the angles between their respective z µ axes and the momenta of the incoming protons. In the CS frame, those angles become precisely equal for massless hadrons, see figure 1. In the GJ frame, the z axis is aligned with the direction of P a , while the naive case falls in between for Y > 0.
It is easy to see that for q T → 0, the two boosts for any frame collapse into a single boost from the lab frame to the leptonic frame. Since for q T → 0 all frames coincide, we have so the scaling properties in table 1 also hold for the W j . Following our previous analysis, we would now discard all power-suppressed W j structure functions, evaluate the remaining ones (j = −1, 2, 4, 5) at leading power, and dress them with exact leptonic tensor components, (2.124) Since all i = −1, . . . , 7 are now populated by R −1 ij , one might think that the different choice of tensor decomposition amounts to a linear power correction compared to the left-hand side as given in eq. (2.102), because the R ij differ from unity by O(λ), but as indicated it is only of O(λ 2 ).
To show that the induced difference is indeed only O(λ 2 ), first note that rotations around the z axis amount to a trivial shift in ϕ. This induces an O(λ 2 ) difference at cross-section level due to our assumption of azimuthal symmetry at leading power in eq. (2.88). Hence, it is sufficient to consider the SO(2) subgroup of rotations around the y axis parametrized by one remaining Euler angle α, with c α ≡ cos α and s α ≡ sin α. By a straightforward calculation, the W i can be shown to transform under the following representations of SO(2): It is a simple exercise in special relativity to show that the resulting Wigner rotation of a generic frame defined by boost y relative to the CS frame with y = Y is where = q T /Q ∼ λ and γ = √ 1 + 2 . For any Y, y, eq. (2.127) turns out to be bounded by the leading term, |tan α| ≤ |tanh Y −y 2 | ≤ . In particular for the GJ frame we have tan α GJ = , with which we recover the well-known result for the relation between W −1,0,1,2 in the CS and GJ frames [107]. Eq.
(2.127) shows explicitly that α ∼ λ. To see that eq. (2.124) indeed holds up to quadratic power corrections, we invert eq. (2.126) by taking α → −α and expanding in α ∼ λ, to find (2. 128) We see that under the rotation, the leading structure functions mix into themselves and into each other only by an O(λ 2 ) amount. The subleading structure functions are populated precisely by an amount commensurate with their intrinsic scaling, see table 1. Combined with the scaling of the corresponding leptonic structure functions as in figure 2, we find that the effect of any O(λ) rotation on the cross section is indeed only of O(λ 2 ) for LPazimuthally symmetric observables. It is natural to ask whether a specific frame choice should be preferred in order to also capture an optimal set of terms at O(λ 2 ). It is well known that leading-logarithmic terms are absent in W 1 in the CS frame [107]. This is natural from the point of view that the CS z axis does not receive power corrections, reducing "spill-over" of terms from the leading-power currents. A reduced size of power corrections by a symmetric choice of frame has also been found for 0-jettiness T 0 [160][161][162], albeit for a somewhat different physical reason. In section 4.4 we will find some numerical evidence that the CS decomposition also reduces the size of power corrections in the φ * spectrum. Taken together, this suggests that the CS frame might indeed be the optimal choice, although the size of the frame-dependent O(λ 2 ) corrections should still be assessed.

Relation to the literature
Early approaches to resummation effects on the Drell-Yan cross section differential in the lepton kinematics [65][66][67] typically picked the CS angles to parametrize the decay phasespace integral, but did not discuss the ambiguity inherent in this choice or the relative size of power corrections. In these approaches, the CS frame primarily serves as a tool to enable generic lepton-differential observables. In ref. [107], the logarithmic structure of W −1,0,1,2 at low q T in collinear factorization was discussed in detail, but the implications for the structure of power corrections to e.g. the fiducial cross section were not explored.
An implementation of generic lepton-differential observables in q T resummation was presented in ref. [73] and more recently in ref. [74] based on a parton-level Monte-Carlo generation of the leptonic final state. There, the choice of rest frame used in earlier results was recast as a q T -recoil prescription for how to distribute the nonzero q T between the colliding partons in the rest frame where the leptonic decay is evaluated. They also showed that the ambiguity in the q T -recoil prescription is in one-to-one correspondence with the ambiguity of the rest-frame choice and that it vanishes for q T → 0. They thus argue that the recoil effects are O(q T /Q) effects that cannot be unambiguously computed through the q T resummation, but some recoil prescription is nevertheless required for practical purposes to satisfy transverse momentum conservation in the parton-level generation of the leptonic decay. (Note that similar recoil prescriptions to preserve momentum conservation are also commonly used in parton-shower Monte Carlos.) Very recently, refs. [86,87] used N 3 LL perturbative baselines to fit nonperturbative models for the rapidity anomalous dimension and TMDPDFs using fiducial Z q T spectra among other data. They also retain the exact dependence of the fiducial phase space on q T , with the analytic leptonic decay matrix element contracted against the LP hadronic tensor ∝ g µν ⊥ . 9 This is essentially equivalent to an exact treatment of L −1 in our notation (up to an overall O(λ 2 ) difference in the projection itself), while L 4 does not contribute to the observables they consider. They also do not provide formal arguments for the exact treatment of the leptonic contributions.
If the fully-differential and flavor-channel dependent leptonic decay matrix element is evaluated in the vector-boson rest frame, including exact leptonic cuts, and contracted against the resummed LP hadronic tensor, either explicitly at the analytic level or during the parton-level generation via the recoil prescription, then this essentially amounts to retaining the exact q T dependence of the leptonic tensor. To the best of our understanding this is the case in refs. [73,74,86,87]. Our analysis thus provides formal justification for doing so, showing that the ambiguity is only of O(λ 2 ) and that for a large class of common measurements (those that are azimuthally symmetric at LP) it actually unambiguously predicts all linear power corrections along with their resummation. In addition, it is formally required in phase-space regions that exhibit leptonic power corrections.
In ref. [83], fiducial lepton cuts are implemented in the resummed cross section strictly on Born kinematics at q T = 0, while fiducial power corrections are obtained through the fixed-order matching. Large power corrections from the fixed-order matching were observed in the fiducial case compared to the inclusive case. From our analysis, this is explained by the linear power corrections induced by the fiducial cuts.
Sometimes a multiplicative fixed-order matching procedure is employed, as e.g. in refs. [83,88], where in order to Sudakov-suppress the fixed-order matching corrections at small q T they are multiplied by the ratio of the LP resummed contribution to its fixed-order expansion. While this procedure is unlikely to produce the correct Sudakov suppression for genuine hadronic power corrections, one might ask if it correctly "dresses" the fiducial power corrections with the LP resummation to achieve their resummation. For this to be the case, the multiplicative matching at minimum has to reproduce eq. (2.97) for the linear power corrections. Clearly, this can only happen if the multiplicative matching only involves a single (effective) hadronic structure function at a time and if it is performed fully differentially in q 2 , Y , and q 2 T . This is typically not the case. For example, the multiplicative matching in refs. [83,88] is performed at the cumulant level, and thus does not satisfy this requirement.

Resummation of leading-power hadronic structure functions
In this section, we discuss our specific resummation setup for the leading-power hadronic structure functions W (0) −1,4 in eq. (2.103). The setup follows standard procedures and is deliberately kept simple, e.g., by ignoring nonperturbative corrections or quark-mass effects [89] at small q T , allowing us to focus on the effect of resumming the fiducial power corrections in the following sections. As discussed in the previous section, the fact that their resummation can be obtained in terms of the leading-power hadronic structure functions, and that this captures all linear as well as leptonic power corrections holds independently of how precisely the LP resummation is performed.

Renormalization group evolution
Directly resumming the logarithms of q T /Q in momentum space is challenging due to the vectorial nature of q T , though by now approaches for doing so exist [76,77]. As shown in ref. [77], the correct solution of the RGE system in momentum space in terms of distributions is actually equivalent to the canonical solution in conjugate (b T ) space modulo different boundary conditions. We therefore bypass this issue, as is commonly done, by solving the RGEs in conjugate (b T ) space. Using that the beam and soft functions only depend on the magnitude b T = | b T |, eq. (2.24b) can be written as where J 0 (x) is the zeroth-order Bessel function of the first kind. We use the framework of the rapidity renormalization group [62] in conjunction with the exponential regulator of ref. [64], where the beam and soft functions in Fourier space obey the coupled RGEs where the µ anomalous dimensions on the first line have the all-order expressions Here, Γ q cusp (α s ) is the cusp anomalous dimension in the fundamental representation, which is known to four loops [163][164][165][166][167][168][169][170][171][172] (see ref. [170] for a complete list of earlier references). At N 3 LL we also require the QCD beta function to four loops [173][174][175][176]. The beam and soft noncusp anomalous dimensionsγ q B (α s ) andγ q S (α s ) are related to those of the 0-jettiness beam and soft functions through consistency relations [177] and are known to three loops [131,[178][179][180]. The all-order form of the rapidity anomalous dimension γ q ν reads This form follows from the commutativity of the µ and ν derivatives in eq. (3.2). The fixed-order boundary termγ q ν,FO is known to three loops [131,134,181]. The RG evolution factors forB q andS follow by solving the coupled systems of equations in eq. (3.2). One possible way of writing the solution is given bỹ Here, we evolve first in ν and then in µ and the final scale µ is set to the scale µ H at which the hard function is evaluated. Any other path in the two-dimensional (µ, ν) space connecting (µ H , µ B , ν B , µ S , ν S ) is viable as well, and the path independence is ensured by exactly satisfying the RG consistency relations between all anomalous dimensions, in particular by using eq. (3.4). By choosing appropriate boundary scales µ H,B,S and ν B,S , the hard function, as well as the beam and soft functions on the right-hand side of eq. (3.5) are free of large logarithms and can be evaluated in fixed-order perturbation theory, with all large logarithms resummed through the evolution kernels. Different methods to numerically evaluate the evolution kernels and integrals over anomalous dimensions entering in eqs. (3.4) and (3.5) were analyzed in detail in ref. [92], which found that their numerical precision becomes relevant at N 3 LL. Here we use the approximate unexpanded analytic results, which provide sufficient numerical precision for our purposes, with the explicit expressions up to N 3 LL given in ref. [92]. The anomalous dimensions and fixed-order boundary conditions required for the resummation at N 3 LL are collected in our notation in refs. [177,182].

Canonical scales and nonperturbative prescription
The canonical scales of the hard, beam, and soft functions, and the rapidity anomalous dimension in b T space are given by where b 0 ≡ 2e −γ E ≈ 1.12291. By evaluating the functions in the factorization theorem at their canonical scales and evolving them to common scales (µ, ν), all logarithms of In ref. [77] it was shown that the resummation in b T space with the canonical scales in eq. (3.6) is in fact equivalent to the exact solution of the RGE in momentum space, i.e., it reproduces the canonical, distributional logarithms in q T , except for the fact that one effectively uses a shifted set of finite terms in the boundary conditions (similar to the difference between renormalization schemes). We exploit this and require that for q T Q, eq. (3.6) is exactly satisfied, such that the resummed q T spectrum in this region is obtained from the numerical inverse Fourier transform of the canonical b T -space result.
With the choice in eq. (3.6), the beam and soft functions and the rapidity anomalous dimension become sensitive to nonperturbative effects at large b T Λ −1 QCD . The Landau pole can be avoided by choosing µ 0 (and µ B , µ S ) as a function of b T such that it asymptotes to a perturbative value at large b T , see e.g. ref. [80] for a concrete implementation. Alternatively, a global replacement of b T by a suitable function b * (b T ) may be performed at the level of the cross section, where b * (b T ) itself is bounded by some perturbative value b max 1/Λ QCD [56,57]. In either case, the mismatch to the full result is then in principle captured by a corresponding nonperturbative contribution. Recently, it was shown that the fullγ i ν may also be determined from lattice calculations [183][184][185][186][187], and that estimates of the first subleading power in b T Λ QCD can also be related to the gluon condensate [188].
Since nonperturbative effects in the region q T ∼ Λ QCD are not the main focus of this work, we use a simpler prescription to ensure that α s remains perturbative. Specifically, we freeze out both the running coupling and the PDFs entering the hadronic structure functions W We choose the smooth function µ fr (µ) governing the freeze-out as (3.8) In practice, we pick Λ fr = 1 GeV for our central results. The behavior of α fr s at low scales is illustrated in the left panel of figure 4. This choice constitutes a (fairly crude) model for the large b T behavior ofγ i ν that is sufficient to regulate the large b T region, and formally amounts to a power correction in Λ QCD q T . We similarly ignore contributions of O(b T Λ QCD ) in the beam and soft function boundary conditions, beyond the ones encoded in our global choice of µ fr (µ). This is also consistent with neglecting the hadronic structure functions W 2,5 involving Boer-Mulders functions in the regime Λ QCD q T Q altogether, see the discussion below eq. (2.103).

Fixed-order matching and profile scales
We extend the description of the cross section to the fixed-order region q T ∼ Q by an additive matching to the fixed-order result via profile scales, dσ = dσ (0+L) (µ res ) + dσ FO (µ FO ) − dσ (0+L) (µ FO ) . Here the µ res argument of the first term indicates that we evaluate the resummed LP hadronic structure functions in dσ (0+L) using the resummation profile scales given below. The argument µ FO in the last term indicates that they are instead evaluated at common fixed-order scales µ FO , which can be done directly in momentum space. The last term effectively acts as a differential subtraction term for the full fixed-order cross section dσ FO in the second term, such that the difference in square brackets is a nonsingular power correction. We will refer to the outcome of eq. (3.9) as N 3 LL (0+L) +NNLO 0 when the resummed LP hadronic structure functions at N 3 LL are combined with the exact leptonic tensor as discussed in section 2.5, and matched to the fixed O(α 2 s ) NNLL 0 result. When instead evaluating the leptonic tensor at strict LP, we refer to it as N 3 LL (0) +NNLO 0 , in which case the fiducial power corrections are only included through the fixed-order matching. The analogous notation is used at lower orders.
Approaching q T ∼ Q, the q T resummation must be turned off to ensure the delicate cancellations between singular and nonsingular contributions and to properly recover the correct fixed-order result for the spectrum, i.e., such that the difference between the first and last terms in eq. (3.9) smoothly vanishes for q T ∼ Q. To achieve this, we use the method of profile scales [189,190]. We use the hybrid profile scales constructed in ref. [80], which depend on both b T and q T , and undergo a continuous deformation away from the canonical b T scales in eq. (3.6) as a function of the target q T value, schematically, (3.10) We note that µ 0 does not need to asymptote to µ FO towards large q T because its effect on the matched result is already turned off as ν S → ν B . In this limit, the first and last term in eq. (3.9) exactly cancel, leaving the fixed-order result dσ FO . We choose central scales as where f run is a hybrid profile function given by f run (x, y) = 1 + g run (x)(y − 1) . (3.12) It controls the amount of resummation by adjusting the slope of the scales in b T space as a function of q T /Q via the function (3.13) As a result, for q T ≤ x 1 Q, the slope is unity yielding the canonical resummation, while for q T ≥ x 3 Q, the slope vanishes so the resummation is fully turned off. In between, the slope smoothly transitions from one to zero, which transitions the resummation from being canonical to being turned off. This is illustrated in the right panel of figure 4. Note that in contrast to ref.
[80], we do not require a deformation away from the canonical scales to regulate the Landau pole at large b T , but instead rely on the replacement in eq. (3.7). For our central results, we use the transition points (x 1 , x 2 , x 3 ) = (0.3, 0.6, 0.9).

Estimate of perturbative uncertainties
We identify several sources of perturbative uncertainties, which we estimate as follows.
In the limit q T Q, the perturbative uncertainty is driven by the combined uncertainty from truncating the expansion of the soft, beam, and rapidity anomalous dimensions. To estimate them, we adopt the set of SCET II profile scale variations introduced in ref. [191] for a SCET II -like jet veto. They are given by where each of the four variation exponents can be v i = {+1, 0, −1}, and the amount of variation is governed by i.e., as q T → Q, this source of perturbative uncertainty is turned off together with the resummation as a whole. The central scale choice corresponds to (v µ S , v ν S , v µ B , v ν B ) = (0, 0, 0, 0), and a priori there are 80 possible different combinations of the v i . Since the arguments of the resummed logarithms are ratios of scales, some combinations of scale variations will lead to variations of these arguments that are larger than a factor of two, and are therefore excluded. After dropping these combinations we are left with 36 different scale variations. By taking their maximum envelope, we obtain an estimate of the resummation uncertainty ∆ res . Note that independent variations of µ H need not be considered as part of ∆ res because the corresponding change in the argument of the resummed logarithms is already covered by eq. (3.14). Second, we estimate the fixed-order perturbative uncertainty ∆ FO from the maximum envelope of overall variations of µ FO by a factor of two. These variations are inherited by all the resummation scales in eq. (3.11), so they leave the resummed logarithms invariant. Third, we estimate the inherent uncertainty ∆ match in our matching procedure eq. (3.16) Finally, we consider two independent variations of Λ fr = {0.8, 1.5} GeV away from our central choice Λ fr = 1 GeV as a rough estimate of the uncertainty ∆ Λ in our nonperturbative prescription. Combining all sources of uncertainty in quadrature, we take as an estimate of the total (perturbative) uncertainty on our results.

Resumming fiducial power corrections
As discussed in section 2.5, fiducial power corrections arise entirely from the leptonic tensors L i (q, O, Θ), and accordingly can be treated exactly in the factorization by keeping the L i exact. In this section, we consider three applications to discuss this mechanism in more detail, namely the q T spectrum with fiducial cuts (section 4.2), the lepton transverse momentum distribution (p T ) (section 4.3), and the φ * observable (section 4.4). In all cases, we consider our primary examples of Z → + − or W → ν.

Numerical inputs and computational setup
All our numerical results in this and the following sections are obtained using the following setup. We use the CT18NNLO [192] PDF set, and correspondingly use the three-loop running to obtain the numerical value of α s at any required scale with α s (91.1870) = 0.118 as starting point. We use the same PDF also at lower orders, which is consistent and allows us to exhibit the genuine size of perturbative corrections.
For the resonant W and Z propagators, we work in the fixed-width pole scheme. We use the following electroweak parameters [193]   For the electroweak couplings, we use the G µ scheme, with the values for α em and sin 2 θ w obtained from m W , m Z , and G F as All factorized cross sections, both at fixed order and including resummation up to N 3 LL accuracy as described in section 3 are obtained from the C++ library SCETlib [114]. By default we use the CS tensor decomposition, and LP cross sections including fiducial power corrections are denoted as σ (0+L) , while those at strict LP without fiducial power corrections are denoted by σ (0) . We have also implemented alternative tensor decompositions using eq. (2.124), in particular the one that corresponds to the GJ frame, and will denote cross sections evaluated using this choice as σ (0+L) GJ . SCETlib uses the Cuba 4.2 library [194,195] for adaptive multi-dimensional integration over Q and Y , combined with q T integrals whenever they cannot be performed analytically. To perform oscillatory Bessel integrals for inverse Fourier transforms we use a double-exponential method for oscillatory integrals [196][197][198].
The integral over the leptonic phase space appearing in eq. (2.60), is carried out semi-analytically as follows. We focus on binned observablesΘ(q, θ, ϕ) and assume that the differential measurement O on the decay products is being integrated over. We assume that the measurement cuts and binsΘ(q, θ, ϕ) evaluate to 1 when all cuts are passed and 0 otherwise, i.e., we take it to be a product of θ functions. An explicit dependence on θ or ϕ, e.g. to apply angular projections, can also easily be accommodated, but this is not needed for the results obtained in this paper. For given values of q and ϕ, we 10 The pole-scheme values are converted from the on-shell ones using first determine all intervals in cos θ that pass the given cuts. Notably, for all observables considered here (p T 1,2 , η 1,2 , φ * , and any of their combinations), the interval boundaries can be evaluated analytically even for nonzero q T . The integral over cos θ over these intervals is then carried out analytically. The remaining integral over ϕ is performed by adaptive numerical integration. In practice, the sum over hadronic structure functions i, see eq. (2.102), can be pulled under the integral in eq. (4.4) because the structure functions only depend on the given value of q, so the decay phase-space boundaries only have to be determined once. Typical evaluation times even for complicated phase-space volumes are in the few-millisecond range on a single Intel R Core TM i5-7200U CPU @ 2.50 GHz for a target relative numerical precision of 10 −7 . The algorithm is not restricted to the leading-power structure functions, but can also be used standalone with generic hadronic structure functions that are provided. Fixed-order results for q T and φ * at LO 1 and NLO 1 for the relevant Born+1-parton cross sections are obtained from MCFM8 [199][200][201]. These results are used in the fixedorder matching. In addition, they are used to obtain q T (or φ * ) integrated cross sections at NLO 0 and NNLO 0 by combining them with q T (or φ * ) subtractions including fiducial power corrections supplied by SCETlib. This setup is discussed in more detail in section 5. We stress that the inclusion of fiducial power corrections in the subtractions is essential to obtain numerically stable results down to very small q T and φ * and for p T near the Jacobian peak.

q T spectrum with fiducial cuts
We first discuss the impact of fiducial cuts on the Drell-Yan q T spectrum. We consider the standard kinematic selection cuts of requiring a minimum transverse momentum p min T and maximum rapidity η max of the final-state leptons, where p T,i and η i with i = 1, 2 are the transverse momentum and pseudorapidity of the two leptons.

Origin of power corrections
The p min T cut was already discussed in detail in ref. [104]. Here, we briefly review the key steps and results, and in addition discuss the rapidity cut. A useful parametrization of the total momentum q µ and the lepton momenta p µ 1,2 in the lab frame is where m T = (Q 2 + q 2 T ) 1/2 . As before, we neglect the lepton masses and align the total transverse momentum q T with the x-axis. We denote the azimuthal angle of the first lepton in the lab frame as ψ to distinguish it from the CS angle ϕ, and write its rapidity as y 1 = Y + ∆y. Momentum conservation determines p µ 2 , and fixes the transverse momenta and rapidities of the leptons to (4.7) For compactness, here we do not substitute p T,1 in the expressions for p T,2 and η 2 . The two-particle phase space defined in eq. (2.10) then takes the form The integrated phase space with the cuts in eq. (4.5) can now be written as The integrand in eq. (4.9) depends on q T only through the combinations q 2 T and q T cos ψ. Thus, the expansion of the integrand in the limit q T Q can only yield linear fiducial corrections if the ψ integral does not average out. This is equivalent to requiring that the cuts break azimuthal symmetry, as otherwise the ψ integral can always be trivially carried out, such that all odd powers of q T cos ψ integrate to zero and only quadratic corrections in (q T /Q) 2 arise. Inclusive measurements are a special case, as without cuts Φ L (q) can only depend on q 2 .
To see this mechanism explicitly for cuts on p T and η, we expand the lepton transverse momenta and rapidities in eq. (4.7) in q T /Q ∼ λ, (4.10) All observables in eq. (4.10) have a well-defined, nonvanishing LP limit as q T → 0, and the first correction is proportional to q T cos ψ. Since at q T = 0, ψ and ϕ coincide, we immediately find that the fiducial q T spectrum obeys azimuthal symmetry at leading power, so the discussion in section 2.5.1 applies. Naively, one might also expect that all linear fiducial corrections vanish upon integration over ψ. However, this is spoiled by the minimum and maximum in the θ functions in eq. (4.9), as can be easily seen for the p min T cut. For cos ψ > 0, one has p T,1 > p T,2 , and thus the θ function in eq. (4.9) only restricts p T,2 . Vice versa, for cos ψ < 0 it is p T,1 that is constrained. This leads to two different integrands of the ψ integral in the two integration regions, leading to a nonvanishing ψ integral. This shows that the azimuthal symmetry is explicitly broken at O(λ) leading to linear fiducial power corrections. However, it also shows that if one were to only apply cuts to one of the leptons while being fully inclusive in the other, no linear power corrections from the cuts would arise, since the ψ integral would average out when integrated against the g −1,2,4,5 (θ, ϕ) (using again that the difference between ψ and ϕ is itself of order q T ).
The situation is more complicated for the rapidity cut. Determining the transition point of the maximum in the corresponding θ function in eq. (4.9), i.e. the value ψ tp for which |η 1 | = |η 2 |, we find that If |cos ψ tp | ≥ 1, then the θ function in eq. (4.9) only restricts either |η 1 | or |η 2 | but not both for the whole ψ range. In this case, the rapidity cut does not break azimuthal symmetry. For small but nonvanishing values of q T , the Q/q T scaling in eq. (4.11) can be overcome by a sufficiently small value of the vector-boson rapidity Y . To be precise, eq. (4.11) has a solution in the physical range |cos ψ tp | < 1 when . (4.12) Note that the η 1 constraint always requires that |Y + ∆y| ≤ η max . Furthermore, we are only interested in the q T Q limit, which implies that eq. (4.12) only becomes important when |Y | 1. Hence, linear fiducial corrections will only arise in the region while for q T q tp T only quadratic power corrections arise. Note that in the region q T ∼ q tp T Q the standard power counting breaks down, as one has to simultaneously expand |Y | ∼ q T /Q 1. This is an example of a leptonic fiducial power correction as discussed in section 2.5.2, where it is crucial to keep the lepton phase space exact to correctly account for both small scales q T /Q and |Y |. In practice, the size of this region is of O(|Y |) and thus small by construction, and hence its contribution to the cross section when integrated over or binned in Y is suppressed as well.
To illustrate and validate our observations, we have numerically implemented the exact phase space with cuts in eq. (4.9). In figure 5, we show the relative difference of the dilepton phase space as a function of q T /Q for the cut η max = 2.4 (left) and the combined cut p min correction has no simple scaling behavior. Finally, for the relatively large value Y = 1 one has quadratic corrections essentially throughout the whole q T spectrum (green dotted line). Also note that, in general, the corrections to the phase space are rather small, as even for Y = 0.1 they do not exceed 1% for q T /Q < 0.1. In the right plot in figure 5, we apply both η max and p min T cuts. For small values of Y , the p min T cut is a stronger constraint on the phase space than the η max cut. Hence, the two curves for Y = 0 (red solid) and Y = 1 (blue dashed) are equal, as the p min T constraint is independent of rapidity. For large Y ∼ η max , the rapidity cut dominates over the p T cut, as illustrated for Y = 2 (green dotted), and one only has quadratic power corrections.

Numerical results
Having explicitly demonstrated that linear power corrections to the q T spectrum arise from fiducial cuts on the final-state leptons, we now verify that they can indeed be captured in the factorization theorem by keeping the lepton kinematics exact, as discussed in section 2.5.1. In the following, we will always consider the fiducial cut p T ≥ 25 GeV , |η | ≤ 2.4 , (4.14) as employed in the CMS Drell-Yan measurement at 13 TeV [15].
In figure 6, we show the q T spectrum for Q = m Z at NLO 0 , both without (left) and with (right) fiducial cuts. In both figures, the red points illustrate the full NLO 0 result, while the solid blue line shows the result at σ (0) (left) and at σ (0+L) (right). The various dotted and dashed lines show the differences between the full result and different singular limits. In the inclusive case (left panel), there are no linear power corrections, and thus σ − σ (0) (green, dashed) scales quadratically in q T , as expected. With fiducial cuts (right panel), σ −σ (0) (gray, dotted) clearly suffers from linear power corrections, and as explained before, these linear corrections can be accounted for by keeping the leptonic tensor exact. This is illustrated by the green-dashed line, which shows the difference σ − σ (0+L) between the exact and NLP result, and only depends quadratically on q T . In particular, the size of these corrections is comparable to the quadratic corrections in the inclusive case. The orange, dot-dashed curve shows the difference σ (0+L) GJ − σ (0+L) between two choices of the tensor decomposition, corresponding to the CS frame and the GJ frame. This difference scales quadratically in q T , confirming that the ambiguity from the choice of tensor decomposition is quadratically suppressed. Moreover, we observe that this ambiguity is numerically much smaller than σ (0+L) itself, indicating that it may be completely negligible in practice.
It is also interesting to study the impact of the NLP corrections on the resummed q T spectrum. In the top-left panel of figure 7, we show the difference between the LP and the exact q T spectrum at NLO 0 (blue, short-dashed) and NNLO 0 (red, long-dashed). For reference, the gray line shows our best prediction σ * at N 3 LL (0+L) +NNLO 0 , scaled down to 5% of its original size. At both NLO 0 and NNLO 0 , the power corrections diverge as q T → 0 due to the overall 1/q T behavior (compared to 1/q 2 T at LP). The opposite signs at small q T also illustrates the poor perturbative convergence in this regime.
In the top-right panel of figure 7, we show the difference between the LP q T spectrum and the resummed and matched q T spectrum, at NLL (0+L) (green, dotted), NNLL (0+L) +NLO 0 (blue, dashed) and N 3 LL (0+L) +NNLO 0 (red, solid). Since the resummation includes the linear power corrections, the divergence as q T → 0 is cured, and we observe very good perturbative convergence between the different resummed predictions.
Finally, in the bottom panel of figure 7 we show the difference between the NLP and the exact q T spectrum at NLO 0 (blue, short-dashed) and NNLO 0 (red, long-dashed), again including our best prediction σ * for reference. Since all terms diverging as 1/q 2 T or 1/q T are included in σ (0+L) , this difference is finite as q T → 0, and the overall size is much smaller compared to the top left panel, indicating that corrections beyond the linear NLP limit can be safely included at fixed order.
Given the large effect that the resummation of fiducial power corrections has on the q T spectrum at small q T , it would be very interesting to investigate whether a net resummation effect on the total fiducial cross section remains after integration over q T , as would be Our best prediction dσ * for the total spectrum at N 3 LL (0+L) +NNLO 0 is indicated as a light gray line for reference, scaled down to 5% of its original size. In the bottom row we restrict to the remaining nonsingular (quadratic) power corrections, which are finite for q T → 0 (note the difference in vertical scale). relevant for PDF determinations from fiducial Z and W rapidity spectra. It has been argued that recoil ambiguities might prevent a first-principle calculation of these effects [202], but since the ambiguity from the choice of tensor decomposition is fairly small, and vanishes as soon as the resummation is fully turned off, it is not unlikely that a surviving net effect is unambiguous and can be calculated.

Lepton p T spectrum
We next study the distribution in the lepton transverse momentum p T . To be concrete, we consider W ± → ± ν , for which p T is an essential observable. For simplicity, we do not consider any additional fiducial cuts. This serves as a prototypical example of the appearance of leptonic power corrections near a radiation-sensitive edge of Born phase space as discussed in section 2.5, which in this case happens near the Jacobian peak at p T ∼ Q/2.

Origin of power corrections
Using the parametrization of the lepton momenta in terms of CS angles in eq. (2.55), the lepton p T can be expressed as where s θ ≡ sin θ, s ϕ ≡ sin ϕ, c ϕ ≡ cos ϕ. We also introduced the variable κ = 2p T /Q to parameterize the distance of p T from the Jacobian peak at p T = Q/2, which will be useful in the following. Eq. (4.15) can be easily solved for s θ , with physical solutions constrained by 0 ≤ s θ ≤ 1. In the following, we restrict ourselves to the case κ > , which will be the relevant region to describe large p T . In this case, the only physical solution for s θ is given by where for brevity we introduced the abbreviations (4.17) The leptonic structure functions defined in eq. (2.60) are then given by (4.18) Note that eq. (4.18) still shows the full dependence on ϕ and s θ , such that one could easily reinstate fiducial cuts. Eq. (4.16) yields two solutions for θ ϕ , related by θ → π − θ. Since the g i (θ, ϕ) for i = 1, 4, 6 are odd under this transformation, eq. (4.18) immediately vanishes for these cases. In all other cases, the g i are even under this transformation, such that we obtain where θ ϕ can be either of the two physical solutions defined by eq. (4. 16).
It is now straightforward to expand in to study the q T → 0 limit, which yields (4.20) The constraint κ ≤ 1 reflects the strict bound p T ≤ Q/2 in the Born limit. We also recover that at LP only the i = −1, 0 contributions survive. The i = 4 contribution, which in principle can contribute at LP and gives rise to the forward-backward asymmetry, vanishes due to the symmetry of p T under θ → π − θ. Before proceeding, it is instructive to illustrate the phase space differential in p T , which is closely related to L i (q, p T ) and provides a bound on the L i (q, p T ) since all g i (θ, ϕ) are bounded. It can be evaluated more easily using the parametrization of the lepton momenta given in eq. (4.6). After some effort, one obtains Here, α = √ 1 + 2 − , and K(x) is the complete elliptic integral of the first kind. The appearance of three distinct regions can easily be understood from eq. (4.19): For κ > 1/α, the θ function in eq. (4.19) becomes incompatible with c ϕ ≤ 1, while for κ < α it imposes no constraint in addition to c ϕ ≥ −1.
In figure 8, we show the differential phase space in p T in eq. (4.21) for Q = m W . In the left panel, we show the exact phase space as a function of p T for fixed q T = 8 GeV (red solid), with the gray vertical lines indicating the edges of the different regions in eq. (4.21). In the q T = 0 LP limit (blue dashed), they collapse to the kinematic Born limit p T ≤ Q/2. In the right panel, we fix p T = 40.1 GeV very close to the Born edge, and show the phase space as a function of q T . The exact result is shown by the red curve, the q T = 0 LP limit by the blue dashed line, and their difference by the green dotted line. The thick vertical line at q T = Q − 2p T ≡ p L shows the transition to the second region of eq. (4. 21). For sufficiently small q T p L , we see a clear (quadratic) power suppression, while near and above this value the power corrections become O(1). (The sharp dip in the green line is just an artefact of the logarithmic scale and the green line changing its sign.) Clearly, expanding in ∼ λ is only well-defined if κ 1, i.e., away from the Jacobian peak p T Q/2. This is already evident from the divergence of eq. (4.20) as κ → 1. As long as κ 1, is the only small scale in the problem, which justifies expanding in and leads to at most linear fiducial power corrections. On the other hand, close to the Jacobian peak, the distance p L = Q − 2p T emerges as an additional small scale, and the naive expansion in q T is actually an expansion in q T /p L , which is only allowed for q T p L but breaks down for p L ∼ q T . To illustrate this explicitly, we can expand in the regime p L ∼ q T by simultaneously counting both scales as small. To do so, we take where as before = q T /Q ∼ λ, such that formally q T /(Q−2p T ) ∼ 1. With this replacement, we can expand eq. (4. 19) in λ, This vanishes for all g i odd in ϕ, which only leaves i = −1, 0, 2, 3. This should be contrasted with the naive LP result in eq. (4.20), which only receives contributions from i = −1, 0. The i = 2 contribution is proportional to the double-Boer-Mulders effect, which we can neglect, see the discussion below eq. (2.103). For i = 3 we have W 3 ∼ O(λ), see table 1, which thus yields a linear power correction. Hence, we find the interesting effect that the proximity to the Jacobian peak induces sensitivity to new hadronic structure functions at O(λ), which do not contribute at O(λ) away from the peak region. From eq. (4.23) it is evident that naively expanding in q T near the Jacobian peak would amount to expanding in q T /p L , which is not allowed. However, eq. (4.23) is only valid near the peak, because by counting p L /Q ∼ λ we have expanded away the dependence on κ = 1 + O(λ), which is not allowed away from the peak. Hence, to cover the full range of p T , we must not expand in p L , while near the peak we must count q T ∼ p L to avoid inducing uncontrolled leptonic power corrections in q T /p L . Clearly, the simplest way to satisfy both requirements is to not expand at all and keep the exact result corresponding to eq. (4. 19).
Finally, note that the breakdown of the naive power expansion around p T = Q/2 does not immediately affect the leptonic tensor if we only consider a fiducial cut p T ≥ p min T , since we can evaluate it as  Figure 9. Lepton transverse momentum spectrum for on-resonance W + production at the LHC at fixed order (left) and including the resummation of fiducial power corrections to N 3 LL (right). The horizontal axes shows the distance to the Jacobian peak at p T = m W /2. Thus, the leptonic power corrections in this case scale as q T /(Q − 2p min T ), and so as long as p min T Q/2, the effect of p min T can be treated as a linear fiducial power correction as discussed for the q T spectrum with fiducial cuts in section 4.2.

Numerical results
There are two key insights from our analysis of the differential p T phase space. First, the p T spectrum near the Jacobian peak is directly sensitive to the small transverse momentum q T of the decaying vector boson. This causes fixed-order predictions to become unreliable in this region, which is a well-known effect. Second, the strict q T → 0 limit by itself cannot describe the p T spectrum in this region, which means the strict LP q T resummation is also insufficient. Both problems are cured simultaneously by combining the exact leptonic tensor, which encodes the exact decay kinematics and automatically retains all leptonic power corrections, with the q T -resummed hadronic tensor, thus allowing us to obtain physical predictions around the Jacobian peak.
We illustrate this in figure 9 for the p T spectrum in W + → + ν decays, where we show the spectrum both at fixed order (left) and after resummation including fiducial power corrections (right). In both panels, the horizontal axis shows the distance of p T to the Jacobian peak at p T = m W /2, and to avoid smearing out the peak we consider the spectrum at a fixed point Q = m W . The fixed-order spectrum (left) is shown at LO 0 (green dotted), NLO 0 (blue dashed), and NNLO 0 (red solid). The LO 0 result corresponds to Born kinematics and clearly shows the kinematic edge at p T = Q/2. Starting at NLO 0 , the W boson can have nonvanishing q T , which opens up the phase space beyond the edge. However, in the vicinity of the edge, the fixed-order predictions become unstable due to the sensitivity to small q T , which is clearly visible by the diverging NLO 0 and NNLO 0 curves, and in particular by the sign change between NLO 0 and NNLO 0 at p T ≈ Q/2.
In the right panel in figure 9, we show the resummed p T spectrum at NLL (0+L) (green dotted), NNLL (0+L) +NLO 0 (blue dashed), and N 3 LL (0+L) +NNLO 0 (red solid). The resummation including leptonic power corrections cures the unphysical behaviour of the fixed-order results, yielding a well-behaved spectrum in the full p T range, with a resummed Sudakov shoulder at p T ≈ m W /2. Note that the cross section beyond the edge is already populated at NLL (0+L) without any fixed-order matching. We stress that without including the exact leptonic tensor, the resummation would only affect the region p T < m W /2, and not cure the peak region. In fact, the results with strict LP resummation would look very similar to the pure fixed-order results, with the N 3 LL (0) +NNLO 0 essentially indistinguishable from the pure NNLO 0 result. This is the first time that resummed N 3 LL results for the p T spectrum are presented, and we observe extremely good perturbative convergence, with the results at NNLL (0+L) +NLO 0 and N 3 LL (0+L) +NNLO 0 falling on top of each other. We leave a more detailed phenomenological analysis of the p T spectrum to future work.

φ * spectrum
The φ * observable was first proposed in ref. [203], extending earlier work on the a T observable [204,205]. Both observables are sensitive to small q T , but promise better experimental resolution than q T itself due to being based on angular measurements, which are less susceptible to the momentum resolution of the individual lepton momenta than q T itself.
The factorization and resummation for a T was first studied in ref. [206] at NLL, and extended to both a T and φ * at NNLL+NLO in refs. [70,207,208], and also in refs. [209] including a study of nonperturbative contributions, and also more recently at N 3 LL+NNLO in refs. [81,83]. None of these calculations incorporate the finite recoil of the dilepton system in the calculation of the φ * observable, i.e. the employed definition of φ * is only an approximation of the actually measured observable, as discussed below. The resummation with the exact definition of φ * was considered in refs. [73,84] at NNLL+NNLO and NNLL+NLO via parton-level MC integration of the leptonic final state.
Ref. [203] defines the two closely related observables φ * CS and φ * η . We only consider the latter, as it is more commonly used in experiments. It is defined as where the acoplanarity angle is φ acop = π − ∆ϕ, with ∆ϕ being the azimuthal opening angle between the leptons in the lab frame, and cos θ * η = tanh where η 1,2 are the two lepton rapidities.

Origin of power corrections
Using the parametrization of the lepton momenta in the Collins-Soper frame as given in eq. (2.55), eq. (4.25) can be written as (φ * ) 2 = 8κ( sin ϕ sin θ) 2 (κ − 2 + α 2 ) 2 (κ + 2 − α 2 + 2) , κ 2 = ( 2 − α 2 ) 2 + 4 2 sin 2 θ sin 2 ϕ , α 2 = 1 + 2 cos 2 ϕ sin 2 θ . (4.27) Note that φ * is boost invariant and thus independent of Y , and can depend on q T only through the dimensionless ratio = q T /Q. From eq. (4.27), one easily finds the special values (4.28) The singularity arises from the case where both momenta are parallel to each other in the transverse plane, such that φ acop = π and eq. (4.25) becomes ill-defined. Numerically, we have also tested that φ * monotonically decreases with |sin θ|, such that φ * can be uniquely inverted on the intervals θ ∈ [0, π/2] and θ ∈ [π/2, π], with the two solutions trivially related by θ 1 = π − θ 2 . Expanding eq. (4.27) in 1, one obtains the commonly employed approximation φ * (0) = |sin ϕ| . (4.29) The monotonicity of φ * with |sin θ| implies that this is a lower bound to φ * , 11 but from eq. (4.28) it follows that this bound is only saturated for θ = π/2. It is thus natural to ask whether there is a better approximation of φ * . From eq. (4.27), it is easy to see that φ * only vanishes if either = 0 or sin ϕ = 0. Expanding eq. (4.27) in these limits, we find Note that the expansion in small in the second line can be recovered by reexpanding the small-sin ϕ limit, because in eq. (4.27) each term in either multiplies sin ϕ or is enhanced relative to sin ϕ. The first line in eq. (4.30) is ill-defined for |tan θ| < already at the leading O(sin 2 ϕ), while this singularity only appears at the second order in . Eq. (4. 30) suggests that the fundamentally small quantity to be power counted is |sin ϕ|, not itself. In particular, any given value in φ * can in principle receive contributions from arbitrarily large .
In figure 10, we compare the exact result for φ * (red solid) to its leading expansions in small sin ϕ (blue dashed) and small (green dotted). We fix ϕ = π/8 and show results for = 0.25 (left panel) and = 0.5 (right panel). The gray vertical line shows the breakdown of the small-sin ϕ approximation at cos θ = (1 + 2 ) −1/2 . Since φ * is fairly insensitive to cos θ in a rather large range of cos θ, the small-expansion is a fairly good approximation in that region. However, it quickly deteriorates for moderate to large cos θ. In contrast, the small-sin ϕ expansions follows the exact curve remarkably well, almost up to the point where it becomes ill-defined. At large cos θ, both approximations break down, as φ * is driven by the small-θ behavior φ * ∼ 1/θ 2 .
We find that the region |tan θ| < cannot be described by the expansion in small-sin ϕ, which breaks down, or by the expansion small-, which assigns an artificially small value φ * (0) in this region. However, this does not invalidate the LP description of φ * , as this region of phase space is suppressed as O( 2 ), Another important property of φ * is that it is not azimuthally symmetric at LP due to its explicit dependence on ϕ. To identify which W i contribute to φ * at LP, we evaluate eq. (2.60) with the approximate observable in eq. (4.29), which yields dσ dQ 2 dY dcos θ dφ * (0) = 3Q 16πE 2 cm i=−1,0,2,4 (4.32) Here, theW i are the hadronic structure functions in Fourier space. In eq. (4.32), we have already carried out the integral over ϕ, which gives rise to kernels K i , while we are still differential in cos θ. One can easily incorporate any LP fiducial cuts that depend on θ, but are independent of ϕ into eq. (4.32), which holds for most common cuts such as eq. (4.5). The nonvanishing kernels entering eq. (4.32) are given by where si(β) = − ∞ β dt sin t/t is the sine integral. The kernels for i = 1, 3, 5, 6, 7 vanish because the corresponding g i (θ, ϕ) are odd in ϕ or under ϕ → ϕ + π. Since g −1,0,4 (θ, ϕ) are independent of ϕ, they give rise to the same kernel K −1,0,4 . In contrast, W 2 is dressed with a different kernel K 2 due to the nontrivial ϕ dependence of g 2 (θ, ϕ) ∝ cos(2ϕ). In particular, for β → ∞ one has K 2 (β) ≈ − cos β = K −1,0,4 (β), and thus there is a relative phase shift of π.
Eq. (4.32) is convenient, as it effectively only reweights the (resummed) hadronic tensor in Fourier space with K i (b T Qφ * ), compared to J 0 (b T q T ) appearing in q T resummation. In momentum space, this is equivalent to the fact that the spectrum for the LP φ * (0) in eq. (4.29) can be obtained by reweighting the (resummed) q T distribution with the angle to the dilepton system. The convenient form of eq. (4.32) was first noticed in refs. [206,207], where it was also noted that the cos(b T Qφ * ) gives rise to a plateau in the resummed φ * spectrum, in contrast to the Sudakov peak encountered in q T resummation. However, the form in eq. (4.32) has the distinct disadvantage that it does not allow one any longer to include fiducial power corrections due to additional fiducial cuts beyond the strict LP. The above previous works did not consider the contribution from W (0) 2 , which involves the double-Boer-Mulders contribution, see eq. (2.103). Comparing to table 1, W 2 does not contribute to azimuthally symmetric observables at LP, and thus is not encountered in the LP q T resummation. In practice, we expect this contribution to be rather small, as h ⊥ 1 only matches onto subleading twist-3 collinear PDFs, and thus we will not consider it in our numerical study. Nevertheless, it is interesting to note that the φ * spectrum may give direct access to the double Boer-Mulders effect, and we leave a more detailed study of this for future work. At leading twist-2, W 2 is suppressed as O(λ 2 ), such that our resummed spectrum σ (0+L) still fully captures all linear power corrections arising from small q T with leading-twist collinear PDFs, and this also holds when including additional fiducial cuts.

Numerical results
We now turn to the numerical study of power corrections to φ * at one loop. For simplicity, we work at fixed Q = m Z , and normalize all results to the tree-level cross section. The results for inclusive and fiducial Drell-Yan are shown in the left and right panel of figure 11, respectively. In both plots, the exact results σ are shown by the red points, the factorized prediction including fiducial power corrections σ (0+L) by the blue line, and their difference by the green dot-dashed curve. This is contrasted by the gray-dashed curve which shows the difference σ − σ (0) between the exact and strict LP (i.e. employing φ * (0) ) result. The orange, dot-dashed curve shows the difference σ (0+L) GJ − σ (0+L) between our default CS tensor decomposition, and an alternative choice corresponding to the GJ frame.
In the inclusive case (left panel), σ (0+L) and σ (0) only differ by whether φ * is imple-mented exactly or using φ * (0) . In this case, we observe large linear corrections to σ (0) , whereas corrections to σ (0+L) appear to be quadratically suppressed. Interestingly, the σ (0+L) GJ seems to have linear corrections as can be seen from the linear scaling of the difference σ (0+L) GJ −σ (0+L) . Hence, σ (0+L) for a generic frame receives linear corrections, although they are roughly an order of magnitude suppressed compared to σ (0) .
For fiducial Drell-Yan production (right panel), we observe linear power corrections for both σ (0) and σ (0+L) , which are larger than the power corrections in the inclusive case, especially for σ (0) . Nevertheless, we again see see that σ (0+L) has significantly smaller corrections than σ (0) , despite having the same linear scaling in φ * . The ambiguity between the two choices of tensor decomposition is again a linear effect, but again at much smaller overall magnitude than the corrections beyond σ (0) .
Overall, we find that φ * generically receives linear power corrections. In addition to the common fiducial corrections, φ * is affected by corrections from expanding the observable itself, and by the fact that even very small φ * receives contributions from large q T , as is apparent from eq. (4.29). Hence, a priori there is no reason to expect that corrections to φ * are quadratically suppressed. Nevertheless, σ (0+L) includes all linear power corrections from small-q T , which is reflected by the corrections to σ (0+L) being significantly reduced compared to σ (0) . We also note that the choice of tensor decomposition strongly affects how well contributions from large q T are captured. Empirically, we find that our default choice corresponding to the CS frame minimizes the size of power corrections, but we have not been able to identify an underlying reason for this observation.
We conclude this section by studying the impact of the fiducial power corrections on the resummed φ * spectrum with fiducial cuts. In the top-left panel of figure 12, we show the difference between the strict LP and the exact φ * spectrum at NLO 0 (blue, short-dashed) and NNLO 0 (red, long-dashed). For reference, the gray line shows our best prediction σ * at N 3 LL (L) +NNLO 0 , scaled down to 5% of its original size. We note a large discrepancy between NLO 0 and NNLO 0 as φ * → 0, indicating large, unresummed logarithms in the power corrections and consequently poor perturbative convergence in this regime.
In the top-right panel of figure 12, we show the difference between the LP φ * spectrum and the resummed and matched φ * spectrum, at NLL (0+L) (green, dotted), NNLL (0+L) +NLO 0 (blue, dashed) and N 3 LL (0+L) +NNLO 0 (red, solid), which corresponds to the left panel but with the linear power corrections resummed. As a result, the divergence as φ * → 0 is cured, and we observe very good perturbative convergence between the different resummed predictions. Note that in contrast to q T , see figure 7, the resummed power corrections do not vanish as φ * → 0, due to the different weighting of the Sudakov factor with cos(b T Qφ * ) rather than J 0 (b T q T ), cf. eq. (4.32).
Finally, in the bottom panel of figure 12 we show the remaining fixed-order power corrections after including the fiducial power corrections in the resummation at NLO 0 (blue, short-dashed) and NNLO 0 (red solid), again including our best prediction σ * for reference. Since all terms diverging as φ * → 0 are included in σ (0+L) , the remaining power corrections are well-behaved as φ * → 0, which makes their overall size almost negligible. Our best prediction dσ * for the total spectrum at N 3 LL (0+L) +NNLO 0 is indicated as a light gray line for reference, scaled down to 5% of its original size. The bottom panel shows the remaining fixed-order power corrections, which are finite for φ * → 0 (note the difference in vertical scale).

Applications in fixed-order subtractions
The inclusion of fiducial power corrections in the q T factorization theorem by treating the leptonic tensor exactly can be immediately applied to improve the q T subtraction method. In fact, we have employed this in section 4.3 to obtain NNLO predictions for the p T spectrum. Here, we elaborate in more detail on this, illustrating that this approach can significantly and systematically improve the q T subtraction. In particular, it becomes crucial in regions of phase space sensitive to leptonic power corrections, such as p T close to the Jacobian peak, where the strict LP expansion breaks down, and thereby also the subtractions based on the strict LP limit.

q T subtraction with fiducial power corrections
The q T subtraction method has been originally proposed in ref. [115] for the calculation of color-singlet processes at NNLO, and its application at N 3 LO was discussed in refs. [177,210]. Here, we briefly review this method, and discuss how to extend it to include fiducial power corrections. Here, we only discuss the simplest implementation of q T subtractions as a phase-space slicing. As discussed in ref. [104], fiducial cuts can also be explicitly accounted for in a differential implementation of q T subtractions, which can also be further improved with the methods we discuss here, which we leave for future work.
We denote the cross section to be calculated as σ(X), where X summarizes all observables to be differential in, such as Q, Y , and p T , as well as possible fiducial cuts. The cross section can be written as where the cumulative cross section as a function of q cut T is defined as The q T subtraction is typically implemented as a slicing method by adding and subtracting a suitable subtraction term, For color-singlet processes, q T vanishes by construction in the Born limit. Hence, the integral in eq. (5.3) necessarily involves at least one resolved real emission, and can thus be calculated from the corresponding Born+1-parton process at one order lower than σ(X) itself. The cancellation of virtual and real divergences occurs only in the limit q T → 0. By constructing σ sub such that it fully describes this limit, the cancellation of IR divergences only occurs within σ sub , and ∆σ is a power correction that vanishes as q cut T → 0, and thus can be neglected for sufficiently small values of q cut T . To construct σ sub and study the size of ∆σ, it is useful to expand the differential cross section and its cumulant for q T Q and q cut T Q, respectively, where the different contributions scale as The dσ (0) /dq T diverges as 1/q T for q T → 0, and is in fact precisely given by the factorization theorem in eq. (2.103). Note that it contains δ and plus distributions to regulate this divergence, which precisely encodes the cancellation of virtual and real IR divergences. The dσ (m) /dq T with m > 0 are integrable as q T → 0, and consequently σ (m) (q cut T → 0) → 0. The subtraction term σ sub (X, q cut T ) must contain all singular terms, and has the form It follows that the correction term in eq. (5.3) scales as where m corresponds to the first term in the sum in eq. (5.4) that is not contained in σ sub . So far, q T subtractions have always been implemented in this fashion by choosing σ sub = σ (0) . As discussed in section 2, for inclusive processes the sum in eq. (5.4) starts with m = 2, i.e. ∆σ is suppressed as O(q 2 T /Q 2 ). These terms were explicitly calculated in ref. [160] for Higgs and Drell-Yan production differential in Q and Y , but inclusive in their decay. The corresponding results for inclusive Higgs and Drell-Yan production, i.e. integrated over Y , were also calculated later in refs. [139,211].
When applying fiducial cuts, one expects corrections to be only suppressed as O(q T /Q), i.e. m = 1, as shown in ref. [104] and discussed more generally in sections 2 and 4. 12 In section 4, we also demonstrated that more complicated observables such as p T or φ * can also have much larger corrections than expected.
In section 2, we showed that the linear power corrections in q T /Q are unique for any observables that are azimuthally symmetric in the q T → 0 limit, in which case they only arise from the leptonic tensor. Since the subtraction method is inherently tied to a calculation at leading-twist in collinear factorization, this also holds for more complicated observables such as φ * . This means that in many relevant cases the dσ (1) /dq T term in eq. (5.4) can be unambiguously predicted and included in the subtraction term. Hence, instead of eq. (5.6) we use where σ (0+L) includes all fiducial power corrections using our default CS tensor decomposition. In principle, other choices are possible as well. From our discussion of φ * in section 4.4, we have some indication that the CS decomposition tends to minimize the power corrections in more complicated cases, which motivates using it as our default choice also in the subtraction context.

Application to the lepton p T spectrum
We illustrate the advantage of including the leptonic corrections in eq. (5.8) for the case of the p T spectrum. As discussed in section 4.3 (see e.g. figure 8), we can distinguish three regions in the spectrum according to their sensitivity to small q T . The expected dependence on q cut T and the linear subtraction terms are as follows: 1. Sufficiently below the Jacobian peak, p T (m W − q cut T )/2, only quadratic power corrections arise, and thus we expect no significant impact from including the linear terms in the subtraction.
2. Around the Jacobian peak, |2m W − p T | q cut T , the standard power counting breaks down, and one is sensitive to the singular behavior as q T → 0. In this region, keeping fiducial power corrections becomes strictly necessary, since the strict LP expansion breaks down due to leptonic power corrections q T /(2m W − p T ). In particular, the region m W /2 < p T < (m W + q cut T )/2 requires to take these corrections into account, as otherwise the singular prediction is limited to p T < m W /2, and the region beyond the peak will contain leftover singular behavior.
3. Far above the Jacobian peak, p T (m W + q cut T )/2, there are no contributions from the subtraction term, and only the above-cut integral in eq. (5.1) is relevant.
We illustrate these effects numerically at NLO. Reference results are obtained by evaluating pp → W + → + ν at NLO using MCFM, which employs Catani-Seymour dipole subtractions [212,213]. The corresponding results using q T subtraction are obtained by combining the below-cut contribution dσ sub /dp T obtained from SCETlib, both for σ sub = σ (0) and σ sub = σ (0+L) , with the above-cut contribution obtained by evaluating pp → (W + → + ν ) + 1 parton at LO 1 using MCFM. We work in the narrow-width approximation such that Q = m W , which avoids smearing out the Jacobian peak, as is necessary to clearly see the effects described above. To keep numerical uncertainties from the Monte Carlo integration small, we bin the p T spectrum in bins of width 1 GeV, with p T = m W /2 being aligned with a bin boundary to cleanly separate below-peak and above-peak bins.
In figure 13, we show the relative difference for the p T spectrum obtained with q T subtraction to the exact NLO result, so smaller values imply that the q T subtractions work better with smaller missing power corrections. We show results for the two values q cut T = 0.4 GeV (left) and q cut T = 1.6 GeV (right). In the top row, we are inclusive in the lepton, while in the bottom row we apply the cuts p T ≥ 25 GeV, |η | < 2.4. The blue points use the strict LP subtractions, while the red points include the fiducial power corrections in the subtraction. The error bands indicate the MC integration uncertainties. For reference, we also show the exact NLO spectrum in gray in arbitrary units.
Due to the rather large bin size of 1 GeV, only the two bins adjacent to the Jacobian peak at p T = m W /2 are significantly affected by the choice of q cut T . Below the first bin left of the peak, power corrections are essentially negligible for all choices of q cut T , irrespective of whether or not the fiducial power corrections are included and independent of the cuts. Beyond the first bin right to the peak, there is no contribution from the subtraction term, and the spectrum is entirely given by the above-cut contribution. However, as expected, the bins adjacent to the Jacobian peak are significantly affected by the large leptonic power corrections. At strict LP (blue points), these are of O(80%) for q cut T = 0.4 GeV, and even of O(10) for q cut T = 1.6 GeV, indicating the breakdown of the subtractions. Including the fiducial power corrections (red points) automatically retains all leptonic power corrections and as a result the power corrections are of the order of a few percent, leading to stable predictions even very close to the peak. The results with and without addition fiducial cuts are very similar, as expected by the absence of linear fiducial power corrections when being inclusive in the neutrino. Figure 14 shows the q cut T dependence of the lepton p T spectrum, for the bin −1 ≤ p T − m W 2 ≤ − 1 2 GeV (left), and the bin − 1 2 ≤ p T − m W 2 ≤ 0 GeV (right), i.e. the two bins left to the Jacobian peak. (Note the finer bin width of 0.5 GeV, compared to 1 GeV above). The result for the first bin right of the peak is very similar to the right figure, while bins further away have negligible q cut T dependence that is entirely driven by MC integration uncertainties, and hence we do not explicitly show these. We also only show the inclusive case. In figure 14, we show the relative difference of the result obtained using q T subtractions to the MCFM result using Catani-Seymour subtractions. As before, the blue and red points correspond to implementing the q T subtractions without and with fiducial power corrections, and the error bars indicate MC integration uncertainties.
In the lower bin (left panel), the q cut T dependence of the q T subtraction including fiducial power corrections (red points) is consistent with zero within the MC uncertainties, illustrating the extremely good convergence of the subtractions. In contrast, at strict LP (blue points) there is a sizable difference to the exact result. While it is falling with q cut T , one needs a very small value q cut T ∼ 0.2 − 0.4 GeV to achieve an accuracy of O(10 −2 ). In the bin adjacent to the peak (right figure; note the different plot range), the difference between Catani-Seymour and strict-LP q T -subtraction results is consistently larger by a factor ∼ 10 2 , indicating the much stronger sensitivity to small q T ∼ p T − m W /2. As before, at strict LP (blue points) there is a strong dependence on q cut T , but even at the extremely small value of q cut T = 0.1 GeV the result has an O(1) uncertainty due to missing power corrections, showing the breakdown of the method. After including the fiducial power corrections (red points) the q cut T dependence is negligible within the MC integration uncertainties.
Note that there seems to be a systematic difference compared to the MCFM result obtained with Catani-Seymour subtractions. However, in this extreme region of phase space it is not a priori clear which of the results is more reliable and whether the reported MC uncertainties are still trustworthy. Thus, we refrain from a definite statement which subtraction method performs better in this scenario.
Our numerical checks confirm the theoretical observation that close to the peak, where one is sensitive to small q T , reliable predictions can only be obtained when taking the fiducial power corrections into account. In bins sufficiently far away from the peak, this is not as crucial, but can still significantly improve convergence of the method, i.e. it allows one to use much larger values of q cut T than the strict LP subtraction whenever linear fiducial corrections are present. This becomes particularly important at higher orders, where it becomes numerically challenging to push q cut T to such small values. In our analysis in section 4.3, we require the p T spectrum at NNLO 0 , which we obtain by combining the NLO 1 calculation of pp → W + +1 parton from MCFM with a q T subtraction including fiducial power corrections supplied by SCETlib, choosing q cut T = 0.8 GeV. We verify that this yields a sufficiently precise result in the first bin left of the Jacobian edge, −1 ≤ m W /2 ≤ 0 GeV. Compared to above we have increased the bin width from 0.5 GeV to 1 GeV to reduce numerical uncertainties. In the left panel in figure 15, we show the absolute cross section in this bin as a function of q cut T for σ (0+L) (red points) and σ (0) (blue points). Clearly, only the result using σ (0+L) shows a stable convergent behavior by giving a result insensitive to q cut T , whereas the result for σ (0) does not converge and is quite far from the true result. We then use a simple quadratic extrapolation of σ (0+L) (q cut T ) to q cut T → 0 as a proxy for the exact result in this bin. In the right figure, we show the relative difference of σ (0+L) (red points) and σ (0) (blue points) to this extrapolation. As expected, σ (0) shows O(1) deviations due to the clear lack of convergence towards small q cut T due to uncontrolled leptonic power corrections. In contrast, σ (0+L) is in good agreement with the extrapolation, with deviations of O(10 −3 ) that are largely driven by numerical uncertainties. This justifies our choice of q cut T = 0.8 GeV for the results presented in section 4.3. It also confirms the expectation that the systematic pattern of improvement from including fiducial power corrections in the subtractions observed at NLO 0 will carry over to NNLO 0 . We leave a more detailed investigation at NNLO 0 to future work.

Application to Drell-Yan with symmetric cuts
As a second application to illustrate the q T subtraction including fiducial power corrections, we consider Drell-Yan pp → Z/γ * → + − with fiducial cuts on the final-state leptons. In particular, we focus on enhanced power corrections that arise when one applies symmetric cuts, where both leptons are required to have p T ≥ p min T . For this case, it was already remarked in ref. [117] that the q T subtractions become very sensitive to q cut T leading to instabilities in the perturbative calculation. In fact, these instabilities have already been pointed out long ago in ref. [116], and traced back to a sensitivity to small transverse momenta, which may require resummation to obtain well-behaved results. In addition, as explained in ref. [104] and in section 4.2, the lepton cuts induce linear power corrections.
Here, we will illustrate in more detail that the enhanced linear power corrections can be avoided by employing q T subtraction with fiducial power corrections. Since this allows one to reach the precision of a LP subtraction at larger values of q cut T , it also reduces the region of phase space where fixed-order calculations need to be evaluated, and thus should also reduce the effect of the aforementioned instabilities. In the following we consider the total cross section for the Drell-Yan process, studying the effect of different p T cuts on the two leptons, namely p T h ≥ 25 GeV , 25.1 GeV , 25.5 GeV , 26 GeV , p T s ≥ 25 GeV . (5.9) Here, p T h (p T s ) is the harder (softer) transverse momentum. We only vary the cut on the harder lepton, while the cut on the softer lepton is fixed to p T s ≥ 25 GeV, and we always apply the same rapidity cut |η | ≤ 2.4 for both leptons. As in section 5.2, we consider the cross section at NLO, combining a one-loop below-cut calculation of pp → Z/γ * → + − from SCETlib, with a tree level above-cut calculation of pp → (Z/γ * → + − ) + j from MCFM, always working in the narrow-width approximation. We consider the following values of the q cut T subtraction cutoff: q cut T ∈ {0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4} GeV . (5.10) First, we study the q T spectrum, comparing the exact result with σ (0) and σ (0+L) . This is shown in figure 16, where the four panels correspond to the four lepton cuts in eq. (5.9), and the corresponding plot without fiducial cuts is shown in figure 6. In each panel, the red points show the full result σ and the blue line the LP prediction including fiducial power corrections σ (0+L) . In all cases, including the fiducial power corrections substantially improves the accuracy, with the leftover power corrections σ − σ (0+L) (dashed green) being quadratically suppressed and of similar size as in the inclusive case. On the other hand, using the strict LP limit, the leftover power corrections σ − σ (0) (dotted gray) are around an order of magnitude larger and only linearly suppressed. Varying the cut on the harder lepton has almost no effect on the green curves, but changes the gray curve, which has a dip from changing its sign, whose position is very sensitive to the difference of the two p T cuts.
To show the precise effect on the q T subtraction, we show in figure 17 the relative difference between the total cross section obtained via q T subtraction to the reference NLO result obtained from MCFM using Catani-Seymour subtractions as a function of q cut T . The blue points show the result from q T subtractions at strict LP, σ sub = σ (0) , while the red points show the result including fiducial power corrections in the subtractions, σ sub = σ (0+L) . The blue points are slightly shifted for better visibility. The uncertainties show the MC integration uncertainties.
For reference, the top-left panel shows the result without any fiducial cuts, in which case there is no difference between the two subtractions. The remaining power corrections are quadratically suppressed, and thus one observes a very good convergence, with deviations to the exact result of only 0.1% already for q cut all three cases, the power corrections to the subtraction with fiducial power corrections (red points) are of similar size as in the inclusive case. Even though the q cut T dependence becomes somewhat smaller in magnitude with increasing difference p min T h − p min T s , reflecting a sensitivity on q T ∼ |p min T h − p min T s |, this is mostly accounted for by including the fiducial power corrections, and as before choosing q cut T ∼ 3.2 GeV suffices for permille accuracy. In contrast, when implementing the q T subtraction at strict LP (blue points), the size of missing power corrections increases by one to two orders of magnitude, reflecting the enhanced power corrections from the fiducial cuts. To reach the 1 precision goal requires a substantially smaller q cut T ≈ 0.1 GeV. Furthermore, we observe again that the shape of the q cut T dependence of σ (0) (q cut T ) strongly depends on the precise difference |p min T h − p min T s | between the lepton cuts. In particular, for p T h ≥ 25.1 GeV (bottom left), we see that the q cut T dependence accidentally vanishes around q cut T = 0.8 GeV, while it is quite flat in the region q cut T ∼ 0.1 − 0.8 GeV. This region roughly corresponds to the the sign change visible in the corresponding q T spectrum, see figure 16 (top right). This behavior is very dangerous, as it can easily be confused with an apparent convergence if the q cut T dependence is not properly extrapolated. 10 -1 Figure 17. Relative difference between the exact NLO Drell-Yan cross section and its approximation with q T subtractions as a function of q cut T for different fiducial cuts, as indicated in the labels. The red and blue points correspond to implementing the q T subtraction with σ sub = σ (0) and σ sub = σ (0+L) , respectively. The blue points are slightly shifted for clarity, and the error bars indicate the MC integration uncertainties.
In summary, we again find that including the fiducial power corrections in the q T subtraction significantly improves the convergence in the q cut T → 0 limit, and thus makes the method much more reliable. In particular, for the case of symmetric lepton p T cuts, it recovers the same level of power corrections as for the inclusive case.

Comparison to data
In this section, we compare our N 3 LL+NNLO 0 resummed predictions for q T and φ * for Drell-Yan, pp → Z/γ * → + − ( = e, µ), with the following precision LHC measurements: We consider two sets of predictions: The strict LP resummation with fiducial power corrections only included via the fixed-order matching is denoted as N 3 LL (0) +NNLO 0 , and analogously at lower orders. The resummation including fiducial power corrections is denoted as N 3 LL (0+L) +NNLO 0 , and analogously at lower orders. In this case, the fixed-order matching only adds the remaining genuine (non-fiducial) power corrections.

q T spectrum
In figure 18, we compare our results to the CMS 13 TeV measurements for q T . The top panels show the q T spectrum at NLL (green), NNLL+NLO 0 (blue), and N 3 LL+NNLO 0 (orange). The bands show the estimated perturbative uncertainties as discussed in section 3.4. Note that the theory predictions are obtained for the measured binning, but are drawn as smooth curves for clearer visibility. The lower panels show the same results normalized to the respective highest-order prediction. The results using strict LP resummation are shown on the left (lighter shading), while those including the resummation of fiducial power corrections are shown on the right (darker shading). In both cases, we observe good convergence of the resummed predictions, with substantially reduced perturbative uncertainties at subsequent higher orders, as well as good agreement with the data. Nevertheless, resumming the fiducial power corrections on the right further improves the perturbative convergence and also yields a systematically better agreement with the data. The data agreement deteriorates in the first two bins, which can be attributed to small-q T nonperturbative effects. These are expected to become important for q T 2 GeV, but the nonperturbative ingredients necessary to account for these effects are not included in our predictions. This is also reflected in the substantially increased perturbative uncertainties in this region.
To further illustrate the importance and impact of the fiducial power corrections, in figure 19 we show the analog of the bottom panel of figure 18 but comparing to the pure resummed results, i.e., without including the fixed-order matching corrections. The strict LP resummation (left) completely fails to describe the data, showing that in this case the fixed-order matching corrections that supply the fiducial power corrections at fixed order are essential. On the other hand, upon resumming the fiducial power corrections (right), the excellent data agreement is restored even without the fixed-order matching. In other words, with the fiducial power corrections included in the resummation, the fixed-order matching becomes essentially unimportant for q T 40 GeV, both at NNLL and N 3 LL.
In figure 20, we show the analogous comparison for the ATLAS 8 TeV measurements [4] in the electron channel, with the top panel showing the q T spectrum itself, while the bottom panel shows the relative differences to the respective highest-order prediction. The ATLAS measurements are normalized to the total cross section. We correspondingly normalize our predictions to the integrated cross section at the corresponding order. 13 As before, we  Figure 18. Predictions for the Drell-Yan fiducial q T spectrum without (left) and with (right) resummed fiducial power corrections compared to CMS 13 TeV measurements [15]. The top panels show the spectrum, with the theory predictions drawn as smooth curves for better visibility. The bottom panels show the percent differences to the respective highest-order prediction central value.  Figure 19. Same as the bottom row of figure 18, but without including power corrections from the fixed-order matching.  Figure 20. Predictions for the normalized Drell-Yan fiducial q T spectrum without (left) and with (right) resummed fiducial power corrections compared to ATLAS 8 TeV measurements [15] in the e + e − channel. The top panels show the normalized spectrum. The bottom panels show the percent differences to the respective highest-order prediction central value. The analogous results for the µ + µ − channel can be seen in appendix B in figure 24. see good perturbative convergence of the predictions, as well as good agreement with the data. The data agreement again improves when resumming the fiducial power corrections on the right, leading to an overall flatter shape and reduced size in the difference between predictions and measurement. The results in the muon channel are practically identical, and are provided for completeness in figure 24 in appendix B.

φ * distribution
In figure 21, we compare our results for the normalized φ * spectrum to the ATLAS 8 TeV measurements [4] in the electron channel. The analogous results in the muon channel are provided in appendix B in figure 25. The top panel shows the predictions at NLL (green), NNLL+NLO 0 (blue), and N 3 LL+NNLO 0 (orange), with the bands showing the estimated perturbative uncertainties as discussed in section 3.4. The predictions are obtained with corrections. Since the uncertainties for the total cross section are much smaller than for the spectrum, they are practically irrelevant for this purpose. We have also checked that treating the µFO variation in a correlated fashion or treating it fully uncorrelated between spectrum and normalization and adding it in quadrature leads to essentially identical estimates of the total perturbative uncertainty.  Figure 21. Predictions for the normalized Drell-Yan fiducial φ * spectrum using the LP resummation (left) and including the resummation of fiducial power corrections (right) compared to ATLAS 8 TeV measurements [4] in the e + e − channel. The top panels show the spectrum, with the predictions drawn as smooth curves for better visibility. The bottom panels show the percent differences to the respective highest-order prediction central value. The analogous result for the µ + µ − channel can be seen in appendix B in figure 25. the experimental binning but are drawn as smooth curves for better visibility. The lower panels show the same results normalized to the respective highest-order predictions. The results using strict LP resummation for both the observable itself and the fiducial cuts are shown on the left (lighter shading), while those including the resummation of fiducial power corrections for observable and cuts are shown on the right (darker shading). In both cases we observe good convergence of the resummed predictions. For large φ * 0.5, the spectrum enters the fixed-order region and consequently the NLL (green) results start to deviate substantially, and to lesser extent also the NNLL+NLO 0 (blue) results. The first one or two bins are again sensitive to small-q T nonperturbative effects, which is reflected in their increased perturbative uncertainties. As for the q T spectrum, we find excellent agreement with the data, which is further improved on the right by resumming the fiducial power corrections, especially at NNLL+NLO 0 where the shape improves significantly.
In figure 22, we show the analogous comparison for the CMS 13 TeV φ * measurements [15]. The top panels show the spectrum itself, and the bottom panels the relative difference to the respective highest-order prediction. The predictions show the same be-  Figure 22. Predictions for the Drell-Yan fiducial φ * spectrum using the LP resummation (left) and including the resummation of fiducial power corrections (right) compared to CMS 13 TeV measurements [15]. The top panels show the spectrum, and the bottom panels show the percent differences to the respective highest-order prediction central value. haviour as at 8 TeV, and we again find good agreement with the data. Here, the improvements from resumming the fiducial power corrections are even more striking. While the strict LP resummation on the left shows a clear trend of overshooting the data at small φ * , we find nigh-perfect agreement across the spectrum with resummed fiducial power corrections. To further highlight this, in figure 23 we show the analog of the bottom panel of figure 22 but comparing to the pure resummed results only, i.e., without including fixedorder matching corrections. As we already saw for the q T spectrum, the LP resummation alone (left) basically fails to describe the data, showing that in this case the fixed-order matching is necessary in order to supply the fiducial power corrections at least at fixed order. In contrast, when resumming the fiducial power corrections (right), we find the same excellent data agreement as before up to φ * 0.5. This shows that the φ * spectrum has rather large sensitivity to power corrections throughout its spectrum and profits enormously from including them in the resummation. At the same time, the remaining fixed-order power corrections become negligible in this range. Beyond φ * 0.5, we enter the fixed-order region and as expected, the pure resummed results quickly deteriorate and matching to the fixed-order results becomes strictly necessary.

Conclusions
We have studied the impact of fiducial cuts and other generic leptonic measurements on the factorization of the Drell-Yan process at small transverse momentum q T Q. They generically induce fiducial power corrections in q T /Q relative to the well-studied leadingpower terms predicted by q T (equivalently TMD) factorization, which are significantly larger than the quadratic power corrections arising for the inclusive q T spectrum.
Using a Lorentz-covariant tensor decomposition of the leptonic and hadronic tensors combined with formal power-counting arguments in SCET, we have shown that for a large class of observables (those that are azimuthally symmetric at leading power), the fiducial power corrections are the only source of linear power corrections. Furthermore, by retaining the exact leptonic structure functions, the fiducial power corrections are unambiguously predicted from factorization and are correctly resummed to the same order as the leadingpower corrections.
We have also shown that the naive power expansion in q T /Q 1 can break down near the edge of Born phase space due to uncontrolled leptonic power corrections ∼ q T /p L , where p L is the distance from the edge of Born phase space. In such regions, it is strictly required to keep all leptonic power corrections ∼ q T /p L to correctly describe the actual leading-power limit. An important example is the p T spectrum near the Jacobian peak p T = Q/2 with p L = Q − 2p T . This provides another formal reason to keep the exact leptonic structure functions, because doing so guarantees that all required leptonic power corrections are automatically retained. The kinematic recoil prescriptions used in practical implementations usually yield an exact description of the leptonic decay and measurements. Our analysis shows for the first time that this is not only justified, but even necessary to obtain a description that is formally valid across the entire leptonic phase space. These conclusions also immediately apply to scalar processes such as Higgs production.
The tensor decomposition can be interpreted as a specific choice of vector-boson rest frame, which naturally emerged to be the Collins-Soper frame as defined by boosting from the lab frame, even when keeping nonzero masses of the initial state hadrons. The CS tensor decomposition yields nine Lorentz-scalar hadronic structure functions, which are defined for an arbitrary leptonic final state, and for Z/γ * → or W → ν decays directly map onto the commonly used angular coefficients for the cross section in the CS angles. We also discussed that Born leptons can be theoretically well defined in terms of an IRsafe Born projection of the full leptonic final state, including QED final-state radiation. We have shown that the cross section in the CS angles of the so-defined Born leptons admits a LO-like complete angular decomposition in terms of spherical harmonics, with the corresponding generalized angular coefficients modified by QED corrections.
We have presented resummed predictions with and without the resummation of fiducial power corrections at N 3 LL. The comparison of our predictions to precision Drell-Yan q T and φ * measurements from ATLAS and CMS confirms the importance of fiducial power corrections and their resummation in a striking way: While the strict LP resummation is able to describe the data within (theory) uncertainties, it fails to do so without including the sizeable fixed-order corrections. On the other hand, including the fiducial power corrections in the resummation systematically improves the agreement with the data, particularly at very small q T and φ * . Furthermore, the fixed-order matching corrections now become very small and essentially do not play a role below q T 40 GeV and φ * 0.5. Since the importance of unresummed, fixed-order power corrections is significantly reduced by resumming the fiducial power corrections, this has important implications for all phenomenological applications at small q T , for example the extraction of nonperturbative inputs to TMD distributions from data.
Furthermore, it is also much cheaper computationally to predict or obtain the fiducial power corrections via factorization instead of including (or extracting) them numerically from the full fixed-order calculations, since the latter become expensive quickly toward small q T or may not even be available. This is also reflected in a significantly improved performance for the q T subtraction method when the fiducial power corrections are included in the subtraction term predicted from factorization. In particular, this extends the applicability of the subtraction method to phase-space regions where it would otherwise break down due to uncontrolled leptonic power corrections. We have explicitly demonstrated this numerically for the p T spectrum near the Jacobian edge p T = m W /2. As a second example, we have studied the q T subtraction for Drell-Yan production with symmetric cuts on the final-state leptons, where we observed that enhanced power corrections from choosing symmetric cuts are fully accounted for by our formalism.
We have also considered the resummed p T spectrum, which plays a crucial role in the precision m W measurement at the LHC, and have obtained the first N 3 LL predictions for it. We have demonstrated that a reliable prediction of the physical p T spectrum near the Jacobian peak is possible and that it relies in an essential way on the interplay between small-q T resummation effects and the exact treatment of leptonic power corrections that describe the recoil of the leptonic system. The effective enhancement from leptonic power corrections also induces an enhanced sensitivity to additional hadronic structure functions.
We stress that these effects are suppressed in the q T spectrum itself, where they only enter indirectly through fiducial cuts. In contrast, they are directly exposed in the p T spectrum near the Jacobian peak. We therefore encourage precise, unfolded measurements of the p T spectrum for both W and Z in this region. Other interesting measurements would be the double-differential spectrum in the lepton momenta p T 1 , p T 2 , or equivalently p T 1 ± p T 2 , which can be expected to provide further sensitivity to these effects.
Here we used that the projected hard function can only depend on the Lorentz scalar q 2 = ω a ω b .
The case of an incoming antiquark in the a direction follows from a ↔ b, H µν V V q q (n a , n b ; ω a , ω b ) = H µν V V qq (n b , n a ; ω b , ω a ) .

(A.3)
This implies as expected from parity. In the inclusive case discussed in section 2.2, we used the shorthand Inserting the expression in eq. (2.114) into eq. (A.1), and suppressing the q 2 argument for brevity, we find where denotes the real part, the vector (v f ) and axial (a f ) couplings of a quark of flavor f to the Z boson were given in eq. (2.3). The terms denoted as O(α 4,6 s ) are proportional to the square of the singlet matching coefficients, which only start at C af = O(α 2 s ) and C vf = O(α 3 s ), whereas the nonsinglet coefficient C q = 1 + O(α s ). In the parity-odd case, we have H 4 ZZ qq = 16πα em N c δ qq 2v q a q |C q | 2 + 2 f v q a f C * q C af + a q v f C * q C vf + O(α 5 s ) , H 4 γγ qq = 0 , For W ± exchange, we have where V qq denotes the CKM-matrix element for q ∈ {u, c, t} and q ∈ {d, s, b} (and we take it to vanish in all other cases). The overall relative factor of 2 between H −1 and H 4 is due to the conventional normalization of g 4 (θ, ϕ) = cos θ rather than 2 cos θ in eq. (2.61). The renormalized matching coefficients C q (q 2 , µ) and C vf,af (q 2 , µ) can be extracted from the IR-finite parts of the qq vector and axial-vector form factors, which admit the same flavor decomposition as eq. (2.114). Explicit expressions for C q (q 2 , µ) in our notation can be found in ref. [182]. The two-loop results which enter in our analysis follow from the two-loop quark form factors [214][215][216][217]. In principle, there is an O(α 2 s ) contribution to the axial-vector singlet coefficient if the top quark is taken to be massive at the hard scale [147]. These contributions have however been found to be small at the level of the total cross section [17,123], and we neglect them in this paper.

B Additional comparisons to data
In section 6, we compared our predictions to the ATLAS 8 TeV measurements [4] in the pp → Z/γ * → e + e − channel. For completeness, here we provide the analogous results in the pp → Z/γ * → µ + µ − channel, which differ from the electron measurements by the lack of the lepton-rapidity exclusion region 1.37 < |η| < 1.52. In figure 24, we show the results for the q T spectrum in the muon channel, which corresponds to figure 20 in the electron channel. In figure 25, we show the results for the φ * spectrum in the muon channel, which corresponds to figure 21 in the electron channel. The results for the muon channel are essentially identical and confirm the conclusions drawn in section 6.