The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT subtraction method: electroweak corrections and power suppressed contributions

Building upon the formulation of transverse-momentum resummation for heavy-quark hadroproduction, we present the first application of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT subtraction formalism to the computation of electroweak corrections to massive lepton pairs through the Drell–Yan mechanism. We then study the power suppressed contributions to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT subtraction formula in the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\mathrm{cut}}$$\end{document}rcut, defined as the minimum transverse momentum of the lepton pair normalised to its invariant mass. We analytically compute the leading power correction from initial and final-state radiation to the inclusive cross section. In the case of initial-state radiation the power correction is quadratic in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\mathrm{cut}}$$\end{document}rcut and our analytic result is consistent with results previously obtained in the literature. Final-state radiation produces linear contributions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\mathrm{cut}}$$\end{document}rcut that may challenge the efficiency of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_T$$\end{document}qT subtraction procedure. We explicitly compute the linear power correction in the case of the inclusive cross section and we discuss the extension of our calculation to differential distributions.


Introduction
The q T subtraction formalism [1] is a method to handle and cancel the IR divergences appearing in QCD computations at next-to-next-to-leading order (NNLO) and beyond. In its original formulation it has been applied to carry out a variety of NNLO QCD computations for the production of colourless final states in hadronic collisions [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Most of the above computations are now publicly available in Matrix [17]. A first application of q T subtraction to the computation of the approximate next-to-next-to-next-to-leading order (N 3 LO) QCD corrections to Higgs boson production through gluon fusion has been presented recently [18]. In the last few years, thanks to the formulation of transverse-momentum resummation for heavy-quark production [19][20][21][22][23] the method has been extended and applied to the production of top-quark pairs [24][25][26]. The q T subtraction counterterm is constructed by exploiting the universal behavior of the associated transverse-momentum (q T )d i stribution. Therefore, the subtraction is intrinsically non local and in practice the computation is carried out by introducing a cut, r cut on the tranverse momentum of the colourless system normalised to its invariant mass. When evaluated at finite r cut both the contribution of the real emission and the one of the counterterm exibit logarithmically divergent terms plus additional power suppressed contributions that vanish as r cut → 0. The efficiency of the subtraction procedure crucially depends on the size of such power suppressed contributions.
In the inclusive production of a colourless final state the power suppressed contributions are known to be quadratic in r cut (modulo logarithmic enhancements) [12]. This allows us to obtain precise predictions by either evaluating the cross section at sufficiently small r cut , or carrying out the r cut → 0 extrapolation 1 [17]. The power suppressed contributions to the next-to-leading order (NLO) total cross section have been explicitly evaluated in Refs. [28,29]. In the case of heavyquark production the r cut dependence is found to be linear [25,26,30], and it is an interesting question to investigate the origin of this peculiar behavior.
Up to now the q T subtraction formalism has been applied only to higher-order QCD computations. The formulation of the method for heavy-quark production can be straightforwardly extended to the computation of NLO electroweak (EW) corrections to the Drell-Yan process. The purpose of the present paper is twofold. We first present and discuss the first application of the q T subtraction formalism to the computation of the NLO EW corrections to the production of massive lepton pairs. Our computation consistently includes initial-state radiation, final-state radiation from the massive leptons and their interference, and our results are compared with the ones obtained with the well established public generator Sanc [31] and to an independent computation that we carry out with the dipole subtraction formalism [32]. Then, we present the analytic computation of the power suppressed contributions, and we confirm the linear r cut behaviour by computing its NLO coefficient. We also extend our results to the case in which cuts are applied.
The paper is organised as follows. In Sect. 2 we review the q T -subtraction formalism, by detailing its implementation up to NLO in the case of heavy-quark production. In Sect. 3 we apply the formalism to the computation of NLO EW corrections to the Drell-Yan process. In Sect. 4 we study the power suppressed contributions and we explicitly compute the leading power corrections in the case of final-state and initial-state radiation. In Sect. 5 we summarise our results.

2T h eq T subtraction formalism
The q T subtraction formalism [1] is a method to handle and cancel the IR divergences appearing in higher-order QCD computations. The method uses IR subtraction counterterms that are constructed by considering and explicitly computing the transverse-momentum q T distribution of the produced final-state system. At Born level such distribution is proportional to δ(q 2 T ). At higher perturbative orders multiple radiation of soft and collinear partons makes the distribution divergent in the q T → 0 limit. If the produced final-state system is composed of non-QCD (colourless) partons (e.g., leptons, vector bosons or Higgs bosons), the small-q T behaviour has a universal (process-independent) structure that is explicitly known up to the NNLO level (and, in part, at N 3 LO [18,33]) through the formalism of transverse-momentum resummation [34]. These results on transverse-momentum resummation are sufficient to fully specify the q T subtraction formalism for this entire class of processes. By using the formulation of transverse-momentum resummation for heavy-quark production [19][20][21][22][23], the q T subtraction formalism has been recently extended to this class of processes [24][25][26].
According to the q T subtraction method [1], the parton level differential cross section dσ QQ (N )NLO for the inclusive production process pp → QQ + X can be written as where dσ QQ+jet (N )LO is the QQ+jet cross section at (N)LO accuracy. The square bracket term of Eq. (1) is IR finite in the limit q T → 0, but its individual contributions, dσ QQ+jet (N )LO and dσ QQ, CT (N )NLO , are separately divergent. The IR subtraction counterterm dσ QQ, CT (N )NLO is obtained from the (N)NLO perturbative expansion (see, e.g., Refs. [24,35]) of the resummation formula of the logarithmically-enhanced contributions to the q T distribution of the QQ pair [19][20][21]: the explicit form of dσ QQ, CT (N )NLO can be completely worked out up to NNLO accuracy.
In the following we will limit ourselves to consider Eq. (1) up to NLO accuracy. The explicit expression of dσ QQ, CT NLO in the partonic channel ab → QQ + X reads [24] where M is the invariant mass of the QQ pair and the symbol ⊗ denotes convolutions with respect to the longitudinalmomentum fractions z 1 and z 2 of the colliding partons. The functions Σ (1) cc←ab in Eq. (2) can be written as where r = q T /M, and the coefficients Σ (1,k) The coefficients A (1) c and B (1) c are the first-order coefficients for transverse-momentum resummation ( A The functions P (1) ab (z) are the lowest-order DGLAP kernels (the overall normalisation is specified according to the notation in Eq. (41) of Ref. [35]). The functionsĨ k (r ) (k = 1, 2), which appear in Eq. (3), encapsulate the singular behavior at small q T , and they read (see Appendix B of Ref. [35]) where b 0 = 2e −γ E and we have introduced the modified Bessel function of imaginary argument The coefficient Σ (1,2) cc←ab (z 1 , z 2 ) in Eq. (4) controls the leading logarithmic contribution at small q T , while the coefficient Σ (1,1) cc←ab (z 1 , z 2 ) in Eq. (5) controls the next-to-leading logarithmic term. The latter has a first term (first line in Eq. 5) which is identical to what we have in the case of the production of a colour singlet. The second term (second line in Eq. 5) is due to soft radiation and it is an additional term that is specific of the q T subtraction method for the case of heavy-quark pair production [24]. Here Ŵ (1) t is the first-order contribution to the soft anomalous dimension for transverse-momentum resummation in heavy-quark production and its explicit expression is given in Eq. (33) of Ref. [21]. The soft anomalous dimension is a matrix acting on the colour indeces of the four hard partons in the Born level scattering amplitude |M cc→QQ . At this perturbative order the soft anomalous dimension is expressed in terms of colour correlators T i · T j with definite kinematic dependence, where the indices i and j refer to the hard-scattering partons.
The first-order hard-collinear coefficients H QQ NLO in Eq. (1) are also completely known [19][20][21]. In the next Section we apply the method to the computation of EW corrections to the Drell-Yan process.

NLO EW corrections to the Drell-Yan process
We consider the hadroproduction of a dilepton pair through the Drell-Yan mechanism. NLO EW corrections to this process have been considered in Refs. [36][37][38][39]. A tuned comparison of various Monte Carlo codes is presented in Ref. [40].
The NLO q T subtraction formalism for heavy-quark production reviewed in Sect. 2 can be straightforwardly extended to the computation of the NLO EW corrections to the Drell-Yan process. In this case the heavy-quark pair is replaced by a massive lepton pair and the abelian limit is carried out along the lines of Ref. [41]. The partonic cross section up to NLO EW can be evaluated by using where dσ R is the real emission cross section and the functions H NLO and dσ CT are obtained from the corresponding functions appearing in Eq. (1) with the replacements As is well known, at LO (i.e. O(α 2 )) both the qq and the γγ partonic channels contribute and we can write for the hadronic cross section where σ qq LO and σ γγ LO are the Born level cross sections in the qq and γγ channels, respectively. At NLO EW we can write where we have introduced the O(α 3 ) correction in the qq channel, ∆σ qq , the corresponding correction in the q(q)γ channel, ∆σ qγ , and the correction in the γγ channel, ∆σ γγ . Since the γγ channel provides only a very small contribution to the Drell-Yan cross section, ∆σ γγ will be neglected in the following discussion.
Our calculation is carried out by using an extension of the numerical program of Ref. [2]. All the required tree level matrix elements are computed analytically while the virtual EW corrections for qq → l + l − , which include vertex and box diagrams, are obtained by using Gosam [42,43]. We use the setup of Ref. [44], and, in particular, we work in the G µ scheme with and use the complex-mass scheme [45] throughout. More precisely, real and virtual photons emissions are controlled by α(0), while the α 2 in the LO cross section is derived from G F , m Z and m W . Following Ref. [44], the MRST2004qed [46] parton distribution functions (PDFs) are used setting the factorization scale to µ F = M Z . The following set of cuts are applied To validate our implementation, we have performed a tuned comparison with the public generator Sanc setting the mass of the heavy lepton to the muon mass m l ≡ m µ = 105.658369 MeV. In Table 1 we report our result for the lowest order cross sections σ qq LO and σ γγ LO , and the NLO EW corrections in the qq and qγ channels, ∆σ qq and ∆σ qγ .T h e NLO correction ∆σ qq is obtained performing the calculation at different values of r cut and extrapolating to r cut → 0 Table 1 Comparison of NLO EW corrections to Drell-Yan dimuon production (m l ≡ m µ = 105.658369 MeV) computed with q T subtraction and with the Sanc generator. In the qq channel the q T result is obtained with a linear extrapolation in the r cut → 0 limit (see Fig. 2), while in the q(q)γ channel it is obtained at r cut = 0.01%. The LO result in the qq and γγ channels is also reported for reference. The small discrepancy in the σ qq LO cross section, below 0.5 per mille, can be ascribed to the use of the complex mass scheme in our computation through a linear fit. This procedure will be fully motivated by an analytic computation presented in the next section. We find that the agreement is quite good, at few per mille on the NLO correction. The small discrepancy in the σ qq LO cross section, below 0.5 per mille, is consistent with the use of different schemes to handle unstable particles. In particular, at variance with what is done in Sanc, complex couplings are used within our computations.
The comparison with Sanc can be extended to differential distributions. In Fig. 1 we compare the invariant mass (m ll ) and rapidity (y ll ) distribution of the dilepton pair, the transverse momentum and rapidity of the positively charged lepton ( p T,l + and y l + ) obtained with our implementation with the corresponding results obtained with Sanc. For ease of reference, we plot also the LO result in black. We see that the agreement is good, the differences typically being well below the percent level. In addition, we have repeated our calculation by using the dipole subtraction method 2 [32] and the independent matrix-element generator Recola [48,49] for the virtual corrections, finding the same level of agreement as reported in Table 1.

Power corrections
In this section we are going to analytically study the power suppressed contributions to the q T subtraction formula in Eq. (8). In particular, we will show that the leading power suppressed contributions have a soft origin. In order to avoid complications due to the small lepton mass, which are sensitive to collinear logarithmic contributions and could obscure the r cut behavior, in the following we will consider a pair of heavy leptons of mass m l = 10 GeV. By using the setup of Sect. 3 (see Eqs. 12, 13), we have studied the dependence of the NLO corrections for the fiducial cross section on r cut . We have varied r cut in the range 0.01% ≤ r cut ≤ 1% and we have used the r cut -independent cross section computed with our inhouse implementation of the dipole subtraction method as normalisation. The results for the r cut dependent correction δ q T = ∆σ/σ qq LO in the qq and qγ channels are shown in Fig. 2. A distictive linear behavior in the dominant qq-annihilation channel emerges. Nonetheless, as reported in Ref. [17], it is known that symmetric cuts on the transverse momenta of the final state leptons challenge the convergence of q T -subtraction leading to a stronger dependence on r cut even in the case in which a charge-neutral final state is produced. In Fig. 3 we show the dependence of the NLO corrections for the inclusive cross section on r cut when no cuts are applied. Again a distinct linear behavior in the dominant qq-annihilation channel emerges, in agreement with what has already been observed for the case of the tt cross section [25], which can be clearly interpreted as a genuine new effect due to the emission of radiation off the massive final-state leptons.
In the following we analytically study the behavior of NLO cross sections computed with q T subtraction in the r cut → 0 limit. We are interested in determining the structure of the leading power correction to the inclusive cross section, and in identifying the origin of the linear behavior observed in Fig. 3. We recall that when applying the q T subtraction formula the second term on the right hand side of Eq. (8)is computed by introducing a lower limit r cut on the q T /M ratio. With such a cutoff we can treat separately the real contribution dσ R and the counterterm dσ CT . We start our discussion from the contribution of the counterterm. From Eq. (2)w e have The NLO coefficient Σ (1) cc←ab depends on r = q T /M only through the functionsĨ i (r ). Therefore we have In the small r limit the integralsĨ 1 (r ) andĨ 2 (r ) read , i.e., they depend quadratically on r modulo logarithmic terms. This results holds also at NNLO and beyond. It follows that the leading power corrections from the counterterm are always quadratic in r cut , independently on the perturbative order. As a consequence, the linear behavior with r cut that we observe in heavy-quark production and in the EW corrections to dilepton production must be due to the real emission.
In the following we analytically compute the real-emission contribution at small values of r cut .
We consider the production of a massive lepton pair in pure QED in the diagonal channel with p 2 3 = p 2 4 = m 2 . We define the variables Since there is a lower limit on the ratio r = q T /M we can safely work in d = 4 dimensions. The differential cross reads and the angular integral is defined in the centre-of-mass frame of the final-state leptons. By integrating Eq. (20) over q 2 T and M 2 and keeping into account the phase space constraints we obtain dσ qq dr 2 cut =− 1 32 where z min = 4m 2 s z max = 1 − 2r cut 1 + r 2 cut + 2r 2 cut .
The matrix element squared |M| 2 can be divided into three separate gauge invariant contributions: final state radiation, initial state radiation and interference. The interference contribution is odd under the exchange p 3 ↔ p 4 and therefore vanishes after angular integration. We now discuss the finaland initial-state contributions in turn.

Final-state radiation
The integration of the matrix element squared corresponding to final state radiation over the angular variables can be carried out along the lines of Ref. [50]. After partial fractioning, the required angular integrals have the form in terms of two coefficient functions, K 1 and K 2 , which are regular at z = 1 (soft limit) and do not depend on the cut-off parameter r cut : and In the small-r cut limit the integral in Eq. (24) can be computed by using the expansions and we obtain for the r cut dependence of the partonic cross section ≡σ FS LP (s; r cut ) +σ FS NLP (s; r cut ) + O(r 2 cut ) (28) where we have dropped terms which do not depend on r cut and we have introduced the Born cross section with β = 1 − 4m 2 s . Equation (28) shows that the final-state contribution to the NLO cross section, integrated down to r cut , contains the expected single logarithmic term in r cut , which is due to soft emission and will be cancelled by the subtraction counterterm [more precisely, by the abelian limit of the term in the second line in Eq. (5)]. The next-to-leading power contributionσ FS NLP (s; r cut ) is linear in r cut and it is responsible for the behavior observed in Fig. 3.

Initial-state radiation
The integration of the matrix element squared corresponding to initial-state radiation over the angular variables is straightforward and we obtain where the coefficient functions K 3 and K 4 now read The coefficient function K 3 (z; z min ) controls the most singular term, and is proportional to the Altarelli-Parisi splitting function. In order to evaluate the integral in Eq. (30)wehave to expand the distribution in the small r cut limit. Since we know that linear terms in r cut are absent, we have to expand up to O(r 2 cut ). Such expansion is tricky, because the functions K 3 (z; z min ) and K 4 (z; z min ) contain a square root which vanishes at z = z min . At variance with the case of final-state radiation, this leads to spurious singularities when the distributions appearing in the expansion involve derivatives at z = z min . We found it convenient to split the integration over z as follows The integral from z min to a can be safely computed by expanding in r cut . The integral from a to z max can be computed by using the expansion where we have defined the distributions δ (n) (z − b), 1 1−z a , D (1) (z, a) and D (2) (z, a) through their action on a test function f (z) as By combining the two contributions z min < z < a and a < z < z max the dependence on a cancels out and we obtain for the r cut dependence of the partonic cross section where we have dropped terms which do not depend on r cut and the dots stand for terms that vanish faster than r 2 cut as r cut → 0. At variance with Eq. (28), Eq. (39) contains a double and a single logarithmic term in r cut , which will be cancelled by the subtraction counterterm. As expected, the next-to-leading power contributionσ IS NLP (s; r cut ) is quadratic in r cut , modulo logarithmic enhancements.
As a byproduct of our calculation, we can reobtain the power suppressed terms in the case of the production of a vector boson of mass M. To get rid of the decay it is enough to carry out the m → 0 limit, while the constraint on the mass of the vector boson M eliminates the integration over the z variable. By using the expansion in Eq. (34) with a = 0the contributions from the functions K 3 (z; z min ) and K 4 (z; z min ) agree with the results in Eq. (4.7) and (4.8) in Ref. [29].

Numerical validation
In order to check the results presented in Sects. 4.1, 4.2 we have numerically implemented the exact real emission contribution to the cross section and the expansions in Eqs. (28) and (39). Fig. 4 Subtracted partonic cross section for final-state radiation (left panel) and initial-state radiation (right panel). The solid lines represent the subtraction of the leading-power term, while the red solid line is obtained by subtracting also the next-to-leading power terms in Eqs. (28) and (39), respectively. The upper panels show the result normalised to the Born cross section, while the lower panels show the result normalised to the r cut → 0 limit. The computation is carried out at fixed β = 0.6 In Fig. 4 we report the exact real emission partonic cross section in the qq channel for β = 0.6 as a function of r cut from which we have subtracted the leading-power contribution (black curve) and both the leading and next-to-leading power contributions (red curve). The numerical computation is separately carried out for the final-state radiation (left panel) and initial-state radiation (right panel) contributions. Both for final-state radiation and initial-state radiation the leading-power contribution exactly matches the divergent behavior of the real emission cross section which is finite in the small-r cut limit. The subtraction of the leading-power contribution exactly corresponds (up to quadratic terms in r cut ,s e eE q .16) to the second term on the right hand side of Eq. (8) and it is thus what is usually done in the standard q T subtraction procedure. In the case of final-state radiation (left panel) the subtracted cross section exibits the expected linear behavior, while for initial-state radiation (right panel) the subtracted cross section scales quadratically with r cut . When besides the leading-power contribution, also the nextto-leading power (linear) term is subtracted the final-state subtracted cross section (red curve) behaves quadratically with r cut , consistently with the result in Eq. (28). In the case of initial-state radiation, the additional subtraction of the next-to-leading power (quadratic) term makes the subtracted cross section almost independent on r cut .

Hadronic cross section
We now briefly comment upon the behavior of the hadronic cross section. In the following, we will show that when the fully inclusive cross section is considered, the convolution with the PDFs potentially introduces an additional linear term in r cut . In the case of final-state radiation such contribution could modify the parton level result. In the case of initial-state radiation such contribution could potentially change the power counting, by making the power correction linear. However, we will argue that, both for final-state and initial-state radiation, such additional term vanishes.
The real contribution to the hadronic cross section reads where S is the hadronic CM-energy. The presence of a finite r cut implies that where z max , defined in Eq. (22), behaves linearly with r cut The hadronic cross section in Eq. (40) can be rewritten as 1 0 dx 1 appearance of a further linear term through integration. We thus conclude that, as anticipated, in the case of final-state radiation the linear term in r cut is completely driven by the parton level result, while for initial-state radiation the convolution with PDFs will not produce linear terms in r cut .

Final-state radiation at next-to-leading power: beyond inclusive observables
In Sect. 4.1, we have established by means of an analytical computation that in the case of final-state radiation off massive emitters the leading power corrections are linear in r cut and we have explicitly evaluated the coefficient of the linear term. The feasibility of our computation and the methods involved crucially rely on the fact we consider a sufficiently inclusive observable as the total cross section. In the following we propose a way to promote the calculation of the final-state radiation contribution to differential level. Our starting point is the expansion of the real contribution to the differential cross section in the soft limit. According to the discussion in Appendix A, the leading soft contribution allows us to obtain the leading-power term in r cut ,whilethe next-to-soft contribution will allow us to obtain the next-toleading power. By inspection of Eq. (27) which involves the derivative of the δ-distribution, we indeed expect that higher order terms in the soft expansion contribute to the next-toleading power.
In the following we propose a strategy to numerically prove the above result, which in turns provides a procedure to compute the next-to-leading power in a fully differential way. From the soft contributions we can construct a local counterterm which cancels the singularities of the real cross section but does not contribute to the next-to-leading power. Then the subtracted cross section is finite in four dimensions and can be integrated numerically in the unresolved region r < r cut . By construction, since the standard q T subtraction counterterm does not lead to linear contributions in r cut (see Eq. 16), the combination of the standard q T subtraction formula for r > r cut with such an additional subtracted term for r < r cut will be free of linear terms in r cut .
To construct the additional counterterm we need a mapping which reabsorbs the radiation into a Born-like configuration. Among the available mappings at NLO, we choose the mapping proposed in Ref. [51]. It is a massive FKS [52] mapping dedicated to the case of the radiative emission off final state massive emitters and present some peculiar features: -the radiation is reabsorbed in such a way not to modify the partonic CM energy; -the energy of the radiation (in the CM frame), which controls the way the soft limit is approached, appears explicitly among the variables of the mapping.
In the above mapping, we identify an emitter and a radiated parton. The radiation variables are given by the radiation energy fraction ξ = 2E rad / √ s (s is the partonic CM energy), the cosine of the angle between the emitter and the radiated partons y and an azimuthal angle φ (we refer to Ref. [51]for more details). The phase space reads dΦ R = dΦ B × J (ξ, y,φ)dξ dydφ (46) where dΦ B is the Born phase space element and the jacobian J is given in Eq. (49) of Ref. [51], and reduces to where β = 1 − 4m 2 /s and y phy is the cosine of the physical angle between the emitted photon and the leptons in the Born configuration (in practice we have either y phy = y or y phy = −y [51]). The next-to-leading power correction as function of r cut is obtained by subtracting the new soft counterterm from the real emission contribution in the unresolved region r < r cut where in the argument of the second theta function, we take the soft limit M(φ R ) → M(Φ B ) = √ s. The expression in Eq. (48) is fully differential, so that it can be used also in the case in which cuts on the final state are applied. The contribution in Eq. (48) can be combined with the standard q T subtraction formula in Eq. (8) to obtain an improved subtraction procedure. In Fig. 5 we study the r cut dependence of the NLO EW correction to the complete Drell-Yan process when the q T subtraction formula is supplemented with the next-to-leading power term in Eq. (48). We consider pp collisions at √ S = 7 TeV and we compute the r cut dependent correction δ q T (r cut ) in the case in which no cuts are applied (Fig. 5 (left)) and when asymmet-ric cuts on the transverse momenta and rapidities are applied: p T,l − > 25 GeV, p T,l + > 20 GeV and |y l | < 2.5 ( Fig. 5  (right)). We see that in both cases the linear dependence with r cut is nicely cancelled. 3 The crucial point for this additional subtraction to be effective is that the additional counterterm in Eq. (47) scales like dξ/ξ, thereby leading to purely logarithmic contributions in r cut . We have checked that alternative local subtractions which do not fulfill this property do not lead to a cancellation of the linear term.
We conclude this section with few comments on the above results. The subtraction of the linear r cut behavior through Eq. (48) does not require any analytic integration. It just requires an appropriate phase space mapping. The reader may of course argue that there is no need to introduce the modification of Eq. (48) to achieve a smooth cancellation of the soft singularity. Indeed, at NLO one can simply use a local subtraction scheme like FKS or dipole subtraction to carry out the fully differential computation. Nonetheless, the strategy adopted here could prove itself useful when extending the computation to the mixed QCD-EW corrections with the q T subtraction formalism. In this case, given that we aim at the computation of an effect of the order of few per mille,having a quadratic instead of linear r cut behavior could dramatically improve the numerical control of the O(αα S ) contribution.

Summary
In this paper we have considered an application of the q T subtraction formalism to the production of massive lepton pairs in hadron collisions, and we have used this process as a case study to investigate the power suppressed contributions in the parameter r cut . We have shown that q T subtraction can be applied to evaluate NLO EW corrections to this process through a straightforward abelianisation procedure from heavy-quark production in QCD. The computation can be carried out for lepton masses as small as the one of the muon without substantial complications, and we have been able to successfully reproduce inclusive and differential results obtained with the numerical program Sanc. Our calculation paves the way to possible applications to the computation of mixed QCD-EW corrections [41,[53][54][55][56] and to NNLO QED corrections [41] to the Drell-Yan process.
We have then studied the power suppressed contributions to the q T subtraction formula in the parameter r cut .A si s known, in the case of the inclusive cross section, initial state radiation leads to power corrections quadratic in r cut , and we have explicitly evaluated the corresponding NLO coefficient. Generally speaking, linear power suppressed terms arise when cuts on the lepton transverse momenta are applied. We have shown that, even in the case of the inclusive cross section, final state radiation leads to a linear power correction in r cut . We have explicitly computed the coefficient of such linear term in the case of the inclusive cross section. By exploiting the purely soft nature of final state singularities, we have also discussed how the result can be extended to differential distributions. The method used to carry out such extension could be used in future applications of the q T subtraction formalism at O(αα S ).

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The research presented in this manuscript is of purely theoretical nature and does not have associated data which need to be deposited.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Soft Power Counting
In this Appendix, we discuss the soft power counting for final-state radiation. In the strictly soft limit, the phase space of the emitted photon with momentum k exactly factorizes (A.1) The leading power constribution to final-state radiation is given by the soft-factorisation formula where The power counting is more easily understood if we consider light-cone coordinates Then, the 1-body phase space volume has the form with k − = k 2 ⊥ /2k + . Considering for example the contribution from S 34 , the leading power unconstrained soft integral is given by In the above formula, we have enforced the soft kinematic with two back-to-back massive leptons. The azimuthal average is straightforward, after disentangling the product occuring in the denominator by means of the partial fractioning relation 1 (a 3 − cos θ)(a 4 + cos θ) = 1 a 3 + a 4 × 1 a 3 − cos θ + 1 a 4 +cos θ , (A. 8) and it gives To make the scaling with the transverse momentum manifest, we apply the following change of variables at fixed k ⊥ : (A.10) The soft integral becomes where s is the partonic CM energy. We can complete the calculation of the leading power contribution by performing the integration over the x variable. The relevant integrals are of the form which exactly matches the coefficient of the leading logarithmic divergence proportional to the charge product e 3 e 4 =−1 in Eq. (28). The contributions from I soft 33 and I soft 44 can be obtained in a similar way and reproduce the remaining term in Eq. (28). The power counting for the linear power correction follows now easily observing that the energy of the radiation scales with the transverse momentum 15) This implies that corrections to the soft approximation will produce linear terms in r cut .