Infrared Renormalons in Kinematic Distributions for Hadron Collider Processes

Infrared renormalons in Quantum Chromodynamics are associated with non-perturbative corrections to short distance observables. Linear renormalons, i.e. such that the associated non-perturbative corrections scale like one inverse power of the hard scale, can affect at a non-negligible level even the very high-energy phenomena studied at the Large Hadron Collider. Using an Abelian model, we study the presence of linear renormalons in the transverse momentum distribution of a neutral vector boson $Z$ produced in hadronic collisions. We consider a process where the $Z$ transverse momentum is balanced by a sizable recoil against a coloured final state particle. One may worry that such a colour configuration, not being azimuthally symmetric, could generate unbalanced soft radiation, associated in turn with linear infrared renormalons affecting the transverse momentum distribution of the vector boson. We performed a numerical calculation of the renormalon effects for this process in the so-called large $b_0$ limit. We found no evidence of linear renormalons in the transverse momentum distribution of the $Z$ in the large transverse-momentum region, irrespective of rapidity cuts.


Introduction
Due to the absence of clear new physics signals at the Large Hadron Collider (LHC), precision calculations for LHC-physics processes have become increasingly important in recent years, in order to aid in the search for new physics phenomena that manifest themselves as modest deviations of the production and decay properties of the Standard Model particles. Improving the precision of QCD calculations is particularly demanding, given the size of the strong coupling constant, and a considerable effort is under way to improve the precision of QCD calculations of collider processes with the inclusion of two-loop, and even three-loop corrections.
The transverse momentum spectrum of the Z boson is one of the most precise observables measured at the LHC. More specifically, the normalized distributions are measured with a precision that reaches the sub-percent level in the low-intermediate values of the transverse momentum [1][2][3][4]. The uncertainties associated with theoretical calculations are, however, still at the percent level. The process of Z+jet production in hadronic collisions has been computed at NNLO QCD [5][6][7]. The current state of the art is given by NNLO+N 3 LL [8], and the resummation effects have been proven to have a large effect in the region where the Z boson transverse momentum is small, which is also plagued by important non perturbative effects [9][10][11]. The ATLAS measurement at 13 TeV [4] is in excellent agreement with the NNLO+N 3 LL result for p Z T < 30 GeV, while some tension at the level of few percent is found for larger values, where the impact of resummation is negligible. This tension is similar in size to the residual scale uncertainty, and, furthermore, it is not observed in the 8 TeV data [12]. These measurements may have important implications for constraining the strong coupling constant and the PDFs at the LHC (see e.g. Ref. [13]).
In this article we investigate the presence of non-perturbative effects in the Z-boson transverse-momentum distribution in the moderately large transverse momentum region, where resummation effects should not play an important role. We consider both the single and double differential cross sections where p T and y are the transverse momentum and rapidity of the Z boson, compute explicitly their cumulants σ(p T > p cut ), σ(p T > p cut , 0 < y < y cut ), and look for linear renormalons in the result. If we assume that this distribution is affected by linear non-perturbative corrections, their natural size would be of order Λ/p T , where Λ is a typical hadronic scale and p T is the Z transverse momentum, thus leading to values that may easily reach the few percent level, larger than the present experimental error and the current theoretical accuracy. A further concern stems from the fact that the soft-radiation pattern in the Z+jet process is not azimuthally symmetric, and infrared renormalons are related to soft radiation. If we believe that we can model renormalon effects as being equivalent to the emission of a soft particle with transverse momentum of order Λ, it would be reasonable to assume that they may affect linearly the transverse momentum of the recoil system, and thus also the transverse momentum of the Z boson.
There is a well-known relation between non-perturbative corrections and the factorial growth of the coefficients of the perturbative expansion in field theories. An observable R, in a generic renormalizable quantum field theory, can be expressed as a series in the renormalized coupling α R(α) = n r n α n . (1.2) Generally, the coefficients r n grow factorially for large orders, and one must worry about the associated ambiguities in the definition of the sum in (1.2). There are three known sources of factorial growth: ultraviolet (UV) renormalons, infrared (IR) renormalons and instantons. In QCD, IR renormalons lead to a factorial growth of the form where p is a positive integer, b 0 is defined as usual and n l is the number of light flavours. Eq. (1.3) implies that the terms of the perturbative expansion grow by a factor n(2b 0 /p)α s as the order increases, and thus reach a minimum when n ≈ p/(2b 0 α s ). Using Stirling's formula (n! ≈ n n exp(−n)), the size of the minimal term is easily estimated to be where µ is the renormalization scale. This is the well known connection between IR renormalons and power corrections. It tells us that the size of the minimal term of the perturbative expansion plays the role of a power suppressed ambiguity in the value of the corresponding observables. In our case, as already mentioned earlier, we worry about the case when p = 1, that will be referred to in the following as linear renormalons. We also stress that, due to the relatively large value of the QCD coupling constant, the minimal term in the perturbative expansion may be reached quite early, and one may worry that, in some cases, already at the N 3 LO level linear renormalons may become relevant. UV renormalons also lead to a factorial growth of the same form as in eq. (1.3), with p assuming negative integer values starting with −2. In this case, formula (1.5) still holds if we replace p with |p|, and the minimal term is, at worse, of order Λ 2 /µ 2 . Furthermore, the perturbative expansion is alternating in sign, and can in principle be resummed with Borel techniques. Instantons lead instead to a much stronger power suppression, and need not be considered here. 1 General arguments based upon field theory assure us that linear renormalons cannot be present in certain quantities. This is the case for observables that admit an Operator Product Expansion (OPE), and are such that there are no operators of dimension higher by one power with respect to the leading contribution. In fact, IR renormalons in hard processes originate from subgraphs with low momenta, and, if the process in question admits an OPE, it must be possible to organize the Feynman graphs in terms of expectation values of local operators [19,21]. In case of processes that do not have a known OPE, we lack a safe guideline for the classification of non-perturbative corrections, and thus also of renormalon effects. In such cases, one has to resort to assumptions and approximations in order to gain some insight into the problem.
A very valuable method for the study of IR renormalons is the so-called large-b 0 approximation (see [20] and references therein). In essence, the method considers an alternative theory such that the renormalon structure can be fully computed, i.e. an Abelian gauge theory in the limit of a large number of fermions n f . In this theory the leading terms, of order (α s n f ) k , are fully calculable, and the theory exhibits IR and UV renormalons. In order to have an estimate of the renormalon effects in the full original theory, at the end of the calculation one replaces the b 0 in the Abelian theory (that is proportional to −n f ) with the full QCD b 0 .
We recall that the authors of ref. [22] showed that, in the large-b 0 approximation, the Drell-Yan cross section does not have linear renormalons. This result was used to argue that, contrary to previous claims, soft gluon resummation cannot give solid evidence of the presence of linear renormalons in the total Drell-Yan cross section. Later on, in ref. [23], it was shown that under certain assumptions the rapidity distribution of Drell-Yan pairs is free of linear power corrections. This result can be shown to imply that in the large-b 0 approximation there are no linear renormalons in the rapidity distribution of Drell-Yan pairs. Figure 1: Typical diagram contributing to the transverse momentum distribution of a Z boson, which is represented by a blue zigzag line. The soft radiation pattern (dashed line) is associated with the dipoles formed by the outgoing hard gluon (red curly line) and the initial state quarks, and is not azimuthally symmetric.
In this work we study in the large-b 0 approximation the transverse momentum distribution of a vector boson, and in particular, we want to determine whether it is affected by linear renormalons. Our interest in this process stems from the fact that, unlike the inclusive Drell-Yan case, the soft emission pattern in the production of a vector boson in association with a hard jet is not azimuthally symmetric, as shown schematically in fig. 1. Under these circumstances, we can expect that soft gluons may induce linear renormalon corrections to the transverse momentum distribution of the vector boson, since they are not emitted according to an azimuthally symmetric pattern. We also notice that in ref. [24], it was found that certain leptonic observables in top production and decay are affected by linear renormalons, thus casting doubts on the common assumption that leptonic observables should be less affected by non-perturbative corrections with respect to hadronic ones.

On the small p T distribution of vector bosons
The transverse momentum distribution of vector bosons in the region of small transverse momenta has been the subject of considerable theoretical activity since the early days of perturbative QCD [9][10][11] up to the present days [12,[25][26][27][28][29][30]. It has a very peculiar perturbative expansion in QCD, starting with a δ 2 (p T ) at leading order, and receiving singular contributions due to soft and collinear gluon emissions at higher perturbative orders. It is dealt with by resumming the singular contributions to all orders in the perturbative expansion. This turns out to yield (for very large masses) finite results down to zero transverse momentum, according to a mechanism first exposed by Parisi and Petronzio in ref. [10]. It was however immediately recognized that some non-perturbative input is required in these calculations, since the resummation of soft-collinear gluons must be cut-off at transverse momenta of the order of typical hadronic scales. Also, by intuitive reasoning, one expects that Fermi motion effects should smear the δ 2 (p T ) leading order contribution (and the higher order singular ones) to a transverse size of the order of a Fermi. 2 Considering for the sake of argument the smearing due to Fermi motion, let us call f ( q T )d 2 q T the primordial transverse momentum distribution of the quarks in a hadron. If f ( q T ) is dominated by values of q T of the order of a hadronic scale, f can be expanded in moments where Λ is a typical hadronic scale, corresponding in parameter space to the behaviour Notice that linear terms in ∂ p T (or b) are naturally absent (at least as long as one does not consider spin structures in the distribution), and in fact the behaviour given in eq. (1.7) has been advocated as early as ref. [10], where the formf (b) = exp(−Λ 2 b 2 ) was proposed. The same form also appears in several subsequent works. 3 When convoluted with perturbative mechanisms giving rise to the transverse momentum of the vector boson, the form of eq. (1.6) leads to corrections of order Λ 2 /p 2 T . The absence of linear corrections in this context has also a rather simple intuitive explanation. The primordial transverse momentum smearing gives a transverse kick, of the order of typical hadronic scales, to the perturbative distribution. However, it is azimuthally symmetric. Thus, its first-order effects cancel out, leaving only quadratic corrections. On the other hand, the non-perturbative corrections introduced in this context are all that is needed in order to regularise the ill-behaved perturbative expansion in the small p T region. They cannot be claimed to be the only non-perturbative corrections that are relevant to the problem, especially if we move to regions of phase space where the QCD perturbative expansion is totally well-behaved.
At variance with the case of the small transverse momentum calculations, we deal with a process that has a well behaved perturbative expansion in the strong coupling constant. We do not make any assumption regarding the origin of non-perturbative corrections. We compute the full cross section, and look for renormalon effects in the full result. We can only do this, however, by considering a simpler theory, i.e. the large n f limit of QCD, where we can compute the cross section exactly, at all orders in the perturbative expansion.
We expect that linear effects may appear only if the soft radiation pattern is not azimuthally symmetric, and thus we considered a process where such asymmetric configuration is realized. For example, we did not consider the process where two incoming partons produce a Z plus a photon with large transverse momentum. Such process does yield an azimuthally symmetric pattern for gluon radiation, and thus we do not expect linear renormalons to arise there.
We also remark that the thrust of a jet does receive linear power corrections associated with linear renormalons (see for example eq. (5.56) in ref. [20]). Since in the process we consider the Z is recoiling against a jet, the worry that such linear corrections may affect the Z transverse momentum distribution is not without ground. Figure 2: A Born diagram for Z production in photon-quark collision. The incoming photon is represented by the blue wavy line. This process has an azimuthal asymmetric colour pattern of soft emission, and is suitable for testing for the presence of linear renormalons in the Z transverse momentum distribution.

Our calculation
A calculation of the large-b 0 corrections to the process depicted in fig. 1 is too demanding, and in fact, at the moment, no large-b 0 corrections to processes already involving a gluon emission or exchange has been ever carried out. Such calculation would inevitably involve dressed gluon propagators joining into a three-gluon vertex, also including a vertex correction formed by a quark triangle graph (such correction is in fact of order g s (g 2 s n f )). Thus, we will instead compute the process depicted in fig. 2, i.e. the production of a Z boson in photon-quark collisions. More specifically, we assume that the incoming photon couples to a single quark flavour (that we will call for definiteness a d quark) but there is a large number n f of light quark flavours that we denote with the symbol q. We can thus compute this process in the large n f limit, and turn to the large b 0 limit in the usual way, i.e. by replacing the b 0 in the Abelian theory with the full QCD b 0 . In this process (as in the realistic QCD one) the soft emission pattern is not azimuthally symmetric, and we can investigate whether linear IR renormalons are present in the Z transverse momentum distribution due to recoil against soft emissions. The radiative corrections that one needs to compute are represented schematically in fig. 3. The solid blob insertion in the gluon propagator represents the inclusion of all corrections given by a fermion loop, as represented by the recursive graphic equation . (1.8) The inclusion of all corrections embodied in eq. (1.8) amounts to considering α s n f to be of order 1. Because of this reason, the gluon splitting into fermion pairs must also be included, since, when squared, it also gives rise to a factor α s n f . We also notice that we do not need to worry about the interference of the fermions arising from gluon splitting with the fermion coming from the initial hadron, since this term is suppressed by a factor of n f . We have developed a parton-level generator for the computation of this process that can be used for the computation of any IR finite differential distribution, provided a cut in the Z transverse momentum is included. The rest of the paper is organized as follows. In Section 2 we describe the procedure that we adopt in order to perform the calculation. In essence, the result is expressed in terms of a corresponding next-to-leading correction to the Born process in a fictitious theory where the gluon has a finite mass, plus a "remainder" correction that accounts for the difference of the calculation in the massive case with respect to the one performed in the full theory [34][35][36]. This difference does not contribute to quantities that do not depend upon the kinematics of coloured objects in the final state, and thus does not contribute to the transverse momentum distribution of the Z boson. The procedure for the removal of the initial state collinear singularities is discussed here.
In Section 3 we give further details on the setup of the calculation. In particular, it is shown that the real cross section must be partitioned into three terms, each of them appropriate to one of the singular regions of the real cross section, and each term (requiring an appropriate integration importance sampling) is computed independently. Finally, in section 4 we present the results of our calculation. In sec. 5 we give our conclusions.

The method
We assume we are calculating a cross section with a given set of cuts, that we represent with a function of the kinematic configuration Θ(Φ), that takes the value 1 if the cuts are satisfied, and zero otherwise. The cross section for our process (borrowing from the notation of ref. [37]) is given by where B represents the Born term for the process dγ → dZ (see fig. 2), V refers to the virtual correction to the dγ → dZ process, R g refers to the process dγ → dZg, and R qq to the process dγ → dZqq (see the diagrams labelled as (v), (g) and (qq) in fig. 3). The C ⊕ term is the subtraction of the collinear singularities arising in both R g and R qq when the gluon or the qq pair becomes collinear with the incoming quark. Both these collinear configurations are associated with the same underlying Born dγ → dZ, and thus can be represented by a single term. The C g and C qq terms represent the collinear counterterms associated with the collinear singularity due to the photon splitting into a dd pair in the R g and R qq contributions respectively. We have defined where p ⊕ , p are the momenta of the incoming positive and negative rapidity hadrons, x ⊕ , x are the momentum fractions of the incoming light quark and photon, and k Z and k d are the momenta of final state Z and light quark d. Furthermore where B is the squared amplitude for the process divided by the flux factor, and f d (f γ ) are the incoming quark (photon) distribution functions. The definitions of dΦ g and dΦ qq , as well as dΦ g , dΦ qq and R g , R qq are fully analogous to eqs. and Thus, the two C collinear subtractions correspond to diagrams (g) and (qq) of fig. 3, when the outgoing quark connected to the incoming quark line becomes collinear to the incoming photon. The D g and D qq amplitudes correspond to processes initiated by a qq collision with a Zg and a Zqq final state respectively, represented in fig. 4, and Φ g D , Φ qq D are the corresponding phase spaces. The C dd and C dγ functions are the universal collinear divergent factors for the d → d + X and γ → d + X splittings.
D g D qq Figure 4: Representative graphs of the dd initiated processes arising as subprocesses in the collinear limit for the incoming photon splitting, contributing respectively to the D g and D qq squared amplitudes.
We remark that no interference arises in our approximation from the final state down quark connected to the incoming quark line and the final state quarks arising from gluon splitting. In fact, such interference term would be down by a factor of n f .
The diagrams of fig. 3 are affected by both ultraviolet and infrared divergences, that we regulate using dimensional regularization. It turns out that ultraviolet divergences cancel, since they are all associated with vertex and propagator corrections in an Abelian model. For the collinear subtraction C , C dγ is given by where α is the electromagnetic coupling and C d is the down-quark charge. This must be accompanied by the rescaling µ 2 F → µ 2 F exp(−γ E )/(4π) according to the MS subtraction prescription. The MS subtraction term for the C ⊕ collinear singularity is much more complicated, since it must include all corrections of order α s n f , i.e. all vacuum polarization insertions in the emitted gluon, and the gluon splitting into a qq pair. Here we avoid this problem (following ref [22]), and perform our collinear subtraction using the so-called DIS scheme, that is defined by requiring that the structure function F 2 has the expression (where C i are the electric charges of the species i, and i runs over all quarks and antiquarks) to all orders in perturbation theory. In this way we simply need to compute F 2 (x, Q 2 ) in the same approximation that we have used for our process, i.e. by including all fermion loop insertions in the gluon propagator, and then express our cross section in terms of the q i in the DIS scheme. The diagrams involved are shown in fig. (5). When translating the DIS scheme cross section into the MS scheme no new linear renormalons arise. This is a consequence of the fact that F 2 obeys an operator product expansion where power corrections are controlled by the twist of the operator, and the dominant power corrections, corresponding to twist 4, are quadratic.
The computation of the large n f cross section can be technically performed by computing a similar process with a finite gluon mass λ. The result is identical to the one in Figure 5: Representative graphs for the calculation of F 2 in the large n f limit. eq. (3.22) of Ref. [24], with the inclusion of the collinear remnants, that must be included if the process involves initial-state singularities. More specifically, we have With g * we denote a gluon with mass λ, and with the superscript (λ) applied to previously defined objects we denote the same objects computed for a massive gluon with mass λ. Thus, for example, V (λ) are the virtual corrections depicted schematically in fig. 3 and labelled with (v), where the gluon propagator with the red blob is replaced by a massive gluon propagator with gluon mass λ, Φ g * denotes the phase space for the process dγ → Zdg * , and so on. In eqs. (2.20) and (2.22) we have defined k qq = k q + kq. The fact that the large n f calculation can be given in terms of cross sections involving a massive gluon is well known. The appendix of ref. [24] summarises all the technical steps that lead to this equivalence, but this is not the reference where this equivalence has been found first. It has been used in ref. [22] in order to compute the Drell-Yan and DIS process in the large n f limit. The need of the ∆ terms has been first pointed out in ref. [34], and is the only correction that is needed to connect the massive gluon calculation to the realistic case, where the virtual gluon decays into a qq pair.
For the purpose of the present work, the ∆ terms do not contribute. In fact, our cuts only depend upon the kinematics of the Z boson, and thus are insensitive to whether we replace the qq pair arising from gluon splitting with an undecayed massive gluon with the mass equal to the virtuality of the pair. Thus, the differences of theta functions in the square brackets are always zero in our case.
The C (λ) dd (z) term was computed in eq. (3.10) of ref. [22], and in terms of the quantities defined there we have The z variable in the f (1),real term is limited to the range As a consequence of the Adler sum rule, we have Formula (2.13) relates the massive calculation with the full leading n f perturbative expansion. We remark that the integrand in (2.13) represents formally a power expansion in b 0 α s with factorially growing coefficients. In fact, it turns out that T (λ) is finite as λ → 0, and it has an expansion in λ of the form with the possible presence of logarithmic enhancements in the quadratic remainder. If T (0) = 0 then the perturbative expansion of the result has a linear renormalon. Its structure is examined in detail in appendix A. The computation of the collinear subtraction arising from the collinear splitting of the incoming photon is a standard NLO subtraction, and can be handled by standard means, as those already coded in the POWHEG BOX package [38]. It is useful, however, to clarify from the beginning that we do not expect any linear renormalon corrections to the Z transverse-momentum distribution from this collinear region. In fact, as stated earlier, the T ∆ (λ) term does not contribute, and if the final state quark is collinear to the incoming photon, the transverse momentum of the Z must be balanced by the massive gluon. Since we are only considering distributions where the Z has a sizable transverse momentum, the massive gluon will also have a large transverse momentum, and under these conditions the mass correction to the process is quadratic in λ.

Details of the calculation
As illustrated in the previous section, the computation of the IR renormalon contribution requires us to perform the calculation in the context of an Abelian theory with a massive gluon. The Born diagrams for our process are represented in fig. 6. For simplicity, we have taken the Z boson to be vectorially coupled and stable. The Born and Virtual contributions are computed analytically at fixed external momenta. The computation was performed using the symbolic manipulation program MAXIMA [39], using a package for the computation of Feynman Diagrams developed by one of the authors [40]. One-loop scalar integrals have been evaluated using the COLLIER library [41]. The virtual corrections are obtained by attaching a (massive) gluon propagator to two points along the fermion line in all possible ways, avoiding however to attach both ends to the same external line. Virtual corrections to the external lines are instead dealt with according to the usual LSZ prescription, i.e. one multiplies the Born diagram by the wave function renormalization correction. The virtual contribution to the cross section is computed by multiplying the sum of the virtual amplitudes with the conjugate of the sum of Born amplitudes, and by taking twice the real part of the result. The virtual corrections are infrared finite, since the gluon mass regulates infrared divergences. However, each individual contribution is affected by ultraviolet divergences, that we regulate using standard dimensional regularization, so that gauge invariance is preserved at every step of the calculation. The sum of all contributions is ultraviolet finite, and thus one can take its limit to 4 dimensions.
Upon integration over the external momenta, the Born and Virtual contributions are finite as long as we require a lower bound on the Z transverse momentum. The numerical treatment of this stage of the calculation will be detailed further on.
The real corrections are obtained by attaching one external massive gluon in all possible ways to the Born graph. Upon phase space integration, real corrections are not infrared finite even if one requires a finite transverse momentum of the Z-boson, due to the collinear singularities arising from the photon splitting into two massless quarks, as shown in fig. 7. This kind of singularities should also be regulated in dimensional regularization, and an MS subtraction should be performed as required by the factorization formalism. In fact, this procedure is automatically implemented in the POWHEG BOX, and does not require any further work on our part.
Singularities associated with soft or collinear gluon emissions are instead regulated by the gluon mass. These are illustrated in fig. 8, corresponding to an initial and a final-  state collinear singularity, associated with a soft one. Real contributions thus behave as log 2 λ for small λ. The same log 2 λ term, but with opposite sign, comes from the virtual corrections, and cancels the real one. After the cancellation of the double logarithmic terms, the only remaining singularity is a single logarithm of λ, that cancels against an opposite contribution in the collinear counterterm (see eq. (2.19)). At this point we get a result that has a finite limit as λ → 0: the cross section goes to a constant and our task is to determine whether this constant is approached with a slope linear in λ.
It is clear that, in view of the massive cancellations involved, in order to get a convincing numerical evidence of the small λ behaviour of the cross section, the numerical integration should be performed with an appropriate importance sampling near the regions that are singular in the λ → 0 limit. Furthermore, a direct calculation in the massless case, i.e. at λ = 0, is also necessary, since it gives us a point with negligible error, that would be difficult to obtain with small values of λ. In order to perform the calculation, the real contribution is separated into three terms: where p T,d is the transverse momentum of the final-state quark, and m T,g is the transverse mass of the final-state gluon, defined as m dg is the mass of the quark-gluon final-state system, and E d , E g are the energies of the quark and gluon in the partonic center of mass. Notice that the expression is of the order of the transverse mass of the gluon relative to the final state quark. The three superscripts (1), (2) and (3) label the three regions where each of the three contributions is singular, i.e. region (1) is singular when the final state quark is collinear to the incoming photon; region (2) when the gluon is collinear to the incoming quark, and region (3) when the final state gluon is collinear to the outgoing quark. The real contributions for regions (2) and (3) are evaluated independently, and with a different parametrisation of the phase space. In case of region (2), the phase space is factorized as the product of the two body phase space formed by the final state gluon recoiling against the quark-Z system, and the two body phase space for the quark-Z system itself. In the case of region (3), the phase space is factorized in terms of the two body phase space for the system comprising the Z, and the quark-gluon system recoiling against it, with the quark-gluon system itself parametrised as a two body phase space. With these parametrisations it is easy to perform importance sampling integration near the singular configurations. We notice that both the real and virtual contributions, as well as the Born terms, are singular when the transverse momentum of the Z vanishes. However, since we are interested in the transverse momentum distribution of the Z for transverse momenta comparable to the Z mass, when computing the integral we can suppress this region with a factor proportional to the Z transverse momentum. We choose where p T,cut is a parameter close to the transverse momentum cut that we would like to apply to our cross section. The adaptive Monte Carlo integration is performed in a standard way, multiplying the Real, Born, Virtual and collinear subtraction contributions by this suppression factor, in order to obtain a convergent result. This factor is divided out when computing the cross section with cuts, in order to obtain the correct and finite result. In fact, since we always require a transverse momentum of the Z larger than a given cut, F supp never vanishes in the region of interest.

Region (1)
The contribution of region (1) can be evaluated using directly the POWHEG BOX framework, taking the real process dγ → Zdg, and the Born process dd → Zg. The real cross section is taken equal to R (1) . The Born contribution itself is absent, since we assume that our incoming hadrons have only a down quark (for the positive rapidity incoming hadron), and a photon (for the negative rapidity one) content. The Born subprocess dd → Zg, however, enters in the collinear subtraction, that is automatically performed by the POWHEG BOX to implement the factorization of collinear singularities, as well as in the collinear remnant, that is also automatically computed by the POWHEG BOX framework.
The region of phase space dominated by soft and collinear gluons is expected to be highly suppressed in region (1), to such an extent that we do not expect any linear renormalons here. In fact, the collinear and soft-collinear region associated with an ISR gluon emitted by the quark is suppressed by a m 2 T,g factor (see eq. (3.3)). In this singular limit, for small λ, R (1) behaves like where E g is the gluon energy, θ g its angle with respect to the incoming (or outgoing) quark. The gluon mass λ provides a lower cutoff on the E g and θ g integrations, so we see that this term cannot develop a linear sensitivity to λ.

Computation of the T (λ) functions
For every observable that we consider, we compute the associated T (λ) function. The calculation is performed as follows: • The λ = 0 contribution is computed by implementing our process in the POWHEG BOX framework. The contributions of regions 2 and 3 are computed together. They have singular regions associated with initial state radiation of a gluon from the (positive rapidity) incoming quark, and final state radiation of a gluon from the final state quark. The collinear singularity associated with a vanishing transverse momentum of the outgoing quark is absent in regions 2 and 3. Care must be taken to perform the subtraction of the initial state collinear singularity, that, for consistency, has to be performed in the DIS scheme. This affects the collinear remnant, and we suitably modified the POWHEG BOX code in order to comply with this request. The factorization scale, corresponding to the scale Q of the DIS scheme subtraction, is taken equal to the Z mass. The virtual corrections to the process dγ → Zd enter this calculation, and are computed directly for λ = 0 using dimensional regularization. The cancellation of the associated soft and collinear divergences with those arising from the real contributions is already implemented in the POWHEG BOX in a general way, and requires no further action.
• The calculation for λ = 0 for the regions 2 and 3 cannot be simply implemented in the POWHEG BOX framework, that was not designed to handle singularities regulated by a mass. They were thus implemented in a dedicated Fortran code. In particular, for the real contribution care was taken to adopt phase space parametrisations that are suitable to handle each region with adequate importance sampling. The virtual contribution was also evaluated independently. Rather than computing separately the real and virtual contributions of the collinear subtraction, we implemented a local cancellation of the associated soft divergence by making use of relation (2.25).
• The contribution from region 1, for either λ = 0 or λ = 0, is also implemented in the POWHEG BOX. In this case the only singular region present is the one associated with the collinear splitting of the initial state photon into a dd pair. The underlying Born for this singularity is given by the dd → Zg process, where the gluon has mass λ, different from the case of regions 2 and 3. The collinear singularity is treated in the MS scheme in this case, and the collinear remnant is automatically provided by the POWHEG BOX.

Results
As our benchmark set-up, we have taken two colliding particles with center-of-mass (CM) energy of 300 GeV. The positive rapidity incoming particle (labelled as (1)) has a parton density consisting only of down quarks, while the negative rapidity particle (labelled as (2)) has a parton density consisting only of photons, distributed as This totally arbitrary choice is only dictated by simplicity, and is adequate for our purposes. We compute the cross section for the production of a stable vector boson Z of mass M Z = 91.188 GeV, that is only vectorially coupled. The Z and γ couplings are both given by g 2 Z/γ = 4π, and the down quark is taken to have charge −1/3. 4 The Born diagrams have been computed supplying the correct colour factor (that is 1), and in the calculation of the virtual and real corrections we have included the appropriate QCD colour factor C F . We have chosen the factorization scale µ F = M Z . The renormalization scale choice does not affect T (λ).
To begin with, we show in fig. 9 the result for the T (λ) function defined in eq. (2.16), associated to the cross section for the production of a Z boson with a transverse momentum larger than 20 GeV ( fig. 9 on the left) and 40 GeV ( fig. 9 on the right) as a function of the gluon mass λ. In order to extract the slope around λ = 0, which is responsible for linear renormalons (see eq. (2.26)), we fit T (λ) using the function where the inclusion of the single and double logarithmic terms are motivated by the findings in the Drell-Yan case [22,23]. We neglected the point for λ = 5 GeV in our fitting procedure, in order to increase the quality of the fit near λ = 0. We performed two fits, one including b as fit parameter, and the other fixing it to 0, in order to assess its impact on T (λ). In Tab. 1 we report the results of the fits. We observe that the linear   coefficient has a negligible impact on the fitting functions, its size is at least an order of magnitude smaller than the coefficient of the dominant quadratic term, and its value is consistent with zero. Thus, we find no evidence of the presence of a linear renormalon in the Z boson transverse momentum distribution, and furthermore we find that the value of the corresponding coefficient, if non-vanishing, is much smaller than the coefficients of the quadratic terms.  Table 1: Results of the fit of the T (λ) function, defined in eq. (2.16) and illustrated in Fig. 9. The fit function is given in eq. (4.2). In the first fit, corresponding to the blue lines in the figures, b in unconstrained, while in the second fit, corresponding to the red lines, b has been set to 0. The last line corresponds to the associated reduced χ 2 .
We also performed a more exclusive analysis, imposing an additional cut over the rapidity of the Z boson y Z , besides the one over the transverse momentum. The results are shown in fig. 10 and in Tab. 2. Again we do not find numerical evidence of a linear sensitivity to λ, implying that the doubly differential distribution in rapidity and transverse p T c = 20 GeV, 0 < y Z < 0.6  Figure 10: As in fig. 9, supplemented with a cut on the Z rapidity 0 < y Z < y c , with y c = 0.6.  Table 2: Results of the fit of the T (λ) function defined in eq. (2.16) illustrated in Fig. 10. The fit function is given in eq. 4.2. The first fit corresponds to the blue lines, while in the second fit the linear coefficient has been set to 0 and corresponds to the red lines. The last line corresponds to the associated reduced χ 2 .
momentum is free from linear renormalons. By looking at the coefficients reported in Tabs. 1 and 2, we notice that when we set p c T = 40 GeV instead of 20 GeV we encounter larger errors in the determination of the coefficients c and d, since the corresponding contributions are suppressed by two powers of λ/p c T , and thus their relative importance diminishes for higher cuts.

Conclusions
The current LHC physics demands high precision theoretical predictions, and has promoted an unprecedented theoretical effort in pushing perturbative calculations beyond next-toleading order, and in some cases even beyond the next-to-next-to-leading order. At the current level of precision, possible non-perturbative effects that are suppressed by a single power of the hard scale can sometimes be comparable or larger in size than the current theoretical uncertainties. Unfortunately, for collider physics observables we lack a theory of even the most important non-perturbative corrections. This is unlike others frameworks, including also heavy flavour physics, where the existence of an operator product expansion allows us to classify and parametrise non-perturbative effects.
Calculations of renormalon effects in Abelian models, using the so called large-b 0 approximation, can provide a way to explore the structure of non-perturbative effects in collider physics. These model calculations have helped in the past to investigate the structure of renormalons in Drell-Yan processes [22,23], and have been used recently in getting some insights on issues regarding the precision measurements of the top mass [24]. Related methods have also contributed to a clarification of the structure of linear renormalon effects in jet physics (see [20] and references therein).
In this paper, we have used the large-b 0 approximation in order to understand whether linear renormalon effects can yield linear power corrections to the inclusive differential distribution for the production of a vector boson at the LHC, in the regime where the transverse momentum is safely in the perturbative region, and where the resummation of transverse momentum logarithms is not needed. The process we are considering is sufficiently complex, since it involves gluon radiation both from the initial and final state, and entails single and double logarithmic singularities that cancel when combining the virtual and real corrections with the factorization of initial state singularities. We have chosen a process that mimics the Z production in an Abelian model, but is such that, as in the full QCD case, the soft radiation pattern is not azimuthally symmetric, and thus, on intuitive ground, may be associated with non-perturbative recoil effects that affect linearly the Z transverse momentum. We find no evidence of linear power suppressed effects, and find instead that the renormalon structure is well represented by quadratic terms associated with some logarithmic enhancement.
Our numerical evidence gives us a useful indication, but, by its own nature, can never be considered a solid proof of the absence of linear renormalons. Nevertheless, since we found that the coefficient of the linear renormalon is much smaller than the coefficient of the quadratic one, and is consistent with zero, we can conjecture that linear renormalons are absent for the observables that we have considered.
Currently it seems that in all collider processes considered so far that involve massless flavours, the large-b 0 approximation shows that linear renormalons are absent in observables that are inclusive in the production of coloured partons. This conclusion does not hold if massive quarks are present, as reported in ref. [24]. Further analytical work is needed in order to put these conjectures on more solid ground, perhaps also shedding some light on the underlying mechanisms that lead to the formation and cancellation of renormalon effects in hadron collider physics. work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 788223, PanScales) and by the UK Science and Technology Facilities Council (grant number ST/P001246/1).

A Renormalon structure
As far as the region of λ < µ C , that is the region relevant for IR renormalons, a linear term in T (λ) will lead to the contribution (see eq. (2.13)) where P indicates that the principal part of the integral should be taken. The first two terms in eq. (A.3) are obviously analytic in a neighbourhood of a = 0. The third term is a typical Borel representation of a power expansion with factorially growing coefficients, resummed using a principal value prescription for handling the pole on the real axis. If the pole is instead handled by moving the integration contour above or below the real axis by a tiny amount, the imaginary part of the resulting integral is equal (in absolute value) to the last term in the formula. When replacing a = b 0 α s = 1/ log µ 2 C /Λ 2 , it leads to a linear power correction: We thus see that the integral on the left-hand side of eq. (A.3) cannot be cast in the form of a Borel sum with principal value prescription, because of the presence of the last term, whose origin can be traced back to the discontinuity of the arctangent in the left-hand side integral when λ = λ L , with In practical applications where one is interested in the value of the resummed expansion using the principal value prescription, a suitable θ function is added to the arctangent integrand to make it continuous, see ref. [42].