One-loop evolution of parton pseudo-distribution functions on the lattice

We incorporate recent calculations of one-loop corrections for the reduced Ioffe-time pseudo-distribution ${\mathfrak M} (\nu,z_3^2)$ to extend the leading-logarithm analysis of lattice data obtained by Orginos et al. We observe that the one-loop corrections contain a large term reflecting the fact that effective distances involved in the most important diagrams are much smaller than the nominal distance $z_3$. The large correction in this case may be absorbed into the evolution term, and the perturbative expansion used for extraction of parton densities at the $\mu \approx 2$ GeV scale is under control. The extracted parton distribution is rather close to global fits in the $x>0.1$ region, but deviates from them for $x<0.1$.


INTRODUCTION
Feynman's parton distribution functions (PDFs) [1] f (x) are the crucial building blocks in the description of hard inclusive processes in quantum chromodynamics (QCD). Accumulating nonperturbative information about the hadron structure, the PDFs are a natural subject for a lattice study. However, straightforward definitions of PDFs refer to matrix elements of bilocal operators on the light cone z 2 = 0, the intervals inaccessible on the Euclidean lattice.
The ideas of how to get information from space-like intervals date to the pioneering paper of W. Detmold and D. Lin [2] who proposed a lattice study of the deepinelastic-type Euclidean correlators of heavy-light currents. Later, V. Braun and D. Müller [3] proposed to use Euclidean correlators to extract the pion distribution amplitude, another function [4] playing a fundamental role in perturbative QCD studies of hard exclusive processes. The use of correlators in the form of "lattice cross sections" was more recently advocated in the papers by Qiu and collaborators [5,6].
The current correlators involve a quark propagator connecting the current vertices. This factor is avoided in the proposal by X. Ji [7] to study the quasi-PDFs Q(y, p 3 ) that describe the distribution of the spatial z 3 -component of the hadron momentum p 3 . While being different from the Feynman PDFs f (y) describing the distribution of the hadron's "plus"-momentum p + = p 0 + p 3 , they coincide with f (y) in the infinite momentum limit p 3 → ∞.
Both on the lattice and in the usual continuum space, the basic object for all types of PDFs is the matrix element M (z, p) generically (i.e. ignoring the inessential spin complications) written as p|φ(0)φ(z)|p . By Lorentz invariance, it is a function of the Ioffe time (pz) ≡ −ν [8] and the interval z 2 , M (z, p) ≡ M(ν, −z 2 ).
In the (formal) light-cone limit z 2 = 0, the Fourier transform of M(ν, 0) with respect to ν gives f (x). In this sense, the ν-dependence of the Ioffe-time distribution (ITD) M(ν, −z 2 ) reflects the longitudinal structure of the PDFs. As shown in Ref. [9], the z 2 -dependence of M(ν, −z 2 ) determines the k ⊥ -dependence of the (straight-link in the case of QCD) transverse momentum dependent parton distributions (TMDs) F(x, k ⊥ ).
Since the quasi-PDFs Q(y, p 3 ) are given by the Fourier transform of M(z 3 p 3 , z 2 3 ) with respect to z 3 , their shape is distorted by nonperturbative transverse momentum effects entering through the second argument of M(z 3 p 3 , z 2 3 ). While, in a general perspective, the k ⊥ -dependence of F(x, k ⊥ ) provides information about the three-dimensional structure of hadrons, in the case of the quasi-PDFs it is a nuisance responsible for the unwanted difference between Q(y, p 3 ) and f (y) that is very strong at momenta reached in existing lattice calculations of quasi-PDFs.
To decrease the impact of the z 2 -dependence of the ITD M(ν, −z 2 ), it was proposed [10] to consider the reduced ITD M(ν, −z 2 ) given by the ratio of M(ν, −z 2 ) and the rest-frame distribution M(0, −z 2 ). Though there are no first-principle grounds that the nonperturbative part of the z 2 -dependence disappears in this ratio, it is natural to expect that it is strongly reduced.
The ideal case when M(ν, −z 2 ) is just a function of ν corresponds to factorization of the x and k ⊥ dependencies of the TMD F(x, k ⊥ ). In fact, the idea that 1 GeV 2 is a standard assumption of the TMD practitioners (see, e.g., Ref. [11]), with a Gaussian being the most popular form for K(k ⊥ ).
An exploratory lattice study of the reduced ITD was performed in Ref. [12] (and also described in Ref. [13]). The results show that M(ν, z 2 3 ) is basically a universal function of ν, with small deviations from the common curve for the points corresponding to the smallest values of z 3 .
As demonstrated in Ref. [12], these deviations may be explained by perturbative evolution. While the leading logarithmic approximation (LLA) used in Ref. [12] is sufficient to analyze the ln z 2 3 dependence, one needs to go beyond it to specify the scale µ which should be attributed to the extracted scale-dependent PDFs f (x, µ 2 ).
To this end, one needs complete expressions for one-  [14,15]. Our goal here is to give a more detailed discussion of the LLA treatment of the evolution, and also to extend the analysis beyond the LLA. As we will show, the one-loop correction contains a large contribution that considerably changes the results obtained in the LLA.
To make this article self-contained, we outline in Sec. II the basics of the Ioffe-time distributions and pseudo-PDFs. In Sec. III, we discuss the structure of one-loop corrections. In Sec. IV, we describe the evolution effects revealed in the lattice study of Ref. [12], and convert the data for the reduced ITD M(ν, z 2 3 ) into the standard parton densities f (x, µ 2 ) defined in the MS scheme. The summary of the paper and conclusions are given in Sec. V.

II. IOFFE-TIME DISTRIBUTIONS AND PSEUDO-PDFS
The basic object for defining parton distributions is a matrix element of a bilocal operator that (skipping inessential details of its spin structure) may be written generically like p|φ(0)φ(z)|p . Due to invariance under Lorentz transformations, it is given by a function of two scalars, the Ioffe time (pz) [8] (which will be denoted by −ν) and the interval z 2 (again, the sign for the second argument is chosen so as to have a positive value for spacelike z). One can demonstrate [9,16] that, for all relevant Feynman diagrams, its Fourier transform P(x, −z 2 ) with respect to (pz) has −1 ≤ x ≤ 1 as support, i.e., In this covariant definition of x, one does not need to assume that z is on the light cone z 2 = 0 or that p is light-like p 2 = 0. On the light cone z 2 = 0, we formally have P(x, 0) = f (x). Hence, the function P(x, −z 2 ) may be treated as a generalization of the concept of PDFs onto non-lightlike intervals z 2 , and following [10] , we will refer to it as the pseudo-PDF. In view of lattice applications, we will take the separation z = {0, 0, 0, z 3 } oriented in the direction specified by the hadron momentum p = {E, 0, 0, P }.
In renormalizable theories (including QCD), the function M(ν, −z 2 ) has logarithmic ∼ ln(−z 2 ) singularities. In deep inelastic scattering (DIS), they result in a logarithmic scaling violation with respect to the photon virtuality Q 2 . A wide-spread statement is that the Q 2 -dependent DIS structure functions W (x B , Q 2 ) probe the hadron structure at distances ∼ 1/Q. In the case of the pseudo-PDFs P(x, z 2 3 ), one may say that they literally describe the hadron structure at the distance z 3 .
Just like the DIS form factors W (x B , Q 2 ) are written in terms of the universal parton densities f (x, Q 2 ), the pseudo-PDFs obtained from lattice calculations may be expressed through the usual parton distributions. The latter are defined by the operators on the light cone z 2 = 0, i.e., in a logarithmically singular limit. In the approach based on the operator product expansion (OPE), the standard procedure is to remove these singularities with the help of some prescription.
The most popular of them is the MS scheme based on the dimensional regularization. Consequently, the resulting PDFs have a dependence on the renormalization scale µ, and therefore one should write the PDFs as f (x, µ 2 ).
Switching from x to the Ioffe time ν gives the functions introduced in Ref. [17] and called there the Ioffe-time distributions. In this context, the functions M(ν, −z 2 ) that are the Fourier transforms of pseudo-PDFs, should be called the Ioffe-time pseudo-distributions or pseudo-ITDs.
To get a relation between the pseudo-PDFs P(x, z 2 3 ) and the MS parton densities f (x, µ 2 ), one can use the nonlocal light-cone OPE [18,19] (see also [15]) for the matrix element defining P(x, z 2 3 ), i.e., for the pseudo-ITD. The result has the structure similar to that of the usual OPE for the DIS structure functions W (x, Q 2 ). In this expression, the twist-2 coefficient functions C i are given by an expansion in the strong coupling constant α s , while O(z 2 ) symbolizes higher-twist terms. However, the application of the OPE to the pseudo-ITDs and pseudo-PDFs in QCD faces complications related to the gauge link. Namely, when z is off the light cone, the link generates linear ∼ z 3 /a and logarithmic ∼ ln(1 + z 2 3 /a 2 ) ultraviolet (UV) divergences, where a is an UV regulator with the dimension of length (it may be a finite lattice spacing). Though disappearing for z 3 = 0, these divergences require an additional UV regularization when z 3 is finite.
As stated in Ref. [12], for small spacelike intervals z 2 = −z 2 3 , and at the leading logarithm level, the reduced pseudo-PDFs are related to the MS distributions by a simple rescaling of their second arguments, namely, where γ E is the Euler's constant (a more detailed discussion will be given later on). This rescaling factor is very close to 1, since 2e −γ E = 1.12. However, this factor may be changed by the O(α s ) terms present in the coefficient function.

A. Confinement and infrared cut-offs
There are several standard techniques to calculate gluon radiative corrections in QCD. Most of them are oriented to work in the region of absolute perturbative QCD (pQCD) applicability. A straightforward use of such methods, however, may need some care in applications involving energy scales that are not very large. For this reason, let us discuss some features of calculations on the border of applicability of perturbative methods.
To begin with, one should remember that quarks and gluons are confined, i.e. the propagators of all diagrams (even in a continuum case) are embedded in a finite volume whose size is determined by the hadron's radius. The confinement effects lead, in particular, to a rapid decrease of correlators like ITDs or pseudo-PDFs at distances z 3 larger than the hadronic radius R. Still, at short distances one can use asymptotic freedom and obtain, in particular, the ln z 2 3 singularities. Thus, it makes sense to treat pseudo-ITDs and pseudo-PDFs as sums of the soft and hard parts. The soft part basically reflects the size of the system and is assumed to be finite for z 3 = 0. The hard part is singular for z 3 → 0, and is produced by perturbative interactions. The hard part may be visualized then as generated from the soft part through a hard exchange kernel H(0, z; z 1 , z 2 ), In the standard pQCD factorization approaches, the soft part is mimicked by on-shell parton states, and the ln z 2 3 -singularities appear either as ln(z 2 3 m 2 ), where m is the parton mass or ln(z 2 3 µ 2 IR ), where µ IR is the scale used in dimensional regularization of infrared singularities in the case of massless partons.
Since M soft (z 1 , z 2 ) in Eq. (7) rapidly decreases for large separations |z 1 −z 2 |, the hadronic size R provides an infrared cut-off for the integral, even when the quarks are massless. While at short distances one gets the ln(z 2 3 /R 2 ) behavior, the logarithmic form is just an approximation valid for z 3 R. Such a restriction may be hard to implement on the lattice.
Of course, the exact form of the IR regularization imposed by confinement is not known. To get a feeling, let us take an infrared regularization by a mass term. A typical integral producing the ln z 2 3 singularity then has the form where α is the Schwinger's α-parameter and m is the infrared regulator. One can see that where K 0 (mz 3 ) is the modified Bessel function. Its expansion for small z 3 explicitly shows the expected ln(z 2 3 m 2 ) singularity. The usual pQCD factorization procedure is to split ln(z 3 /R) into the short-distance part ln(z 3 µ) that is attributed to the coefficient function and the long-distance part ln(1/µR) that is absorbed into the "renormalized" PDF f (x, µ 2 ). Given the commonly used lattice spacing a ∼ 0.1 fm and the hadron size R 1 fm, the question is whether there is enough interval for the logarithmic part of the z 3 -dependence to be visible in the data at all.
An important feature of the Bessel function K 0 (mz 3 ) is that it exponentially decreases when z 3 exceeds the infrared cut-off 1/m. Thus, if instead of the short-distance approximation of I K (z 2 3 ) by ln(1/z 2 3 ), one would use the "exact" I K (z 2 3 ) function for the evolution term, there will be no evolution corrections for large z 3 . In other words, the logarithmic evolution disappears at large distances.

B. Rescaling relation
To fix a relation between the pseudo-PDF scale z 3 and the MS scale µ, one should take into account constant terms, like 2 ln(2e −γ E ) in Eq. (9). In the MS-OPE approach, one takes z 2 = 0 and then applies the dimensional regularization which adds the α factor into the integral (8) making it convergent. After that, one uses the MS-prescription, which is arranged to produce exactly ln(µ 2 /m 2 ) as the result in this case, Thus, the constant term in Eq. (9) provides the leadinglogarithm rescaling coefficient 2e −γ E between the pseudo-PDFs and MS parton distributions expressed by Eq. (6).
One may ask what happens if one uses another type of the IR regularization. In particular, the Gaussian models for TMDs suggest that the decrease for large z 3 is also Gaussian. One may expect that the hard correction should resemble, for large z 3 , the behavior of the soft part. Thus, the exponential e −m|z3| fall-off of the modified Bessel function may look too slow. A Gaussian decrease can be easily provided by a sharp IR cut-off applied to Eq. (8). For small z 2 3 , the incomplete gammafunction Γ(0, z 2 3 /z 2 3 ) has a logarithmic singularity while for large z 2 3 , the function I G (z 2 3 ) has a Gaussian e −z 2 3 /z 2 0 fall-off. Again, we can calculate the z 3 = 0 version of Eq. (11) using the MS-scheme to obtain One can see that the pseudo-PDF/PDF rescaling (6) remains intact. This is a natural result, because the relation between the finite-z 3 and MS cut-offs concerns only the short-distance properties of the bilocal operator.

C. One-loop correction
The discussion given in the previous section addresses only the overall rescaling between two regularization schemes (just like the relation between the values of the QCD scale Λ in, say, MOM and MS schemes). To establish a connection between the pseudo-PDFs and the MS-PDFs, we need, in addition, the constant part of the oneloop coefficient function in the nonlocal OPE of Eq. (4). It was given in Refs. [14] and [15], with some differences between them. After rechecking our calculation and fixing typos, we present our result in the form Turning to the PDF counterpart, we take z 2 = 0 and using the MS scheme for the UV divergence, obtain The logarithmic part here involves a convolution that may be symbolically written as B ⊗ M (ν) where is the Altarelli-Parisi kernel [28]. Combining Eqs. (14) and (15), we obtain the relation which is in agreement with a recent result of Ref. [29] (see also Ref. [30]). Eq. (17) allows one to convert the data points for M(ν, z 2 3 ) into the "data" for I(ν, µ 2 ). The first contribution in the second line is an obvious term reflecting the general multiplicative scale difference between the z 2 and MS cut-offs. If all the further terms are neglected, then the only difference between M(ν, z 2 3 ) and I(ν, µ 2 ) is just the rescaling µ 2 = 4e −2γ E /z 2 3 . In that case, one can evolve the M(ν, z 2 3 ) data to a particular z 3 value z 0 , and treat (in this approximation) the resulting function M(ν, z 2 0 ) as the MS ITD corresponding to the scale µ = 2e −γ E /z 0 , which is numerically close to 1/z 0 . This simple rescaling relation (used in Ref. [12]) is modified when the further terms of Eq. (17) are included. In particular, the term proportional to the Altarelli-Parisi kernel B(w) may be absorbed into the ln z 2 3 term, which would just change the rescaling relation into µ = 2e −1/2−γ E /z 0 .
The term with [ln(1 − w)]/(1 − w) produces a large negative contribution. In Feynman gauge, according to Ref. [15], it comes from the evolution part of the vertex diagrams involving the gauge link (see Fig.1). The key point is that the gluon is attached there to a running tz 3 position on the link. After integration over t, etc., the net outcome is that the z 3 -dependence of these diagrams is generated by an effective scale smaller than z 3 . Indeed, let us combine the [ln(1 − w)]/(1 − w) term with the ln z 2 logarithm by rewriting Eq. (17) as We see that z 3 enters now through a running (1 − w)z 3 location. The remaining (w + 1) ln(1 − w) term is much less singular than B(w) for w = 1, and does not produce large contributions.
Thus, the magnitude of the one-loop correction is governed by the combined evolution logarithm. It cannot be made zero by a particular choice of µ because it depends on the integration variable w. Still, the w-integrated contribution will vanish for some µ that we may write as where 1 − w is the "average" value of 1 − w. Since B(w) is strongly enhanced for w = 1, we should expect that 1 − w is numerically small, leading to a µ ∼ k/z 3 rescaling with a rather large coefficient k. As we will see, k ∼ 4 in this case.
Again, one may ask if the perturbative formula (17) involving the ln z 2 3 logarithm may be applied to actual lattice data. In particular, our exercise with the massterm IR regularization and the resulting Bessel function shows that the logarithmic behavior ln z 2 3 of the hard term is valid only for z 3 values well below the IR cut-off R, which is given by the hadron size in our case. Hence, a practical question is whether the data really show a logarithmic evolution behavior in some region of small z 3 .

A. General features
An exploratory lattice study of the reduced pseudo-ITD M(ν, z 2 3 ) for the valence u v − d v parton distribution in the nucleon has been reported in Ref. [12]. An amazing observation made there was that, when plotted as functions of ν, the data both for real and imaginary parts lie close to respective universal curves. The data show no polynomial z 3 -dependence for large z 3 . Given that z 2 3 /a 2 changes in the explored range from 1 to about 200, we interpret this result as the total absence of higher-twist terms in the reduced pseudo-ITD.
As explained in Refs. [10,12] and in the Introduction, such an outcome corresponds to a factorization of the νand z 2 3 -dependences of the soft part of the Ioffe-time distribution M(ν, z 2 3 ) = M (ν)M(0, z 2 3 ). In terms of TMD F(x, k 2 ⊥ ), this corresponds to factorization of its xand k 2 ⊥ -dependences in the region of soft k ⊥ . However, as observed in Ref. [12], there is quite visible z 3 -dependence for small values of z 3 , namely, z 3 6a, that may be explained by perturbative evolution.
Let us consider first the real part. It corresponds to the cosine Fourier transform of the function q v (x) corresponding to the valence combination, i.e., the difference q v (x) = q(x) −q(x) of quark and antiquark distributions. In our case, q = u − d. In Ref. [12], it was found that the data for the real part are very close (see Fig. 2) to the curve R f (ν) generated by the function  This shape was obtained by forming cosine Fourier transforms of the normalized x a (1 − x) b -type functions and fixing the parameters a, b through fitting the data. While all the data points have been used in the fit, the shape of the curve is obviously dominated by the points with smaller values of Re M(ν, z 2 3 ). To give a more detailed illustration, we show in Fig. 3 the points corresponding to z 3 values in the range 7a ≤ z 3 ≤ 13a. As one can see, there is some scatter for the points with the largest values of ν in the region ν 10, where the finitevolume effects become important. Otherwise, practically all the points lie on the universal curve based on f (x). In this sense, there is no z 3 -evolution visible in the large-z 3 data.
In Fig. 4, we show the points in the region a ≤ z 3 ≤ 6a (note that, on the lattice, z 3 = 0 means that also ν = 0, and M(0, 0) = 1 by definition). In this case, all the points lie higher than the universal curve. We recall that the perturbative evolution increases the real part of the pseudo-ITD when z 3 decreases. Thus, one may conjec- ture that the observed higher values of R for smaller-z 3 points may be a consequence of the evolution. A typical pattern of the z 3 -dependence of the lattice points is shown in Fig. 5 for a "magic" Ioffe-time value ν = 3π/4 that may be obtained from five different combinations of z 3 and P values used in Ref. [12]. The shape of the eye-ball fit line is given by the incomplete gammafunction Γ(0, z 2 3 /30a 2 ). This function entirely conforms to the expectation that the z 3 -dependence has a "perturbative" logarithmic ln(1/z 2 3 ) behaviour for small z 3 , and rapidly vanishes for z 3 larger than 6a.
As expected, R(ν, z 2 3 ) decreases when z 3 increases. We also see that the evolution "stops" for large z 3 . In this context, the overall curve based on Eq. (21) corresponds to the "low normalization point", i.e., to the region, where the perturbative evolution is absent.

B. Building MS ITD
Thus, we see that the data of Fig. 5 show a logarithmic evolution behavior in the small z 3 region. Still, the z 3 -behavior starts to visibly deviate from a pure logarithmic ln z 2 3 pattern for z 3 5a. This sets the boundary z 3 ≤ 4a on the "logarithmic region". So, let us try to use Eq. (17) in that region to construct the MS ITD.
It is instructive to split the contributions in Eq. (17), where we will denote Re I(ν, µ 2 ) ≡ I R (ν, µ 2 ). The first, "evolution" part, given by (recall that R(ν, z 2 3 ) ≡ Re M(ν, z 2 3 )) corresponds to the leading logarithm approximation used in Ref. [12]. For z 3 = 2e −γ E /µ, the logarithm vanishes, and we have This happens, of course, only if, for an appropriately chosen α s , the ln z 2 3 -dependence of the one-loop correction cancels the actual z 2 3 -dependence of the data, visible as scatter in the data points in Fig. 4. In Ref. [12], it was found that this happens when α s /π ≈ 0.1. Thus, Eq. (17) is accurate only in the region, where the data show a logarithmic dependence on z 3 , i.e., z 3 ≤ 4a in our case.
As we have discussed, the L ⊗ R f term reflects the fact that the actual scale in the evolution part of the vertex diagrams is less than z 3 . To illustrate its impact, we show, in Fig. 6, the functions B ⊗ R f and L ⊗ R f . One can see that the last one is negative and rather large. Its ν-dependence is similar to that of the B⊗R f function. In fact, in the ν < 5 region, we have L ⊗ R f ≈ −3.5B ⊗ R v . Thus, the combined effect of these two terms is close to that of −2.5B ⊗ R f . As a result, the inclusion of these terms may be approximately treated as a LLA evolution ⌫ I R (⌫) FIG. 7. Function IR(ν, µ 2 ) for µ = 1/a calculated using the data with z3 from a to 4a. The upper curve corresponds to the ITD of the CJ15 global fit PDF.
with a modified rescaling factor. Specifically, we may write Thus, the rescaling factor has changed by a factor of 4 compared to the original LLA value! We may use µ ≈ 4/z 3 as a guide, but the actual numerical calculations should, of course, be done using the "exact" Eq. (17). To proceed, we choose the value µ = 1/a which, at the lattice spacing of 0.093 fm used in Ref. [12] is approximately 2.15 GeV. The estimate (25) tells us that the ITD I R (ν, µ 2 ) at this scale should be close to the pseudo-ITD R(ν, z 2 3 ) for z 3 ≈ 4a, a distance that is on the border of the z 3 ≤ 4a region. Taking the value α s /π = 0.1 used in Ref. [12] and applying the full oneloop relation (17) to the data with z 3 ≤ 4a, we generate the points for I R (ν, (1/a) 2 ).
As seen from Fig. 7, all the points are close to some universal curve with a rather small scatter. The curve itself was obtained by fitting the points by the cosine transform of a normalized N x a (1−x) b distribution, which gave a = 0.35 and b = 3. The magnitude of the scatter illustrates the error of the fit for the ITD in the ν ≤ 4 region. In Fig. 7, we compare our µ = 1/a ITD with the ITD obtained from the global fit PDFs corresponding to the CJ15 [31] global fit. One can see that our ITD is systematically below the curve based on the global fit PDFs.
The reason for the discrepancy may be understood from Fig. 8, where we compare the normalized N x 0.35 (1 − x) 3 ≡ q v (x, µ = 2.15 GeV) distribution to CJ15 [31] and MMHT 2014 [32] global fit PDFs, taken at the scale µ = 2.15 GeV. Unlike the ∼ x 0.35 function, these PDFs are singular for small x, which leads to the enhancement of ITDs for large and moderate values of ν.
To fit the points for I R (ν, µ 2 ), we have used the same simplest N x a (1 − x) b Ansatz for the PDF as in Ref. [12].
15 GeV built from the data shown in Fig. 7 and compared to CJ15 and MMHT global fits.
In principle, one may use more complicated models for PDFs and get practically the same fitted curve for the ITD in the ν ≤ 4 region, while a somewhat different curve for PDF q v (x). The reason is simple: the inverse cosine Fourier transform is unique only when one exactly knows the ITD in the whole 0 ≤ ν < ∞ region. Performing such a transform from a limited ν ≤ 4 region, one needs to add some assumptions either about the behavior of the ITD outside this region or about a functional form of the PDF q v (x). We fixed our choice by taking The study of how the shape of q v (x) varies if one uses more complicated forms, in particular, those used in the global fits [31,32] is an interesting problem that, however, goes beyond the scope of the present paper.
Comparing to the LLA results of Ref. [12], we observe that the large negative one-loop correction in Eq. (17) has visibly changed the extracted PDF, which is now further from the global fit PDFs. The main reason is that the z 0 = 2a pseudo-ITD constructed in Ref. [12] was treated there as corresponding to the µ ≈ 1 GeV scale, while according to the modified rescaling relation (25), it should correspond to µ ≈ 4 GeV. Hence, to get the µ ≈ 2 GeV curve, one needs to evolve it down in µ.
Still, the guiding idea of Ref. [12], that the MS ITDs I R (ν, µ 2 ) can be obtained from the reduced pseudo-ITDs R(ν, z 2 3 ) by an appropriate rescaling µ = k/z 3 , works with a rather good accuracy for all z 3 ≤ 6a if one takes k ≈ 4. By this rescaling relation, the µ = 1/2a ≈ 1 GeV ITD corresponds to the z 3 ≈ 8a reduced pseudo-ITD. As we discussed, a boundary point beyond which the evolution stops, is z 3 ≈ 6a. Hence, the pseudo-ITD at this distance is given by the ITD R f (ν) corresponding to the universal fit function f (x) of Eq. (21). This result may be also obtained by a direct numerical calculation based on Eq. (17).
Using Eq. (17) one may also evolve the MS ITD below µ = 1/2a, and the resulting functions will be changing with µ. On the other hand, the pseudo-ITDs do not change with z 3 when z 3 6a. Hence, the rescaling connection I R (ν, µ 2 ) ≈ R(ν, (4/µ) 2 ) in this region becomes less and less accurate when µ decreases, and eventually makes no sense.

C. Imaginary part
Imaginary part of the pseudo-ITD may be considered in a similar way. It corresponds to the sine Fourier transform of the function given by the sum q(x) +q(x) of quark and antiquark distributions. This function differs from the valence combination Fig. 9, we show the data for large z 3 values z 3 ≥ 7a. Just like in the case of the real part (see Fig. 3), the points with ν 10 are close to a universal curve. Representing q(x) +q(x) = q v (x) + 2q(x) and taking f (x) of Eq. (21) as q v (x), we find Note that in Ref. [12], the fit was made for all the z 3 points (i.e. the points with z 3 ≤ 6a have been also included), and the overall coefficient forq(x) was obtained to be 0.07 rather than 0.1. In Fig. 10, we show data with z 3 ≤ 4a. As one can see, all these points are below the curve obtained by fitting the z 3 ≥ 7a data. This is in agreement with the fact that, in the region ν 6, the perturbative evolution decreases . Function II (ν, µ 2 ) for µ = 1/a calculated using the data with z3 from a to 4a. The curve is described in the text. the imaginary part of the pseudo-ITD when z 3 decreases. Note that the 1-loop relation holds for the whole function M = Re M + i Im M. So, we should just separate there real and imaginary parts, and the construction of the MS function Im I(ν, µ 2 ) ≡ I I (ν, µ 2 ) proceeds in the same way as for the real part.
The results are shown in Fig. 11. Again, all the points are rather close to a universal curve with a rather small scatter. The curve shown corresponds to the sine Fourier transform of the sum of the valence distribution q v (x, µ = 1/a) = N x 0.35 (1 − x) 3 obtained from the study of the real part, and the antiquark contribution 2q(x, µ = 1/a). The latter was found from the fit to be given byq(x, µ = 1/a = 2.15 GeV) = 0.07[20x(1 − x) 3 ].

V. SUMMARY AND CONCLUSIONS
In this paper, we have extended the leading-logarithm analysis of lattice data for parton pseudo-distributions and reduced pseudo-ITDs performed in Ref. [12]. To this end, we incorporated recent results for the reduced pseudo-ITDs at the one-loop level [15] (see also [14,29]).
It was found that the correction contains a large term resulting in essential numerical changes compared to the LLA. The large correction appears since effective distances involved in the most important diagrams are much smaller than the nominal distance z 3 . This leads to a change (from k LLA ≈ 1 to k ≈ 4 in the case of our particular ITDs) of the coefficient k in the rescaling relation µ = k/z 3 that allows to (approximately) convert the pseudo-PDFs P(x, z 2 3 ) into the MS PDFs f (x, µ 2 ). While the rescaling relation serves as an instructive guide for quick estimates and semi-quantitative analysis, the MS ITDs may be directly constructed applying the exact one-loop formula. Using it, we have obtained the ITD I(ν, µ 2 ) at the µ = 1/a ≈ 2.15 GeV MS scale using the data in the 0 ≤ z 3 ≤ 4a region.
We found that I(ν, µ 2 ) at this scale is close to the reduced pseudo-ITD M(ν, z 2 3 ) for z 3 ∼ 4a. Since all the data in the a ≤ z 3 ≤ 4a region do not differ much from the z 3 = 4a ones (see Fig. 4), the conversion of the M(ν, z 2 3 ) data into I(ν, 1/a 2 ) does not involve large changes, i.e., the perturbative expansion for MS ITD I(ν, µ 2 ) in terms of the reduced pseudo-ITDs M(ν, z 2 3 ) is under control. A formal reason is that the large correction in this case can be absorbed into the z 2 3 -dependent evolution term, with remaining corrections being small.
Phenomenologically, the PDF extracted in this way follows the trend of those given by the global fits in the x > 0.1 region, but does not reproduce their singular behavior in the x < 0.1 region. The latter is usually related to the x −0.5 pattern of the ρ-meson Regge trajectory. Since the ρ-meson is essentially a rather narrow resonance in the ππ system, one should not expect to accurately reproduce the ρ-meson properties in a lattice simulation in which the pions are as heavy as 600 MeV. Thus, one may hope that using simulations at physical pion mass would produce a better agreement with the global fits in the small-x region. This hope is supported by recent extractions [33,34] of q v (x) using the quasi-PDF lattice simulations at physical pion mass. ACKNOWLEDGMENTS I thank J. Karpie, K. Orginos and S. Zafeiropoulos, my collaborators on Ref. [12], who performed the lattice simulations, the results of which were analyzed in that paper and also used in the present work. I am especially grateful to K. Orginos for collaboration on the pseudo-PDF evolution in Ref. [12] and further discussions of this subject. I thank N. Sato for providing the code generating the global fit PDFs and Y. Zhao for discussions