Towards fully non-perturbative computation of inelastic $\ell N$ scattering cross sections from lattice QCD

We propose a fully non-perturbative method to compute inelastic lepton-nucleon ($\ell N$) scattering cross sections using lattice QCD. The method is applicable even at low energies, such as the energy region relevant for the recent and future neutrino-nucleon ($\nu N$) scattering experiments, for which perturbative analysis is invalidated. The basic building block is the forward Compton-scattering amplitude, or the hadronic tensor, computed on a Euclidean lattice. Total cross section is constructed from the hadronic tensor by multiplying a phase space factor and integrating over the energy and momentum of final hadronic states. The energy integral that induces a sum over all possible final states is performed implicitly by promoting the phase space factor to an operator written in terms of the transfer matrix on the lattice. The formalism is imported from that of the inclusive semileptonic B meson decay [P. Gambino, S. Hashimoto, arXiv:2005.13730]. It can be generalized to compute the $\ell N$ scattering cross sections and their moments, as well as the virtual correction to the nuclear $\beta$-decay. Necessary quark-line contractions for two current insertions corresponding to the Compton amplitude to be computed on the lattice are summarized.


I. INTRODUCTION
Deep inelastic scattering (DIS) played an important role on the emergence of the parton picture of nucleon and the discovery of the asymptotic freedom, which lead to the fundamental theory of strong interaction, Quantum Chromodynamics (QCD).The key finding was that the structure function W µν may well be described by a sum of parton distributions, and the partons inside nucleon behave as if they are free particles despite the strong force that binds them together.And, one does not have to take account of the details of the individual hadronic final states when calculating the total cross section, which is called the quark-hadron duality.Perturbative analysis of DIS is based on these observations.Theoretically, the DIS process is factorized into perturbative scattering amplitudes of elementary partons, i.e. quarks and gluons, and non-perturbative functions called parton distribution functions (PDFs), which have to be determined by fitting the experimental data.
(For a review of the early days of the development, see, for instance, [2].)However, the separation of perturbative and non-perturbative contributions, so-called the factorization [3], is not obvious, especially for higher twist contributions which become relevant when one tries to calculate beyond the leading order of the operator product expansion (OPE).
There is even an intrinsic difficulty in the scale separation using OPE due to the renormalon ambiguity (see, for instance, [4,5]).Therefore, the theoretical analysis of the lepton-nucleon ( N ) scattering has been limited in the high-energy regime, where the effects of higher twist operators are negligible.
The other limit, i.e. the low-energy limit, is given by the elastic scattering of nucleon, whose form factors can be calculated using lattice QCD, and the theoretical computation is based on a solid ground.In the intermediate energy region, where extra pion(s) can be generated but the process is still highly non-perturbative, the details of the final state hadrons also become relevant, and one finds resonances and other structures in the differential cross sections, i.e. the duality is violated.The theoretical analysis of such processes is far more complicated and quantitative computation based on QCD has remained impractical.
In this work we construct a formalism to compute total inelastic cross-section of the N scattering fully non-perturbatively using lattice QCD.The formalism does not involve the parton distribution functions.Rather, we directly compute the (weighted) integral of the hadronic structure function, that defines the total cross section.Since the method does not rely on the separation of perturbative and non-perturbative contributions, there is no fundamental limitation on the range of Q 2 that can be treated.Namely, the analysis is useful even when the momentum transfer Q 2 is small, where the standard perturbative analysis is not applicable.This is an advantage for the analysis of the neutrino-nucleon scattering cross section at low energy, i.e. the initial neutrino energy in the range from several hundred MeV to a few GeV, which is the energy range relevant for the neutrino experiments such as T2K, NOνA and DUNE [6].In fact, in this energy region, the cross section is neither dominated by the quasi-elastic scattering, nor described by deep inelastic scattering.The contributions of inelastic processes including a few pions is significant, and there has so far been no established theoretical method that can reliably treat this region.
We propose to use the hadronic tensor, or equivalently the forward Compton-scattering amplitude, computed using lattice QCD.The hadronic tensor is defined as a matrix element of two weak-currents inserted between nucleon states that represent the initial target nucleon.
On the Euclidean lattice, the two currents are placed away from each other in space and imaginary time directions, so that the matrix element obtained after Fourier (or inverse-Laplace) transform is limited in a unphysical kinematical region where the final states never become on-shell.The physically relevant structure function corresponds to the imaginary part of the hadronic tensor on the physical cut, which is hard to obtain with ordinary lattice computation methods.(One has to solve the inverse problem.See below.)Our proposal is to use the method developed to calculate the total semileptonic decay rate of the B meson [1], that enables us to compute the the sum over all possible states between the two weak currents.With an appropriate weight of the energy, so-called the phase space, the sum corresponds to the total scattering cross section after integrating over spatial momentum.
An approximation is introduced for the phase space factor that appears as a weight in the energy integral (or the sum over all possible final states).As we describe later, the approximation is very precise when the weight is a smooth function of the energy.This is not the case in practice because of the sharp upper limit on the energy transfer from the initial lepton to the hadronic system, and the weight factor is actually a discontinuous function of energy.Therefore, the associated error has to be carefully examined.The experience in the study of the B meson semileptonic decay suggests that it is not substantial [1], and we study the size of the potential systematic effect for the case of the N scattering assuming a form of spectrum of the final states.
Recently, some attempts have been proposed to compute the inclusive processes using lattice QCD [7][8][9][10][11][12][13]. In the context of the N scattering, they correspond to the calculation of the forward Compton-scattering amplitude, which is the same as what we treat.The kinematical point accessible on the lattice is apart from the physical cut, and the methods are devised to relate the lattice data to the physical amplitude by solving the inverse problem, which is extremely difficult and there is no satisfactory solution that provides reliable quantitative results for the physical amplitudes.The difference of our proposal is that we do not try to solve the inverse problem, but we advocate to compute only the energy integral of the physical amplitude.In that way, the difficulty of the inverse problem is circumvented, and the physical quantity of interest is accessible.
The method of energy integral with two current insertions may potentially be applied also to the study of the Cottingham formula that relates the electromagnetic contribution to the proton-neutron mass difference to the forward Compton-scattering amplitude [14].
(See also [15][16][17] and references therein.)The formula has the form of an integral of the hadronic tensor in terms of the inserted energy and momentum, which is the same structure as the total cross section, but there is an additional complexity due to the ultraviolet divergence and some dedicated analysis would be necessary.It is also related to the two-photon exchange diagram in the N scattering, which is relevant to the precise computation of the electromagnetic radius of proton [18].
Another interesting application is the γW exchange contribution to the nuclear β decay.
At O(α), the nucleon γW box diagram gives rise to a nuclear-structure independent correction to the super-allowed nuclear β decays.Its hadronic uncertainty limits the accuracy of the determination of the Cabibbo-Kobayashi-Maskawa matrix element |V ud |, for which recent phenomenological estimates suggest a tension with the CKM unitarity [19,20].In the γW box diagram involving a lepton and a nucleon in the initial/final states connected by a photon propagator and a W boson propagator, the hadronic states between the weak and electromagnetic currents can be any excited states (with a corresponding quantum number), and their contributions have to be taken into account.The integral over the internal momentum resembles that of the total N cross section, and the method developed in this work is applicable.
The lattice computation of the necessary four-point functions including two current insertions is a major challenge especially when flavor-changing currents are involved.The calculation of the forward Compton-scattering amplitude has been performed so far using the Feynman-Hellmann technique [12,13].In this work we figure out all the necessary quark-line contractions including the cases of the flavor-diagonal, flavor-changing, as well as those for the β-decay for future computations This paper is organized as follows.In Sec.II we outline the method to perform the integral over all possible intermediate states with an appropriate weight.The method is essentially the same as those in [21] for current two-point functions and in [1] for B meson inclusive semileptonic decays.The kinematics of the N scattering is summarized in Sec.III and the master formula of the total cross section is given.An application to the β decay is discussed in Sec.IV.Explicit formula and some examples of the energy integral are then described in Sec.V. Potential errors due to approximation are also investigated.We provide some details of the quark-line contractions, which are necessary to compute the hadronic tensor.They are especially complicated when the charged current, which involves the change of flavors, as described in Sec.VI.Further discussions and future prospect are given in Sec.VII, and our conclusions are in Sec.VIII.

II. LATTICE CORRELATORS AND SPECTRAL FUNCTIONS
We first outline the basic idea of the method we are proposing.As in the standard analysis, the inelastic scattering cross section can be written in terms of a product of the leptonic tensor and hadronic tensor.Using the spectral decomposition, the hadronic tensor can be viewed as a spectral function; the total cross section is its integral over final-state energy with a weight factor determined by the leptonic tensor.Therefore, once the spectral function is extracted from the lattice data, the cross section can be obtained as emphasized in [7].In practice, the extraction of the spectral function needs a solution of the inverse problem, which is a well-known example of ill-posed problems; there have been no practical methods developed so far that allow sufficiently accurate quantitative estimates (see, for instance, [7,8,22]).The problem can be overcome by combining the energy integral with the computation of the forward Compton-scattering matrix element [1,21], as outlined below.The idea was developed from an analysis to relate the different kinematical regions of the Compton amplitude using analytic continuation [23].
Let us consider a matrix element of a nucleon with two current insertions For the moment, we ignore the Lorentz and flavor indices of the current J as well as the spin of the nucleon state |N ; a more concrete definition will be given in the following sections.The matrix element of the form (1) can be computed on the lattice from four-point correlation functions including the source operators to create and annihilate the external state |N .On the second line, we introduce a Fourier transform of the current J(q) ≡ x e −iq•x J(x), and V is the spatial volume of the lattice.We assume that the initial state |N is at rest, so that (1) describes the process where a momentum q is injected at t = 0 and taken out at t. (An extension to the case of non-zero initial momentum of nucleon is straightforward.)The time separation between the two currents is imaginary since the calculation is performed on the Euclidean lattice.The time evolution is then described by a transfer matrix e − Ĥ with Ĥ the Hamiltonian of the system.Here we use the lattice unit, i.e. the lattice spacing is a = 1.In the analysis of the lattice data, we do not need the explicit form of the lattice Hamiltonian Ĥ; it is introduced to remind us that the evolution of individual intermediate eigenstates with energy ω is given by e −ωt .The eigenvalues of z ≡ e − Ĥ are limited in the range [0, 1].
The correlator (2) can be formally decomposed into the contributions of individual energy eigenstates, with the spectral function where the sum runs over all possible states X(q) with a specified momentum q.The δfunction in (4) picks the states of a certain energy ω among all possible states with energy Since the spectral function ρ(ω; q) describes the transition rate of the initial state |N to the states with a certain energy ω and momentum q, the total cross section can be written as an integral of the spectral function with an appropriate weight function K(ω; q): The weight function is determined by the details of the process of interest.For the inelastic N scattering, the spectral function corresponds to the hadronic tensor and the weight function K(ω; q) represents the phase space of the scattering specified by the leptonic tensor.
Since the energy integral in (5) extends up to infinity, the kinematical upper limit for the energy ω is also encoded in the weight function by a step function.The momentum q is also integrated over to obtain the total rate.
One approach to obtain the total rate Γ from ( 5) would be to first extract the spectral function ρ(ω; q) from the lattice correlator by solving the inverse problem (3) and then to use it in (5).The problem is, however, that the inverse problem is extremely difficult in practice since the lattice data C(t; q) are known only at limited values of t with nonnegligible statistical noise.In the context of the nucleon structure, several methods have been proposed to solve the inverse problem, e.g. a reconstruction through moments [12,13], the Backus-Gilbert method [7].More extensive tests of various methods in the market, including the Maximum Entropy Method, Bayesian reconstruction, and neural network methods are found in [11,22].Unfortunately, none of them allows fully quantitative computation of the spectral function ρ(ω; q) as a function of ω.This is only natural because of the complicated structure of the spectrum including resonances and scattering states with interactions.To circumvent the problem, smeared spectral function has been considered.It is defined as ρ(ω) = dω ∆(ω, ω )ρ(ω ) with a certain smearing kernel ∆(ω, ω ) that typically has a peak at ω ≈ ω and rapidly decreases for larger |ω − ω |.The smeared spectrum ρ(ω) is easier to reconstruct with limited numerical data [8,21] when the smearing width is sufficiently large.The computation of the integral (5) from ρ(ω) would then become another non-trivial problem as it needs to take the limit of vanishing smearing width.
A practical method to actually compute the smeared spectral function was proposed in [21], and it has been applied for the inclusive decay rate of B meson [1].The key idea was to identify the weight function K(ω; q) in ( 5) as a smearing kernel, and then to perform the ω-integral using the correlator computed on the Euclidean lattice.The ω-integral can be carried out formally using the relation dω K(ω; q)ρ(ω; q) ∝ dω K(ω; q) = N | J(−q)K( Ĥ; q) J(q)|N (0) .
On the first line, the definition of the spectral function (4) and are inserted.The sum over all possible states X(q) |X(q) X(q)| is performed with the δ-function δ(ω − E X(q) ), replacing the energy of individual states E X(q) by the Hamiltonian Ĥ.
The intermediate form (7) can be viewed as a smeared spectral function.In fact, a spectral function N | J(−q)δ( Ĥ − ω) J(q)|N corresponding to the N scattering process is integrated over the energy with a smearing factor K(ω, q).This smearing kernel does not have a peak structure around ω, but the mathematical form is equivalent.On the last line, the ω-integral is carried out by introducing an operator K( Ĥ; q).The remaining task is then to find an expression of K( Ĥ; q) that can be practically implemented in the lattice calculation.
The use of the smeared spectral function has been extended towards a different direction, i.e. to compute the scattering amplitude [9,10], where the smearing kernel is identified as a factor that appears in the LSZ reduction formula, so that the necessary scattering amplitude is directly obtained.The expression (8) should also be applicable in such a case.
Comparing ( 8) with (2), we notice that the integral can be evaluated using the lattice correlators if the operator K( Ĥ; q) is approximated by a polynomial of the form because the matrix element of the right-hand side is nothing but C(t; q)'s up to a normalization factor.There may be various ways to construct this approximation, and the method introduced in [21] uses the Chebyshev approximation, which is described in the following.
The Chebyshev approximation is not an expansion in terms of some small parameters.
Rather, it attempts to approximate the whole function consider the kernel as a function of z, which can take values between 0 and 1, and expand K(− ln z, q) using an orthonormal set of functions {T * j (z)} with j = 0, 1, 2, • • • , N .Here, T * j (x) represents the shifted Chebyshev polynomials, which are related to the standard Chebyshev polynomials of the first kind T j (x) as T * j (x) ≡ T j (2x−1).The shifted Chebyshev polynomials T * j (x)'s are defined in 0 ≤ x ≤ 1.The polynomial order N controls the precision of the approximation.
The Chebyshev polynomials are obtained by a recurrence relation: so that the shifted ones are derived as The first few of the shifted Chebyshev polynomials are then , and so on.They are shown in Fig. 1 as a function of x (or z) as well as ω through z = e −ω .
The Chebyshev polynomials are designed such that they evenly oscillate between −1 and +1 with j the number of nodes.The Chebyshev approximation of a function f (x) defined in [0, 1] has the form f (x) , and each order represents (a sort of) frequency component of the function.As the higher order terms are added, the approximation can reproduce finer details of the original function.
The Chebyshev approximation of the matrix element of the operator K( Ĥ; q) and thus of the matrix element of ẑ ≡ e − Ĥ can be obtained as follows.We define a state |ψ(q) as |ψ(q) ≡ e − Ĥt 0 J(q)|N with t 0 a small time separation introduced to avoid any divergence when evaluating ψ(q)|ψ(q) .The approximation may then be written as which is based on the corresponding approximation formula with z = e −ω .
The approximation evaluated on the state |ψ(q) , as in (12), may be obtained from that for a c-number (13) by considering a decomposition into energy eigenstates |i of energy ω i , i.e. |ψ = i a i |i .Then, the matrix element of T * j (ẑ) may be written as a sum i |a i | 2 T * j (e −ω i ), while the matrix element of the kernel operator is also given as , each term of which may be approximated using (13).
Each term of the right-hand side of (12), the matrix element of the Chebyshev polynomials T * j (ẑ) may be constructed from those of ẑ = e − Ĥ using which is immediately obtained from the lattice data (2).Then, the right-hand side of ( 12) is nothing but a linear combination of ( 14) with different t's.
The coefficients c * j (q) in ( 12) are obtained by an integral according to the general formula of the Chebyshev approximation.Since K(ω; q) = K(− ln z; q) is a known function, the coefficients can be obtained easily using numerical integration.
The Chebyshev approximation provides the best possible approximation at a given order of the polynomials of the form (12) and thus of any polynomials of that order.It is the best in the sense that the maximum deviation from the true function is minimal in the range 0 ≤ z (= e −ω ) ≤ 1, which covers all positive energy eigenvalues of the final states.
For smooth kernel functions K(ω; q), the coefficients c * j (q) rapidly decrease (often exponentially) for larger j, so that the contributions from higher order terms are suppressed, since the Chebyshev polynomials are bounded: |T * j (x)| ≤ 1.When the kernel function K(ω; q) is non-smooth or even discontinuous, the approximation becomes highly oscillatory near the discontinuity and higher order terms become necessary to suppress them.Examples relevant for the computation of the inelastic scattering cross sections are discussed in Sec.V.
One may also consider to use the Chebyshev approximation to obtain more information of the spectral function ρ(ω; q).In fact, it is possible to introduce a kernel function that has a finite value only in a small bin of ω and vanishes otherwise.With such a filtering function, one can stochastically count the number of energy eigenvalues in that bin, if the Chebyshev approximation works.Unfortunately, the filtering function is highly discontinuous and thus needs higher order Chebyshev polynomial terms to achieve good approximation, which is not practical for this particular application.The Chebyshev eigenvalue filtering technique has been used in the context of lattice QCD computation for the calculation of the Dirac operator eigenvalue spectram, through which one can extract the chiral condensate of QCD [24] as well as the spectral function in the full energy range [25].
Our proposal is to combine the operator representation (8) of the ω-integral with the Chebyshev polynomials in order to write it using the correlators computed on the lattice.
The remaining integral over q, see ( 5), has to be carried out with the lattice data obtained at several values of q.This program has been demonstrated for the B meson inclusive semileptonic decays in [1].The formulation for the N scattering is described in the following sections.
In (5) the integral over the energy ω of the hadronic final state corresponds to the sum over all possible final states with a given spatial momentum q 2 .Many of them are multiparticle states such as N π, N ππ, etc., which have continuous spectra in the infinite volume limit.On the lattice of finite spatial extent, they are discretized to satisfy the periodic boundary condition, so that the ω-integral is actually a sum over various allowed states.
The energy of each state receives power corrections of the form 1/V [26] and the limit of V → ∞ has to be taken.We expect that the power-like finite volume effect is marginalized by the integrals over ω and then q 2 , because the finite-volume correction should be most significant for the low-energy and low-momentum states while the total cross section receives more contributions from higher energy regions due to the phase-space enhancement.The problem remains severe when the upper limit of the energy integral is relatively low so that only a limited number of states are kinematically allowed due to the finite volume effect.
Studies of individual states, e.g.N π states, in the finite volume would be more useful in such cases.

III. νN SCATTERING: KINEMATICS
In this section, we summarize the kinematics of the inelastic N scattering partly to establish our notations and to identify the phase space factor that plays the role of the weight function of the energy integral.We are particularly interested in the νN scattering, which is relevant to the recent and future neutrino experiments, but the formulation can also be applied for electromagnetic scattering of electron (or muon) with a slight modification.
The diagram for the νN scattering is shown in Fig. 2. We assign the energy-momentum as p µ = (E, p), p µ = (E , p ) for the incoming (ν) and outgoing ( ) leptons, and P µ = (M N , 0), P µ X = (ω, P X ) for the target nucleon (N ) and outgoing hadronic system (X), respectively.The momentum transfer is then q µ = (p − p ) µ = (E − E , p − p ).The rest frame of the target nucleon is assumed for simplicity, and M N denotes the nucleon mass.The weak current is denoted as For a lepton scattering off unpolarized nucleon through a W -boson exchange, the differential cross section is given by where dΩ represents a solid angle of the final lepton measured with respect to the direction of the incoming lepton.The leptonic tensor L µν is for the weak current.The last term with the totally anti-symmetric tensor µναβ arises from the cross term between the vector and axial-vector currents.For the electromagnetic interaction the W -boson propagator 1/(q 2 − M 2 W ) in (16) has to be replaced by the photon propagator 1/q 2 .Also, the weak coupling g 2 /2 needs to be replaced by the electric charge e 2 , and the last term in the leptonic tensor (17) has to be omitted.
The hadronic tensor contains the contribution from various hadronic states: Here, |X(P X ) represents arbitrary hadronic final states with total spatial momentum P X .
The states are normalized such that d 3 P X X(P X ) |X(P X ) X(P X )| = 1.Since the nucleon is unpolarized, its polarizations (pol.) are averaged.The mass dimension of W µν is −1.
The total cross section may be obtained by integrating out the differential cross section (16).We choose ω and q 2 = (p − p ) 2 as the kinematical variables as they are convenient for the lattice calculation.Using dE dΩ = (π/EE )dωdq 2 , which assumes an integral over a cylindrical angle, we obtain the total cross section In terms of these variables, ω and q 2 , the standard kinematical variables are written as and the Bjorken scaling variables are The lower and upper limits of the ω-integral in (19) correspond to the kinematical limits of x = 1 and x = 0, respectively.
Figure 3 shows the integral paths of ω and q 2 in the plane of x and Q 2 .For a fixed |q|, the ω integral forms a trajectory shown in the plot from x = 1 down to x = 0. Since the upper limit of |q| is fixed by the initial lepton energy E, the integral region is the range below a line of |q| = E.
Here we write down the explicit forms of the leptonic tensor L µν /2.The metric is chosen as g µν = diag(1, −1, −1, −1).We set the direction of the initial neutrino on the z-axis, i.e.
p µ = (E, 0, 0, E).In the following, we label the z-direction by µ = k.Other two spatial directions are denoted as i and j, which are either x-or y-axis .Each component of L µν is written in terms of ω and q 2 as where q 0 = ω − M N and q k = q 0 + (q 2 − q 2 0 )/2E.Here, L 0i will be combined with W 0i , which supplies another factor of q i , so that the sum L 0i W 0i + L 0j W 0j is proportional to k , which is written as a function of ω and q 2 .The νN scattering cross section is thus obtained from the master formula (19).The dependence on the initial neutrino energy appears only through the upper limit of the q 2 integral.The neutrino energy in the recent neutrino experiments are typically in the range of several hundred MeV to a few GeV.(For a review, see [6].)When the initial neutrino energy cannot be controlled event by event, an weighted average of σ(E) with respect to the incoming neutrino energy distribution is actually observed.Such an average can be easily obtained once the total cross section is calculated as a function of the initial neutrino energy E.

IV. W γ EXCHANGE CONTRIBUTION TO β-DECAY
Another application of our formulation is about the higher order corrections to the neutron β decay.It could also be considered as a nuclear-structure independent correction to the nuclear β decay, which is relevant for the precise determination of |V ud |.
At the leading order the β decay occurs through a virtual W exchange.At the next order in α, there may be another exchange of a virtual photon between the final state proton and electron, that makes a box diagram.An estimate of such diagram can be obtained easily using the electromagnetic form factor of proton if the intermediate state can be assumed to be a ground-state proton, but actually the contribution from excited states has to be taken into account, and it is a source of significant uncertainty.In our formulation, the integral over the inner-loop momentum can be carried out with the effects of all possible intermediate states included.
The correction to the tree-level amplitude of the nucleon β-decay may be written as [19] V where is the first Nachtmann moment of the structure function F (0) 3 [27,28]: Again, we may rewrite the integrals over Q 2 and x by those of ω and q 2 .The Jacobian is given by dQ 2 dx = 2M N x/νdωdq 2 .Then, the formula ( 24) can be written in the form with r and Q 2 also rewritten using ω and q 2 .(See ( 20) and (22).) The structure function is defined as the part including the -tensor i µναβ in W µν as and ) By looking at the component of (µ, ν) = (i, j) = (1, 2) we obtain which is to be combined with (26).
There is a recent lattice computation of the same quantity but for pion [29], which has been used to estimate the contribution for nucleon [30].They used the coordinate-space integral instead of the momentum space integral given above.Another method to use the Feynman-Hellmann theorem has also been proposed [31].

V. ENERGY INTEGRAL
We apply the method outlined in Sec.II to the computation of the total N cross section (19) or the loop correction (26).
On the lattice, one can calculate the forward Compton-scattering amplitude from a four-point function that includes the operators to create and annihilate the initialstate nucleon |N (0) .We omit the spin index of the nucleon state, but an avarage over the nucleon spin is assumed as indicated by 1 2 pol. .On the second line, the correlation function is rewritten using the transfer matrix e − Ĥt , and a Fourier transform of the current is introduced: Let us consider the ω-integral in (19).The factor L µν /(q 2 − M 2 W ) 2 in front of W µν has the form ω l with l either 0, 1, or 2, depending on the components µ and ν, up to an ω-independent factor.The W -boson propagator can be approximated by a constant, 1/M 4 W , for low-energy scatterings.For the electromagnetic scattering, on the other hand, 1/(q 2 ) 2 = 1/((ω − M N ) 2 − q 2 ) 2 has to be multiplied.We write the factor in front of W µν collectively as K(ω, q), so that the integral to be performed is written as dω K(ω)W (ω), where the indices other than ω are omitted for simplicity.(The definition of the integral kernel K(ω) will be slightly modified.See below.) Along the line of ( 8), the ω-integral at a fixed q can be rewritten in the form The indices µ, ν as well as q are omitted on the left-hand side.The upper and lower limits of the ω-integral has to be included in the definition of the kernel operator K( Ĥ).For the total cross section (19), the kernel K(ω) is proportional to ω l (l = 0, 1, or 2) until it hits the upper limit M N + |q| where it vanishes discontinuously.The power l depends on the form of the leptonic tensor L µν , so that (30) corresponds to each term appearing in (19) and we have to add them together in the end.
As we have already mentioned, the approximation of K(ω) using the Chebyshev polynomials becomes more difficult when the target function is discontinuous.In order to avoid this problem, we propose to smooth out the discontinuity.For example, we can replace the Heaviside step function θ(M N + |q| − ω) to realize the upper limit by is introduced to specify the range of the smoothing (or smearing).To be explicit, we may take the kernel function of the form where the factor e 2ωt 0 is introduced to compensate the small time evolution that appears when we define |ψ µ (q) ≡ e − Ĥt 0 Jµ (q)|N (0) .The kernel function (31) does not reflect the lower limit of the ω-integral at M 2 N + q 2 .That is because there is no state contributing to the integral below the lower limit, which corresponds to the elastic scattering.In fact, the forward Compton-scattering amplitude C JJ µν (t; q) in ( 29) behaves as e − √ M 2 N +q 2 t at large time separations.Therefore, we can safely extend the lower limit of the ω-integral to zero.
The kernel function should be adjusted to treat the electromagnetic scattering with a photon exchange or the γW -box correction to the β-decay (26), as they have complicated prefactors.But, the basic strategy is unchanged.
In order to demonstrate how the Chebyshev approximation of K(ω) works, we show some examples of K(ω) and its approximations in the following.We set the lattice cutoff 1/a = 2.4 GeV and the lattice size L = 32, which are typical in today's lattice QCD simulations.
The nucleon mass is taken at the physical value M N = 0.96 GeV.The small time duration t 0 is taken to be minimum, t 0 = 1, in the lattice unit.The plots in the following are all in the lattice unit.
We choose the kernel function K(ω) in (31) with the smoothing parameter σ = 0.2, 0.1 and 0.05 (lattice unit).|q| is taken as the smallest non-zero momentum on the lattice: |q| = 2π/L, so that the upper limit is at ω = 0.60.
Even when the approximation for the kernel function is not very accurate, the integral over would then contain larger errors when |q| ∼ M π , as we discuss in more details below.
Fortunately, the single pion threshold is relatively easier to treat theoretically using baryon chiral perturbation theory, and the potential error may be corrected when higher precision is required.The thresholds of more than one pions would be less problematic since such contribution is added on the spectrum of fewer pions and thus their impact is less significant.
Once the approximation is constructed, it is straightforward to evaluate the ω-integral.
A practical problem is that the Chebyshev polynomials T * j (x) involve huge cancellations among different orders of x k (k = 0, 1, ..., j), since the coefficients may grow as fast as 4 k with alternating signs.As a consequence, ψ µ (q)|T * j (e − Ĥ )|ψ ν (q) / ψ µ (q)|ψ ν (q) with large j may have an error larger than 1.Since the Chebyshev polynomials are constructed such that |T j (x)| ≤ 1 is satisfied, the terms whose magnitude is greater than 1 due to statistical fluctuations can easily destroy the whole approximation.In order to avoid such a problem, we should add constraints | ψ µ (q)|T j (e − Ĥ )|ψ ν (q) / ψ µ (q)|ψ ν (q) | ≤ 1, when we determine their numerical values from the lattice data.This can be done using a constrained fit.
Namely, we extract the Chebyshev matrix elements ψ µ (q)|T * j (e − Ĥ )|ψ ν (q) / ψ µ (q)|ψ ν (q) from the correlators ψ µ (q)|e − Ĥt |ψ ν (q) / ψ µ (q)|ψ ν (q) by a fit with the constraints.Due to the statistical error of the correlators, the higher order Chebyshev polynomials are not well determined but limited within the range between −1 and +1.Such poorly determined higher order terms are still useful to estimate potential errors due to the truncation of the polynomial.(See [21] for details).
Systematic error due to ignored higher order terms in the Chebyshev approximation can be estimated from their coefficients.Since the shifted Chebyshev polynomials, T * j (x), form a orthonormal basis of functions in 0 ≤ x ≤ 1 and those of higher j represents rapid variation of the true function, we expect that the high-j coefficients c * j are suppressed when the true function varies only slowly.This can be confirmed with the examples of the approximations given above.Figure 8 shows how c * j decreases for higher polynomial orders j.The plots are presented for different l's and momenta |q|'s.For all the cases, the coefficient c * j basically falls exponentially as j gets higher.The rate of decrease is faster when the smearing width σ is larger (σ = 0.2), and |c * j | becomes O(10 −3 ) or even smaller already at j = 10.This size represents the relative error in the final result, because the Chebyshev matrix elements ψ µ (q)|T * j (e − Ĥ )|ψ ν (q) / ψ µ (q)|ψ ν (q) are constrained between −1 and +1 by construction.The decrease of c * j becomes slower for smaller σ.Thus, for more rapidly varying kernel functions, much larger polynomial orders are necessary to achieve the same precision.When the lattice data are available only in a limited range of the time separation between two currents, this sets the limit of the method, and we have to choose sufficiently large σ such that the truncation error is under control.The results should then be extrapolated to the limit of σ → 0.
An example of the extrapolation σ → 0 is shown in Fig. 2 of [1] for the inclusive semileptonic B decays.It is for the same type of function with l = 2, and the data show that the dependence on the smearing parameter σ is mild and becomes essentially flat for some small values of σ ( 0.1 in the lattice unit).
For the N scattering, for which lattice computation of the relevant amplitude is not available in the form useful to perform this analysis, we consider a simple model that describes an elastic scattering plus a single pion production processes.The elastic channel corresponds to a delta function δ(ω − M 2 N + q 2 ) in the spectral function W (ω) in (30).The single pion production begins at s = (M N + m π ) 2 for s = ω 2 − q 2 .The spectral function typically has a shape between the elastic contribution (a peak at ω = M 2 N + q 2 ) and the N π continuum is unknown, so that the height of the peak is also arbitrary.The task is then to obtain a convolution integral of W (ω) with the (smeared) kernel K(ω).Apparently, the contribution of the elastic channel is underestimated by the smearing, while the N π contribution is likely overestimated because the spectrum is an increasing function.
We estimate the error due to the smoothing with the model described above.This The N π contribution depends on the details on the kernel function, thus on l, and the error increases with l (filled circles).It appeared that the error in this particular case is quite significant for the N π continuum contribution, and one probably needs σ = 0.1 or smaller to control the extrapolation, which seems to be well described by a linear dependence on σ 2 .The error would probably cancel between the elastic and N π contributions.When the elastic contribution is relatively large, the cancellation becomes stronger and the total error might not be as substantial as the estimate for the N π contributions suggests.
The actual error depends on the details of the spectrum.The value of σ can be chosen at the analysis stage of the lattice data, and thus the analysis can be repeated for various σ without extra computational cost.How much one can reduce σ depends on how large time separations the lattice data exist at without overwhelmed by the statistical noise.

VI. QUARK-LINE CONTRACTIONS WITH TWO CURRENT INSERTIONS
The challenging part in the lattice computation of the inelastic N cross section is the calculation of the forward Compton-scattering amplitude (29), which involves two current insertions.It may be obtained utilizing the Feynman-Hellmann technique as applied in [12,13], but here we consider the more conventional approach to contract the quark lines.
In order to extract the matrix element (29) we need to compute the following two-point and four-point functions on the lattice.For completeness, we also define the three-point function: with t sep ≡ t snk − t src , τ 1 ≡ t 1 − t src and τ 2 ≡ t 2 − t src .The spatial momentum components are discretized on the lattice as p i = 2πn i /L with an integer n i for a lattice of the spatial extent L. In the end we take the initial nucleon momentum to zero, p = 0, and the source position x can be fixed, e.g.x = 0, without loss of generality.The operators O 1 and O 2 are the weak or electromagnetic current with some Lorentz structure.

A. Preparation
The nucleon interpolating operators at source and sink may be defined as for P (proton) and N (neutron).Roman letters a, b, ... stand for the color index, while the Greek indices α, β, ... distinguish the spinor components.S is a diquark spin matrix, which we choose S = Cγ 5 = γ 1 γ 3 and S = −S.
We define the quark propagator / D −1 of a flavor q as ( / D −1 q ) ab αβ (y, x) = q a α (x)q b β (y) .The flavor q may either represent up (u) or down (d) quark.We call it the propagator, because we assume it describes a propagation in the same direction in time as the nucleon propagator.In the following, we use a notation F q (y, x) ≡ ( / D −1 q ) ab αβ (y, x).We will define the backward propagator as well, shortly.
The nucleon two-point function is obtained by taking all possible contractions as where T is a projection matrix and the repeated indices are summed.The projection matrix is set to extract desired nucleon state.For instance, to sum over the nucleon spin it is set as T = diag(1, 1, 0, 0).
To compute the three-point function, our strategy is to build a sequential propagator and then to contract at the location of the operator.To do so, we define the backward (sinksequential) propagators B u,d , which depend on the flavor of the valence sequential backward quark propagator: where S u,d 's are diquark propagators that reflect the structure of the nucleon interpolating operator at the sink.They are composed of u-and d-or two u-quark forward propagators as The connected three-point functions with a neutral current insertion is then constructed as where the current operator has a γ-matrix structure Γ.

B. Four-point functions
Here we describe a scheme for the computation of the general four-point nucleon correlation functions for both flavor-diagonal (or neutral) and flavor-changing (or charged) currents.

Neutral current
For the neutral currents the current insertions and the corresponding quark lines from the source to sink nucleon operators are shown in Fig. 10.
In order to contract at the location of the current J (1) , z 1 , we need to compute three j X a 7 N q 9 7 F 6 x l y w o y 1 R l k U d 9 P y 9 2 0 F H O F Y s 8 L n j q H m t U K p W q t X C U 5 l g i Z X M r k / 0 A m n x M R W 6 5 d T y B d z 6 p 2 t f 7 j 0 1 z z 3 K e 5 Z l 7 V 5 k 9 x m h S 1 H d K b H f Q W h W 9 r d i q 8 P M R d q t P z r E 4 X c M Q q 6 5 9 O D B v r c I F k 6 Z Y R i P + o c a c j w f L P q a I o / e s b p E H S g U f X t G P g P 4 P g Z f L q y s 3 j m x S 0 Y a k c J e n Y U T z 2 G M T l A 9 5 1 k w R a b n Y Y 1 S l i 1 5 c E E P q Y s d c 3 p y A K 9 + q 6 G q d 7 W N 2 + / J d l p H N 9 F b q A I t 0 0 a 3 0 Q f o D u o g q 3 C z 8 G 5 h v / B h s V / 8 v P h l 8 a s M e m V N c t 5 A c 6 P 4 z T 8 v z k w R < / l a t e x i t > (2) d < l a t e x i t s h a 1 _ b a s e 6 4 = " t k i A f A N 5 q 9 m l f R E 8 e l A F B E 6 u J (2) (2) u < l a t e x i t s h a 1 _ b a s e 6 4 = " p G l P 1 u t w w D D 3 u J (2) d J (2) u J (2) d , and different current sequential propagators.One is the current-sequential propagator (C): and two others are the current-sink sequential and sink-current sequential propagators called G and E, which were introduced in [32], where we should replace the CP-violating operator with the current operator J (2) for our purpose (See Fig. 11).Then, the sink-current sequential propagators E are defined as The current-sink sequential propagators G can be expressed using the diquark propagators S as In what follows, we omit the arguments (x, p, ...) in the propagators for simplicity.Then, the connected four-point functions for the neutral currents are given using F, C, E, and G < l a t e x i t s h a 1 _ b a s e 6 4 = " s q 1 J r g 3

C
< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 0 C 2 0 j V j t u 9 p P t V 9 q v x q / 1 y / W F + t X p t A L C y X n C j q 0 6 h / / C 2 m M P 5 A = < / l a t e x i t > / z e 9 8 c 4 q e H t w P w Z c P P S d n T m e J 5 h E 5 5 O o 0 n J 6 H I R q j t Z J n z 11. Propagators for computing the connected nucleon four-point functions. as (2) (2) where Tr denotes a trace over color, spinor indices as well as the location of the current z = z 1 , which is limited on the time-slice t 1 .The first term in either case corresponds to the diagram given on the left in Fig. 10.The G-type contribution does not appear for J d J d .

Charged current
The quark-line contractions are more complicated for the charged current, for which the quark flavor changes.All the diagrams for the charged current insertions J − J + and J − are shown in Fig. 12.The first and third diagrams in both cases resemble the first and third ones given in the J u J u current in Fig. 10, where we can use the backward sequential propagators B q and E q for these diagrams.On the other hand, for the second diagram the Wick contractions of the sink and source quarks are different from that of the J u J u or J d J d current (See Fig. 13), so that the backward sequential propagators G defined above can not be used.
To evaluate these crossing diagrams, it is convenient to define generalized nucleon twopoint functions: k S qb j ]q a i , where we introduce fictitious valence quarks with three different flavors q i , (i, j, k = 1 , 2, or 3).Therefore it is understood that the Wick contraction for the function where S J − J + [F 1 , F 2 ] is a diquark propagator for crossing diagrams, and the concrete form will be given below.C + is a current-sequential forward propagator defined as Obviously, C + = C u,d in the isospin symmetric limit.We note that such diquark propagator for crossing diagrams can be obtained by taking the functional derivative of the generalized nucleon two-point functions (N ijk [F q 1 , F q 2 , F q 3 ]).For example, the diquark propagator corresponding to the first diagram in Fig. 13 is given as where the function S q l ijk [F q l , F q l ] is defined as For example, S q 3 132 [F q 1 , F q 2 ] is The full contribution of the diquark propagator S J − J + is then obtained by the sum of four diagrams, For the other two diagrams, they are evaluated using the backward propagators B u and E, since the diagrammatic structure is the same as the flavor-diagonal case.However the quark flavors in the third (E-type) diagram in Fig. 12 should be different from the original E-type propagators, since the quark flavor for backward propagators will change in the case of the charged current.For this purpose we introduce generalized functions E q 1 q 2 defined as + J (2) u < l a t e x i t s h a 1 _ b a s e 6 4 = " O F T 6 e H C e + k 5 Q C j 8 g n + (2) + < l a t e x i t s h a 1 _ b a s e 6 4 = " F 0 W K S Q 6 v h a X S o r s 4 2 e 6 0 s K 5 + J (2) u J (1) + , J + J d , and J (1) + .
Then the first diagram in Fig. 12 is Tr[B u •Γ (1) •C + ], and the third diagram is Tr[ Next, we consider the J + J − insertions.A similar argument to the previous analysis for J (1) + can apply.As shown on the right-hand side of Fig. 13, there are four crossing diagrams like those of J + .Another diquark propagator for crossing diagrams S J + J − is given as where C − is the current-sequential propagator Then the backward current-sink sequential propagator for the crossing diagrams In summary, the full connected contribution of the nucleon J − J + and J + J − correlation functions are written as where the first term in either case corresponds to the diagram on the left in Fig. 12.

β-decay amplitude
Finally, we also consider the proton to neutron transition amplitudes that are given by four-point functions of J + J u , J u J + , J + J d , and J d J + .Again, the crossing diagrams exist for these correlation functions.As shown in Fig. 14, we have two types of diagrams for each operator combination.We also note that there are four different types of the contractions for each diagram, which are common to all diagrams, since the flavor changing current is commonly given by J + .
As for the diagram on the left of each operator, it is convenient to introduce the following diquark propagator for the proton to neutron transition For example, the diagrams given in Fig. 15 are obtained by where B N P d is the sequential backward propagator of d quark, Using B N P d we also obtain the backward sink-current sequential propagators for the left diagrams of J u J + , and J + J (2) d in Fig. 14, For the diagrams on the right for each operator combination in Fig. 14, we use the following four different types of the backward current-sink sequential propagators, where each of diquark propagators is defined as Using these propagators, the full connected contribution of the nucleon J + J u , J u J + , J + J (2) d , and + correlation functions are written as We note that the neutron amplitudes of the Compton scattering as well as the β − decay are also obtained from these functions by interchanging u and d flavor indices.Therefore, all the necessary quark-line contractions for the β-decay are given by these formulae.
As implemented in [12,13], some of these matrix elements may be obtained utilizing the Feynman-Hellmann technique.However, such numerical implementation would be more demanding.For example, the four-point functions with an insertion of J − J + are obtained from a second derivative of the two-point function in the presence of two external source + .When the currents are flavor-changing, it corresponds to a modification of the Dirac operator to a flavor-dependent one / D(m q ) → / D(mu , and one has to compute its inverse at various values of 1,2 to take a numerical derivative and then an extrapolation 1,2 → 0. Furthermore, in order to control the time separation between the two currents one has to repeat the method with the source terms only at a given time slice, so that the total cost may be as high as the direct computation of the four-point function. We use this technique to verify the quark-line contraction codes with the flavor changing currents.We place the external field corresponding to J (1,2) multiplied by a small parameter 1,2 at a given time slice τ 1,2 and compute the two-point function.By taking a numerical derivative we confirm that the results agree with what we obtain with the four-point functions including two current insertions.The details of the analyses and numerical results are shown in Appendix A.

VII. DISCUSSIONS
This work proposes a method to compute the inelastic N scattering cross section using lattice QCD.The key steps for computing the total cross section are to decompose the forward Compton-scattering amplitude into the contributions of different energies and then to integrate with an appropriate weight factor that represents the phase space.Instead of literally performing this program, we consider the spectral decomposition only virtually and realize the energy integral by identifying the weight factor as an operator constructed from the Hamiltonian.This weight operator between the two currents in the Compton amplitude can be reconstructed from the corresponding lattice computation.The essential point is that we can avoid the explicit spectral decomposition, which is a well-known example of the ill-posed inverse problem.
In the standard analyses of DIS, the cross section is measured depending on Q 2 and x, from which one can determine the structure functions F i (x, Q 2 ).(The subscript distinguishes distinct kinematical structures, but the details are not important in this discussion.)In the parton model, the structure functions are further decomposed into the contributions of partons (quarks and gluons) and written in terms of the parton distribution functions (PDFs).The basic assumption here is that the power corrections of the form 1/Q 2 and higher can be neglected, which is justified only above Q 2 ∼ several GeV 2 .In fact, in the low Q 2 region, there are resonance structures in F i (x, Q 2 ) near x ∼ 1 due to low-lying energy states.Such resonances may not be treated by perturbation theory.Even using the lattice QCD calculation, treating the individual excited states is a very challenging task, since they are actually multi-particle states like N π, N ππ, • • • , whose spectrum becomes dense on large volumes.Our proposal is to consider only a sum over such states, which is an well-defined quantity and the correspondence with the lattice observable can be established.In other words, instead of the structure functions F i (x, Q 2 ), we only analyze their weighted integrals.
The problem of the standard analysis at low Q 2 is related to the assumption of quarkhadron duality.In the perturbative QCD calculation one takes the quark and gluon external states, which are unphysical, instead of hadronic states.Such an assumption can be justified when one sums over a certain range of kinematical variables because the non-perturbative effects which are enhanced near the resonances can be avoided [33].In the analysis of DIS, the duality has been observed to be satisfied after averaging over some appropriate range of x [34][35][36] (see also [37] for a review).Conversely, some sort of smearing (or averaging) of the experimental data is necessary to compare with perturbative QCD, and it becomes more prominent in the low Q 2 region.So far, no quantitative measure on how much smearing has to be introduced for a desired precision has been known.Our proposal is one way to define such smearing on a theoretically solid foundation.In this paper, we have focused on the calculation of the total cross section, but it can be easily extended to the cases of partially integrated cross section, which plays the role of the smearing.Fully non-perturbative computation is possible in our framework, and no assumption of duality is necessary.
One of the main themes in the study of nucleon structure is the determination of PDFs.
With our method, the x-dependence of PDFs is not accessible as we need an integral over the energy of final hadronic system.Still, some moments of the cross section can be computed by choosing an appropriate weight function including some power of x, for instance.Once various moments are obtained, they can be used to constrain the overall shape of PDFs.
We have to be careful, though, because the integral over ω while fixing q 2 corresponds to an integral along a curve on a (x, Q 2 ) plane as shown in Figure 3.They cannot be simply written using the moments conventionally defined as 1 0 dxx n F i (x, Q 2 ) at a fixed Q 2 .One may consider instead integrated moments of the form 0 dxx n F i (x, Q 2 ), which can be calculated using the method introduced in this work.Given the fact that the Q 2 dependence of F i (x, Q 2 ) is a subdominant effect and can be understood using evolution equations, one could still gain some information from such analyses.
Towards realistic computation of the inelastic N scattering cross section, there are a couple of challenges we would be confronted.One of the major limitations, which is actually common in lattice QCD calculations, is that the spatial momentum |q| one can reach would be less than a few GeV/c.The momentum is of course limited due to the discretization effects, but in practice the limitation rather comes from large statistical noise of higher momentum correlators, for which the signal is rapidly overwhelmed by the noise.As a consequence, it will be very challenging for the lattice computations to reach the DIS kinematical region, a few GeV/c.The energy region of interest for the neutrino experiments, T2K and DUNE, on the other hand, is below a few GeV, which would be within reach.
Even at small momenta, the computation of the forward Compton amplitude is a challenging task, because it requires a saturation of the ground state nucleon on the both ends to prepare N | and |N .The signal-to-noise ratio for nucleon decreases as exp[−(m N −3m π /2)t] for large time separations where the ground state would dominate [38].Even for a nucleon two-point function, it requires a lot of effort to make sure that the ground state has been reached, and there have been many studies to investigate the ground-state saturation for three-point functions, which are relevant for various nucleon charges.(See, for instance, [39][40][41].More details and a full list of references may be found in [42].)The up-to-date simulations have a time separation between the nucleon source and sink operators about 1 fm, which is not ideal to sufficiently suppress the excited states of mass gap about the pion mass but is a necessary compromise.The insertion of two currents separated from each other in the time direction requires even larger time separation than the computation of the three-point functions.The computations carried out so far [11][12][13] do not include extensive tests of the ground-state saturation.We emphasize that the signal-to-noise problem is common for all lattice calculations, especially for those of nucleon properties, and various methods are being studied to improve the situation.

VIII. CONCLUSIONS
This paper describes a new method that enables us to explore a class of new applications of lattice QCD.It is about inclusive processes such as the lepton-nucleon scattering without specifying the final hadronic states.On the lattice, the corresponding quantity is the forward Compton-scattering amplitude calculated at various Euclidean time separations between the two inserted currents.The necessary quark-line contractions are summarized.The lattice observable can be related to the physical cross section.The method has been already successfully tested for inclusive B meson decays [1].
The method opens a new possibility to study the νN scattering in the low-energy region, which is relevant to the neutrino oscillation experiments, such as T2K and DUNE.So far, theoretically solid analysis has been possible only for elastic scattering, for which the form factors computed on the lattice can be used, and for deep inelastic scattering, which is described by perturbation theory.The region between these lowest and high energy scales can be treated fully non-perturbatively using the technique proposed in this work.
Although this paper describes only the total cross section in order to be explicit, the proposed method can be utilized to compute other related quantities, such as moments of x or other variables.The change of the analysis is only to modify the smearing kernel, and the extension is straightforward as long as it does not introduce more discontinuities.
Through such moments, one can extract more information about the process and the nuclear structure.
The actual computation is yet to be carried out.The computational cost is significantly more demanding for the forward Compton-scattering compared to that of the three-point function used to extract the form factors. Furthermore, the signal-to-noise problem is severer, so that the realistic calculation would be a substantial challenge in lattice QCD in the next decade.N 321 [F 2 , F, F 1 ], where the first term corresponds to the sum of the first and third diagrams in Fig. 12, and each of the other four terms corresponds each of the crossing diagrams in Fig. 13.
, it is useful to consider the following combination Since in the numerical calculations we only have the data at some discrete values of 1 2 , we fit the data to obtain the linear coefficient of the 1 2 term in CJ (1) J (2)  2pt−conn = f ( 1 , 2 ).The fits are shown in Fig. 16, where an excellent linear dependence on 1 2 may be observed.The comparison between the direct and the background field method is shown in Table I.We can confirm a precise correspondence.

4 FIG. 3 .
FIG.3.Integral paths in terms of q 2 and ω on the plane of x and Q 2 .The spatial momentum is chosen as |q| = 2π/La × k with k = 1, 2, 3, 4 (from bottom to top).With a lattice cutoff 1/a = 2.4 GeV and the lattice extent L/a = 32, the physical lattice extent is L = 2.62 fm.The momentum |q| drawn in this plot corresponds to the range between 0.47 GeV and 1.88 GeV from bottom to top.

Fig. 4 FIG. 4 .
FIG. 4. Kernel functions K(ω) with l = 0 (top), 1 (middle), 2 (bottom) and their modification due to the smearing.The modified functions are shown with σ = 0.2, 0.1 and 0.05 in the lattice unit.Other parameters are described in the text.The solid purple curve represents a mock data for the spectrum used to test the method (in an arbitrary unit).

20 FIG. 5 .
FIG. 5. Kernel function K(ω) and its Chebyshev approximations plotted in the lattice unit.Top panels are for l = 0 while the middle and bottom panels show them for l = 1 and 2. From left to right, the smearing width gets narrower (σ = 0.2, 0.1 and 0.05).The approximations are written with the polynomial order N = 5, 10, and 20.The energy integral is truncated at ω = M N + |q|;
which comes from the imaginary part of an one-loop diagram describing a creation and annihilation of a N π pair in the s channel.(Strictly speaking, this is valid only for two scalar particles.To be more realistic, one should use for instance the heavy-baryon chiral perturbation theory.The qualitative feature near the threshold is expected to be unchanged, though.)It is shown in Fig.4(purple curve) with an arbitrary unit.The relative weight

FIG. 9 .
FIG.9.Extrapolation of the integral to the of σ → 0. The estimate is given with a mock spectrum.The integral at finite σ divided by that of σ = 0 is shown.The contribution from the elastic channel (gray), which is independent of l, extrapolates from below.The N π contribution depends on l (l = 0 (black), 1 (blue), 2 (red)) and is overestimated at finite σ.The horizontal axis is σ 2 .The data points correspond to σ = 0.15, 0.1, 0.05 as well as 0.02.The last point is already very close to σ = 0.

< l a t e x i t s h a 1 _ b a s e 6 4 =
" 5 6 a V i R k 3 m A c c i g W d b V g Z / d i g / S o = " > A A A R 3 n i c 7 V d N b x t F G J 6 U j 4 Y F 6 r p c k H q Z E h v Z k m 3 N r p 0 4 I Y p U 4 I L C p R + Y V P I a a z / G 9 s j 7 x e 4 4 I V n t k Q s S 4 g g S J 5 A 4 I K 7 8 A y 7 8 A Q 7 9 C Y h j k b h w 4 J 3 d c b 1 2 7 N B U E S p V Z m V 7 9 p 3 n e T 9 m 3 3 0 0 N g O H R Z y Q h 2 t X X n j x p Z e v r r + i v P r a 6 9 c K 1 4 s 3 P o 7 8 S W j R j u U 7 f v j A N C L q M I 9 2 O O M O f R C E 1 H B N h x 6 Y 4 / f F + s E hD S P m e x / x 4 4 D 2 X G P o s Q G z D A 6 m f v H a z 0 q 5 0 y k r u k m H z I s H b D g J a d I d c T P o K b p F P U 5 D 5 g 2 n 6 5 y N T w J m c Q F 6 z K H H 3 s j w b D B E F B L w h n w U 6 1 O r 7 f O x b c q u 6 O 2 x 3 b Y w f K S r t 8 7 7 n 3 n q q 6 d d K G b 7 O Q E / J i 4 9 p r r 7 / x 5 l u b h e J W 6 e 1 3 r t / Y v v k 4 9 M a B S Y 9 N z / a C p 4 Y e U p u 5 9 J g z b t O n f k B 1 x 7 D p E 2 P 0 i f A / O a V B y D z 3 c 3 7 m 0 5 6 j D 1 z W Z 6 b O w X S y X b j T N e i A u V G f D c Y B n R T L 9 + 6 V i 4 m N s 9 G 5 z 0 w u H S m O n r l D 3 b X A E F I o 5 w 7 4 M O q m V s v j I T u n k 0 h V N M e Z j x F U E 5 D P o a B N + 7 z Y B Z K c f o k r u l r F O s e V O l G 0 G q k e Y j m m b i 3 j V p V 2 9 b D Y 9 Q P P 1 w e d P g 1 q O K B W L 0 7 C P d z x x v x I 1 U g N M / d I a 5 K e j D 8 s l h 1 m W k 3 P e 2 N z Z R F 1 n I Q y Y a I w d R 5 C I O c x v p K I R P B 6 m I I B 9 sP R S B L Y A Z k 3 6 K J q g I s W N A U U D o Y B 3 B c w C / O o n V h d 8 i Z y i j T a h i w 1 8 A k R i V y W / k e / K S / E x + I H + Q f x b m i m Q O w e U M vo 0 4 l v o n 1 7 9 6 7 9 H f r 4 x y 4 J u j 4 T R q K W e O + q g t u T L g 7 k u L W I U Z x 5 + e f / 3 y 0 Z 2 H 5 e g j 8 i 3 5 E / g / J y / I j 7 A C 9 / Q v 8 7 s H 9 O E 3 S / j 0 4 X k G u R z w L N 4 5 s U u 6 r E h h r 5 b h x D m M 0 D m q Z z K L S M F 0 n a i h j L K S O P h A D 6 n 5 j r k 8 e Q y S 1 F R a D 7 S d u x 8 n 3 b S J 3 k c f o g p 0 z D 6 6 i z 5 F 9 9 E x M g v P C z 8 V f i n 8 u v V 7 a a t 0 s 3 Q r h l 7 b S G L e R T O j 9 M G / e w C B L w = = < / l a t e x i t > u u d u u d J L b p X t d h 5 0 m L A b a c q K X b 5 1 z r 2 3 b l 0 f x Y 4 8 l n B C n m x c e+ n l V 1 5 9 b X N L e b 3 2 x p v X b 2 y / 9 T A J 0 9 i h X S f 0 w v i x b S X U Y w H t c s Y 9 + j i K q e X b H n 1 k T z 4 R 6 4 9 O a Z y w M P i M n 0 e 0 7 1 u j g A 2 Z Y 3 E w D b a 3 v j Z t O m J B N m S j N K Z T p d 6 9 f b t b V w o r Z 5 O L i D l c L p V I e h 6 M r c A F Q 0 I h Y D D i 4 8 w s r W 7 I E 3 Z B p 5 m u G b 7 / 7 x i R b A G K O I T 0 6 J A r J q T J 6 Z e 4 Y e l N b H H c 2 C O a o Z L m M S 5 X e k 4 c J g k E 6 A u Q c Y l q q 7 r W a u J s W k E b l j H j R N c 6 z W P F H M V W N M a Z D L C 3 1 4 s 8 i w U q j q m r 4 j P G x 9 i K 4 / D s h G i H B y o O U 3 6 i G 0 T F L D g x D k g / D 7 j A E q D D A t M m f R l 1 e q z Uf e a 6 H q 1 2 Z B c 7 I n I 7 c l R r R r l W Z B n F Y W S N e k M a y y j 9 n M 5 D 3 B P h j o p w h y I c c C F a z E b j m

FIG. 10 .
FIG.10.Connected diagrams for the nucleon 4 point functions with J D q d t p P u z m 5 2 p y D d 9 G 7 8 A h 4 8 a e L B + D G 8 + A U 8 8 B G M B w + Y e P H g m + 0 i F A T j b H b 3 z X u / 3 3 u / e f O q g S s j T e n u w O D Q m c z Z c 8 P n r Q s X L 1 0 e G R 2 7 8 i T y 2 y E X 6 9 x 3 / f B p l U X C l U q s a 6 l d 8 T 5 z C t 7 j b 3 V e 7 z 2 6 t 5 a N 7 9 B 3 9 B v q f 0 t 3 6 S c 8 g d r 6 w d + v i r U 3 p + i p 4 3 c H c 3 k Y O b l z p k s s q S i w V 6 f h z D 2 0 o A N T h z I b p l H 6 P 6 x m w q q l P H x w h u y j E 3 P c 2 C g V 7 Z m i b a + W J p b u p + M 0 D D f g N u R w Z O 7 C E j y A F V g H P v h 9 6 N r Q z a F b m e n M s w z L 8 B 5 0 c C D l X I W + l X F / A + a O g T A = < / l a t e x i t > z J(2) e X X K l d W X v 8 s C U c x 4 7 s s l G F 8 z 6 U J l y L g u 0 o o y e 9 F M a e + K / m e O / x Y 9+ / t 8 z g R Y f C p O o x 4 x 6 f 9 Q P Q E o w p M 3 Z W l L w 2 j d q N m O C 7 v i y D t i f 4 o 5 l l 7 o N y o Y z i M B 4 r H I u h P + 5 U Y j i P B l A Y 9 4 v D D Y E A D D w w J B w F B X w 1 S Z 2 r 1 Q p W I M c 9 S y 7 R 9 / 5 8 x W n w B i l R m 1 C T v K c M B 2 Y p / g e v U a m C q c H 2 N m P Z 6 k z S 2 c X 7 N + u 1 y v 2 V u N r Y N J 4 r D i P b b P R 4 3 s S t H v D P x o 0 L c D k d q x y K k i U W w Y 3 9 A O r m H b a P m C 8 + T f B b X L e K S P O b U 3 G Z x m C Q w q o 5 G z C D m e g O n 2 Q x X d + 1 p X y G p H 9 No g N P c 7 9 p a O 5 J U B E 0 c c 6 + J s d a 0 V U j a I F P X M 9 S B U A N M 4 z g 8 2 I F R t g r a H E u H z G A g s e g P S v l j j 1 S a k D 5 z q 4 1 w o F 3 9 i P e 9 V e a n D N W 7 N l 3 n v e d 9 3 3 m n S e r u y b 3 B S H P V 2 7 c f O 3 1 N 9 5 c f a t U f r t y 6 / a d t X c e + 0 7 g U X Z E H d P x n u q a z 0 x u s y P B h c m e u h 7 T L N 1 k T g 5 B L Q N P B 2 e w o u c 4 v g A j u O s B 0 R X 1 6 8 a U x 5 J 3 m p T f U q z u T K j s m m y n T K h Q o H 3 D Q N T z v p y M m B b m p 0 v L 5 D 1 u P 2 k / a C q Z d c N M x M k 7 s + k 1 u 0 q I V B u O A 8 l T a 1 6 n M i B m y 9 T d Z P R l y w 8 5 A 5 W y / r s + V R R 7 6 r U R Y q r U X S e V V 5 v b S 6 X i C Z M + o 7 X 0 9 n x G W R n h Z a P t P K W E 4 X X A v t l S 7 P K + n p X A a m x H P q 3 8 u U f u 4 t 1 s 8 p O W i 1 5 l W t t q 6 l 8 V o a r 6V x W d Q 8 l 2 1 5 E n t x V U U G 7 v e D L 8 O a U o / w / f 5 H M I M G 2 Z R F M t v I v 7 f K x + K r b Q x I X 4 A L o 3 9 n g z R J P P D s R E k n G y g d h 8 7 a y i r q I g M 5 i K I A W Y g h G w m Y m 0 h D P n w 6 S E E E u W D r o R B s H s x 4 v M 5 Q h E r g G w C K A U I D 6 x i + h / D U S a 0 2 P M u Y f u x N Y R c T / j z w x K h K / i A / k R f k G f m F / E n + W x g r j G P I X E 7 h V 0 9 8 m d u / / c 2 7 j 1 5 e 6 G X B r 0 C j i d f S n A U a o H a c K 4 f c 3 d g i q 6 C J / / H Z t y 8 e f f K w G n 5 I f i B / Q f 7 f k + f k d 6 j A P v 6 H / v i A P f x u S T 4 D + D 6 F W B a s L G Z O s q T F O z L g a h l O n s M Y n a G t X G T p K T O 9 j N c o 9 j J S P / h A D y n T H T M 7 e Q w 3 d a f Z e q B u 3 P s 0 7 a Z V 9 B 7 6 A N W g Y 3 b R P f Q 5 O k R H i J a / L v 9 c / r X 8 W + V Z 5 e / K v 5 W X C f T G S u p z F x X G L f Q / i b A R E A = = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " l x + V R t B T h c y k b e O c 8 p p 0 m n Y I a q Q = " > A A A O T 3 i c 7 V b N b x t F F J + U r 3 R N 0 g Y u S L 1 M S F z Z w r F m 1 3 H t N L J U w Q X 1 l L a k r W Q b a z / G 9 s j 7 p d 1 x 0 m S 1 / 0 A l z h w 4 g c Q B c e Q K J y 5 c O X D o n 4 A 4 F q k X h H i z H 9 l d f 6 R N B I d K j B V 7 9 s 3 v v X n v N 2 9 + W c0 1 m c 8 J e b Z y 5 Y 0 3 3 3 r 7 n d W r U u n d t f V r 1 z f e e + g 7 U 0 + n h 7 p j O t 5 j T f W p y W x 6 y B k 3 6 W P X o 6 q l m f S R N v l E r D 8 6 o p 7 P H P s z f u L S v q W O b D Z k u s r B N N g o 3 e x p d M T sY M h G U 4 + G U v n A k H p j 3 1 V 1 G j Q t K 5 S S d c 4 m p y 7 T e Q R K f e i J P V Z t A w w + h a 3 t E R 8 H v d R q O N x n p z Q M 5 L o S B V q A E W k n I J f D 5 i Y d c i z 1 I G N O n + C K K l e x y n F l h 9 R 3 a 6 S 6 j 9 O V r u 4 5 v g 8 7 9 A V I y a H k e q O K g 3 A Z t B C w 3 i x C K 2 o W q Q G R 2 t V 9 q T f y V H e M g y i X n Z 2 u a 6 r M r m G P G j W M n S n v y A q p Y W Z 3 l A Z J k w H Y k H o 1 r J l T W o t A S i s G 7 Z E 0 i y x U h r q V R I p R S j X c l 8 o W M w y T Z p R o S Q W k w A e Y M z s w s I + j c b a s p x z F F D X z r n p h q V C z t q T m v S T R l k h U K 5 Z z z P g Y q 5 7 n H O M O h l y a i e O c n 5 5 Q l Y 8 + D 4 p Z 8 N h o z K V c O V n K B R p 0 O b 8 S E Z G j K C 0 0 Y S C 3 l F t p p + T B y J j Q F x z Z o p J e R k V 2 2 D M c K o t b o o i K y Z B 6 t m P Q s 0 a V a 6 K b q x g H P S i G R 5 o Q A J t h s D 3 d D s E h g 5 M L Q J P A j f Y M X O Q W w g 1 2 H G 4 7 P L q 9 e N u Y 8 Y 7 z U u r K Z Z z J p R 3 j T e X Z l A s V D p l p G p 5 6 3 B W T j m a q + m R z l 2 x G 7 S f s B V M / v m i Y m i Z z f S q 2 a O o W B u W C 8 5 T b u l V d E H F K N 9 t k 8 3 j M O D 0 L m b P 1 0 z 4 7 P 2 q i w X J z m X Z e V l 9 f Q V 4 X C O W l t D c v q H P i s k x Q C y 2 f a m W k p 4 s l Q n 2l u / P 6 y 2 m H 1 G / 9 1 2 I K 5 y Q v O e o Z a Z 3 9 1 z k n o f + a d k L d u 4 3 8 u c p C S P 8 X x d d O F F s L R L F 1 M V H M c 9 k W J 7 E X V V V k 4 O 7 g o 8 + D i l w N 8 d 2 B A T N o k G 1 R J L W N / C u r e C y 8 1 Z Y T S P I m X B i D 6 1 u k T q K B 5 y d y M t l C y T h w N 0 e m X z x / c v l 8 O b p J v y B + Q / 9 f k G f k Z K r C P / t S / v U f v f 3 V O P k P 4 P o F Y F q w s Z 0 6 w p E Y 7 U u D q P J w 4 h w k 6 R T u 5 yM J T Z H o R r 3 H k Z S R + 8 I E e k m c 7 Z n 7 y E O 7 q b r 1 5 T 9 m 6 8 3 H S T a v o B v o Q V a B j W u g O + h Q d o E O k l 7 4 o / V D 6 s f T T 2 q 9 r L 9 b + X k + g V 1 a S y f u o M N a v / g P Y s x R / < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " 2 k a h t 2 8 l a t O s q e 2 L b 6 1 9 t P O Y D R o = " > A A A N / H i c 7 V Z P b x t F F J + U P w 1 r T B u 4 I H G Z k B j Z q m P N O j F x G 0 W q 4 I J 6 S l v S V r K N 2 Z 0 d e 6 f e f 9 o d J y S r 5 Q P w B T h w A s E B V V z 5 A l y 4 I w 7 9 C I h j k b h w 4 M 3 s b r 3 r 2 O k f i Q N S Z u X d 2 f d + v z f v v X n z v G b g 8 E g Q 8 n j l 0 i u v v v b 6 5 d U 3 t M q b 1 b e u X F 1 7 + 1 7 k T 0 P K D q n v + O E D 0 4 i Y w z 1 2 K L h w 2 I M g Z I Z r O u y + O f l Y 6 u 8 f s T D i v v e p O A n Y w D X G H h 9 x a g g Q D d e 0 3 / o m G 3 M v H v H x N G S J V r M O t E w k + O Q 0 4 F Q o e Q 5 j J 5 5 t e B Y I I g a r e W N h x / 1 c a v k i 4 q c s i f V W 2 3 U X Y 6 S n G S g Q s J 7 D R k L r g 4 + C f Y H r h t 7 A h s D 1 L d L a a Z L G H s 4 1 P R r 6 U Q Q L D C S o h G p 1 G j h O l k H b B a j e 2 i 5 D 6 0 Z Z 3 W 3 s a f 1 x a A Q 2 j p U v W 1 u 9 w D G 4 1 8 Q h s 5 o Y + 1 O x f 5 0 0 M f f 2 2 7 s k 9 w V Q I x Y 2 s e l M W X M Rp l 2 0 t A z V b i R 7 W s 3 l l u W w W U b M L F Z S S g e I Z 3 K I a g + r 8 V R N 8 7 j T u D p F K i 2 p S i G b z x O y q

FIG. 14 .
FIG. 14. Diagrams for proton to neutron transitions for J t e x i t s h a 1 _ b a s e 6 4 = " v W H t a z B N i l M O 3 w w s m b a s a o g j + X d 2 f d 9 3 p n 3 m Z 1 5 Z r S q Y 3 B P E P J y a f m N a 2 + + 9 f b K d e n G O + + + d / P W 6 v u P P L v r U n Z M b c N 2 n 6 i K x w x u s W P B h c G e O C 5 T T N V g j 9 X O Z 4 H / 8 S l z P W 5 b X 4 y T X N Y P 0 h U W O y J D M e Y O 7 b g d Y B D s u F m y b E I 2 K 1 d C j N u D K c 1 W k 4 q 9 N w p l n O Y 0 A R Z Z e 3 9 d Q k o H I / v w x n K q c 9 I e v U e C S s Y r o p V 8 q z n 4 w U l D 5 t K g + / n 1 G U R s B G j M 4 0 L 1 u N m E s N y 9 Z Y M g V 3 5 H I w E U s Y + w 3 I X I Q L 3 4 e h 6 / m b 3 c 0 e B P T h J A c 0 b n h n f w A e 5 N a D p W n b w r J

r e 2 T 9 q c 4 F
u 2 g y Z W s m k 2 p y q 8 z S 0 s o Z P G b F 9 a u u o k 0 r v M U A e C X E N 7 M k B 5 V 3 e F H O Q 3 n z L M Z J m 8 L M q j y b I A / I 0 9 V V 5 c H 9 t 5 a S 3 z n p b 6 7 N N k p j o c E L D c 6 n w Q 3 d c x T K f L k 2 T m r n K M d D G j y 7 0 F y 6 5 p Z x T p W Z T l i z 2 0 p K Y y / 1 7 D 8 v 1 c 0 l S X M 5 D I / b d h d H 4 I X 8 X k H 5 v c Q j s M 5 c 9 j q U d + g M 9 P 8 T 4 J z n / 7 l I b y 4 p m t e J d 6 r P T a / t D L w Q 4 I U A x 5 b 4 o 7 A E 5 e T W B q m

FIG. 16 .
FIG.16.Results for C2pt−conn .The linear fits to the lattice data are also shown.

TABLE I .
2) − ¥hline Re(C 4pt−conn ) −3.908 0.9488 4.114 5.098 Comparison of the results from the direct and the background field propagator methods.