Lattice calculation of hadronic tensor of the nucleon

We report an attempt to calculate the deep inelastic scattering structure functions from the hadronic tensor calculated on the lattice. We used the Backus-Gilbert reconstruction method to address the inverse Laplace transformation for the analytic continuation from the Euclidean to the Minkowski space.


Hadronic tensor on the lattice
There has been a lot of interest and developments in calculating the structure functions and parton distribution functions on the lattice in recent years [1]. The Euclidean hadronic tensor has been formulated in the path-integral formalism [2][3][4]. The hadronic tensor in deep inelastic scattering, from which the parton distribution functions are obtained through the factorization theorem, can be calculated from its Euclidean counterpart via an inverse Laplace transform [1,2,4]. It has been revealed that, in addition to the valence partons, there are two types of sea partons -connected sea (CS) and disconnected sea (DS) in three topologically distinct path-integral diagrams. The extended evolution equations to accommodate both the connected and disconnected sea partons are derived [5]. It is essential to have separately evolved CS and DS partons so that comparison with lattice calculations of unpolarized and polarized moments of PDF can be made. Only with the extended evolution equations will the CS and DS partons remain separated at different Q 2 to facilitate global fitting of PDF with separated CS and DS partons.
where |p is the nucleon state and J µ is the vector current. In the Euclidean path integral formalism, the hadronic tensor related nucleon matrix element can be expressed as the ratio of the four-point function and the two-point function. To be specific, they are where E N and m N are the energy and mass of the nucleon and Γ e = 1+γ 4 2 is the unpolarized spin projector. After inserting the complete set of intermediate states, we havẽ We see that in Eq. (5) there is an exponential dependence on the Euclidean τ. It will exponentially decay when the lowest E n is heavier than E N and it will exponentially grow when the lowest E n is lighter than E N . To analytic continue to the Minkowski space, an inverse Laplace transform is needed, i.e.
with c > 0. This is basically doing the anti-Wick rotation back to the Minkowski space to recover the delta function in energy as shown in Eq. (1). The topologically distinct insertions are shown in Figure 1. The first three involve leading-twist contributions from the valence +CS partons in (a), the CS antipartons in (b), and the DS partons in (c). The last two insertions ((d) and (e)) are higher-twist contributions which are suppressed by O(1/Q 2 ) and will be ignored.
We will focus on the first two insertions in this work. It has been pointed out that the second insertion are the connected seaū andd contribution, which are responsible for the Gottfried sum rule violation [2,4]. It is noted that this approach has the advantage that no renormalization is needed when the conserved vector current is used and only finite renormalization is needed for local vector current. Furthermore, the structure functions from the hadronic tensor are frame independent so that they can be calculated with the nucleon at rest or at low lattice momenta. However, the inverse Laplace transform, which is needed to convert the hadronic tensor from the Euclidean space to Minkowski space, is a challenge. We shall test a numerical approach to address the inverse problem.

Backus-Gilbert reconstruction
The most challenging part of our calculation is to convert the Euclidean hadronic tensorW(q 2 , τ) as a function in τ to the Minkowski W(q 2 , ν) in ν through certain numerical approximation. This is a typical case of the inverse problem. We shall consider the Backus-Gilbert reconstruction method as suggested in [6].
A general form of the inverse problem is where c i are the discrete data points that we have, k i (x) are the corresponding integral kernels and w(x) is the desired continuum function. In principle, one cannot get a unique solution for w(x) since the number of input data points is always finite. In other words, we do not have all the information needed to reconstruct every detail of the continuum function. Even though the continuum function w(x) can be discretized into w j , the number of points of w j we need will in general be larger than that of c i . So the inverse problem is, in principle, an ill-defined problem. However, we are still able to find the most probable solution to the problem by using some algorithms, such as the maximum entropy method [7] and the Backus-Gilbert method [8,9]. In our case, the input data are τ-dependent Euclidean correlation function, the kernel is the Laplace kernel and the function ω(ν) is the spectral density Another reason to use the Backus-Gilbert method is when considering the exclusive multi-hadron final states on the lattice, the finite volume effects must be handled correctly. Based on [6], the regularized delta-function introduced in the Backus-Gilbert method will ensure that the reconstructed spectrum density function has a well-defined infinite volume limit.

Numerical details
The central idea of the Backus-Gilbert method is that a linear combination of the kernel functions can be used to approximate the delta function, Assuming the equation holds, we will have If we can compute a(ν 0 ) ( a means the coefficients vector (a 0 , a 2 , ...) T ) for each ν 0 and for each ν 0 Eq. (9) always holds, then we can have the value of the function ω(ν) at arbitrary ν, which means that the inverse problem is solved. Practically, the finite number of kernel functions cannot span a complete function base and, therefore, the delta function approximated form the linear combination cannot be a true delta function, but a broadened Gaussian-like function instead (the regularized delta function). The more data points we have, the more sharp the Gaussian-like function will be. The sharpness of the Gaussian-like function actually embodies the resolution of the reconstruction. This is consistent with our intuition, the more data points we have, the more information we can extract.
On the other hand, to solve the coefficients a(ν 0 ), we can first define a function K(ν 0 , a) of a as the "deltaness" of the Gaussian-like function, and the coefficients a which minimize the value of K will serve as the desired a(ν 0 ). Therefore to solve for the coefficients a(ν 0 ) we need to solve the linear equations Moreover, to remove the trivial solution of a = 0, a Lagrange multiplier is added into the equations. The error of the Backus-Gilbert result σ(ν 0 ) is defined as where (Cov) is the covariance matrix of the input data. Since in general the errors of the data are not negligible, the very details of the data are not accurate (altered randomly by the noise), one needs to introduce an SVD cutoff to remove these imprecise small eigenvalues when solving the linear equations. This SVD cutoff is the only tunable parameter when using the Backus-Gilbert method. When the cutoff is set to be small, more eigenvalues will be kept and thus better resolution one can have. However, the small eigenvalues above the cut will result in larger a t values and, consequently, larger final errors from Eq. (13). On the other hand, when the cutoff is set to be larger, less eigenvalues will be kept and thus the resolution will be worse, but the errors will be smaller and the results will be more stable. So this leads to a trade-off. In practice, the SVD cutoff is chosen according to the errors of our data. We try to find a stable window of the cutoff which gives consistent results and then choose the one producing the smallest errors in that window.

Toy model tests
To test the Backus-Gilbert method, first, we generate a mock two-point correlation function containing three mass terms without noise by We set A i = 20, 40 and 60 and m i = 0.2, 0.6 and 1.2 respectively. The number of time slices used is 25, which is close to that of the preliminary lattice calculation in the following section. The Backus-Gilbert reconstructed results are shown in the left panel of Figure 2. We can see the three peaks of the spectral density function are indeed located close to the input positions of m i (vertical lines in the plot), which means the Backus-Gilbert method extracts the correct spectral information from the toy two-point function. We can also estimate the spectral weight A i from the spectral function. To do so, we use Gaussian functions to approximate the peaks. We know that the integrals of Gaussian For convenience sake, we choose to use the width-at-half-maximum w H i = 2.773 σ i to replace σ i , then   Next, we change the setup to a more relevant one: one isolated state (m 0 = 0.6 and A 0 = 200) plus a section of dense spectrum (simulating the continuum spectrum in the DIS region) and also some noises. The dense spectrum is form ν = 1.5 to ν = 10 with interval δν = 0.01 and weight factor A = 0.5. The signal-to-noise ratio is set to be 100 at all time slices. In this case, since the noise will cover some details of the correlator, we need to use an SVD cutoff of 10 −6 to remove the inaccurate small eigenvalues. As show in the right panel of Figure 3, the peak shows itself in the right position; the near constant tail is stable and the value is commendably consistent with the expected value that is A/δν = 50. Again, for the spectrum weight of the isolated state, we can read A 0 ∼ 280 and w H 0 ∼ 0.7 from the plot and the estimated A 0 ∼ 208 is consistent with the input. These two toy model tests enable us be confident in the feasibility of applying the Backus-Gilbert method in our hadronic tensor calculation. We will show the Backus-Gilbert results calculated from our preliminary lattice data in the next section. Moreover, we will also consider other inverse methods such as the maximum entropy method and χ 2 fitting of model spectral function in the future to see whether we can have consistent results. This will help to estimate the systematic error of the inverse Laplace transform.

Preliminary results
Our exploratory calculation is carried out on a small 12 3 × 128 anisotropic clover lattice generated by the CLQCD collaboration [10] with m π ∼ 640 MeV, a s = 0.1785(53) fm and ξ = 5. The number of configurations used is 500.  We choose the direction of the currents µ = ν = 1, so what we calculate is the one-one component of the Euclidean hadronic tensorW 11 ( p, q, τ). We use two sequential sources (one from t 0 through t 1 to t 2 while the other from t 0 through t f to t 2 ) to calculate the four-point functions. Therefore, 3 times of inversions are needed for each p, q, t 1 and t f . The Chroma software [11] is used to implement our calculation. The results ofW 11 are shown in Figure 3 as a function of the Euclidean time τ. The left plot is the connected sea anti-parton case from Fig. 1 (b) with p = 0. We see that both curves ofū andd go down as τ increases since the energy of the lowest intermediate state (i.e. nucleon energy with a momentum of p + q) is larger than that of the nucleon at rest. The right one is the valence plus connected sea case with p 0 and q = − p. Both curves go up as τ increases since the lowest intermediate state is the nucleon at rest whose energy is lower than that of the initial state with a moving nucleon (i.e. m N < E p .) The Minkowski hadronic tensor W 11 (q 2 , ν), up to a constant, is just the spectral density ω(ν) in the expression of the Backus-Gilbert method. After doing the Backus-Gilbert reconstruction using the data in the p = 0 case, we get the Minkowski hadronic tensor W 11 as shown in the left plot of Figure 4.
The lowest peaks are the elastic peaks at ν = E p+ q − E p as expected. The second peaks are related to the quasi-elastic peaks of the γN reaction which are the excited nucleon states which includes Roper state, S 11 , πN, etc. When the deep inelastic region sets in as in the experimental γp total cross section σ γp tot in the right panel of Fig. 4, the density of the spectral function is expected to be smooth,  reflecting the parton nature of DIS. For the DIS, the kinematics requires ν < | q| so that Q 2 > 0 and 0 ≤ x = Q 2 2m p ν ≤ 1. Furthermore, ν should be several GeV above the quasi-elastic peaks to be in the DIS scaling region. When we examine Fig. 4 (a), we find that when the quasi-elastic peak tapers off at ν ∼ 1 which is ∼ 5.5 GeV on this lattice, | q| = 3 corresponds to 1.7 GeV. Thus, Q 2 < 0 in this case, there is no room for the physical DIS region.
We also analyze the q = − p 0 case , the corresponding Euclidean data onW 11 ( p, q, τ) are on the right panel of Figure 3. The Backus-Gilbert reconstructed results are shown in Figure 5. The elastic peaks are at ν < 0 since E p+ q < E p . However, again the kinematics of the setup still cannot bring us to the physical region of Q 2 and x, the Q 2 at the tail of the inelastic peak at ν ∼ 0.5 is still negative.
It is clear from these preliminary studies that one needs to have a large | q| so that there is enough room for the energy transfer ν to be less than | q| and yet a few GeV above the inelastic peak. To achieve this, one needs lattices with smaller spatial lattice spacings in order to have a relatively large