Transverse Momentum Measurements with Jets at Next-to-Leading Power

In view of the increasing precision of theoretical calculations and experimental measurements, power corrections to transverse-momentum-dependent observables are highly important. We study the next-to-leading power corrections for transverse momentum measurements in $e^+ e^- \rightarrow 2$ jets. We obtain a factorized expression for the cross section, which involve twist-2 and twist-3 operators, and identify the new jet functions that appear in it. We calculate these jet functions at order {\alpha}s for a family of recoil-free schemes, and provide the corresponding anomalous dimensions at leading order. Additionally, we show that the (endpoint) divergences that typically arise in sub-leading-power factorization can be subtracted and cancel for our case. By working with jets, everything is perturbatively calculable and there are substantial simplifications compared to the general next-to-leading power framework. Importantly, our analysis with jets can be extended to semi-inclusive deep-inelastic scattering, with the future Electron-Ion Collider as key application.


Introduction
The momentum distributions of quarks and gluons (collectively called partons) inside hadrons is one of the most interesting subjects in Quantum Chromodynamics (QCD).These distributions are crucial for any prediction when hadrons are collided, such as the Large Hadron Collider (LHC) and the Electron Ion Collider (EIC).The fundamental theoretical tool for accessing these distributions is provided by factorization, which allows one to express (factorize) a cross section in terms of calculable coefficients and nonperturbative parton distributions.The universality of these distributions ensures that they can be measured in one process and then applied to another one.
Parton distributions are correlation functions of partons located at two different space-time points in a hadron.If these two points lie on a light-cone direction, one gets the usual parton distribution functions (PDFs).When there is also a transverse separation, the correlation functions are transverse-momentum-dependent distributions (TMDs).In principle one can also consider more general structures, such as Wigner distributions, generalized TMDs (GTMDs), etc., depending on more kinematic variables and involving different hadronic states, but they will not be the subject of this work.Furthermore, one can also consider the different functions that account for the spin of partons or/and hadrons, so as to obtain the most general information about partonic entanglement [1].There are also corresponding nonperturbative distributions that describe (the fragmentation of) hadrons in the final state from energetic partons.We will focus on TMDs in the final state and work with jets rather than hadrons, because they can be described in perturbation theory.
The factorization of the cross section for TMDs has a long history, starting with the seminal works of refs.[2][3][4].It was observed that perturbative calculations involve a new type of divergences, nowadays called rapidity divergences.These divergences need special treatment and point to a more involved form of factorization than for the pure collinear case involving PDFs.TMD factorization has been achieved in more recent times, using special regulators for rapidity divergences and effective field theory [5][6][7] or the background field method [8].In these works, TMD factorization results from expanding the cross section in powers of λ, which is the ratio of a small transverse momentum q T and the scale Q of the hard scattering process.In Drell-Yan (DY), semi-inclusive deep-inelastic scattering (SIDIS) and e + e − → 2 jets (or 2 hadrons), q T is the small transverse momentum of the di-lepton pair and Q their invariant mass.
Restricting to the leading power in this expansion, and performing perturbative calculations for the renormalization group energy evolution, one has arrived at the extraction of unpolarized TMDs up to next-to-next-to-next-to-next-to-leading logarithmic (N 4 LL) order [9,10].This is currently the highest logarithmic order that has been used to determine a nonperturbative distribution.The TMD formalism has also been used to make a theoretical prediction for the Z-boson transverse momentum distribution at the same order [11].In principle, the same accuracy can also be achieved for the other, spin-dependent distributions described by TMD factorization, because the TMD evolution kernel [12][13][14][15][16], cusp anomalous dimension [17,18] and hard factors [19][20][21][22] don't depend on spin.
Given this high theoretical precision for the cross section at leading power, there have recently been advancements in the exploration of hadron structure at higher powers of λ.These power corrections are expected to introduce significant insights into our understanding of hadrons and their properties [8,[23][24][25][26][27][28].The need for power corrections is also clear from the per mille-level accuracy of recent DY measurements at the LHC.Additionally, SIDIS and e + e − experiments (like EIC [29], EIcC [30] and Belle) are typically conducted at low energy, where power corrections are typically larger [31].Summarizing, power corrections will enable higher precision in many experimental situations.The process that we consider here can be important to establish the impact of NLP effects.
Recent progress has opened the possibility to study power corrections within the TMD formalism.In particular, we currently have at our disposal a complete basis of operators up next-toleading power (NLP) [8,24], whose evolution properties have been studied [8,27].These studies have demonstrated that the factorization of the DY and SIDIS processes remains valid up to NLP precision.However, additional non-perturbative operator matrix elements enter in the cross section, indicating the presence of new non-perturbative physics.Furthermore, the factorized formula for the cross section contains additional divergences characteristic of sub-leading-power factorization, which require careful treatment.
We find it desirable to have access to cross sections where non-perturbative QCD effects are minimal, to enable direct tests of NLP operators and shed light on the formalism.To this end, we investigate NLP effects using jets as hadronic final states.This is because jets are infraredsafe quantities whose properties can be largely fixed perturbatively, i.e. they can be calculated in perturbation theory 1 unlike the T-odd jets in refs.[32,33] that require a nonperturbative hypothesis.Specifically, we are considering the impact of NLP corrections on di-jet production at e + e − colliders.This offers the possibility to test NLP factorization at e + e − colliders such as LEP and Belle, and its extension to SIDIS will be particularly useful at the EIC [34].
The possibility to use jets in TMD factorization has been explored in several works, see f.i.[35][36][37][38][39][40][41][42][43][44][45][46].Refs.[38,39] focused mainly on the consistency of the TMD factorization theorems with jet definitions and algorithms, finding that the definition of the jet axis and radius are essential to establish the factorization.The standard jet definition suffers from non-global logarithms, due to its sensitivity to whether soft radiation gets clustered into the jet or not, substantially limiting the accuracy.By contrast, recoil-insensitive schemes [47] such as the Winner-Takes-All (WTA) jet axis [48] are insensitive to this effect, enabling high precision calculations.
In this work, we address two fundamental questions: First, we want to identify the new jet functions that arise from the additional operators at NLP.Second, we want to obtain the factorized expression for the cross section in terms of jet functions and hard matching coefficients.The answers to these questions that we obtain, enable us to study how NLP corrections impact transversemomentum-dependent measurements, including angular asymmetries.We find that in order to get a non-vanishing contribution from NLP operators, one has to break the symmetry between the energetic radiation going into the two different directions.For the case of di-jet production, this can be achieved by e.g. using two different (recoil-free) recombination schemes.
Despite the complexity of the subject we have made an effort to present our material in a didactic way.The choice of the process and the jet definitions that we use, greatly reduces the number of independent distributions and simplifies the evolution with respect to the original papers [8,28,49,50], see eq. (4.24).In this sense, the present work can also be viewed as an introductory overview to the study of NLP effects in TMD factorization.
The outline of this paper is as follows: We start in sec. 2 by introducing the kinematic variables needed to describe e + e − → 2 jets, and writing the cross section in terms of the hadronic and leptonic tensor.In sec.3, the hadronic tensor is factorized and rewritten, with Fierz identities relegated to app. A. Our final expression for the cross section in terms of jet functions is given in eq.(3.36).In sec.4, we discuss the cancellation of rapidity divergences and additional divergences that show up at sub-leading-power.This section also covers the renormalization and evolution of the jet functions.We calculate and present all ingredients of the cross section to first order in a s in sec. 5. Our conclusions are given in sec.6.

Kinematics
The process that we consider in this paper is the annihilation of an electron and positron, with momenta ℓ and ℓ ′ , into two jets with momenta P 1 and P 2 , e − (ℓ) e + (ℓ ′ ) → j 1 (P 1 ) j 2 (P 2 ) X . (2.1) We will assume a large jet radius, such that all collinear final-state particles are clustered in one of the two jets, and any radiation outside the jets X is soft2 .At lowest order in the electromagnetic coupling, this process proceeds by the electron and positron annihilating into an intermediate photon with momentum q = ℓ + ℓ ′ , which subsequently creates a quark-antiquark pair that produces two jets.In the following, we assume that the jet algorithm yields massless jet momenta, P 2 1 = P 2 2 = 0 (as in the WTA recombination scheme), and we introduce light-like vectors n and n such that In order to introduce our observables, we define two planes, which we will refer to as the transverse and perpendicular plane.The transverse plane is orthogonal to the two outgoing jets and the perpendicular plane is orthogonal to the intermediate photon momentum q and one of the ℓ ℓ′ outgoing jets.These planes can be described in terms of the metric tensors that project onto the corresponding sub-spaces, and we define the transverse and perpendicular components of a vector v as Note that, by construction, We will sometimes write two-dimensional vectors in boldface notation, where it is implied that the two-dimensional metric has positive signature.As we mostly encounter transverse vectors, we will suppress the subscript T , e.g.b = b T .We now list the variables that we use to describe the cross section, which are the same as in ref. [51] for the massless case: as well as an azimuthal angle The angles ϕ 1 and α are not independent.It is natural to use the former in the center-of-mass frame and as depicted in fig. 1, while it is convenient to use the latter in the frame where the jets are exactly back-to-back.They are related by All dot products of four-vectors can be expressed in terms of the observables listed above.
While the above expressions are Lorentz invariant, the interpretation we now give is valid for the (di-lepton) center-of-mass frame: First of all, the z i (i = 1, 2) are the energy fractions of the two jets.The variable y corresponds to the angle θ 2 between the momentum of the jet and the lepton, y = (1 − cos θ 2 )/2.Next, ϕ 1 is the azimuthal angle between the perpendicular component of jet 1 and the electron momentum ℓ.Finally, where θ 12 is the angular decorrelation of the two jets.These angles are shown in fig. 1.Note that the invariant mass Q of the collision is not an observable, but it is determined by the initial state of the experiment.
In this paper we consider the large jet radius limit, where z 1 = z 2 = 1 up to power corrections, so we will not be differential in them.This leads us to consider the following cross section where L µν and W µν are respectively the leptonic and hadronic tensor.The leptonic tensor reads and in this work we will neglect the lepton helicity contribution proportional to λ e .Our main focus is thus on the hadronic tensor, which contains all QCD effects.For this, the recent work of ref. [8] opens the door to a more systematic study of the hadronic tensor up to NLP, and this is the subject of the next section.The case of two hadrons (instead of jets) has been studied before, especially in relation to Belle experiment, see the reviews in refs.[31,[52][53][54].The distributions in this case were classified in ref. [51], and mass corrections to kinematics have been included in ref. [55].

Factorization at NLP
The goal of this section is to obtain a factorized expression for the cross section in terms of a set of jet functions and hard functions.The expression we obtain holds at the bare and unsubtracted level, and we postpone the discussion of renormalization and (overlap) subtraction for the next section.
To obtain a factorized expression for the cross section we first expand the hadronic tensor up to next-to-leading power using the TMD operator product expansion method.The complete and un-expanded hadronic tensor for the process under consideration is given by, where J µ (y) is the electromagnetic current Since we assume a large jet radius, all energetic final-state particles are clustered into one of the two jets and X only contains soft radiation.The explicit notation that we use for the power expansion of the hadronic tensor is In order to arrive at our final expressions for W LP,NLP , we proceed in several steps, starting from the TMD operator expansion method and the basis of operators in sec.3.1.In particular, we list a set of operator building blocks that have definite twist, which allow for a clean separation in terms of power counting.We use these definitions in sec.3.2, to obtain the power expansion of the product of currents in eq.(3.1) in position space (indicated by the tilde) Next, we include the external jet states and Fourier transform back to momentum space, to obtain an expression for the hadronic tensor in terms of a set of jet functions with open spinor indices.In sec.3.3, we apply Fierz relations to obtain a set of jet functions with closed spinor indices, and use discrete symmetries and Lorentz decomposition to simplify and reduce the basis of jet functions, resulting in just one new jet function at NLP! Sec.3.4 contains the final expression for the hadronic tensor in terms of jet functions, as well as its contraction with the leptonic tensor.

Operator basis
In this work we make use of the TMD operator product expansion [8].In this method, background fields are introduced for the two collinear directions The background fields q n , q n, A µ n , A µ n will be used to construct our operators, while the fields χ, B µ are purely dynamical.As usual, we decompose the quark background fields into "good" and "bad" components, where such that each field has a definite power counting (ξ ∼ λ, η ∼ λ 2 ).The collinear and anti-collinear regions overlap in the soft region.To avoid the double counting of this region its contribution must be carefully subtracted.We discuss the subtraction of the overlap region in sec.4.1.
Using the background fields one can construct a set of building-block operators of definite twist (dimension minus spin), Here, the first subscript on the operator denotes the twist of the operator and the second subscript denotes the collinear sector to which the operator belongs.These operators contain Wilson lines, which we denote by The gauge field appearing in the Wilson line matches the collinear direction under consideration, e.g.A n (A n) for the Wilson line in U 2,n (U 2,n ).The operators in eqs.(3.8) and (3.9) will enter in the matrix-element definition of the jet functions.Note that these operators are process dependent via their dependence on L and L, with L = −∞ for incoming partons and L = +∞ for outgoing partons.Since we consider jet production, we set L = L = +∞ from now on.In the Wilson lines we will also introduce a δ-regulator whenever necessary (see refs.[8,56] for more details).

Operator basis at LP
At lowest order in the power expansion of the product of currents in eq.(3.4), for which the result will be shown in sec.3.2, the following combinations of operators appear The (±) superscript denote causal and anti-causal fields from the Keldysh formalism [57].The color indices of the fields are contracted, which is made explicit for O 11 by the color trace tr c , and i, j are spin indices.The gauge invariance in the expressions above is restored using the building-block operators in eqs.(3.8) and (3.9), by (3.12)

Operator basis at NLP
At next-to-leading power in the expansion of eq.(3.4), one finds contributions from the following operators To write these in terms of the building-block operators we need to do a few manipulations.First we replace the transverse gauge fields that appear here by a field strength via the following relations, which introduces a δ-regulator for rapidity divergences.Subsequently, we apply the following identity where a positive δ ± E = δ ± /q ± ensures that the integral is well defined and shows how δ regulates rapidity divergences.Note that in this identity the sign of ξ may be chosen freely as it is integrated over all values.It turns out that it is convenient to choose +ξ for the O 21 operators and −ξ for the O 12 operators, as the jet functions that will be built from these operators then have support for ξ ∈ [0, 1].This convention differs from earlier work, and allows for several simplifications in our case.
The above manipulations allow us to express the operators in eq.(3.13) in the NLP hadronic tensor in terms of the building-block operators in eqs.(3.8) and (3.9) as follows, We have omitted the space-time arguments on the operators on the left-hand side as they are the same as in eq.(3.13).The sign in front of iδ differs between O 12 and O 21 .As we will see, their matrix elements are related by complex conjugation.

Expansion of the hadronic tensor
The aim of this section is to write the hadronic tensor in terms of matrix elements of the buildingblock operators in sec.3.1.For this, we first consider the power expansion of the product of currents in eq.(3.4) in position space.Following [8], the leading power contribution can be written at leading perturbative order as For brevity we suppressed the space-time arguments of the operators, but they are the same as in eq.(3.12).At next-to-leading power more operators appear and the expression for the hadronic tensor is a bit more involved, where again we have considered the leading perturbative order and have suppressed the spacetime arguments of the operators for brevity.The non-trivial Wilson coefficients that enter beyond leading order in perturbation theory involve a convolution in the light-cone positions, so we prefer to include them only when switching to momentum space in eq.(3.19).As explained in more detail in ref. [8], two different types of power corrections can be distinguished in eq.(3.17), those involving derivatives of twist-2 operators (kinematic power corrections) and genuine twist-3 contributions.
At this point we are ready to write the hadronic tensor in momentum space, including the Wilson coefficients and external states, allowing us to pass from operators to jet functions.To switch to momentum-space, we note that the inverse derivatives that appear in eq.(3.17) can be treated using the following identity, Since the background fields of the different collinear sectors are not coupled at this order in the power counting, the matrix elements factorize and the result for the hadronic tensor can be written in terms of products of jet functions, Here, C 1 and C 2 are matching coefficients that can be calculated in perturbation theory.These coefficients arise when the current operator is written in terms of the background fields, and we follow the definitions of [8], and the symbol J denotes jet functions with open spinor indices.At leading order in perturbation theory C i = 1, consistent with eqs.(3.16) and (3.17).The coefficient C 1 is known up to three-loops [21] (see also [18,58]), while the coefficient C 2 is calculated up to NLO in [8], section 6.For definiteness, the C 2 in eq.(3.19) is related to the one in ref. [8] by The bounds of the ξ integral arise from the range over which the jet functions have support.
The jet functions that appear above are defined in terms of matrix elements of the building-block operators in eqs.(3.8) and (3.9).The twist-2 jet functions are given by and the twist-3 jet functions are defined as and we have introduced for convenience Note that whether the operator appears in the amplitude or conjugate amplitude is determined by whether the operator is causal or anti-causal.From now on we will drop these superscripts on the operators.

Jet function definitions
The next step in the factorization is to define a minimal set of (bare) jet functions that have no open spinor or Lorentz indices, as the set of jet functions in sec.3.2 contains redundancy.First, to get rid off the open spinor indices, we apply the Fierz relations in app. A. This gives rise to four twist-2 jet functions with no Lorentz index and eight twist-3 jet functions with one open Lorentz index, The factor 2N c accounts for averaging over the spin and color of the initial quark field.The different overall sign between J 21,n and J 21,n is chosen because ϵ µν T = ϵ µνρσ n ρ nσ changes sign when n ↔ n.Discrete symmetries provide an additional reduction of independent jet functions.Jet algorithms are spin-independent and respect the C,P and T symmetries of QCD.This allows us to derive that the LP jet functions satisfy, reducing the number of independent twist-2 jet functions down to two, one for n and one for the n direction.For the twist-3 jet functions we have and similarly for J q 12 .This reduces the number of independent twist-3 jet functions from eight down to two.Since the quark and anti-quark distributions are identical we will drop the corresponding label from now on.The minus signs appearing in the complex conjugation cancel when performing the decomposition in eq.(3.28).Thus the signs in eq.(3.25) were chosen precisely such that these symmetry relations do not contain minus signs.
Both the NLP jet functions and the derivatives of the LP jet functions that appear in the hadronic tensor carry an open (transverse) Lorentz index.To achieve a fully factorized formula for the hadronic tensor we must decompose these objects into the available structures b µ and ϵ µν T b ν .Parity implies that only b µ contributes, leading us to define This final step removes all open Lorentz indices from the jet definitions.Since all functions now only depend on the magnitude of b, we can perform the angular integration using the following identities, Applying the above decomposition to the definitions presented in eqs.(3.21) and (3.22), the final expressions for the jet functions read, where the trace is over both the color and Dirac indices.The jet functions for the n-collinear sector can be obtained by replacing n ↔ n.

Factorization
Starting from eq. (3.19), we now apply the relations and definitions of the previous section to obtain a factorized expression for the hadronic tensor, Here, H 1 and H 2 are defined in terms of the Wilson coefficients and are referred to as the hard functions There are two things to note here.First, note that the hard function H 2 is complex valued for the process we consider here.Second, we have replaced the arguments 2q + q − of the Wilson coefficients by Q 2 , as this replacement can be made up to power corrections (these are kinematic power corrections, see ref. [59]).For this latter point, it has been shown that higher-power corrections in fact lead to the replacement 2q + q − → Q 2 in the Wilson coefficients, as expected by Lorentz invariance.Next, we combine the hadronic tensor with the leptonic tensor to obtain the cross section for e + e − → 2 jets.In the above expression, the hadronic tensor is decomposed into several Lorentz structures.As an intermediate step, we record the following results for the contractions between these Lorentz structures and the leptonic tensor, Note that the power-suppressed contributions coming from the two Lorentz structures are identical, making it experimentally impossible to distinguish between the different types of power corrections.Combining the full expression for the hadronic tensor of eq.(3.33) with the leptonic tensor in eq.(2.12), we find that the cross section is given by Note that in the above expression the hard functions and jet functions are bare and unsubtracted quantities.

Subtraction and renormalization
In the previous section we arrived at a factorized expression for the cross section in terms of a set of jet functions and Wilson coefficients.We still need to treat the soft overlap between the collinear and anti-collinear region and renormalize all ingredients.Furthermore, the above result for the cross section contains some divergences, known in the literature as special rapidity divergences and endpoint divergences, that require treatment.In this section we discuss the overlap subtraction and the resulting cancellation of these divergences, as well as the renormalization of the jet functions.For other cases a scheme needs to be developed for treating endpoint divergences [60][61][62].First, in sec.4.1, we discuss the subtraction of the overlap region.In summary, the lack of a clear boundary between the collinear and anti-collinear regions gives rise to rapidity divergences, which cancel once the overlap region has been properly subtracted.In the method employed in this work, not all divergences cancel after the overlap region has been subtracted.These remaining special rapidity divergences already enter at lowest order in a s = α s /(4π), and are contained in the kinematic and higher-twist corrections.As a consequence, ill-defined convolutions between plus distributions and logarithms appear in the finite terms at higher orders in perturbation theory, which are commonly referred to as endpoint divergences.We discuss the subtraction and cancellation of these divergences in sec.4.2.The treatment of the overlap region is one of the main differences between the approach in SCET and the TMD operator expansion, which we compare in sec.4.3.Finally, the renormalization of the jet functions and the corresponding evolution equations are discussed in sec.4.4.

Overlap subtraction
As mentioned in sec.3.1, the two collinear regions overlap in the soft region, which is doubly counted.The presence of this overlap has been extensively discussed in the literature about TMD factorization [5][6][7], and it has a similar treatment in the background field method [8].
In the TMD operator expansion method, the overlap between the collinear and anti-collinear regions can be subtracted for bare operators by making the following replacements with similar replacements for the other combinations of operators that appear in eq.(3.17).Using the δ-regulator, the overlap factor is identical for both the leading twist and sub-leading twist operators, and is given by the soft function The overlap subtraction as presented in eq.(4.1) happens at the level of the product of two operators from the different collinear sectors, and the rapidity divergences cancel between the O n , O n and S.However, it is desirable to express the cross section in terms of ingredients that are each individually free of rapidity divergences.To define a rapidity-finite jet function, it is necessary to implement this subtraction by absorbing a square root of the soft function into each of the jet functions.To do this we use the fact that the soft function can be separated into two pieces, where and ζ is referred to as the rapidity scale.Using the above identity, the overlap subtraction procedure can be implemented by making the following replacements and similarly for the n-collinear direction.After the subtraction procedure, the leading-power jet function J bare,sub 11 is free of rapidity divergences but depends on the rapidity scale ζ.The dependence on ζ will be discussed in sec.4.4.

Special rapidity and endpoint divergences
The subtraction procedure of the TMD operator expansion method does not cancel all rapidity divergences.Specifically, J ′ 11 and J 21 , which respectively correspond to the the kinematical and higher-twist corrections, still contain divergences.These divergences can be subtracted and cancel in the cross section between the n and the n sectors.We now discuss the cancellation of these divergences separately for the kinematic and higher-twist corrections.

Higher-twist correction
For the twist-3 jet functions J 21 (ξ, b 2 ) divergences appear in the convolution with H 2 in eq.(3.36) as ξ → 0. In most of the literature this is referred to as the problem of endpoint divergences, while in the TMD operator expansion literature these divergences are considered part of the special rapidity divergences.In the case of TMD factorization, different methods for treating these divergences exist and we compare the two in sec.4.3.In this work, we subtract the singular behavior of the jet functions as ξ → 0 by means of a sub-leading-power soft function.This ensures that the integral of the hard function and the twist-3 jet function is well defined.This was first discussed in refs.[63][64][65], which focused on SIDIS.
The key insight is that in the limit ξ → 0 the gluon field strength becomes soft.In the soft limit, the gluon decouples from the collinear sector and the field strength can be factored out of the collinear matrix element.As a result, the limit ξ → 0 is captured by where . . .denotes terms which are regular as ξ → 0 and a similar expression holds for the opposite collinear sector.In this expression the jet functions on the right-hand-side are the familiar twist-2 jet functions and S 21,n is an NLP soft function.The NLP soft functions for the two collinear sectors are defined as This definition follows from replacing the collinear field strength operator in the O 21 operators in eq.(3.17) by a soft field strength, and following the same steps as in sec.3. Note that, strictly speaking, this derivation requires the introduction of a background field for the soft modes, although the soft function drops out in the final result for the NLP case we consider here.Using eq. ( 4.8), we can subtract the singular behavior from the jet functions in the ξ → 0 limit Note that the subtracted jet function is regular as ξ → 0. As a result, its convolution with the hard function is well defined and free of any endpoint divergences.
The divergent ln δ ± E and singular 1/ξ + parts of the jet function are now contained in the subtraction terms, which cancel in the hadronic tensor between the two collinear sectors as the jet functions only appear in the combination For this cancellation to work, one needs the two soft functions in eq. ( 4.9) to be identical, which is only the case if see ref. [24].Note that this constraint arises beyond leading order in a s , where the hard function contains terms involving ln ξ.Summarizing, the special rapidity and endpoint divergences in the higher-twist correction can be subtracted and cancel between the n and n sector on the condition that the rapidity regulators for the two sectors are related.

Kinematical correction
For the kinematical correction, only the following combination of jet functions appears after the subtraction procedure of sec.4.1 has been performed, Special rapidity divergences also appear in this expression because the derivative that defines J ′ 11 is taken before the soft function is divided out.
As for the higher-twist correction, the divergences in the above expression can be cancelled between the two collinear sectors.To see this we first define a subtracted jet function J ′ 11 , that is free of any rapidity divergences, because we take the derivative after subtracting.Using this definition, one can write where the special rapidity divergences are contained in the second term.Inserting this in eq.(4.13), we see that the rapidity divergent terms combine as follows, where the expression on the right-hand-side can be obtained using the exponentiation of the soft factor with D(b) the Collins-Soper kernel.We thus see that the special rapidity divergences for the kinematical corrections cancel.In fact if one takes δ + E = δ − E , as required by the cancellation of divergences for the higher-twist corrections, the subtraction terms cancel altogether.

Discussion of special rapidity and endpoint divergences
The subtraction of the overlap region and the cancellation of special rapidity and endpoint divergences is treated differently in SCET [24,[63][64][65] and the TMD operator expansion literature [8,25,50].Here we compare the two approaches.
In SCET, a soft mode is included explicitly, and its contribution at next-to-leading power leads to a family of soft functions [24].At the level of bare matrix elements however, the contributions of these soft functions cancel in the cross section and as a result no NLP soft factor needs to be introduced.Nevertheless, these soft functions can be used to subtract the overlap regions and define subtracted TMD distributions at NLP [63][64][65].On the other hand, in the TMD operator expansion method, no soft mode is included and the overlap between the regions is subtracted by the same soft function that appears at leading power.The main consequence is that in this method, at NLP, not all rapidity divergences are cancelled by the subtraction procedure.The divergences that remain after the subtraction procedure are referred to as special rapidity divergences.This includes the singular behavior in the higher-twist corrections that appear as ξ → 0, which in the literature are commonly referred to as endpoint divergences.
For both approaches, the remaining divergences cancel in the cross section between the n and n sectors.However, the subtraction and cancellation of the divergences is formulated differently between the two methods.In this work and in the SCET literature, it is formulated at all orders in perturbation theory by means of a soft matrix element.After the singular behavior of the jet functions has been subtracted, the divergent subtraction terms cancel in the combination Since the hard function contains logarithms of ξ at higher orders in perturbation theory, the only way for the divergences to cancel in the above expression is to have δ + E = δ − E .The above cancellation holds at the bare and renormalized level, however for the latter the distributions for the n and n sectors have to be evolved together from the same initial to the same final scale.
In refs.[9,25,50], on the other hand, the cancellation of special rapidity divergences is formulated order-by-order in perturbation theory.In ref. [25] a subtracted distribution is defined by where R is to be constructed order-by-order in α s such that the convolution between J 21 and H 2 is finite.Note that we have converted their expression to our notation.At leading order in perturbation theory the hard function contains no terms involving ln ξ, and so no constraint on the rapidity regulators arises.However, at higher orders in perturbation theory one would run into additional divergences coming from the integral of ln ξ in the hard function with the 1/(ξ − iδ ± E ) in the jet functions, leading to ln 2 δ ± E .It is the cancellation of these additional divergences that constrains δ + E = δ − E .We therefore conjecture that if the approach of ref. [25] is extended to higher orders, one would find such a constraint.

Renormalization and evolution
The subtracted jet functions defined above are at the bare level and require renormalization.Since renormalization and evolution takes place at the operator level, it is independent of the recoil-free jet algorithm that is employed.
For the twist-2 jet functions the renormalization is multiplicative, and we define the renormalized jet function as The renormalized jet functions depend on a renormalization scale µ and their evolution is given by Here, Γ cusp is the cusp anomalous dimension and γ V is the anomalous dimension of the quark vector form factor, which to one-loop are given by, The same equations apply to J ′ 11 as the renormalization and evolution kernels do not depend on b.The renormalization and evolution of the twist-3 operators is well known [25,49,[66][67][68] and involves a convolution in the momentum fraction ξ.A few simplifications allow us to recast the evolution presented in refs.[25,49] in a more elegant form.Since we consider large radii jets, our twist-3 jet functions depend on a single momentum fraction ξ.In our case ξ is restricted to the interval [0, 1], which can be shown by inserting a complete set of states in between each field that enters the matrix-element definition of the jet function.Because of this constraint, the mixing between different regions of momentum fraction space and mixing between T-even and Todd contributions, as discussed in ref. [25], is absent in our case.Following this, the renormalized twist-3 jet functions can be expressed as and the evolution reads The anomalous dimension for the twist-3 jet functions has been calculated to one-loop order [49] and is given by The subtracted and renormalized jet functions also depend on the rapidity scale ζ, in addition to their renormalization scale dependence.For all jet functions the rapidity scale dependence is described by the Collins-Soper kernel, but for J In principle, one could worry about the role of the subtraction term in the evolution.Specifically, one could argue that the introduction of subtraction terms alters the anomalous dimension.However, the contribution of the subtraction term cancels between the two jet functions in the cross section.Therefore, as long as one evolves the jet functions from the two collinear sectors together, the presence of the subtraction term will not have an effect on phenomenology.

Ingredients at order a s
In the previous sections we obtained a factorization formula for the hadronic tensor in terms of a set of jet functions.While the factorization relies on a recoil-free jet axis, the form of the factorization and the above definitions are independent of the precise choice.In order to make a theoretical prediction for the cross section, however, the jet functions must be calculated for a specific jet axis.We discuss the jet axes we consider in sec.5.1, and we will present the calculation of the corresponding jet functions in sec.5.2.

Jet algorithms
In this work we consider the jet axes defined through the E n -recombination schemes, with the special case of the winner-take-all (WTA) axis corresponding to the n → ∞ limit.This is a direct generalization of the p n T -scheme [69,70] to e + e − collisions.In principle, these axes also depend on the clustering, e.g.anti-k T vs. Cambridge/Aachen, but this distinction is not relevant at the order we are working, where the final state in the fixed-order jet function contains at most two partons 3 .The resummed jet function will of course contain the dominant effect of multiple emissions.
For all axes, we consider the large radius limit, such that all final-state particles will be clustered into the jet.For this is the reason, we write the final state of the jet function as |J alg ⟩.The jet axis enters our calculation as follows Here X is the collinear final-state produced by the fields, and the prefactor cancels that of the phase-space integral over the total jet momentum.Our transverse momentum measurement is encoded in the transverse position of the fields, by fixing our coordinate system such that the jet transverse momentum P alg J vanishes.This allowed us to use the same operators as in other transverse momentum measurements.
We now present the specific formulae to determine the transverse momentum of the jet for the E n scheme, which depends on the number of particles in the final state.In the case where the final state consists of a single parton with momentum k, we simply have P J = k.When the final state contains two partons with on-shell momenta k 1 and k 2 , respectively, the 3-momentum of the jet in the E n -scheme is given by where the complicated prefactor arises from the condition that the resulting 4-momentum is massless.However, up to NLP we can approximate ⃗ k 1 • ⃗ k 2 ≈ k 0 1 k 0 2 , since their angle is small due to being in the same collinear direction.This leads to (5.3) We will restrict ourselves to n > 1 for which the resulting jet axis is recoil-free and our factorization analysis is justified.For n = 1, P E n J is simply the total transverse momentum of particles in the jet, and thus not recoil free.The WTA recombination scheme corresponds to the n → ∞ limit, in which case eq. ( 5.3) simplifies to (5.4)

Calculation
We now calculate all jet functions to first order in a s .For convenience, we summarize their definitions here.The bare and unsubtracted jet functions are defined as with U 1 and U 2 defined in eqs.(3.8) and (3.9) and similar definitions hold for the jet functions in the n direction.The subtracted jet functions are given by with Finally, the renormalized jet functions are given by (5.8)The diagrams that give a vanishing result are not shown.

Bare jet functions
The (leading-order) matrix elements that are required in this calculation are given by with the relevant diagrams shown in fig. 2. The index a on the polarization vector indicates the color of the gluon (though technically is not part of the polarization).We have omitted the virtual corrections, as the virtual correction to J 11 vanishes in dimensional regularization, and the virtual correction to J 21 is exactly zero.Inserting the expressions in eq.(5.9) in eq. ( 5.5), we find to order a s These expressions only depend on the details of the jet algorithm through P alg J .From here, one performs the integral over k using δ (d−2) (P alg J ), and the lightcone components of p and k.For the E n scheme, leading to (5.14)The remaining integral over p in J 11 can then be performed by using and the integral that appears for J 21 can be obtained by differentiating with respect to b.Finally, we perform the integral over ξ for J 11,n and expand the result in ϵ ) with (5.18)

Subtraction and renormalization
To obtain the subtracted and renormalized jet functions at first order in a s , we need the LP and NLP soft functions and the renormalization factor Z J11 to order a s .We do not need the one-loop expression for Z J21 at this order.The LP soft function with δ regulator has been calculated to one-loop order in ref. [6], and the result reads The renormalization factor for J 11 is given by The integrals over k + and k − can be performed using the delta functions and the remaining integral over k can be performed using eq.(5.15).The final result for the NLP soft function reads

Final results
We now present the final results for all the ingredients of the cross section to first order in a s .The subtracted and renormalized jet functions read (5.23) Note that since J 21 and J ′ 11 only start at order a s , the renormalization does not effect them at this order in perturbation theory.The WTA axis corresponds to the n → ∞ limit.In this case, The twist-3 jet function J 21 is shown in fig. 4 for n = 2, 3 and ∞ (which corresponds to WTA).By dividing J 21 by a s , we removed dependence on the renormalization scale.Note that these jet functions have a singularity at ξ = 1 2 arising from the E n scheme.However, since the hard function is smooth at this point, there is a cancellation between ξ = 1 2 + δξ and ξ = 1 2 − δξ to yield a finite result (principle value).For completeness, we also present the one-loop results for the hard functions, (5.25)

Conclusions
Due to the very high precision of the Large Hadron Collider and the advent of the Electron-Ion Colliders, there has been an upsurge in interest surrounding TMDs, as this presents a unique opportunity to determine these distributions with unprecedented precision.This necessitates a deeper understanding of the factorization theorem for transverse momentum dependent cross sections.
In this study, we examine select aspects of the recently developed operator basis that contribute to the differential cross sections at next-to-leading power.We focus on jets, yielding an infrared and collinear safe observable, that can be calculated in perturbation theory.In order to further simplify the description, we consider only large radii such that all energetic radiation is inside jets.By using a recoil-free jet recombination schemes, we avoid non-global logarithms.It is natural to examine how higher twist operators that appear at NLP enter in the description of jets and generate new measurable effects.We illustrate all these aspects on a prototypical process, e + e − → 2 jets, which is theoretically simple and experimentally relevant.The final factorized cross section in eq.(3.36) was obtained by starting from the power expansion of the hadronic tensor, decomposing in spin structures and applying discrete symmetries.It involves only one independent twist-2 and twist-3 jet operator matrix element for each collinear direction.To our knowledge, this is the first time that the contribution of a twist-3 operator to a jet has been determined.We have found that NLP effects exist in e + e − → 2 jets only when the two jets are different, e.g. when two different recoil-free jet algorithms are used or a jet is replaced by an identified hadron (TMD fragmentation function).
It is well known that factorization theorems at sub-leading power contain additional divergences that require careful treatment.In sec.4.2 we have shown that these divergences can be subtracted by means of a soft function and cancel in the cross section between the two collinear sectors.Moreover, we demonstrated that the cancellation of these divergences requires that the jet functions of the two collinear sectors are evaluated at the same rapidity and renormalization scales, and they need to be evolved together when considering resummation.The evolution equations necessary to perform resummation were presented in sec.4.4.
Having arrived at a finite and factorized expression for the cross section, we computed the twist-3 jet functions using WTA and E n recombination schemes at leading order in perturbation theory.We discussed how these jet algorithms can be implemented in perturbative calculations.All results that are needed to obtain a resummed expression for the cross section are presented to order a s in sec.5.2.
There are several aspects of this study that we consider relevant for future applications.First of all, the measurement we consider allows for a substantial simplification of the NLP formalism, making it more palatable.The process that we describe can be measured at the Belle experiment, but we leave a detailed phenomenological analysis to future work.It would be intriguing to consider how the nonperturbative effects that we study impact QCD research at the Future Circular Collider (FCC).More importantly, measurements of hadron structure at the EIC will benefit from our NLP jet analysis, which we consider a priority.Our framework can be naturally extended to also take into account polarized targets at the EIC.

Figure 1 .
Figure 1.Kinematics of the di-jet production in e + e − .

Figure 2 .
Figure 2. Diagrams for the matrix elements that enter the calculation of the twist-2 jet function at NLO (left and middle) and the twist-3 jet function at LO (right).The double line represents the Wilson line, and the ⊗ denote insertions of the quark field (at the corners) and gluon field strength (on the double line).The diagrams that give a vanishing result are not shown.

Figure 3 .
Figure 3.The only diagram contributing to the NLP soft function at leading order.The ⊗ denotes the field strength operator and the double lines represent the Wilson lines.