Precision Study of Excited State Effects in Nucleon Matrix Elements

We present a dedicated analysis of the influence of excited states on the calculation of nucleon matrix elements. This calculation is performed at a fixed value of the lattice spacing, volume and pion mass that are typical of contemporary lattice computations. We focus on the nucleon axial charge, g_A, for which we use about 7,500 measurements, and on the average momentum of the unpolarized isovector parton distribution,_{u-d}, for which we use about 23,000 measurements. All computations are done employing N_f=2+1+1 maximally-twisted-mass Wilson fermions and using non-perturbatively calculated renormalization factors. Excited state effects are shown to be negligible for g_A whereas they lead to an O(10%) downward shift for_{u-d}.


Introduction
Understanding nucleon structure is one of the fundamental goals of lattice QCD.Such an endeavor is becoming more realistic as present day calculations are being performed closer to the limit of physical quark masses, small lattice spacings and large volumes [1].Thus a direct comparison of results from lattice calculations and experimental measurements is becoming feasible, allowing us to probe QCD as the underlying theory of the strong interactions.
In order to establish that lattice QCD can indeed provide results that address this challenge, we are focusing, in this work, on a number of benchmark observables for which QCD is expected to produce the right results.However, even with the significant advances of the past few years, there is presently an unexpected discrepancy between lattice calculations and experimental measurements of physical observables such as the nucleon axial charge, g A , the average momentum of the unpolarized isovector parton distribution, x u−d , and the charge radius of the nucleon [2][3][4][5].Clearly, these discrepancies need an explanation.Naturally, this demands a careful study of the systematic effects that play an important role in lattice calculations, i.e. lattice artifacts, finite volume effects and non-perturbative renormalization.Indeed, these systematic effects are currently being addressed by many lattice QCD collaborations.In this work, we focus on the investigation for two key observables, namely g A and x u−d .These quantities are determined on the lattice by computing suitable ratios of 3-point and 2-point correlation functions that reach a constant plateau value for asymptotically large Euclidean time separations.
There are several studies of g A and x u−d at different values of the lattice spacing and various volumes [6][7][8][9].For the values of the pion masses considered, these results reveal that taking into account lattice spacing and finite volume effects is probably not sufficient to reconcile the lattice calculations with the current experimental value of g A and the phenomenological extractions of x u−d .In addition, there is a recent study indicating that the discrepancy persists even for pion masses quite close to the physical value [9].One might still argue that calculations performed at precisely the physical point might eliminate these discrepancies but that would require a strong pion mass dependence making such an explanation increasingly less plausible.
An issue that is currently under scrutiny within the lattice community concerns excited state contamination.These are sub-leading contributions in Euclidean time correlation functions that can cause a systematic effect in determining the desired nucleon-nucleon matrix element.In order to understand whether the physical plateau appears at larger Euclidean times than used in present day calculations, we need to examine the relevant 3-point functions at Euclidean times that are as large as possible.The difficulty associated with probing such sub-leading contributions to these ratios is that the signal-to-noise ratio decreases exponentially fast with Euclidean time.
To tackle this problem, we have performed a dedicated high precision cal-culation of the 3-point functions for g A and x u−d .In particular, we have used about 7500 measurements for g A and about 23000 for x u−d .With such large numbers of measurements, we are able to calculate the correlation functions for g A and x u−d for large Euclidean times to sufficient accuracy.This allows us to study carefully possible excited state contributions, which is the goal of this work.
This paper is organized as follows.In section 2, we describe briefly how the matrix elements for the observables of interest in this work are calculated in Euclidean field theory and how the excited state contributions enter the calculation.In section 3, we give the details of the calculation and explain the open sink method used in this dedicated analysis.In section 4 we present our results and in section 5 we summarize and conclude.

Nucleon Matrix Elements and Excited State Contributions
In order to make the paper self-contained, we introduce the basic definitions of the quantities studied here.In the following discussion, we assume an infinite volume, while in all our practical calculations, periodic or antiperiodic boundary conditions are taken as needed.The effects of a finite space-time lattice are felt in the standard finite-size effects of the matrix elements and the so-called thermal contributions that distort the Euclidean time dependence of correlation functions.These systematic effects are both addressed by the finite-size studies discussed previously and are henceforth ignored.
The zero momentum nucleon 2-point correlation function on the lattice is defined as where JN is a nucleon interpolating field that creates a state on the lattice with the same quantum numbers as the nucleon, Γ is a matrix acting in Dirac space, and the sum over the Dirac indices α is implicitly understood.The Euclidean time t denotes the separation between the creation time t source and annihilation time t sink of the nucleon and is often referred to as the source-sink separation, t = t sink − t source .We use translational invariance to set t source = 0 resulting in t = t sink .The nucleon interpolating field is where C = iγ 0 γ 2 is the charge conjugation operator.The transfer matrix formalism allows us to relate lattice correlation functions to matrix elements of operators.Application of the standard methods gives the spectral representation of the 2-point function in terms of the eigenstates of the transfer matrix or equivalently the Hamiltonian H.The resulting expression is which in the limit t → ∞, will be dominated by the nucleon ground state, In these expressions, m k labels the masses in the nucleon channel and m N denotes specifically the nucleon mass.Additionally, we have introduced the symbols J N and J(k) N for the overlap of the interpolating fields with the k th eigenstate of H. Strictly speaking, the limit t → ∞ cannot be realized on a finite lattice; in practice however, it suffices to take t large enough so that the correction coming from the lowest lying excited state can be neglected.For the evaluation of nucleon matrix elements that we are interested in, we additionally need to calculate 3-point correlation functions.They are defined as where O is a local field corresponding to the operator Ô of interest and Γ is an appropriately defined matrix acting in Dirac space.We denote by t the insertion time of the operator under consideration.Like the 2-point function, there is a spectral representation that can be derived from the transfer matrix formalism and reads The asymptotic limit of large Euclidean time again isolates the nucleon contribution as follows lim The desired nucleon matrix element 0| Ô |0 is then obtained from the asymptotic Euclidean time limit of the ratio of the 3-point and the 2-point function It is the main goal of this paper to investigate how large t − t and t should be so that the contribution of the lowest lying excited state -understood as a systematic error -becomes negligible within the desired precision.
Let us begin by discussing in what ways the excited states of the nucleon contribute to the nucleon matrix element calculated from the ratio of the 3-point to the 2-point function.The expressions given in Eqs. ( 1) and ( 3) illustrate how calculations of Euclidean time correlation functions can be used to determine matrix elements in the limit of large Euclidean time separations.However, in practice the finite time extension of the lattice prevents us from taking the asymptotic limits and therefore one has to carefully examine the sub-leading contributions usually ignored.In the following we assume that at the values of t we use, excited state contributions in the 2-point function can be ignored.This is plausible since excited state contributions to the 2-point functions are generically more suppressed than those contributions to the 3-point function.The reasoning for this is that the fields in the 2-point function are always separated by a distance t and those in the 3-point function are separated by t − t or t , both of which are smaller than t in practice.In particular, in the 3-point function we have a double limit requiring both t − t and t to be asymptotic.Let us now consider the leading contributions to the ratio of the 3-point and the 2-point functions originating from taking into account the contribution of the first excited state: where ∆M is the mass gap between the nucleon ground state and the first excited state.As can be seen, there are two additional time dependent contributions to leading order.

Lattice fermion action
For this work, we employ maximally-twisted-mass Wilson fermions [10].We use the gauge field configurations generated by the European Twisted Mass Collaboration (ETMC) with N f = 2 + 1 + 1 flavors, thus fully accounting for the first two quark generations.We refer to Ref. [11] for the details of our lattice formulation.Since we are aiming at a precise result, we concentrate on only one ensemble with a pion mass of m π ≈ 380 MeV and a lattice spacing of a ≈ 0.078 fm.This pion mass is chosen small enough to be relatively close to the physical pion mass but still large enough to ensure that finite size effects are suppressed.In addition, the propagators can be calculated with moderate computational cost which allows us to analyze a large statistical ensemble in order to obtain an accurate result.We emphasize that maximally-twisted-mass fermions realize an automatic O(a)-improvement for which no additional operator specific improvements are needed.Therefore, at the lattice spacing employed here, one expects that discretization effects are also suppressed.This is confirmed by direct calculations of these matrix elements at three different lattice spacings smaller than 0.1 fm in Refs.[2,6,12].
Although in this analysis we use twisted-mass fermions, the most important aspects of excited state contributions are expected to be universal and independent of the particular lattice discretization used.Thus the conclusions obtained here are of direct relevance to the calculations of other groups including those using different lattice actions.

Fixed Sink Method
An efficient computation of the connected piece of the 3-point function is possible by means of the so-called sequential source method.This technique requires two sets of propagators.The first are the forward quark propagators of the light flavors that are also used to compute the 2-point function and are independent of the operator and hadronic state as long as the same source is used.Those propagators are then used to build a sequential source for a second generalized propagator, again for each of the light flavors, that is specific for the hadron state we are interested in.The desired matrix element is obtained by contracting the corresponding operator with the free ends of these two propagators as illustrated in Fig. 1.This method has the advantage that we can use the same set of propagators for any choice of the operator insertion and hence it is the method of choice for the calculation of generalized form factors of a particular hadron.A disadvantage of this method is that the source-sink separation time t in Eq. ( 3) must be fixed before the sequential propagator can be computed.Obviously, changing t would require another set of propagator computations.Thus, a prudent choice for the value of t is mandatory.This is because, on the one hand, according to Eq. ( 3), a large source-sink separation is desirable for the suppression of excited state contributions.On the other hand, the signalto-noise ratio drops exponentially fast with the source-sink separation.A reasonable choice is therefore a source-sink separation at which the contribution from excited states becomes negligible compared to the statistical error.However, this can only be determined a posteriori and only after repeating the entire calculation for several values of t.Thus in practice, having chosen a reasonable value for t, one looks at the time dependence of the right hand side of Eq. ( 3) as a function of the operator insertion time t .If a plateau, as a function of t , is observed, then it is assumed that the excited states have been sufficiently suppressed and the plateau value is identified as the matrix element of interest.However, there still remains the possibility that the asymptotic plateau value has not been reached.Therefore in this paper, we carry out a thorough investigation of the excited state contributions using a more appropriate approach as described in the next section.

Open Sink Method
In the fixed sink method the summation over x in the 3-point function of Eq. ( 2) is done by the sequential inversion.One can instead perform the summation over y through a sequential inversion.Therefore we need to fix the particular operator that we are interested in and also the time slice t where it is inserted.The sequential source is then constructed at the operator rather than the sink and hence, the sequential propagator is now operator-dependent but state and t independent.For the example of an operator that involves only one lattice site, the sequential propagator is given by where the Roman superscripts are color indices and the Greek subscripts are Dirac indices.The generalization of this expression for operators that involve more than one lattice site and gauge links, such as derivative operators, is straightforward.The sequential propagator is obtained for all source sink separations t, thereby allowing us to study the effects of excited states.The open sink or fixed current method is illustrated in Fig. 1.Clearly, this method is not practical when a large number of matrix elements of different operators must be computed since for each new operator extra inversions must be performed.However, for our dedicated study of excited state effects in g A and x u−d it is clearly very useful.

Observables
As stated above, in this work, we concentrate on two relatively simple but nonetheless phenomenologically very relevant quantities.The first is the nucleon axial charge, g A , which plays an important role in the beta decay of the neutron and appears as a low energy constant in effective chiral Lagrangians.It has been precisely measured and it is straightforward to calculate in lattice QCD using the techniques described in the previous sections.However, the values obtained from various lattice QCD calculations are typically 5% to 10% lower while having themselves a statistical accuracy of the order of 1%, see Fig. 2 for the example of our own calculations.
The second observable is the lowest non-trivial moment of the unpolarized parton distribution function in isovector flavor combination, x u−d .This quantity is determined phenomenologically from a global analysis of deep inelastic scattering data, and the discrepancy between the phenomenological and lattice values is even larger, roughly 50% to 60%, see Fig. 2. For the precise definitions of the corresponding 3-point functions, we refer to Refs.[4,6].Though this does not affect excited state contamination, it is important to note that we use a non-perturbative renormalization of our bare matrix elements.The corresponding renormalization factors are calculated in the RI MOM scheme and are matched to the MS scheme at a scale of (2 GeV) 2 .For more details we refer to Refs.[13,14].The values of the renormalization constants used in this work are Z A = 0.774 for the renormalization of the bare g A and Z x = 0.998 for the renormalization of x u−d .

Results
In order to have a reference value, we first perform a calculation of the nucleon axial charge g A using the sequential fixed sink method with a fixed source-sink separation of t = 12a ≈ 0.94 fm.Gauge invariant Gaussian smearing of the quark fields, including APE-smeared gauge links, is used in order to improve the overlap with the nucleon ground state.In Fig. 2 we compare the value obtained for the N f = 2 + 1 + 1 ensemble to our previous N f = 2 results at various pion masses.As can be seen, the value we find for g A is in good agreement with the N f = 2 results obtained at nearby pion masses.
We then perform an analysis on the same N f = 2 + 1 + 1 ensemble using the open sink method.The time slice of the operator insertion was fixed to t = 9a.This was chosen to safely suppress excited state contributions from the source, as can be verified from the 2-point function.The result of the analysis using the open sink method is shown in Fig. 3. Clearly, the value of g A does not show any statistically significant dependence on the source-sink separation t.Hence, the plot demonstrates that there is no contribution from excited states detectable within the statistical accuracy of 2.5%.Note that, although t = 9a, a plateau for g A is already reached at t = 11a.It is worth mentioning that, in order to reach a comparable statistical accuracy as the one obtained when using the fixed sink method with t = 12a with 500 measurements, we had to perform roughly 7500 measurements when we take for example t = 18a.This was achieved by choosing randomly located source points with typically 2 sources per configuration.
As a second benchmark quantity, we have examined x u−d .As for the case of g A , we first determine the value of this observable using the fixed sink method increasing the statistics of the calculation presented in [14].For the open sink method, we have chosen the operator insertion time to be t = 11a.We expect that for this choice excited state effects from the source are sufficiently suppressed for this operator.We perform in total about 23, 000 measurements for x u−d using randomly distributed source points with 5 sources per configuration.With such statistics and at a sourcesink separation of t = 18a, we could equal the precision of the fixed sink method that was done with a source-sink separation of t = 12a using 1300 measurements.In Fig. 4 In the left panel, we show the relative deviation of ETMC lattice results for gA from the experimental value [15].In the right panel, we show the relative deviation of ETMC lattice results for x u−d from a result obtained from a phenomenological analysis [16].The lattice values for N f = 2 at the various pion masses are from Refs.[6,12].The filled (magenta) diamonds show the results using the N f = 2 + 1 + 1 ensembles.separation t.We also indicate the value obtained from the fixed sink method analysis as well as the experimental result extracted from a recent global analysis [16].
As can be seen, for this observable there is a shift of the value of x u−d and a plateau is reached at larger values of the source-sink separation than what we have used in the fixed sink method.Despite the fact that the results for larger values of t decrease they clearly do not reach the phenomenological value.In order to estimate any residual dependence on t, we determined the value of x u−d for an infinite source-sink separation by fitting the expected exponential behavior, to the lattice results with a fixed t = 11a.The result of this fit is x u−d = 0.22(1), which is 12% lower than the result of x u−d = 0.250(6), obtained using t = 12a in the fixed sink method.The error of the fit is estimated by varying the fit range and by comparing the use of a fixed parameter ∆M versus fitting ∆M directly.

Summary and Conclusions
In this letter we have performed precision calculations of g A and x u−d for a single ensemble of gauge field configurations with N f = 2 + 1 + 1 dynamical fermions employing a non-perturbative renormalization.We have investigated the behavior of these benchmark quantities as a function of the  source-sink separation in order to assess the influence of excited states on the current lattice results for g A and x u−d .This is particularly important given that excited states may play a role in explaining the presently observed discrepancy between lattice computations and phenomenological evaluations of several important nucleon observables.
We find that for the here considered pion mass of about 380 MeV and lattice spacing of a ≈ 0.078 fm, the contamination of excited states is negligible for g A , but for x u−d , the effect is of the order of 10% compared to our previous calculations, where the source-sink separation has been set to about 1 fm.This is an effect larger than the finite volume and lattice spacing effects we observe at this value of the pion mass, volume and lattice spacing.Moreover, this demonstrates that contributions from excited states are operator dependent and should be investigated separately for each operator.
One way to better control excited state effects is to use a variational method such as the generalized eigenvalue method [17,18].Recently, a new approach to deal with excited state contamination of hadronic matrix  elements has been developed and applied for the B * Bπ coupling in ref. [19].Whether the generalized eigenvalue method can improve the calculation of matrix elements of the nucleon needs still to be tested, though.
However, if the 10% shift for x u−d as we found here persists at smaller pion masses, excited state effects can not be the single dominating systematic effect responsible for the tension between lattice and phenomenology.Of course, we cannot exclude that at smaller values of the pion mass excited state effects might become significantly larger.Therefore, in order to clarify the deviation between lattice calculations and experimental determinations of nucleon matrix elements, a very careful and accurate analysis of systematic errors will be needed, taking into account the possible contamination of excited states as observed in this work.

Figure 1 :
Figure 1: Diagrammatic illustration of the sequential method through the sink (left) and the open sink method (right).
we plot x u−d as a function of the source-sink

Figure 3 :
Figure 3: Results for gA for a range of source-sink separations obtained from the open sink analysis on one N f = 2 + 1 + 1 ensemble.The light grey band indicates the result obtained from the fixed sink method using a source-sink separation of 12a and the dark grey band shows the experimental value.
fixed sink method, t sink = 12a ABMK fit open sink method

Figure 4 :
Figure 4: Results for x u−d for a range of source-sink separations obtained by means of the open sink method.The operator insertion was at a temporal separation from the source of t = 11a.The value (including errors) obtained from the fixed sink method using a source-sink separation of 12a is indicated by the light grey band.The phenomenologically extracted value is shown with the dark grey band.The blue solid line corresponds to a fit described in the text.