Towards the nucleon hadronic tensor from lattice QCD

TL;DR

This paper presents the first lattice QCD calculation of the nucleon hadronic tensor, enabling insights into structure functions and neutrino scattering, using novel algorithms to solve an inverse problem with promising results.

Contribution

It introduces a lattice QCD approach to compute the nucleon hadronic tensor and compares three algorithms for the inverse problem, demonstrating feasibility and initial results.

Findings

Successfully reconstructed the vector charge from the hadronic tensor.

Observed nucleon resonances and inelastic scattering contributions at certain energies.

Identified the need for smaller lattice spacings to access higher energy transfers.

Abstract

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.

Figures

Tables

Table 1. Table 1 : The kinematic setup of nucleon momentum p → → 𝑝 \vec{p} , three-momentum transfer q → → 𝑞 \vec{q} , proton energy E p subscript 𝐸 𝑝 E_{p} , the energy of the lowest intermediate state E n = 0 subscript 𝐸 𝑛 0 E_{n=0} , the modulus of the three-momentum transfer | q → | → 𝑞 |\vec{q}| , the range of energy transfer ν 𝜈 \nu , the range of four-momentum transfer Q 2 superscript 𝑄 2 Q^{2} and the corresponding Bjorken x 𝑥 x .

$\overrightarrow{p}$ ($2\pi /L_s$)	$\overrightarrow{q}$ ($2\pi /L_s$)	$E_p$ (GeV)	$E_{n=0}$ (GeV)	$\left|\overrightarrow{q}\right|$ (GeV)	$\nu$ (GeV)	$Q^2$ (GeV2)	$x$
(0,3,3)	(0,-6,-6)	2.15	2.15	3.57	[2.96, 3.68]	[4, 2]	[0.16, 0.07]

Equations

$$W_{\mu \nu }=\frac{1}{4\pi }\int d^{4}z{e}^{iq\cdot z}⟨p,s\left|\left[J_{\mu }^{†}(z)J_{\nu }(0)\right]\right|p,s⟩.$$

$$W_{\mu \nu }=\left(-g_{\mu \nu }+\frac{q_{\mu }{q}_{\nu }}{q^2}\right)F_1(x,Q^2)+\frac{{\hat{p}}_{\mu }{\hat{p}}_{\nu }}{p\cdot q}{F}_2(x,Q^2)$$

$$W_{\mu \nu }^M=\frac{1}{2}\sum _n\int \prod _i^n\left[\frac{d^3{\vec{p}}_i}{(2\pi {)}^{3}2E_i}\right]⟨p,s|J_{\mu }^{†}(0)|n⟩⟨n|J_{\nu }(0)|p,s⟩(2\pi {)}^3{\delta }^4(q-p_n+p),$$

$$W_{\mu \nu }^{\prime }=\frac{1}{4\pi }\sum _n\int dt{e}^{\left(\nu -(E_n-E_p)\right)t}\int d^3\vec{z}{e}^{i\vec{q}\cdot \vec{z}}⟨p,s|J_{\mu }^{†}(\vec{z})|n⟩⟨n|J_{\nu }(0)|p,s⟩,$$

$$W_{\mu \nu }^{\prime }=\frac{1}{4\pi }\sum _n\frac{e^{\left(\nu -(E_n-E_p)\right)T}-1}{\nu -(E_n-E_p)}\int d^3\vec{z}{e}^{i\vec{q}\cdot \vec{z}}⟨p,s|J_{\mu }^{†}(\vec{z})|n⟩⟨n|J_{\nu }(0)|p,s⟩,$$

$$C_4(t_f,t_2,t_1)=\sum _{{\vec{x}}_f}{e}^{-i\vec{p}\cdot {\vec{x}}_f}\sum _{{\vec{x}}_1{\vec{x}}_2}{e}^{-i\vec{q}\cdot ({\vec{x}}_2-{\vec{x}}_1)}⟨{\chi }_N({\vec{x}}_f,t_f)J_{\mu }^{†}({\vec{x}}_2,t_2)J_{\nu }({\vec{x}}_1,t_1){\bar{\chi }}_N(\vec{0},t_0)⟩,$$

$$C_2(t_f)=\sum _{{\vec{x}}_f}{e}^{-i\vec{p}\cdot {\vec{x}}_f}⟨{\chi }_N({\vec{x}}_f,t_f){\bar{\chi }}_N(\vec{0},t_0)⟩.$$

$$W_{\mu \nu }^E(\vec{p},\vec{q},\tau )=\frac{E_p}{m_p}\frac{\mathrm{Tr}[{\Gamma }_e{C}_4]}{\mathrm{Tr}[{\Gamma }_e{C}_2]}\to \sum _{{\vec{x}}_1{\vec{x}}_2}{e}^{-i\vec{q}\cdot ({\vec{x}}_2-{\vec{x}}_1)}⟨p,s|J_{\mu }^{†}({\vec{x}}_2,t_2)J_{\nu }({\vec{x}}_1,t_1)|p,s⟩,$$

$$W_{\mu \nu }^E=\sum _n{A}_n{e}^{-(E_n-E_p)\tau },$$

$$A_n\equiv \sum _{{\vec{x}}_1{\vec{x}}_2}{e}^{-i\vec{q}\cdot ({\vec{x}}_2-{\vec{x}}_1)}⟨p,s|J_{\mu }^{†}({\vec{x}}_2,0)|n⟩⟨n|J_{\nu }({\vec{x}}_1,0)|p,s⟩$$

$$W_{\mu \nu }^M(\vec{p},\vec{q},\nu )=\frac{1}{i}{\int }_{c-i\infty }^{c+i\infty }d\tau e^{\nu \tau }{W}_{\mu \nu }^E(\vec{p},\vec{q},\tau ).$$

$$W_{\mu \nu }^E(\vec{p},\vec{q},\tau )=\int d\nu W_{\mu \nu }^M(\vec{p},\vec{q},\nu )e^{-\nu \tau }.$$

$$c({\tau }_i)=\int k({\tau }_i,\nu )\omega (\nu )d\nu ,$$

$$c({\tau }_i)=\sum _{j}k({\tau }_i,{\nu }_j)\omega ({\nu }_j)\Delta {\nu }_j;$$

$$\sum _{i}a({\tau }_j,{\nu }_0)k({\tau }_i,\nu )\sim \delta (\nu -{\nu }_0),$$

$$\sum _{i}a({\tau }_i,{\nu }_0)c({\tau }_i)\sim \int \delta (\nu -{\nu }_0)\omega (\nu )d\nu =\omega ({\nu }_0).$$

$$P[\omega |D,\alpha ,m]\propto \frac{1}{Z_S{Z}_L}{e}^{Q(\omega )},$$

$$S=\sum _j\left[\omega ({\nu }_j)-m({\nu }_j)-\omega ({\nu }_j)\log \left(\frac{\omega ({\nu }_j)}{m({\nu }_j)}\right)\right]\Delta {\nu }_j,$$

$$L=\frac{1}{2}\sum _{i,j}\left(c({\tau }_i)-c^{\omega }({\tau }_i)\right)C_{ij}^{-1}\left(c({\tau }_j)-c^{\omega }({\tau }_j)\right)$$

$$P[\omega |D,\alpha ,m]\propto e^{Q^{\prime }(\omega )},$$

$$S^{\prime }=\sum _j\left[1-\frac{\omega ({\nu }_j)}{m({\nu }_j)}+\log \left(\frac{\omega ({\nu }_j)}{m({\nu }_j)}\right)\right]\Delta {\nu }_j,$$

$$P[\omega |D,m]=\frac{P[D|\omega ,I]}{P[D|m]}\int d\alpha P[\alpha |D,m],$$

$$W_{44}^E(\vec{0},\vec{0},\tau )\sim \sum _n⟨p,s|\bar{\psi }{\gamma }_4\psi |n⟩⟨n|\bar{\psi }{\gamma }_4\psi |p,s⟩e^{-(E_n-M_p)\tau }.$$

$$W_{44}^M(q^2,\nu )=\delta (q^2+2m_p\nu )\frac{2m_p}{1-q^2/{4m_p}^2}\left(G_E^2(q^2)-\frac{q^2}{4M_N^2}{G}_M^2(q^2)\right).$$

$$W_{\mu \nu }^{\prime }(T)$$

$$W^2=(q+p{)}^2=m_p^2-Q^2+2m_p\nu .$$

Full text

I Introduction

In scattering processes involving nucleons such as deep inelastic scattering (DIS) and neutrino-nucleon scattering at low energies, the hadronic tensor $W_{\mu\nu}$ is used to characterize the nonperturbative nature of the nucleon structure. It is the imaginary part of the forward virtual Compton scattering amplitude W_{\mu\nu}=\frac{1}{2\pi}{\rm Im}T_{\text{\mu\nu}} and can be expressed as a nucleon matrix element with two current operators inserted

[TABLE]

The hadronic tensor of the vector currents can be further decomposed, according to its Lorentz structure, into structure functions, i.e.,

[TABLE]

for the unpolarized case, where $\hat{p}_{\mu}=p_{\mu}-\frac{p\cdot q}{q^{2}}q_{\mu}$. $p_{\mu}$ and $q_{\mu}$ are the nucleon 4-momentum and momentum transfer, respectively. $x$ is the Bjorken $x=\frac{Q^{2}}{2p\cdot 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}=\sum_{a}c_{i}^{a}\otimes f_{a}$, where the convolution kernel $c_{i}^{a}$ 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 (Liu and Dong, 1994; Liu, 2000), it has only recently become feasible to compute the necessary 4-point correlation functions, thanks to the increases in computing power Liu (2016); Liang et al. (2018). 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 (Ji, 2013; Lin et al., 2018), Compton amplitude (Chambers et al., 2017), Pseudo-PDFs (Radyushkin, 2017; Orginos et al., 2017) and Lattice cross sections (Ma and Qiu, 2018; Sufian et al., 2019). 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 Liu (2016) 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 (Liu and Dong, 1994; Liu et al., 2012) that the Gottfried sum rule violation (i.e., the $\bar{u}$ and $\bar{d}$ 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 (Liang et al., 2019), 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 (Acciarri et al., 2015) 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 Alvarez-Ruso et al. (2018). 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 (Formaggio and Zeller, 2012). 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 neutrino-nucleus 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 Kronfeld et al. (2019)333We thank A. Kronfeld for bringing to our attention the relevance of the hadronic tensor to the neutrino-nucleon scattering..

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 Liu (2016):

[TABLE]

where $q$ is the momentum transfer, $p$ is the nucleon momentum and $p_{n}$ is the momentum of the $n$’th intermediate state. The 4-dimensional Dirac delta function ensures the conservation of 4-momentum and picks out the contribution of a particular momentum transfer. However, it is noted in Liu (2016) that if one carries out the integral in Eq. (1) in the Euclidean case, the Fourier transform in the time direction becomes a Laplace transform

[TABLE]

and after the integration one has

[TABLE]

where $\nu$ is the energy transfer, $E_{p}$ and $E_{n}$ are the energies of the external nucleon and the $n$’th intermediate state, and $T$ is the integration length in the time direction. The factor of $\frac{e^{\left(\nu-(E_{n}-E_{p})\right)T}-1}{\nu-(E_{n}-E_{p})}$ is problematic since it does not converge if $\nu-(E_{n}-E_{p})>0$. This happens when $\nu$ is greater that the energy gap between the nucleon and the intermediate state, such as $\Delta$, Roper, … . Besides, even if the numerator converges, $\frac{1}{\nu-(E_{n}-E_{p})}$ is not a good approximation to the delta function and one cannot pick out the clean contribution for a specific $\nu$ as nearby states will mix.

Instead, we construct the following 4-point correlation function with only a 3-dimensional Fourier transform

[TABLE]

and the normal nucleon 2-point function as

[TABLE]

Then, the Euclidean hadronic tensor ${W}^{E}_{\mu\nu}(\vec{p},\vec{q},\tau)$ is defined by the ratio of the 4-point function to the 2-point function Liu and Dong (1994); Liu (2000); Liu et al. (2012); Liu (2016); Liang et al. (2018)

[TABLE]

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

[TABLE]

where

[TABLE]

and $\tau=t_{2}-t_{1}$.

The computation of 4-point functions consumes the most computer resources in our calculation. Doing the Wick contraction for the 4-point function leads to several topologically distinct diagrams in the Euclidean path-integral formulation (shown in Fig. 1). We have not specified the flavor of the quark lines in the figure. In practice where the flavor is taken into account, the contractions can be more complicated according to the types of the currents (neutral or charged). For the calculation with small momentum transfers, all the diagrams contribute and one needs to combine them to have physical results. However, for the case with large momentum and energy transfer as in the DIS, Figs. 1(a), 1(b) and 1(c) are dominated by the leading twist while Figs. 1(d), 1(e) and 1(f) are suppressed since they involve only high twists. The respective leading twist parton degrees of freedom are classified by the first three diagrams, namely the valence and connected sea (CS) partons $q^{\rm v+cs}$ (Fig. 1(a)), the CS anti-partons $\bar{q}^{\rm cs}$ (Fig. 1(b)) and the disconnected sea (DS) partons and anti-partons $q^{\rm ds}+\bar{q}^{\rm ds}$ (Fig. 1(e)) Liu and Dong (1994); Liu (2000). 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 Liu and Dong (1994); Liu (2000); Liu et al. (2012).

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:

[TABLE]

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}_{\mu\nu}$ Liu (2016); Liang et al. (2018),

[TABLE]

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

[TABLE]

where $c(\tau_{i})$ denotes discrete lattice data with finite number of points (usually $O(10)$), $k(\tau_{i},\nu)$ is the integral kernel that is a function of both $\tau_{i}$ and $\nu$, and $\omega(\nu)$ is the target function which is usually continuous with respect to $\nu$. 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 $\omega(\nu)$:

[TABLE]

however, the number of $\nu_{j}$ one needs to reproduce the structures of $\omega(\nu)$ 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.

III.1 Backus-Gilbert method

The Backus-Gilbert (BG) method Backus and Gilbert (1968, 1970); Hansen et al. (2017) utilizes the fact that the kernel functions can be linearly combined to approximate the Delta function, if they span a complete function basis

[TABLE]

where $a(\tau_{i},\nu_{0})$ are the coefficients for the $i$’th kernel function at a specific point $\nu_{0}$, which can be calculated by assuming a criterion of “deltaness” and solving the linear equations. Having $a(\tau_{i},\nu_{0})$, the value of the target function at $\nu_{0}$ is

[TABLE]

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 Hansen et al. (2017). 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 Kim et al. (2018), so careful checks are alway needed.

III.2 Maximum Entropy method

The Maximum Entropy (ME) method Rietsch (1977); Asakawa et al. (2001) makes use of the Bayesian probability with prior information about the target function to find the most probable solution

[TABLE]

where $P[\omega|D,\alpha,m]$ denotes the probability that $\omega$ is the solution given lattice data $D$, prior information $m$ and a hyper parameter $\alpha$. $Q=\alpha S-L$ is a combination of the Shannon entropy

[TABLE]

which entails the constraint from the prior information, and the likelihood function

[TABLE]

with $c^{\omega}(\tau_{i})$ being the correlator reconstructed using Eq. (14) and $C$ the covariance matrix, which embodies the constraint from the data. $m(\nu_{j})$ is the default model (the prior information we plug in) and $\alpha$ is the weight that balances the two constraints. If $\alpha$ is zero, ME reduces to the normal $\chi^{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 $\alpha$ and the results of different $\alpha$’s are averaged in a range of $\alpha$ based on certain assumptions Asakawa et al. (2001). Practically, the parameter space for finding the maximum probability $P[\omega|D,\alpha,m]$, i.e. the maximum value of $Q$, is reduced to a smaller one instead of the whole $\nu$ space by employing singular value decomposition Bryan (1990); Asakawa et al. (2001), which makes the maximum search easier while the resolution may be affected. Improved ME with an extended search space Rothkopf (2013) is also proposed and we will check whether it produces better results in our future study.

III.3 Bayesian Reconstruction

Bayesian Reconstruction (BR) Burnier and Rothkopf (2013) is an improved Bayesian method. The Bayesian probability is

[TABLE]

where $Q^{\prime}=\alpha S^{\prime}-L-\gamma(L-N_{\tau})^{2}$ and

[TABLE]

which is an alternative way to encode the constraint from the prior. $\gamma$ is a numerically large number such that the term $\gamma(L-N_{\tau})^{2}$ helps to prevent over-fitting. The hyper parameter $\alpha$ here is integrated over as

[TABLE]

where $P[\alpha|D,m]$ is the probability of $\alpha$ in the present of $D$ and $m$. The search space of the maximum is also the whole parameter space which enhances the ability of finding the real maximum while the price to pay is the need of high-precision architecture (e.g., 512-bit floating point numbers). For both Bayesian-type methods, ME and BR, the choice of the default model is in principle arbitrary, but a reliable reconstruction should not depend on the default model. In all the following calculations of the paper, the default model is chosen to be a constant.

To test the three methods, mock data are generated. The first set of mock data is a 2-point function with three states of mass 1.0, 1.5 and 1.8 GeV respectively and the spectral weights are all unity. The lattice spacing is tuned to $0.1$ fm and the number of time slices available is 20. Noises are added by assuming normal distributions around the central values and the signal-to-noise ratio is set to be 100. This is to check the resolution for peak structures. The results of the reconstruction can be found in Fig. 2. The left panel shows the reconstructed spectral functions and the right panel shows the comparison between the original input data and the 2-point functions computed using the spectral functions. For the spectral functions we see that the result of BG shows basically no peak structures, reflecting its poor resolution for this setup. The result of ME is much better as two peaks around 1 and 2 GeV are clearly seen. Considering the broad widths of the peaks, it is consistent with the input data. BR gives the best reconstruction in this case as two sharp peaks appear at $\sim 1$ and $\sim 1.6$ GeV. Although the second and the third states are not separated, the resolution of BR is much better. From the right panel we see that the regenerated 2-point functions from the spectral functions of ME and BR are well consistent with the input data, but for the BG case, differences occur at large $t$. This exhibits the fact that, as pointed out before, BG does not guarantee that the reconstructed result reproduces the input data.

The second set of mock data is a 2-point function with two isolated states of mass 1.0, 1.5 and a dense spectrum (simulating the continuous spectrum) from 2 GeV to 15 GeV. This is like the energy dependence of the neutrino-nucleon scattering cross section. The lattice spacing is set to $0.02$ fm in this case and the number of time slices is 100. The signal-to-noise ratio is set to be 100 too. This is to check the case of a flat spectral function up to higher energies. The results of the this test are shown in Fig. 3. Similarly, the left panel shows the reconstructed spectral functions and the right panel shows the comparison between the original input data and the 2-point functions recomputed using the spectral functions. From the left panel we see that both ME and BR reconstruct a peak at $\sim 1$ GeV (the BR one is much sharper) while BG shows no peak in that energy region. This agrees with what we observed in the previous test. However, in the region 2 to 8 GeV, BG presents more consistent results with the input (the red dashed line) while both ME and BR show unphysical oscillations. These unphysical oscillations are called “ringing” and are artifacts of the reconstruction Kim et al. (2018). 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 Fischer et al. (2018); Kim et al. (2018) 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 Karpie et al. (2019). 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.

IV Preliminary results

IV.1 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 Blum et al. (2016) with clover fermions as valence quarks. The configuration is preprocessed by HYP smearing and the tadpole improved clover coefficient $C_{\rm sw}=1.033$ is used to generate the clover term of the clover action for the valence quarks. The pion mass is tuned to be close to the unitary point $\sim 371$ MeV. The lattice spacing is about 0.06 fm which we expect to be fine enough such that the inverse algorithms can give reasonably reliable results.

For this case, we choose $\mu=\nu=4$ and $\vec{p}=\vec{q}=0$ for ${W}^{E}_{\mu\nu}$. So Eq. (8) becomes

[TABLE]

For simplicity, the two currents are both inserted on the $d$ quark line so only Fig. 1(a) contributes. And for $\tau\gg 0$ only the ground state survives, so ${W}^{E}_{4,4}=g_{V}^{2}=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. 1(a)). 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 $\tau=t_{2}-t_{1}$ is from $1$ to $9$. The result of ${W}^{E}_{44}(\tau)$ is plotted in Fig. 4(a). It shows that within errors, ${W}^{E}_{4,4}(\tau)$ is indeed a constant of value $1$. The drop at $\tau=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$.

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 Cheng and Li (1984) is

[TABLE]

In our case $W^{M}_{44}(\nu)\overset{q^{2}=0}{=}\delta({\nu})G^{2}_{E}(0)=\delta(\nu)$. This is easy to understand since the spectral function should be a Dirac delta function at $\delta(\nu-\nu_{0})$ when the Euclidean correlator is a single exponential $\sim e^{-\nu_{0}t}$ and here the constant is a special case of an exponential with $\nu_{0}=0$.

The converted results of $W^{M}_{44}(\nu)$ using ME and BG are plotted in Fig. 4(b) and Fig. 4(c), respectively. They both give a peak around $\nu=0$ and in this sense the results are consistent with the theoretical prediction of $\delta(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 $\nu=0$, which is because BG has difficulties in resolving the results of negative $\nu$. An important check is that the area under the peaks should be $g_{V}^{2}=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.

IV.2 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. $\sim 0.1$). For this calculation, we use an anisotropic clover lattice Lin et al. (2009) with $a_{t}\sim 0.035$ fm. The pion mass is about 380 MeV and the momentum unit is $\frac{2\pi}{L_{s}}\sim 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. 1. We choose $\vec{p}=(0,3,3)$ and $\vec{q}=(0,-6,-6)$ in lattice unit and $\mu=\nu=1$, such that only the $F_{1}$ structure function survives; thus $W^{M}_{11}=F_{1}(x,Q^{2})$. Since the energy transfer $\nu$ is not fixed by the lattice 3-dimensional Fourier transform, we can choose a range of $\nu\in[2.96,3.68]$ GeV such that the corresponding $Q^{2}$ is in a range of 2 to 4 GeV2. 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 $\vec{p}+\vec{q}=-\vec{p}$, therefore the energy of the lowest intermediate state $E_{n=0}=E_{p}$ and for large enough $\tau$, ${W}^{E}_{11}\propto e^{-(E_{n=0}-E_{p})\tau}$ is a constant.

The results of the Euclidean hadronic tensor ${W}^{E}_{11}(\tau)$ for both $d$ and $u$ quarks are shown in Fig. 5. In this calculation, we set $t_{0}=0$, $t_{1}=8$ and $t_{f}=28$ in lattice unit and still we focus on the diagram in Fig. 1(a) only for the moment. The figures show that for small $\tau$, the Euclidean hadronic tensor ${W}^{E}_{11}(\tau)$ decays exponentially since the energy of the intermediate states $E_{n}$ are larger than $E_{p}$. For larger $\tau$, it gets flatter for the $d$ quark case (left panel) which is consistent with what we expect, i.e. $E_{n=0}=E_{p}$. The tail of the Euclidean hadronic tensor for the $u$ quark case (right panel) goes up after $\tau\sim 15$, which we believe is due to the contamination of the sink nucleon excited states.

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 $\nu=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\leavevmode\nobreak\ \rm{GeV}^{2}$ in this case, the hadronic tensor is suppressed by a factor of $(1+Q^{2}/0.71)^{-4}\sim 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 $\nu$ 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 $\tau$ 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 $\nu=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.

From these two examples, we see that the elastic and the resonance contributions, where the energy transfer is not too large, can be extracted for the current setup. However, it seems that there are no contribution in the DIS region. We will discuss the possible reason and solution in the next section.

It is suggested by A.J. Chambers et al. Chambers et al. (2017) that the structure functions can be extracted through the time-ordered current-current correlator on the lattice through the operator product expansion of the forward Compton amplitude. This is carried out by setting $\nu=0$ in the integral in Eq. (4) to avoid the divergence when $\nu-(E_{n}-E_{p})>0$ as discussed in Sec. II after Eq. (5). However, this leads $W_{\mu\nu}^{\prime}(p,q,T)$ in Eq. (5) (considered to be the Compton amplitude $T_{\mu\nu}(p,q)$ in  Chambers et al. (2017)) to

[TABLE]

where the state label [math] is the nucleon state with momentum $\vec{p}+\vec{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\rightarrow\infty$, the $n\geq 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 $\rm{GeV}^{2}$ in this case), the hadronic tensor for the elastic scattering is highly suppressed (by a factor of $\sim 10^{-5}$). Whereas, the resonance contribution around $1-2$ GeV at $Q^{2}\sim 9-11\,\rm{GeV}^{2}$ is much larger as shown in Fig. 6. To estimate how large a $\nu$ is needed for DIS, we can look at $W$, the total invariant mass of the hadronic final state

[TABLE]

The global fitting of PDF usually take a cut with $W^{2}>10\,\rm{GeV}^{2}$. When we take $Q^{2}=4\,\rm{GeV}^{2}$, this gives $\nu>6.5$ GeV. Therefore, taking $\nu=0$ in Eq. (4) will not yield PDF in the DIS region which needs both $Q^{2}$ and $\nu$ to be large.

V Discussion and Summary

To explore the reason why there is no contribution for $\nu\gtrsim 2$ GeV in Fig. 6 and Fig. 7, we calculate the effective mass from 2-point functions which is a quick way to check the highest energy of intermediate states that our Euclidean hadronic tensor contains. The results are plotted in the left panel of 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 $\rho$ meson with different lattice setups (right panel of Fig. 8). The reason we choose to check $\rho$ meson is because the hadronic tensor involves two vector currents inserted between the nucleon states and the correlator of $\rho$ can be treated as two vector currents inserted between the vacuum. Although the exact value of how high we can reach in the $\rho$ 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 Blum et al. (2016), although the pion masses and volumes are not the same, the highest effective masses are similar, around 3 GeV. For 24J and 16J Lin et al. (2009), the spacial lattice spacings are similar to the ones of 24I and 48I ($\sim$ 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 Bazavov et al. (2010) and the two setups of 32If Blum et al. (2016), despite their different fermion actions and volumes, their highest effective masses are all about 5.5 GeV and their lattice spacings are $\sim$0.06 fm. For the lattice with lattice spacing $\sim$0.045 fm (64H Bazavov et al. (2010)), 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 $\sim$0.045 fm can be a better choice to reach $\nu>2$ GeV.

In this paper, we formulate our approach of calculating the hadronic tensor on the lattice. 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 $\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 for the lattice calculation. The HISQ lattice with lattice spacing $\sim$0.045 fm should be suitable to study the neutrino nucleus scattering at DUNE where the beam energy is between $\sim 1$ to $\sim 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.

Acknowledgements.

KFL thanks X. Feng, A. Kronfeld, J. C. Peng, J. Qiu and Y. Hatta for illuminating and useful discussions. We also thank the RBC and UKQCD Collaborations for providing their DWF gauge configurations. This work is supported in part by the U.S. DOE Grant No. DE-SC0013065 and DOE Grant No. DE-AC05-06OR23177 which is within the framework of the TMD Topical Collaboration. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This work used Stampede time under the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1053575. We also thank the National Energy Research Scientific Computing Center (NERSC) for providing HPC resources that have contributed to the research results reported within this paper. We acknowledge the facilities of the USQCD Collaboration used for this research in part, which are funded by the Office of Science of the U.S. Department of Energy. Y.Y. is also supported by the CAS Pioneer Hundred Talents Program.