Parton distribution functions from lattice QCD using Bayes-Gauss-Fourier transforms

We present a new method, based on Gaussian process regression, for reconstructing the continuous $x$-dependence of parton distribution functions (PDFs) from quasi-PDFs computed using lattice QCD. We examine the origin of the unphysical oscillations seen in current lattice calculations of quasi-PDFs and develop a non-parametric fitting approach to take the required Fourier transform. The method is tested on one ensemble of maximally twisted mass fermions with two light quarks. We find that with our approach oscillations of the quasi-PDF are drastically reduced. However, the final effect on the light-cone PDFs is small. This finding suggests that the deviation seen between current lattice QCD results and phenomenological determinations cannot be attributed solely on the Fourier transform.


I. INTRODUCTION
Parton distribution functions (PDFs) are fundamental objects that describe the structure of hadrons probing the distribution of momentum and spin among their constituent partons. They also serve as an essential input for collider experiments to obtain the cross-section of a given process. PDFs are extracted from global QCD analyses of worldwide experimental data assuming certain parametrizations that thus imply an intrinsic model dependence. Since PDFs are inherently non-perturbative observables, an ab-initio calculation using lattice field theory methods based solely on QCD is highly desirable. However, such an effort was not possible for a long time due to the fact that PDFs are defined on the light cone and thus they are not accessible by Euclidean lattice QCD simulations.
A way to circumvent this difficulty was suggested by Ji [1] who proposed to use matrix elements from purely spatial correlations that are accessible in lattice QCD. In this way, quasi-PDFs can be connected to the true PDFs through a matching procedure. The method replies on the so-called large momentum effective theory (LaMET), which requires hadron states boosted to large enough momentum. Several works implemented this proposal within lattice QCD with very promising results [2][3][4][5][6] demonstrating the applicability of this approach. A major step forward was the development of perturbative [7] and non-perturbative [8,9] renormalization of the lattice matrix elements. An overview of the current status of lattice QCD calculation of PDFs can be found in Refs. [10,11].
In Ref. [12] a detailed analysis of systematic effects that enter in the computation of quasi-PDFs and need to be investigated was presented. One of the open issues that was identified was the unphysical oscillations seen in the x-dependence of quasi-PDFs, which also affect the final PDFs. A possible origin of these oscillations is the periodicity of Fourier transformation and the fact that a truncation is implemented due to the finite Wilson line. Another possible explanation is that not large enough momentum boosts are currently feasible. It is the goal of this paper to examine whether the discrete Fourier transformation is responsible for the oscillatory behavior and to replace it by evaluating a continuous Fourier transform on the function resulting from a Gaussian process regression (GPR). Although our approach is reminiscent to the one described in Ref. [13], it differs from it in a number of aspects, as discussed in Sec. II.
The rest of the paper is organized as follows: In Sec. III we present the N f = 2 lattice ensemble under study in this work analyzing the main features of phenomenological data for the unpolarized PDFs. In Sec. IV we analyze in more detail the discretization procedure related to the computation of the quasi-PDFs and the Discrete Fourier Transform (DFT). In Sec. V we introduce the main concepts of GPR arXiv:2007.13800v1 [hep-lat] 27 Jul 2020 and explain how to use it in the context of Fourier transforms (FT). In Sec. VI we test the proposed method on a mock data-set generated from a function whose FT is known in closed form. Finally, in Sec. VII we apply the Bayes-Gauss-Fourier transform to a set of lattice data, and in Sec. VIII we draw our conclusions.
In Eq. (1) the state |h(P 3 ) is the boosted hadron with fourmomentum P = (E, 0, 0, P 3 ) and W (0, z) is the Wilson line taken along the boost direction z. The structure of the matrix Γ, acting in Dirac space, determines the three types of PDF, namely Γ = γ 3 is for the unpolarized, Γ = γ 5 γ 3 for the polarized and Γ = σ 3ν for the transversity. From finite momentum lattice QCD measurements ofq(x, P 3 , µ) from Eq. (1), the physical PDF can be obtained through a matching procedure based on Large Momentum Effective Theory (LaMET) [4,[14][15][16][17][18]. This allows to obtain, after a proper renormalization, the physical PDF, in the limit of large enough momenta. The matching formula reads where q(x, µ) is the physical PDF and C(ξ, η) is the matching kernel that depends on the type of PDF. A detailed description of the matching procedure, including the analytic expression of the kernel for unpolarized PDF, is reported in Appendix A. However, on the lattice, renormalized matrix element is only available for a limited range of z-values, where z is the length of the Wilson line entering the non-local operators. Since the theory is defined on a discrete set of lattice sites, the matrix element is known only for discrete values of z. Thus, the integral defining the quasi-PDF becomes a finite sum where a is the lattice spacing In particular, the transformation given in Eq. (4) means that, instead of computing the Fourier transform (FT), we compute an analytic continuation of the discrete Fourier transform (DFT) defined for continuous x values. The DFT frequencies are with N = (2k max /a+1) and ak max = z max /a. However, the FT frequencies relevant for the computation of quasi-PDFs are where the momentum P 3 = 2πP/L, L is the spatial extent of the lattice, P ∈ N and x assumes continuous values in the interval [−1, 1]. The discretization procedure described in Eq. (4) introduces a systematic bias in the quasi-PDF [19]: 1. Knowing the matrix element only at discrete z-values limits the high frequency components to |x| < π aP3 . Therefore, higher frequency components cannot be resolved and they are wrongly measured as lower frequency components below the threshold. This phenomenon is also known as aliasing.
2. The finite spatial extent of the lattice introduces a cut-off on z that is limited by z max ≤ L 2 . A limitation in the number of points reduces the frequency resolution of the discretized Fourier transform defined in Eq. (4). When frequency components of the signal do not correspond to the discrete frequencies, the discretized Fourier transform suffers from a distortive effect known as frequency leakage. Cut-off effects become significant for aP 3 ∼ 1.
These two problems could be solved if it were possible to reconstruct a continuous form of the renormalized lattice matrix element, before evaluating the FT. In this work we use the formalism of Gaussian Process Regression (GPR) [20] to perform both the interpolation and the extrapolation required in order to tackle, respectively the two issues mentioned above.
The choice of GPR for the continuous reconstruction is based on the following reasons: • the interpolation is non-parametric, so it has the flexibility to adapt to any dataset without being restricted to a specific parametrized function; • since GPR is based on Bayesian inference, the information about the behaviour of the function towards infinity (extrapolation) can be incorporated into the prior distribution, and taken into account for the interpolation; • the uncertainties of the measurements are incorporated into the interpolation through Bayes theorem; • it is possible to impose a chosen level of smoothness to the interpolating function; • the result of the interpolation is continuous, defined over whole domain of interest and its Fourier transform is computable in closed form.
An application of GPR to Fourier analysis has already appeared in Ref. [13], where the Fourier transform of the interpolating function is referred to as Bayes-Gauss-Fourier Transform (BGFT). Even though the main concepts remain the same, the procedure that we propose for computing the BGFT is different from the one described in Ref. [13]. In that paper, the prior covariance function of the Gaussian process is determined by the inverse FT of a spectral density model that fits the DFT data. In our case, we tune the hyper-parameters of the prior covariance function using the standard Maximum Likelihood Estimation of type II (MLE-II) procedure described in Ref. [20], since, within this approach, it is easier to implement the prior function the available information about the behaviour of the renormalized matrix element.

III. PROPERTIES OF THE LATTICE MATRIX ELEMENT
We apply the BGFT approach to the MMS renormalized matrix element of the unpolarized PDF. We use the results of Ref. [12] that were computed using a gauge ensemble with two dynamical mass degenerate light twisted mass quarks (N f = 2) generated by the Extended Twisted Mass Collaboration (ETMC) [21]. The lattice spacing is a = 0.0938(2)(3), the lattice size 48 3 × 96 and the sourcesink time separation is fixed at t s = 12a ≈ 1.1 fm. All the relevant parameters regarding the employed gauge ensemble are reported in Table I. This ensemble is referred to as the cA2.09.48 ensemble.  Table I. Simulation parameters of the cA2.09.48 ensemble used to extract the unpolarized quasi-PDF. The nucleon mass (mN ), the pion mass (mπ) and the lattice spacing (a) were determined in Ref. [22].
We show the results on the renormalized nucleon matrix element of the unpolarized operator in Fig. 1 taken from ref. [12]. The momentum is P 3 = 10π/L 1.38 GeV. A total of N conf = 811 gauge configurations were analyzed, with a total number of measurements N measures = 73000. The errors were obtained using the jackknife resampling method. Before presenting the methods used to determine the quasi-PDFs, it is necessary to analyze some important aspects regarding the properties of these lattice matrix element. In a recent paper [12], based on the data from the cA2.09.48 ensemble, the authors identified possible sources of systematic uncertainties in the reconstruction of PDFs from lattice QCD simulations, including the dependence on the cutoff z max .
The results showed that the PDFs extracted are contaminated by a larger noise and stronger oscillations as z max increases above the value of z max /a 10. In particular, the higher the value of z max , the bigger the errors prohibiting the PDFs from reaching a zero value in the large positive x region. However, all the light-cone PDFs extracted using different values of the cut-off z max /a are compatible within errors [12]. We chose to compare the results of the proposed non-parametric regression procedure with the DFT quasi-PDF obtained with z max /a = 10. Let us examine first some important properties exhibited by the renormalized matrix element that will be useful to better design our algorithm. Firstly, they are hermitian functions of z, namely which means that the real part is even and the imaginary part is odd in z. This property will be exploited by our BGFT method in order to make the algorithm more stable, as described in Sec. V D. Secondly, both the real and imaginary parts are expected to decay to zero, with a rate that increases with the nu-cleon boost. This crucial aspect will be reflected in the choice of the prior mean of the GPR, which encodes all the relevant information about the renormalized matrix element, as illustrated in Sec. V C. Thirdly, the real and imaginary parts of the matrix element can be written in polar coordinates, i.e. at fixed z, the complex function h(z) = Re h(z) + i Im h(z) can be written as The function arctan2(y, x) reads In Fig. 2 we show the functions ρ(z) and φ(z) for the renormalized nucleon matrix element of the unpolarized operator. The modulus ρ(z) of the complex function h(z) is an even function that decays to zero with a rate that depends on the nucleon boost. However, due to the increasing errors in the renormalization function, the function ρ(z) deviates from zero as z increases. On the other hand, the argument φ(z) shows a linear behavior for all relevant values of z. Such a property means that the matrix element can be written as or, equivalently This property will be used in the BGFT to improve the results as described in Sec. V E. The linear dependence of φ(z) on z can also be observed directly from the phenomenological PDFs, as discussed in the following section.

A. Analysis of phenomenological data
To illustrate the relation between the real and imaginary parts of the underlying matrix element, we consider the NNPDF3.1 phenomenological determination of the unpolarized PDF [23]. This data set is shown in Fig. 3 and we apply the inverse matching procedure to derive the corresponding matrix element. The inverse matching can be interpreted as the inverse of the operation given in Eq. (3), which allows us to obtain the quasi-PDF from the light cone PDF bỹ The matching kernelC(ξ, η) is reported in Appendix A, while the quasi-PDF for P 3 = 1.38 GeV is shown in Fig. 4.
Having the quasi-PDF, the lattice matrix element can be computed through the inverse FT As can be seen in Fig. 5, the matrix elements extracted from the NNPDF3.1 phenomenological unpolarized PDF show qualitatively the same behavior as that seen in the lattice QCD data. In particular, they decay to zero at sufficient large z/a with a trend that depends on the nucleon boost P 3 . It is interesting to observe that one needs to boost to P 3 = 2.75 GeV so that the matrix element decay to zero for z/a-values larger than 10. Moreover, the z-dependent phase φ(z) behaves linearly within the range of z/a-values for which the real matrix elements is nonzero, as shown in Fig. 6. In order to measure the deviation from the linear behavior of the function φ(z), we present in Fig. 7 the difference between the argument of the matrix element computed at P 3 = 0.83 GeV and a linear fit performed in the interval z/a ∈ [− 11,11]. The obtained curve is compatible with zero, therefore we have elements to presume that the linearity of φ(z) is indeed an intrinsic property of the matrix element that can be exploited to design an efficient regression algorithm. We will provide more details on how we exploit this property in Sec. VII.

IV. DISCRETE FOURIER TRANSFORM
The quasi-PDF given in Eq. (1) is defined as the Fourier transform of the renormalized lattice matrix element of the unpolarized operator. As mentioned in Sec. I, the physical PDF can be obtained from the quasi-PDF by applying a perturbative matching procedure. In particular, at oneloop order, the integral of Eq. (3) consists of a convolution of the quasi-PDF with a function possessing a singularity at x/ξ = 1. Therefore, to compute the light-cone PDF it is crucial to obtain a trustworthy reconstruction of the quasi-PDF for continuous x-values and, in particular, in the region x ∈ [−1, 1]. In what follows we give more details on the issues that arise in the reconstruction of a continuous quasi-PDF originating from the renormalized matrix element. In particular, we show the difficulties in accessing the small-x region having available a limited amount of  Table I.
The problem of reconstructing a continuous momentumspace function starting from a discrete and finite set of momentum-space points is mathematically ill-posed. Indeed, a well-defined transformation is the one that maps . Such a transformation is the Discrete Fourier Transform (DFT), that can be defined as follows with where N = (2z max /a + 1) and ak max = z max /a. However, the DFT assumes periodicity of the matrix element h(z, P 3 ) with a period T that is a function of z max and the DFT frequencies given in Eq. (13) have a meaning only under this property. However, the matrix element is not periodic and Eq. (14) does not hold. In particular, this implies that we cannot attribute to the DFT frequencies of Eq. (13) any special meaning. However, as it will become clear later on, it is still interesting to examine the DFT of the renormalized matrix element. In Fig. (8) we show the DFT corresponding to the matrix element of the unpolarized operator (introduced in Sec. III) with a boost P 3 = 10π/L and z max /a = 10.
As already pointed out, to compute the physical PDF a continuous momentum-space function q(x, P 3 ) is required. A possible solution is to compute the sum in the righthand side of Eq. (4) for continuous x values. The resulting transformation can be written as follows 1 with x ∈ IR. Such a transformation maps an N -point discrete position-space sequence h(z, P 3 ) defined on the interval z ∈ [−z max , z max ] into a continuous momentum-space functionq(x, P 3 ). However, since the matrix element is available up to z max , the definition in Eq. (15) cannot be applied. What we compute instead is the convolution of the matrix element with a function χ I (z) defined by that restricts z in the interval I = [−z max , z max ]. This procedure is equivalent to computing the sum in Eq. (15) within the interval [−z max , z max ], where the matrix element is known. In what follows, we will show that the transform obtained through this procedure is an interpolation of the DFT and, for this reason, will be referred to as interpolated-DFT (iDFT). The results of applying iDFT to the matrix element as compared to the DFT are shown in Fig. (8). The iDFT is the transformation that is being employed to compute the quasi-PDF starting from the renormalized matrix element and, for this reason, it is interesting to analyze in depth the features characterizing this quantity. The matrix element can be expressed as the inverse DFT 1 In signal processing, this kind of transformation is referred to as the Discrete-Time Fourier Transform (DTFT).
whereq(ω k , P 3 ) is the k-th DFT coefficient, and k max = z max /a. Substituting Eq. (17) into Eq. (15), and computing the sum over z we get where ∆ω k = ω − ω k and is the so-called Dirichlet kernel. In particular, in the limit x → 0, the Dirichlet kernel converges to 1, thus the iDFT is equivalent to the DFT for ω = ω k , k ∈ [−k max , k max ].
In summary, given an N −point discrete position space sequence h(z, P 3 , µ), it is not possible to define an appropriate discrete transformation returning an N -point discrete momentum-space sequenceq(ω k , µ). The most straightforward solution is to use analytical continuation for real values of the variable x of the DFT. In particular, the iDFT is the convolution of a step function χ I (z) defined in the interval [−z max , z max ] with the matrix element. This is equivalent to setting the value of the matrix element to zero for |z| > z max and, as it will be shown in Sec. VI, this procedure introduces nonphysical oscillations in the quasi-PDF.

A. Gaussian process regression
A Gaussian process (GP) [20] f (z) is a collection of random variables, labeled by the continuous real index z, that follows a joint multivariate Gaussian distribution with mean function µ(z) and covariance function k(z, z ). This means that, given any number n of domain points z 1 , ..., z n , the random variables f 1 , ..., f n , where f i ≡ f (z i ), are distributed according to: Given the values of the matrix element h(z 1 ), ..., h(z n ), the aim of GPR is to find the most suitable µ(z) and k(z, z ) such that µ(z) is an estimate of h(z), and the covariance k(z, z ) describes the deviations of µ(z) from h(z). If the matrix element h(z 1 ), ..., h(z n ) has zero error, then an appropriate GPR would return a mean function such that µ(z i ) = h(z i ) and k(z i , z i ) = 0 ∀i. For this reason, the resulting GP is also called a surrogate model for h(z).

B. Bayesian inference
The regression is usually performed with Bayesian inference, for which it is needed to specify a prior GP that is then adapted to the measured data computing the conditional probability and using Bayes theorem. The resulting posterior GP is the surrogate model.
In order to specify the prior GP f P (z), it is sufficient to define a prior mean function µ P (z) and a prior covariance function k P (z, z ). The prior mean function represents our belief about the behaviour of h(z). The prior covariance function quantifies the amount of expected deviation from the mean function, and the correlation between those deviations at different values of z.
Denoting by z 1 , ..., z n the z values in which the matrix element measures h 1 , ..., h n are available with uncertainties ∆h 1 , ..., ∆h n and assuming the errors to be Gaussian, it is possible to analytically compute the mean function µ(z) of the posterior Gaussian process f (z) in terms of the prior GP [20]: where K ij = k P (z i , z j ) + ∆h 2 i δ ij .

C. Choice of the prior
It is possible to verify from Eq. (20) that, if the errors ∆h i are zero, then µ(z) passes through all measured points. This property, together with smoothness conditions that can be imposed on GP, as it will be detailed below, ensures that the dependence of the mean posterior on the mean prior becomes weaker as we get closer to measured points. Therefore, the choice of GP prior does not affect the domain areas with high density of measured points. On the other hand, for z values that are far from the measured values, the choice of the prior plays a decisive role.
With the choice of the prior GP, it is possible to impose different levels of smoothness to a GP [24,25]. In particular, choosing µ P (z), k P (z, z ) ∈ C ∞ guarantees that both the GP prior and the GP posterior are infinitely meansquare differentiable, which we expect to be the most appropriate description for our problem at hand.
In order to not compromise the performance of the method and to reduce over-fitting, the prior covariance function is commonly chosen to be stationary and symmetric so that k P (z, z ) = k P (|z − z |, 0). Another common choice is to use a covariance function monotonically decreasing with distance: k P (a, 0) > k P (b, 0) ∀a, b such that 0 < a < b. This last property states that the correlation between predictions at different z 1 and z 2 decreases with their distance |z 1 − z 2 |, which means that the value of µ(z) will be mostly determined by the value of the neighboring measured points.
We opted for the squared exponential covariance function, which is the standard choice for a C ∞ function that satisfy all the properties listed above. We thus consider where the real values σ and , also called hyper-parameters, are fixed using the maximum likelihood estimation of type II described in Ref. [20].
Since the behaviour of the mean function tends to be independent of the measured values at long distances from them, the asymptotic behaviour depends only on the choice of the prior mean function. As mentioned in Section III, we know that renormalized matrix element should tend to zero in the limit of |z/a| → ∞. It is then possible to guarantee this limit for the posterior mean function by choosing a prior mean function that tends to zero at infinity.

D. Strategy for complex Hermitian data
The GPR described until now is defined on real-valued data. There are many possible ways to extend the procedure to complex data. The approach that we choose consists in performing two independent fits: a GPR for the absolute values and a minimum χ 2 linear regression for the complex argument.
Since our target functions are Hermitian, we can restrict the fit procedure to the positive semi axis and then use the Hermitian symmetry to obtain the results on the negative semi axis. With this strategy there is a reduction by half of the number of points used for the non-parametric regression, which improves the stability and the performance of the algorithm. Since the complex argument is an odd function of z, the linear regression is performed with the intercept fixed to zero.
Denoting by µ(z) the result of the GPR for the absolute value and θ the coefficient resulting from the linear regression, the surrogate model obtained for the renormalized matrix element in the positive semi axis reads: The corresponding Hermitian function defined over the full real domain is then: A side effect of this procedure is the loss of continuity of the derivative at z/a = 0. However this is not an issue for our procedure because no subsequent passages require this property to hold.

E. Analytic Fourier Transform
A useful feature of the GPR is the possibility to perform analytically the FT of the posterior mean, obtaining an improved stability and performance compared to what is achievable with a numerical integration.
The FT definition that we adopt is the following: If T is the integral transform defined by Eq. (1), it is possible to write T and T −1 in terms of the FT of Eq. (23a) as follows Thus, after performing the fit, it is possible to estimate the quasi-PDF by computing the FT of the fit function of Eq. (22) using the convention of Eq. (23a), and then by evaluating it using Eq. (24a). When computing the FT, the phase of Eq. 22 simply corresponds to a shift in the FT: In order to compute the FT of µ(|z|), it is useful to observe that Eq. (20) is just a linear combination of the covariance function reported in Eq. (21): where w i ≡ j K −1 ij (h j − µ P (z j )) and µ P (z) is a generic prior mean assigned to the absolute value of the renormalized matrix element.
The FT of k P (|z|, z i ) is available in closed form given by Thus, the quasi-PDF transform of the fit is If the chosen prior mean ρ P (|z|) has a known analytical FT, then the quasi-PDF transform of the fit is available in closed form.

VI. TESTING BGFT ON A MOCK DATA SET
To gain an insight on the artifacts that may be introduced by the discrete Fourier transform and on the effectiveness of the proposed method, we produced a mock dataset that mimics the behaviour of the matrix element shown in Fig. 1. Given the rescaled Gaussian g(x; µ, σ, c) z max =5 z max =7 z max =18 Figure 11. Dependence of the DFT on the cutoff compared to the shifted Gaussian g(x, µ, σ, c) from which we generated the mock data-set.
its inverse quasi-PDF transform of Eq. (24b) reads To be consistent with the results reported in sec. VII, we choose P = 5 and L = 48. The complex function of Eq. (27) is then sampled in the interval z ∈ [−25, 25], z ∈ N, with c = 2.22, µ = 0.315, σ = 0.230. The employed coefficients µ, σ and c correspond to the best fit performed with the function of Eq. (26) on the data obtained by evaluating the iDFT using the discrete set of grid points. The resulting fit is shown in Fig. 9. In order to mimic the behaviour of lattice data, we generated N = 100 numbers at each fixed integer z from a Gaussian distribution centered in T −1 [g(x)](z) with variance increasing linearly with z 2 , obtaining a sample of N mock matrix element. The average and the Jackknife standard deviation of this sample are shown in Fig. 10. The dependence of the discrete Fourier transform on the cut-off z max is investigated. In Fig. 11 we show the DFT computed with four different values of z max , together with the shifted Gaussian g(x, µ, σ, c) from which we generate the mock data-set. In particular, setting z max = 5, a huge bias is introduced in the DFT and big oscillations afflict the final result. Moreover, the bias becomes negligible only with z max = 18, where the iDFT coincides with the analytical FT within the error. The observed behavior of the iDFT is due to the fact that, if the z cut-off is too small, then the frequency resolution of the DFT is not fine enough to capture the behaviour of the analytical FT in the small x region. Considering the case of z max = 7, we apply the regression described in Sec. V to the mock data set, using a constant zero function as the prior mean function. The results are illustrated in Fig. 12. The conclusion is that, given the mock matrix element up to z max = 7 and their asymptotic behavior specified by the prior as a zero con-    stant function, the non-parametric regression is able to reproduce the data outside the fit range. As a consequence, the results after applying the Bayes-Gauss-Fourier transform shown in Fig. 13, are compatible with the analytical transform of the mock matrix-element within error.

VII. APPLICATION OF BGFT TO THE RENORMALIZED LATTICE MATRIX ELEMENT
We apply the method detailed in Sec. V on the case of interest, namely the renormalized nucleon matrix element of the unpolarized operator, for which the data are given in Sec. III. As we pointed out in Sec. V C, in the region with high density of data-points, the posterior mean is strongly dependent on the renormalized matrix element rather than on the prior mean. For this reason, we exclude from the fit the region of the nonphysical negative values of Re h(z), performing the non-parametric regression up to z max /a = 7.
As mentioned before, it is possible to incorporate information about the objective function into the prior mean function, and this choice is particularly relevant to determine the asymptotic behaviour. There is not much a priori knowledge about the function reproducing h(z, P 3 ). For this reason, we test our method with two different prior mean functions: the uniformly zero function and the function of Eq. (27) obtained through applying the inverse transform to the Gaussian fit of the DFT. Using different prior distributions is a method to cross-check that the final conclusions are independent of the prior choice.
Let us here summarize the key steps of the procedure: 1. at a fixed z-value we rewrite the complex number h(z) = Re h(z) + i Im h(z) in the polar complex plane as with ρ(z) = Re h(z) 2 + Im h(z) 2 φ(z) = arg(h(z)) = arctan2 (Im h(z), Re h(z)) .
The function arctan2(y, x) is defined in Sec. III 2. as pointed out in Sec. III, the function ρ(z) is asymptotically zero, while φ(z) can be taken as a linear function of z. After choosing a prior mean function, we perform a non parametric regression of the function ρ(z), while a linear fit is sufficient to reproduce φ(z), as shown in Fig. 14; 3. in order to check the result of the fit to the renormalized matrix element, we can go back to the Cartesian coordinates, as shown in Fig. 15; 4. employ the formula in Eq. (25) to compute the quasi-PDF.
In the upper panel of Fig. 16 we compare the iDFT quasi-PDF to the BGFT quasi-PDF. While in the physical region x ∈ [−1, 1] the two results are compatible, for larger |x| the non-physical oscillations due to the periodicity of the discrete FT are strongly suppressed. However, the physical meaning of quasi-PDF can be made explicit only after having performed the matching procedure. In Fig. 18 we display the light-cone PDF reconstructions obtained via iDFT and BGFT. As can seen, although the the non-physical oscillations in the quasi-PDF are suppressed, in the physical PDF the effect is small. The nonphysical negative PDF in the antiquark region x ∼ −0.1 remains, as well as a mild oscillatory behaviour in the large |x| region. This means that this behavior does not appear to be caused by the cut-off in z and the discrete FT. Finally, it is interesting to investigate how the nonparametric regression curves for the real and imaginary parts of the matrix element go to zero. In particular, both the tail of the real and imaginary parts of the nonparametric regression curves can be modeled by the func- The parameters a and b have been computed by minimizing χ 2 . We obtained a = 25(16), b = 9.4(9) fm −1 for the real part and a = 270(80), b = 13.8(7) fm −1 for the imaginary part.
As previously stated, we use as an alternative prior mean the function of Eq. (27) with µ = 0.315, σ = 0.23 and c = 2.22 in order to cross-check our results. The outcome of the non-parametric regression with non-zero mean prior is shown in Fig. 19, while the lower panel of Fig. 16 shows the FT of this function together with the resulting BGFT. As pointed out in Sec. V C, the choice of the prior mean function modifies the result of the GPR in the region where there is a low density of data points. In this specific case, it slightly modifies the decay rate in the large z region, further reducing the amplitude of the remnant oscillations present in the BGFT. However, the effect of the prior mean is not observable in the light-cone PDF that is still compatible with the reconstruction obtained with iDFT, as shown in the lower panel of Fig. 18.

VIII. CONCLUSIONS
In this work we address the nonphysical oscillations, which appear in the computation of PDFs from lattice QCD simulations. Due to the lattice discretization, the continuous Fourier transform in Eq. (1) cannot be computed. Moreover, the continuous FT cannot be replaced by the discrete Fourier Transform, since this would require periodicity of the matrix element. To obtain a continuous reconstruction of the FT, we can employ the analytical continuation of the DFT, defined for continuous real values of x. However, this transform, would require the knowledge of the discretized matrix element for arbitrary large values of the Wilson line length |z|. The commonly adopted solution in the literature is to assume that the matrix element goes to zero for |z| > z max . This assumption leads to the definition of the interpolated DFT (iDFT) that, in contrast to the DFT, provides a continuous frequency domain function that is suitable for computing the light-cone PDF. However, the iDFT is only an interpolation of the DFT, consisting of a linear combination of Dirichlet kernels. As a consequence, despite being a continuous function, it is still afflicted by the same problems that appear in DFT (aliasing and leakage) that hinder the evaluation of the transform. These considerations are in accordance with the results obtained in a recent paper [19], and suggest that the non-physical oscillations observed in the PDFs computed from lattice QCD matrix element may be due to the discretization of the Fourier Transform. The problems afflicting the iDFT can be solved by reconstructing a continuous form of the renormalized matrix element defined over the whole domain before evaluating the continuous FT. In order to obtain such a continuous reconstruction we employ Gaussian Process Regression. It consists of a Bayesian non-parametric regression that leverages on the smoothness properties of the renormalized matrix element function, and its asymptotic behaviour. Another property that makes GPR a useful tool for quasi-PDF computations is that the Fourier transform of the regression result is analytically computable. We demonstrate the applicability of this approach in Sec. VI using a mock data set generated from a function whose FT is known in closed form. Even though this mock data set possesses the same limitations characterizing lattice QCD data, the Bayes-Gauss-Fourier Transform is able to capture all the relevant feature of the analytical FT.
The method is applied to the MMS renormalized matrix element of the unpolarized PDF, computed on the cA2.09.48 ensemble with N f = 2 flavours of quarks, lattice size 48 3 × 96 and a source-sink time separation t s = 12a ≈ 1.1 fm. The BGFT shows a significant reduction of non-physical oscillations in the large |x| region, while it is compatible within error bars with the DFT transform for x ∈ [−1, 1]. However, the improved behaviour of the quasi-PDF is small for the physical PDF, where no substantial deviation can be detected as compared to the PDF obtained with the iDFT.
This finding suggests that the presence of the nonphysical negative values in the light-cone PDF for x < 0 cannot be ascribed to the discrete Fourier transform, or at least, this cannot be the sole cause of this behaviour. Thus, it seems to become mandatory to either reach higher nucleon boosts in lattice calculations or compute higher perturbative orders for the matching and the conversion between renormalization schemes. Also, higher twist effects need to be understood. Nevertheless, in this paper we have presented an alternative way to analyze the renormalized matrix element, which lead to quasi-PDFs with suppressed oscillation. This approach can thus provide a valuable cross-check of lattice computations of parton distribution functions.
Moreover, in Sec. III A we apply the inverse matching procedure to obtain the matrix elements starting from the phenomenological PDFs. Given the light-cone PDF, the quasi-PDF operator can be computed as q (y, µ, P 3 ) =