Towards the nucleon hadronic tensor from lattice QCD

We present the first calculation of the hadronic tensor on the lattice for the nucleon. The hadronic tensor can be used to extract the structure functions in deep inelastic scatterings and also provide information for the neutrino-nucleon scattering which is crucial to the neutrino-nucleus scattering experiments at low energies. The most challenging part in the calculation is to solve an inverse problem. We have implemented and tested three algorithms using mock data, showing that the Bayesian Reconstruction method has the best resolution in extracting peak structures while the Backus-Gilbert and Maximum Entropy methods are somewhat more stable for the flat spectral function. Numerical results are presented for both the elastic case (clover fermions on domain wall configuration with $m_\pi\sim$ 370 MeV and $a\sim$ 0.06 fm) and a case (anisotropic clover lattice with $m_\pi\sim$ 380 MeV and $a_t\sim$ 0.035 fm) with large momentum transfer. For the former case, the reconstructed Minkowski hadronic tensor gives precisely the vector charge which proves the feasibility of the approach. While for the latter case, the nucleon resonances and possibly shallow inelastic scattering contributions around $\nu=1$ GeV are clearly observed but no information is obtained for higher excited states with $\nu>2$ GeV. A check of the effective masses of $\rho$ meson with different lattice setups indicates that, in order to reach higher energy transfers, using lattices with smaller lattice spacings is essential.


I. INTRODUCTION
In scattering processes involving nucleons such as deep inelastic scattering (DIS) and neutrino-nucleon scattering at low energies, the hadronic tensor W µν is used to characterize the nonperturbative nature of the nucleon structure. It is the imaginary part of the forward virtual Compton scattering amplitude W µν = 1 2π ImT µν and can be expressed as a nucleon matrix element with two current operators inserted The hadronic tensor of the vector currents can be further decomposed, according to its Lorentz structure, into structure functions, i.e., for the unpolarized case, wherep µ = p µ − p·q q 2 q µ . p µ and q µ are the nucleon 4-momentum and momentum transfer, respectively. x is the Bjorken x = Q 2 2p·q . The structure functions are valuable quantities which reveal the inner structure of the nucleon. They can be used to extract parton distribution functions (PDF's) through the QCD factorization theorem F i = a c a i ⊗ f a , where the convolution kernel c a i is perturbatively calculable. Due to their nonperturbative nature and importance, it is natural to explore the possibility of calculating the hadronic tensor and structure functions with Lattice QCD, which is a first-principle nonperturbative method of solving the strong interaction. While the proposal for a lattice QCD based evaluation has been put forward more than 20 years ago [1,2], it has only recently become feasible to compute the necessary 4-point correlation functions, thanks to the increases in computing power [3,4]. In recent years, there has been a lot of effort in the lattice community focusing on the computation of x-dependent PDF's. Examples are Quasi-PDFs and LaMET [5,6], Compton amplitude [7], Pseudo-PDFs [8,9] and Lattice cross sections [10,11]. Each of the approaches has its own advantages and difficulties. Since the hadronic tensor is scale-independent and the structure functions are frame-independent, the lattice calculation of the hadronic tensor has its unique advantages in that no renormalization nor large nucleon momentum are needed. However, to convert the hadronic tenor from Euclidean space to Minkowski space involves an inverse problem [3] which presents a substantial numerical challenge. In case that reliable lattice results for the Minkowski hadronic tensor are obtained for a set of kinematic setups, they can be used to obtain parton distribution functions via the factorization theorem, as are carried out in global fittings of experiments.
Another feature of calculating the hadronic tensor on the lattice is that it reveals explicitly the connected-sea anti-parton contribution. It has been pointed out [1,12] that the Gottfried sum rule violation (i.e., theū andd difference of PDF's) can be explained by the existence of the connected-sea anti-partons. Recently, the ratio of strange to u/d momentum fraction in disconnected insertions was calculated [13], which helps to separate the connected and disconnected-sea parton distributions in global fittings. The calculation of hadronic tensor using Euclidean 4-point correlation functions will provide a direct proof of the degrees of freedom of connected-sea anti-partons and can finally resolve the puzzle of Gottfried sum rule violation.
In addition to deep inelastic scattering, the hadronic tensor plays an important role too for scatterings at lower energies. One example is the experiments of neutrino-nucleus scattering, e.g. LBNF/DUNE [14] at Fermilab, which aims to study the neutrino properties.
These experiments face several challenges like the reconstruction of the neutrino beam energy and flux and the consideration of the nuclear effects and models [15]. In view of this, the input of accurate determination of the neutrino-nucleon scattering is vital to investigating the nuclear effects of neutrino-nucleus scattering. However, it is not trivial to study the neutrino-nucleon scattering since at different beam energies, different contributions (elastic (EL), resonance (RES), shallow inelastic (SIS) and DIS) dominate the total cross section [16]. Nevertheless, the hadronic tensor is useful in all the energy regions. For example, in the EL region of neutrino-nucleon scattering which is relevant to the quasi-elastic neutrinonucleus scattering, the hadronic tensor is actually the square of the elastic form factors of the nucleon and as a result, the cross section of the neutrino-nucleus scattering can be calculated by combining the nucleon form factors and nuclear models about the nucleon distribution inside a nucleus. In the RES, SIS and DIS regions, inelastic neutrino-nucleon scatterings emerge and one will need to have the hadronic tensor to cover all the inclusive contributions.
In this sense, calculating the hadronic tensor is so far the only way we know that Lattice QCD can serve the neutrino experiments in the whole energy range [17] 1 .
Since Lattice QCD is formulated with Euclidean time and the hadronic tensor involves a 4-dimensional Fourier transform, one cannot calculate the hadronic tensor directly on the lattice. Instead, we calculate its counter part in Euclidean space and then convert it back to Minkowski space. The formalism of Euclidean hadronic tensor is discussed in Sec. II. The conversion to the Minkowski space is implemented by solving the inverse problem which is the most challenging part of our calculation. We will discussed three methods and two examples in Sec. III. Sec. IV presents numerical results for both the elastic case and a case with large momentum transfer. Discussion on the results comes in Sec. V.

II. LATTICE FORMALISM OF HADRONIC TENSOR
After inserting a complete set of intermediate states in Eq. (1) and carrying out the integral, one comes to the following expression of the hadronic tensor [3]: where q is the momentum transfer, p is the nucleon momentum and p n is the momentum of the and after the integration one has where ν is the energy transfer, E p and E n are the energies of the external nucleon and the is not a good approximation to the delta function and one cannot pick out the clean contribution for a specific ν as nearby states will mix. Instead, we construct the following 4-point correlation function with only a 3-dimensional Fourier transform and the normal nucleon 2-point function as Then, the Euclidean hadronic tensor W E µν ( p, q, τ ) is defined by the ratio of the 4-point function to the 2-point function [1][2][3][4]12] where E p and m p are the energy and mass of the nucleon. We can insert the intermediate states again between the two currents and we have where and τ = t 2 − t 1 .  (Fig. 1a), the CS anti-partonsq cs (Fig. 1b) and the disconnected sea (DS) partons and anti-partons q ds +q ds (Fig. 1e) [1,2]. In our approach, they can be calculated separately which is a great feature especially for the CS anti-partons that are responsible for the Gottfried sum rule violation [1,2,12]. After the Euclidean hadronic tensor is calculated, we need to convert it back to Minkowski space to obtain physical results. Formally, the inverse Laplace transform fulfills this objective: However in practice, the Euclidean hadronic tensor is a function of Euclidean time, which is real, so that the integral in the inverse Laplace transform along the imaginary time axis is not possible. Numerically, one can try to solve the inverse problem of the Laplace transform to determine an estimation of W M µν [3,4], Details about solving the inverse problem are discussed in the next section.

III. SOLVING THE INVERSE PROBLEM
A general form of the inverse problem reads where c(τ i ) denotes discrete lattice data with finite number of points (usually O(10)), k(τ i , ν) is the integral kernel that is a function of both τ i and ν, and ω(ν) is the target function which is usually continuous with respect to ν. In principle, determining every detail of a totally unknown continuous function with finite input information is not possible; videlicet, more than one solution can be found to match the input data. Numerically, we can discretize however, the number of ν j one needs to reproduce the structures of ω(ν) is, in many cases, much larger than the number of input points, so the problem is still ill-posed. Nevertheless, many algorithms are available to extract the most probable solution at a certain resolution.
Actually, this is a common problem, not only in physics, and the algorithms have been kept updated and improved.
In this section, we will briefly introduce three methods of solving the inverse problem, discuss their features and use some mock data to test their resolutions and robustness.

A. Backus-Gilbert method
The Backus-Gilbert (BG) method [18][19][20] utilizes the fact that the kernel functions can be linearly combined to approximate the Delta function, if they span a complete function where a(τ i , ν 0 ) are the coefficients for the i'th kernel function at a specific point ν 0 , which can be calculated by assuming a criterion of "deltaness" and solving the linear equations.
Having a(τ i , ν 0 ), the value of the target function at ν 0 is It is worthwhile noting that the number of independent kernel functions is equal to the number of the discrete lattice data points so usually the function basis is far from complete, leading to a coarse resolution. In some sense, this is a feature instead of a disadvantage, since the broadened delta function ("regulated delta function") can be treated as a smoothing procedure and ensures a well-defined infinite volume limit [20]. This is because the lattice spectrum is discrete and the volume correction for multi-particle states is not negligible.
However in our case, we are only concerned about the inclusive contribution rather than the contribution from one particular state, the infinite volume limit is well-defined. And in fact, it is nearly impossible to isolate the multi-particle-state contributions even with other methods of better resolution from the present day lattice data. Another feature of BG is that, unlike other Bayesian-type methods, it solves the problem point by point and dose not guarantee that the reconstructed result reproduces the input data [21], so careful checks are alway needed.

B. Maximum Entropy method
The Maximum Entropy (ME) method [22,23] makes use of the Bayesian probability with prior information about the target function to find the most probable solution where P [ω|D, α, m] denotes the probability that ω is the solution given lattice data D, prior information m and a hyper parameter α. Q = αS − L is a combination of the Shannon which entails the constraint from the prior information, and the likelihood function with c ω (τ i ) being the correlator reconstructed using Eq. (14) and C the covariance matrix, which embodies the constraint from the data. m(ν j ) is the default model (the prior information we plug in) and α is the weight that balances the two constraints. If α is zero, ME reduces to the normal χ 2 fit which has no unique solution for the inverse problem, since the number of parameters is larger than the number of input data. The uniqueness is guaranteed for finite α and the results of different α's are averaged in a range of α based on certain assumptions [23]. Practically, the parameter space for finding the maximum probability P [ω|D, α, m], i.e. the maximum value of Q, is reduced to a smaller one instead of the whole ν space by employing singular value decomposition [23,24], which makes the maximum search easier while the resolution may be affected. Improved ME with an extended search space [25] is also proposed and we will check whether it produces better results in our future study.

C. Bayesian Reconstruction
Bayesian Reconstruction (BR) [26] is an improved Bayesian method. The Bayesian probability is where which is an alternative way to encode the constraint from the prior. γ is a numerically large number such that the term γ(L − N τ ) 2 helps to prevent over-fitting. The hyper parameter α here is integrated over as   [21]. The ringing of the ME method seems weaker, which is due to the fact that SVD contains an additional smoothing that suppresses the ringing but also leads to the significantly larger width of peak structures. Both the BR and ME methods suffer from this disadvantage, so they are not the optimal method for reconstruction of flat spectral functions. Actually, there is an update of BR aiming to address this problem [21,27] and we will also try to test this in our future study. Similar to the first case, the right panel of Fig. 3 shows that the regenerated 2-point functions from the spectral functions of ME and BR are well consistent while that of BG is not.
Solving the inverse problem is the most challenging part of our calculation. More detailed studies on how the lattice spacing, the error of the correlation functions and the number of time slices affect the reconstruction and more inverse methods regarding this physics problem are needed. Recently, several inverse algorithms have been applied to evaluate the efficiency of obtaining x-dependent PDF's from mock Euclidean time correlators [28]. From the above two tests we know that BR has the best resolution for peak structures while BG and ME are more stable for flat spectral functions. Knowing the different methods' advantages, one can combine them to resolve different parts of the spectrum, e.g., BR for the sharp peak structures and ME for the flat region.

A. Elastic case
Having discussed how to solve the inverse problem, we now apply the algorithms to realistic lattice data. The first example is to check the vector charge. This is a calculation that is relevant to neutrino-nucleon scattering and also serves as a benchmark of the whole approach. The calculation is done on RBC/UKQCD domain wall lattice 32Ifine [29] with clover fermions as valence quarks. The configuration is preprocessed by HYP smearing and the tadpole improved clover coefficient C sw = 1.033 is used to generate the clover term of For simplicity, the two currents are both inserted on the d quark line so only Fig. 1a contributes. And for τ 0 only the ground state survives, so W E 4,4 = g 2 V = 1, given proper normalization factor Z V . Two sequential propagators are used for constructing the 4-point function with one starting from t 0 through t 1 to t 2 and the other starting from t 0 through t f to t 2 (Fig. 1a). Therefore, for each calculation, the source point t 0 and the two sequential points t 1 and t f are fixed while all values of t 2 are available. In this particular example, we choose t 0 = 0, t f = 15 and t 1 = 5 in lattice unit so t 2 should be in the range of [6,14] to exclude the contact points and the corresponding τ = t 2 − t 1 is from 1 to 9. The result of Fig. 4a. It shows that within errors, W E 4,4 (τ ) is indeed a constant of value 1. The drop at τ = 9 is likely due to the fact that it is too close to the nucleon sink.
The errors are around 0.4% and the number of configurations used is 100.  Figure 4: The results of the elastic case. Fig. 4a: the Euclidean hadronic tensor W E 44 as a function of τ . Fig. 4b: the Minkowski hadronic tensor W M 44 as a function of energy transfer ν from the ME method. Fig. 4c: the Minkowski hadronic tensor as a function of ν from the BG method.
To convert the results to the Minkowski space, the ME and BG methods are employed.
Actually, the exact form of hadronic tensor for elastic scatterings [30] is In our case W M 44 (ν) . This is easy to understand since the spectral function should be a Dirac delta function at δ(ν − ν 0 ) when the Euclidean correlator is a single exponential ∼ e −ν 0 t and here the constant is a special case of an exponential with ν 0 = 0.
The converted results of W M 44 (ν) using ME and BG are plotted in Fig. 4b and Fig. 4c, respectively. They both give a peak around ν = 0 and in this sense the results are consistent with the theoretical prediction of δ(v). However, similar to the cases of the mock data, ME shows much better resolution than BG. Another problem of the result of BG is that it is not symmetric about ν = 0, which is because BG has difficulties in resolving the results of negative ν. An important check is that the area under the peaks should be g 2 V = 1. The values of numerical integral of the results from BG and ME are 1.18(6) and 1.001 (7) respectively. Again, ME shows a more precise result. Although it is not necessary and cumbersome to calculate the vector charge by constructing the 4-point function and solving the inverse problem, it nevertheless shows the feasibility of our approach. The vector charge can be obtained reliably. For more complicated cases such as non-zero momentum transfers or charged currents, this approach may show its advantages and provide the inclusive contribution of all intermediate states.

B. Non-zero nucleon momentum and momentum transfer
As pointed out in the introduction, another important motivation of calculating the hadronic tensor is to have the lattice results of structure functions in the DIS region which can be used together with experimental inputs to better pin down the parton distribution functions. To this end, we need to have large momentum transfers to make the scattering "deep" enough to access the parton degrees of freedom. Meanwhile, the momentum of the external proton cannot be too small (and it is better to be in the opposite direction of the momentum transfer) if one wants to reach small x (e.g. ∼ 0.1). For this calculation, we use an anisotropic clover lattice [31] with a t ∼ 0.035 fm. The pion mass is about 380 MeV and the momentum unit is 2π Ls ∼ 0.42 GeV. The reason we switch to this lattice is that the signal-to-noise ratio will be much worse than the previous elastic case when the momentum transfer or the nucleon momentum is large. Thus, having more data points in the t direction helps the inverse algorithms to have more stable results.
The detailed kinetic setup is listed in Table. I. We choose p = (0, 3, 3) and q = (0, −6, −6) in lattice unit and µ = ν = 1, such that only the F 1 structure function survives; thus . Since the energy transfer ν is not fixed by the lattice 3-dimensional Fourier transform, we can choose a range of ν ∈ [2.96, 3.68] GeV such that the corresponding Q 2 is in a range of 2 to 4 GeV 2 . The Bjorken x that can be accessed is between 0.07 and 0.16 for this setup. An interesting point of this setup is that p + q = − p, therefore the energy of the lowest intermediate state E n=0 = E p and for large enough τ ,   Similarly, we need to solve the inverse problem to obtain results in Minkowski space.
The results from the ME method are shown in Fig. 6. The error bands are mainly from the average of different default models. The behaviors of d and u quarks are similar. We do not observe a peak around the elastic point ν = 0 which is because at this point the hadronic tensor is the square of the electromagnetic form factor and the form factor for elastic scattering is highly suppressed. Taking a dipole form for the form factor with Q 2 = 12.7 GeV 2 in this case, the hadronic tensor is suppressed by a factor of (1+Q 2 /0.71) −4 ∼ 10 −5 as compared to the charge at Q 2 = 0. We do observe a broad structure shows at about 1 GeV, which should be the combined contribution of nucleon resonances and possibly SIS.
As discussed above, the preferred ν range that can lead us to the parton structure functions is from 2.96 GeV to 3.68 GeV, however, our results show that it is basically zero within error in that region. To check whether this is a resolution issue of ME, we also use the BR method that shows better resolution for discrete structures in our mock data test to handle the same data. The results are shown in Fig. 7. And this time, to show exactly the effect of different default models, we plot the results with different default models separately in log scale. We also check the effect of including or excluding the data points of large τ since they can have large excited-state contaminations. This time, except that the peaks around 1 GeV are much shaper than the ME case, the basic conclusion we learn is the same. No elastic contribution shows at ν = 0 and the only structure is around 1 GeV. When the energy transfer goes above, say, 2 GeV, the reconstructed results approach to those of the default model values which means the data have no constraint in that region.  Compton amplitude T µν (p, q) in [7]) to where the state label 0 is the nucleon state with momentum p + q so that the first term is the elastic scattering which diverges as T and reflects the elastic scattering pole in Eq. (5).
As T → ∞, the n ≥ 1 state contributions are suppressed by 1/T . However, as shown in Fig. 6 when T is finite (0.7 fm in our case), the excited states including nucleon resonances and those in the SIS and DIS regions, all contribute. When Q 2 is large (e.g., 12.7 GeV 2 in this case), the hadronic tensor for the elastic scattering is highly suppressed (by a factor of ∼ 10 −5 ). Whereas, the resonance contribution around 1 − 2 GeV at Q 2 ∼ 9 − 11 GeV 2 is much larger as shown in Fig. 6. To estimate how large a ν is needed for DIS, we can look at W , the total invariant mass of the hadronic final state The global fitting of PDF usually take a cut with W 2 > 10 GeV 2 . When we take Q 2 = 4 GeV 2 , this gives ν > 6.5 GeV. Therefore, taking ν = 0 in Eq. (4) will not yield PDF in the

V. DISCUSSION AND SUMMARY
To explore the reason why there is no contribution for ν 2 GeV in Fig. 6 and Fig. 7 Fig. 8. We see that for either u or d quark, the highest effective mass is around 1 GeV, which means that there is simply no information of higher excitations for this particular case. This should be due to lattice artifacts, since the lattice we are using has finite volume (resulting in discrete momenta and discrete spectrum), finite lattice spacing (an UV cutoff) and unphysical pion mass (unphysical multi-particle states).
To sort out the most important factor, we calculate the effective mass of the ρ meson with different lattice setups (right panel of Fig. 8). The reason we choose to check ρ meson is because the hadronic tensor involves two vector currents inserted between the nucleon states and the correlator of ρ can be treated as two vector currents inserted between the vacuum.
Although the exact value of how high we can reach in the ρ meson case may not have much to do with the hadronic tensor case, the fact that how lattice artifacts affect the effective mass should be relevant. The legend of the figure shows the features of the setups. Each label in the legend has four parts: valence quark type/sea quark type ("O" denotes overlap, "D" for domain wall, "C" for clover and "H" for HISQ), spacial size plus a suffix which serves as an identifier, spacial lattice spacing, and sea pion mass. It is easy to see that, for 24I and 48I [29], although the pion masses and volumes are not the same, the highest effective masses are similar, around 3 GeV. For 24J and 16J [31], the spacial lattice spacings are similar to the ones of 24I and 48I (∼ 0.1 fm), but the highest effective mass can be higher than 5 GeV, which is because these two lattices are anisotropic and their temporal lattice spacings are about 0.035 fm. Then, for 48H [32] and the two setups of 32If [29], despite their different fermion actions and volumes, their highest effective masses are all about 5.5 GeV and their lattice spacings are ∼0.06 fm. For the lattice with lattice spacing ∼0.045 fm (64H [32]), the highest effective mass is close to 8 GeV. This test shows that the lattice spacing is the most important factor in order to have the information of higher excitations.
In view of this comparison, the HISQ lattice with lattice spacing ∼0.045 fm can be a better choice to reach ν > 2 GeV.
In this paper, we formulate our approach of calculating the hadronic tensor on the lattice. 24J and 16J are anisotropic lattices with a s /a t = 3.7, so their highest effective masses are higher than those of 24I and 48I that have similar spacial lattice spacings. To increase visibility, some points in the right panel are shifted slightly in the horizontal direction.
We point out that this is the approach that includes the inclusive contribution of all the intermediate states which is crucial to providing information for the neutrino scattering experiments at low energies. It is also promising to calculate the structure function in the DIS region which can be used in the global fittings of parton distribution functions.
However, solving the inverse problem is the most challenging part. We have implemented and tested three algorithms using mock data, showing that the BR method has the best resolution in extracting peak structures while BG and ME are more stable for the flat spectral function. Realistic lattice results are presented for both the elastic case and a case with large momentum transfer. For the elastic case, the reconstructed Minkowski hadronic tensor from the ME method gives precisely the vector charge which shows the feasibility of this approach. For the latter case, the RES and possibly SIS contributions around 1 GeV are observed but no information is obtained for higher excited states with ν > 2 GeV.
A check of the effective masses of ρ meson with different lattice setups indicates that, in order to reach higher energy transfers, using lattices with smaller lattice spacings is essential for the lattice calculation. The HISQ lattice with lattice spacing ∼0.045 fm should be suitable to study the neutrino nucleus scattering at DUNE where the beam energy is between ∼ 1 to ∼ 7 GeV. In the future, working on lattices with lattice spacing of 0.3 fm or smaller would be desirable for studying the parton physics.