Transverse momentum dependent distributions in e+e− and semi-inclusive deep-inelastic scattering using jets

The extraction of transverse momentum dependent distributions (TMDs) in semi-inclusive deep inelastic scattering (SIDIS) is complicated by the presence of both initial- and final-state nonperturbative physics. We recently proposed measuring jets (in- stead of hadrons) as a solution, showing that for the Winner-Take-All jet axis the same factorization formulae valid for hadrons applied to jets of arbitrary size. This amounts to simply replacing TMD fragmentation functions by our TMD jet functions. In this paper we present the calculation of these jet functions at one loop. We obtain phenomenological results for e+e−→ dijet (Belle II, LEP) and SIDIS (HERA, EIC) with a jet, building on the arTeMiDe code. Surprisingly, we find that the limit of large jet radius describes the full R results extremely well, and we extract the two-loop jet function in this limit using Event2, allowing us to achieve N3LL accuracy. We demonstrate the perturbative convergence of our predictions and explore the kinematic dependence of the cross section. Finally, we investigate the sensitivity to nonperturbative physics, demonstrating that jets are a promising probe of proton structure.


Introduction
Since the early days of the parton model, the structure of the proton has been a major focus of the nuclear and particle physics communities.In addition to being of intrinsic interest, it is of direct relevance for describing the initial state at hadron colliders such as the LHC, and therefore important in the search for new short-distance physics.The essential theoretical ingredient is factorization, which allows one to separate the cross section into a hard scattering, that can be calculated in perturbation theory, and process-independent parton distribution functions (PDFs).The PDFs parameterize the proton structure, describing the momentum fraction of partons in the proton along the direction of motion.
We will focus on transverse momentum dependent PDFs, where in addition the transverse momentum of partons in the proton is probed.Since a transverse momentum measurement can also be thought of as the measurement of an angle, it is natural that TMD factorization theorems generically involve two TMD distributions.Traditionally, the relative transverse momentum of two hadrons in e + e − , the transverse momentum of a hadron semi-inclusive deep-inelastic scattering (SIDIS), ep → ehX, and the transverse momentum of a γ * /Z boson in pp collisions have been considered.
We recently proposed replacing individual final-state hadrons by jets in the above measurements [1].Jets are collimated sprays of hadrons, that appear in high-energy collisions because of the collinear singularity of quantum chromodynamics (QCD).Practically they are identified by clustering particles according to a specified algorithm.On the theoretical side, we demonstrated that one can simply replace the TMD fragmentation functions entering factorization theorems with our TMD jet functions, for which the use of the Winner-Take-All (WTA) recombination scheme [2] played a key role.The advantage of our approach is that such functions are perturbatively calculable, thus removing an important source of uncertainty.Specifically, the intrinsically nonperturbative distribution of the momentum fraction of individual hadrons is removed by using jets.
The advent of the electron-ion collider (EIC) will enable the extraction of PDFs with unmatched precision, with SIDIS experiments playing an important role [3].In this context, replacing the nonperturbative TMD fragmentation functions with calculable jet functions would allow one to increase the sensitivity to initial-state nonperturbative physics.Of course, for small transverse momenta, the jet functions themselves will also receive nonperturbative corrections.However, this can be addressed by exploiting the universality of the nonperturbative structure of the TMD jet function, with e + e − → dijet providing a useful testing ground.Explicitly, data from e + e − collisions could be used to fit a model for nonperturbative corrections to the jet function to be later applied to SIDIS.
A number of other jet observables that account for transverse momentum dependence have recently been considered.The main focus has been on the transverse momentum of hadrons fragmenting in jets, in both inclusive [4] and semi-inclusive [5,6] processes.In the same context, refs.[7,8] used soft-drop jet grooming [9] to reduce sensitivity to soft radiation within the jet.These studies consider the transverse momentum with respect to the standard The horizontal direction represents the beam axis.For dijets the relevant quantities q and θ are the transverse momentum decorrelation and angular decorrelation of the system, defined with respect to the relative orientation of the two jets.We consider almost back-to-back jets, θ 1, and study different hierarchies between θ and the jet radius R. In SIDIS q represents the transverse mometum of the jet, and the corresponding angle is measured with respect to the beam axis.We work in the Breit frame, where the jet recoils almost in the direction of the incoming proton, θ 1.
jet axis (SJA); instead, as an alternative way to reduce sensitivity to soft radiation, refs.[10,11] performed a similar analysis for the transverse momentum with respect to the Winner-Take-All (WTA) axis.The transverse momentum of the jet itself was also recently considered in photon + jet production [12] and lepton-jet correlation in deep-inelastic scattering [13].
Besides showing a full derivation of the results presented in ref. [1], the main purpose of this paper is performing a numerical analysis of e + e − → dijet and semi-inclusive deepinelastic scattering (SIDIS) using arTeMiDe [14,15], to study the phenomenology of TMDs with jets.
In the case of e + e − → dijet, the main physical quantity we consider is the transverse momentum decorrelation.It is defined as where p i are the jet transverse momenta measured with respect to a common direction and z i = 2E i / √ s are their energy fractions, √ s is the center-of-mass energy of the collision.Since factorization requires a small transverse momentum decorrelation, we will always assume

.2)
A related quantity is the angular decorrelation, shown in the left panel of fig. 1, where the final expression exploits eq.(1.2).This makes it explicit that we consider configurations where jets are almost back to back. 2 The angular decorrelation is similar to the azimuthal decorrelation in hadronic collisions, calculated at next-to-leading logarithmic accuracy in refs.[16][17][18][19].
In principle, the definitions in eqs.(1.1) and (1.3) depend on the choice of axis with respect to which the jet transverse momenta are measured.However, differences induced by this choice are suppressed by powers of q 2 T /s.Of course, the definition is sensitive to the details of the jet algorithm: our default throughout the paper will be the WTA axis with anti-k T [20], but we will also consider the SJA and other clustering algorithms of the k T family.We will also explore the dependence on the jet radius R.
In SIDIS, shown in the right panel of fig. 1, we choose to work in the Breit frame and define the transverse momentum as where P J is the transverse momentum of the jet with respect to the beam axis and z = 2E J /Q is the jet energy normalized to the di-lepton invariant mass Q.In analogy with eq. ( 1.3) we define a corresponding angle θ and require We use the same symbols q T and θ for analog quantities in different processes since they play the same role in factorization formulae, and their meaning should be clear from the context.To summarize our main findings: when using the WTA axis, the same factorization formulae valid for hadrons hold for jets, independently of the hierarchy between the angle θ and the jet radius parameter R. Because the factorization theorem ensures that hadronization effects in the jets are universal, they can be estimated in e + e − and then used in the analysis of SIDIS experiments.We anticipate that the main nonperturbative effects come from the evolution factor.These effects are universal (i.e. the same in e + e − , SIDIS, and Drell-Yan experiments and independent of the polarization of the hadrons) and their estimation is one of the major goals of TMD analyses.In this context we note the vital role played by the ζprescription [15], which ensures that the nonperturbative contribution to the evolution factor (that is responsible for the resummation) is uncorrelated with other nonperturbative effects.
Another observation has lead us to focus on the large radius regime of the jets.In fact, at one-loop order we notice that our jet function is well described by its large-R limit.In this limit the jet functions simplify considerably, and are determined by renormalization group evolution (RGE) up to a constant.We exploit this fact to numerically extract the two-loop, large-radius jet function from Event2 and push the accuracy of the calculation to N 3 LL in this case.Surprisingly, the validity of this regime extends down to fairly small values of the jet radius, allowing us to get precise results across the whole range in transverse momentum.This brings the perturbative precision of TMDs with jets on par with TMDs with final-state hadrons.
The paper is structured as follows: In sec. 2 we discuss the factorization formulae, considering different hierarchies between R and θ, illustrated in fig. 1.We present expressions for both the transverse momentum decorrelation in e + e − → dijet, as well as the transverse momentum of the jet in SIDIS.In addition to our default choice of using the WTA axis, we also consider the standard jet axis (SJA), for which we show that the factorization is significantly more complicated when θ R. In sec. 3 we explicitly compute the quark jet functions at one-loop order, performing the calculation in both transverse-momentum and impact-parameter space.The renormalization and resummation is discussed in sec.4, and the two-loop jet function for θ R is extracted from Event2 in sec. 5.In sec.6 we present our numerical results for e + e − → dijet and SIDIS, and we conclude in sec.7. A summary of our conventions and perturbative ingredients are collected in the appendices.

Factorization of the cross section and definition of the jet functions
The factorization of the cross section depends on the quantity For small values, R is just the jet radius parameter R, but in general the parameter R allows us to capture some power corrections.In the following we will use R when considering transverse momenta, while we use R when considering angles.In this section we review the factorization formulae of ref. [1] for all possible hierarchies, while in the remainder of this paper we concentrate on the ones that play a role in our phenomenological results.We start here by introducing the jet function, which is the main new ingredient of our analysis, providing its definition and briefly discussing its renormalization.Our factorization analysis is carried out using Soft-Collinear Effective Theory (SCET) [21][22][23][24], in which the jets are described by collinear modes and the radiation outside the jets is described by a soft mode.The typical momentum scaling of the momenta of these modes are summarized in table 1, in terms of light-cone coordinates Here n µ and nµ are light-like vectors along the directions of the jets, with n • n = 2.The jet function, that enters the factorization theorem for θ ∼ R, is written in b-space as the following collinear matrix element Here, z is the light-cone momentum fraction of the jet with respect to the initiating quark, E is the energy of the initiating quark, and P is the momentum operator.The trace in eq.(2.3) is over Dirac indices, and χ n (y) = W † n (y)ξ n (y), where ξ n is the collinear quark field in the light-like direction n µ and W n is a collinear Wilson line, ensuring collinear gauge invariance.

Mode
Table 1.The parametric scaling of the momenta (p − , p + , p) corresponding to the modes in SCET, for the various hierarchies between θ and R. For θ R the modes differ between the Winner-Take-All and standard jet axis.
Gluon-initiated jets do not enter for e + e − and SIDIS, but we give the corresponding definition for completeness, where is the collinear gluon field, with F αβ n the collinear field strength tensor.We will also perform the calculation in momentum space, which simply involves replacing The above definitions are for the bare jet functions, as indicated by the absence of renormalization scales.A perturbative calculation shows that both ultraviolet (UV) and rapidity divergences affect these distributions, so that one should consider the renormalized quantities and similarly for J g .Here Z q is the UV renormalization factor, R q is the rapidity renormalization factor, and rapidity divergences are removed first, as in ref. [25].A key observation is that these renormalization factors are the same as in the case of TMDs, as we discuss in sec.4.

R ∼ θ 1
We now turn to the factorization analysis, starting with dijet production in e + e − scattering at a center-of-mass energy √ s, where θ ≈ 2q T / √ s ∼ R 1.This is the simplest case since there are only two scales, √ s and q T .The cross section differential in the momentum decorrelation q and the jet energy fractions z i = 2E J,i / √ s factorizes as3 The hard function H e + e − encodes the hard scattering process, in which a quark-anti-quark pair is produced.It contains virtual corrections, but no real radiation because that would result in q T ∼ √ s.For convenience we have extracted the tree-level cross section σ e + e − 0 , which contains a sum over quark flavors.The jet functions describe the fraction z i of energy of the initial (anti)-quark that goes into the jet, as well as their transverse momentum through the impact parameter b (its Fourier conjugate).They depend on the jet algorithm, as indicated by the argument √ s 2 R, but this does not affect its anomalous dimension, as required by RG consistency.Soft radiation does not resolve the jet because its typical angle is order 1, whereas R 1. Consequently, we do not have to consider clustering soft radiation in the jet algorithm, and we can simply include its effect as an overall recoil of the system, as indicated in eq.(2.8).The soft function has been absorbed into the jet functions in the above expression, as we will discuss in sec.4.There we will also show that the RG evolution between the hard scale µ H ∼ √ s and jet scale µ J ∼ q T in eq.(2.8) resums invariant mass logarithms of µ H /µ J ∼ √ s/q T , and similarly that ζ is related to the resummation of invariant rapidity logarithms of √ s/q T [14,15], see also refs.[26][27][28][29][30].
The corresponding factorization theorem for the cross section of semi-inclusive deepinelastic scattering is given by dσ ep→eJX dQ 2 dx dz dq = which is differential in the di-lepton invariant mass Q 2 , Bjorken x, the energy fraction z of the jet generated by the splitting of the quark, and the jet transverse momentum q T .We work in the Breit frame, where z = 2E J /Q, and apply an e + e − jet algorithm.The modification to the factorization theorem compared to eq. (2.8) is fairly modest: the hard function is replaced by the one for SIDIS, one of the jet functions is replaced by a TMD PDF, and the sum over quark flavors must be explicitly included because both σ DIS 0,q and F q depend on it (J q does not, since we only consider light quarks, which we treat as massless).The hard function is slightly different, where C V is the Wilson coefficient for the hard matching, l Q 2 = ln(µ 2 /Q 2 ) and a s = g 2 /(4π) 2 .
The NNLO and NNNLO expression can be found in ref. [31], taking into account that H e + e − (Q 2 , µ) is the same as for the Drell-Yan process.The two loop expressions are provided in eqs.(A.19, A.20) of the appendix, to make the paper self-contained.

R θ 1
We now consider the case where we have an additional hierarchy due to the small size of the jet radius, R θ 1.While this regime will be of limited phenomenological interest to us, it allows us to make contact between our framework and TMD measurements with final state hadrons, corresponding to the R → 0 limit.The modes are again listed in table 1, and involve additional collinear modes whose scaling is set by R.
The factorization in this case is an extension of eqs.(2.8) and (2.9).The jet function contains two scales √ sR q T , which can be separated through a further collinear factorization, (2.11) Only collinear radiation at angular scales θ, encoded in C i→j , can affect q T .However, subsequent splittings down to angles of order R will change the parton j with momentum fraction z into a jet with momentum fraction z.This is described by the semi-inclusive jet function J j , which has been calculated to O(α s ) in refs.[32,33] (our notation matches that of ref. [32]).The distinction between WTA vs. standard jet axis is irrelevant, since θ R. The additional RG evolution between µ J ∼ q T and µ J ∼ ER sums single logarithms of The (z ) 2 in front of C i→j was chosen to ensure that these matching coefficients coincide with those for TMD fragmentation, given to O(α 2 s ) in refs.[25,34].That these same matching coefficients enter here is is not surprising, since for R → 0 the semi-inclusive jet function becomes the fragmentation function (summed over hadron species) [35].Thus in this limit we reproduce the known results for TMD fragmentation to hadrons.For convenience, we collect the relevant one-loop expressions for the matching coefficients and semi-inclusive jet function in eqs.(A.31) and (A.32) of the appendix.

θ R for the Winner-Take-All axis
We now consider θ R for the Winner-Take-All axis.For R ∼ 1, the modes in table 1 are expected.The factorization takes on a rather simple form, since soft radiation does not affect the position of the jet axis, due to the WTA recombination scheme.In particular, there is no distinction between soft radiation inside and outside the jet.The only effect of soft radiation is its total recoil, which is therefore described by the standard TMD soft function.Since θ R, the collinear modes do not resolve the jet boundary, so z = 1 and the ER dependence drops out, (2.12) For completeness we also provide a definition of and a similar formula can be written for the gluon case.For θ R 1, one would expect the same modes as are listed for the standard jet axis in table 1.In this case the soft function does not resolve the jet boundary, because R 1, but collinear-soft modes with scaling resolve the jet boundary and contribute to q T .However, by the same reasoning as before, their only effect is a total recoil on the system, independent of whether emissions are inside or outside the jet.Consequently, these additional modes do not need to be considered, since they will simply be removed by the zero-bin subtraction [36], due to their overlap with the soft mode.This leads to the interesting conclusion that, for the WTA axis, the cross section for θ R is independent of R.

θ R for the standard jet axis
For completeness we also discuss θ R for the standard jet axis.We do not present any numerical results for this case, and therefore limit our discussion to the dijet momentum decorrelation in e + e − collisions.First we consider the case θ R ∼ 1, for which the modes are given in table 1. Energetic emissions outside the jet are not allowed because these would lead to θ ∼ R. Because the standard jet axis is along the total momentum of the jet, momentum conservation implies that q T is simply determined by the transverse momentum of soft radiation outside the jets.In particular, the angle of energetic emissions inside the jet is unrestricted.Since R ∼ 1, these emissions are hard, explaining the absence of a collinear mode.Each of these hard emissions induces a soft Wilson line, implying the presence of nonglobal logarithms (NGLs) [37] of √ sR/q T .The corresponding cross section can be described using the framework of refs.[38,39] (see also refs.[40,41]) (2.15) We have eliminated the measurement of the momentum fractions of the jets, since z i = 1 in this limit.H m denotes the hard function with m real emissions inside the jets, along the light-like directions n i .The soft function S m describes the transverse momentum q T of soft radiation outside the jets, produced by the Wilson lines along the directions n i .The color indices describing the representation of the hard emissions/Wilson lines connects the hard and soft function, and Tr c denotes the trace over these color indices.Finally, ⊗ denotes integrals over the light-like directions n i .
Moving on to θ R 1, we have collinear modes whose angular size is set by R, and additional collinear-soft modes with scaling which are fixed by the requirement that they resolve the jet boundary and contribute to q T .Because R 1, no hard real emissions are allowed, and the soft function does not resolve the jet.However, each collinear emission produces a collinear-soft Wilson line, in direct analogy to the soft Wilson lines generated by hard emissions for R ∼ 1.Using again the framework of refs.[38,39], the corresponding cross section is given by dσ (2.17) The hard and soft function are the same as for θ ∼ R. The jet function J m describes m collinear emissions inside a jet along light-like directions n i , and the collinear-soft function U m describes the resulting q T from collinear-soft emissions of these Wilson lines.

Quark jet function at one loop
In this section we present a detailed calculation of the one-loop quark jet function that enters the factorization formula in eqs.(2.8) and (2.9).We use dimensional regularization with d = 4 − 2ε to handle UV divergences, and the modified δ-regulator for the rapidity divergences [34,42], The Feynman diagrams and measurement are discussed in sec.3.1.We present a detailed calculation in momentum space in sec.3.2 and in impact-parameter space in sec.3.3, thus providing a cross check of our results.The advantage of performing the calculation in momentum space is that this is the space in which the jet algorithm is defined.
On the other hand, the renormalization and resummation are simpler in impact-parameter space.Here ⊗ represents the collinear quark field χ n , which contains a collinear Wilson line that can emit gluons.A sum over cuts is understood, where cuts through loops describe real emissions, while only cutting the quark line corresponds to virtual corrections.The latter vanish for our choice of regulators.

Feynman diagrams and measurement
The one-loop diagrams that contribute to the quark function are given in fig. 2. This leads to the following expression for the bare jet function up to one loop, Here E is the energy of the quark field initiating the jet, and its small light-cone component + (and thus virtuality) is integrated over.The phase-space of the outgoing gluon, with momentum k µ , and quark, with momentum µ − k µ is integrated over.The δ + (k 2 ) ≡ δ(k 2 )θ(k 0 ) and δ + [(k − ) 2 ] denote the corresponding on-shell conditions.The coupling has been replaced by the renormalized one in the MS scheme, leading to the prefactor (. . . ) ε .There are three different cases we need to consider: (a) both partons are inside the jet, (b) the gluon is outside the jet, (c) the quark is outside the jet.
These cases are identified by Θ case , and the transverse momentum q 2 T case and jet energy E J,case depend on the case and jet algorithm, and are given in table 2 in terms of the energy fraction of the quark and of the jet size R, defined in eq.(2.1).
Table 2.The Θ case encodes the various regions of phase space, and the corresponding jet energy E J case and transverse momentum q 2 T case .At this order the only difference between jet algorithms is the recombination scheme, i.e. standard jet axis vs. Winner-Take-All.
At one loop, there are only two partons, so every distance measure gives the same clustering condition (as we will see in sec.5, this is no longer true at two loops).There are differences between the standard and Winner-Take-All (WTA) recombination scheme.This distinction is only relevant when both partons are inside the jet, in which case the standard jet axis is along their total momentum while the WTA axis is along the most energetic one.
Switching from k − to the quark energy fraction x, using the on-shell conditions, and exploiting azimuthal symmetry, we rewrite the one-loop term of eq.(3.1) as Here we replaced the δ − regulator by its dimensionless counterpart After similar manipulations, the corresponding one-loop gluon jet function is From this expression one can obtain the one-loop result for the gluon jet function presented in ref. [1], following step by step the calculation of the quark function detailed below.

One loop results in momentum space
In order to perform the calculation in transverse momentum space we directly solve the two integrals in eq.(3.3), inserting the measurements for the various cases in table 2. We start with the case of both partons inside the jet.
In the case of the standard jet axis, the dependence on the transverse momentum is trivial and the calculation reduces to the one performed in ref. [32] for the semi-inclusive quark jet function.After integration over the transverse momentum, Here we set the rapidity regulator δ − E to zero because the endpoint x = 1 is already regulated by dimensional regularization.The remaining integral over the energy fraction is a combination of Euler Beta functions, whose expansion up to O(ε 0 ) yields where For the WTA axis, the transverse momentum dependence becomes nontrivial.The condition min(x, 1 − x) reduces to x > 1  2 if we symmetrize the integrand, (3.9) Performing the remaining integral requires to treat the integrand as a two-dimensional distribution, see eq. (A.10), and yields the result Finally we consider the cases where only one particle is inside the jet, that are independent of the jet algorithm.We use x → 1 − x to combine the case where the gluon is outside the jet, with the case where the quark is outside.Both the integrals over transverse momentum and energy fraction are fixed by the δ functions enforcing the measurement, resulting in Expanding the result in ε and δ − E requires again some algebra with distributions, that is performed explicitly in app.A.2.We obtain We now combine the expressions in eqs.(3.7) and (3.10) with (3.12), to obtain the bare quark jet function at one loop The dependence on the algorithm occurs via the functions ∆ axis q , that explicitly read The expression for the WTA axis is more involved because it introduces the threshold z > 1 2 .We notice that This implies that the dependence on the jet algorithm vanishes in the regime R θ, as predicted from the factorization formula in eq.(2.11) (the semi-inclusive jet function J that enters there is independent of the jet axis).

One-loop results in impact-parameter space
The calculation of the quark jet function at one loop can also directly be performed in impactparameter space.This calculation provides a check of the results in the previous section.We perform the same two integrals of eq.(3.3) with the cases shown in the table 2 as in the momentum-space calculation, but first carry out the Fourier transform of the jet function J alg [1]   q (z, b, ER) = dq e ib•q J alg [1]   q (z, q, ER). (3.17) The case with both partons inside the jet is the only one that depends on the choice of axis.
The result for SJA has a trivial dependence on the transverse momentum and can be written as Note that for this calculation the IR divergences are regulated by ε and we can safely neglect the δ − E regulator.
The WTA axis choice introduces a non trivial dependence on the transverse momentum of the jet function.Symmetrizing the integral over x, as in eq.(3.9), we rewrite the jet function as Integration over the transverse momentum allows us to rewrite eq.(3.19) as The jet function depends on the transverse position in terms of the dimensionless combination The remaining step is the integration over x.The integral in the first term (first two lines) is straightforward to perform analytically.On the other hand, the second integral has a part for which we were unable to obtain a closed analytical expression.The result of this second integral is given by the function G(B ER ), whose explicit expression is given in eq.(A.28).This leads to Note that the only difference between the SJA and WTA results is G.When B ER 1 the function G is zero, as required by the axis independence in this limit.
Next we consider the case when only one parton is inside the jet.By using x → 1 − x, we can combine case (b) and (c).As we now have an explicit dependence on the momentum fraction of the jet, the rapidity regulator δ − E needs to be kept.We find, The terms which a divergent behavior in the limit z → 1 should be understood as regulated under the +-prescription.For clarity we have split the δ(1 − z) contribution into three pieces: the first term will be eliminated after the renormalization of rapidity divergences, and the third term is exactly cancelled by the corresponding part of the case with both particles inside the jet, removing IR divergences presented here as double poles in ε.
The final result for the quark jet function for both choices of axis is obtained summing eq.(3.18) (SJA) or (3.22) (WTA) with (3.23), where We have checked that these expressions agree with those obtained in sec.3.2, which is partially numerical for the WTA axis.For the numerical implementation, we find it most convenient to take a numerical Fourier transform of the momentum space result, rather than using the above expressions.

Rapidity renormalization
The jet function in eq. ( 2.3) has the same renormalization as in the case of TMDs.Here we summarize the main points of rapidity renormalization, referring to e.g.ref. [25] for further details.The rapidity renormalization factor R q in eq.(2.7) can be extracted from the soft function taking into account the zero-bin Zb q , that accounts for the overlap with collinear modes.The soft function for SIDIS is given by the following vacuum matrix element of soft Wilson lines where the coordinates in brackets indicate the position of both Wilson lines, and T ( T ) denotes (anti-)time ordering.The Wilson lines are defined as usual In eq. ( 4.2) we did not include the transverse gauge links which are necessary to preserve the gauge invariance in singular gauges [43][44][45], because we do not use them in our computations.For e + e − all Wilson lines are future pointing, which corresponds to Sn → S n in eq. ( 4.2), but this fact has no practical consequences as in the case of TMDs [28,[46][47][48].The overlap of collinear and soft modes depends in general on the rapidity regulator used in the perturbative calculation and for the modified δ-regulator used in the present work one finds that in fact S q (b, ζ, µ) = Zb q , so that the rapidity renormalization factor has the simple form R q (b, ζ, µ) = 1/ S q (b, ζ, µ).The parameter ζ in R q is a scale that comes from splitting the soft function in two factors, each one of which is absorbed in one of the jet functions or TMDs.Specifically, the perturbative calculation of the soft function reveals that it depends linearly on ln(µ 2 /(δ + δ − )), where the δ ± are the rapidity regulators for each of the collinear modes in the factorization theorem.To separate them, ζ ± are introduced with 4 , and 2E is the hard scale of the process under consideration.In the calculation of a jet function along the direction n one can effectively replace δ − E = δ + E ζ, so that the subscripts ± for the variable ζ can be omitted.While the rapidity renormalization factor is simply multiplicative in b-space, the jet function can also be calculated in momentum space, as we have shown in the previous section.

One-loop renormalization of the jet function and small and large R limits
Our bare jet function in eqs.(3.13) and (3.24) is still affected by divergencies.As discussed in eq.(2.7) and sec.4.1, its renormalization is particularly easy to implement in impactparameter space, where it is purely multiplicative and takes the same form as for hadron TMDs.The explicit one-loop UV and rapidity renormalization factors are Z [1]  q (ζ, µ) = − ) leading to the renormalized expression The corresponding momentum-space result is presented in eq.(A.24).From eq. (4.7) (or equivalently from eq. (A.24)) one can take the limits R → 0 and R → ∞, to approach the factorization regimes described respectively in sec.2.2 and in secs.2.3 and 2.4.In the small-R limit the two axes give the same result, and we explicitly checked that the jet function factorizes further as in eq.(2.11).The perturbative ingredients in which the jet function factorizes are listed in app.A.4.The large-R limit is particularly interesting for the WTA axis, where the jet function simplifies as in eq.(2.12).We verified that the dependence on the jet radius drops out in this limit, obtaining

Resummation and ζ-prescription
The renormalization group equations (RGEs) of the TMD jet function are the same as for the standard hadronic TMD, where D q and γ q are the rapidity and UV anomalous dimension, respectively.We only consider the quark jet function, because the gluon does not enter in our phenomenological results.As in the hadronic TMD case we have where | f.p. denotes the finite parts.
Since the order of derivatives can be interchanged, one obtains [29,49], where Γ cusp q is the quark cusp anomalous dimension.Consequently, where and γ V is the finite part of the renormalization of the vector form factor.Both γ V and D are known up to O(a 3 s ) [50][51][52][53][54], and a numerical computation of the fourth-order cusp anomalous dimension was recently presented in ref. [55].All these anomalous dimensions are collected in app.A.3.
The high energy scale value for µ is always set at the hard scale, i.e. √ s for e + e − and Q for SIDIS.As for the TMD case, the evolution of the jet function in the plane (µ, ζ) is governed by eq.(4.9).A systematic treatment of this case has been provided in ref. [15], and in our results we have implemented the optimal solution suggested in that work.Summarizing the main points: The solution of eq.(4.9) is in principle path independent, when the anomalous dimensions are known to all orders.This means that the evolution is effectively provided by an evolution potential: in the plane (µ, ζ) one can identify null-evolution curves corresponding to equipotential lines and the only true evolution occurs only between jet functions belonging to different equipotential lines.When the perturbative expansion of the anomalous dimension is truncated, it is possible to recover a path independent result through e.g. the improved-γ scheme of ref. [15], which only affects terms in the perturbative expansion beyond the order that one is working at.At this point we are left to choose an initial equipotential line ζ µ (b), which is known as the ζ-prescription.A special line is provided by the saddle point of the evolution potential.This line exists for all values of b (at least for b T < 1/Λ QCD ) and covers all the ranges on µ and ζ, providing the optimal solution Explicitly, at two-loop order The evolution of the optimal distribution to a generic set of scales (µ, ζ) is then simply given by where (µ 0 , ζ µ 0 (b)) is a point on the special line and U q R is the TMD evolution factor Choosing the simplest possible line which connects the initial and final point of the evolution in the improved-γ scheme, eq. ( 4.17) reduces to4 which is convenient for numerical calculations.The rapidity anomalous dimension D q has a nonperturbative part, which is independent of other nonperturbative inputs of the jet distribution and should be estimated by itself.The ζ-prescription (unlike e.g. the b * -prescription) allows this separation theoretically.At the moment, the only extraction of the nonperturbative part of the evolution factor from data within this prescription has been carried out in ref. [56], so that in our phenomenological analysis we use their parametrization for the nonperturbative contribution to the rapidity anomalous dimension, Here D res q is the resummed perturbative part of D q , and where the constants B NP and c 0 parametrize the nonperturbative effects.The perturbative expansion of the resummed rapidity anomalous dimension D res q is where X = β 0 a s (µ) ln(µ 2 b 2 T e 2γ E /4), β 0 is the leading coefficients of the QCD beta function and a s = g 2 /(4π) 2 .The leading term reads and we have used this expansion up to third order in a s , which incorporates the four-loop anomalous dimension.The complete expression up this order can be found in ref. [15,57].The unresummed expression for the rapidity anomalous dimension is reported in eq.(A.16).

Numerical implementation of evolution
We use arTeMiDe to run the double scale evolution from the initial scale of the TMD jet function/PDF where µ 0 is frozen at 2 GeV to avoid the Landau pole and (µ 0 , ζ 0 ) belongs to the special line, to the hard scale Since the rapidity resummation is the dominant source of uncertainty and to consistently use the nonperturbative parameters extracted in ref. [56], we will always use the highest known order in the evolution, even though the jet function for generic R is only calculated at one-loop order.The nonperturbative parameters of the evolution kernel in eqs.(4.19) and (4.20) are set to 5 Quark jet function for large R at two loops As we will see in our numerical analysis, the large-R limit captures the dominant part of the perturbative corrections.This justifies focussing on the quark jet function in the large-R limit, J WTA q , which is completely determined at two loops by known anomalous dimensions, except for a constant j [2] .Explicitly, We extract this constant using the Event2 generator [58], which we run with n f = 5 and an infrared cutoff ρ = 10 −12 , generating about a trillion events.Specifically, we consider the difference at O(a 2 s ) between the the cross section with a cut on the angular decorrelation θ ≤ θ cut obtained from Event2 and our factorization theorem, extracting the overall factor of a 2 s = α 2 s /(4π) 2 .This is shown in fig.3, where the different panels correspond to the (e + e − version of) anti-k T [20], Cambridge/Aachen [59,60] and k T [61] jet algorithm.The different curves in each panel correspond to the C 2 F , C F C A and C F T F color structure, with the bands indicating the statistical uncertainty.From varying the infrared cutoff we conclude that the cross section obtained from Event2 can be trusted for log 10 θ cut > −3, corresponding to the plotted range.
The clear plateau at small values for θ cut shows that our factorization theorem predicts the singular part of the cross section correctly.The value of the plateau corresponds to the missing two-loop constant j [2] (the overall factor of 1/2 was chosen to cancel the factor of 2 from the two jet functions in the factorization theorem).The decomposition of j [2] in terms of the C 2 F , C F C A and C F T F color structures is given by i.e. the color structures are inside the constants.The result was extracted for n f = 5, and the generalization to arbitrary number of flavors only involves rescaling j [2] T F .While the power corrections, corresponding to contributions to the cross section not included our factorization theorem, decrease as θ cut becomes smaller, the statistical uncertainties increase.This motivates us to fit the plateau to a constant, considering as fit range −3 ≤ log 10 θ cut ≤ log 10 θ cut max , where we vary log 10 θ cut max between −2.9 and −2.We take the spread between the values obtained in this manner as an estimate of our uncertainty, finding anti-k T : j T F = −13.0± 0.3 . ( While these constants are remarkably similar for anti-k T and Cambridge/Aachen, they differ substantially for k T .

Results
The region of interest for TMDs is small q T , for which the regimes θ ∼ R and θ R are most relevant.This leads us to exclusively focus on the WTA axis, which is well behaved in the large-R limit.We start by considering the transverse momentum decorrelation in e + e − collisions, obtaining numerical predictions for the Belle II and LEP experiments.We use e + e − to test the perturbative convergence, and explore the dependence on the jet radius R and cut on the jet energy fraction z.In the case of SIDIS we provide numerical predictions for the EIC, and investigate the sensitivity of our cross section to nonperturbative effects.
In our numerical implementation we build on the arTeMiDe code [14,15] to obtain resummed predictions for TMD cross sections.The original version of arTeMiDe [62] provides cross sections for Drell-Yan and SIDIS with fragmentation into hadrons.However, its modular structure allowed us to extend it to processes involving jets with a modest amount of modification.Specifically, we have added e + e − → dijet and jet-SIDIS high-level modules, and a jet TMD low-level module that provides our perturbative input for the quark jet functions in b-space at the initial scale.

Momentum decorrelation in e + e − collisions
In our analysis of the e + e − cross section, differential in the transverse momentum decorrelation, we consider two experiments: • Belle II: √ s = 10.52 GeV, 4 quark flavors.
We account for both the photon and Z-boson contribution, and restrict the plotted q T range to a region where the power corrections to the factorization theorem can be neglected.In the Belle analysis we omit b-jets, since we do not include quark mass effects in our calculation of the jet function.(Experimentally, these are of course relatively easy to distinguish from light quark jets.)We start our analysis by studying the dependence on the jet radius parameter R in fig. 4 for LEP.The cross section is shown for various jet radii, ranging from R = 0.1 to 0.7, using the factorization formulae for θ ∼ R in sec.2.1.We consider two representative cuts on the jet energy fraction: z > 0.25 (left panel) and z > 0.75 (right panel).For comparison we also show the large-R limit, discussed in sec.2.3.We use the one-loop jet function (since we only have the one-loop result for θ ∼ R), but include the hard function at two-loop order and perform the resummation at N 3 LL accuracy.
As expected, as R increases the results approach the R → ∞ limit.In both cases, the cross section for R = 0.7 is indistinguishable from the large-R result, and for z > 0.25 the difference is even minimal for R = 0.5.This observation will be used in the rest of our analysis, to justify including the two-loop jet function in the large-R-limit, as this will capture the dominant two-loop contribution.Explicitly, we will combine results according to dσ dq T where NLO and NNLO indicate the order of the jet function.In each term we use the NNLO hard function and include the resummation at N 3 LL accuracy.Therefore, our approximate N 3 LL result returns the exact two-loop expression in the large-R limit, but will miss some (small) O(θ/R) corrections.
Next we study the perturbative convergence of the TMD cross section in fig. 5. We take R = 0.5, z > 0.25 and show results for the cross section for Belle II (left panel) and LEP (right panel) at NLL, NNLL and N 3 LL.The ingredients that enter in the various perturbative The N 3 LL result is obtained with the prescription in eq. ( 6.1).The bands encode the perturbative uncertainty, as described in the text.Figure 6.Dependence of the transverse momentum decorrelation distribution on the cut on jet energy fraction z, for Belle II with R = 0.7 (left panel) and LEP with R = 0.3 (right panel).The dependence on this cut is larger for smaller R, as discussed in the text.In both cases, the results for z > 0.5 (solid red curve) exactly coincide with the large-R limit, see footnote.
orders are summarized in table 3. The perturbative uncertainty is estimated by varying the scales µ i in eqs.(4.23) and (4.24) up and down by a factor 2 around their central value and taking the envelope.The band obtained by this procedure at LL is artificially small and not shown.As expected, the N 3 LL correction is small compared to the NNLL one, and the uncertainty bands overlap and are reduced at higher order.
In fig.6 we investigate the dependence of the cross section on the cut on the jet energy fraction z > z cut for a fixed value of the jet radius, which provides a complementary picture to fig. .Estimate of the sensitivity of the TMD to nonperturbative effects in the rapidity resummation at Belle II (left) and LEP (right).We vary the parameter c 0 in the range of its statistical uncertainty, testing both the fixed and variable B N P schemes of ref. [56].
functions.For R = 0.7 the dependence on the cut on z is relatively mild, which reflects the fact that in the large-R limit the jet function is proportional to δ(1−z), and thus independent of this cut.For R = 0.3 there is a stronger dependence.This is not surprising, since the cross section diverges as z cut → 0 and has large logarithms of 1 − z cut for z cut → 1.We found that, regardless of the jet radius, for z cut = 0.5 the cross section coincides with the large-R result.This is due to a one-loop accident. 5s a next step, we study how sensitive these cross sections are to B N P and c 0 that parametrize the nonperturbative contribution to the rapidity evolution, see eqs.(4.19) and (4.20).We considered both the "fixed B N P " and "variable B N P " schemes used in the recent fit in ref. [56], and varied the parameters within the statistical errors listed in their table 4. In practice, we found that the B N P variation is subdominant, so in fig.7 we simply plot variations of c 0 .As one would expect, the sensitivity to nonperturbative effects is much larger at Belle, commensurate with its smaller center-of-mass energy, and increases at low transverse momenta.The conclusions obtained from the two schemes are compatible with each other.The situation is similar for LEP, though the relative variation is substantially lower (below 1% for most of the range in q T ).
Finally, we have investigated the impact of the choice of jet algorithm, specifically, the impact of the different two-loop constants in eq.(5.3).We found the difference with respect to anti-k T to be negligible for Cambridge-Aachen (< 0.1%) and very small for the k T algorithm (< 1%).

Transverse momentum dependent distributions in SIDIS
In this section we show representative results for TMD measurements with jets in SIDIS.
We focus on the EIC, which is a future facility for the study of TMD distributions, and consider some representative values for the experimental kinematic parameters.Specifically, we assume a center-of-mass energy of √ s = 100 GeV.We take 10 ≤ Q ≤ 25 GeV and study the transverse momentum distribution for q T ≤ 3 GeV, ensuring that power corrections of order q 2 T /Q 2 to the factorization theorem can be neglected.In this kinematic range we expect quark mass effects to be negligible, so we ignore them.We work in the Breit frame, impose a cut on the jet energy fraction z > 0.25 and set the jet radius to R = 0.5.Our e + e − analysis in the left panel of fig. 4 shows that in this case the large-R approximation works extremely well, so we again include the two-loop, large-R jet function of sec.5, using eq.(6.1).
We use the quark TMD PDFs obtained in ref. [56].In this fit the matching of the TMDs onto PDF is incorporated at NNLO, using the NNPDF 3.1 PDFs [63] with α s (M Z ) = 0.118.The additional nonperturbative component of the TMD PDFs is modeled with the ansatz where the values for λ i were fit in ref. [56].
Our results are shown in fig.8, for which we consider different intervals in the elasticity y in the range 0.01 < y < 0.95.In each case, we obtained the uncertainty band by independently varying the scales µ H and µ 0 up and down by a factor of 2 around their central values, and taking the envelope.We find that roughly half of the contribution to the cross section comes +5% 5% 0 Figure 9.The sensitivity of the cross section to nonperturbative effects is estimated by varying the parameter c 0 , that controls the nonperturbative contribution to the evolution kernel, within its current statistical uncertainty [56].from low elasticity (y < 0.2).The variation in shape between the different elasticity intervals is modest; at high elasticity the peak of the distribution is shifted towards larger transverse momenta.
The EIC will shed light on the nonperturbative structure of TMDs, which motivates us to investigate the sensitivity of our observable to nonperturbative hadronic physics.A rough impression of the sensitivity of our cross section to nonperturbative corrections can be obtained by varying the parameters B N P , c 0 and λ i (see eqs. (4.19), (4.20) and (6.2)) that enter our nonperturbative model.In principle, these parameters are highly correlated and a full error estimate would require taking data with a large number of replicas, along the lines of the original analysis in ref. [56].In practice, we observe that the nonperturbative uncertainty is dominated by the variation of the single parameter c 0 .Therefore, we obtain a realistic estimate of the size of NP effects by simply varying c 0 within its statistical uncertainty, which we show in fig.9.The effect of varying c 0 is not large (below 5% in the whole range), but nonnegligible, and grows for small q T .This plot suggests that this measurement can likely be used to improve our knowledge of the nonperturbative part of the evolution kernel, parametrized by c 0 , which is very relevant because it is universal.We have explored the dependence on R, z cut and the range in Q and y, finding similar sensitivity to nonperturbative effects.

Conclusions
The study of the transverse momentum distribution of the proton can benefit from using jets (instead of hadrons) as final state.A clear advantage is that the jet momentum can be calculated in perturbation theory, while the fragmentation of hadrons is an intrinsically nonperturbative process.We provided an initial formulation of this idea, using a modern definition of jets, in ref. [1].There we observed, for the first time, that the cross section for dijet production in e + e − collisions and SIDIS with a jet in the final state can have the same factorization as for hadronic TMD measurements, simply replacing a TMD fragmentation function by our TMD jet function.This factorization depends on the jet radius R and recombination scheme, holding only for all values of R if the Winner-Take-All axis is used.In particular, in the regime of small q T , which is interesting for extracting the intrinsic transverse momentum of partons in the proton, the cross section for the standard jet axis does not satisfy the usual TMD factorization.
To explore the ramifications of these ideas, we presented numerical results in this paper for Belle II and LEP (e + e − collisions), and the EIC (SIDIS), building on the existing arTeMiDe code.We reported the details of the NLO calculations of the TMD jet function, and have also numerically evaluated the NNLO contribution in the large-radius limit with Event2.This was motivated by the observation that the NLO result is well described using the large-R jet function, for all experimental cases we consider.Consequently we can achieve the same N 3 LL accuracy as in the corresponding hadronic TMD cases.
We have verified the perturbative convergence of our numerical predictions, achieving perturbative uncertainties of order 5% in the peak of the distribution at N 3 LL.We also find that our cross sections have similar sensitivity to nonperturbative effects as the corresponding hadronic case, without the burden of additional nonperturbative effects from fragmentation.Specifically, we have investigated how the cross section changes when varying the nonperturbative parameters within the errors provided in ref. [56], concluding that in principle these experiments can provide important constraints on these parameters.Here we benefit from using the ζ-prescription, which ensures that the nonperturbative parts of the evolution kernel and the rest of the TMD are uncorrelated.
The nonperturbative effects to the jet TMD have not been estimated in this work.However our factorization theorems ensures that these effects can be included in the definition of the jet functions and are therefore universal, i.e. the same in e + e − collisions and SIDIS.In this respect, the hadronization of jets can be treated in the same way as the nonperturbative part of a hadron TMD, and is therefore expected to be subdominant compared to the nonperturbative part of the evolution.Consequently, jet measurements may provide one of the best ways to constrain the nonperturbative part of the evolution kernel.To reduce the sensitivity to hadronization effects one can consider grooming, which will be investigated in a forthcoming publication [64].
Related "cut" distributions are defined as L cut n (q T , q 0 ) = L n (q T , q 0 )θ(q 0 − q T ) .(A.8) Plus distributions naturally arise in the expansion of logarithmically-singular terms in dimensional regularization, T ) + L 0 (q T , µ) + O(ε) .(A.9) In the calculation of the jet function in momentum space, one encounters terms where the above expansion cannot be used because the divergences in the limits ε → 0 and δ − E → 0 are mixed by a step function.In particular, expanding eq.(3.11) involves the identity where the last term involves a genuine two-dimensional distribution.This identity was obtained by switching to cumulative distributions in both variables, then expanding in δ − E , and finally expanding in ε, Every term can now be identified as the cumulative of a distribution, resulting in eq.(A.10).We note that the last term in the expansion, defined by is only well-defined by the prescription in eq.(A.5) if one first carries out the integral over q 2 T before the integral over z.
Table 3. Various orders in resummed perturbation theory, and the fixed-order (F.O.) and resummation ingredients they involve.The fixed-order ingredients are the perturbative expansion of the hard function, jet function and the coefficients in the matching of the TMD PDFs onto collinear PDFs.We also use the PDFs extracted at this order as well, and use the corresponding running of the coupling.

A.3 Anomalous dimensions
We now list the anomalous dimensions that enter the double-scale evolution described in sec.4. Our predictions use N 3 LL resummation by default, corresponding to the first row in table 3.An exception is fig.5, where we compare different orders to test the convergence of resummed perturbation theory.We only need the anomalous dimensions for quarks, which we expand as a n s D [n]  q . (A.13) The coefficients in the expansion of the cusp anomalous dimension are given by Γ [0] q = 4C F , Γ [1]  q = C F C A 268 9 − 4π 2 3 − 80 9 n f T F , Γ [2]  q = C F C The fourth-order result was computed numerically in ref. [55], which also provides a full decomposition in terms of color structures.The non-cusp anomalous dimension is given by, V,q = C F C F − 3 + 4π The rapidity anomalous dimension can be conveniently expressed in terms of Γ cusp q in eq.(A.14) as D [1]  q = Γ [0] q 2 L µ , L µ + D [2]  q (0) , q β 1 L 2 µ + 2β 0 D [1]  q (0) + Γ [2]   q L µ + D [3]  q (0) .(A.16) The first two coefficients of the QCD beta function, that enter here, are given by and the constant terms read D [3]  q (0) = C F C

Figure 1 .
Figure1.Geometry of the event for e + e − → dijet (left) and SIDIS (right).The horizontal direction represents the beam axis.For dijets the relevant quantities q and θ are the transverse momentum decorrelation and angular decorrelation of the system, defined with respect to the relative orientation of the two jets.We consider almost back-to-back jets, θ 1, and study different hierarchies between θ and the jet radius R. In SIDIS q represents the transverse mometum of the jet, and the corresponding angle is measured with respect to the beam axis.We work in the Breit frame, where the jet recoils almost in the direction of the incoming proton, θ 1.

Figure 2 .
Figure 2. Cut diagrams that contribute to the one-loop quark jet function in SCET.Here ⊗ represents the collinear quark field χ n , which contains a collinear Wilson line that can emit gluons.A sum over cuts is understood, where cuts through loops describe real emissions, while only cutting the quark line corresponds to virtual corrections.The latter vanish for our choice of regulators.

Figure 3 .
Figure 3.The difference between the O(α 2 s ) contribution to e + e − cross section with a cut on the angular decorrelation θ ≤ θ cut , obtained from Event2 and from our factorization theorem.The panels correspond to the (e + e − version of) anti-k T , Cambridge/Aachen and k T jet algorithm, and the curves correspond to the different color structures, see eq. (5.2).The uncertainty bands indicate the statistical uncertainty.The missing two-loop constant in the quark jet function is the value of the plateau at small θ cut .

Figure 4 .
Figure 4. Dependence of the cross section differential in the transverse momentum decorrelation on the jet radius parameter R, for cuts on jet energy fraction z > 0.25 (left) and z > 0.75 (right).We use the NLO jet function computed in the regime R ∼ θ, and show the large-R result (red solid) for comparison.

Figure 5 .
Figure 5. Perturbative convergence of the cross section differential in transverse momentum decorrelation, for Belle II (left) and LEP (right), for jet radius R = 0.5 and jet energy fraction z > 0.25.The N 3 LL result is obtained with the prescription in eq.(6.1).The bands encode the perturbative uncertainty, as described in the text.

4 .Figure 7
Figure 7. Estimate of the sensitivity of the TMD to nonperturbative effects in the rapidity resummation at Belle II (left) and LEP (right).We vary the parameter c 0 in the range of its statistical uncertainty, testing both the fixed and variable B N P schemes of ref.[56].

Figure 8 .
Figure 8. TMD cross section for SIDIS with jets with √ s = 100 GeV, 10 < Q < 25 GeV and different intervals in elasticity within the range 0.01 < y < 0.95.