QCD factorization for twist-three axial-vector parton quasidistributions

The transverse component of the axial-vector correlation function of quark fields is a natural starting object for lattice calculations of twist-3 nucleon parton distribution functions. In this work we derive the corresponding factorization expression in terms of twist-2 and twist-3 collinear distributions to one-loop accuracy. The results are presented both in position space, as the factorization theorem for Ioffe-time distributions, and in momentum space, for the axial-vector quasi- and pseudodistributions.


Introduction
Twist-three effects originate from the quantum-mechanical interference between a single parton and a gluon-parton pair. Their investigation is exceptionally important for the understanding of the nature of strong interactions, but it is also a very challenging task. The main complication is that twist-three contributions are usually of subleading power in the hard scale, and either contaminated by the leading-twist contributions or statistically suppressed. Additionally, a typical twist-three observable is sensitive only to a certain projection of the underlying quark-antiquark-gluon correlation function and a combination of many observables is required to unravel its structure.
More recently, the interest in twist-three effects was fueled by studies of transverse momentum dependent (TMD) parton distributions. In this case, the twist-three distributions arise in the collinear limit for many TMD distributions already at the leading power [34][35][36][37]. In particular, the single transverse spin asymmetry is sensitive to the Sivers function [38], which matches the Qiu-Sterman function [39] in the collinear limit. The latter is a certain integral of the same quark-antiquark-gluon correlation function that contributes to g 2 (x, Q 2 ), see, e.g., [40][41][42]. In this way, one can extract twist-three distributions from the transverse momentum dependence, which has been done recently for the Qiu-Sterman function [43,44].
Lattice QCD simulations have the potential to explore the plethora of twist-three distributions more directly by the calculation of Euclidean correlation functions that are specially tailored to the twist-three effects. The lattice calculations of the leading-twist parton distributions (PDFs) have attracted a lot of attention, see [45,46] for a review. In this case, a convenient object is an equal-time correlation function of two quark (or gluon) fields connected by the straight Wilson line, which can be factorized [47] in terms of the corresponding PDF convoluted with a perturbatively calculable coefficient function. The lattice "data" can be analyzed either directly in position space in terms of Ioffe-time distributions [48], or Fourier-transformed to momentum space where they are referred to as quasidistributions (qPDFs) [49] or pseudodistributions (pPDFs) [50].
In Refs. [51,52] the same approach was suggested to extract the twist-three PDF g T (x) [53] from a lattice calculation of the correlation function N (p)|q(z)γ µ γ 5 [z, 0]q(0)|N (p) , (1.1) where z µ is a space-like vector and [z, 0] is the Wilson line. We refer to this object as a "quasidistribution" in a broad sense, to distinguish from the parton distributions defined for the light-like separation z 2 = 0. In this work, we derive the factorization formula for the correlation function in (1.1) to twist-three accuracy at next-to-leading perturbative order (NLO). The results are presented both in position space, as a factorization theorem for Ioffe-time distributions, and in momentum space, for the axial-vector quasi-and pseudodistributions.
Our main results can be summarized as follows: • The factorization theorem to twist-three accuracy has a more complicated structure as compared to the leading twist [47]. In particular the quasidistribution (1.1) cannot be factorized in terms of the parton distribution g T (x) [53], as conjectured in [51,52].
• We have checked that the factorization scale dependence of the coefficient functions in our expressions agrees with the known results on the renormalization of twist-three operators, cf. [5,7,12,41,54].
• The simple LO relation between the twist-two contributions to the quasidistributions (1.1) with longitudinal and transverse projections of the vector index is violated at NLO. Hence the Wandzura-Wilczek relation is violated for quasidistributions. This is different from DIS, in which case the Wandzura-Wilczek relation for twist-two contributions holds on the level of structure functions.
• In the limit of large number of colors, the logarithmic terms in the coefficient functions are simplified [10,18] and can be combined to express the result in terms of the quarkantiquark transverse spin distribution g T (x). However, this simplification does not occur for finite terms.
• We consider a short-distance expansion of (1.1) which may provide a method to calculate the matrix element of the twist-three operator of the lowest dimension avoiding a complicated procedure of nonperturbative renormalization in the presence of power divergences, cf. [31].
The presentation is organized as follows. Sect. 2 and Sect. 3 are introductory and contain main definitions and notations. In Sect. 4 we formulate the factorization theorem and derive the NLO expressions at the operator level and for the position-space (Ioffe-time) distributions. In Sect. 5 the corresponding results are presented for the qPDFs and pPDfs. The final Sect. 6 contains a discussion and outlook. Technical details are delegated to the Appendices.

Definitions
In the present work we study the twist expansion of a product of quark and antiquark fields where z µ is a four-vector, q are quark fields, and [z, 0] is the straight Wilson line in the fundamental representation of the gauge group For space-like separations z 2 < 0 which are relevant for lattice calculations, the timeordering in Eq. (2.1) is redundant. The Dirac matrix γ 5 is defined (in d = 4 dimensions) as where is the Levi-Civita tensor with 0123 = − 0123 = 1. Flavor indices of the quark fields are omitted for brevity.
If not stated otherwise, we tacitly assume that the nonlocal operator (2.1) is renormalized, The renormalization constant Z O in the MS scheme is known to three-loop accuracy [55,56]. It is the same for all Dirac structures. Both sides of Eq. (2.4) depend on the renormalization scale µ R . This dependence is not shown explicitly in what follows, unless it is important for understanding. In renormalization schemes with an explicit regularization scale, the Wilson line in Eq. (2.1) suffers from an additional linear ultraviolet divergence [57] which has to be removed. This can be done by the renormalization of a residual mass term, similarly to the heavy quark effective theory, or, alternatively, by forming a suitable ratio of matrix elements involving the same operator [58,59]. This issue is well known and does not require an extra elaboration.
In this work, we study only the forward matrix elements. Thus, the global positioning of the operator is unimportant, and only the distance z between the fields plays a role. Without loss of generality, we can put the quark position at the origin. The nucleon matrix element of the operator in Eq. (2.1) can be parametrized in terms of three invariant functions, Here, M 2 = p 2 is the nucleon mass, s µ is the spin vector, (s · p) = 0, s 2 = −1, and We refer the functions G as Ioffe-time quasidistributions (qITDs) [58,60,61], and the variable as the Ioffe time [60]. The subscript T indicates the projection onto the transverse plane, which is orthogonal to z µ and p µ : (2.9) The transverse projection of the metric tensor in four dimensions satisfies For further use we also introduce the anti-symmetric tensor If the four-vectors z µ and p µ are positioned in the (t, z)-plane, then 12 T = − 21 T = 1 (and other components vanish). Let us note that the terms ∼ z 2 p 2 can be omitted, since in the following we study the limit z 2 p 2 /p 2 z → 0. However, we keep these terms in (2.9) and (2.11) to ensure exact orthogonality conditions.
The three qITDs G 1 , G T and G 3 at small z 2 can be matched to collinear PDFs. This matching is studied in great detail for G 1 , where only the twist-two collinear PDF contributes at accuracy O(z 2 ). The corresponding coefficient function is known to two loops (NNLO) [62][63][64]. The factorized expressions for G T and G 3 are more complicated. At the same power accuracy, they receive contributions from collinear PDFs of twists 2,3 and twists 2,3,4, respectively. In this paper we derive the leading-power factorization theorem for the qITD G T distribution to NLO (one-loop) accuracy. It incorporates twist-2 and twist-3 collinear distributions. Starting from this result, the coefficient functions for qPDFs and pPDFs can be obtained by the appropriate Fourier transform. The necessary definitions are given in the corresponding sections.

Light-ray operators
Operator product expansion (OPE) can be conveniently organized in terms of the generating functions of (renormalized) local operators [7,[65][66][67][68][69] The superscript µ F indicates the renormalization scale for the operator. In what follows, we often suppress the scale dependence not to overload the notation. Importantly, the operator O γ µ γ 5 (z, 0) in Eq. (3.1) and O γ µ γ 5 (z, 0) in Eq. (2.1) are different beyond the tree level. Local operators on the r.h.s. of Eq. (3.1) can be decomposed into a sum of contributions with different geometric twist (dimension minus spin). In particular, the twist-two operators are obtained by symmetrization over all Lorentz indices and subtraction of traces. We define the twist-two projection of the nonlocal operator (3.1) as the generating function for (renormalized) local twist-two operators where n µ is an auxiliary light-like vector, n 2 = 0, and O and can be written more explicitly in several different representations [7,[67][68][69][70]. Similarly, it is possible to construct the projectors onto twist-three and higher-twists operators.
In this work we consider the transverse component O γ µ T γ 5 (z, 0) = g µν T O γν γ 5 (z, 0) and neglect twist-four contributions. Hence, we can omit terms O(z 2 ) corresponding to the subtraction of traces. To this accuracy it is sufficient to use the following decomposition where [68] [ Here and below we use the notation T µ + (az, bz, cz) = igq(az) γ ρ / zγ µ γ 5 F ρz (bz)q(cz) (3.8) where g is the QCD coupling constant, F µν is the gluon field-strength tensor, F µz = F µν z ν , and F µν = 1 2 µναβ F αβ . The Wilson lines connecting the three fields in the quark-antiquarkgluon operators (3.8), (3.9) are implied. Let us emphasize that Eqs. (3.6) and (3.7) are exact operator identities for renormalized operators. We also note that the T µ ± operators involve the QCD coupling, hence [O γ µ γ 5 (z, 0)] tw3 vanishes in the free theory. So, the twist-three effects arise due to quark-gluon interactions.
The position-space PDFs S ± are related to the PDFs in momentum fraction space by the Fourier transformation where The functions S ± obey the symmetry relation [11,18] and similar for S ± . It follows that S − (ζ, βζ, αζ) = S + (ᾱζ,βζ, 0), and, therefore, the two terms in (3.13) are equal to each other. Thus we write simply and, going over to the momentum fraction representation, obtain (3. 18) In this expression one can get rid of two integrations using the delta-functions, but the resulting expression is unwieldy due to a multitude of integration regions. For comparison, the corresponding expression in position space (3.17) is rather compact. This situation is rather general. For this reason, we carry out the major part of the calculations in position space, and go over to momentum fractions only in the very end.
The relation x 1 + x 2 + x 3 = 0 (3.15) implies that the twist-three PDFs are functions of two variables. However, the three-variable notation used here is more convenient for many reasons. First, it simplifies the symmetry relations (3.16). Second, twist-three correlation functions have different parton interpretation in each kinematic domain x i ≶ 0 [71] and can most naturally be presented using three-component barycentric coordinates, see [41,42].
There is no single established notation for twist-three PDFs. Our definition is common in the literature on the DIS structure function g 2 (x, Q 2 ), e.g., [4,18,72]. In applications to semi-inclusive reactions [36,37,[40][41][42]73], it is customary to use a different pair of twist-three PDFs defined as The corresponding momentum fraction distributions T and ∆T (defined as in (3.14)) are related to S ± as The T and ∆T functions satisfy the following symmetry relations A detailed comparison between different notations can be found in Refs. [41,42]. The correlation functions S ± (or, equivalently, T and ∆T ) are scale-dependent. Their evolution is autonomous, in the sense that it does not involve other distributions (apart from the gluon twist-3 distribution in the singlet case, which we do not discuss in this work). The evolution equation takes the form where K is a kernel which, in general, depends on six variables x = {x 1 , x 2 , x 3 } and y = {y 1 , y 2 , y 3 }. The leading-order (LO) expression for the kernel K can be found in [41,54]. This expression is, again, relatively compact in position space but rather lengthy in terms of the momentum fractions. For completeness we present the position-space kernel [41] in App. A. On the contrary, the evolution of the function g T is not autonomous, since the projection of variables (x 1 ,x 2 ,x 3 ) onto a single momentum fraction x is not an eigentransformation of the evolution kernel (3.23), except for the large-N c regime (and in LO only), where the function g T obeys the DGLAP-type evolution equation [10,18].

QCD factorization for quasidistributions
At tree level O γ µ γ 5 (z, 0) = O γ µ γ 5 (z, 0) and therefore where the ellipses stand for higher-twist O(z 2 ) corrections. Including higher-order perturbative corrections these expressions can be generalized to the following factorization theorems: where the coefficient functions C k are given by a series expansion in the QCD coupling We write with tree-level expressions The validity of QCD factorization for qITDs G 1 , G T is a direct consequence of the existence of the operator product expansion. The factorization theorems for quasidistributions in momentum space are obtained by the Fourier transform of above expressions and do not require any additional justification. The corresponding expressions for qPDFs and pPDFs are derived in Sec. 5.
The main purpose of this work is the calculation of coefficient functions C k to one-loop accuracy. The calculation of higher-twist contributions can be conveniently done using the background field technique pioneered by Schwinger [74] and adapted to the applications in nonabelian gauge theories in Refs. [7,42,68,[75][76][77][78][79]. The basic idea is explained briefly in what follows 1 . Explicit examples of calculations using two different versions of this approach are presented in the Appendix C.
Thanks to the exact relations in (3.12), (3.13), the twist-two contribution to G T in Eq. (4.5) can be rewritten in terms of g tw2 T , and a part of twist-three contributions in terms of the two-particle distribution g tw3 T : One should have in mind that the separation of the two-particle and three-particle contributions of twist three is not unique. We discuss this possibility and its limitations in Sect. 4.4.

Background field technique
The separation of coefficient functions and operator matrix elements in a certain amplitude can be understood in the spirit of Wilson's approach to the renormalization group as integrating out the high-momentum degrees of freedom. We introduce the scale µ and define "fast" and "slow" fields as modes with momenta p > µ and p < µ, respectively: where ψ and B are the "fast", and q and A are the "slow" components. For any gaugeinvariant operator O one can integrate over the "fast" fields giving rise to an effective operator that only depends on "slow" degrees of freedom where S is the QCD action. In this expression q and A can be considered as given fields which satisfy classical QCD equations of motion (EOM). The background field technique [75,76] is the method to evaluate such integrals paying due attention to gauge invariance. The presence of background field modifies the structure of the gauge-fixing term in the action. The analog of the conventional covariant gauge fixing condition for the "fast" fields is [75,76] ( is the covariant derivative in the background field. The expression for the action S[q + ψ, A + B] − S[q, A] in the background field gauge can be found in Ref. [75]. The major advantage of the background field method in our context is that the functional integral (4.10) is invariant under the local gauge transformations of "slow" (classical) fields. This observation greatly simplifies the calculation as it implies that one can use any suitable gauge for the background fields, and in particular a "physical" gauge where the gluon field can be expressed directly in terms of the strength tensor.
There are several methods to evaluate functional integrals in background fields, which have their advantages and drawbacks. We have performed the whole calculation using two techniques. The first approach, hereafter referred to as method A, uses traditional perturbation theory for the background field action [75] and axial gauge z µ A µ (x) = 0 for the classical field. The second approach, method B, uses the expressions for the lightcone expansion of the quark and gluon propagators in the background field [7] and Fock-Schwinger gauge x µ A µ (x) = 0. Intermediate expressions in these two methods have a very different structure. The final results are in agreement, which provides a strong check of their correctness. For pedagogical purposes we present a detailed calculation of a certain gauge-invariant subset of diagrams using both approaches in App. C.

Renormalization factors and treatment of γ 5
We use the dimensional regularization with d = 4 − 2 and (modified) minimal subtraction scheme. Calculating the relevant Feynman diagrams in the presence of background fields we obtain the expression for the bare qPDF operator (2.1) in terms of the bare light-ray operators. The result has the following schematic structure where the bare coefficient functions depend on and are singular at → 0. The 1/ terms in the coefficient functions are due to ultraviolet (UV) and infrared (IR) singularities and are removed by the renormalization procedure. In the present case, it is not necessary to distinguish UV and IR poles during the calculation. The UV singularity is removed by the (multiplicative) renormalization of the qPDF operator, Eq. (2.4). To one-loop accuracy [81] where C F = (N 2 c − 1)/2N c is the quadratic Casimir operator in the fundamental representation.
The IR singularities in the coefficient functions are removed by the renormalization of light-ray operators. For twist-two [7] 14) The coefficient of 1/ in this expression is the renown DGLAP kernel in the position-space representation. The renormalization factors of twist-3 operators T µ ± have similar structure. The necessary expressions can be found in [41]. For readers' convenience we collect them in App. A, see Eq. (A.2).
between the coefficient functions and the operators in each term in Eq. (4.12), we take the limit → 0 and obtain the final result Here Z i stands for the renormalization kernel for the appropriate operator. This expression involves two scales. The dependence on the renormalization scale µ R is governed by the anomalous dimension (4.13). The dependence on the factorization scale µ F is canceled between the coefficient functions and light-ray operators. In the following we set The final comment concerns the treatment of γ 5 . The usual MS scheme is defined in such a way (see, e.g., [82][83][84]) that the renormalization of flavor-nonsinglet vector and axialvector operators is assumed to be the same. This is achieved by starting with a suitable γ 5 definition in d = 4 − 2 dimensions [85][86][87] aided by a finite renormalization that effectively restores the anticommutation property {γ µ , γ 5 } = 0. This procedure is mandatory for flavor-singlet operators, but whether it has to be followed for flavor-nonsinglet operators depends on the method how the calculation is done. In our case we can simply assume {γ 5 , γ µ } = 0 as γ 5 does not appear in traces in which case using anticommutation property leads to algebraic inconsistencies.

One-loop coefficient functions for light-ray operators
In this section we present the results for bare coefficient functions of twist-two and twistthree operators. We have done the calculations using two versions of the background field technique and got the same result, see App. C for the details. Below we write the results for the coefficient functions separating an overall factor We also omit the "bare" superscript in order not to overload the notation. Twist-two contributions: Twist-three contributions: where [O γ µ T γ 5 (z, 0)] tw3 is given by Eq. (3.7). Note that the contributions in the first two and the last two lines in Eq. (4.20) correspond to different ordering of the fields on the light cone: In the first two lines the gluon field is in between the quark and the antiquark, and in the last two lines the gluon is outside (to the left or to the right). The "wrongly ordered" contributions are suppressed in the large-N c limit.
Renormalized coefficient functions are obtained as explained above, applying the renormalization factors for the light-ray operators and the overall renormalization factor Z O for the qPDF operator. The resulting expression is finite at → 0, as it should. Cancellation of the 1/ poles provides one with a further check of the calculation.

qITDs at NLO
The qITDs are obtained by taking the matrix element of the OPE (4.15) such that the structure of the expansion is essentially retained. In the following expressions 21) and the plus-distribution is defined as usual, We obtain where This result is in agreement with earlier calculations [47,59,88]. Note, that the coefficient function (4.24) is nothing but the renormalized C 1 given in Eq. (4.19).
The transverse qITD G T receives twist-two and twist-three contributions where the twist-two term is given in Eq. (4.26) and the twist-three term is given in Eq. (4.31). First, we present the twist-two part, which reads where C (1) and we remind that (3.12) 1 (α, L z ). As a consequence, the Wandzura-Wilczek relation does not hold for quasidistributions, This is different from the "classical" Wandzura-Wilczek relation for polarized DIS [2], which is exact at the level of structure functions, cf. [18]. The twist-three part deserves more attention. The first term in Eq. (4.20) is written in terms of the two-particle twist-three operator [O γ µ T γ 5 (z, 0)] tw3 with the same coefficient function as the twist-two contribution. If this "two-particle" term were the only one, the qITD G T would be factorizable in terms of the "full" PDF g T (ζ) = g tw2 T (ζ) + g tw3 T (ζ), as conjectured in Refs. [51,52]. To clarify the role of the remaining three-particle contributions it is instructive to consider the large-N c limit such that the contributions in the last three lines in (4.20)) can be dropped and the expressions become much simpler. We have  To simplify this expression we have used Eq. (3.17) and the symmetry relation (3.16) which equalizes the contribution of T µ + and T µ − . We see that the three-particle contributions involving the weight factorβ in the integral over the gluon position can be rewritten in terms of the two-particle twist-three PDF g tw3 T (αζ; µ). However: 1. With this addition, the coefficient function of g tw3 T (αζ; µ) becomes different from the coefficient function of g tw2 T (αζ; µ). This is already true for the logarithmic term ∼ L z , i.e. for the evolution kernel in the large-N c limit [10,18]. 3 The main effect to the LO accuracy is the constant shift by N c of the anomalous dimensions of the twist-three operators as compared to the leading twist ones [10].
2. A genuine three-particle contribution involving ln β appears, see the last line of Eq. (4.30). It cannot be rewritten in terms of g tw3 T (αζ; µ).
The 1/N c corrections have a more complicated structure and do not allow a separation of two-particle contributions in a natural way. Collecting all terms we obtain the NLO expression for the twist-three part of the qITD G T (ζ, z 2 ) (cf. (4.8)) where the logarithmic part is given by 5 Quasi-and pseudo-distributions qPDFs and pPDFs [47,49,90] are defined as functions of the parton momentum fraction so that their interpretation is more close to the traditional PDFs. There seems to be no established notation in the literature for all their variants. Traditionally, the distributions and structure functions related to axial-vector operators are denoted by the letter g with various subscripts, cf. [3,4,53]. To follow this practice, on the one side, and to distinguish various types of distributions, on the other side, we use different fonts: The qPDFs are denoted by the typewriter font letter g, and pPDFs are denoted by the blackletter font g, and similarly the corresponding coefficient functions. Following [47,49] we introduce qPDFs as Fourier transforms of the qITDs with respect to the distance z. The orientation of the vector z µ is kept fixed. We define where v µ = z µ /|z| is the unit vector along z µ , and p v = (p · v). In turn, pPDFs [50] are defined as Fourier transforms of qITDs with respect to the momentum p keeping its orientation fixed, These distributions have different properties. In particular, pseudodistributions have a natural "partonic" support |x| < 1 [90], whereas qPDFs have |x| < ∞.
In this section, we present the NLO expressions for axial-vector pPDFs and qPDFs. Going over from the qITDs (4.23), (4.31) to the pseudodistributions is relatively straightforward, since all Fourier integrals are reduced to Dirac delta-functions and are easily taken. The derivation of qPDFs is less trivial due to Fourier transformation of logarithmic contributions. These terms have "unnatural" support properties, which makes the direct calculation cumbersome. To avoid this complication we use the identity from which simple relations between the coefficient functions for qPDF and pPDFs can be derived, see App. B. Using (B.9) (B.10) we are able to obtain the NLO expressions for qPDFs from the corresponding pPDFs with relatively little effort.

pPDFs at NLO
Using Eqs. (4.23), (4.31) and performing the Fourier transformation (5.3), (5.4) we obtain where C (1) and with the logarithmic part being (5.10) One can show that the integrands in these expressions are finite at x i → 0 and α → 1. Moreover, the first and the second moments of (5.9) and (5.10) vanish. All expressions are defined for −1 < x < 1.

qPDFs at NLO
The qPDFs can be most easily derived from the corresponding pPDFs (5.6), (5.7) applying the double-Fourier transformation (5.5). This transformation is non-trivial only for the terms involving L z (B.10), which are also the ones responsible for extending the support property of qPDFs beyond the partonic region |x| < 1 to |x| < ∞. In the expressions given below we use the following notation: The first moment of both distributions is zero.
We obtain (5.15) and the coefficient functions are given by where Note that in the large-N c limit the genuine three-particle contributions for qPDFs and pPDFs are the same.
Our results for the coefficient functions C 1 and C T coincide with known expressions [47,88]. The remaining expressions are new results. In general, the structure of the threeparticle contributions to qPDFs is rather unwieldy due to multitude of integration regions. For two-particle contributions, as well known, three domains x > 1, 0 < x < 1 and x < 0 are distinguished. For three-particle contributions one ends up with 30 domains for the variables (x 1 , x 2 , x 3 ). This structure is not explicit in Eq. (5.19) but is revealed once the integrations over the delta-functions are performed.

Discussion and outlook
We have formulated the factorization theorem for the space-like axial-vector correlation function (2.1) in terms of parton distributions to twist-three accuracy and calculated the corresponding coefficient functions to NLO accuracy. In this section we discuss the results in connection with possible lattice calculations.
Since the twist-two contributions to the "transverse" part of all versions of the quasiparton distributions are given by the helicity PDF ∆q(x, µ), the utility of the lattice approach crucially depends on the possibility to subtract (or at least minimize) such terms and reveal the twist-three contributions of interest. It is well-known that this subtraction can be implemented exactly for the correlation function of two vector currents by virtue of the Wandzura-Wilczek relation. For the quasidistributions the cancellation is not complete, as explicitly demonstrated by our calculation. We obtain for the qPDFs and pPDFs, respectively. The twist-two remainder on the r.h.s. of this relation for the qPDF case is unfortunately rather large.
To illustrate this point, consider the first nontrivial moment of g T which can be accessed by considering the small-distance expansion of the qITD. Using Eqs. (4.29) and (4.31) one obtains where (we use the notations of Ref. [25]) For an estimate, we take a 2 0.05 for the u-quarks in the proton at µ 2 = 4 GeV 2 [91,92]. Then (20/9)C F a s a 2 ∼ 3.6 · 10 −3 which is of the same order as the expected size of the twistthree matrix element |d 2 | ∼ (1÷5) · 10 −3 [25,[29][30][31]. Reducing this "twist-two pollution" can pose a serious problem for the qPDF approach in the studies of twist-three effects 4 . As far as the twist-three contribution itself is concerned, constraining the quarkantiquark-gluon correlation function in its full complexity from present-day lattice calculations is probably unrealistic. Thus trying to reduce the nonperturbative input to a function of one variable, g T (x), as attempted in [51,52], is certainly logical. However, the shortcomings of such a reduction have to be clearly understood. Any approximation of this kind is theoretically self-consistent if and only if it is maintained at all scales, in other words if g T (x) does not mix with the "genuine" three-particle contributions that are neglected. This condition is, indeed, satisfied to LO accuracy in the large N c limit [10,18], which can be sufficient at the current stage. This decoupling does not mean, however, that the coefficient functions of g T (x) to logarithmic accuracy can be calculated from the two-particle quark-antiquark matrix elements. The quark-antiquark-gluon matrix elements must be considered and contribute to the splitting functions (and to finite terms). As the result, the coefficient functions of the twist-two and twist-three contributions to g T (x) in the factorization theorem for the qPDFs are different already in the large N c limit, see (5.18). This difference is missed in [51,52]. Another issue is that at NLO finite corrections ∼ N c appear that cannot be reduced to g tw3 T . Whether such terms can be minimized in some way, remains to be studied.
To conclude, we have presented the first NLO analysis of axial-vector quasidistributions of the nucleon to the twist-three accuracy. The same method can be extended in a straightforward manner to chiral-odd twist-three quasidistributions that are of particular interest, cf. [93]. We plan to consider them in a separate publication.

Acknowledgments
This study was supported by Deutsche Forschungsgemeinschaft (DFG) through the Research Unit FOR 2926, "Next Generation pQCD for Hadron Structure: Preparing for the EIC", project number 40824754. Y.J. also acknowledges the support of DFG grant SFB TRR 257.

Appendices A Evolution kernel for twist-3 distributions
The evolution equations for twist-three quark-antiquark-gluon distributions can be found in Ref. [41,54]. For the readers' convenience, we collect the relevant expressions in this appendix.
The evolution equation for the function S − has the form The evolution equation for S + is obtained trivially, using the symmetry relation (3.16). The corresponding expression for momentum space distributions is lengthy as the evolution kernels in different sectors x i ≶ 0 are not the same. Explicit expression can be found in [54]. The distribution g T is proportional to the integral (3.17). Application of the operator H to this expression yields The ∼ N c part of this expression that includesβf (α) in the integrand, can be rewritten as a convolution of a two-particle kernel with g T with the proper rescaling of ζ. The 1/N c contributions cannot be simplified in this way and contribute to the "genuine" three-point part, see Sect. 4.4.

B Relation between the coefficient functions for qPDFs and pPDFs
The qPDF g(x, p v ) is related to the pPDF g(x, z 2 ) by the double Fourier transformation The pPDF on the r.h.s. of this equation can usually be presented in the form of the Mellin convolution where Q(x) (with −1 < x < 1) is some parton distribution function. Applying the transformation (B.1) to (B.2) one obtains where the coefficient function C(x, L p ) is related to C(α, L z ) as In these formulas we denote (B.5) Let us emphasize the factor |y| in the L p that is present in the quasidistribution (B.3). It appears due to the change of variables. Also we observe that the quasi-distribution is nonvanishing for −∞ < x < ∞.
The coefficient function C(α, L z ) depends on ζ only via the logarithms which appear in increasing powers in higher orders of the perturbative expansion. The nontrivial part of the transformation (B.4) consists of the evaluation of the Fourier transform of ln n ζ. Using that L z = l p + ln(ζ 2 e 2γ E ) we write the coefficient function C as a series C (α, L z ) = ∞ n=0 ln n (ζ 2 e 2γ E )C n (α, l p ). (B.6) Here C n is given by a perturbative series starting from a n s . The coefficient function for the qPDF is then The integral (B.8) can be evaluated explicitly for any given n. To this end one should first move the integration contour over α to the complex plane, evaluate the integral over ζ, and finally close the α-integration contour on the branch cuts of the integrand. The result for arbitrary n is simple but lengthy. To the NLO accuracy we only need n = 0 and n = 1: and Derivation of Eq. (B.10) is straightforward under the assumption that the integrand is regular for y ∈ [0, 1]. However, C(y) has a singularity at y → 1 in the form of plus distributions. These singularities result in additional pole terms ∼ δ(x) that are easy to miss in a direct evaluation. To bypass this complication, we have used that the first Mellin moment of plus distributions vanishes and thus these additional contributions can be found by enforcing vanishing of the first moment of C(x) by adding suitable ∼ δ(x) terms.

C Sample calculations
For pedagogical purposes we present in detail a part of the calculation. We used two independent methods which we call Method A and Method B. Method A is based on explicit expansion of the background field Lagrangian with subsequent evaluation of the diagrams. Method B uses partially integrated action with the propagators in the background field. Although both methods are based on the same concept, the actual calculation is rather different and in particular the diagrammatic decomposition of the relevant contributions is not the same, although there are correspondences between classes of diagrams. For illustration, in Method A and Method B we also employ a different gauge fixing condition for the background field, leading to essentially different algebra and intermediate expressions.
Naturally, the final results coincide. In this appendix, we present the calculation of the subset of diagrams which involve gluon exchange between the quark and the antiquark, see Fig. 1. This proves to be the most cumbersome part. In what follows we calculate the relevant contributions using both methods and find the same result.

C.1 Method A
Interaction with the background field is described by the action where q and A are classical background fields, ψ and B are quantum fields, f ABC are the SU (N c ) structure constants, and v µαβγ = 2g µβ g αγ − g µα g βγ − 2g µγ g αβ . (C. 2) The dots indicate terms that are not needed in the present computation.
Using this effective Lagrangian we derive the expressions for the diagrams shown in Fig. 1 where we have already done the color algebra, and are the quark and gluon propagators in position space (with d = 4 − 2 ) in Feynman gauge. We use the axial gauge for the background gluon field With this choice that there is no Wilson line, [z, 0] = 1l, and Note that here we use the retarded prescription to fix the residual gauge dependence, A µ (t → −∞) = 0. In the case of the advanced prescription, A µ (t → +∞) = 0, the lower limit of the integration would become +∞z 0 . In the following we tacitly assume z 0 > 0.
The diagrams are to be evaluated in the limit z 2 → 0 neglecting terms O(z 2 ). A similar computation is made in Ref. [42] for a TMD operator which is different from the present case only by z µ → z µ T . A very detailed description of the calculation of the (analog of) diagram D (d) can be found in App. B of Ref. [42]. Here we consider D (a) , which is algebraically simpler but incorporates all the same principal steps. In Method B, the same contribution is given by a certain mixture of D (a) and D (b) .
Thus we start with (cf. Eq. (C.19)) The expansion of the integral in powers of z 2 can be done with the following trick. We join the denominators using auxiliary integrations with Feynman parameters, and perform a shift of integration variables such that the denominator acquires the form A + z 2 B where A does not depend on z. We obtain where λ = αβ + αγ + βγ and [dαdβdγ] = dαdβdγδ(α + β + γ − 1). Next, we expand the quark fields at x, y → 0. The resulting integrals are of the form x µ 1 ...y µn /[A + z 2 B] 5−3 and can be easily taken.
To the required twist-three accuracy the diagram should be evaluated up to terms ∼q∂q so that we need to expand the fields in Eq. (C.8) to first order. Computing the (tadpole) loop integrals one ends up with (C.10) To present the expression (C.9) in this simple form we have used that z µ T = 0 and {γ 5 , γ µ } = 0 to simplify the algebra, and also changed to the dual Feynman variables (see Eq. (B.8) in Ref. [42]).
The expression in the parenthesis in (C.9) gives rise to a genuine twist-three contribution. It can be rewritten in terms of quark-antiquark-gluon operators using QCD equations of motion. For example, where T is defined in Eq. (3.8).
The remaining contributions in Eq. (C.3) are computed in the same way. We obtain dβ αᾱ T + (z, βz, αz) + T − (z, βz, αz) (C.14) where we have made a total shift of the operator position and performed a series of changes of variables to present the expression in a simpler form. Summing up everything, we observe that the gauge-dependence due to the choice of the "retarded" integration limit −∞ cancels, and the final result for this set of diagrams reads This expression coincides with Eq. (C.27) after appropriate change of variables. We have checked that the same result is obtained using the advanced axial gauge. Altogether, these different versions of the calculation provide an excessive check of the result.

C.2 Method B
In this appendix we explain the calculation of the same contribution, with a gluon exchange between the quarks, using the technique of Ref. [7]. In this case we have only one Feynman diagram where the Wick contractions correspond to the quark and gluon propagators in the background field [7] q a (x)q b (y) = i where ∆ = x − y, σF = σ αβ F A αβ t A ab , and For bookkeeping purposes it is convenient to split the calculation in the contributions with and without the background gluon field in the propagators, similarly to Fig. 1. The contribution (a) is by far the most complicated one; let us consider it in some detail.
The basic idea is that the remaining (classical) fields in the integral can be expanded as q(y + ux) = q(ux) + y ξ ∂ ξ q(ux) + . . . (C.20) after which the x-integration becomes trivial, and it is easy to convince oneself that the terms with more than one derivative produces corrections ∼ O(z 2 ) and can be dropped. Wilson lines also have to be expanded and to this end using Fock-Schwinger gauge for the classical gluon field x ξ A ξ (x) = 0, A ξ (0) = 0, proves to be very convenient. In this gauge [0, y + ux] = 1l and Taking the x-integral, we repeat the same procedure to do the y-integration: combine the remaining two propagators introducing another Feynman parameter, shift the integration variable, and expand the classical fields along the z-direction (the antiquark field and the The last term can be handled by rewriting ∂ ξ = D ξ + igA ξ , using (C.21) and the operator identity Here F zξ = z η F ηξ . In this way we obtain (d = 4 − 2 ) dt tq(vz) γ ρ γ µ T − γ µ T γ ρ F zρ (tz)/ zγ 5 q(uz) ,