Higher-order power corrections in a transverse-momentum cut for colour-singlet production at NLO

We consider the production of a colourless system at next-to-leading order in the strong coupling constant αS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{S}}$$\end{document}. We impose a transverse-momentum cutoff, qTcut\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{T}}^{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{cut}}$$\end{document}, on the colourless final state and we compute the power corrections for the inclusive cross section in the cutoff, up to the fourth power. The study of the dependence of the cross section on qTcut\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{T}}^{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{cut}}$$\end{document} allows for an understanding of its behaviour at the boundaries of the phase space, giving hints on the structure at all orders in αS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{S}}$$\end{document} and on the identification of universal patterns. The knowledge of such power corrections is also a required ingredient in order to reduce the dependence on the transverse-momentum cutoff of the QCD cross sections at higher orders, when the qT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{\mathrm{T}}$$\end{document}-subtraction method is applied. We present analytic results for both Drell–Yan vector boson and Higgs boson production in gluon fusion and we illustrate a process-independent procedure for the calculation of the all-order power corrections in the cutoff. In order to show the impact of the power-correction terms, we present selected numerical results and discuss how the residual dependence on qTcut\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{T}}^{{\displaystyle } {\scriptstyle } {\scriptscriptstyle } {\scriptscriptstyle } \mathrm{cut}}$$\end{document} affects the total cross section for Drell–Yan Z production and Higgs boson production via gluon fusion at the LHC.


Introduction
The current precision-physics program at the Large Hadron Collider (LHC) requires Standard Model (SM) theoretical predictions at the highest accuracy. Data belonging to "benchmark" processes, which are measured with the utmost precision at the LHC, need to be tested against theoretical results at the same level of accuracy. This is not only important for the extraction of SM parameters per se, but also for searches of signals of new physics, that can appear as small deviations in kinematic distributions with respect to the SM predictions. Reaching the highest possible level of precision is then the main goal and the calculation of perturbative QCD corrections plays a dominant role in this context. Until a few years ago, the standard for such calculations was next-to-leading order (NLO) accuracy. In recent years, a continuously-growing number of next-to-next-to-leading order (NNLO) results for many important processes has appeared in the literature, giving birth to the so called "NNLO revolution". For several "standard candles" processes, the first steps towards the calculation of differential cross sections at N 3 LO have also been taken (see e.g. [1,2]).
The computation of higher-order terms in the perturbative series becomes more involved due to the technical difficulties arising in the evaluation of virtual contributions and to the increasing complexity of the infrared (IR) structure of the real contributions. In order to expose the cancellation of the IR divergences between real and virtual contributions, the knowledge of the behaviour of the scattering ampli-tudes at the boundaries of the phase space is then a crucial ingredient and it is indeed what is used by the subtraction methods in order to work. These methods can be roughly divided into local and slicing. Among the first, the most extensively used at NLO were proposed in Refs. [3,4]. As far as the NNLO subtraction methods are concerned, the past few years have witnessed a great activity in their development: the transverse-momentum (q T ) subtraction method [5][6][7][8], the N -jettiness subtraction [9,10], the projection-to-Born [11], the residue subtraction [12,13] and the antenna subtraction method [14][15][16] have all been successfully applied to LHC phenomenology. The first application of the q T -subtraction method to differential cross sections at N 3 LO was recently proposed in Ref. [1], in the calculation of the rapidity distribution of the Higgs boson.
While a local subtraction is independent of any regularising parameter, slicing methods require the use of a cutoff to separate the different IR regions. Such separation of the phase space introduces instabilities in the numerical evaluation of cross sections and differential distributions [17][18][19][20], and some care has to be taken in order to obtain stable and reliable results.
The knowledge of logarithmic and power-correction terms in the cutoff plays a relevant role in the identification of universal structures, in the development of regularisation prescriptions and in resummation programs [21][22][23][24][25][26][27][28][29][30][31]. According to their behaviour in the zero limit, the cutoff-dependent terms can be classified into logarithmically-divergent or finite contributions, and power-correction terms, that vanish in that limit. In particular, the terms that are singular in the smallcutoff limit are universal and are cancelled by the application of the subtraction methods, while finite and vanishing terms are, in general, process dependent. However, after the subtraction procedure, a residual dependence on the cutoff remains as power corrections. While these terms formally vanish in the null cutoff limit, they give a non-zero numerical contribution for any finite choice of the cutoff.
From a theoretical point of view, the knowledge of the power corrections greatly increases our understanding of the perturbative behaviour of the QCD cross sections, since more non-trivial (universal and non-universal) terms appear. The origin of these terms can be traced back both to the scattering amplitudes, evaluated at phase-space boundaries, and to the phase space itself. Thus, several papers have tackled the study of power corrections in the soft and collinear limits [32][33][34], while studies in the general framework of fixed-order and threshold-resummed computations have also been performed [35][36][37][38][39][40][41][42].
From a practical point of view, the knowledge of the power corrections makes the numerical implementation of a subtraction method more robust, since the power terms weaken the dependence of the final result on the arbitrary cutoff. This is not only valid when the subtraction method is applied to NLO computations, but it is numerically more relevant when applied to higher-order calculation, as pointed out, for example, in the evaluation of NNLO cross sections in Refs. [19,20].
In this paper we consider the production of a colourless system at next-to-leading order in the strong coupling constant α S . In particular, we discuss Drell-Yan (DY) V production and Higgs boson production in gluon fusion at NLO, in the infinite top-mass limit. We impose a transversemomentum cutoff, q cut T , on the colourless system and we compute the power corrections in the cutoff, up to order (q cut T ) 4 , for the inclusive cross sections. The knowledge of these terms will shed light upon the non-trivial behaviour of cross sections at the boundaries of the phase space, and upon the resummation structure at subleading orders. In addition, it allows to have a better control on the cutoff-dependent terms, when the q T -subtraction method of Ref. [5] is applied to the numerical calculation of cross sections, allowing for a use of larger values of the cutoff. We also describe a processindependent procedure that can be used to compute the allorder power corrections in the cutoff.
The outline of this paper is as follows. In Sect. 2 we introduce our notation, and we briefly summarize the expressions of the partonic and hadronic cross sections, in a form that is suitable for what follows. In Sect. 3 we outline the calculation we have done and in Sect. 4 we present and discuss our analytic results for V and H production, along with a study of their numerical impact. We draw our conclusions in Sect. 5. We leave to the appendixes all the technical details of our calculation.

Kinematics and notation
We briefly introduce the notation used in our theoretical framework and we recall some kinematic details of the calculations presented in this paper.

Hadronic cross sections
We consider the production of a colourless system F of squared invariant mass Q 2 plus a coloured system X at a hadron collider (2.1) We call S the hadronic squared center-of-mass energy and we write the hadronic differential cross section for this process as f a/b are the parton densities of the partons a and b, in the hadron h 1 and h 2 respectively, and dσ ab is the partonic cross section for the process a + b → F + X . The dependence on the renormalisation and factorisation scales and on the other kinematic invariants of the process are implicitly assumed.
In Appendix A we have collected all the formulae for the calculation of the partonic cross section. Using Eqs. (A. 19) and (A.20), we write the hadronic cross section as where s is the partonic center-of-mass energy, equal to We have also made explicit the dependence on z, the ratio between the squared invariant mass of the system F and the partonic center-of-mass energy, and on q T , the transverse momentum of the system F with respect to the hadronic beams. Using Eqs. (2.5) and (2.3) and integrating over x 2 we obtain We then introduce the parton luminosity L ab (y) defined by so that we can finally write (2.8)

Partonic differential cross sections
In this section we recall the formulae for the first-order real corrections to the Drell-Yan production of a weak boson V (W or Z ) and to the Higgs boson production in gluon fusion, in the infinite top-mass limit.
The partonic cross sections dσ ab (q T , z)/dq 2 T in Eq. (2.8) are computable in perturbative QCD as power series in α S The Born contribution dσ (0) (q T , z)/dq 2 T and the virtual contributions to dσ (1) ab (q T , z)/dq 2 T are proportional to δ(q T ). Applying the formulae detailed in Appendix A, with a little abuse of notation, 1 we can write the partonic differential cross sections for the real corrections to V and H production as where is the Born-level cross section for the process qq → V , with c W the cosine of the weak angle and g, g v , g a the weak, the vector and the axial coupling, respectively. With N c we denote the number of colours, C F = N 2 c − 1 /2N c = 4/3 and T R = 1/2. For W production, the flavours of the quarks in Eqs. (2.10)-(2.12) are different, and the corresponding Cabibbo-Kobayashi-Maskawa matrix element has 1 In the rest of the paper we deal only with the real corrections to V and H production. We then use dσ (1) ab (q T , z)/dq 2 T to indicate them.
to be included in Eq. (2.12). The expression of the Altarelli-Parisi splitting functions p qg (z) andp qq (z) are given in Appendix D.
H production (2.14) where is the Born-level cross section for the process gg → H , with v the Higgs vacuum expectation value and C A = N c = 3. The expression of the Altarelli-Parisi splitting functions p gq (z) andp gg (z) are given in Appendix D.
We notice that the terms proportional to the Altarelli-Parisi splitting functions in Eqs. (2.10)-(2.14) embody in a single expression the whole infrared behaviour of the amplitudes, i.e. their soft and collinear limits. The structure of these terms was derived in a completely general form, from the universal behaviour of the scattering amplitudes in those limits, in Ref. [59].

Description of the calculation
In the small-q T region, i.e. q T Q, the real contribution to the perturbative cross sections of Eqs. (2.10)-(2.14) contains well-known logarithmically-enhanced terms that are singular in the q T → 0 limit [21][22][23][24][25][26][27][28][29][30]. In the context of inclusive NLO fixed-order calculations, the logarithmic terms are cancelled when using the subtraction prescriptions. For more exclusive quantities, such as the transverse-momentum distribution of the colourless system, the same logarithmic terms need to be resummed at all orders in the strong coupling constant to produce reliable results. Although our studies are of value in the context of the transverse-momentum resummation, here we limit ourselves to the case of inclusive fixed-order predictions at NLO, leaving the resummation program to future investigations. In this paper we compute power-correction terms to the cross section that, although vanishing in the small-q T limit, may give a sizable numerical contribution when using a slicing subtraction method.
To explicitly present the perturbative structure of these terms at small q T , it is customary in the literature [47,59] to compute the following cumulative partonic cross section, integrating the differential cross section in the range 0 ≤ q T ≤ q cut T , The cross section in Eq. (3.1) receives contributions from the Born and the virtual terms, both proportional to δ(q T ), and from the part of the real amplitude that describes the production of the F system with transverse momentum less than q cut T . The virtual and real contributions are separately divergent and are typically regularised in dimensional regularisation. Since the total partonic cross section is finite and analytically known for the processes under study, following what was done in Refs. [60,61], we compute the above integral aŝ where q max T is the maximum transverse momentum allowed by the kinematics,σ tot ab (z) is the total partonic cross section andσ > ab (z) is the partonic cross section integrated above q cut T . The advantage of using Eq. (3.2) is that the partonic cross section integrated in the range 0 ≤ q T ≤ q cut T is obtained as difference of the total cross section (formally free from any dependence on q cut T ) and the partonic cross section integrated in the range above q cut T of Eq. (3.4). Since q T > q cut T > 0, the last integration can be performed in four space-time dimensions, with no further use of dimensional regularisation. In Refs. [60,61] the computation of the cumulative cross section was performed in the limit q cut T Q, neglecting terms of O (q cut T ) 2 on the right-hand side of Eq. (3.2). In this paper, we compute these terms up to O (q cut T ) 4 included.

q T -integrated partonic cross sections
In this section we present the results for the partonic cross section in Eq. (3.4), integrated in q T , from an arbitrary value q cut T up to the maximum transverse momentum q max T allowed by the kinematics of the event, given by at a fixed value of z. The integrations are straightforward and do not need any dedicated comment. To lighten up the notation, we introduce the dimensionless quantity 2 that will be our expansion parameter in the rest of the paper, and we define that will allow us to write the upcoming differential cross sections in a more compact form.
V production We do not consider the process qq → Hg since it is not singular in the limit q T → 0 and the corresponding analytic/numeric integration in the transverse momentum can be performed setting explicitly q cut T = 0. A couple of further comments about the above expressions are also in order. In first place, the part of the cross sections proportional to the Altarelli-Parisi splitting functions in Eqs. (3.8)-(3.11) has a universal origin, due to the factorisation of the collinear singularities on the underlying Born. The rest of the above cross sections is, in general, not universal. In addition, for Higgs boson production, the NLO cumulative cross sections that we have computed coincide exactly with the jet-vetoed cross sections σ veto ( p veto T ) of Ref. [62], provided we identify q cut T = p veto T .

Extending the integration in z
According to Eq. (2.8), in order to compute the hadronic cross section we need to integrate the partonic cross sections convoluted with the corresponding luminosities. In the calculation of the total cross sections, the upper limit in the z integration is unrestricted and is equal to 1. When a cut on the transverse momentum q T is applied, the reality of Eqs. (3.8)-(3.11) imposes the non negativity of the argument of the square roots, i.e.
1 − π 2 T ≥ 0, (3.12) that in turn gives Since our aim is to make contact with the transversemomentum subtraction formulae, that describe the behaviour of the cross sections in the soft and collinear limits, we need to extend the integration range of the z variable up to 1, i.e. the upper integration limit of z in a Born-like kinematics. In fact, only in the z → 1 limit we recover the logarithmic structure from the soft region of the emission. In order to obtain explicitly all the logarithmic-enhanced terms in the small-q cut T limit, we have then to expand our results in powers of a. Since both the integrand and the upper limit of the integral depend on a, the naïve approach of expanding only the integrand does not work, due to the appearance of divergent terms in the z → 1 limit, that have to be handled with the introduction of plus distributions.
Using the notation of Ref. [60], we first introduce the func-tionR ab (z), defined by 14) where the upper integration limit in z in the last integral is exactly 1 andσ (0) is the partonic Born-level cross section for the production of the colourless system F. The function R ab (z) can be written as a perturbative expansion in α Ŝ where the δ(1 − z) term is the Born-level contribution, and δ B = 1 when partons a and b are such that a + b → F is a possible Born-like process, otherwise its value is 0. The coefficient functionsR ab (z) can be computed as power series in a. It is in fact well known in the literature [60] that the NLO coefficientR (1) ab (z) has the following form 3 (3.16) 3 The notation for the expansion of R (1) ab (z) follows from the number of powers of α S , log(a) and a 1 2 , i.e. and the aim of this paper is to compute the first unknown terms in Eq. (3.16) that were neglected in Refs. [60] and [61], namely R (1,m,r ) ab (z), for r up to 4 and for any m. In a way similar to what was done in Eq. (3.14) forR ab (z), we introduce the functionĜ ab (z) defined by and σ tot ab is independent of a, the coefficients of the terms that vanish in the small-q T limit in the series expansion in a of R ab (z) andĜ ab (z) are equal but with opposite sign, at any order in α S . We recall thatR ab (z) contains terms of the form δ(1 − z), coming from the Born and the virtual contributions, that are independent of a and are obviously absent inĜ ab (z).
In the rest of the paper we compute the first terms of the expansion in a ofĜ (1) ab (z), that will be obtained from the following identity We have elaborated a process-independent formula to transform an integral of the form of the first one in Eq. (3.19) into the form of the second one, producing the series expansion ofĜ (1) ab (z) in a. The application of our formula reorganizes the divergent terms in the z → 1 limit into terms that are integrable up to z = 1 and logarithmic terms in a. Since this is a very technical procedure, we have collected all the details in Appendix B, and we refer the interested reader to that Appendix for the description of the method.

Results
In this section we summarize our findings. We present in Sect. 4.1 the analytic results for theĜ (1) ab (z) functions we have In Refs. [6,60], the leading-logarithmic R computed. In the calculation of these functions, we kept trace of all the terms originating from the manipulation of the contributions proportional to the Altarelli-Parisi splitting functions, in the partonic cross sections of Eqs. (2.10)-(2.14).
These terms constitute what we call the "universal part" of our results, as detailed in Sects. 2.2 and 3.1. We will indicate these terms with the superscript "U", while the remaining terms will have a superscript "R". We stress here that the distinction between universal and non-universal part is purely formal, and it does not have a physical implication. The reason of this separation is to have hints on the general structure of the q cut T dependence of inclusive cross sections for the production of arbitrary colorless systems. We comment on the results that we have obtained in Sect. 4.2.
In Sect. 4.3 we study the numerical significance of the power-correction terms we have computed, discussing first their impact on the different production channels for Drell-Yan Z boson and Higgs boson production in gluon fusion. Then we present their overall effect, normalising the results with respect to the total NLO cross section, in order to have a better grasp on the size of these contributions.

Results for theĜ (1) ab (z) functions
We indicate withĝ ab (z) the remaining part, stripped off of a common colour factor. Our expressions forĜ (1) ab (z) contain derivatives of the Dirac δ function, δ (n) (z), up to n = 5, and plus distributions up to order 5. We report here the definition of a plus distribution of order n where g(z) has a pole of order n for z = 1, and l(z) is a continuous function in z = 1, together with all its derivatives up to order (n − 1). For completeness, we collect in Appendix E more details on the plus distributions, and the identities we have used to simplify our results. and +2 δ (2) (1 − z) a log(a) and In Eq. (4.7) we have written the (1 + z 2 ) terms coming from the numerator of thep qq (z) splitting function as (4.10) in order to keep track of the universal origin of those terms.
H production and (4.14) and In Eq. (4.16) we have written the 2(z 2 − z + 1) 2 /z terms coming from the numerator of thep gg (z) splitting function as (4.19) in order to keep track of the universal origin of those terms.

Comments onĜ (1) ab (z)
The leading-logarithmic (LL) and next-to-leadinglogarithmic (NLL) coefficients of theĜ (1) ab (z) functions that we have computed agree with the ones in the literature, along with the finite term. Their values have been known for a while [26,29] and are related to the perturbative coefficients of the transverse-momentum subtraction/resummation formulae for V [28] and Higgs boson production [63], as pointed out in Sect. 3.2. The coefficients of the terms of order a log(a) and a, and of order a 2 log(a) and a 2 are instead the new results computed in this paper. The general form of theĜ all the other coefficients being zero. We will refer to the terms in the first and in the second line of Eq. (4.20) as leading terms (LT). These terms are either logarithmically divergent or finite in the a → 0 limit. We name the terms in the sum in the third line of Eq. (4.20) as next-to-leading terms (NLT), and the terms in the fourth line as next-to-next-to-leading terms (N 2 LT), and so forth.
We notice that the NLT and N 2 LT terms are at most linearly dependent on log(a), consistently with the fact that the LL contribution is a squared logarithm. In addition, no odd-power corrections of √ a = q cut T /Q appear in the NLT and N 2 LT terms. This behaviour is in agreement with what found, for example, in Ref. [48], i.e. that, at NLO, the power expansion of the differential cross section for colour-singlet production is in (q cut T ) 2 . We do not expect this to be true in general when cuts are applied to the final state.

Soft behaviour of the universal part
The origin of some of the terms in the diagonal channels, i.e. the qq channel for V production and the gg channel for H production, can be traced back to the behaviour of the Altarelli-Parisi splitting functions in the soft limit, i.e. z → 1. In fact, in this limit, 21) 4 The notation for the expansion of G  so that Insertingp qq (z) andp gg (z) in Eqs. (3.9) and (3.11), respectively, they give rise to a contribution of the form where, following the notation of Appendix D, we have defined The details for the derivation of Eq.

The non-universal part
It is also interesting to notice that the non-universal part of theĜ (1) ab (z) functions contains terms proportional to log(a), multiplied by powers of a. These powers are controlled by the form of the non-universal parts in Eqs. (3.8)-(3.11), and to the way they enter in our generating procedure described in Appendix B. In fact, by inspecting Eq. (B.10), we see that they contribute toĝ R(1) ab with terms of the form +a 2 log(a) + a log(a) + · · · n = 1 −2 a 2 log(a) + · · · n = 2 +a 2 log(a) + · · · n = 3 −6 a 3 log(a) + · · · n = 4 +2 a 3 log(a) + · · · n = 5 (4.25) where the dots stand for power terms in a with no logarithms attached. This also explains why, for V production,ĝ qq (z) in Eq. (4.8) contain both terms a log(a) and a 2 log(a): they receive contributions from all the terms in Eq. (4.25) starting from n = 1, since Eqs. (3.8) and (3.9) contain a term proportional to (1 − z) 1 − π 2 T . Instead,ĝ R(1) gq (z) in Eq. (4.13) andĝ R(1) gg (z) in v (4.17) contain only the term a 2 log(a), since they receive contributions from the terms in Eq. (4.25) starting from n = 2, due to the fact that Eqs. (3.10) and (3.11) contain a term proportional to (1 − z) 2 1 − π 2 T and (1 − z) 3 1 − π 2 T , respectively. As far as the finite term in the diagonal channels is concerned, we notice that, in the qq channel of DY production, the first term in Eq. (4.8) happens to correspond to the first-order collinear coefficient function defined in the "hardresummation scheme", introduced in Refs. [64,65] within the q T -subtraction formalism. Instead, the first term in the gg channel of H production in Eq. (4.17) has no connection with the first-order collinear coefficient function, that is zero for this production channel. In conclusion, the structure of the terms in the non-universal part depends on the peculiar form of the differential cross sections.

Higher-order soft behaviour of the squared amplitudes
In this section, we extend the study performed in Sect We have then applied the algorithm described in this paper to each of the terms of the expansions that we have calculated, in order to compute their behaviour as a function of a. This has allowed us to trace the origin of the a log(a) and a 2 log(a) terms. Our findings are collected in the following: We reproduce the a log(a) behaviour ofĝ  Collecting our result in a table, we have: In summary, the next-to-leading-soft approximation of the exact amplitudes reproduces the a log(a) term only for V production and for the non-diagonal channel of H production. For the diagonal channel of H production, only the expansion up to next-to-next-to-leading-soft order reproduces the a log(a) term. We would like to point out that, of the three terms contributing to the N 2 LS of Eq. (A.42), only the constant one, i.e. the number 16, is needed to reproduce the a log(a) coefficient. The u/t and t/u terms do not give rise to any a log(a) contribution.
Moreover, only the expansion up to next-to-next-toleading-soft order in the diagonal channel for V production and in the non-diagonal channel for H production reproduces the a 2 log(a) coefficient. Higher orders in the expansion in the softness of the final-state parton are needed for the nondiagonal channel of V production and for the diagonal channel of H production.

q T -subtraction method
In the original paper on the q T -subtraction method [5], the expansion in α S of the transverse-momentum resummation formula generates exactly the three terms in Eq. (3.16), plus extra power-correction terms.
In the formula forR (1) ab (z) that we can build from our expression ofĜ (1) ab (z), by changing the overall sign and adding the δ(1 − z) contribution from the virtual correction, the power-correction terms are exactly those produced by the expansion of the real amplitudes. If one is interested in using our formula forR (1) ab (z) to reduce the dependence on the transverse-momentum cutoff, within the q T -subtraction method, the aforementioned extra terms need then to be subtracted from our expression ofR (1) ab (z).

Numerical results
As previously pointed out, NLO (and NNLO) cross sections computed with the q T -subtracted formalism exhibit a residual dependence on q cut T , i.e. the parameter a we have introduced in Eq. (3.6). This residual dependence is due to power terms which remain after the subtraction of the IR singular contributions, and vanish only in the limit a → 0 (limit which is unattainable in a numerical computation). In this section we discuss the residual systematic dependence on q cut T due to terms beyond LT, NLT and N 2 LT accuracy.
We present our results for Z and H production at the LHC, at a center-of-mass energy of √ S = 13 TeV. In our NLO calculations we have set the renormalisation and factorisation scales equal to the mass of the corresponding produced boson, and we have used the MSTW2008nlo partondistribution function set [66]. The mass of the Z boson m Z and of the Higgs boson m H have been set to the values 91.1876 GeV and 125 GeV, respectively.
As an overall check of our calculation, we compared the results obtained with the analytically q T -integrated cross sections in Eqs. Then, in order to study the residual q cut T dependence of the NLO cross sections for all the partonic subprocesses, we insert the expansion in Eq. (4.20) into Eq. (3.19), and we introduce the following definitions × a 2 log(a)Ĝ where we have dropped the > and (1) superscripts for ease of notation, since there is no possibility of misunderstanding in this section, because we present only the NLO results we have computed for theĜ (1) ab (z) functions. [pb] qq → Zg @ 13 TeV   .7), for the channels that contribute to Z and H production at NLO. We have expanded the luminosity functions on the basis of the Chebyshev polynomials up to order 30. In this way, the computation of the derivatives of the luminosity functions can be performed in a fast and sound way.
In the forthcoming figures, we plot the following quantities as a function of q cut T (the corresponding value of a is given on top of each figure): where σ >(1) ab is the cumulative cross section defined on the left-hand side of Eq. (3.19), obtained by integrating the exact differential cross sections of Eqs. (3.8)-(3.11). We expect that, by adding higher-power terms in a, these differences tend to zero more and more quickly when q cut T → 0. And in fact, the results shown in the following figures confirm this behaviour.
We first present our findings separated according to the partonic production channels. In all the figures presented in this paper, the statistical errors of the integration procedure are also displayed, but they are always totally negligible on the scales of the figures.
In Fig. 1 we collect the results for the aforementioned cross-section differences, as a function of q cut T , for the qg → Zq (left) and qq → Zg (right) channels, and in Fig. 2 we collect similar results for the gq → Hq (left) and gg → Hg (right) channels. As expected, NLT and N 2 LT contributions increase the accuracy of the expanded cross section, with respect to the exact one.
To give a more quantitative estimation of the powersuppressed corrections, we present results for the total   Fig. 3 Results for 1 − σ >(1) − σ /σ NLO as a function of q cut T , for Z boson production, in pp → Z j. Same legend as in Fig. 1. In the left pane, the low-q cut T region is displayed, while, in the right pane, a larger region in q cut T is shown. The total cross section at NLO for Z production, σ NLO , has been taken equal to 55668.1 pb  Fig. 4 Results for 1 − σ >(1) − σ /σ NLO as a function of q cut T , for H boson production, in pp → H j. Same legend as in Fig. 1. In the left pane, the low-q cut T region is displayed, while, in the right pane, a larger region in q cut T is shown. The total cross section at NLO for H production, σ NLO , has been taken equal to 31.52 pb hadronic cross section, normalised with respect to the corresponding exact NLO cross section σ NLO (i.e. including also the virtual contributions). The results are shown in Figs. 3 and 4, where we have used σ NLO = 55668.1 pb for Z production and 31.52 pb for H production. On the left panes we plot results in a smaller q cut T region, while, on the right panes, we extend the q cut T interval to higher values.
These plots show exactly how the residual cutoff dependence of the cross sections changes when the q T -subtraction counterterm is corrected by the NLT and N 2 LT power terms. For example, for Z production and for q cut T = 10 GeV, corresponding to a = 0.012, the LT cross section gives an estimate of the exact cross section within the 5‰, that reduces to below the 1‰ when the NLT contribution is added and becomes less than 0.01‰ when also the N 2 LT is present.  For Higgs boson production, the residual cutoff dependence is even more pronounced: in fact, at q cut T = 10 GeV, corresponding to a = 0.0064, the LT is precise within the 1% level. When the NLT is added, the precision reaches the 0.2‰, and is below 0.001‰ with the addition of the N 2 LT.
An interesting question is to estimate the impact of the universal parts of theĜ  Fig. 5, for Z production, and in Fig. 6, for H production.
Comparing these figures with the corresponding ones with the fullĜ (1,n,m) ab (z) functions, i.e. Figs. 3 and 4 , we see that the non-universal contributions play a crucial role for Z production, while their role is minor in Higgs boson production. This is due to the fact that the non-universal part in Eqs. (3.10) and (3.11) is suppressed by higher powers of (1 − z), with respect to the corresponding expression for Z production, in Eqs. (3.8) and (3.9), confirming the conclusions drawn in Sect. 4.2.2.

Conclusions
In this paper we considered the production of a colourless system at next-to-leading order in the strong coupling constant α S . We imposed a transverse-momentum cutoff, q cut T , on the colour-singlet final state and we computed the power corrections for the inclusive cross section in the cutoff, up to the fourth power. Although we studied Drell-Yan vector boson production and Higgs boson production in gluon fusion, the procedure we followed is general and can be applied to other similar cases, up to any order in the powers of q cut T . We presented analytic results, reproducing the known logarithmic terms from collinear and soft regions of the phase space, along with the finite contribution, and adding new terms as power corrections in q cut T . We found that the logarithmic terms in q cut T show up at most linearly in the powercorrection contributions, consistently with the fact that the LL contribution is a squared logarithm. In addition, no oddpower corrections in q cut T appeared in our calculation, in agreement with known results in the literature for the NLO differential cross section in colour-singlet production. We do not expect this to be true in general when cuts are applied to the final state.
Along the calculation we kept track of the origin of the newly-computed terms, so that we were able to separate them into a universal part, and a part that depends on the process at stake. In particular we derived and identified the contribution to the universal part coming from soft radiation, present in the diagonal partonic channels for Z and H production. We could also explain some features about the presence of power-suppressed logarithmic terms, appearing in the non-universal part of the power corrections. Furthermore we showed that the knowledge of the squared amplitudes at the next-to-leading-soft approximation is not enough to predict the (q cut T ) 2 log q cut T behaviour of the power corrections. The same conclusion can be drawn for the knowledge of the next-to-next-to-leading-soft approximation in predicting the (q cut T ) 4 log q cut T power correction. We also studied the numerical impact of the power terms in the hadronic cross sections for Z and H production at the LHC at 13 TeV, both by keeping track of the different par- We plotted the behaviours of the cross sections while adding more and more orders of the power-correction terms, as a function of q cut T , and comparing them with the exact cross sections. For example, in Z production and for a value of q cut T = 10 GeV, the sensitivity on the cutoff can be reduced from 5 to 0.01‰, when adding the (q cut T ) 4 contributions to the (q cut T ) 2 ones. Higgs boson production suffers from a larger sensitivity on the cutoff, and the dependence goes from 1% to 0.2‰, when all the power corrections we computed are added. By performing the same numerical comparisons for just the universal part of the power corrections, we showed that the non-universal contributions play a crucial role for Z production, while their role is minor in Higgs boson production.
The knowledge of the power terms is crucial for understanding both the non-trivial behaviour of cross sections at the boundaries of the phase space, and the resummation structure at subleading orders. Within the q T -subtraction method, the knowledge of the power terms helps in reducing the cutoff dependence of the cross sections. While the application of the q T -subtraction method in NLO calculations is superseded by well-known local subtraction methods, at NNLO it still plays a major role, also in view of the fact that, as shown in Refs. [1,19], the sensitivity to the numerical value of the cutoff increases at higher orders.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The manuscript has no associated data.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A Partonic phase space and partonic cross sections at NLO
In this section we collect some basic formulae used to compute the cross section of the partonic process where a, b and c are quarks or gluons, in a combination compatible with the production process of the colourless system F. In parentheses, the four momenta of the particles are given.

A.1 Partonic phase space
The standard Mandelstam relativistic invariants are given by and we also define the threshold variable z The phase-space volume with the appropriate flux factor is given by Since the colourless system recoils against the final coloured parton, their transverse momenta are equal. Calling θ the angle between p 1 and k, we can write Inverting the system, we find the relations which lead to an expression of the phase-space volume in terms of q T and t On the other hand, using the identity t u = s q 2 T , (A. 10) we can write the argument of the δ function in Eq. (A.4) as As a consequence, it is possible to write We then add a dummy integration over the z variable, that allows us to rewrite the phase-space volume as

A.2 Partonic cross sections
We can write the partonic cross sections for a 2 → 2 process as where M is the amplitude for the partonic process, that in general can be written as a function of the Mandelstam variables s, t and u. From Eqs. (A.3) and (A.10) we can express s and u as functions of z, q T and t, and using Eq. (A.16) we can write dσ = 1 16π where we have defined (A.20) A.3 Squared amplitudes and their soft limit In this section, for completeness, we collect the squared amplitudes |M(s, t, u)| 2 for V +jet and H +jet, at the lowest order in α S , stripped off of trivial coupling, colour and spin factors. The normalization of the following amplitudes is such that is exactly the numerator of Eqs. (2.10)-(2.14).
Together with the exact squared matrix elements, we give also the soft behaviour of the amplitudes, using the energy k 0 of the final-state parton as expansion parameter. We have computed the soft expansion adopting the following procedure: we first got rid of s in favour of Q 2 , t and u using the identity (A. 22) In this way, the only dependence on the energy of the finalstate parton is through t and u, that are linearly dependent on k 0 (see Eq. (A.2)). Then we perform a Laurent expansion in k 0 , and define leading soft (LS) the term proportional to the highest negative power of k 0 , next-to-leading soft (N 1 LS) the subsequent term, and so on. Finally we re-express all the soft-expansion contributions in terms of t and u. As a result, at each order of the expansion, all the terms proportional to a given power of k 0 are included, and only them. This is an unambiguous way to define the softness order of the expansion.
V production The exact squared amplitude is given by with soft-expansion terms The exact squared amplitude is given by with soft-expansion terms We notice that, for q(q) + g → V + q(q) production, the LS term of Eq. (A.24) has only one negative power of k 0 , while in q +q → V + g the LS term of Eq. (A.30) has two negative powers of k 0 , in agreement with the eikonal approximation for soft-gluon emission.
H production The exact squared amplitude is given by with soft-expansion terms The exact squared amplitude is given by with soft-expansion terms Similar conclusions to V production can be drawn for H production: one negative power of k 0 in g + q(q) → H + q(q), see Eq. (A.35), and two negative powers for g + g → H + g, see Eq. (A.40).

B Process-independent procedure for extending the z integration
In this section we describe the procedure followed to extend the integration range of the z variable up to 1, as displayed in Eq. (3.19), performing an expansion in a. We consider an integral of the following form , l(z) well behaved for τ ≤ z ≤ 1 and l(z) = 0 for z < τ. We also assume that l(z) is C ∞ , so that we can derive it as many times as necessary.
We then suppose that g(z) can be written as an expansion in (negative) powers of (1 − z) where, in z = 1, the g i (z, a) are not singular or have an integrable singularity. If not identically zero everywhere, the g i (z, a) are different from 0 for z = 1 and i ≥ 1, and, in general, the g i (z, a) functions contain growing powers of a as i increases. The point we would like to make here is that the righthand side of Eq. (B.2) is convergent only for z ≤ 1 − f (a), and it converges to g(z). For z > 1 − f (a) the series does not converge to g(z), otherwise g(z) would be defined in this region too.
We also assume that we can exchange the order of integration and summation of the series, and we write Eq. (B.1) as where we have extended the z-integration down to 0, since l(z) = 0 for z < τ. Each term of the series can now be manipulated as shown in the following. is finite and poses no problems.
where we have added and subtracted the first term of the Taylor expansion of l(z) around the point z = 1, and performed straightforward manipulations of the integration limits. The first and second integrands in the above equation are well behaved when z → 1, since the numerator goes to zero at least as fast as (1−z), cancelling the divergence of the denominator.

I 2
In a similar way, we can manipulate I 2 to have where we have added and subtracted the first two terms of the Taylor expansion of l(z) around z = 1. Again the first two integrands are finite when z → 1, since the numerator

Final expression
The same procedure can be applied to all the integrals I n and leads to the final result (1 − z) 2 + · · · (B.9) (B.10) Notice that inĨ 2 the sum of the terms of the series add up to give back g(z), since the upper integration limit is 1 − f (a), so that we are within the region of convergence of the series. The integrals inĨ 2 have to be evaluated exactly analytically, and this is the harsh part of the calculation. The integrals inĨ 3 can instead be computed by performing an expansion in a, and this part of the calculation poses no problems. Examples of resolution of these integrals are given in Appendix C.
Finally, by using of the plus distributions defined in Appendix E, we can writeĨ 1 in a more compact form This completes our process-independent procedure for the manipulation of the integral in Eq. (B.1).

C Detailed derivation of the results for V and H production
In this section we present in detail how we applied the method described in Appendix B to perform the series expansion in a for every production channel of the processes at stake. In particular, we specify for each channel which functions are assumed to be the l(z) and g(z) functions of Eq. (B.1). For ease of notation, in the following sections, the subscripts of the parton luminosities are suppressed, since any misunderstanding is prevented by the title of the section itself.
Also, in the summary of each of the following sections, a distinction is made while separating the final result in a universal and a non-universal part. As detailed in Sects. 2.2 and 3.1, the contributions proportional to the Altarelli-Parisi splitting functions constitute what we call the universal part of the results. The remaining ones constitute the nonuniversal one.
C.1 V production: q g channel The relevant integral, corresponding to that in Eq. (B.1), is We can express I as the sum of three integrals where and for each of the three integrals we apply the procedure detailed in Appendix B.

Integral I a
We define and expanding g(z) according to Eq. (B.2), we have  Then, writing I U and I R in the form We can express I as the sum of three integrals where (C.27) (C.28) and for each of the three integrals we apply the procedure detailed in Appendix B.

Integral I a
We define and expanding g(z) according to Eq. (B.2), we have

Integral I b
We start by separating I b into two further integrals, writing and for each of them we follow our integration and expansion procedure.
• Integral I b1 We define and we deal with this case as with a case with g 0 = 0, g 1 (z, a) = 1 and all the other g i functions equal to 0. Then, we perform the integrations in Eqs. (B.9)-(B.11).

• Integral I b2
We define and we deal with this case as with a case with g 0 = 0, g 1 (z, a) = log(1 − z) and all the other g i functions equal to 0. Then, we perform the integrations in Eqs. (B.9)-(B.11).

Integral I c
We define

Integral I b2
We define and we treat this case as the case with g 0 (z, a) = 0 and  (C.122)

D Altarelli-Parisi splitting functions
The zero-order Altarelli-Parisi splitting functions are defined as The unregularised Altarelli-Parisi splitting functions are given bŷ

E Plus distributions
We define a plus distribution of order n as where g(z) has a pole of order n for z = 1, and l(z) is a continuous function around z = 1, together with all its derivatives up to order (n − 1). For example, the first three plus distributions read