A stochastic method for computing hadronic matrix elements

We present a stochastic method for the calculation of baryon three-point functions that is more versatile than the typically used sequential method. We analyze the scaling of the error of the stochastically evaluated three-point function with the lattice volume, and we found a favorable signal-to-noise ratio suggesting that our stochastic method can be used efficiently at large volumes to compute hadronic matrix elements.


Introduction
Hadron structure calculations in lattice QCD have emerged as a powerful tool for providing comparison to and guidance for experiments, see e.g.refs.[1][2][3][4].Examples include moments of generalized parton distribution functions, as well as form factors. Lattice computations of such quantities have been carried out at several values of the lattice spacing allowing the continuum limit to be taken.In addition, small, sometimes even physical, values of the pion mass have been employed in the calculation of these quantities leading to an improved understanding of their quark mass dependence and how they approach the physical point.Unfortunately, studies of excited state contributions [5][6][7][8][9] suggest that for some quantities these effects can play an important role as a systematic uncertainty that affects hadronic three-point function computations.Safely accounting for these effects requires large statistics, hence methods to speed-up these calculations are highly desirable.
The progress in nucleon matrix element calculations on the lattice has prompted an effort to go beyond the simplest observables and to pursue a larger variety of interesting hadronic quantities.The evaluation of more observables will deepen our knowledge of hadron structure and provide a more comprehensive test of QCD.However, using the conventional sequential method [10] to calculate these matrix elements, it is necessary to perform a new computation of the needed quark propagators for each observable of interest 1 .This then leads to a high computational demand if many physical quantities are being sought.
In this paper, we describe an alternative approach, based on a stochastic method, that allows us to obtain a large class of observables with only a single computation of the propagators.To this end, we employ stochastically computed all-to-all propagators.Since a calculation based on a stochastic evaluation of propagators may lead to very noisy results, we perform a detailed study to determine whether the stochastic noise can be controlled with a moderate number of stochastic sources.We determine the signal-to-noise ratio as a function of the lattice size to test whether our stochastic method can be used in large volumes, such as 48 3 × 96, that are used in contemporary lattice computations.These questions are addressed specifically for the example of the axial charge g A of the nucleon.We would like to emphasize that we do not perform an analysis of systematic effects, since our goal is solely to test the stochastic method.
During the course of our work, a similar stochastic approach was employed in the calculation of meson three-point functions [11], where a (heavy) all-to-all propagator is estimated stochastically.There it was found that the stochastic method is competitive and in some cases even superior to the sequential one.Of course, it is not guaranteed that the same conclusions hold for baryon matrix elements, since those are subject to a stronger exponentially decreasing signal-to-noise ratio [12].
This paper is organized as follows: in Sec. 2 we outline our stochastic method, in Sec. 3 we present the results of this method and compare to those obtained by the fixed sink method, and we summarize our findings in Sec. 4.

Stochastic method for baryon three-point functions
Quantities that are needed for investigating hadron structure can be extracted by computing matrix elements of baryons with local operators.In lattice QCD these baryon matrix elements are obtained from baryon threepoint functions in Euclidean space-time that are of the form O( y, τ ) = q( y, τ )Γq( y, τ ) .
Γ represents a combination of γ-matrices and covariant derivatives and we have used translational invariance to set the source point to zero.Naively one would need an all-to-all propagator from every lattice point ( y, τ ) to all points ( x, t) for the evaluation of the above three-point function, which is, of course, prohibitively expensive to calculate.Such a demanding computation can be circumvented by applying the sequential method to perform the summation over the spatial coordinates of either the sink or the current [10].
For the example of the fixed sink method, the momentum as well as the time slice of the sink are fixed and an additional inversion for each flavour is needed.An alternative approach is to estimate the all-to-all propagator stochastically, which is the method that we explore in this paper.
A generic three-point function of a baryon B is defined as AA .The insertion time of the operator is denoted by τ .For illustration let us now consider the three-point function of a proton and the operator dΓd.We use the interpolating field I p possessing the quantum numbers of the proton, namely where C = iγ 0 γ 2 is the charge conjugation operator.In terms of quark propagators, the connected piece of this three-point function reads where x = ( x, t) and y = ( y, τ ) and where the up (down) quark propagator is denoted by S (u) (S (d) ).P is an appropriate spin projector, which we will specify later.
The sequential method with fixed sink makes use of the fact that we can perform the sum over x in Eq. ( 1) by an additional inversion.Then a generalized propagator for fixed time slice, projector P and sink momentum p is obtained, as indicated by the shaded area in the left diagram of Fig. 1.This renders the explicit calculation of an all-to-all propagator unnecessary.
Our alternative method uses a stochastic estimate of the all-to-all propagator appearing in three-point functions like Eq. (3).Such an estimate is obtained via where M is the Dirac matrix.In the above equations we have suppressed Dirac and color indices.η r is a random source obeying The stochastic method is diagrammatically illustrated in the right diagram of Fig. 1.Using the stochastic estimate, we can decompose the double sum in Eq. ( 3) into the product of two single sums, which is significantly computationally cheaper than the naive double sum and reads δd where x = ( x, t) and y = ( y, τ ).We have suppressed the average over the number of stochastic samples, cf.Eq. ( 4), and used the γ 5 Hermiticity property of the Dirac matrix to obtain the above expression.As before, we use superscripts to denote the quark flavor.The drawback of the stochastic method is that we have to average over a sample of N S stochastic sources.This requires N S inversions compared to just twelve (one inversion per Dirac and color index) in the sequential method.However, a major benefit of this method is its flexibility, since we do not need to fix the spin projector or the sink momentum, nor even the sink time slice, in principle.

Assessment of the stochastic method
To test the applicability of the stochastic method, we need to determine how large N S must be in order to keep the stochastic noise under control.This depends on the observable of interest.To be concrete, we compute a relatively simple benchmark observable of nucleon structure, namely the nucleon axial charge g A , using our stochastic method.This quantity can be obtained from a nucleon matrix element of the isovector axial current.For the evaluation of matrix elements of the nucleon, we need to introduce the zero momentum nucleon two-point function, We then examine the ratio R g A of the nucleon three-point and two-point correlation functions, In the limit of large Euclidean time separations, R g A converges to g A up to the renormalization factor Z A , In this paper we use the value of the renormalization constant Z A = 0.757(3) determined non-perturbatively [13,14].
In order to demonstrate that the stochastic method can indeed produce results with a reasonable computational effort and is potentially competitive with the sequential method, we performed a benchmark calculation with N f = 2 + 1 + 1 flavors of quarks.We employed twisted mass fermions at maximal twist with a lattice spacing of a ≈ 0.082 fm determined from the nucleon mass [14], a pion mass of m π ∼ 370 MeV and a volume of L 3 ≈ (2.6 fm) 3 .In Fig. 2, we show R g A (t, τ ) obtained using the stochastic method as a function of the insertion time τ for a fixed source-sink separation t = 12a.We compare to the value R g A (t = 12a, τ = 6a) obtained using the sequential method.For different values of τ close to the middle of the plateau the picture is similar.We use spin-color diluted random Z(4) vectors as stochastic sources, We have used a fixed number of gauge field configurations N gauge = 460 in both our stochastic and sequential approaches.Our observation is that using at most N S = 4 spin-color diluted stochastic noise vectors per configuration for the estimate of the all-to-all propagator is sufficient to reach the same statistical accuracy as with the sequential method.
In terms of inversions, N S = 1 corresponds to using the same number of inversions as in the sequential method i.e. (4 × 12) per gauge field configuration, for every additional set of stochastic sources we need (2×12) additional inversions to obtain the forward and back propagators, respectively.Thus this method would require four times more inversions.We would like to remark, however, that in the stochastic approach we can compute the correlation functions of proton and neutron without additional inversions, in contrast to the sequential method.In addition we used 3 operators k = 1, 2, 3 in Eq. ( 5), and correspondingly 3 spin projectors for the stochastic method, again without the need of additional inversions when using the stochastic method.Our observation is that this procedure reduces the statistical error by about a factor √ 6, corresponding to a factor of about 6 in statistics.We would like to remark however that this effective gain in statistics is rather specific for g A and will change for other observables.
Having demonstrated that it is in fact possible to compute g A using the stochastic method with a reasonable computational effort compared to the sequential method, we would like to know how the situation changes when the volume is varied.A potential danger of our stochastic method is that the number of stochastic sources required to reach the same precision as the sequential method may increase rapidly as the volume increases.
To study the volume effects, we use N f = 2 flavors of maximally twisted mass fermions, instead of the N f = 2 + 1 + 1.We expect that the stochastic noise should not noticeably depend on the number of dynamical flavors, and for the N f = 2 case, there exists a series of four different volumes at the same value of the lattice spacing, a ≈ 0.082 fm [16] and a pion mass, m π ≈ 370 MeV [17].These volumes are V = L 3 × T , where L/a = 16, 20, 24, 32 and the temporal extent of the lattice is T = 2L in all cases.This enables us to thoroughly study the volume dependence of the stochastic method over a relatively large range of spatial volumes, from about (1.4 fm) 3 to (2.8 fm) 3 .This corresponds to 1.95 m π L 3.9.
We performed an analysis using the stochastic method on a fixed number of gauge field configurations N gauge = 200 at each of the four volumes.The source-sink separation is fixed to 12a in all cases, which corresponds to about 1.067 fm.In the left panel we show results for a lattice extent of L/a = 24 and in the right panel we show results for L/a = 32.The dashed lines represent the error when using the sequential method.A fixed number of gauge field configurations Ngauge = 200 was used, with a pion mass of mπ ≈ 370 MeV and a lattice spacing of a ≈ 0.082 fm.
In Fig. 3 we show the relative error of g A as a function of the number of stochastic sources for two of the four volumes L/a = 24 and L/a = 32, where also the sequential method has been applied.For both volumes a convergence towards the error of the sequential method can be seen, which looks better for the larger volume, where the error is close to the error of the sequential method for N S = 4, to be conservative.This confirms the observation in the N f = 2 + 1 + 1 calculation mentioned above, which was done at about the same physical volume.The error of the error is not shown, but it is roughly of order 10% of the error.
We show the combined statistical and stochastic error of g A obtained using the stochastic method with N S = 8 for four different volumes in Fig. 4. In order not to be subject to a systematic error we do not fit the ratio given in Eq. ( 5) using an estimated plateau range τ between source and sink.Instead, we take the renormalized ratio R g A (t = 12a, τ = 6a) and its error in the middle between source and sink as our estimate for g A , where contributions from excited states are expected to be the smallest.For the larger volumes the error scales like the one of the sequential method, V −0.5 [18].Therefore the plot is consistent with the error being dominated by the gauge noise for larger volumes, which we have demonstrated above, see Fig. 3. Thus this method appears to work as well or even better at the larger volumes typically employed in current calculations, an observation which, of course, needs to be verified for other quantities and different physical situations.q q q q 0.1 0.2 0.5 The suitability of the method is demonstrated in Fig. 5 where we compare  ) and Gp(Q 2 ) computed with the stochastic and standard sequential method (refered as exact in the legend).For the stochastic estimate we use a fixed number of stochastic sources NS = 4 and we use one projector for the proton correlators.The computation is performed using N f = 2 + 1 + 1 configurations fixed to Ngauge = 500 with a pion mass of mπ ≈ 370 MeV and a lattice spacing of a ≈ 0.082 fm.
the nucleon axial form factors G A (Q 2 ) and G p (Q 2 ) [14] computed using the standard approach and our new method.For this comparison we use the same number of configurations and one projector.Since we use four stochatic noises diluted in colour and spin the computation is four times more expensive producing errors that are comparable to the exact case.However, the new method allows to compute the three point functions of the neutron as well as for three different projectors for free thus compensating for the increases cost.This compared with the fact that one can consider different final states e.g. a nucleon carrying momentum for free makes the new method more versatile.

Conclusion
We have applied a stochastic method for the calculation of nucleon matrix elements using spin-color diluted time slice sources.We have taken the case of the nucleon axial charge as a typical example of a three-point function to explore the method.Our conclusions for our test case is that the error is comparable to the error of the sequential method already at a moderate number of stochastic sources, namely four spin and color diluted timeslice stochastic sources, when using the same number of gauge field configurations.In this particular case we have effectively increased the statistics used in the stochastic method by a factor 6 through averaging over neutron and proton correlators and using 3 different currents, which does not require the computation of new propagators.Moreover, our results indicate that the convergence behaviour in the number of stochastic sources, N S , appears to improve when the volume is increased.
Since the stochastic method needs of O(10) more inversions (we needed N S = 4 for g A but this maybe different for other matrix elements) it is competitive with the sequential method when one computes O(10) more types of matrix elements.In the case of g A the increase in computational effort is easily compensated by computing 6 different matrix elements with 3 different spin projections for the proton and the neutron.Thus, even with the additional computational overhead, the great versatility of the stochastic method explained in this paper outweighs the sequential method when many baryon matrix elements or form factors are computed.

I
B is the baryon interpolating field and A and A summarize the indices depending on the quantum numbers of the baryon B, which are appropriately contracted with the function ζ (B)

Figure 1 :
Figure 1: Diagrammatic illustration of the sequential method through the sink (left) and the stochastic method (right).

Figure 2 :
Figure2: Plateau region of the ratio Rg A (t = 12a, τ ) obtained from the stochastic method on one N f = 2 + 1 + 1 ensemble with a pion mass of mπ ≈ 370 MeV and a lattice spacing a ≈ 0.082 fm.We use NS = 2, 4 and 6 spin-color diluted stochastic time-slice sources, with the source located at the sink.The source-sink separation is 12a and we compare to the standard sequential method, of which we show the ratio Rg A (t = 12a, τ = 6a) with the light gray band.The dark gray band indicates the PDG value[15].

Figure 3 :
Figure 3: The filled circles show the relative error on the value obtained for gA as a function of the number of spin-color diluted stochastic sources NS per gauge field configuration.In the left panel we show results for a lattice extent of L/a = 24 and in the right panel we show results for L/a = 32.The dashed lines represent the error when using the sequential method.A fixed number of gauge field configurations Ngauge = 200 was used, with a pion mass of mπ ≈ 370 MeV and a lattice spacing of a ≈ 0.082 fm.

Figure 4 :
Figure4: The filled circles show the error of gA using the stochastic method with a fixed number of stochastic sources NS = 8 on a fixed number of gauge fields Ngauge = 200 for four different lattice sizes.The pion mass is mπ ≈ 370 MeV and the lattice spacing is a ≈ 0.082 fm.The dashed lines indicate a scaling proportional to V −1/2 and V −1 and are solely meant to guide the eye.The error of the error is not much bigger than the symbol size, hence we do not show it.Please note the double logarithmic scale.

Figure 5 :
Figure 5: Comparison of GA(Q 2) and Gp(Q 2 ) computed with the stochastic and standard sequential method (refered as exact in the legend).For the stochastic estimate we use a fixed number of stochastic sources NS = 4 and we use one projector for the proton correlators.The computation is performed using N f = 2 + 1 + 1 configurations fixed to Ngauge = 500 with a pion mass of mπ ≈ 370 MeV and a lattice spacing of a ≈ 0.082 fm.