Target Mass Effects in Parton Quasi-Distributions

We study the impact of non-zero (and apparently large) value of the nucleon mass $M$ on the shape of parton quasi-distributions $Q(y,p_3)$, in particular on its change with the change of the nucleon momentum $p_3$. We observe that the usual target-mass corrections induced by the $M$-dependence of the twist-2 operators are rather small. Moreover, we show that within the framework based on parametrizations by transverse momentum dependent distribution functions (TMDs) these corrections are canceled by higher-twist contributions. We identify a novel source of kinematic target-mass dependence of TMDs and build models corrected for such dependence. We find that resulting changes may be safely neglected for $p_3 \gtrsim 2M$.


Introduction
The parton quasi-distributions (PQDs) Q(y, p 3 ) recently proposed by X. Ji [1] convert into usual twist-2 parton distribution functions (PDFs) f (y) when the hadron momentum p 3 tends to infinity. Unlike PDFs that are defined through a correlator of quark fields separated by a light-like interval z, the definition of QPDs refers to the interval z that has only a space z 3 component. This opens a possibility to extract PQDs from Euclidean lattice gauge calculations.
It is expected that, for a finite p 3 ≡ P, the difference between Q(y, P) and f (y) is explained by the higher-twist and targetmass corrections in powers of Λ 2 /P 2 and M 2 /P 2 , respectively.
The target-mass dependence of the twist-2 matrix elements is well-known since mid 70's [2,3]. In Ref. [4] (see also [5,6,7]), this information was used to connect x n moments of PDFs f (x) and y n moments of the twist-2 part of QPDs Q(y, P). In Ref. [8], this connection was converted into a direct relation between the functional forms of Q twist−2 (y, P) and f (y). In the present paper, we give our derivation of this relation and emphasize that for y > 0 and y < 0 components of f (y), it reduces to a simple rescaling by factors depending on the ratio M 2 /P 2 .
We also observe that the P-evolution pattern exhibited by the corresponding components of Q twist−2 (y, P) is rather different from the nonperturbative evolution of PQDs Q(y, P) in the models considered in our recent paper [9]. Furthermore, the comparison of the two cases indicates that the M 2 /P 2 target-mass corrections in Q twist−2 (y, P) are much smaller than the Λ 2 /P 2 higher-twist corrections in our model PQDs Q(y, P).
According to Ref. [9], the PQDs are completely determined by the transverse momentum dependent distributions F (x, k 2 ⊥ ). Thus, our next goal is to find the twist-2 part F twist−2 (x, k 2 ⊥ ) of the total TMD. Using the formalism of virtuality distribution functions (VDFs) [10,11] we find the explicit form of such a TMD [it coincides with the results of earlier studies [12,13] based on a particular on-mass-shell Ansatz for the parton-hadron blob χ(k, p)]. The form of F twist−2 (x, k 2 ⊥ ) is fully specified by the PDF f (x), and its k 2 ⊥ -support is limited by k 2 ⊥ ≤ x(1 − x)M 2 . As a consequence, the average transverse momentum induced by such a TMD is rather small. In particular, for a toy PDF f (x) = (1 − x) 3 , it is given by k 2 ⊥ = M 2 /30, which is about (170 MeV) 2 in case of the nucleon, that is considerably smaller than a folklore value of (300 MeV) 2 .
Our further study shows that the twist-2 part is not the only source of kinematic target-mass corrections: they also come from the higher-twist contributions. After incorporating the analysis of target-mass dependence for Feynman diagrams in the α-representation and studying equations of motion for the full TMD F (x, k 2 ⊥ ), we conclude that F (x, k 2 ⊥ ) should depend on k 2 ⊥ through the combination k 2 ⊥ + x 2 M 2 , and that this is the only "kinematically required" target-mass effect for the full TMD. Making this modification in the models used in Ref. [9], we observe that these M 2 /P 2 -corrections may be neglected well before the PQDs closely approach the limiting PDF form.
The paper is organized as follows. In Section 2, we start with recalling the definition of PQDs and their relation to VDFs established in Ref. [9]. Then, using the α-representation, we analyze the target-mass dependence of Feynman diagrams. In Section 3, we investigate the M 2 -dependence of the twist-2 part of the PQD Q(y, P). Using the VDF formalism, we also find the twist-2 parts of the relevant VDF and TMD. In section 4, we study the M 2 -dependence of higher-twist contributions and equations of motion for TMDs. Since the basic relations between various types of parton distributions are rather insensitive to complications brought in by spin, in Sections 2 -4 we refer to a simple scalar model. In Section 5, we discuss modifications related to quark spin and gauge nature of gluons in quantum chromodynamics (QCD). In Section 6, we discuss the k 2 ⊥ → k 2 ⊥ + x 2 M 2 modification of models for soft TMDs used in Ref. [9], and present numerical results for nonperturbative evolution of PQDs obtained in this modeling. Summary of the paper and our conclusions are given in Section 7.

Definition of PQDs and their relation to VDFs
The parton quasi-distributions originate from equal-time bilocal operator formed from two fields φ(0)φ(z) separated in space only [1], which corresponds to z = (0, 0, 0, z 3 ) [or, for brevity, z = z 3 ]. Then the PQDs are defined by (2.1) In our paper [9] we have analyzed the PQDs in the context of a general VDF representation [10,11] (where M 2 = p 2 ) that basically reflects the fact that the matrix element p|φ(0)φ(z)|p depends on z through (pz) and z 2 , and may be treated as a double Fourier representation with respect to these variables. The VDF representation holds for any p and z, but it is convenient to take the frame in which p = {E, 0 ⊥ , p 3 = P}. When z has only the minus component z − , the matrix element is parameterized by the parton distribution function (PDF) f (x) that depends on the fraction x of the target momentum component p + carried by the parton. The relation between the VDF Φ(x, σ) and the collinear twist-2 PDF f (x) is formally given by The σ-integral diverges when Φ(x, σ) has a ∼ 1/σ hard part generating perturbative evolution of PDFs. Our primary concern is nonperturbative evolution, so we will always imply the soft part of Φ(x, σ) for which the σ-integral converges. If we take z having just the third component, z = z 3 , we have This gives a relation between PQDs and VDFs, For large P, we have and Q(y, P → ∞) tends to the integral (2.4) producing f (y). The deviation of Q(y, P) from f (y) for large P may be described by higher-twist corrections in powers of Λ 2 /P 2 (where Λ is a scale like average primordial transverse momentum) and target mass corrections in powers of M 2 /P 2 .
As shown in our paper [9] , PQDs are completely determined by TMDs, so building models for TMDs we generate evolution patterns showing how Q(y, P) may depend on P due to the transverse-momentum effects.

Target mass dependence of VDFs
To discuss the origin of the target-mass dependence of VDFs it is convenient to switch to the momentum space description of the bilocal matrix element in terms of the function χ(k, p) (see Fig. 1) which is an analog of the Bethe-Salpeter amplitude [14]. A crucial observation is that the contribution of any (uncut) diagram to χ(k, p) may be written as (see, e.g., [15]) , where P(c.c.) is the relevant product of coupling constants, L is the number of loops of the diagram, and l is the number of its lines. The functions A(α), B s (α), B u (α), C(α), D(α) are sums of products of the nonnegative α j -parameters. Using Eq. (2.11) we get the representation with a function F(x, λ; p 2 ) specific for each diagram (2.11) Transforming Eq. (2.10) to the coordinate representation and changing λ = 1/σ gives (2.12) Note that the quadratic dependence on x in the exponential was produced by the k → z Fourier transformation: originally all terms in the exponential of Eq. (2.10) have linear dependence on x. Basically, one gets −x 2 M 2 after manipulating Absorbing the factor exp[−ix 2 M 2 /σ] into F(x, 1/σ) and defining the Virtuality Distribution Function gives the VDF representation (2.2).
Taking z that has z − and z ⊥ components only, i.e., projecting on the light front z + = 0, we define the Transverse Momentum Dependent Distribution in the usual way as a Fourier transform with respect to remaining coordinates z − and z ⊥ . The TMD may be written in terms of VDF as . In addition, Φ(x, σ; M 2 ), and hence also F (x, k 2 ⊥ ) have a "dynamical" or "kinematically unpredictable" M 2 -dependence contained in F(x, 1/σ; M 2 ) that comes from the last line in the α-representation (2.11).

PQD for twist-2 part
Another (and well-known) example of the kinematical target mass dependence is given by the M 2 -structure of the matrix elements of the twist-2 local operators.
To get the twist-2 part of the bilocal operator φ(0)φ(z), one should start with the Taylor expansion in z and then change the product of derivatives ∂ µ 1 . . . ∂ µ n into its traceless part {∂ µ 1 . . . ∂ µ n }. In a short-hand notation (z∂) n → {z∂} n , and (zp) n → {zp} n , so that Note that for z = z − we have {zp} n = (zp) n , which reproduces Eq. (2.3). To proceed in a situation with z 2 0, we use the fact that the structure of {zp} n is related to the Gegenbauer polynomials C 1 n (cosh θ) equal to Chebyshev polynomials U n (cosh θ) = sinh((n + 1)θ)/ sinh θ. As a result, where r = 1 − z 2 p 2 /(zp) 2 (see, e.g., Ref. [16]). Using p = (E, 0 ⊥ , P) and taking z = z 3 , we have (zp) = −z 3 P, z 2 = −z 2 3 and p 2 = M 2 . Thus, we have This gives and we get the twist-2 part of PQD in the form . This result was originally obtained (in somewhat different way and notations) in Ref. [8]. As noticed there, the integral over y is preserved One can check that the momentum sum rule also holds, Since the PQD Q(y, P) for negative y may come both from the y > 0 and y < 0 parts of the PDF f (y), it makes sense to split f (y) in these two parts and analyze PQDs coming from each of them separately. For illustration, we take the same model as in Ref. [9], namely, the function f ( According to Eq. (3.5), the PQD for positive y is obtained from the original f (y) by stretching it by factor (1 + ∆) in the horizontal direction and squeezing by factor (1 + 2∆) in the vertical one (see Fig. 2). For negative y, one should take f (−y) and contract it by factor ∆ in the horizontal direction, with the same squeeze by (1 + 2∆) in the vertical one.
Thus, if the twist-2 target mass corrections were the only ones here, it would be very easy to reconstruct such a PDF from the PQD at positive y: one should just perform the (1 + ∆) and (1 + 2∆) rescaling mentioned above.
For comparison, we show in Fig. 3 the P-dependence of PQD due to the nonperturbative evolution in the Gaussian model of Ref. [9]. Notice that the curve for P = 10 Λ is close in height to the P = M curve of Fig. 2. We expect that the scale Λ is about 300 to 500 MeV, or from 1/3 to 1/2 of the nucleon mass. Hence, 10Λ corresponds to about 3 -5 M. One can see that already the P = 2M curve from Fig. 2 is very close to the limiting curve (in this case ∆ = 0.06), which means that the target mass corrections in this comparison are visibly smaller than the higher-twist corrections governed by Λ (despite the fact that Λ was taken to be 2 -3 times smaller than M).

VDF for twist-2 part
The P-evolution patterns in Figs. 2 and 3 are rather different. It is interesting to find a physical reason for this difference. As shown in Ref. [9], the PQDs are completely determined by the TMDs, The Gaussian model mentioned above corresponds to a factorized Ansatz So, let us find out what kind of TMD corresponds to the twist-2 part of the matrix element. The first step is to find the VDF corresponding to the twist-2 contribution (3.1). To this end, we start with the decomposition of the traceless combinations over the usual ones, that follows from the ξ k expansion of the Gegenbauer polynomials C 1 n (ξ). This gives a double expansion in (pz) and z 2 for the sum in Eq.
At this stage, it is convenient to treat x > 0 and x < 0 parts of f (x) separately. For definiteness, we take x > 0. Notice now that where the functions f k (x) are defined by the recurrence relation with f 0 (x) = f (x). As a result, Comparing with the VDF representation (2.2), we find

TMD for twist-2 part
To proceed with the formula (2.14 ) producing the TMD we use which results in the δ (n) (k 2 ⊥ ) expansion that is equivalent to the following expression for the (k 2 ⊥ ) n mo- It is easy to check that the moment relation (3.20) is satisfied by the function In the M = 0 limit, we have Thus, no transverse momentum is generated in the case of a massless target. Our illustration model f ( wherex ≡ 1 − x. One can check that using the TMD (3.23) in the TMD→PQD conversion formula (3.8) one obtains the PQDs dictated by Eq. (3.5) and shown in Fig. 2.
The interpretation of the twist-2 approximation in terms of the transverse momentum dependent function given by Eq. (3.21) is known [12,13] from the early days of the ξ-scaling approach [3]. It was derived by imposing the k 2 = 0 condition on the hadron-parton blob χ(k, p) through the Ansatz while keeping the target mass finite p 2 = M 2 , see, e.g., Ref. [12]. In a similar context, Eq. (3.21) was obtained in Refs. [17,18] (see also [19]).
Our VDF-based derivation shows that the twist-2 TMD (3.21) can be obtained without additional assumptions.

Comparing TMDs
Note that, because the support of f (x) is 0 ≤ x ≤ 1, the twist-2 TMDs (3.21) vanish for k 2 ⊥ ≥ xxM 2 . This should be contrasted with the usual expectation (incorporated into our TMD models in Ref. [9] ) that TMDs are smooth functions of k ⊥ with a support extending to k 2 ⊥ = ∞. From a physical point of view, the twist-2 part F twist−2 (x, k 2 ⊥ ) describes a situation when a free massless quark happens somehow to be bound within a system with a total mass M. This results in a kinematic transverse momentum described by a rather artificially-looking TMD of Eq. (3.23) type. Clearly, this is just a model construction mimicking a hadron by a combination of free quarks with the total invariant mass M.
Comparing TMDs, it is instructive to calculate the average transverse momentum For comparison, the Gaussian TMD (3.9) gives k 2 ⊥ G = Λ 2 . Thus, taking Λ = M/3 we should expect that Λ 2 /P 2 corrections for PQD in the Gaussian model are about 3 times larger than the M 2 /P 2 corrections in the twist-2 part of the PQD. This observation explains the difference between Figs. 2 and 3.
Note that for more realistic valence PDFs f (x) that are singular for x = 0, the value of k 2 ⊥ twist−2 is even smaller. In particular, for f ( x it equals to M 2 /66, resulting in k 2 ⊥ twist−2 ≈ (116 MeV) 2 , i.e. factor of 8 smaller than the expected folklore value of 0.1 GeV 2 .
A rather exotic form of the twist-2 part of the TMD contradicts a natural expectation that TMDs should be smooth functions of k ⊥ with an unlimited support. To produce such a smooth TMD (having, moreover, a much larger k 2 ⊥ ), the higher-twist terms should literally wipe out the features brought in by the twist-2 term. This is only possible if the higher-twist terms also have the M 2 -dependence.

Twist decomposition
The twist-2 contribution appears as the first term in the twist decomposition of the original bilocal operator (see, e.g., Ref. [16]). The operators containing powers of ∂ 2 have higher twists, and their contribution to the light-cone expansion is accompanied by powers of z 2 . For PQDs, z 2 would result in a 1/P 2 suppression factor, just like for the target-mass corrections in twist-2 contribution. To analyze the interplay between the twist-2 and twist-4 terms, let us take the terms bilinear in z, As we discussed, {zp} 2 = (zp) 2 − z 2 M 2 /4 contains the z 2 M 2 target-mass correction term.
Since a VDF contains all information about the z-dependence of the original matrix element, it should provide the VDF representation for the twist-4 matrix element p|φ(0)∂ 2 φ(0)|p as well. To this end, we calculate z B(z, p) in the VDF representation (2.2) involving the x > 0 part of the VDF, and get (We remind that p in the VDF representation (2.2) is the actual hadron momentum, with p 2 = M 2 ). Assuming a soft Φ(x, σ) and taking z = 0, we get the twist-4 matrix element As one can see, it contains a term which a) is proportional to M 2 and b) is completely specified by the twist-2 PDF f soft (x). This means that the kinematical target-mass correction terms z 2 M 2 are contained not only in the twist-2 part of the original matrix element B(z, p), but also in its higher-twist parts. Most importantly, when substituted in Eq. (4.3), this term cancels the z 2 M 2 term coming from the twist-2 part. As a result, we have the expression free of the M 2 z 2 terms.

TMD parametrization
A similar result may be easily obtained in general case if one expands the exp[−iσz 2 /4] factor in the VDF representation (2.2) and uses the relation Then one obtains the representation of the matrix element in terms of the TMD F (x, k 2 ⊥ ). The sum over l gives the Bessel function J 0 , so we may also write The TMD parametrizations (4.8) and (4.9) provide another form of the z 2 -expansion, alternative to the twist decomposition (4.1). Its advantage is that the (pz)-dependence comes through the plane waves e −ix(pz) producing simple powers (pz) n rather than complicated traceless combinations {pz} n containing z 2 M 2 target-mass dependent terms that are simply artifacts of the twist decomposition. The TMD representation (4.9) is especially convenient in applications to PQDs. In particular, it directly leads to the TMD→PQD conversion formula (3.8).
One may argue that, due to equations of motion, like ∂ 2 φ = λψφ in a scalar λφ 2 ψ theory, one may write p|φ(0)∂ 2 φ(0)|p as p|φ (0) Still, it is an interesting question of how to incorporate equations of motions in the VDF/TMD parametrizations of the bilocal matrix element.

VDF parametrization for off-shell quarks
Since quarks in the nucleon are virtual, the matrix element B(z, p) does not satisfy the free-quark equation of motion z B(z, p) = 0. Keeping nonzero z and integrating by parts in Eq. (4.4), we obtain By equations of motion, this should be equal to the 3-body quark-quark-gluon contribution. For example, in a λφ 2 ψ scalar model, this should be equal to p|φ(0) λψ(z)φ(z)|p . Thus, building the VDF parametrization for the matrix element of the 3-body φψφ operator in a situation when ψ and one of the φ's are at the same point (and may be treated as one field) we should impose the condition reflecting equations of motion. For the TMDs constructed from Φ's using Eq. (2.14) (with k 2 ⊥ substituted by κ 2 to avoid too clumsy notations below) this gives 12) or, differentiating with respect to κ 2 , xF (x, κ 2 ) . (4.13) For the twist-2 part, when the l.h.s. of Eq. (4.13) vanishes, we have seen in Eq. (3.21) that the function xF (x, κ 2 ) depends on x and κ 2 through the combination η ≡ x + κ 2 /xM 2 . (4.14) Noticing that we can rewrite Eq. (4.13) in terms of x and η variables, Now, treating xF (x, κ 2 ) as a function G(x, η) of x and η, and introducing G 3 (x, η) ≡ ∂F φ(ψφ) (x, κ 2 )/∂κ 2 we have and finally If G 3 (x, η) vanishes, then we conclude that G(x, η) must be a function of η, in agreement with Eq. (3.21). If G 3 (x, η) does not vanish, the only restriction imposed by the equation of motion is Eq. (4.18). Thus, we may take any reasonable model for the two-body function G(x, η) and then just incorporate Eq. (4.18) [or original Eq. (4.11)] as a restriction that should be satisfied by the three-body function G 3 (x, η), when the qGq contribution is included, say, in a DIS calculation.
Of course, choosing a model for G(x, η) one should take care that the resulting G 3 (x, η) is also reasonable. In other words, if one has some information/expectations about the form of G 3 (x, η), one should make an effort to find a form of G(x, η) that would lead to the desired (or close) form of G 3 (x, η).
An important lesson is that, in the context of equations of motion, it is natural to build models of TMDs F (x, k 2 ⊥ ) in the form of functions of x and k 2 ⊥ + x 2 M 2 . This observation is in full accord with the general conclusion made at the end of the Section 2 that TMDs F (x, k 2 ⊥ ) must depend on k 2 ⊥ through the k 2 ⊥ + x 2 M 2 combination.

Equations of motion for spinor quarks
In spinor case, one deals with the matrix element of a B α (z, p) ≡ p|ψ(0)γ α ψ(z)|p (5.1) type. It may be decomposed into p α and z α parts: B α (z, p) = 2p α B p (z, p) + z α B z (z, p). These parts are not completely independent, since there are restrictions imposed by equations of motion.
Consider the handbag contribution for the virtual Compton amplitude, whose imaginary part gives the deep inelastic scattering (DIS) cross section. It may be written as where z β /z 4 comes from the spinor massless propagator , p), and s µναβ ≡ −g µν g αβ + g µα g νβ + g να g µβ .
To check the electromagnetic gauge invariance, we calculate The antisymmetric term will be eliminated if one takes B α (z, p) to be a derivative B α (z, p) = ∂ α B(z, p) of some "generating" scalar function B(z, p). After that, B α ,α = z B(z, p), and the equation of motion for B α ,α = p|ψ(0)/ ∂ψ(z)|p brings us to a study of z B(z, p), which completely parallels that performed in the previous section.
As for the remaining violation of the EM gauge invariance for the DIS handbag, it is proportional to zB (z, p), i.e., we still have it, as it is caused by the virtuality of the active quarks. In a Yukawa gluon model, we have / ∂ψ(z) = igφ(z)ψ(z), and this violation will be compensated when one includes terms coming from the 3-bodyψφψ diagrams, provided that one imposes the restriction (4.18).

Equations of motion in QCD
In QCD, one should take the operator O α q (z 2 , z 1 ; A) ≡ψ(z 2 ) γ αÊ (z 2 , z 1 ; A)ψ(z 1 ) (5.4) involving the gauge linkÊ(z 2 , z 1 ; A) along the straight line connecting z 1 and z 2 . The equation of motion, applied to the relative coordinate z, takes the form As a result, we have Taking again B α (z, p) = (∂/∂z α )B(z, p) reduces the equation of motion to the equation for z B(z, p) involving a scalar function B(z, p), and we can use all the results of Section 3, since the explicit form of Bψ Gψ (z, p) was not essential there.

M 2 -dependence of TMDs
Thus, if one uses the VDF/TMD representations (2.2), (4.9) for matrix elements, there are no kinematic z 2 M 2 -corrections that are artifacts of expansion over traceless {pz} n combinations. Furthermore, the PQDs are given by the conversion formula (3.8), and the target-mass dependence of Q(y, P) may only come from that of F (x, k 2 ⊥ ). According to the general statement made at the end of Section 2, the TMDs F (x, k 2 ⊥ ) must depend on k 2 ⊥ through the k 2 ⊥ + x 2 M 2 combination. This is a "predictable" or "kinematical" target-mass dependence.
We also noted there that F (x, k 2 ⊥ ) may have a "dynamical" M 2 -dependence due to the M 2 -dependence of the underlying function F(x, 1/σ ; M 2 ) of Eq. (2.12). This kind of M 2 -dependence cannot be derived from kinematics, and in this sense it is "unpredictable". In principle, there is nothing special in the fact that F(x, 1/σ ; M 2 ) depends on the hadron mass, just like there is no wonder that the shape of a PDF f (x) may be different if the hadron mass would be different. This is to say that some part of the M 2 -dependence of F(x, 1/σ ; M 2 ) may be absorbed into the form of the PDF f (x), and would not lead to M 2 /P 2 corrections describing the difference between a QPD Q(y, P) and its PDF f (y). Still, some part of the unpredictable M 2 -dependence may lead to the M 2 /P 2 corrections, and it is a challenge to build VDF models that would "realistically" reflect that part of the M 2 -dependence.
Leaving this problem for future studies, in what follows we will investigate the consequences of the "mandatory" change k 2 ⊥ → k 2 ⊥ + x 2 M 2 in the TMD models that have been used for generating nonperturbative evolution of PQDs in our paper [9].

Gaussian model
Adding the M 2 -dependence into our Gaussian model (3.9) by the k 2 ⊥ → k 2 ⊥ + x 2 M 2 prescription, we get wheref (x) = f (x) e −x 2 M 2 /Λ 2 . Thus, we have a simple change in the form of the PDF, f (x) →f (x), that would not be reflected by M 2 /P 2 terms in the difference betweenQ(x, P) andf (x).

Simple non-Gaussian model
Another VDF model proposed in Ref. [9], intends to reproduce the large-|z| exponential ∼ e −|z|m fall-off of the perturbative propagator D c (z.m) of a particle with mass m, while removing its 1/z 2 singularity at small z 2 by a "confinement" factor e iσ/Λ 2 reflecting the finite size of a hadron. This model corresponds to the TMD given by Using the prescription k 2 ⊥ → k 2 ⊥ + x 2 M 2 amounts to the change m 2 → m 2 + x 2 M 2 in this model. To avoid a two-parameter (Λ and m) modeling, in our paper [9] we took m = 0. Let us do the same here. In the context of the m-model (6.2), the resulting TMD model (6.4) corresponds to assuming that the parton mass m is a fraction xM of the nucleon mass M. This assumption does not look absolutely unnatural in view of the fact that the VDF representation (2.2) involves the plane wave factor e −ix(pz) in which p is the actual hadron momentum p satisfying p 2 = M 2 .
For the quasi-distribution, the model (6.4) gives

Numerical results
Now we have two parameters, the nucleon mass M and the transverse momentum scale Λ, and we need to decide what is their ratio. To this end, we calculate the average transverse momentum in the model of Eq. (6.4) with f (x) = (1 − x) 3 , and find To illustrate the impact of the M 2 terms in Eq. (6.5) on the shape of quasi-distributions, we take again f (y) = (1 − y) 3 θ(y), and compare curves for M/Λ = 3 and M/Λ = 0 at P/Λ = 3 (i.e. P/M = 1). As one can see from Fig. 4, the two curves are very close to each other. At the same time, they are still very far from the limiting (1 − y) 3 shape. Increasing P to 2M, we get the curves that practically coincide (see Fig. 5), still being rather far from the asymptotic P/Λ → ∞ shape.
Thus, in this scenario, when one reaches the momentum P that is sufficiently large to stop the nonperturbative evolution of the PQD Q(y, P), there is no need to bother about target mass corrections. Given the expected accuracy of lattice gauge calculations, they may be safely neglected starting with P ∼ 2M.

Summary and conclusions
In this paper, we have studied the target-mass dependence of the parton virtuality distributions. Our main result is that if one uses the VDF/TMD representations (2.2), (4.9) for matrix elements, there are no kinematic z 2 M 2 -corrections that are an inherent feature/artifact of expansions over traceless {pz} n combinations that appear in the twist decomposition. In our approach, the PQDs are given by the TMD→PQD conversion formula (3.8). In the P → ∞ limit of the latter, the PQD Q(y, P) tends to the twist-2 PDF f (y) irrespectively of the fact that the VDF/TMD representation does not involve the twist decomposition.
We have established that TMDs F (x, k 2 ⊥ ; M 2 ) must depend on k 2 ⊥ through the k 2 ⊥ + x 2 M 2 combination. Hence, the x 2 M 2 addition here may be considered as a kinematic targetmass correction. Furthermore, TMDs may have a dynamic M 2 -dependence that cannot be predicted from kinematical considerations. Just like the form of the k 2 ⊥ -dependence of the TMDs, this part of the M 2 -dependence can only be modeled in our approach.
We have studied the effect of the k 2 ⊥ → k 2 ⊥ + x 2 M 2 modification of the TMD models used in our paper [9], and found that the M 2 /P 2 corrections become negligible well before the PQD curves Q(y, P) become close enough to the corresponding PDF f (y). Thus, we see no need to correct the lattice gauge calculations of PQDs for M 2 -effects.
A similar analysis of the target-mass effects can be made for the pion quasi-distribution amplitude studied recently on the lattice in Ref. [7] and in the VDF approach in Ref. [20]. Since the pion mass m π is much smaller than the nucleon mass M (even when m π is taken in its lattice version m π ∼ 310 MeV [7]), while the pion size scale ∼ 1/Λ is not very different from that of the nucleon, the target-mass effects in that case may be completely ignored.
A possible future extension of our findings is an application of the VDF/TMD approach to inclusive DIS, with the goal to investigate if the target-mass corrections described there by the Nachtmann ξ variable [2] are a genuine feature of the process or just an artifact of the twist decomposition.