On linear power corrections in certain collider observables

We study linear power corrections O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(ΛQCD/Q) to certain collider observables. We present arguments that prove that such corrections cannot appear in observables that are inclusive with respect to QCD radiation, such as total cross sections as well as rapidity and transverse momentum distributions of color-neutral particles. Although our calculations are carried out in a simplified framework, our arguments and conclusions are applicable, with some reservations, to processes both at lepton and hadron colliders. We also show how an improved understanding of the origin of linear power corrections allows us to simplify their calculation. As an application, we compute the leading non-perturbative corrections to the C-parameter and the thrust in e+e− annihilation in a generic three-jet configuration.


Introduction
An important part of the LHC physics program consists in the exploration of phenomena that occur at distances that are between a hundred and a thousand times smaller than the size of the proton. Thanks to the celebrated properties of Quantum Chromodynamics (QCD) such as asymptotic freedom and factorization, physics at such distances can be described using perturbation theory, where elusive quarks and gluons play the role of fundamental physical degrees of freedom.

JHEP01(2022)093
Corrections to this perturbative picture are expected to be small, suppressed by ratios of the non-perturbative QCD parameter Λ QCD ∼ 0.3 GeV and the typical energy scale Q of the process (or observable) under consideration. This hard scale Q typically ranges from a few tens to a few hundred GeV. It follows that these non-perturbative effects may change perturbative predictions by a relative amount proportional to (Λ QCD /Q) n ∼ (0.01) n − (0.001) n which, depending on the value of Q and of the exponent n, varies from a percent for n = 1 and Q ∼ 30 GeV to a permille and even smaller values for larger values of Q and n.
Perturbative predictions for cross sections and distributions, on the other hand, are controlled by powers of the strong coupling constant α s (Q) ∼ 0.1 for Q ∼ 30 − 100 GeV. Currently, perturbative computations are often performed to second or even third order in the expansion in α s , leading to theoretical predictions for hard processes which are typically accurate to within one to ten percent [1]. Further development of methods for perturbative calculations in QCD may improve the precision of such theoretical predictions, perhaps by about an order of magnitude. If this happens even for a few selected processes and observables, perturbative predictions at this precision will have to be supplemented with non-perturbative corrections provided that O((Λ QCD /Q) n ) contributions with n = 1 exist for a particular process or observable. On the other hand, non-perturbative corrections with n > 1 are too small to be of any relevance for most hard processes at the LHC.
For color-singlet decay rates, deep-inelastic scattering structure functions and inclusive decays of heavy quarks [2], it is well-known that operator product expansion techniques allow one to conclude that fully inclusive observables do not receive linear power corrections. However, it is currently not known how to generalize these results to more differential observables and to the case of hadron-hadron collisions. Indeed, at present there is no full theory of non-perturbative corrections to short-distance processes at lepton and hadron colliders. It is then not possible to predict the exponent n for a generic process or observable, let alone compute the contribution of O((Λ QCD /Q) n ) terms precisely. However, it is wellunderstood that one source of non-perturbative corrections is present within perturbation theory itself. Indeed, the appearance of the Landau pole in the strong coupling constant leads to an intrinsic ambiguity when integrating over soft momenta. Since such an ambiguity will have to cancel with contributions that arise from physics beyond perturbation theory, it can be used as an estimate of, at least, some non-perturbative contributions.
It is well-known that the ambiguity related to the appearance of the Landau pole can be studied within the approximation of a large (and negative) number of massless fermion species (see ref. [3] for a review). This approach is particularly simple if no gluons appear in a given process at leading order. Indeed, in such cases the appearance of linear power corrections can be investigated by computing O(α s ) corrections to the process (and observable) under consideration that originate from virtual exchanges and real emissions of massive gluons, in the limit of a small gluon mass λ, see ref. [3]. 1 The presence of terms that are linear in λ implies that a particular observable receives leading power corrections that are of the type O(Λ QCD /Q), while their absence can be interpreted as an indication that non-perturbative corrections are further suppressed.

JHEP01(2022)093
In the context of high-energy collider physics, early studies of linear power corrections were mostly focused on studying shape variables in electron-positron collisions [4][5][6][7][8][9][10][11][12][13][14], 2 on the heavy-quark mass definition [16,17], on the Drell-Yan process [18][19][20] and on jets [21,22]. Recently, first attempts were made to extend such studies to more complicated processes that could be considered as proxies for realistic processes at hadron colliders. In particular, appearances of linear power corrections in top production and decay processes and in the transverse momentum distribution of the Z bosons produced in photon-hadron collisions were studied in refs. [23,24], respectively. In both cases, calculations were performed numerically for finite values of the gluon mass λ. The presence or absence of linear power corrections was established by a numerical extrapolation to vanishing values of λ.
In the case of the Z transverse momentum distribution studied in ref. [24], no evidence of linear power corrections was found. Although this result is fully sufficient for phenomenological purposes, it is interesting to understand if the presence or absence of such linear corrections in hard processes can be deduced on more general grounds. This is what we set out to do in this paper. Unfortunately, we are not yet able to perform an analysis of fully realistic processes since we have to restrict ourselves to cases where there are no gluons at leading order. Apart from this rather substantial restriction, we keep our discussion general. Whenever we are interested in a process that does contain gluons at leading order (e.g. the Z transverse momentum distribution, or e + e − event shapes in the three-jet region) we follow the approach of ref. [24] and use photons as proxies for hard gluons. Our main findings can be stated as follows: • no linear powers of λ arise from virtual corrections in generic hard processes with massless partons; • observables that are inclusive with respect to momenta of colored final state particles do not receive O(λ) and, therefore, O(Λ QCD /Q) power corrections. Observables of this type include e.g. total cross sections as well as kinematic distributions of colorless particles.
From these findings it immediately follows that no linear power corrections appear in the inclusive Drell-Yan cross section [18] and in the rapidity distribution of the Drell-Yan pair [19], at least away from the kinematic boundaries. Similarly, they also imply that the Z transverse momentum distribution computed in the simplified model of photon-hadron collisions studied in ref. [24] also does not receive linear power corrections, even if rapidity cuts on the Z boson are applied. Our results have interesting implications for event-shape studies in e + e − annihilation. Recently, this topic has received renewed attention in relation to the extractions of the strong coupling constant α s from event shapes. Indeed, it was argued in ref. [25] that a better control on non-perturbative corrections is crucial for a reliable determination of α s . 3

JHEP01(2022)093
In particular, in ref. [25] the standard approach to computing power corrections, that consists in extrapolating them from the two-jet to the three-jet region, was criticized. By considering shape variables like the C-parameter that exhibit two Sudakov regions (one near the two-jet limit and the other at the three-jet symmetric point), the authors of ref. [25] argued that the coefficient of the linear power correction near the two-jet region cannot be reliably extrapolated to the three-jet one. With our formalism, we can compute the coefficient of the power corrections in the three-jet region for several shape variables, irrespective of the presence of Sudakov regions. 4 This paper is organized as follows. In section 2 we study non perturbative corrections in a toy model, namely the production of two scalar color-charged particles in the decay of a massive vector boson. Within this toy model, we argue that the decay rate in this case is free of linear O(Λ QCD /Q) corrections. Rather than presenting new results, the purpose of this section is to illustrate basic features of our approach and to provide arguments that can be generalized to more complex cases. Such generalization is discussed in section 3, which is devoted to more complex processes with additional hard particles in both the initial and final states. There, we generalize arguments given in section 2 and argue that also in more complex cases linear O(Λ QCD /Q) terms are not present for observables that are inclusive with respect to QCD radiation.
In sections 4 and 5 we consider the implications of our result for the calculation of shape variables in e + e − annihilation. In particular, in section 4 we discuss a specific observable, namely the C-parameter, and show how our formalism can be applied to compute nonperturbative corrections to it in an approximation where the splitting of a massive gluon into a qq pair is neglected. In section 5 we present a general framework for dealing with a broader class of shape variables. Using this framework, we compute linear power corrections to both the thrust and the C-parameter distributions and compare these results with a numerical calculation at finite λ, extrapolated to λ → 0. We find consistent results, confirming our analytical findings. We conclude in section 6. This paper also contains several appendices. In appendix A we describe the computation of integrals relevant for our study. In appendix B, we detail the analytic calculation of the various integrals that we use in our analysis of the C-parameter in section 4. In appendix C, we report technical details of the calculation of shape variables in the largen f limit that we discuss in section 5. Finally, in appendix D we study non-perturbative corrections to the C-parameter in the two-jet limit.

A toy model: vector-boson decay to scalars
In this section, we consider the decay of a spin-one boson into two colored charged scalars φ. Our goal is to understand linear power corrections in this model and present arguments that, on the one hand, can be easily verified in this simple case and, on the other hand, are sufficiently general to be applicable in more complex situations. Because of this, we refrain 4 We stress however that at the present stage we are only able to obtain robust results for processes of the form e + e − → qq + γ, i.e. using photons as proxies for gluons. We will speculate on the full QCD generalization in section 5.

JHEP01(2022)093
as much as possible from using the exact form of the various matrix elements relevant for this calculation, and instead focus on their general structure.
We investigate power corrections to the process Following the discussion in the introduction, we do this by computing O(α s ) corrections to this process in a QCD-like theory where the gluon has a small mass λ, and by checking whether or not such corrections contain terms that are linear in λ. To keep our analysis as simple as possible, in this section we only consider the total decay rate. 5 We will discuss more complicated processes and observables in section 3. As stated earlier, we only consider the case of massless scalars, p 2 1 = p 2 2 = 0. We begin with the analysis of virtual corrections. There are two contributions that need to be studied -the wave-function renormalization constant for the external φ particles and the one-loop matrix element. We start with the former. We work in dimensional regularization, define the space-time dimension as d = 4 − 2 , and use the Feynman gauge for simplicity. The scalar's self-energy reads and σ 1,2 are two constants whose specific form is irrelevant for our discussion. Using Feynman parameters, B(p, λ) can be written as (

2.4)
This representation makes it apparent that which in turn implies that the wave-function renormalization constant Z φ does not contain terms that are linear in λ.
We then move to the one-loop matrix element. There are three diagrams that contribute to it. Using the Passarino-Veltman reduction [28], one can express this matrix element through four scalar integrals. Omitting color indices for simplicity, we schematically write It is well-known [3] that the total decay rate does not receive linear power corrections. However, as we have stressed in the introduction, we study this process as a first step towards establishing more general results.

JHEP01(2022)093
where B(p, λ) is given in eq. (2.3) and the other loop integrals are defined as 6 . (2.7) As indicated in eq. (2.6), the coefficient functions c, b 1,2,3 and a are rational functions of λ 2 ; this is a direct consequence of how the Passarino-Veltman reduction proceeds. It remains to consider the scalar integrals. Eq. (2.5) implies that B(p 1,2 , λ) ∼ λ −2 . Also, dimensional analysis dictates that A ∼ λ 2−2 . These integrals cannot then generate odd powers of λ upon expansion in both and λ. The case of the scalar triangle C is less trivial. However, it is easy to see that it admits the following representation Although it is straightforward to complete the integration over x and express C in terms of polylogarithmic functions, this is not necessary as the above representation makes it obvious that the function C does not contain odd terms in the small-λ expansion. We then conclude that the behavior of the renormalized one-loop virtual amplitude in the limit of small λ is described by the following formula 7 The logarithmic divergences in λ 2 are the usual soft-collinear singularities, that are canceled by analogous contributions in the real emission terms. Apart from that, eq. (2.9) implies that the dependence on λ in virtual corrections starts at O(λ 2 ). We conclude that for the process V → φφ virtual corrections do not induce any linear sensitivity to infrared physics. As the next step, we need to analyze the real-emission contribution with k 2 = λ 2 . The amplitude of this process can be written as (2.11) The two first terms on the right hand side describe emissions by external particles while the third one describes "structure-dependent" radiation. In eq. (2.11), g s is the strong coupling, T a 12 is the generator of the SU(3) color algebra in the fundamental representation and is the polarization vector of the gluon. Also, M 0 (p 1 , p 2 ) is the color-stripped matrix element with no extra emissions, while the structure of M 3 will be discussed in the following. 8 6 We note that since C(p1, p2, λ) is finite we can evaluate it directly in d = 4. 7 For brevity, we refer to all possible O(λ n ln m λ) terms as O(λ n ) terms. 8 A straightforward calculation shows that M0(p1, p2) = e(p2 −p1)µ µ V where V is the polarization vector of the decaying vector boson and e is its coupling to the scalars. Also, M µ 3 = 2 µ V .

JHEP01(2022)093
It is obvious that if the emitted gluon is resolved, the amplitude squared and the phase space can be expanded in powers of λ 2 . The situation is, however, more delicate in the soft and collinear regions where a) the amplitude develops singularities in the λ → 0 limit and b) one becomes sensitive to restrictions on the phase space induced by the gluon mass that can be linear in λ. These regions are a potential source of linear power corrections and we now study them in detail.
We begin by discussing the emission of a soft gluon. Simple power counting arguments show that this region could give rise to linear power corrections. Indeed, consider a situation where the gluon energy ω is comparable to λ, ω ∼ λ. The phase space is proportional to ω dω β θ(ω − λ), with β = 1 − λ 2 /ω 2 . Since for small λ the real emission amplitude is expandable in powers of ω starting with M ∼ 1/ω, linear terms in λ could potentially be generated. To see whether this is the case, we need to study both the matrix element and the phase space in more detail. The power counting argument implies that linear power corrections can only originate from next-to-leading terms in the small-ω expansion of both the matrix element squared and the phase space. We now discuss how to compute them.
We consider first the matrix element, and construct an expansion of the real-emission amplitude in powers of k [29,30]. We write (2.12) We can determine M µ 3 if we require that the above expression satisfies the Ward identity which means that upon replacing µ with k µ in eq. (2.12) we should get zero. We find that this is achieved if the following condition is satisfied It follows that (2.14) The amplitude expanded to first subleading order in k then reads We now need to square this amplitude and integrate it over the phase space of the three final state particles, working through next-to-leading approximation in the soft limit. Although this can be done by choosing a particular parametrization of the three-particle phase space, we do this in a way that can be generalized to more complex cases. Consider the three-particle phase space

JHEP01(2022)093
To expose its dependence on λ, we follow ref. [31] and introduce a Lorentz transformation Λ that boosts the vector q − k to the vector κq where κ = (q − k) 2 /q 2 . Specifically, We also find it convenient to definel for a generic momentum l. It is then straightforward to obtaiñ Since Λ is a Lorentz transformation, it follows that This, together with allows us to re-write the phase space as We note that in eq. (2.22) the dependencies on the gluon momentum and its mass are separated from the rest of the phase space.
We can now explicitly check whether or not soft gluons lead to contributions proportional to λ in the total rate. To this end, we need to consider 24) where N is an irrelevant normalization factor, expand the integrand through next-toleading order in small k and integrate over k. To facilitate this, we perform the Lorentz transformation discussed above. We write For soft gluons, the Lorentz transformation is small. To first order in k, it is straightforward to obtain We use q =p 1 +p 2 and rewrite eq. (2.27) as where the antisymmetric tensor B µν reads We now consider the matrix element M(p 1 , p 2 , k). From the discussion above it follows that it is sufficient to consider the approximation eq. (2.15), that we now write in terms of the momentap i . To this end, we note that and We also introduce the soft current and use eq. (2.15) to write the expansion of the matrix elements in powers of k as (2.33) Computing the square of this amplitude is straightforward. One can replace the sum over the gluon polarizations with −g µν since the Ward identity is satisfied. Contracting the soft current J µ with the various structures that appear in eq. (2.33), we obtain Using these results, we find that all derivatives of the leading order amplitude drop out from the matrix element squared and we obtain col,pol (2.35)

JHEP01(2022)093
In eq. (2.35), we need to sum over the gluon polarizations but it is not necessary to sum over the polarizations of the decaying particle.
We are now ready to ascertain whether soft gluon emission leads to linear power corrections. To do this, we study the ratio In principle, we need to integrate this formula over soft gluon momenta but it is actually easier to integrate it over all possible values of k. Since the region of large gluon momenta only gives rise to quadratic terms in λ, it is safe to extend the integration region. The computation of the integrals in eq. (2.37) is described in appendix A where it is proven that they do not contain terms which are linear in λ. We conclude that Eq. (2.38) implies that the emission of soft gluons does not generate O(λ) contributions to the total decay rate.
Since we work at next-to-leading order in the soft expansion, one may wonder whether soft scalars give rise to linear terms in λ. To study this situation, consider the matrix element in eq. (2.11). We are interesting in its behavior in the limit p 1 → 0. Momentum conservation then implies that both p 2 and k are large. Since µ k µ = 0 and since M 3 is non-singular in the p 1 → 0 limit, we conclude that the matrix element is not singular for p 1 → 0. This, together with the fact that the phase space is proportional to E 1 dE 1 , allows us to conclude that also this kinematic region cannot produce terms that are linear in λ.
The last potentially problematic region is the hard-collinear one. For definiteness we will consider the kinematic configuration where the gluon momentum is collinear to the outgoing momentum p 1 . To study this region, we employ the Sudakov decomposition to parametrize the gluon momentum k and write Then, introducing s 12 = 2(p 1 p 2 ), we write (2.40) where in the last step we have integrated over β to remove the δ-function. The z → 0, 1 limits correspond to the soft-gluon and soft-scalar cases that we have already discussed. In the hard-collinear region, k ⊥ should be integrated from zero to some value which is large compared to the gluon mass and z should be integrated between some minimal value that JHEP01(2022)093 is much larger than λ/q to z ∼ 1. Inspection of propagators shows that they are quadratic in λ. Indeed, (2.41) We conclude that the contribution of the hard-collinear region leads to an expansion in powers of λ 2 /q 2 and cannot produce terms that are linear in λ unless linear terms in k ⊥ appear.
The presence of such linear terms in k ⊥ is, in principle, possible in processes where more particles are involved. For example, suppose that there is another particle with momentum p 3 in the process. Then considering the collinear region k||p 1 and writing the Sudakov decomposition for p 3 = z 3 p 1 + β 3 p 2 + p 3⊥ , we find (2.42) Since in the hard-collinear region s 12 zβ 3 ∼ 1 and k ⊥ ∼ λ, we can expand the above propagator in powers of k ⊥ . If the integration over the directions of k ⊥ is not restricted, all odd powers of k ⊥ disappear. 9 Even powers of k ⊥ , on the other hand, correspond to even powers of λ. Hence also in this case the hard-collinear region does not give rise to linear power corrections. Finally, we note that, since collinear radiation is local in momentum space, this conclusion is general and applies to any process, regardless of its complexity. To conclude, in this section we have studied the process V → φφ(+g) and provided general arguments showing that no contributions that are linear in the gluon mass can appear. While this result is neither surprising nor new (see e.g. [3] for a generic discussion of decay rates), we avoided using explicit formulas for the matrix elements and for the phase space in the hope that the above arguments can be then easily generalized to more complex processes and observables. Indeed, the next section is devoted to the discussion of such generalization, for a broader class of processes which occur at both leptonic and hadronic colliders.

General case
We now turn to the discussion of more general cases and study processes at lepton and hadron colliders, with the usual caveat that we do not consider processes that involve gluons at leading order. 10 Although our discussion is general, for simplicity we focus on cases where only two massless color-charged partons are present in the Born amplitude (i.e. where there is only one emitting QCD dipole), while keeping the number of colorless particles arbitrary. One of our motivations for studying this case is the analysis of the Z transverse momentum distribution in photon-proton collisions in ref. [24] that we would like to understand analytically. Also, we wish to develop a general formalism that allows one to deal with non-perturbative corrections to a large class of event-shape variables. 9 We note that this cancellation is not guaranteed if observables have a non-trivial dependence on k ⊥ /| k ⊥ | in the collinear limit. We will come back on this issue when discussing event-shape variables in section 5.
10 As a consequence, the "hadron collider" may become a photon-proton collider as it is done in ref. [24].

JHEP01(2022)093
Compared to the discussion of section 2, if we consider "hadron" colliders we should also study the renormalization of parton distribution functions (PDFs). However, at least as long as the collinear factorization framework holds, PDFs renormalization is processindependent, and can be then studied in deep-inelastic scattering. But there an operator product expansion allows one to conclude on very general grounds that power corrections start at O((Λ QCD /Q) 2 ), and no linear terms are present. 11 We then only need to discuss virtual and real corrections, in analogy with what we did in section 2. We devote the next two subsections to this.

Virtual corrections
To study O(λ) contributions arising from virtual corrections, we need to consider both the wave-function renormalization constant and the one-loop amplitude. However, a simple calculation shows that even in the case of quarks the former can be expanded in powers of λ 2 and ln λ, so that it cannot give rise to linear power corrections. Hence, we only need to investigate one-loop amplitudes. 12 We have discussed virtual corrections in the toy model of the previous section where we argued that Passarino-Veltman reduction combined with explicit formulas for scalar integrals that contribute to the one-loop matrix element for the V →φφ process makes it obvious that virtual corrections possess an expansion in powers of λ 2 and logarithms of λ. To generalize this discussion, we note that the Passarino-Veltman reduction argument remains valid also for more complex processes, but the scalar integrals that one obtains are more complicated.
The analysis in the previous section was based on an explicit computation of the threepoint function C(p 1 , p 2 , λ). To understand what happens in the case of more complex integrals, it is useful to go back to that computation and ask whether an expansion of the three-point function in powers of λ can be constructed directly in the momentum representation. The answer to this question is known [32]. To obtain such an expansion, one writes the following identity 13 where the three operators produce particular Taylor expansions of various propagators. 14 11 An explicit calculation of the collinear counterterms that shows that this is indeed the case can be found e.g. in ref. [18]. 12 The situation is more delicate for heavy quarks. Indeed, in this case it is well-known that the mass renormalization counterterm can receive linear power corrections. see e.g. [3] for a review. 13 We consider C(p1, −p2, λ) rather than C(p1, p2, λ) for convenience. Indeed, C(p1, −p2, λ) is symmetric under 1 ↔ 2 exchange. 14 Although this integral is finite in four dimensions, individual expansion terms may exhibit divergences that we regularize dimensionally. In fact, it is known in this case that dimensional regularization is not sufficient to regularize T (1) k 2 separately. This subtlety is not important to us since it only affects terms that contain logarithms of λ.

JHEP01(2022)093
Specifically, It is obvious that the operator T γ λ produces terms that only contain even powers of λ. To see that this is also the case for T (1,2) k 2 , consider the j-th term in the expansion generated by T Using the Sudakov decomposition it is easy to see that upon rescaling k ⊥ → λk ⊥ and α → λ 2 α, the j-th term in the sum scales as (λ 2 ) j , modulo logarithmic corrections. We conclude that the triangle C(p 1 , −p 2 , λ) can be expanded in powers of λ 2 and no linear corrections can be generated, in agreement with the explicit result of section 2. We note that the reason why the three operators T γ λ 2 , T (i) k 2 , i = 1, 2, are needed to expand the three-point function in powers of λ is as follows. Starting from the Sudakov decomposition, it is possible to recognize [32] that only three kinematic configurations may contribute to the expansion of the three-point function C in power of λ. They are where s 12 = 2(p 1 p 2 ). These three regimes correspond to the T γ λ 2 , T k 2 operators, respectively (see ref. [32] and references therein for more details).
We continue with the discussion of more complex scalar integrals. A typical case that arises e.g. in the computation of corrections to the Z-boson transverse momentum distribution is the four-point function that, in addition to the three propagators that appear in C(p 1 , −p 2 , λ), contains a further propagator that does not go on the mass-shell in any of the singular limits (i.e. when k becomes soft or collinear to external particles). We write where The expansion of D in powers of λ proceeds in the same way as for the three-point function. We find The way the three operators act on the off-shell propagator follows from the scalings of various Sudakov parameters in relevant regions, cf. eq. (3.5). Then, the operator T γ λ does nothing to the last propagator whereas the operators T (1,2) k 2 produce its expansion in powers JHEP01(2022)093 of α, k 2 and k ⊥,µ q µ ⊥ or in powers of β, k 2 and k ⊥,µ q µ ⊥ , respectively. It follows that the operator T γ λ generates an expansion in powers of λ 2 . The action of the operator T (1) k 2 leads to integrals of the following type where we made use of the fact that upon averaging over directions of k µ ⊥ all odd powers of q ⊥,µ k µ ⊥ disappear. Also, numerator terms such as (k ⊥,µ q µ ⊥ ) 2n can be rewritten upon azimuthal integration in terms of k 2 ⊥ and q 2 ⊥ . Upon rescaling α → λ 2 α, k ⊥ → λk ⊥ , we observe that the integral in eq. (3.8) is proportional to (λ 2 ) j . The analysis of how the operator T (2) k 2 acts on the integrand leads to the same conclusion. It follows that also the box integral D(p 1 , p 2 , p 3 , λ) can be expanded in powers of λ 2 .
Using the Passarino-Veltman procedure, higher-point integrals can be reduced to boxes, triangles, bubbles and tadpoles. The latter two can be straightforwardly expanded in powers of λ 2 and ln λ, as shown in the previous section. Box and triangle integrals that do not develop infrared and collinear singularities in the λ → 0 limit can be expanded in powers of λ 2 in a straightforward way. In fact, for such integrals, the first correction to the λ → 0 limit scales as O(λ 2 ). Integrals that do develop infrared and collinear singularities, on the other hand, can be related to the box and triangle cases discussed above. We therefore conclude that virtual corrections for generic processes with massless particles do not generate linear power corrections.

Real radiation
We now discuss real corrections. Specifically, we consider a generic process I → F , where I and F are short-hand notations for the collection of initial and final state particles, respectively, and study the real-emission corrections I → F + g where g is a gluon with mass λ. We imagine that there are two and only two massless partons with QCD charges each of which can be either in the initial or in the final state. We do not consider cases when one of these partons is a gluon. On the other hand, we allow for an arbitrary number of (massless or massive) QCD-neutral particles.
From the discussion in section 2, it follows that to expose the infrared sensitivity we only need to consider singular kinematic configurations. Furthermore, in section 2 we argued that collinear emissions do not lead to linear power corrections. 15 Since collinear emissions are local in momentum space, the argument of section 2 holds for generic processes. We then only need to consider soft radiation.
When discussing the toy model of section 2, we have argued that soft scalars do not lead to linear power corrections. One may wonder if the same conclusion also holds for soft quarks. To study this, consider a quark with momentum p µ = E(1, n), where the energy E is small, E → 0. The phase space volume element of a massless quark is proportional to EdE and the most singular contribution to any matrix element squared is proportional to JHEP01(2022)093 E/(2(pk) + λ 2 ) 2 where one power of the energy in the numerator comes from the density matrix of a soft quark. Hence, the contribution of the small-energy region arises from energies that are proportional to λ 2 and is then given by where c is independent of λ. We conclude that soft quarks cannot produce contributions that are linear in λ.
As a result, we reach the conclusion that we only need to investigate the emission of soft gluons. From the discussion in section 2 it follows that it is sufficient to expand both the matrix element and the phase space in the soft limit retaining the first subleading (i.e. next-to-eikonal) terms. In the remainder of this section, we discuss how to do this for processes which are more general than the one discussed in section 2. We need to consider three distinct cases: the two QCD partons are in the final state ("final-final" dipole), one parton is in the initial state and the other one is in the final state ("initial-final" dipole) and, finally, both partons are in the initial state ("initial-initial" dipole). In what follows, we study these three cases separately. For definiteness, we will always denote the momenta of the partons that form the dipole as p 1 and p 2 , irrespective of whether they are in the initial or in the final state.

Final-final dipole
We consider a process where colorless particles with a total momentum p I produce final state particles with momenta p 1 , p 2 , p 3 , . . . , p N (3.10) Particles with momenta p 1,2 have QCD charges; 16 all other particles are colorless. We only consider cases where the QCD-charged partons are massless, p 2 1 = p 2 2 = 0. For ease of notation, in this section we will assume that also the non-QCD final-state particles are massless. To explore power corrections in this situation, we consider the emission of a massive gluon with momentum k and mass λ. Momentum conservation then reads We note that since, by construction, there is only one gluon participating in the process, the non-Abelian nature of QCD is immaterial and Ward identities are trivially satisfied. As we discussed in section 2, to study soft-gluon emission it is convenient to construct mappings of hard final state particles that preserve both the on-shell conditionsp 2 i = p 2 i , i = 1 . . . N , and the momentum conservation constraint (3.13)

JHEP01(2022)093
As we argued in section 2, since we are interested in linear power corrections, we only require these mappings to first order in the gluon momentum k. We now discuss how to construct them. Although we only need to find one particular mapping that satisfies the above requirements, we keep our discussion general as different mappings may offer different advantages and disadvantages when used in practical applications. One assumption is that the mappings behave as for small gluon momentum, where the tensors K µν i are constructed using the momentap i and the metric tensor. The momentum-conservation constraint implies We also require our mapping to satisfy the following form of the on-shell condition where λ i are analytic functions of momenta. Using eq. (3.14), we find We are now in position to express the phase-space element for final-state particles in terms of the momentap i . We write . We now make use of the fact that we only need this expression to first order in k. Then, using a relation between the determinant and the trace of a matrix that is nearly the identity matrix and the fact that λ i ∼ k, we find where To proceed further, we need to specify the mapping explicitly. To this end, we focus on the so-called dipole-local mappings, i.e. mappings where the momenta of the particles that do not belong to the radiating dipole are not transformed. By assumption, the dipole in our case is formed by the final state particles with momenta p 1,2 . Therefore, we choose JHEP01(2022)093 K i = 0 for i = 3, 4 . . . , N . Furthermore, we want to construct K µν 1,2 using onlyp 1,2 and the metric tensor. Then, writing the most general form of K µν i=1,2 (3.21) and using eq. (3.15) together with the fact that the coefficients of the tensor do not depend on k, we find the following constraints The requirement thatp i,µ K µν i k ν ∝p 2 i leads to the following set of equations (3.23) In total, eqs. (3.22), (3.23) provide 9 equations for the ten unknowns We decide to express the solution in terms of e 1 . It reads Denoting e 1 = −α, we finally obtain It is now straightforward to finalize the computation of the phase-space transformation. We find With the help of these equations, the Jacobian of the transformation J in eq. (3.20) is found to be The integration over k is restricted by the same condition that we discussed in section 2.
In particular, writing q =p 1 +p 2 and using momentum conservation p 1 + p 2 = q − k, we find that the condition (q − k) 2 > 0 puts an upper bound on the possible values of k. Before continuing, we now briefly comment on the form of the transformation eq. (3.25). First, if we compare it with the analogous transformation in section 2, we immediately see that the mapping used there corresponds to the symmetric case α = 1/2. This case is particularly simple because the phase-space Jacobian does not receive linear corrections. For the sake of generality, however, in this section we will keep α arbitrary. Second, it is interesting to note that the mapping eq. (3.25) automatically satisfies nice infrared JHEP01(2022)093 conditions. The soft limit is not particularly interesting, since for k → 0 one has p i =p i by construction. The collinear limit is however less trivial. In this case, if we formally replace k with η p 1 , we find p 1 = (1 − η)p 1 and p 2 =p 2 , which is exactly what we expect from a collinear-safe mapping. An analogous result holds for the k → η p 2 case. 17 Having studied the phase-space transformation, we need to discuss the matrix element and its integration. Since we only have one QCD dipole, the matrix element squared summed over gluon and quark polarizations can be written in the following way (3.28) The functions A, B 1,2 are polynomials in k. The limiting behavior of the function A follows from the standard soft approximation. Hence, we can write is a four-vector that, in principle, depends on all vectors p i . To understand the contributions proportional to B 1,2 , we note that they can only appear from squares of diagrams where a gluon is emitted and absorbed by the same line. Focusing on the function B 1 for the sake of definiteness, we can write where we have used µ * ν = −g µν to sum over gluon polarizations. 18 A simple computation then gives A similar calculation shows that B 2 admits an analogous decomposition. The term proportional to (p 1 + k) 2 in eq. (3.31) removes the double pole in eq. (3.28); therefore, it can be treated together with the term A/(p 1 + k) 2 /(p 2 + k) 2 in eq. (3.28). Finally, power counting arguments show that contributions of the form λ 2 /[(p i + k) 2 ] 2 are not required 17 Although for our purposes it is sufficient to consider a mapping of the form eq. (3.25), we note that in principle we could have also employed less smooth mappings. For example, assume that K µν 1,2 can be generically written as Although K µν i,⊥ could depend in a non-trivial way on (p1,2k), it is easy to see that this term leads to an odd linear dependence on the transverse momentum component of k, which vanishes after azimuthal integration. Because of this, it is then possible to show that all the arguments presented in this section would apply to this case as well, with no significant modification. 18 The sum over gluon polarizations for massive gluons contains a term kµkν /k 2 . However, this term can be dropped because of the Ward identities that are valid in the (abelian) problem.
What remains to do is to integrate the amplitude squared, expanded through first order in k, over the gluon phase space, after the p →p transformation is performed. To remap the matrix element squared, we use eq. (3.28) but we discard double poles, for the reasons we just explained. Writing the propagators as . (3.32) and expanding them to next-to-leading order in k ∼ λ, we obtain.
We now consider theoretical predictions for an observable that is inclusive with respect to QCD radiation. It follows that For such observables, we can write , (3.35) where N is an irrelevant normalization factor and the vector v µ depends on momentap i ; its exact form is not important for our purposes. We note that the upper bound on kintegration follows from the constraint θ[(q − k) 2 ]. In principle, since eq. (3.35) refers to the expansion around the soft k ∼ λ region, the integration could have been restricted accordingly. However, since our goal is to understand whether O(λ) terms appear in the differential cross section, we can extend the integration to all values of k since the region where k is hard does not generate linear O(λ) terms. A discussion of the integrals that appear in eq. (3.35) is given in appendix A where we show that they can be written as a power series in λ 2 . We conclude that arbitrary differential cross sections that are inclusive w.r.t. the QCD radiation are free of linear power corrections. On the other hand, if one computes an observable that is sensitive to gluon momenta, linear sensitivity can appear. 19 We discuss this case in details in sections 4 and 5.

Initial-final dipole
In this section, we generalize the discussion of section 3.2.1 to the case where one of the radiating partons is in the initial state and the other one is in the final state. This is relevant for example for the production of a vector boson with non-vanishing transverse momentum in hadronic collisions. At the Born level, we write where we have assigned momenta in such a way that partons with momenta p 1,2 form the dipole and p F stands for the momenta of the colorless particles. Since our formalism in its current form cannot deal with gluons in the Born process, we follow the same approach as in ref. [24] and consider quark-photon collisions We are interested in constructing a local dipole mapping that involves the partons p 1 and p 2 and that can be used to understand linear power corrections in this process. At variance with the discussion of the previous section, when constructing the mapping for the initial-state parton we require that the direction of its momentum does not change. Then, writing the transformation for p 1 as and using momentum conservation Similar to the case of final-final dipoles, we require that the on-shell conditions are not affected by the mapping. This is obviously the case for eq. (3.38) which implies for any κ 1 . The equation for p 2 2 is more informative. We find Hence, to satisfy the condition p 2 2 ∝p 2 2 , we require Since κ 1 is k-independent, it follows that

JHEP01(2022)093
In summary, for an initial-final dipole we find the following momenta mappings We note that also in this case this mapping is well-behaved in the soft and collinear limits. Indeed, by construction, in the soft k → 0 limit one has p i →p i . If we formally replace k with η p 1 , we obtain p 1 = (1 + η)p 1 , p 2 =p 2 . Similarly, for k → η p 2 we obtain p 1 =p 1 , Next, we consider the phase-space transformation. The calculation proceeds exactly as in the case of the final-final dipole in the previous section except that in the current case, we only integrate over the momentum p 2 . Since p 2 2 =p 2 2 , the parameter λ 2 from the previous section should be set to zero. Also, using the result for the Jacobian and the momentum conservation, we obtain Similarly to what we saw in the final-final case, the integral over the gluon momentum is constrained by the requirement (q − k) 2 < 0, where q =p 1 −p 2 , which is assumed in eq. (3.47). To compute hadronic cross sections, we need to convolute the partonic phase space and the matrix element squared with parton distribution functions. We write where P 1,3 are the momenta of the incoming hadrons, s hadr = 2(P 1 P 3 ) is the hadronic center-of-mass energy squared and f q,γ are the quark and photon parton distribution functions. We now interpret eq. (3.45) as a transformation rule for x 1 . Indeed, writing p 1 = x 1 P 1 andp 1 =x 1 P 1 we find through linear order in k We then use the phase space transformation and obtain (3.50) Under the assumption that x 1 is a regular point, the above equation can be expanded in ξ. Since ξ appears in the argument of the quark distribution function f q we write

JHEP01(2022)093
Similarly, the amplitude can be expanded up to next-to-eikonal level in a way that is analogous to what we have discussed in section 3.2.1. The only difference is the expansion of the two singular propagators that now read (3.52) Combining these results, we find that we again need to consider integrals that are identical to the ones for the final-final case. As we have already said, all such integrals are discussed in appendix A where it is shown that they can be expanded in powers of λ 2 . We conclude that also in this case there are no linear power corrections to kinematics distributions of final-state QCD-neutral particles. Among other things, this implies that the transverse momentum of a vector boson does not receive linear power corrections even if rapidity cuts are imposed, at least in our simplified "hadron-photon" setup.

Initial-initial dipole
In this section, we consider the case where both radiating partons are in the initial state. For concreteness, we study the Drell-Yan process Although it is well-known that the cross section of this process does not receive linear power corrections [18], we study it using our formalism for completeness. We begin by considering a suitable phase-space mapping for the process where k 2 = λ 2 . We focus on local mappings. We would like to preserve the directions of both p 1 and p 2 , so we look for mappings of the form (3.55) We note that the above mappings automatically satisfy the momentum conservation condition p 1 + p 2 − (p V + k) =p 1 +p 2 −p V and also preserve the on-shell condition for the incoming partons. Requiring that the vector boson remains on the mass shell Clearly, this equation does not have a unique solution for the two vectors κ 1,2 . However, we can require that the transformation does not change the rapidity of the vector boson Y V in the laboratory frame. Using P 1,2 to denote momenta of the colliding hadrons, we find

JHEP01(2022)093
Hence, if we choose we find It is easy to check that the choice of κ i -vectors in eq. (3.58) satisfies the on-shell conditions eq. (3.56). Once again, eq. (3.58) leads to mappings which are well-behaved in the soft and collinear limits. We now consider the phase space. Since there is no Jacobian factor in this case. However, similar to the initial-final case, we have to consider changes in the momenta of the colliding partons. We interpret them as changes in Bjorken x 1,2 . The corresponding formulas read Then, similar to the initial-final case, we have to expand the parton distribution functions in Taylor series to account for the difference betweenx 1,2 and x 1,2 . The rest of the argument proceeds in full analogy with final-final and initial-final cases. The expansion of the amplitude squared leads to integrals over the gluon momentum k that are identical to the ones discussed in appendix A, where it is shown that they do not contain O(λ) terms. This allows us to conclude that the cross section and rapidity distribution of a color singlet production at hadron colliders is free of linear power correction.

A first application to event-shape variables: the C-parameter
As we have mentioned, an interesting application of the framework developed in the previous sections is the study of non-perturbative corrections to e + e − event-shape variables, for generic kinematic configurations. In this section, we perform a semi-realistic analysis of one of such variables, the so-called C-parameter. In the case of a vector boson with momentum q decaying into N massless final state particles with momenta p 1 , . . . , p N , the C-parameter reads . (4.1) We will also use the same definition for the case of massive final-state particles. We are interested in computing power corrections to this observable in a situation that approximates a three-jet configuration in e + e − annihilations. Since our formalism does not JHEP01(2022)093 allow us to deal with processes that contain gluons at leading order, we follow the same approach as in the previous section and use photons as proxies for hard gluons. We then consider the process V (q) → q(p 1 ) +q(p 2 ) + γ(p 3 ), (4.2) and study O(Λ QCD /Q) power corrections to the C-parameter eq. (4.1) that may arise in this case. To this end, we follow the approach described in the previous sections and study O(α s ) corrections to the process in eq. (4.2) in a theory where gluons are given a small mass λ. We note that this way of computing O(Λ QCD /Q) corrections to the C-parameter is not fully justified since the definition of this observable involves momenta of all particles in the final state. For this reason, if one pursues the standard approach to power corrections that relates computations with large number of massless fermions to calculations with massive gluons, a computation of the C-parameter for the five-particle final state becomes necessary. We discuss such a computation in section 5. In this section we consider a simplified setup where we neglect gluon splitting g * → qq and consider the massive gluon as a final state particle. This allows us to directly apply ideas of the previous sections to a relatively simple but non-trivial example and to provide a connection between the general arguments of section 3 and numerical calculations of event shapes in section 5.
To compute O(α s ) corrections to the process in eq. (4.2) we need to account for virtual and real-emission contributions. We have argued in the previous section that virtual corrections cannot produce linear terms in λ; for this reason we discard them and focus only on the real-emission ones. Hence, we consider the process V (q) → q(p 1 ) +q(p 2 ) + γ(p 3 ) + g(k), To determine the non-perturbative corrections to Σ(c) we need to calculate where dσ R is the differential cross section of the process in eq. (4.3) From our discussion in the previous section it follows that we only need to study kinematic configurations where the gluon is soft. We then proceed as outlined in section 3, remap the momenta p i →p i and expand the matrix element in the soft limit retaining next-to-leading terms. To discuss modifications of the observable, we split C into two contributions C(p 1 , p 2 , p 3 , k; q) ≡ 3 + C 3 (p 1 , p 2 , p 3 ; q) + C k (p 1 , p 2 , p 3 , k; q), (4.7)

JHEP01(2022)093
where . (4.8) Then, it follows from eq. (4.7), that where v = v(p 1 ,p 2 ,p 3 ; q) is a vector whose specific form will not be needed. Finally, to account for changes in the C-parameter due to an emission of a soft gluon, we expand the θ-function to first order in k and write We now combine the required changes in the phase space, the matrix element and the observable and write In the above formula and J is the Jacobian of the transformation discussed in section 3. Using arguments presented in section 3 we conclude that the only potential source of O(λ) corrections to Σ(c) is the term C k . Since C k is proportional to the four-momentum of the soft gluon, the amplitude |M(p 1 , p 2 , p 3 , k)| 2 in the relevant terms can be taken in the leading soft approximation. We find (4.13) where (4.14) and the operator T λ is defined to extract the O(λ) contribution from the function it acts upon.
We will now explain how the function I c can be computed. To this end, we use the definition of C k in eq. (4.7) and write

JHEP01(2022)093
with . (4.16) It is convenient to compute these integrals in the rest frame of q. We find where β = 1 − λ 2 /ω 2 and ω max is an upper integration limit imposed by the condition (q − k) 2 > 0. Since we are only interested in the linear dependence on λ, the explicit form of ω max is irrelevant. Also, we have defined where n i and n are unit vectors that define the directions of the spatial components of the momentump i and of the gluon momentum k, respectively. The functions W (i) can be written as linear combinations of three integrals. They are , , Indeed, using the momentum conservation condition 3 i=1 E i n i = 0, one can write where To compute I (i) c we require the following integrals, see eq. (4.17): Their calculation is described in appendix B, where we show that Hence, the linear-λ dependence of the functions W (i) reads

JHEP01(2022)093
where (4.24) Putting everything together, we find We can use this result in eq. (4.13) to compute the O(λ) correction to Σ(c). Finally, we note that it is customary to present results for the non-perturbative corrections as a shift with respect to the perturbative differential distribution. In our case, this reads (4.26) Eq. (4.26) allows one to immediately compute O(Λ QCD /Q) corrections to the C-parameter for a generic three-jet configuration. However, as we have already said the analysis of this section is only semi-realistic since we are neglecting g * → qq splitting. We deal with the fully realistic case in the next section.

Event-shape variables: the general case
In the previous section, we explained how to simplify a computation of linear power corrections to the C-parameter in a generic three-jet configuration using an improved understanding of potential sources of O(λ) terms. However, the scope of the semi-analytic computation described there is limited since we neglected the splitting of a massive gluon into a qq pair. The goal of this section is to develop a general framework that will utilize findings reported earlier in this paper and will allow us to compute linear power corrections numerically for (almost) any shape variable in general kinematic configurations in a straightforward way.

Shape variables in the γ * → ddγ process in the large-n f approximation
To explain our approach, we consider the process γ * → ddγ. We assume that only d-quarks couple to photons and all other n f quarks only couple to gluons. In the large-n f limit, the dominant corrections arise from the emission of virtual or real gluons dressed with fermion bubbles [3], see figure 1. The solid blob in the gluon propagator in that figure implies that fermion loops have been accounted for to all orders; its exact definition follows from the recursion relation = + . (5.1) (v) (g) (qq) Figure 1. A sample of the radiative corrections that need to be included in order to compute the all-order α s (α s n f ) n corrections to the process γ * → ddγ.
As we already mentioned in the introduction, the large-n f prediction for an observable O can be obtained by computing the NLO QCD corrections to its expectation value; such computation, however, should be performed with a massive gluon. Before discussing how the findings of the previous sections allow us to easily obtain predictions for a wide class of observables, we introduce some notation. We denote the expectation value of O with O and indicate with a superscript the perturbative order at which O is computed. For example, O (0) represents the Born-level result, O (1) represents the O(α s ) correction and so on. Finally, a subscript on O indicates the number of final-state particles that are used for the calculation of the observable. To give a concrete example, we consider the cumulative distribution of the C-parameter, see eq. (4.4). In this notation, the Born-level result is denoted by where O 3 = θ C(p 1 , p 2 , p 3 ; q) − c and N is a normalization factor that we will specify shortly. The calculation of O(α s ) corrections due to a gluon with mass λ is instead denoted by with O 3 defined after eq. (5.2) and O 3+1 = θ C(p 1 , p 2 , p 3 , k; q) − c . In what follows, we choose to normalize our results to the LO rate for the process γ * → dd, i.e. 20 Finally, we introduce the short-hand notation dΦ 3 = dLips(q; p 1 , p 2 , p 3 ), dΦ 3+1 = dLips(q; p 1 , p 2 , p 3 , k), (5.5) 20 As shown earlier, radiative corrections to the total cross section do not lead to linear terms.

JHEP01(2022)093
where k denotes the massive gluon with k 2 = λ 2 and The connection between the massive-gluon calculation and non-perturbative corrections is well-known. One can write [3] where α s = α s (µ), β 0 is the first coefficient of the β-function and κ is an overall normalization which depends on the resummation prescription. Ellipses in eq. (5.7) stand for higher-order corrections, both perturbative and non-perturbative. A linear power correction is present in eq. (5.7) provided that the derivative of O (1) λ with respect to λ does not vanish for λ → 0. In what follows, we will discuss non-perturbative corrections to various quantities. However, unless stated otherwise we will always present results for λ d O (1) g * /dλ, without multiplying them by the other coefficients in eq. (5.7). So far, we have only considered final states with a massive gluon. However, if the definition of the observable O is sensitive to the presence of final-state quarks, then contributions where a massive gluon splits into a qq pair have to be accounted for. A detailed discussion of how to do this can be found e.g. in refs. [23,24], which use the same notation that we have just introduced. Here, we limit ourselves to quote the final result. To account for g * → qq splitting, one has to replace In eq. (5.11), dΦ 3+2 = dLips(q; p 1 , p 2 , p 3 , l 1 , l 2 ) is the phase space for the process and R qq is the corresponding matrix element squared The only requirement that we impose at this stage is that it should vanish in the two-jet limit where the three-jet calculation diverges. We note that eqs. (5.9) and (5.10) exhibit logarithmic singularities in the λ → 0 limit that, however, cancel in the sum. Since eq. (5.8) should be finite in that limit, eq. (5.11) should be finite as well; technically, this happens because the expression in the square bracket in that equation vanishes in the soft limit. Also, we point out that for observables inclusive with respect to g * → qq splittings, the contribution shown in eq. (5.11) would vanish and in this case the computation with massive gluons and the large-n f calculation would be identical. However, since for observables that we are interested in final states with a massive gluon never appear, the contribution shown in eq. (5.11) can be seen as a required correction to the calculation with a massive gluon, where the difference in the observable computed with the qq pair and with the massive gluon is added. In fact, as we will see in more detail below, the term proportional to O 3+(2) in eq. (5.11) and the one proportional to O 3+1 , eq. (5.10), cancel each other exactly.
A careful reader could have noticed that there are unregulated soft and collinear divergences in contributions that we account for even if we make sure to stay in the 3-jet region by, e.g. imposing a cut on the C-parameter. These singularities arise when the final state photon is collinear to one of the primary quarks or when it is soft and, instead, the radiated gluon (eventually decaying into the qq pair) is hard and at large angle. It should be clear from the results of the previous sections that kinematic regions with a hard gluon cannot produce linear power corrections and we ignore them here. A more extensive discussion of how we treat these contributions is given in the next subsection and in appendix C, where we also provide more details about the full large-n f calculation including its connection to the all-orders perturbative expansion shown in figure 1.

Simplified computation of the linear term
The results presented in the previous subsection provide a simple recipe for studying nonperturbative corrections to event shapes. Indeed, one has to perform a calculation with a massive gluon, which eventually splits into a qq pair, and then extrapolate the result of such computation to small values of λ. To do so, one requires the full matrix element for the γ * → ddγ(qq) process as well as virtual corrections to the γ * → ddγ process. For this reason, such a calculation is as complicated as any computation with a multiparticle final state can be.
In this subsection we will explain how this procedure can be dramatically simplified provided that the goal is to determine O(λ) terms only. Indeed, we will show that in order to determine linear power correction to a generic infrared-safe observable all that needs to be known is the amplitude of the Born process γ * → ddγ, the eikonal current for the emission of an off-shell soft gluon and the matrix element for its splitting into a quark-antiquark pair.
To prove this statement, we proceed as follows. First we notice that the 5-body phase space for the final state ddγqq can be factorized into a product of a 4-body phase space for JHEP01(2022)093 the production of a virtual gluon together with a ddγ final state and a 2-body phase space that describes the decay of this gluon into a qq pair. We write with dΦ split = dLips(k; l 1 , l 2 ) . (5.15) For ease of notation we did not indicate that dΦ split depends on l 1 and l 2 . Furthermore, as we discussed in the preceding sections, the 4-body phase space for the ddγg * final state can be factorized into a 3-body phase space for ddγ (the underlying Born configuration) and a radiation phase space for the gluon dLips(q; p 1 , p 2 , p 3 , k) = dLips(q;p 1 ,p 2 ,p 3 ) dΦ rad . . For convenience, we define dΦ rad to include also the Jacobian of this momenta transformation. Finally, we identify dLips(q;p 1 ,p 2 ,p 3 ) = dΦ 3 . (5.17) Using the notation introduced above, we now show that one can write O (1) λ in eq. (5.8) as follows In this equation, M µν is the amplitude squared for the production of the ddγg * final state stripped of the polarization vectors of the virtual gluon g * . Thus g * is defined in eq. (5.6). The P µν split factor in eq. (5.18) is proportional to the matrix element squared for the decay of a virtual gluon with mass λ into a qq pair with momenta l 1 and l 2 . More precisely, we define it to be P µν split = 6π λ 2 tr( / l 1 γ µ / l 2 γ ν ), (5.20) so that the following equation holds

JHEP01(2022)093
Since, on the other hand, qq , (5.22) it also follows that This equation, combined with the normalization condition for P µν split in eq. (5.21), makes it clear that the terms proportional to O 3+1 and O 3+(2) in eqs. (5.10) and (5.11) cancel out and disappear from eq. (5.18).
To proceed further, we rewrite eq. (5.18) as It is clear from the results of the previous sections that no linear power corrections can arise from the second line of the above equation, that includes virtual corrections and real emission contribution integrated over the radiation phase space. We therefore can write where the operator T λ , introduced in the previous section, extracts O(λ) terms from the expression that it acts upon. We note that the second line in eq. (5.24) has a finite λ → 0 limit since virtual and integrated real corrections are combined there. Hence, since the complete result is infrared finite, also the first line in eq. (5.24) should have a finite λ → 0 limit. This implies that the quantity upon which T λ acts in eq. (5.25) starts at O(λ 0 ) and contains higher-order terms in the λ-expansion.
In principle, eq. (5.25) already provides a method for computing linear power correction to shape variables that is much simpler than the full large-n f calculation. However, it can be simplified even further. Indeed, an important simplification arises if we observe that in eq. (5.25) the term in square brackets vanishes when k becomes collinear to the primary quarks, as long as the shape variable is infrared and collinear safe. In fact, in this limit O 3+2 becomes equal to O 3 , and the integral of P µν split becomes equal to −g µν . It seems reasonable to assume that in the hard collinear limit the left-over of the collinear cancellation does not yield terms linear in λ. This is easy to see in the thrust case, where a hard collinear splitting changes the momentum of the splitting parton by an amount that is proportional to the square of the splitting angle, and the sum of the projections of the momenta of the pair onto the thrust axis is equal to the projection of the total. We should, however, worry that this may not be the case for all shape variables. Indeed, a generic shape variable may involve terms

JHEP01(2022)093
where k ⊥ is the transverse momentum of the splitting, ϕ is its azimuthal angle, and f (ϕ) is a function that does not vanish upon azimuthal integration. In this case, hard collinear region may produce O(λ) terms. In what follows, we assume that shape variables that we consider do not give rise to such terms, i.e. that if any term linear in the absolute value of the transverse momentum does arise, it vanishes upon azimuthal integration.
Given the above clarifications about admissible shape variables, we conclude that for them linear power corrections can only arise from the emission of soft massive gluons. However, as we pointed out already, the expression in the square bracket in eq. (5.25) vanishes in the soft limit so that the full integral does not yield O(ln λ) terms. It follows then that linear terms can only arise from the leading soft-singular part of M µν . Therefore we can substitute where B is the Born matrix element, see eq. (5.6), and P µν soft (Φ 3+1 ) is the soft factor that arises from the product of eikonal currents that describe emission of a soft massive gluon in the above process. Eq. (5.25) can then be rewritten as Since the term in the square bracket vanishes in the soft limit, in principle it is not necessary to use an exact phase space to compute O(λ) terms in eq. (5.28). However, for a numerical computation it is often convenient to integrate over the exact phase space and in this case additional issues arise since the soft factor P µν soft (Φ 3+1 ) may develop unwanted singularities as we now explain. Indeed, P µν soft (Φ 3+1 ) can be expressed in terms of the original momenta 29) or the momenta after the mapping If the integration over k is restricted to soft momenta, the two equations are equivalent. However, if one uses the soft approximation outside its range of validity, spurious divergences may appear. We now discuss two examples of this, and how we deal with it. Consider the situation where p 1 becomes soft so that in the rest frame ofp 1 +p 2 (i.e. of p 1 + p 2 + k) the gluon recoils against p 2 and becomes collinear top 1 . Although this is not a singular configuration of the full process, eq. (5.30) develops a collinearp 1 ||k divergence, even if the original p 1 and k are not collinear to each other. This singularity in eq. (5.30) is spurious, and it would be removed by terms in the momentum mapping beyond the soft approximation that we are neglecting. To remedy this situation, it is sufficient to restrict JHEP01(2022)093 the integration over the radiation phase space, to exclude regions where p 1 or p 2 are soft. We do this by inserting a θ-function into dΦ rad in eq. (5.28) with 0 < η < 1. 21 We will not show this θ-function in what follows but it is always assumed to be present in dΦ rad . Similarly, care is needed to deal with kinematic regions where emitted photon is either soft or collinear to one of the primary quarks, but the gluon is hard. This region also contributes to shape variables in the three-jet regions, and its contribution is divergent. However, since in this case the gluon must be hard, no linear terms in λ can arise in this case. 22 We thus suppress this region multiplying the amplitude by the factor that dampens the photon-(anti)quark collinear singularity and approaches one if the gluon is unresolved, so that it does not affect O(λ) terms. Finally, since the integration over k in eq. (5.28) is not restricted to the soft region, there are, in principle, terms associated with hard gluons that contribute at O(λ 0 ). To remove them, we write We now note that, since the observable O is infrared safe, then This allows us to rewrite this equation as Eq. (5.34) is our final result for the calculation of linear power corrections to shape variables. Compared to a full large-n f calculation, the formula in eq. (5.34) is remarkably simple. Indeed, it only requires the knowledge of the matrix element of the Born process and eikonal factors that describe the soft emission of a massive gluon and its splitting into a qq pair. We now proceed with the discussion of how to implement eq. (5.34) in a numerical program. First, we note that the integration over radiation variables can be suitably arranged so that the cancellation among the two terms in the curly brackets occurs locally enhancing 21 In the numerical implementation we use η = 1/2. 22 Including electromagnetic virtual corrections the divergence would cancel. But again these would involve a hard gluon, and thus would not lead to any linear term.

JHEP01(2022)093
the efficiency of the numerical integration. To perform such an integration, we generate random phase-space points in the dΦ 3 dΦ rad phase space. We then compute the weight associated with the phase-space Jacobians and the product of the Born amplitude squared and the contracted soft factor P (soft) µν (−g µν ). Given the ddγg * kinematic configuration, the qq splitting kinematics is instead generated by a hit-and-miss technique, exploiting the normalization of the P µν split factor. The integration weight is then used to fill the histograms of the shape variables, computed both for the 3 + 2 and for the underlying Born phase space.
Although it is obvious that O(λ) power corrections cannot depend on the phasespace mapping, such independence provides a non-trivial check on the implementation of the numerical computation. Hence, we have used the following mappings in our computer program: (i) a mapping that preserves the direction of p 1 , i.e. such thatp 1 ∝ p 1 . This corresponds to the general mapping discussed in section 3.2.1 with α = 0; (ii) a mapping that preserves the direction of the difference p 2 − p 1 ∝ p 2 − p 1 in the dipole rest frame. This corresponds to the general mapping discussed in section 3.2.1 with α = 1/2; (iii) a mapping that preserves the thrust direction of the dipole system in the dipole rest frame. In fact, this mapping is not linear in k for small k. However, the non-linear term cancels after the azimuthal integration over k, and thus also this mapping is acceptable.
The above mappings are all of dipole-local type as defined in section 3. Besides these mappings, we have also considered the so-called global mapping of ref. [33]. All these mappings can be expanded linearly in k for small k, and we have checked that all of them give compatible results when used for the computation of O(λ) terms, as expected.
As a further check, we compare the two alternative formulae for the soft eikonal factors eqs. (5.29), (5.30) and found no significant differences.

Comparison with the result of section 4
As a first non-trivial check of the numerical approach based on eq. (5.34), we compute the C-parameter distribution neglecting the splitting of the virtual gluon into a qq pair and compare it with the result of section 4. To this end, we need to remove the g * → qq splitting from eq. (5.34). This can be done by simply replacing O 3+2 with O 3+(2) there, which corresponds to the computation of the C-parameter using momenta p 1 , p 2 , p 3 , k, where k = l 1 + l 2 , instead of p 1 , p 2 , p 3 , l 1 , l 2 . As in section 4, we adopt eq. (4.1) for the definition of the C-parameter in the massive case. We then compute δ NP defined in eq. (4.26) both from a direct numerical integration of eq. (4.26) and from a general-purpose numerical code based on eq. (5.34). In what follows, we refer to the approach based on eq. (4.26) as semi-analytic, and to the one based on eq. (5.34) as numerical.
Results for δ NP obtained with the two methods are reported in figure 2. While the semi-analytic result is linear in λ by construction, the numerical one also contains higher JHEP01(2022)093 powers of λ so that the linear term only dominates in the λ → 0 limit. This explains the differences between the numerical results obtained with different values of λ, and also the residual differences between the semi-analytic and numerical results. We also notice that, as we approach the endpoint regions c = 0 and c = 3/4, 23 subleading powers of λ become more important, thus explaining the larger differences between the semi-analytic and numerical results there. Overall, we observe good agreement between the results obtained with the two methods.

Comparison with the full large-n f calculation
As a second test of our approach, we compute coefficients of O(λ) terms for various observables using eq. (5.34) and compare them with the result of a full numerical calculation performed in the large-n f approximation. This comparison is shown in figure 3 for the differential distributions of the C-parameter and the thrust. In all cases, we perform the computation for λ = 0.5 GeV and λ = 1 GeV. For the numerical approach based on eq. (5.34), we use the mapping (ii) and eq. (5.29) for the soft amplitude in the numerical implementation of eq. (5.34). Details of the full calculation are reported in appendix C.
Since the results shown in figure 3 are divided by λ, the agreement between the λ = 0.5 GeV and λ = 1 GeV cases indicates that the dependence of the observable on λ is indeed linear and that eq. caveat, figure 3 gives strong evidence that one can use eq. (5.34) to compute linear power corrections to generic shape variables.

Non-perturbative correction as a shift in the shape variable
Having verified a simplified method for computing linear power corrections to generic shape observables, we can now use it to derive non-perturbative corrections to them. We will start with a brief overview of the history of such computations.
Non-perturbative corrections to shape variables in the two-jet limit have been considered in refs. [4, 5, 7-14, 25, 34, 35] (for a review see ref. [3]). These non-perturbative corrections are usually employed together with the perturbative ones, as well as with resummations, to extract the strong coupling constant α s from data on e + e − annihilation into hadrons [36][37][38][39][40]. Non-perturbative corrections are usually fitted in the two-jet region and then extrapolated to the three-jet region, where the value of the strong coupling constant is determined. This approach relies on the assumption that non-perturbative corrections in the three-and two-jet regions are the same.
In a recent paper [25], an attempt has been made to gain some insight into the behaviour of these power corrections away from the two-jet limit. The authors of ref. [25] studied the C-parameter distribution, that, besides the Sudakov region at c = 0, has a second Sudakov region at c = 3/4, corresponding to the symmetric three-jets configuration. The presence of this second region allows for a calculation of non-perturbative effects using techniques identical to the ones used for the two-jet region. It was found [25] that there is significant difference between power corrections in two Sudakov regions. Moreover, it was observed in ref. [25] that power corrections in the region where α s is measured strongly depend on the model used to interpolate between the two Sudakov regions. Clearly these results call for a better understanding of the dependence of non-perturbative corrections on the three-jet kinematics.

JHEP01(2022)093
In the previous sections, we have shown how to compute linear power corrections in the three-jet region in a simplified model with ddγ final state. Hence, we are in the position to compare our findings in this simplified setup with the approximate results of ref. [25]. Conversely, we should be able to reproduce the ratio of non-perturbative corrections in the three-jet symmetric point to the non-perturbative corrections in the two-jet limit obtained in ref. [25]; such a comparison should provide a further test of our numerical approach.
To set up the comparison, we follow the same approach as discussed in section 4 where it was shown that the non-perturbative corrections to a cumulant of the C-parameter can be computed as follows Although this result was derived in section 4 for a massive gluon in the final state, it is clear that it also holds if the g * → qq splitting is accounted for. The non-perturbative correction δ NP (c) defined in eq. (5.35) can be computed directly as a function of c using numerical approach described earlier in this section. However, it is customary to separate it into a normalization factor h that describes non-perturbative corrections in the two-jet region and a c-dependent function ζ qq (c) that parametrizes the dependence of non-perturbative corrections on three-jet kinematics. Hence, we write δ NP (c) = hζ qq (c).
(5. 36) In principle, the non-perturbative correction in the two-jet region can be computed in the same way as the one in the three-jet region. However, since the two-jet cross section is proportional to δ(c) in a fixed-order calculation, it is more convenient to relate h with the average value of the C-parameter computed in the two-jet region. Indeed, since dσ dc = σδ(c + h), (5.37) the average value of the C-parameter computed for two-jet events is just −h. Hence, we can write h = −T λ C (1) λ , (5.38) where because of the two-jet constraints the expectation value has to be computed starting from the Born process γ * → dd. We compute h numerically for a multitude of different values of λ; upon linear extrapolation to λ = 0, we obtain h = −9.21 (1) α s λ Q . (5.39) This result is consistent with the value 15π 2 /16 = 9.253 reported in ref. [41]; we attribute the differences between numerical and analytic results to higher powers of λ that are present in the numerical computation. 24 24 A similar computation for thrust yields 1.9457(8) as a linear slope in λ at λ = 0. We can extract analytic result for this quantity, 5π/8 = 1.9635, from refs. [5,41]. Indeed, in ref. [5] the ratio of the nonperturbative shifts to C-parameter and thrust was computed. This, together with the result of ref. [41] for the C-parameter yields the value for the thrust slope quoted above.  Having determined the normalization coefficient, we can now turn to the discussion of the function ζ qq that parametrizes the dependence of non-perturbative corrections on the C-parameter. We plot ζ qq in figure 4; these results are obtained with λ = 2 GeV and λ = 1 GeV. We note that for small values of C, ζ qq approaches unity. This is explained by the fact that soft emissions factorize independently. So, in the dominant region where both the photon and the gluon are soft, the gluon behaves as if it was radiated by a qq dipole (see appendix D).
Near the three-jet symmetric point, that corresponds to c = 3/4, we find ζ qq (3/4) = 0.226(2) for λ = 0.1 GeV. This value is consistent with the one found in ref. [25], provided that only the radiation of the quark dipole in the abelian limit is considered. We note, however, that since the normalization used here and in ref. [25] differ, we can only compare ratios of ζ-functions computed in [25]; in what follows we will always consider ζ LMS (c) = ζ(c)/ζ(0) when we quote results of ref. [25]. With this clarification, and after setting C A = 0 in eq. (18) of ref. [25], we find ζ LMS (3/4) = 0.224, consistent with our result.

Including radiation from the quark-gluon dipoles
The results described so far have been obtained for the γ * → ddγ process and not for the much more interesting case of γ * → ddg. As we explained in the introduction, this is a wellknown limitation of the large-n f approach to computing non-perturbative corrections since processes with gluons at the Born level cannot be dealt with in the theoretical framework employed in this paper.
Although we do not currently know how to overcome this limitation, the structure of the results that we obtained allows us to speculate that, perhaps, it is straightforward to do so. Indeed, our final result shows that linear power corrections to shape variables are captured by the soft approximation to the full matrix element of the γ * → ddγ + g process. For the cases that we considered so far, the soft approximation originates from a color dipole formed by the dd pair. It is tempting to speculate that for the real three-jet production process γ * → ddg we can compute linear power corrections by simply considering the emission of an additional soft massive gluon by all the QCD dipoles dd, dg anddg that are present in this case. We emphasize that we cannot prove this statement at this point, but we believe that it provides a reasonable conjecture.
Since the contributions of the three dipoles simply add up, we can write where we have exploited the fact that the dg anddg dipoles contribute equally. We have defined ζ qg (x) in the same way as ζ qq (x) discussed in the previous sub-section, except that we now assume that the radiating dipole is dg (ordg). We keep, however, the same color factor and the same normalization h used for the qq case; hence the 1/C F factors in eq. (5.40).
In figure 5 we display the function ζ qg for the C-parameter and the thrust. We observe that in both cases ζ qg approaches 0.5 for small c and t. This is easily understood, since in this limit ζ in eq. (5.40) should be one by angular ordering arguments and, since ζ qq approaches one, it follows that ζ qg approaches 0.5 (see appendix D). In the symmetric three-jet limit ζ qg approaches the same value as ζ qq . This is a consequence of the fact that in the symmetric limit the qq and qg dipoles are geometrically equivalent and, once the color factors are removed, they should give the same results.
We notice that the precision of the numerical result for the qg dipole is inferior to the qq one and also that near the symmetric point it is worse for thrust than for the Cparameter. The first issue is probably related to the fact that the hard emitting gluon is generally softer than the emitting quarks. Thus the effective Q of the emission is smaller in the qg case, leading to larger non-perturbative effects, since they are proportional to λ/Q. Regarding thrust, we recall that it vanishes in the symmetric three-jet configuration at Born level. This is different for the C-parameter, which approaches a constant there.
In figure 6 we plot ζ defined in eq. (5.40) for the C-parameter and for thrust. The results for the C-parameter can be compared to figures 1 and 3 of ref. [25]; we note again that predictions of ref. [25] need to be rescaled so that they assume the value 1 at c = 0. The normalized curves agree at the three-jet symmetric point, c = 0.75, where our result computed for λ = 0.1 is 0.479 (5), and the (re-scaled) result obtained in ref. [25] is 0.476; the difference can be attributed to terms proportional to λ 2 . Among the various extrapolations of the function ζ presented in ref. [25], their ζ b,3 curve seems to be the closest to our result.
As a final comment, we notice that for both the C-parameter and the thrust, the non-perturbative correction that we computed here is smaller than the one obtained by extrapolating it from the two-jet region to a symmetric point, especially in the case of the C-parameter. In ref. [25] a fit to α s using the C-parameter was given under various assumptions about the shape of the function ζ(c). For the function ζ b,3 that, as we said, is closest to our results, the authors of ref. [25] extract the value of the strong coupling constant α s = 0.117 (3). This result is in much better agreement with the world average value α s = 0.118(1) as compared to α s = 0.112(2) obtained in ref. [37] using a more conventional treatment of non-perturbative effects. It would be interesting to see if also for the thrust a similar improvement can be achieved.

Conclusions
Understanding non-perturbative corrections to collider processes is an interesting problem in theoretical particle physics that received surprisingly little attention in the recent past. However, thanks to the rapid development of the precision physics program at the LHC a case for a better control of non-perturbative effects in hadron collisions becomes stronger.
A possible way to investigate them is to make use of the asymptotic nature of QCD perturbation theory and estimate these effects by studying the ambiguities of a purely perturbative treatment. These ambiguities are related to the infrared pole in the running of the coupling constant. For simple enough processes they can be identified by computing O(α s ) corrections in an abelian version of QCD with a massive gluon, and extracting terms that are non-analytic in λ 2 , where λ is the gluon mass. From a phenomenological point of view, of particular importance are linear non-perturbative corrections O(Λ QCD /Q). Their presence is exposed by the appearance of O(λ) terms in the massive gluon calculation.

JHEP01(2022)093
While many explicit computations within the massive gluon framework have been performed in the past, we believe that there is a lack of general understanding of how to approach computations of O(λ) corrections to a generic scattering process or observable. In a certain sense, this is not surprising since understanding of these O(λ) terms requires a theory of soft effects at next-to-leading power which is more complicated than the familiar soft limit of scattering amplitudes and cross sections.
In this paper, we have shown that it is possible to demonstrate on rather general grounds that many terms that arise at next-to-leading power from the expansions of both phase spaces and matrix elements for typical processes and observables do not produce O(λ) terms. This result allows us to argue that certain (simplified) collider processes and observables cannot receive linear power corrections. An interesting example of this is the transverse momentum distribution of vector bosons in proton-photon collisions that even if rapidity cuts are applied does not contain O(λ) terms.
We have also shown that an improved understanding of how O(λ) terms may arise allows us to calculate non-perturbative corrections to shape variables away from Sudakov regions both analytically and numerically. To this end, we have derived a formula that allows us to compute linear power corrections to a generic shape variable in a three-jet configuration. This remarkably simple formula only involves the matrix element of the Born process and soft-radiation eikonal functions of color dipoles. As an application, we have used our formalism to compute non-perturbative corrections to the C-parameter and compared it with the result of ref. [25] which is based on an interpolation between the two-jet limit (c = 0) and the three-jet symmetric point (c = 3/4). As expected, we have found that we can reproduce the results of ref. [25] for c = 0 and c = 3/4. Between these two points, our results are close to one of the interpolations presented in ref. [25] whereas they differ significantly from a few other interpolations provided in that reference. This, of course, is not unexpected since interpolations by their very nature are subject to significant uncertainties.
Our analysis is based on the large-n f approach to the study of non-perturbative corrections; for this reason currently it cannot be applied to processes with gluons at the Born level. It would be interesting to understand how to extend this formalism to deal with these cases as well. Such an extension is, of course, very interesting for hadron collider processes. Also, as we have shown, it may lead to improvements in the description of three-jet events and to more reliable extractions of the strong coupling constant from e + e − data. We look forward to study this interesting problem in the future.

JHEP01(2022)093
which can be re-written as It follows that the integration interval for α is It is then straightforward to see that the various integrals shown in eq. (A.1) can be written in terms of the following ones This representation makes it completely obvious that the integrals in eq. (A.1) are actually functions of λ 2 .

B Computation of the I i integrals
In this appendix, we compute the integrals introduced in section 4, see eq. (4.21). The I i integrals are defined in eq. (4.19), β = 1 − λ 2 ω 2 and ω max is a kinematic bound whose precise form is not important. We are interested in the small-λ expansion of the integrals in eq. (B.1). The I 0 integral does contain linear terms in λ. Indeed We now study I 1 . A simple calculation shows that Integrating by parts, it is then straightforward to obtain which implies that for small λ the integral I 1 can be expanded in powers of λ 2 . To understand whether I 12 contains O(λ) terms in the small-λ expansion, we need to compute I 12 . To this end, we first combine the two propagators using Feynman parameters

JHEP01(2022)093
For the full calculation, we start from eq. (5.8). The required amplitudes have been analytically calculated using symbolic manipulation software MAXIMA [42]. The scalar integrals needed for the computation of the virtual corrections have been calculated with the COLLIER library [43]. The virtual contribution is infrared finite, since the gluon mass regulates infrared singularities. We have used dimensional regularization to regularize and extract ultraviolet divergences.
Integration over the phase space of final-state particles in the process γ * → ddγ diverges in the two-jet limit. To ensure that numerical computations are restricted to a three-jet region, we introduced a suppression factor where C is the C-parameter, that vanishes in the two-jet limit regulating the integral. This factor is then divided out when computing distributions and cross sections with cuts. As a result, we are able to obtain correct results as long as we do not attempt to compute observables sensitive to the two-jet region.
The real corrections are obtained by adding the emission of a massive gluon in all possible ways to the Born diagram. Integrating over its momentum, the real emission corrections are affected by collinear divergences, that arise when the photon becomes collinear to one of the primary quarks. These configurations can contribute to the three-jet cross section if the radiated gluon is hard and not collinear. Such singularities are dealt with routinely by the POWHEG-BOX framework [44], that we adopt for our calculation.
Singularities associated with soft or collinear gluons are regulated by the gluon mass λ and manifest themselves as terms proportional to log λ raised to second or first power. Similar logarithmic contributions, but with the opposite sign, arise from the virtual corrections, such that the sum of real and virtual corrections is free of log λ terms.
A carefully-constructed importance sampling near the singular regions is needed in order to reliably estimate the λ → 0 behaviour. We divide the real contribution into three regions and

JHEP01(2022)093
Here, k i and E i denote the four-momentum and the energy of the particle i. R is a shorthand notation for the R (λ) g * (Φ 3+1 ) appearing in eq. (5.10). The different contributions R (i) in eq. (C.4) correspond to different kinematic configurations of ddγg final state. For example R (1) corresponds to the region where the final state photon becomes collinear to either the d ord quark, whereas R (2) and R (3) project on regions where emitted gluon becomes collinear to d andd, respectively.
The contribution of region (1) is handled in the POWHEG-BOX [44], which implements the required subtractions of IR singularities associated with configurations containing a soft or a collinear photon. The remaining two regions are finite, but require dedicated importance sampling of the region that becomes singular in the limit λ → 0.
Finally, we compute the amplitude for the process γ * → ddγqq. This contribution is IR finite in the λ → 0 limit, but is affected by QED singularity associated with the final state photon. We thus proceed as for region (1) case, by evaluating it within the POWHEG-BOX framework. We also computed the process γ * → ddγ at NLO with a massless gluon and subtracted its result from the λ-dependent one in order to isolate the linear term.
The shape variable distributions are obtained in a standard way, by computing each contribution to sufficient accuracy so that after the cancellation of ln 2 λ, ln λ and λ 0 terms one can extract the λ dependence with enough precision.

D On the two-jet limit of C
In this appendix we elaborate more on the two-jet limit of the cumulant of shape variables within our framework. For simplicity, we focus upon the case of the C-parameter, following the calculation of section 4.
We want to show that, in the two-jet limit, the non-perturbative correction of the cumulant of C becomes proportional to the non-perturbative correction to the average value of C in the γ * → qq process. Three observations are needed to prove this: • As c approaches zero, the Born cross section has two collinear-singular regions, when the photon is collinear to either primary quark; a soft singular region, when the photon is soft; and two soft-collinear regions, when the photon is both collinear and soft.
• The correction to the C-parameter due to the emission of a soft gluon, i.e. the C k function of eq. (4.8), has a smooth limit if any pair of the 1, 2 and 3 particles become JHEP01(2022)093 collinear, as well as if one of them becomes soft. In particular, ifp 3 becomes soft or collinear to eitherp 1 orp 2 we have where we have written collectively {p 1 ,p 2 ,p 3 } = Φ 3 and {p 1 ,p 2 } = Φ 2 .
• The eikonal factor in eq. (4.14) only depends upon the direction of the radiating partons, and not upon the absolute value of their momenta.
For the qqγ final state, we only have to consider the case of emission from the quarkantiquark dipole. Our result for the non-perturbative correction, eq. (4.26) can be written concisely as and T λ is an operator that takes out the term linear in λ from the expression it is applied to. For both the collinear and the soft limits of the Born configuration, the integrand in eq. (D.3) becomes Eq. (D.6) can be immediately interpreted as the non perturbative correction to the average value of C in the two jet case. Note that c → 0 does not imply thatp 3 is either soft or collinear. We could also havep 1 orp 2 soft, or collinear to each other. What is important is that for all dominant singular contributions in the amplitude I c can be taken out of the integral. We now consider the γ * → qqg process, and use the conjecture for non-perturbative corrections that we have described in the text. In comparison to the γ * → qqγ case, we now have to consider qg dipoles as well. For the q(p 1 )g(p 3 ) dipole, the eikonal factor in eq. (D.4) becomes (p 1p3 ) (p 1 k)(p 3 k) C k (Φ 3 , k). (D.7)

JHEP01(2022)093
Whenp 2 is collinear top 3 , it reduces to as before. However, whenp 1 is collinear top 3 it becomes zero. In this case then This result follows because the most enhanced regions when c → 0 are the soft-collinear ones, but only one of the two contributes, hence the factor of one half. The above reasoning is corroborated by an explicit evaluation of the final result of section 4, eq. (4.25). Setting 25 3 (D. 10) in that equation, we find that for both the soft x 3 → 0 and the two collinear z → 0, z → 1 limits. To extend this result to the quark-gluon dipole, we still start from eq. (4.25) but assume that 1 and 2 are the quark and the gluon, and 3 is the antiquark. We now use the parametrization The singular regions of the Born amplitude are then given by x 2 → 0 and z → 0, z → 1.
In this case we find T λ I c ≈ −6π λ q 2 × ξ, (D. 13) with ξ → 2 for z → 0, ξ → 0 for z → 1, and ξ → 4z 2 − 6z + 2 for x 2 → 0. Restricting ourself to the most singular regions, i.e. the soft-collinear ones, we see that indeed only one of the two regions contributes, as we have said before. We further remark that the soft-collinear approximation is enough to get these results, that can thus be considered consequences of angular ordering. In the quark-antiquark dipole case, the result follows also from the full soft factorization that applies in abelian theories. In this case we expect that the limit is reached earlier. This is not the case for the qg dipole, since the Born level gluon and gluon emitted by the qg dipole do not factorize simultaneously in the soft limit.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. 25 We recall that xi = 2(piq)/q 2 , x1 + x2 + x3 = 2.