Transverse momentum dependent operator expansion at next-to-leading power

We develop a method of transverse momentum dependent (TMD) operator expansion that yields the TMD factorization theorem on the operator level. The TMD operators are systematically ordered with respect to TMD-twist, which allows a certain separation of kinematic and genuine power corrections. The process dependence enters via the boundary conditions for the background fields. As a proof of principle, we derive the effective operator for hadronic tensor in TMD factorization up to the next-to-leading power (∼ qT/Q) at the next-to-leading order for any process with two detected states.


Introduction
The transverse momentum dependent (TMD) factorization approach has a long history that started from the "DDT formula" [1], and Collins-Soper-Sterman resummation [2,3]. Nowadays, the TMD factorization approach is a well-developed framework that consistently describes various processes in terms of TMD distribution functions (for a review of the present status and latest development, see [4][5][6]). Even though the recent achievements, the TMD factorization theorem still lacks the systematicness of collinear factorization theorems. The main reason is that in its heart, the collinear factorization is based on the method of the operator product expansion [7], which allows the first-principles treatment, independently of any hadronic model. The absence of a similar systematic approach for JHEP01(2022)110 the TMD factorization raises a series of problems, especially when extending the formalism beyond the leading power approximation. In this work, we propose a formal derivation of the TMD operator expansion, overcoming the present status quo. Structurally, it is similar to the light-cone operator product expansion [8,9] for some specific (and broadly known) cross-sections.
The TMD factorization describes inelastic processes with two-detected states (initial or final). The main examples are the Drell-Yan (DY) process, semi-inclusive deep-inelastic scattering (SIDIS), and semi-inclusive annihilation (SIA). The information about the nonperturbative structure is stored in the TMD distribution functions (it could be TMD parton distributions, TMD fragmentation function, TMD jet-function, etc.). The order parameter of the TMD factorization is q T /Q, where Q is the momentum of the hard probe and q T is its transverse component. The latest global analyses [10][11][12] demonstrate that the TMD factorization is valid for q T /Q 0.25. Beyond that limit, the TMD factorization curve quickly deviates from the experimental measurements. To extend the region of TMD factorization, one needs to incorporate power corrections.
In the past years, computations of power correction for TMD factorized cross-sections have been made by several groups, see [13][14][15][16][17][18][19]. These computations required an update of existing methods. In this work, we go further and present a different approach, which is a natural extension of techniques used for studies of power corrections in the collinear factorization, and which has not been discussed in the framework of TMD factorization yet (as far as we know). Namely, we use the background field method and derive the factorization directly in the space of field functionals (operators). For that reason, we call the technique -TMD operator expansion. Starting from the definition of the QCD Lagrangian, we re-derive some known results such as the leading power (LP) TMD factorization, and we proceed extending it to the next-to-leading power (NLP). Since the computation is made at the level of operators, the results can be applied to other cases (partly extending the present discussion), such as processes with jets or factorization with generalized TMD distributions (GTMDs). Many solutions and more minor results derived in this paper are already known or discussed in the literature, but the operational framework presented here is totally new.
There are several sources of power corrections to the factorization theorem. We classify them as follows: • power suppressed contributions to the cross-section due to convolution between hadronic and leptonic tensors or phase-space factors. Usually, these corrections can be accounted exactly, and in any case, should not be mixed with QCD corrections. We solely concentrate on the hadronic tensor, and we do not consider these corrections in the present work; • corrections in the values of kinematic parameters that flow with the change of the ratio q T /Q. For example, the values of effective momentum fractions in DY process are x 1,2 = q ± /P ± 1,2 = x Bj 1,2 1 + q 2 T /Q 2 . Often, one expands these variables into series, generating a tail of power suppressed contributions. Such expansions are unnecessary and even harmful since the values of kinematic variables are unique for all powers,

JHEP01(2022)110
and their exact accounting supports the frame-independence of the final expression. We make the computation in position space, which allows us to identify values of kinematic parameters unambiguously from the definition of the hadronic tensor; • kinematic power corrections. These are contributions of a higher power with the operator content of lower power terms. The Wandzura-Wilczek relations [20] provide a famous example, which can also be generalized to the TMD case [21]. These corrections inherit the structure of lower power terms. They are essential for the restoration of global properties violated by the factorization theorems, such as electromagnetic (EM) gauge invariance, translation and frame independence, etc. We demonstrate that NLP kinematic corrections restore the EM gauge invariance of hadronic tensor up to q 2 T /Q 2 , as expected; • genuine power corrections incorporate new operators, and hence new TMD distributions.
In our work, we study only QCD specific power corrections, i.e., kinematic and genuine. The separation of kinematic and genuine power corrections is an important and non-trivial task [22]. The most efficient way to solve it consists of ordering the operators with respect to their evolution properties. For the collinear distributions, it implies the separation of operators with respect to the geometrical twist. For TMD distributions, we introduce the notion of TMD-twist, such that TMD distributions with different TMD-twist have separate evolution, and thus their matrix elements are independent observables. The essential feature of the TMD factorization is the appearance of infinite light-like gauge links [23,24]. They are responsible for all distinctive features of TMD distributions, such as rapidity divergences, nonperturbative evolution, T-odd distributions, and process dependence. The reverse of the medal is that one cannot derive TMD factorization from a local operator expansion. This fact leads to the absence of such an important ordering criterium as the twist of the TMD operator. So far, that has not been a problem, since the leading term of the TMD factorization can be derived with the analysis of Feynman diagrams [25], or using the soft-collinear effective theory (SCET) approach [26], and it does not need the systematization of TMD operators. However, dealing with power corrections requires the definition of some ordering for the operator basis. In this work, we introduce the twist of the TMD operator (or TMD-twist) given by a pair of numbers, which are geometrical twists of collinear substructures of the TMD operator. This definition of TMD-twist allows for a certain separation of kinematic and genuine contributions of power corrections.
The main aim of this work is the development of a theoretical basis for TMD factorization beyond the leading power. To keep the discussion as general as possible, we perform all computations at the level of operators without specifying a process or referring to cross-sections. The operator representation allows us to keep expressions relatively compact and avoid long algebraic structures that appear at the TMD distributions level. The expressions for NLP cross-sections will be made in subsequent publication.

JHEP01(2022)110
The background field method is particularly suitable for the computation of power corrections. For applications of this method to the collinear factorization see e.g. [22,[27][28][29]. This method computes the operator expansion directly from the function integral, avoiding any matching procedure typical for many approaches. It allows keeping track of the internal structure and the operators' relations at each evaluation step. The computation is naturally performed in position space (although it can easily be turned to momentum space), simplifying the NLP computations. In the case of TMD factorization, the background field Lagrangian [30,31] must be updated to the case of two independent background fields (we call it a composite background field), which is done in appendix A.
The paper is split into two logical parts. The first part includes sections 2, 3, 4, and it provides a review of the factorization method in general terms. In particular, section 2 introduces the main definitions, such as the definition of hadronic tensors, the composite background field, counting rules, and the notion of effective field operators. Section 3 is devoted to the problem of boundary conditions and gauge fixation for background field. We demonstrate that the choice of boundary conditions and an adequate gauge fixing are related to the analytical properties of generating functions, which in turn are tied to the underlying process. Accounting for these properties leads to the famous process dependence at the level of TMD operators. Section 4 is devoted to the general discussion about the structure of power and perturbative expansion. In this section, we introduce the concept of TMD-twist and TMD operator expansion. The second part consists of sections 5, 6, 7, 8 and 9, and it is devoted to the computation of TMD factorization at NLP and NLO. For pedagogical reasons, the computation is presented with many details. The tree order and definition of all operators is given in section 5. Section 6 presents the NLO computation of the hard coefficient function. The problem of soft overlap and subtraction of overlap region is discussed in section 7. Section 8 considers the properties of TMD operators and it is split into subsections 8.1, 8.2, 8.3 that treat respectively rapidity divergences, ultraviolet (UV) divergences and the renormalization of TMD operators. Finally, in section 9, we demonstrate the cancellation of divergences among elements of TMD factorization, fix the scheme dependence and derive the evolution equations for NLP TMD operators.
The text is supplemented by appendices, which contain additional technical details. Appendix A presents the expression for the QCD Lagrangian in the composite background field. Appendix B demonstrates the technique of computation in the composite background field. Appendix C contains the expressions for evolution kernels in momentum space.
One of the difficulties in writing about power corrections comes from the notation. On the one hand, the topic is complicated and requires accurate and exhaustive formulation. On the other hand, all-inclusive writing conceals the general structure of the expression, which is often simple. For that reason, we accept the following convention: we drop the parts of notation that are not important in the present context, such as, color and spinor indices, arguments, etc. However, we keep all essential elements, and a cautious reader should be able to restore all missed components if needed.

JHEP01(2022)110 2 General structure of TMD factorization
In this section we provide the notation and the basic definitions that are used in this work.
General setup. We study the hadronic tensors for Drell-Yan, h 1 + h 2 → γ * + X, semiinclusive deep inelastic scattering (SIDIS), h 1 +γ * → h 2 +X and semi-inclusive annihilation (SIA), γ * → h 1 + h 2 + X, where J µ (y) is the electro-magnetic (EM) current (for brevity we omit the electric charge) J µ (y) =qγ µ q(y), (2.4) and q(y) is a quark field. At sufficiently high energies, the current should be replaced by the electro-weak (EW) current. The difference between EM and EW currents is only in the tensor structure and does not impact the factorization procedure, so that we stick to the EM case.
The factorization for all three processes is almost identical, since (in our approach) it is derived at operator level. The main differences appear in the boundary conditions for the fields and the sign of the Fourier exponent (section 3). For concreteness we center our discussion on the factorization of the DY reaction, commenting necessary modifications for SIDIS and SIA cases.
The kinematics of the process is defined by the photon momentum q µ , and the hadrons momenta p µ 1 and p µ 2 . To avoid complications related to the target mass corrections we assume that hadrons are massless, Hadrons momenta define two light-cone directions, which we traditionally denote as n µ andn µ with (nn) = 1, We also introduce the usual notation for components of the light-cone decomposition of a vector, and v µ T is a component orthogonal to n µ andn µ . The invariant mass of the virtual photon is In the case, of SIDIS Q 2 = −q 2 . The TMD factorization is derived (as it is demonstrated later) in the limit where Ψ p is the hadron's wave function (formed at the distant past), and S QCD is the QCD's action. The superscript (±) indicates that the element is composed only from causal/anti-causal fields. The functional integration measure incorporates all necessary normalization factors.
In the case of SIDIS or SIA the only change in eq. (2.11) is in the hadron wavefunctions, that should be replaced by the wave-functions of produced hadrons according to the process.

JHEP01(2022)110
The next step is to identify the relevant field modes and the ones that are to be integrated, or kept. We suppose that a fast-moving hadron is composed only of fields with momenta along the hadron's one. So, for a hadron with the momentum along n, the fields which constitute it obey the counting (same for (+) and (−) components) where λ is a generic small scale, λ ∼ Λ/Q. Similarly, a hadron with the momentum along n direction is composed out of fields with the counting In the literature these fields are identified as n andn-collinear or collinear and anti-collinear fields (f.i. [36,37]), or target and projectile fields (f.i. [35]). The background fields, in its own kinematic sector, are ordinary QCD fields and thus satisfy the QCD equation of motion (EOMs). Let us emphasize the " "-sign which is used in these formulas. It indicates that the fields incorporate all possible momenta with the corresponding boundary. This is a principal difference of the background approach from SCET, where the field modes are defined in "boxes" of momentum space, alike {∂ + , ∂ − , ∂ T } qn ∼ Q{1, λ 2 , λ} qn, i.e. their momentum is localized around a given scale (that labels the fields). Another difference is that modes scaling as {∂ + , ∂ − , ∂ T } q ∼ Q{λ, λ, λ} q (called "soft" in SCET nomenclature) are not present in our construction. It is convenient to distinguish "good" and "bad" spinor components of the quark field qn /n (x) = ξn /n (x) + ηn /n (x). (2.14) They are defined as (same for (+) and (−) components) The (massless) quark EOMs / Dq(x) = 0, in the terms of these components are

JHEP01(2022)110
Using tree-level computations (see section 5) one finds the power counting for individual components, The counting rules for components of the gluon field also follows from eq. (2.17), (2.18), Effective operators. In essence, the background-field method consists in splitting QCD fields into dynamical and background components with the subsequent (functional) integration of the former. In many aspect, the procedure resembles the Wilsonian renormalization group formulation and it is the simplest way to obtain the operator product expansion (OPE) with any given power counting. It is also a very explicit method to study the power corrections to any OPE (see e.g. power corrections to DIS [9,27], DVCS [22], quasi-distributions [38], small-b matching for TMD distributions [29]). The distinctive feature of the TMD factorization from cited cases is that the background field has two independent components, collinear and anti-collinear. Thus the causal/anti-causal fields in the functional integral in eq. (2.11) split as where ψ and B are the dynamical components which cover the remaining part of the Hilbert space, see figure 1. Since the fields represent independent Fourier components, 1 the integration measure can be split as The separation of the field modes is done respecting gauge-invariance. The Lagrangian is invariant only under the gauge-transformation of all fields by the same transformation parameter (irrespective of any power counting). Nonetheless, the gauge-transformation for the background field and the dynamical field can be decoupled, if one uses the so-called background gauge [30,31]. In this case the (covariant-gauge-like) gauge fixation condition for the dynamical field takes the form where D µ is the covariant derivative (2.19). This gauge fixation modifies vertices of background-to-gluon (ghost) interaction. The advantage of such a choice is that the background fields transform independently from the rest and their gauge can be fixed in any convenient way. This is one of the most profitable features of the background field method, which essentially simplifies the analysis. Formally, the collinear and anti-collinear background fields transform with the same gauge parameter. In the region where the momenta of the fields do not overlap (i.e. for ∂ µ > λ 2 Q) we treat them as totally independent fields, and thus their gauge transformations are also independent. The problems arise where sectors overlap (see barred part of figure 1). To avoid this region, we momentarily restrict our discussion to the non-small-x domain with p parton λQ. We will return to the discussion of the mode overlap in section 7.
After implementing these definitions in the functional integral in eq. (2.11), the hadronic tensor reads where in the square brackets we indicate the field content of each term (for brevity we omit superscripts (±) on these arguments and indicate similar arguments for currents by dots). The label "unsub." states that this expression has an unsubtracted overlapped region that is discussed in section 7. The gauge fixing terms are included in the S QCD exponents. The cross-modes-interaction term is (2.27)

JHEP01(2022)110
Its explicit form in the background field-gauge is derived in appendix A, eq. (A.9). The derivation takes into account the equations of motion (EOMs) for collinear and anticollinear fields. The cross-modes-interaction term can be split into four terms: S nh (Sn h ) (A.10) describing the interaction of (anti-)collinear fields with hard fields; S nn describing the direct interaction of collinear and anti-collinear fields (A.11); and S nnh (A.11) describing the interaction of all fields simultaneously. Each action S nh and Sn h is equal to the usual QCD action with background field [31], while S nnh and S nn are specific for the composite background case. Let us mention that S nn = O(λ 3 ) due to the power-counting in eq. (2.21), (2.22). Before integrating over the hard modes we specify the content of the hadronic wave functions. The main assumption of the parton model is that the hadron is composed of fields collinear with respect to its momentum, (2.28) Once the wave functions are independent of hard mode, we integrate over it and deduce the expression where J µν eff [qn,q n ,q s , . . .] depends on all background (causal and anti-causal, collinear and anti-collinear) modes and is defined as The effective operator satisfies which are consequences of symmetry and transversality of the hadronic tensor.
In the end of the section let us sketch the further steps of the factorization procedure, which are discussed in detail in section 4. The effective operator is an infinite series of individually gauge-invariant terms. This sum can be ordered with respect to λ, where N represents power order (J µν N,k ∼ λ N +4 ), and k enumerates the operators with the same power counting. Generally, each J µν N,k is a convolution of background fields, and it can be written in the schematic form where C µν  N,(a,b) is a coefficient function, O a are some operators (the possible multi-index and multi-position structure is encoded in the single label a) that also depend on y, and ⊗ represents a convolution in variables and indices.
Let us note, that ordering of operators J µν N,k by power counting implies the expansion of components along "slow"-direction, f.i. qn(x + nx 0 ) = qn(x) + O(λ 2 ) (for non-extreme x 0 ). The resulting derivatives contribute to the operators at higher powers. Expansions that involve y µ need special care. The counting rules for the components of y µ are (yq) ∼ 1, that gives Therefore, the transverse derivatives which involve y have the compensating factor y T ∼ λ −1 . The combination y µ T ∂ µ ∼ 1 is not suppressed, in contrast to other transverse derivatives. Due to it, all dynamics in the transverse plane in the final expression is tied to y.
Changing the (functional) integration and summation order in eq. (2.29) we get and the fields in definition of Φ's are just collinear or anti-collinear. In these expressions the functions Φ unsub (that are unsubtracted TMD distributions) are nonperturbative in the sense that they contain unknown information on the hadronic structure and low-energy QCD interactions. The sketched derivation remains the same for the cases of SIDIS and SIA. The only difference among effective operators in different processes is due to different boundary conditions prescribed to collinear fields for initial and final state hadrons (see the next section). In the rest of the paper, the effective operator is the same for all cases.

Process dependence and gauge fixation
The effective operator is gauge invariant term-by-term in the series in eq. (2.33), as a consequence of the gauge invariant definition in eq. (2.30). Fixing the gauges for collinear and anti-collinear fields in a convenient way, we can simplify calculations in the intermediate

JHEP01(2022)110
steps restoring gauge invariant expressions at the end of the computations. We use the light-cone gauges for background fields, This choice removes O(1) gluon components in eq. (2.22). Therefore, the power counting for operators increases with the number of fields in the operator, which essentially simplifies the computation. In any other gauge, there would be an infinite set of operators of the same order in power counting but with different numbers of A + n or A − n , which eventually sums into Wilson lines but makes the computation cumbersome.
The light-cone gauge conditions in eq. (3.1) alone do not remove all gauge freedom, leaving possible z ∓ -independent transformations (for A ± = 0 gauge). The detailed discussion of this can be found in ref. [23]. To fix the redundant gauge freedom, one imposes specific boundary conditions on the components of the gluon field. The final expression is independent of the choice of boundary conditions. However, an inappropriate choice can lead to unnecessary complications during the computation, and to avoid them we specify convenient boundary conditions for each process.
To deduce the proper boundary conditions, we study a generic expression (for regular graphs) contributing to an effective operator, where fn (n) is a combination of (anti-)collinear fields and their derivatives that are (eventually) integrated along the light-cone positions with some weights. The only important property is that fn (n) depends only on z −(+) . The dependence of fn (n) on z +(−) would violate the counting (there is also dependence on y T , but it is external and does not modify analytical properties). The power α is non-integer in the dimensional regularization. Examples of such expressions can be found in eq. (6.3), (6.4), (6.5), (6.9). The +i0 prescription is the appropriate one for the causal-sector of interactions. For the anti-causal sector it is replaced by −i0, but simultaneously all analytical properties are reverted, leaving the final result unchanged. The denominator of eq. (3.2) possesses a branch cut (for definiteness we place the branch cut of x α along the line (−∞, 0)) along the line This is illustrated in figure 2. The field combinations f in eq. (3.2) give raise to parton distributions and thus we must assume some analytical properties for them to guarantee the existence of Fourier integrals. For the case of incoming partons (i.e. PDFs) f is analytical in the lower half-plane, whereas for the case of outgoing partons (i.e. FFs) f is analytical in the upper half plane. Therefore, depending on the process we deal with different combination of analytical properties, summarized in the following table for DY for SIDIS for SIA fn(z − ) is analytical in lower lower upper half-plane. f n (z + ) is analytical in lower upper upper half-plane. Given these properties we can deform the integration contour placing it on the sides of the branch cut as it is shown in figure 2. Finally, the integral splits into several contributions. The contribution at infinity is troublesome, since it produces ill-defined operators with gluons concentrated at infinity. A properly selected boundary condition nullifies it. We illustrate this discussion for the case of DY reaction. In this case, we can close the integration in the lower half-plane, and deform the contour as shown in figure 2. We get where I k are integrals of fn(z − )/(z − ) α along elements of the contour. I 1 and I 2 are integrals along sides for the branch cut, I 0 is a semi-circle at z − = −i0 and I ∞ is a semi-circle at z − = −∞ − i0. The dimensional regularization regularizes a possible ultraviolet pole at z → 0 (demanding α < 1), and thus I 0 vanishes. However, it also implies that I ∞ contributes to the integral, unless the fields f vanishes at z − → −∞. We can use the freedom to fix the redundant gauge condition such that I ∞ vanishes. 2 The same analysis can be done for z + variable. So, we conclude that for the DY reaction the convenient set of boundary conditions is such that fields fn (n) vanish at z −(+) → −∞.
Repeating the same analysis for the cases of SIDIS and SIA, we arrive to the following set of appropriate boundary conditions for DY: lim for SIA: lim This set guarantees the absence of an I ∞ contribution in the loop integrals. For future convenience, we introduce the variables L andL which can take the values ±∞, and summarize boundary conditions as The value of (L,L) is the only difference between processes at the operator level. Notice that for DY and SIA cases the integrals in eq. (3.2) can be written in the factorized form whereas for SIDIS case both integral representations I DY and I SIA are valid.
Once the boundary conditions are specified, the gauge-fixing condition in eq. (3.1) can be inverted. It gives where F µν is the gluon field-strength tensor. After the computation is complete, and the result is written in terms of F µν 's, one can restore the explicit form of the gauge-invariant operator by multiplying fields with gauge links. The rules are For the anti-collinear fields the expressions are analogous. Finally, we mention that in SCET one often uses the operator A µ ⊥ that (for a collinear field) is defined as [36,41,42] (3.14) In the explicit form this operator reads Comparing with eq. (3.11), (3.12) we find that A µ ⊥ = A µ in light-cone gauge, which gives a map between our expressions and the ones written in the terms of A µ ⊥ .

JHEP01(2022)110 4 General structure of the TMD operator expansion
The computation of the effective operator in the composite background follows the common pattern of computation in a single background field, see, e.g., refs. [22,29,38]. The logical steps following in the computation are: 1. Expansion of eq. (2.30) into monomials of collinear and anti-collinear fields.
3. Rewriting of fields in terms of "good" and "bad" components.
5. Reduction of operators to a given basis, using algebra and EOMs.

Fiertz transformation into TMD operators.
During the evaluation, one should keep in mind that in the very end, the operators are inserted into matrix elements, and some of them vanish (e.g., due to non-zero fermion number or due to non-singlet color representation). Such operators can be eliminated without full consideration. At each step, we find a series of terms with increasing power, and the ordering of power counting is preserved. Thus only terms with the desired power counting are finally kept.
General structure of operators. The starting point for the effective operator expansion is the product of two EM currents separated by a distance y. The fact that y T ∼ (λQ) −1 spoils the usual intuition about the computation of the effective operator because it is not allowed to expand over y T , since (y µ T ∂ µ ) ∼ 1, and the effective operator splits into two parts separated by y T . These parts are independent in the sense that they have separate anomalous dimensions and can be separately expanded in powers (this, however, does not mean that the coefficient function is a product of coefficient functions).
After step 2, a contribution to the effective operator has a general form is a light-cone operator composed from causal collinear fields, that are positioned at the light-ray with transverse coordinate y T with coordinates {z − }, and similar for other V 's. C µν is the differential operator that possibly contains a (loop-)integral to be evaluated at step 4. The symbol ⊗ indicates the integral convolution in positions z ± and contraction in Lorentz and color indices between coefficient function and operators. The illustration for spatial configuration of fields in eq. (4.1) is given in figure 3.
Each operator Vn({z − }, r T ) in eq. (4.1) consists of some number of fields located on the light-ray pointing from r T to Ln (and analogously for V n ({z + }, r T )). The positions {z − } JHEP01(2022)110 are distributed along this light-ray. Each field within V n is accompanied by a semi-infinite Wilson line in eq. (3.12). These Wilson lines recombine with each other and partially cancel. However, in general, a Wilson line does not vanish beyond the position max{z − } and continues till Ln. In this sense, the operator V is semi-compact. It also means that the operator V is not entirely gauge-invariant because, under the gauge-transformation, it receives a gauge-rotation factor at (Ln + ∞ T ). These properties are in contrast to usual DIS-like OPE, where resulting operators are compact (i.e., localized in the finite volume) and entirely gauge-invariant. On the other hand, the semi-compact operator basis is the only difference between DIS-like and TMD operator expansions. Therefore, we can apply the powerful machinery of OPE, correcting it only for semi-compact operators.
Step 3 decomposes operators with respect to power counting of the components of the fields in eq. (2.21), (2.22). The resulting operators have the same form as in eq. (4.1), with the difference that operators V have a "pure" power counting and can not be decomposed further.
Twist-decomposition. The aim of step 5 is to decompose each operator V with respect to a convenient basis. The most convenient and physically motivated decomposition is the decomposition with respect to geometrical twist, which is the "dimension-minus-spin" of the operator. One has where U is an operator with a definite twist (for simplicity, we set it equal N , although there can be several terms with the same twist in the decomposition), c N is an integral-differential operator, and ⊗ is the integral and matrix convolution. An operator U N with a definite geometrical twist belongs to an irreducible Lorentz group and thus does not mix with operators of other geometrical twists. This important property is preserved by perturbation JHEP01(2022)110 theory, and thus operators with definite geometrical twists have an independent evolution, and their matrix elements are separate physical observables. The twist-decomposition in eq. (4.2) is an algebraic procedure (consisting of symmetrization and anti-symmetrization of indices) originally defined for local operators. It can be also generalized for non-local operators by means of generating functionals of local operators (see e.g. [43]), or by implication of differential operators (see e.g. [27]), or by considering the conformal transformation properties (see e.g. [44]). All these methods cannot be applied directly to semi-compact operators V and should be revised. In particular, the method of reconstruction of non-local form of operators via generating functions has been applied to semi-compact operators in ref. [45]. The main idea used in ref. [45] is to set the parameter L finite. It allows to write down a formal local expansion for semi-compact operators, make the twist-decomposition and resum the result into a generating functional. Lastly, the limit L → ∞ is taken. This method correctly reproduces known lower-power properties of semi-compact operators, and it can be applied to operators of any power. In ref. [45], the twist-decomposition for several cases of semi-compact operators has been made, including P-odd operators, which are specific for the semi-compact case.
The lowest twist for semi-compact operators is twist-1. The twist-1 operator is just a single "good" component of quark ξ or gluon field F µ+ (with µ being transverse index) with attached semi-infinite Wilson line. Although such numbering looks unusual, it also follows from the formal counting of geometrical twist by adding up conformal spins [46].

Recombination of divergences. At this stage a contribution to the effective operator has the form
where C is a combination of C from eq. (4.1) and c's from eq. (4.2). The coefficient function C has divergences in (in the dimensional regularization) that are IR. These divergences match the UV divergences of operators U . Let Z N be the renormalization factor for the operator U N , where the convolution is in the position of operators, and µ is the renormalization scale. Then inserting unit factors 1 = Z −1 N ⊗ Z N into eq. (4.3) we obtain an expression of the form where operators U are renormalized at the scale µ, and To demonstrate that eq. (4.5) is finite is a non-trivial task and the object of the factorization theorem.
Let us note that the rapidity divergences do not appear in the effective operator and have no traces in the coefficient function. The origin and the factorization of rapidity divergences are discussed in the section 7.

JHEP01(2022)110
TMD operators. The Fierz transformation at step 7 recouples the color indices and groups operators U into color-neutral TMD operators where we explicitly indicate the color indices A and B (that belong to representation R AB of SU(N c )), and dim(R AB ) is the dimension of their representation. After this operation a contribution to the effective operator has a general form where C contains also derivatives that act on O's. This is the final form, see also eq. (2.34), of the effective operator. For convenience, we define the TMD-twist of the operator O N M equal to (N+M) ("Nplus-M"). Such notation matches the usual jargon. In particular, the leading power TMD distributions that are often referred to as twist-2 TMD distributions (without specification of the meaning of twist for TMD operator) are twist-(1+1) TMD distributions within our formalism. The sub-leading power TMD distributions are referred to as twist-3 and have TMD-twist-(1+2) or TMD-twist-(2+1). In principle, the operators with twist-(1+2) or twist-(2+1) are different and define two separate TMD distributions with different evolution equations, although the C-conjugation relates them. The real profit from this notation comes from the higher powers. So, the twist-4 TMD distributions (in usual terminology) can be twist-(3+1), twist-(2+2) and twist-(1+3). The properties of twist-(2+2) TMD distributions are very different from twist-(3+1) distributions, and they do not relate to each other by any means.
In the limit of small transverse separation y T → 0, TMD operators turn to collinear operators. Consistently, this limit is computed by the light-cone OPE, see e.g. [29,45]. The leading term is equal to the product of U 's at y T = 0, The smallest possible geometrical twist for the operator on r.h.s. is (N + M ), since it is a product of spin N and M tensors. Therefore, at high-q T , where ∼ y T contributions are small, and TMD factorization turns into the resummation approach, TMD distributions of twist-(N+M) match collinear distributions of twist (N + M ) or higher. In this way, computing contributions of only TMD-twist-(1+1) operators (at all powers of TMD operator expansion), one should be able to reconstruct all leading twist terms of collinear factorization, including the fixed order computations, such as in ref. [47].

Effective operator at NLP/LO
Using the scheme depicted in the previous section, we compute the effective operator to leading and next-to-leading power (LP and NLP, respectively) and up NLO in perturbation theory. In this section, we derive the tree order of NLP effective operator and introduce necessary definitions. Although this result is (partially) known, see e.g. [19,24,48], our derivation is novel in many aspects because it is made at the operator level and with no explicit reference to a specific process. For this reason, we give a detailed explanation for each step of the computation. The NLO computation is given in the next section.
Tree order for LP. We start with the decomposition of the EM current Here, the first line provides the leading tree-order contribution. The second line contains fields ψ that are to be contracted with ψ from S QCD/int . These terms contribute to NLP and NLO, and they are considered in the next section. The terms in the third line have two fields from the same collinear sector, and thus they produce disconnected contributions to matrix elements, unless extra fields are taken from S int . The first possible non-vanishing contributions of the terms in the third line happen only at N 4 LP.
To get the LP term, we consider the first line of eq. (5.1), and perform the multipole expansion. We get qnγ µ q n (y) +q n γ µ qn(y) The terms in O(λ 4 ) contains derivatives, such as y +qn ← − ∂ − γ µ q n . Decomposing the field q into "good" and "bad" components as in eq. (2.14) we obtain qnγ µ q n +q n γ µ qn =ξnγ µ T ξ n +ξ n γ µ T ξn (5.3) where we omit the arguments understanding implicitly that each collinear field depends on (y − n + y T ) and each anti-collinear field depends on (y +n + y T ). The first, second, and third lines in eq.
3) has the form of eq. (4.1), with operators having pure counting, in the sense that no further expansion is needed. This term represents the LP term of EM current. The second and the third lines of eq. (5.3) have indefinite twist. They are discussed in the following subsection. Thus, the LP term is a composition of twist-1 and twist-1 operators

JHEP01(2022)110
Let us note that the LP current in eq. (5.4) violates the EM charge conservation. Indeed, However, the expression in r.h.s. of eq. (5.5) is O(λ 3 ), and thus formally the charge is conserved at LP. Combining together the LP currents, eq. (5.4), we get the LP expression for the effective operator The quark fields ξ are operators of twist-1, and thus J µν LP (y) is of (1 + 1) × (1 + 1)-twist in our nomenclature. Combining the fields into TMD operators, eq. (4.6), we get where (i, j, k, l) are spinor indices, and the TMD twist-(1+1) operators have the argument ({y − , 0}, y T ). They are defined as The operator O 11,n and O 11,n are obtained by replacing n ←n everywhere, including the subscripts of the fields. This expression is well known and the basis for the LP TMD factorization, see e.g. [24,48,49]. The matrix elements of operators O 11,n and O 11,n give rise to the quark and anti-quark TMD distributions, respectively.
Tree order for NLP. The NLP part of the effective operator can appear only via the combination of a LP current, eq. (5.4), with NLP part of another EM current. The NLP contribution of order ∼ λ 3 to J µ can be composed in two ways. The first one is combining a "good" and a "bad" component of the quark fields from the second line of eq. (5.3). The second possibility is to have three "good" components of fields, e.g.ξA T ξ. To get such a term, one needs to pull down an interaction term from e iS int and couple it to the first line of eq. (5.4). The diagrams representing the NLP contributions are shown in figure 4. The diagrams A and B correspond to the second line of (eq. (5.3)). The remaining diagrams represent the interaction contribution. The diagrams C, D, E and F are specific for the computation in the composite background field, and would be absent for the case of ordinary background field. The reason is that an ordinary background field does not couple to the dynamical fields via 1PI vertex. In other words, there are no vertices with a single dynamical field and background field(s). Such diagrams (in the sum) compose EOM for the background field, and therefore, vanish (e.g. a diagram C with An replaced by A n is not present). This fact is already taken into account in the construction of the effective background action [30], and thus the corresponding vertices are not present in S int . In the case of the composite background, a vertex with a single dynamical field and different background fields is not forbidden, and actually it is provided by Sn nh in the effective action, eq. (A.11). These vertices do not contribute to the EOM. 3 As a result the diagrammatic expansion is not (explicitly) symmetric with respect to n µ ↔n µ . This symmetry is restored once the operators are rewritten in a unique basis. As an example, we present in detail the evaluation of diagrams A and C, which form a symmetric pair. Diagram C reads (we set the global position of the current to 0 for brevity) where we used the propagators in dimensional regularization, eq.
, and the integral over (d − 2) transverse components can be evaluated, eq. (B.11). The result is where we take into account that only the transverse components of the gluon field contribute to LP term. Note that in this expression one cannot cancel z + between the denominator and the numerator, because it would spoil the analytical properties of the integrand. To evaluate eq. (5.11), we recall the analytical properties of An summarized in eq. (3.4), then we close the contour of the z − -integration in the lower (DY and SIDIS cases) or upper (SIA case) half-plane, shrink it to the pole at z − = i0/(2z + ), and evaluate the residue. We obtain the expression In this expression the anti-collinear fieldsξn / An ,T form an operator of twist-2.

JHEP01(2022)110
The diagram A is given by a simple expression Here, the fieldηn has indefinite geometrical twist, which can be checked e.g. by conformal transformation [28]. Practically, it implies that the field η mixes with a quark-gluon pair during the evolution. To rewrite η in the terms of definite twist operators, we apply EOMs, eq. (2.17). In the present case, we need the EOM forηn, eq. (2.17). In the light-cone gauge it can be rewritten Here, the first term on r.h.s. is the total derivative of twist-1 operator, and the second term is the twist-2 operator. The twist-counting can be confirmed by expanding the operators in a series of local operators (setting L finite), and computing the twist of each term in the series. Inserting EOM into eq. (5.13) we get The second term is complementary to diagram C, eq. (5.12). Together they from a transverse expression. To make it explicit we rewrite diagram A in the form of eq. (5.11) and sum the diagrams together. We obtain Evaluating the rest of the diagrams in the same manner we obtain the expression for the EM current at NLP. We split the result into the following parts The first term contains the derivatives of twist-1 operators The terms J µ 21 and J µ 21,8 contain operators of twist-2 and twist-1 where A is the color index in the adjoint representation and t A is the generator of SU(N c ). n ↔n (also in the subscripts of fields). The contributions to J µ 21,8 and J µ 12,8 arise from the diagrams of E and F (plus n ↔n diagrams). The collinear and anti-collinear parts of J µ 21,8 are in the adjoint representation of the color group. These parts of the EM current do not contribute at NLP, because to form a color neutral TMD operator, eq. (4.6), an another operator in the adjoint representation is required. The first non-zero contribution from J µ 21,8 and J µ 12,8 into effective operator takes place at N 2 LP. Let us define the inverse derivative operator, as With this operator the expressions in eq. (5.18), (5.19) have a simpler form where the positions of fields on r.h.s. and l.h.s. are formally the same. Note that the inverse derivatives in the last line acts on the quark and gluon fields together. The operators J µ 1 1 and J µ 11 combine with the LP current eq. (5.4) and form The terms J µ 1 1 and J µ 11 restore the electric charge-conservation at ∼ λ 3 and ∼ λ 4 orders. Indeed, Clearly, the operators with a particular twist combinations form series where each next term restores the charge conservation to a higher power. These series are known as series of kinematic power corrections. So, the operators J µ 1 1 and J µ 11 are the kinematic power corrections to the LP current. Due to the that fact the kinematic power corrections restore the global properties of current (and hadronic tensor), all operators in the series must have the same coefficient function. We demonstrate it explicitly at NLO in the next section.

JHEP01(2022)110
The effective operator at NLP is obtained by composing LP and NLP terms of EM currents. We obtain where δ ij is the Kronecker delta for spinor indices. All derivatives are with respect to the coordinate y, and act only on the subsequent operator. To derive this expression we have used the identitȳ and similar for anti-collinear fields. To derive this identity we assume that the total derivatives of TMD operators can be eliminated, since they do not contribute to the forward matrix elements. Here, TMD operators of twist-(1+1) have the arguments

Gauge invariant expressions for TMD operators.
The definitions in eq. (5.29)-(5.32) are given in light-cone gauge. They can be written in gauge invariant form using the relations in eq. (3.11), (3.14). There are two usual ways to write the twist-(2+1) operators, the one typical for the SCET literature e.g. [18,42], and the one typical for direct QCD computations e.g. [19,24]. Both have advantages and disadvantages. In particular, the SCET-like notation is convenient for the computation of hard coefficient function, whereas the traditional-like is convenient for computation of the evolution properties. Nicely, they are related by a simple transformation. In the present work we use both kinds of notations, designating them by different fonts.

JHEP01(2022)110
The elementary building blocks are the semi-compact operators of twist-1 and twist-2. They are where index µ is transverse and we omit the transverse links. Both operators have an open spinor and color indices. The operator U µ 2,n has open spinor and vector indices, and is more general than the operator appearing at NLP. The latter reads where both sides are spinors. The semi-compact operator built from anti-collinear fields are obtained from Un with replacement n ↔n. In addition, we define the operators Analogously, the reduced twist-2 operator U 2 is The TMD operators are products of semi-compact operators. So, the twist-(1+1) operators in eq. (5.8), (5.9) are simply The superscripts (±) indicate that semi-compact operators are made out of corresponding causal or anti-causal fields. Notice the enumeration of positions in operators O 21,n and O 12,n . It is adjusted such that the gluon field has uniformly coordinate z 2 . All these operators are matrices in spinor space, and singlets in color space. The operators O are related to O by the charge-conjugation. The matrix elements of these pairs define quark and anti-quark TMD distributions, such as in ref. [48]. The definitions of eq. (5.33)-(5.37) are given in the traditional QCD basis. In the SCET literature instead, one would use
The relations between the operators O and O is 48) where N and M are 1 or 2. The relation (5.48) can be inverted Similar expressions hold for O-type operators.
TMD operators in momentum space. Ordinary, the TMD distributions are defined in the mixed momentum-coordinate representation. Namely, they are Fourier-transformed just in light-cone coordinates. Such a transformation provides a similarity with the parton densities (which are defined in terms of fractions of momentum) and it preserves a simple structure of the TMD evolution (which is diagonal in transverse-position space). We also follow this practice and define where x 1,2 is a shorthand notation for (x 1 , x 2 ), which can be interpreted as the fraction parton's momentum, and p + is the hadron's momentum. The Taking into account the translation invariance, we define where M + N = 3, and

JHEP01(2022)110
At this point, we cannot not specify the domain of integration over x. However, once the matrix elements are taken, the values of x are restricted x ∈ [−1, 1] for parton distributions, and x ∈ (−∞, −1] ∪ [1, ∞) for fragmentation functions. The integral measure in eq. (5.57) takes into account the translation invariance of the TMD operator defined in eq. (5.52). 4 In principle, our computation (since it is done in position space) does not imply any simplification that comes from the translation invariance. Nonetheless, we use it to make notation somewhat lighter. The relation between the operators O and O is simplified in momentum space. According to eq. (5.48), we have where all partons momenta are incoming. Performing the transformation for eq. (5.4), (5.17) we get where x 2,1 = (x 2 , x 1 ) and the repeating argument y T is suppressed for brevity. Note that the arguments x 1 and x 2 are adjusted such that x 1 is the momentum fraction of quark or anti-quark, and x 2 is the momentum fraction of the gluon.
Effective operator in momentum space. Similar expressions can be written for the effective currents eq. (5.7), (5.27). In this case, it is convenient to take into account the

JHEP01(2022)110
Fourier integral that is present in the hadronic tensor, eq. (2.1). Defining In the twist-(1+2) and twist-(2+1) operators the argument is x = (x 1 , x 2 , x 3 ). Note that the operators J µν 1211 and J µν 2111 are J µν 1112 and J µν 1121 with n ↔n. In the case of twist-(1+2) TMD operators, the values of x 1,2,3 are not sign definite. The momentum conservation delta-function fixes the value (and the sign) of only one of them. The other two variables are integrated and can be positive or negative.
To derive these expressions we took into account that the global position of the currents is irrelevant for DY, SIDIS and SIA processes. Also note, that the definition of the

JHEP01(2022)110
Fourier transform in eq. (5.63) is taken similar to the DY case in eq. (2.1). According to eq. (2.2), (2.3) the sign of q is the opposite for SIDIS and SIA.
The expressions for the currents (5.65)-(5.69) are simple to generalize to the case of electro-weak interaction. In this case, one needs to replace γ µ T → γ µ T (g V ± g A γ 5 ), and δ → g V ± g A γ 5 (the sign of g A depends on the term). The general structure and other factors is the same.

NLO perturbative correction to NLP operator
The computation of perturbative corrections to effective operators follows the same path as in section 4. Here we find another difference between TMD factorization and collinear factorization. Namely, the loop integrals can produce extra suppressing factors. That happens due to the unhomogeneous counting rule for the vector y in eq. (2.35). In particular, the diagrams with the interaction of causal and anti-causal fields are suppressed by a factor λ 2 . So, the twist-(1+1)×(1+1) contribution of the exchange configuration is N 2 LP (see refs. [15,18] for a discussion about these contributions). Therefore, at NLP we can consider only the interactions within a single causal sector, i.e. for each EM current independently.
At NLO and NLP we have only two-point and three-point relevant diagrams, that are shown in figure 5 and figure 6, respectively. The two-point diagrams couple operators of twist-1 and operators of twist-1 or 2, while the three point diagrams couple only operators of twist-1 and twist-2.
The computation of the diagrams with the background field is slightly different from the ordinary computation of amplitudes because a part of the loop-variables are also arguments of the background fields. Therefore, the integration over these arguments cannot be computed and we obtain the integral convolution in eq. (4.1). To compute the integrals over the rest of the loop variables we used a technique, which has been developed in ref. [29]. This technique is based on a series of shifts for integration variables and leads to the usual Lorentz-invariant loop-integrals. The technique works equally well for any type of power-suppressed diagrams. For pedagogical purposes in appendix B we present a detailed computation of the three-point diagram 5 in figure 6.
Alternatively, one can present the background field as a Fourier image and perform the computation in momentum space, using standard methods. However, the loop-integral in momentum space becomes complicated once a large number of background fields participate. For example, three-point diagrams in momentum space can contain polylogarithms already at one-loop level, whereas one-loop expressions in position space are polynomial for any number of external fields.
Another complication of a momentum space computation is the necessity to specify the direction of the parton's momentum, which is crucial for the computation of Feynman variable integrals. This is not a problem for two-point diagrams, where one has a unique choice (however, different for DY, SIDIS and SIA cases), but for the three-point diagrams one has already three possible combinations of momenta directions. This is because the individual momenta of the quark-gluon pair can have different sign, and the momentum conservation fixes only the sign of their sum, see second and third lines of eq. (5.64).

JHEP01(2022)110
The computation presented here has been done in position space. It is the first computation of the Sudakov form factor at NLP and the first one made in position space. As a cross-check, we have performed also the NLO computation in momentum (for DY kinematics), and checked that the results coincide with each other. Also we have checked that the LP coefficient coincides with the known one, and NLP coincides with the recent computation 5 in refs. [50,51].
Two point diagrams. The two-point diagrams, shown in figure 5, contribute to LP and NLP effective operator. Note, that in figure 5 we have already split the quark field into ξ and η components for convenience. The diagram 1 is the only diagram that has a LP contribution. Explicitly, the diagram 1 reads where C F = (N 2 c − 1)/2N c is the eigenvalue of the quadratic Casimir operator for the fundamental representation of SU(N c ).
The background fields in eq. (6.1) are expanded along the light-cone, as and similar for other fields. The expansion eq. (6.2) can be safely done under the sign of the loop integration. Each next term increases the counting of the operator. At NLP only the operators with one transverse derivative contributes. Substituting the decomposition (6.2) into eq. (6.1) we obtain a 2d-dimensional integral. Two variables x − and z + are arguments of background fields and thus the integral with respect to them cannot be computed. The integration over the rest (2d − 2) loop-variables is done by the method explained in appendix B. We get
The two-points diagrams contributing to the NLO of the effective operator. The diagram 1 contributes to LP, and the diagrams 2 & 3 to NLP. The diagrams with n ↔n should be added. The blobs indicate the type of background field. The blue lines are the dynamical fields. and dots denote the higher-power terms. Similarly we compute the diagrams 2 and 3, Eq. (6.3), (6.4), (6.5) are structured as in eq. (3.2), and thus they can be rewritten in factorized form as in eq. (3.9) or eq. (3.10). For example, the LP part of the diagram 1 in the case of DY process reads diag LP;DY The prefactor of eq. (6.6) is finite for → 0. However, both integrals are UV divergent at z ± → 0, and regularized by > 0. As it is shown below, these poles exactly reproduce the Sudakov double pole. At z ± → ∞ the loop integrals are regularized by the natural field decay. Other diagrams also can be written in this process-dependent form. However, it is more convenient to keep expression as in eq. (3.2), which makes many properties transparent. The diagrams 2 and 3 are to be rewritten using EOMs eq. (2.17), (2.18). To do so, one should rewrite EOMs in the integral form eq. (5.14), substitute it into the diagrams written as in eq. (6.6), exchange the order of integration, and integrate the remaining integrals. This computation can be simplified using the inverse derivative operators eq. (5.21), and the relations As a result, the sum of the two-point diagrams is

JHEP01(2022)110
The expression for mirror diagrams is equal to eq. (6.8) with n ↔n. Note, that the last two lines can be rewritten with inverse derivatives, reproducing the operator J µ 21 eq. (5.23). The operator in the second line of eq. (6.8) is the J µ 11 current of eq. (5.24). As expected, the coefficient functions for all terms of J µ 11 are the same, such that the current conservation eq. (5.25) is preserved. Since the contribution to NLP part is the result of the combination of several diagrams, it gives a strong check of our computation.
Three-point diagrams. The three-point diagrams are shown in figure 6. The diagrams 7-10 are specific for the composite background field and would be absent in the usual background field computation. Note, that the three-gluon vertex that appears in diagrams 6 and 9 is not equal to three-gluon vertex in QCD but has a modification that comes from the background-gauge gauge-fixing condition.
The computation of these diagrams is straightforward and described in details in appendix B. Here we present the final expression for the sum of three-point diagrams. For convenience we add theξAξ-part of the two-point diagrams, such that the result is the full JHEP01(2022)110 expression for the coefficient function of J µ 21 operator. It reads The terms in eq. (6.9) are grouped such that each line forms a transverse combination.
NLO expressions in momentum space. The passage to momentum space is straightforward. The expression for EM current in eq. (5.59) takes the form where C's are the coefficient functions. Their expression up to NLO are , (6.12) In these expressions the momenta q ± are set in accordance to eq. (5.64). I.e. for the J µ Eqs. (6.12), (6.13) are the bare form of the coefficient functions, that contains the IR poles. These poles are removed by the operation in eq. (4.5). In the present case, this operation can be made for EM currents individually before recombining them into the effective operator. The renormalization constants Z 1 and Z 2 are derived in section 8, and JHEP01(2022)110 they exactly remove the pole part ofC's (see eqs. (9.24), (9.26) and discussion there). So, we obtain in the MS-scheme 6 + O(a s ), (6.14) Here and in eqs. (6.12), (6.13), the notation is adopted such that x 2 is the momentum fraction of the gluon field. The expression for C 1 coincides with earlier computations, see e.g. refs. [25,26,37,52,53]. Nowadays, the coefficient C 1 is known up N 3 LO order [54]. The coefficient functions for the anti-causal sector are obtained from the causal ones with complex conjugation. The expression for J µν in eq. (5.64) became decorated by the coefficient functions, where asterisks denote the complex conjugation, which, in fact, applies only to L Q .

Mode overlap and the soft factor
So far, we have assumed that the collinear and anti-collinear fields are entirely independent. It allows us to impose individual gauge-fixing conditions and separate fields into independent gauge-invariant TMD operators. However, it should be kept in mind that there is JHEP01(2022)110 a part of functional integration phase space where collinear and anti-collinear fields are a single background field. In figure 1, this region is covered by diagonal shading. This is the so-called soft region (or glauber region in the SCET nomenclature). Fields in this region (marked by s) satisfy the counting The soft region is double-counted in the functional integral with the measure of eq. (2.24). Let us stress that the double-counting of the soft region does not effect the TMD factorization procedure described in previous sections. Instead, each TMD operator has an uncompensated rapidity divergence. In fact, the rapidity divergences should cancel between collinear and anti-collinear TMD operators, but they cannot due to the double counting of the soft region. There are several solutions of this problem. Let us list some of them: • The definition of collinear and anti-collinear fields can be modified in the soft region, such that there is no overlap. For example, by introducing a cut in rapidity for each field as it is shown by the red-dashed line in figure 1. Then each TMD operator depend on the cut parameter, such that this dependence is compensated in the product. See e.g. discussion in ref. [49]. The same idea is used in the rapidity factorization approach [13,55].
• The product nature of the functional integral allows to remove double-counting simply dividing by the (functional) integral over soft modes. It is possible if the hadron states do not contain soft fields, which is valid in a non-small-x regime. The resulting factor is known as the soft factor [25,26] or zero-bin subtraction [56]. It is the most popular procedure nowadays. However, there is no general approach to determine the soft factor operators at higher powers. Most plausible, such a simple multiplicative structure does not hold for higher power operators.
• One can ignore problems of overlapping modes entirely, and reconstruct necessary parts (such as rapidity renormalization constants) by demanding that the effective operator is well-defined. This logic is used in the collinear anomaly approach, see [37,57].
In the present NLP computation, we use the second way, because it leads to the correct result at NLP without additional computation. However, we expect that this approach does not hold at higher powers. Let us stress, that all approaches should result into the same final expression, up to some finite terms. The fixation of finite terms is equivalent to the fixation of the scheme, and the definition of the parton distribution. The closest example is the difference between MS and DIS schemes, see e.g. [58]. The only difference of the TMD case from the collinear is, that the coefficient of rapidity divergence is nonperturbative, and thus the finite terms added/subtracted from the physical distribution are also nonperturbative. Consequently,

JHEP01(2022)110
the scheme must be defined by a certain nonperturbative statement. At LP one fixes the scheme defining that cross-section for DY and SIDIS do not have extra nonperturbative functions except TMD distributions. The same definition can be applied to NLP, see section 9.
Determining the soft factors. To determine soft factors for our operators, we use the following procedure. We split soft parts of collinear and anti-collinear fields and similar for other components. Then we isolate the soft fields into a single factor dropping power suppressed contributions. For example, for the first term of the LP effective operator eq. (5.7) where [a, b] s is the Wilson line with the soft gluon field, and we omit transverse links for brevity. The operator for the LP soft factor is where the trace is taken with respect to color indices. The trace and the factor 1/N c appears due to the fact that only gauge-invariant operators have non-zero matrix elements. To derive eq. (7.5) we have used the counting rules for soft fields eq. (7.1), (7.2). The operator S represents the soft part of LP effective operator. Therefore, the soft part to the functional integral eq. (2.11) at LP is given by the vacuum (assuming that the hadrons do not carry soft partons) matrix element of eq. (7.5) The same structure follows from the region-separation method [59]. In SCET literature, this procedure is known as a zero-bin subtraction [56]. We remark that the problem of overlapping modes does not impact TMD factorization, and thus the replacement in eq. (7.7) is valid to all orders in perturbation theory. A similar computation can be done for operators contributing to J µν NLP . We have two principal cases: the operators with derivatives (the first and the second line in eq. (5.27)),

JHEP01(2022)110
and operators with extra field A µ T (other lines in eq. (5.27)). In both cases, we obtain that the soft overlap contribution is equal to the LP soft factor eq. (7.6), since a derivative of a soft Wilson line, or an extra factor A µ s necessarily increase the power counting. Therefore, the subtraction of the soft region for NLP operators has the same form as for LP operator eq. (7.7). Namely, (7.8) and similarly for other terms of effective operator eq. (5.27). The LP soft factor is independent on y ± and thus does modify convolutions in J µν eff at NLP. Having the same soft factors for LP and NLP operators leads to the following consequences: • LP and NLP operators must have the same rapidity divergence and, as the result, the same rapidity anomalous dimension.
• LP and NLP operators must have the same collinear divergent part of the UV renormalization.
Indeed, these divergences arise in the interaction of soft modes, and (in the present approach) they are canceled by the soft factor. In section 8, we independently derive both statements, and we explicitly verify them at NLO. A simple procedure described here allows to determine soft factors for LP and NLP terms. For higher power correction a more systematic procedure should be developed.

δ-regularization.
The soft factor has a complicated combination of divergences. Namely, it has UV divergence, rapidity divergences, and mass divergences (see refs. [60,61] for detailed analysis). The mass divergences cancel in the sum of all diagrams [61], whereas rapidity and UV divergences remain.
An essential feature of rapidity divergences is that they are not regularized by dimensional regularization [62]. Therefore, an additional regularization must be implemented. There are many regularizations of rapidity divergences used in the literature, such astilting of Wilson lines [25], analytic regularization [57,63], exponential regulator [64] and δ-regularization [26,65]. Each of these regularizations has been used in plenty of computations and has advantages and disadvantages. The final result after the recombination of divergences is independent on the rapidity regularizator. In this work we use the δregularization, for the only simple reason that we are experienced in it.
The rapidity divergences arise due to the interaction with the far end of the halfinfinite light-like Wilson line [61]. In the δ-regularization these interactions are regularized by insertion of dumping factor into Wilson line, with δ − > 0 ands = sign(L). Thus, there are two regulator parameters δ + and δ − , which regularize divergences associated with different light-like directions.
Soft factor at NLO. The calculation of the LP TMD soft factor has been performed in many papers, see e.g. NLO calculations [25,26,66,67]. The expressions used here with δ-regularization are taken from ref. [68]. The bare soft factor at NLO reads Here, b is a transverse vector, so −b 2 > 0. This expression contains a product of rapidity divergences associated with different Wilson lines in the form of ln δ + and ln δ − . Note, that some 1/ poles of eq. (7.11) are not the UV divergence, but a part of rapidity divergence. The UV part of the soft factor should be computed separately. In MS-scheme (see footnote 6), it is [65] It contains ln(δ + δ − ) which, in this case, is a remnant of the mass divergence. In ref. [61], it is proven that the rapidity divergences of the LP soft factor can be written as a product of factors. The proof is made using the method of conformal transformation and it is valid to all orders in perturbation theory. Using it we present the LP soft factor in the form where ν ± are some scales, and ν 2 = 2ν + ν − , and S 0 is free of UV and rapidity divergences. The NLO expression for R can be deduced from eq. (7.11). In MS-scheme it reads For future convenience, we introduce UV renormalization factor Z R for R, and write Note, that the identity is valid to all orders in perturbation theory.

JHEP01(2022)110 8 Divergences of TMD operators
The TMD operators are composed of two semi-compact light-cone operators, eq. (4.6), separated by a transverse distance, The hadronic matrix element of the TMD operators in eq. (8.1) defines the TMD distributions. It can be a TMD PDF(s) or TMD FF(s) The TMD distributions F N M and D N M can be decomposed over independent Lorentz components and reveal a plethora of TMD distributions. It is widely known that there are eight quark TMD PDFs defined by the LP operator O 11 in eq. (5.53), see e.g. [48,69]. All these distributions (TMDPDF and TMDFF) obey the same evolution equations, because they are matrix elements of the same TMD operator. In this way, all these distributions are alike from the perspective of TMD operator expansion, despite the fact that they have very different partonic interpretation and are measured with different experimental set-ups.
In this work, we concentrate on the global properties of TMD factorization, and thus we do not systematize NLP TMD distribution (see [24,48]). This systematization as well as, the derivation of cross-section is the object of a subsequent publication. Instead, we study the global properties of LP and NLP TMD operators, and write down the evolution equations for LP and NLP TMD distributions.
The operators U 1 and U 2 in eq. (5.33), (5.35) together with their C-conjugated versions U 1 and U 2 , eq. (5.36), (5.38), set up all TMD operators at NLP eq. (5.39)-(5.44). As we show in the following subsection the singularity and evolution properties of a TMD operator follow from the properties of each U that composes it. Therefore, we concentrate on the studies of U 1 and U 2 rather then on studies of O. In addition, we consider a more general operator U µ 2 eq. (5.34), since it does not complicate the computation but allows us a simpler comparison with known expressions.

Rapidity divergences
Rapidity divergences are specific of TMD operators, eq. (8.1). They arise in the interaction of fields with the distant segments of light-like Wilson lines (accurate definition and properties rapidity divergences can be found in ref. [61]).
Renormalization of rapidity divergences. Rapidity divergences are multiplicatively renormalizable. In ref. [61] this statement is proven for the multi-parton scattering soft JHEP01(2022)110 factors. The proof can be easily generalized to TMD operators O N M,n that appear at LP and NLP. 7 Using the multiplicativity of the rapidity divergence we write the LP operators in eq. (5.8), (5.9) as where we omit the argument {y − , 0} of the TMD operators. The variable ν + is a scale of rapidity divergences renormalization. The operators on the r.h.s. of eq. (8.4) are free from rapidity divergences. We distinguish such operators by explicit indication of rapidity renormalization scale ν + . The rapidity-divergence renormalization factor R is independent on y − , and thus the Fourier transformed TMD operator eq. (5.53), (5.54) is renormalized in the same way. The same renormalization factor absorbs rapidity divergences of NLP operators eq. (5.41)-(5.44), and similar for the Fourier-transformed operators eq. (5.53)-(5.56).
7 One needs to apply the conformal transformation Cn (defined in eq. (5.1) of [61]) to the TMD operator ONM,n. The resulting operator CnONM,n is spatially-compact. It has the form of a Wilson line with two light-like segments that are joined at the origin of the light-cone. The partonic fields are positioned along the Wilson line. The Wilson line has a light-like cusp, whose UV divergence corresponds to the rapidity divergence of the TMD operator before transformation. Using the fact that an UV divergence is multiplicatively renormalizable, and that the conformal invariance of QCD is restored in the Wilson-Fisher critical point [70], one derives that rapidity divergence is multiplicatively renormalizable as well. The derivation is made by iterations order-by-order in perturbative expansion, starting from the LO which respects conformal invariance. Practically, this derivation repeats the one given in section 5.2 of ref. [61] for multi-parton soft factors. The light-like cusp anomalous dimension associated with the cusp of CnONM,n (called the soft anomalous dimension, [71,72]) corresponds to the rapidity anomalous dimension of the TMD operator. This correspondence is a simple equality at LO, but receives modifications beyond the LO due to the breaking of the conformal invariance in QCD. The modification terms can be derived using the same iterative method, which has been done in ref. [73] up to N 2 LO. The result coincides with the three-loop brute force computation [74], which non-trivially validate of the method of conformal transformation used to prove the renormalizability of rapidity divergences.
One of the important consequences of the derivation is that the rapidity divergence of LP operators and NLP operators are the same. It follows from the fact, that cusp divergences of CnO11,n and CnO21,n coincide. It independently confirms the same observation pointed out in section 7 after eq. (7.8). Rapidity divergence at NLO. The factor R that renormalizes the rapidity divergences of TMD operators is the same as for the soft factor eq. (7.13), (7.14) and here we check it at one-loop. For this purpose, we consider a TMD operator composed from U 1,n or U µ 2,n eq. (5.33), (5.34) and any other semi-compact operator U N,n . We denote such operators as O N 1,n = U N,n U 1,n and O N 2,n = U N,n U µ 2,n . These operators are more general than operators appearing at NLP, however, as we demonstrate below, their possibly complicated nature does not impact the rapidity divergences structure.

JHEP01(2022)110
TMD operators O N 1,n and O N 2,n are color-neutral. Thus, the operator U N,n has a Wilson line in the anti-fundamental representation pointing to Ln. We split a far segment of this Wilson line from the rest of the operator, where z 0 is such that max{|z − |} < |z 0 |. The rapidity divergences arise only in the interaction of the Wilson line [z 0 n, ±∞n] with U 1,2 . One-loop diagrams that produce the rapidity divergence are shown in figure 7.
To extract the rapidity divergence, we follow the method developed in ref. [29]. Let us describe it using the diagram r1, as an example. We write the diagram r1 using the background field method with a single background field in the A + = 0 gauge. It reads where P + = γ − γ + /2 is the projector of the "good" component of the quark field, see eq. (2.15). Here we have also used that the TMD operators in eq. (8.1) can be written as a single T-ordered operator (or in terms of the functional integral, the superscripts (±) can be omitted), since all fields in it are separated by light-like or space-like distances. Note, that the spinor indices of U N and the rest of the diagrams are not contracted. Joining the propagators by the Feynman parameter α and making a shift y → y + α(nσ + b), we obtain

JHEP01(2022)110
where we used that P +/ yγ + q = 2y + ξ and b + = 0. The integral over y can be computed in the sense of the generating functional. For that we expand ξ in a Taylor series at y = 0, and integrate this series term-by-term using eq. (B.6). The result reads The n = 0 term is rapidity divergent at α → 0. Indeed, in this limit the field ξ is independent on σ and the integral over σ diverges at σ → L. To reveal the divergence we make a change of variable τ = ασ, and obtain (8.11) where dots indicate the terms finite at α → 0, and thus rapidity-divergence-free. In this expression the divergence is transparent. As eq. (8.11) is independent of z 0 , the operator on U N,n can be promoted to U N,n by limiting z 0 → L.
To regularize rapidity divergences we use the δ-regularization eq. (7.9). It gives the factor e −Lδσ in the integral eq. (8.8), and modifies eq. (8.11) as (8.12) where s = sign(L). Now, the integrals over α and τ can be computed. The result is wherep + ξ = −i∂ + is the momentum of field ξ. The direction of the Wilson line does not impact the rapidity divergent part, but gives rise to a finite term ∼ isπ. The computation of similar diagrams for the C-conjugated operator, gives the same rapidity divergent part.
The most important observation is that the rapidity divergent term in diagram r1 is independent of the second part of the operator. In this sense, we can associate the rapidity divergence in eq. (8.13) with the operator U 1,n . In order to get the complete rapidity divergence of the TMD operator O N 1,n = U N,n U 1,n , we should consider also diagrams where the fields of U N interact with the Wilson line of U 1 .
Next, we study the rapidity divergence associated with U µ 2,n ({z 1 , 0}, 0 T ), eq. (5.34). For that we consider the TMD operator O N 2,n . There are two rapidity divergent diagrams shown in figure 7, r2 and r3. The computation is similar to eq. (8.8)-(8.13). The rapidity divergent part of these diagrams is

JHEP01(2022)110
wherep A is the momentum operator acting on the gluon field. Again the r.h.s. of eq. (8.14), (8.15) reproduces the original operator O N 2,n . Summing together diagrams r2 and r3 we obtain (8.16) wherep is some generic momentum. We observe that the coefficient of the rapidity divergence for LP eq. (8.13) and NLP eq. (8.16) operators coincides, as it is predicted in previous sections. Using eq. (8.13) or eq. (8.16) (together with their corresponding charge-conjugated parts) we deduce the NLO expression for the (unrenormalized) factor R defined in eq. (8.4)-(8.5). In MS-scheme, we have This expression coincides with eq. (7.16), which validates the factorization at one-loop. An important feature of the rapidity divergences is that they have a nonperturbative contribution. The computation presented in this section misses this part. It is clear that at large transverse distances b the gluon propagator is modified by the confinement effects. Therefore, the perturbative result in eq. (8.16) is valid only at small (but finite) values of b 2 . At large values of b 2 the factor R gets power corrections. In dimensional regularization, the presence of these power corrections is indicated by renormalon divergences [75]. The renormalization theorem for rapidity divergences guarantees that the factors R for TMD operators of twist-(1+1) and twist-(2+1) coincides to all powers of small-b 2 expansion. It allows us to assume that rapidity divergences for LP and NLP operators also coincides nonperturbatively.
Rapidity divergences of quasi-partonic operators. The consideration on rapidity divergences and the soft factor presented above allows us a generalisation going far beyond NLP. That is, the rapidity divergences and hence the rapidity anomalous dimensions for quasi-partonic TMD operators coincides with the LP.
We define quasi-partonic TMD operators in analogy to collinear quasi-partonic operators. Namely, a quasi-partonic operator U N is the operator whose geometric twist equals the number of fields 8 (excluding Wilson lines), see e.g. [28,46,76]. Consequently, a quasipartonic TMD operator is composed from two quasi-partonic operators U . It also implies that a quasi-partonic operator consists of "good" components of fields only. E.g. operators U 1,n and U µ 2,n are quasi-partonic. The first non-quasi-partonic semi-compact operator has twist-3.
The equality of rapidity divergent parts for all quasi-partonic TMD operators can be derived in several ways. First of all, one can observe that the soft part of such an operator coincides with a staple of Wilson lines. Indeed, since the operator already contains a JHEP01(2022)110 maximum number of fields (and all of them are "good" fields), the replacement of any field by its soft part as in eq. (7.3) increases the power counting. Therefore, the part of the soft factor, responsible for the cancellation of the rapidity divergence of a quasi-partonic operator, coincides with the LP soft factor.
Another way is to use the method of conformal transformations (see footnote 7). In this case, we utilize that "good" components preserve the projection properties after the conformal transformation Cn [44]. Therefore, the soft anomalous dimension of the transformed operator coincides with the one of LP operator. It leads to the equality of factors R for all quasi-partonic operators. In other words, where R is the same as in eq. (8.4) and the operator on r.h.s. is rapidity-divergence-free. Eq. (8.18) is simple to confirm at one-loop. There are two types of diagrams contributing to the rapidity divergence of a quasi-partonic operator, the interaction with a gluon or a quark field. These diagrams are computed in eq. (8.14), (8.15), and they have the same expressions apart of color factors. In the quasi-partonic case, these color factors should be replaced by [1]t a [2]t a (for interaction with quark) or [1]if abc t b [2]t c (for the interaction with gluon), where [1] and [2] are color matrices of the fields before and after the interacting field. Summing together all such structures, one obtains the original color structure multiplied by t c t c . I.e. rap.div.[U quasi-p.
where C K is the eigenvalue of the quadratic Casimir operator for color representation of U quasi-p. N .

UV divergences
The UV divergence of the TMD operator O N M,n = U N,n U M,n consists of UV divergences of U N,n and U M,n , which are independent of each other. It is obvious since the semi-compact operators are separated by a space-like distance b 2 , and any interaction between them is UV finite. The UV divergence of a semi-compact operator contains also a collinear divergence from the interaction of the distant segment of the light-like Wilson line. Despite some similarity this collinear divergence should not be confused with the rapidity divergences discussed in the previous section. They have different properties (e.g. collinear divergence does not receive power corrections), and they are treated differently within the factorization theorem. On another hand, they are regularized by a common regulator (the δ-regulator in the present computation).
UV renormalization of U 1 . The renormalization factor for U 1,n is equal to the renormalization factor of the "good" component of the quark field in the light-cone gauge. It

JHEP01(2022)110
UV renormalization of U 2 . The UV divergent part of U µ 2,n is computed from the diagrams u2-u10 shown in figure 8. The computation is straightforward, and in fact, it almost coincides with the computation of the 2 → 2 evolution kernels for collinear distributions made in refs. [46,76]. The only difference is the collinear divergence (and constant contributions) due to the half-infinite Wilson line.
The collinear divergences are present in the diagrams u3 and u4. Their expressions are wherep ξ andp A are momenta of quark and gluon fields of the operator. Thus, the collinear singularity ∼ ln δ + is proportional to the cusp anomalous dimension (just alike the U 1,n case), but the regulator δ + is weighted by a different momentum in different diagrams. Different weighting produces a non-trivial contribution to the evolution kernel for TMD operators.
The computation of the remaining diagrams is straightforward. Let us only mention the diagram u2 that should be computed up to a derivative term, which after application of EOMs produces the operator U 2,n . It is convenient to present the sum of diagrams in the form where q + =p A +p ξ . The kernels H 1,2 are quasi-partonic evolution kernels for quark-gluon pair [79]. In the notation of refs. [46,80], they read (1) . (8.27) Here, the H are elementary 2 → 2 kernels that are integral operators acting in position space. In the present case, their explicit expressions are

JHEP01(2022)110
where we use the shorthand notation These elementary kernels are invariant under the SL(2)-conformal transformation. For general expressions of H, see appendix A in ref. [80]. For the corresponding expressions in momentum space, see [76].
Adding the renormalization of the quark wave function, gluon wave function and the coupling constant (the latter two compensate each other) we get the renormalization factor The sign ⊗ indicates the integral convolution between kernels H 1,2 and operator U . Here, we have eliminated ∼ ln(is) terms, due to their cancellation with the charge-conjugated part of the TMD operator. The renormalization of the operatorÛ i 2,n (defined in eq. (5.35)) is diagonal in the spinor indices where Note, that the kernel H 2 vanishes, since γ µ γ ν T γ µ T = 0. The momentum space representation for the kernel H 1 is given in eq. (C.6).
A simple structure of the renormalization factor (8.33) allows us to guess the expressions for LO renormalization of many semi-compact operators, for example gluon or di-quark twist-(2+1) operators. In these cases, one should replace H by the corresponding parton evolution kernel, replace color coefficients in the last line, and take into account a different wave-function renormalization.

Renormalization of unsubtracted TMD operators
Finally, we combine together the rapidity renormalization eq. (8.4), (8.5), and UV renormalization eq. (8.21), (8.34). For the TMD operators of twist-(1+1), eq. (5.8), (5.9) we JHEP01(2022)110 where we omit the common argument ({z 1 , z 2 , z 3 }, b) of all TMD operators. The sign ⊗ indicates that the factor Z U 2 is an integral operator, which acts on the positions the of quark-gluon pair. These expressions also hold for Fourier transformed TMD operators, with the expression for the quasi-parton evolution kernels H given in appendix C. The expressions for the TMD operators composed from anti-collinear fields are analogous with In eq. (8.36)-(8.38), we assume the forward kinematics, and q + being the momentum passing through semi-compact operators. Moreover, we set this momentum equal to the momentum passing though the EM current, since the corresponding combination (p + ξ +p + A for Z U 2 andp + ξ for Z U 1 ) equal to q + in the effective operator, due to the δ-function in eq. (6.17).
Note, that in the eq. (8.36)-(8.38) we are using the renormalized factor R, since its UV divergence is part of the UV renormalization factor for TMD operator. Alternatively, one can use the unrenormalized factor R (7.16). In this case, eq. (8.36) takes the form and similar for other operators. In eq. (8.36)-(8.38), the asterisk denotes the complex conjugation. It affects the phase in the term ∼ ln(iL) only. Due to this conjugation the complex parts of the factors Z cancel. Since we already took this cancellation into account in the definitions eq. (5.33), (5.35) where the terms ∼ ln(iL) are eliminated, the indication of conjugation is obsolete. However, we keep this indication in formulas eq. (8.36), (8.38) to keep track of terms' order.

JHEP01(2022)110 9 Recombination of divergences and scaling of TMD operators
Having computed all elements of TMD factorization (UV and rapidity renormalization factors, hard coefficient functions and the soft factor), we can finally combine them into a divergence-free expression. This procedure defines the scheme of the rapidity and UV renormalization, and thus defines the physical TMD distributions.

Subtracted version of TMD operators.
The renormalized TMD operators eq. (8.36)-(8.38) have the common form Combining the TMD operators with the soft factor, eq. (7.13) we obtain where the factors R are canceled, and , (9.4) and Here, we have introduced the notation for Lorentz invariant combinations, The variables ζ andζ naturally appear as the arguments of the logarithms in subtracted renormalization constants. Explicitly, the renomalization constants are

JHEP01(2022)110
The factor Z sub U 1 is a half of the TMD renormalization constant and coincides with it [65]. The factor Z sub U 2 is a new one, to our best knowledge. The rapidity divergent factor R cancels between the soft factor and TMD operators. 9 The leftover of the rapidity divergences are the scaling parameters ν ± . In the final definition of the TMD operator in eq. (9.5), we replace ν ± by ζ andζ using that The ν 2 = 2ν + ν − is a low-energy parameter related to the definition of soft modes, and it can be hidden in the definition of TMD distribution. The scaling variables ζ andζ satisfy The definition of the subtracted TMD operator incorporates the remnant of the soft factor S 0 , which is a finite number that depends on b 2 . We recall that the function S 0 is nonperturbative. The absorption of S 0 into TMD distributions is a part of a scheme definition for the rapidity renormalization. Effectively it adds finite terms to the "minimal subtraction scheme" factor R. The equation (9.3) serves as the definition of the scheme, i.e. the final expression for the effective operator in DY, SIDIS, and SIA does not contain any extra factors. This statement defines the commonly used physical TMD distributions. It also implies that the TMD factorization formula for processes, which contain a soft factor different from S LP , would have an extra nonperturbative function composed from the remnants of the soft factors. For example, such function occurs in the factorization theorem for quasi-TMD distribution [81].
Evolution equations for TMD operators. The TMD operator O N M,n has two scaling parameters µ and ζ. The evolution equations with respect to these parameters follow from the scaling invariance of bare TMD operator eq. (9.1).
The evolution equation with respect to µ is where we omit the factors Z (which commute with the derivative) for brevity. The second term of this expression is divergent. This divergence is compensated by the end-point divergence of ∼ n µ part of the operator J µν 2111 (5.67) at x2 → 0, which can be easily confirmed at LO. In this way, derivatives of rapidity renormalization factors cancels in-between terms related by EOMs.
The end-point divergence of ∼n µ part of the operator J µν 2111 remains uncompensated. However, this part of the effective operator is a result of the direct interaction with background field. The limit x2 → 0 pushes the gluon field to the overlap region and thus should be subtracted.
The same applies for other terms of the effective operator. The explicit realisation of this procedure will be presented in a different publication.

JHEP01(2022)110
where The symbol ⊗ indicate a possible integral convolution. Each anomalous dimension acts on corresponding arguments. The evolution equation (9.12) is general in the sense that a TMD operator of any twist-(N+M) satisfies it. The value of twist is preserved by the evolution, which is a part of geometrical twist definition. However, if there are several operators U N of the same twist and other quantum numbers they can mix with each other.
In the present case, we have deal on with twist-1 and twist-2 operators. The corresponding anomalous dimensions are (9.20) The function D(b, µ) is the rapidity anomalous dimension. It is also known as the Collins-Soper kernel (K = −2D), originally introduced in ref. [2]. For a relation among different definitions of the CS-kernel, see ref. [84]. The LO expression at small b is where in the second line the limit → 0 is taken. The rapidity anomalous dimension is known up to NNLO [73,74]. The important feature of the rapidity anomalous dimension is its nonperturbative nature. In that respect, eq. (9.21) is valid only at small values of b and gets power corrections once b increases. The ∼ b 2 power correction has been computed in ref. [85] from the analysis of the soft factor. Phenomenologically, the Collins-Soper kernel can be extracted from the analysis of TMD factorization at different scales. The most recent extraction of it can be found in refs. [11,12] (see also [85] for comparison of different extractions). The TMD evolution is given by a pair of equations, eq. (9.12), (9.19). The existence of a common solution is guaranteed by the integrability condition (known also as the Collins-Soper equation [86]). It relates UV and rapidity anomalous dimensions where the last equality follows from eq. (9.17), (9.18). The integrability condition guarantees the path-independence of the evolution, see ref. [84] for an extended discussion.
Cancellation of divergences in the factorized expression. Finally, we should check that the UV poles of the TMD operators cancel the IR poles of the hard coefficient function. For the LP term (the first line of eq. (6.17)), it implies that is finite. Indeed, the pole part ofC 1 , eq. (6.12) is where we use the relation (9.11) between ζ andζ.

JHEP01(2022)110
The relation for NLP part (the second and the third lines of eq. (6.17)) is a bit more complicated, because it involves a convolution. Using the definition of the operators eq. (5.61), (5.62) and the form of the convolution in eq. (6.17), one shows that the cancellation of poles requires that where U is a test function, and x 1 is the momentum fraction associated with the gluon. Inserting the momentum space representation for the Z U 2 renormalization constant eq. (C.6) and performing changes of variables, we confirm that This is an independent check of the computations made in section 6. We also stress that most part of the renormalization constants Z U 2 is equal to the quasi-parton-pair evolution kernel, and it can be cross-checked with the literature e.g. [46,76,79]. Altogether, these checks strongly support the results presented in this work.

Conclusion
In this work, we have developed a method to derive the TMD factorization theorem. It is based on the background field method and has similarities with the high energy expansions [87], SCET, and operator product expansion. The method, to which we refer to as TMD operator expansion, allows a systematic derivation of the TMD factorization formulas at operator level and at any order of power series. As a demonstration, we compute the TMD factorization to NLP at NLO and confirm our computation by comparison with the well-known LP and partially-known NLP expressions. The expression for effective operator is given in eq. (6.17), and evolution equations for NLP TMD operators are given in eqs. (9.12), (9.19). We recover many results found in the literature in one or another form even without reference to TMD factorization. For example, the evolution kernel for TMD operators at NLP incorporates the standard quasi-parton evolution kernels [79], which we successfully reproduce. To our best knowledge, the NLO perturbative correction to the NLP TMD factorization theorem is a new result. The main goal of this work is the formulation of a general approach to power corrections in TMD factorization. Let us remark some of the most important lessons.

JHEP01(2022)110
The TMD operator expansion is derived starting from the definition of QCD in the functional integral form. As a part the of derivation, we suppose that hadron states are built from collinear and anti-collinear fields and perform a functional integration of the remaining components. Due to it, the derivation of a factorized expression is automatic, as all operators and their coefficient functions arise from a single initial definition, avoiding any matching procedure, typical for many other derivations. The computations are done in position space, which is the natural language for power corrections since the momentum space expressions are complicated due to multiple momentum ranges. To our best knowledge, it is the first computation of TMD factorization solely in position space.
The famous process dependence of TMD factorization [24,88] (which consists in the orientation of gauge links) appears due to the boundary conditions that are imposed on the background fields. These boundary conditions follow from the request that (matrix elements of) operators are analytic in a proper part of momentum space in accordance with the process. The purely operator-level derivation of TMD factorization is an important step forward for the future development of the TMD approach. With minimal modifications our expressions can be used for a description of processes with jets, e.g. [89][90][91], or involving generalized TMD distributions (GTMDs).
The TMD operators are built from two semi-compact light-ray operators O M N = U M U N , separated by a transverse distance b. The operators U N can be sorted with respect to geometrical twists, alike ordinary compact operators. They are subject to independent UV renormalization. Therefore, the TMD operators can be classified with respect to TMDtwist, which is given by a pair of numbers (M+N), where M, N are geometrical twists of semi-compact operators U . The TMD distributions derived from operators with different TMD-twist have independent UV scaling, and thus they are independent observables. Let us mention that, in principle, operators of higher TMD-twists can mix with operators of lower TMD-twist via the rapidity scale evolution. This question is to be studied in the future.
In the small-b limit, a TMD operator of twist-(M+N) matches onto collinear operators of twist-(M + N ) or higher. This is an essential property because it ensures that the terms of TMD factorization with TMD twist-(1+1) (that is, the LP term and its kinematic power corrections) turn into the fixed order expression with twist-2 parton distributions in the high-energy limit. This observation gives hope to describe the whole q T spectrum of DY, SIDIS, and SIA processes within the TMD factorization approach.
We demonstrate in three independent ways (from the side of the soft factor, from the side of the rapidity divergences renormalization theorem, and by the direct NLO computation) that the rapidity-scale evolution is the same for LP and NLP operators. We also show that it is the same for a part of N 2 LP TMD operators (namely, for all operators of TMD-twist-(2+2) and TMD-twist-(3+1) quasi-partonic operators). The equality of rapidity evolution for different power TMD operators opens new possibilities to measure the Collins-Soper kernel, which is one of the fundamental functions in QCD [85]. For example, studying the ratios of sub-leading quasi-TMD distributions on the lattice, similar to the existing cases [92][93][94]. The important discussion on the derivatives of rapidity renormalization factors, the end-point divergences and their mutual cancellation is left for the future publication.

JHEP01(2022)110
The derivation allows an accurate separation of sources of power corrections, which are listed in the introduction. In particular, we demonstrate that the twist-(1+1)×(1+1) part of the NLP term extends the LP term and is the kinematic power correction. It restores the EM gauge invariance up the N 2 LP order and has the same hard coefficient function. Another important observation, which was overlooked or ignored in previous studies (except [18]), is that the actual expansion parameter is τ 2 T = q 2 T /|2q + q − | (= q 2 T /(Q 2 + q 2 T ) for DY ans SIA, and (= q 2 T /(Q 2 −q 2 T ) for SIDIS) rather than q 2 T /Q 2 . In fact, the value Q 2 does not appear in any formula, but all kinematic variables are expressed via τ 2 T . Accounting for this simple fact can be important for phenomenology (see some studies in ref. [11]).
In the present work, we consider only the theoretical aspects of the TMD factorization at NLP. The relevant cross-section, systematization and relations to earlier defined NLP TMD distributions [24,48] will be performed in future publications. It is known that many observables are not sensitive to NLP ∼ q T /Q corrections but have the first non-vanishing correction at N 2 LP ∼ q 2 T /Q 2 . Even so, our work will improve the present LP picture due to kinematic corrections and helps to clear up the definitions of some observables linear in q T , such as the Cahn effect [95].

JHEP01(2022)110
where A µ (B µ ) and q (ψ) are background (dynamical) gluon and quark fields. The Lagrangian L[q + ψ, A + B] is invariant under gauge transformations of the dynamical fields and where η is the ghost field, and α is residual gauge parameter. The gluon-field strength tensor is F a The effective action revealed after the integration over the dynamical fields is invariant under the gauge transformation of the background fields.
We split the background field into two non-overlapping components Then the action can be presented as where S int involves all fields, and S QCD(cov.g.) is the QCD action in the covariant gauge. In turn, the interaction term can be conveniently split into S int = S 1h + S 2h + S 12 + S 12h , (A.9) where S 1(2)h contains only fields {A 1(2) , q 1(2) , B, ψ}, S 12 contains only background fields, and S 12h contains combinations of both components of background field and dynamical fields. Since the background fields are external, one can eliminate all terms in the interaction Lagrangians that satisfy EOMs. These are all terms ∼ B 1 in S 1h and S 2h , and some terms in S 12 . The resulting actions are S 1h(2h) = g d 4 x q / Bψ +ψ / Aψ +ψ / Bq

JHEP01(2022)110
here q = q 1(2) and A = A 1(2) , The terms {1 ↔ 2} represent the previous terms in the line with the replacement q 1 ↔ q 2 and A 1 ↔ A 2 . The factor 1/2 in S 12 results from the symmetric elimination of EOMs (with respect to A 1 ↔ A 2 ). The Feynman rules for S 1q and S 2q interaction can be found in ref. [31].

B Computation of the hard coefficient function
In some cases, the computation in coordinate space is simpler than in momentum space. In particular, it is known that computations for higher-twists operators are always simpler in position space (cf. e.g. [76] and [46], section 4 and 5 of ref. [38]). This is due to the fact that momentum fractions are not sign-definite for higher twist observables, which leads to the necessity of a separate evaluation of loop integrals for each domain. The coordinate space does not have this problem, which is reflected in the different relative positioning of fields.
There are many examples of computations with the background field for collinear factorization. The pedagogical introduction into such computation is given in ref. [29], and also in refs. [27,38]. However, the case of composite background fields seems to be completely new. The main problem here is that the presence of two (or more) directions makes loop-integrals very asymmetric and thus cumbersome to evaluate. In this appendix, we present a simple technique that bypasses these problems. The main structure of the computation follows the one described in appendix B of ref. [29]. The technique is universal and can be applied to NLO (and with minimal modification to higher orders) computations of diagrams with an arbitrary number of external fields. We exemplarily show the computation of diagram 4. All other diagrams, including two-point diagrams and the computation of UV divergence in section 8.2, are computed in the same manner.

JHEP01(2022)110
Computing all other diagrams for the coefficient function of the effective operator and summing them together, we obtain eq. (6.9). The same technique, but with one light-cone direction, has been used to compute UV part of the TMD operators in section 8.2.

C Evolution kernels in momentum space
In this appendix, we present the evolution kernels eq. (8.28)-(8.31) in momentum space. In the formulas below the variables x 1 and x 2 are Fourier conjugated to z 1 and z 2 , as it is defined in eq. (5.51). Namely, z 1 and z 2 are the positions of gluon and quark fields, correspondingly. We obtain the following expressions where θ(a, b) = 1, a > 0 and b > 0, 0, a 0 or b 0, δ a = 1, a = 0, 0, a = 0. (C.5) The evolution kernels preserve the total momentum passing though the diagrams (i.e. x 1 + x 2 ). However, they do not preserve the sign of individual components. Summing together these kernels we obtain the expression for the kernels H 1 and H 2 (8.26), (8.27),

JHEP01(2022)110
In the literature, one can find a different set of variables used for the momentum fractions. Namely, x 1 = rx and x 2 = rx. The value of r = x 1 + x 2 is preserved by the evolution kernels, and drops out of the expressions. In these variables, the evolution kernels have the form (for r > 0) x y 2 (θ(x −ȳ,x) − θ(ȳ − x, −x)) U (y) + C A δ x U (0).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.