Variational Approach to the Calculation of gA

We present a lattice QCD calculation of the nucleon axial charge, gA, using a variational approach. After a brief outline of how the variational method is applied to the calculation of form factors, we present results for gA using this method. We find that ground state dominance is rapid, evident in the early onset of a clear plateau in the correlation function ratio proportional to gA. Through a comparison with results obtained via traditional methods, we show how excited state effects can suppress gA by as much as 8% if sources are not properly tuned.


Introduction
In recent years, lattice calculations have taken a tremendous step towards simulating QCD at the physical point. Algorithmic and technological developments have allowed simulations to probe at or near physical quark masses on increasingly larger volumes, with finer lattice spacings and vastly increased statistics. Calculations of the ground state spectrum have yielded results consistent to within a few percent of their physical values with well controlled systematic errors [1,2]. Naturally the next step has been to strive for this level of precision for the matrix elements of these states. Despite the remarkable consistency between lattice and experimental data for the pion form factor F π (Q 2 ), a complete description of other hadronic states, particularly the nucleon, has proven to be remarkably challenging [3,4].
The most notable shortfall is for the nucleon axial charge, g A ≡ G A (Q 2 = 0). In principle g A should be relatively simple to calculate. Being an iso-vector quantity, disconnected loop contributions are absent and as we have direct access to G A (0), we circumvent the need for extrapolations in Q 2 . Unfortunately, the lattice values for g A to date have been consistently lower than the experimental value by as much as 10-15% [5]. In an effort to account for these discrepancies, several studies have carefully examined the systematic errors present in the calculation [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. In this letter we will focus on the role of excited state effects.
Recently there has been an increased effort to understand and reduce the impact of excited states on form factor calculations.
In computing these quantities, it is well understood that to ensure excited state contributions to the correlation function are sufficiently suppressed, one needs large Euclidean time separations between operators. To choose a suitable time separation one should identify the time slices where the correlation functions take on their asymptotic form. For the two commonly used sequential source techniques, this is a relatively simple procedure for the fixed current method. One simply chooses a current insertion time, t C , once the asymptotic behaviour is observed in the two-point correlator. Results are extracted from the data once the asymptotic behaviour is observed in the threepoint correlator.
For the fixed sink method, one requires knowledge of the asymptotic behaviour of the three-point correlator a priori. Unfortunately, the temptation to use earlier sink times in order to obtain more precise results is inescapable. These results can suffer from excited state contaminations, even if a plateau is observed with t C . In refs. [10,16], it was found that for certain matrix elements, eg. x , the source-sink separations often used in the literature were not sufficiently large to suppress excited state effects. Nonetheless, as we move ever closer to the physical point one is naturally forced to choose earlier sink times as the signal degrades much quicker.
To counter this issue, new techniques are being devised to try and control the sub-leading terms to the three point correlator. The use of the summation method [18,20] has shown improvement upon the conventional approach, but the underlying excited states contributions are still present. It is not hard to imagine situations where these still impact significantly and alter the final result.
In this paper we take a somewhat different approach. Rather than reduce the impact of excited states through Euclidean time evolution, we seek to separate them out from the ground state at the source and sink. Drawing upon the techniques developed for excited state spectroscopy calculations, we will use the variational approach to construct interpolating fields that couple with individual energy eigenstates and use these to isolate the desired matrix elements [21,22]. An analogous approach has been presented in [23,24] for the study of B * → Bπ transitions and in [25] for the study of the axial charges of Nucleon excited states. Here we apply it specifically to g A to remove excited state contributions.
This letter is organized as follows. In Section 2 we will examine the variational method in the context of excited state spectroscopy and then outline how this method can be applied to the calculation of hadronic matrix elements. Section 3 outlines the details of this calculation. In Section 4 we present our results and compare our variational method with the traditional, single-operator approach to the calculation of g A . Section 5 is a cost-benefit discussion for the variational method. Finally we provide our concluding remarks in Section 6.

Variational Method for Matrix Elements
The 'variational method' [26,27] is a well established approach for determining the excited state hadron spectrum. It is based on the creation of a matrix of correlation functions in which different superpositions of excited state contributions are linearly combined to isolate the energy eigenstates. A diversity of excited state superpositions is central to the success of this method.
Starting from a basis of operators { χ i (x) | i = 1, . . . , N }, we construct a correlation matrix of two-point correlation functions, Due to the discrete nature of the lattice, we can decompose these correlation functions into a discrete sum over energy eigenstates, where the parameters Z α i ( p) are the coupling strengths of the interpolators χ i (x) with the energy eigenstate of mass m α and Γ projects out the desired parity. We choose new operators φ α (x) to be linear combinations with a suitable choice of coefficients v α i and u α j , such that these interpolators couple to a single energy eigenstate, From Eqs. (2) and (4) we find that the necessary values for v α i and u α j are the solutions of the following eigenvalue equations where the eigenvalue c α = e −m α ∆t . It is important to note that both (5) and (6) are evaluated for a given momentum p and so the diagonalisation condition is only satisfied when we project with the relevant coefficients as follows: v Thus the two-point correlation function for the state | α, p is We can extract the mass m α from G α ( p = 0, t) in the standard way.
To understand how we can utilise the variational method for use in form factor calculations, we must first identify the terms present in the three-point correlation function, where O(x) is the current operator to be inserted. Sandwiching the current between two complete sets of states we end up with three terms, the vertex amplitude, β, p ′ , s ′ | O(0) | α, p, s , and the coupling terms Ω| χ i (0) | β, p ′ , s ′ and α, p, s |χ j (0) | Ω , The coupling parameters take the same form as they did in the calculation of the two-point correlator with two key differences.
The inclusion of a current means that the initial and final momenta need not be the same. Furthermore, there also exists the possibility that the initial and final energy eigenstates are not the same. That is, the current can induce a transition between states. For this calculation the necessary expression is To isolate the matrix element from the three-point function, we construct a ratio in the standard way. In this work we shall use the ratio defined in [28]. For the state α the necessary ratio is, (12) Key to this approach is the utilisation of a basis of operators in which there is diversity in the overlap with various excited states. As there are a limited number of local bilinear operators for a given J PC , a great deal of work has been made by various groups in increasing the number of available operators. Here we choose to use fermion source and sink smearing as a method of extending our operator basis, as outlined in [29,30].

Calculation Details
For this calculation we make use of the PACS-CS (2+1)flavour dynamical-QCD gauge field configurations [31] made available through the ILDG [32]. These configurations are generated using a non-perturbatively O(a)-improved Wilson fermion action and Iwasaki gauge action. The value β = 1.90 results in a lattice spacing a = 0.091 fm, determined via the static quark potential. With dimensions 32 3 × 64, these ensembles correspond to a spatial length of L = 2.9 fm. As the intention of this paper is to examine whether the variational approach is an improvement upon traditional techniques, we will consider only the light quark mass that corresponds to m π ≈ 290 MeV. The resulting value of m π L = 4.26 is comparable to the values used by most groups.
In this work we are primarily interested in isolating the ground state and so have chosen to use a small variational basis upon which to perform our correlation matrix analysis. We use gauge-invariant Gaussian smearing in the spatial dimensions only with smearing fraction α = 0.7 [33]. We consider four levels of smearing with the optimal choice found in [33], these being 16, 35, 100 and 200, applied to the standard, local proton interpolator thus allowing for construction of a correlation matrix of dimension up to 4 × 4. In table 1 we list the rms-radii for our choice of smearing parameters. We choose to use variational parameters t 0 = 18 and ∆t = 2, again taken from [33], where it was found that this choice produced best balance between systematic and statistical uncertainties.
This vertex can be expressed via two independent form factors, the axial form factor G A (Q 2 ) and the induced pseudoscalar form factor G P (Q 2 ), as where q µ = p ′ µ − p µ and Q 2 = − q 2 . Using isospin symmetry, one can show that for the flavour-changing current A ud µ , the matrix element is equivalent to that of the iso-vector current A u−d µ , As we are interested in G A (Q 2 = 0), it suffices to consider the case where the incoming and outgoing momenta are the same, in particular we choose to work in the nucleon rest frame as this will provide the smallest statistical uncertainties. This will mean that the left and right eigenvectors required to project out the three-point function will now correspond to the same momenta.
We choose to insert our fermion source at t 0 = 16. For the calculation of the three-point functions we use a local axial vector current calculated using a sequential source technique with the current held fixed and inserted at t C = 21, at the onset of asymptotic behaviour for the projected two-point function. We choose to use µ = 3 for the current with the corresponding projection matrix being Γ ′ = Γ 3 = Γ 4 γ 5 γ 3 , where Γ 4 ≡ 1 2 (I + γ 0 ). The value for the axial renormalization constant Z A = 0.781 (20) was determined non-perturbatively in [34] using a Schrödinger functional scheme.
The resulting expression from which we extract g A is the ratio of the eigenstate-projected three-point and two-point functions, As a comparison, we also evaluate g A using a single correlator from smeared source to point sink and smeared source to smeared sink. These are indicative of results one would extract from a traditional approach.

Results
In Fig. 1 we present the bare values of g A with increasing sink time t s following the current insertion at t C = 21 for the smeared source to point sink, smeared source to smeared sink (both with 35 sweeps of smearing) and our variational method respectively. Between the traditional approach (upper two plots) and the variational approach (bottom plot), we can see significant differences in the overall shape of the correlation function ratio.
For the smeared source to point sink (upper plot) the Euclidean time suppression of excited state contributions manifests itself as a steady increase in the extracted value. This trend in the data does not have a clear endpoint and so there is no definite plateau. Guided by the χ 2 dof obtained via a covariance matrix analysis, the earliest time slice one could consider is t S = 25, but what is clear is that we are forced to consider fit windows uncomfortably close to regions dominated by noise.
By smearing the sink as well as the source, there is a definite improvement in the quality of the plateau. The excited state behaviour is again present as a steady increase in the value of g A , but somewhat suppressed. In this case there is a definite plateau observed at t S = 24, which is supported by the χ 2 dof . Unfortunately, this is again somewhat close to the region where signal is lost to noise.
In Fig. 1 (c) we see quite a different situation. Our variational approach yields extremely clean results with rapid ground state 21 Figure 2: An overlay of the results from fig. 1. The data sets have been offset from the time slice for clarity -the circles (blue) are the results for the variational approach, the triangles (purple) are the smeared source → smeared sink, while the squares (red) are the smeared source → point sink. The fitted value from the variational approach has been included (blue shaded region) to highlight where the traditional approach is consistent with the improved method.
dominance. The systematic rise in the data is no longer present and the onset of the plateau is within two time slices of the current insertion. In Fig. 2 we have overlaid the three datasets to highlight the excited state behaviour between the traditional and variational approach. If we look carefully at the variational approach, we can see that some excited state contamination is present immediately after the current, but this is short lived. It is worth noting that this is a consequence of the limited size of our variational basis. As is highlighted in [24], an n × n correlation matrix allows one to isolate out the n lightest states in the given channel and so the sub-leading contributions will come from the n th + 1 state. In the case of the ground state, these contributions will be short-lived due to the large mass splitting between the ground state and n th + 1 excited state. If one were to construct a basis whose dimension was the number of states in the given channel, then it would be possible to completely isolate each state.
What is of most concern in Fig. 2 is the lack of overlap between the results of the traditional approach and those of our variational method at t s = 24 and 25 where good fits can be made. In Table 2 we list those fits, for the three data sets with the strict criterion that the χ 2 dof lies between 0.800 and 1.200. In both data sets employing the traditional approach, we can obtain good fits with small uncertainties if we choose to begin fitting around t S = 25 or 26, but find that the results are significantly small. As we move the fit window to later times, the central value increases. Between 25-30 and 28-30 we observe a systematic variation of 6% in the value g A . A consistent result can be extracted from these datasets if we choose to fit at later time slices around t S = 28, but the resulting values have unattractively large uncertainties, as they are close to the onset of noise. It is clear that in this case, we have little control over the excited state systematics. In contrast, for the various fit windows on the data from the variational approach we find the variation in the fits is considerably smaller than the smallest statistical uncertainty.
It is worth considering how the level of smearing affects the extracted value of g A . In Fig. 3, we present the renormalized g A considering each of the four smearings used to construct our variational basis. What we find is a dependence on the level of smearing used in the calculation. It appears that for low levels of smearing the extracted result can be significantly lower, with the smallest level of smearing differing by up to 8% from our improved, variational result. From this evidence, it is clear that if the smearing level is not properly tuned at the source and sink, then excited state effects significantly impact the extracted result for g A .
In principle, one could tune the smearing so that the optimal overlap is observed with the ground state. By using a point source propagator and tuning the smearing through the sink via the two-point correlator, outlined in [35], one removes the need for expensive inversions for the tuning. Unfortunately, the optimal level of smearing depends on the quark mass, β value, momentum or operator. One must tune the smearing for each set of parameters under consideration. Immediately, one can find appeal in the variational approach as there is no longer a need to tediously tune the operators to match the ground state.

Cost-Benefit Discussion
A concern with the correlation matrix approach is the increased cost. For our implementation, we require 2 inversions per configuration for every smearing we include in constructing the correlation matrix. For n smearings we have a total of 2n = 8 inversions per configuration, as opposed to the minimum of 2. In Fig. 2 we can see that, for large Euclidean times, the conventional approach is consistent with the correlation matrix approach, albeit with larger errors. Thus it is worth considering what the required increase in statistical sample would be for the conventional approach to produce results with similar error to that of our correlation matrix method.
Given the error varies with the sample size N as ∆g A ∝ 1 √ N then the relative increase in sample required to obtain an error (∆g A ) desired is given by where (∆g A ) current is the error extracted with the current sample of size N current . Using the leading time-slice of the associated fit-windows as indicative of the uncertainty in g A , which for the smeared-smeared approach is t s = 27 and for the correlation matrix approach t s = 23, we find that Naively we expect a factor 4 increase in statistics, which would require fewer inversions than our correlation matrix method. However, we note that the peak value for the smeared-smeared approach is at time slice 26 and so χ 2 dof analysis would tend to favour earlier points around times 24-25. This is consistent with Table. 2. In the tradition of choosing the earliest possible fit-window to minimise statistical uncertainty, a more appropriate fit window would be times 24-31. Being conscientious of the rapid growth in error bars, the best choice would Table 2: Un-renormalized values of g A from fit windows which give a covariance matrix based χ 2 dof between 0.800 and 1.200. The datasets are identified as (a) Traditional approach with point sink, (b) Traditional approach with smeared sink and (c) Variational approach, where for the traditional approach we have selected the 35 sweeps of smearing. We note how the value of g A increases for the traditional approach as we move the fit window to later times. In contrast, the variational approach is stable across all windows with the desired χ 2 dof .  (3). This result is systematically suppressed, relative to the correct result of g A = 1.47(2), by excited state effects. While one could invest more super-computing resources to reduce statistical error, in this case one will only get the wrong answer very accurately if one does not take care in fine-tuning the source.
To further understand this we note that using the variational approach, ground state domination occurs earlier in Euclidean time thus allowing the current insertion at an earlier time. For this particular ensemble, ground state dominance for the nucleon occurs at time t = 21, so our choice for t C is ideal for the correlation matrix method. For the smeared-smeared approach with 35 sweeps of smearing, ground state dominance does not occur until time t = 23. This is why the peak value is systematically low. The downwards shift for small source smearings is the result of the current being too early, sampling both ground state and excited state contributions to the matrix element. This also gives rise to the smearing dependence illustrated in Fig. 3. Thus for a more comprehensive comparison, one requires a new simulation with t C = 23, two time slices later. However, we can still get some insight from our present analysis into the required increase in statistics. For the the ratio of three-to two-point functions, ground state dominance occurs 6 time slices after the current insertion, so with t C = 23 one would be considering a fit window commencing at t s = 29 as opposed to t s = 27 considered earlier. Here we have a factor 11 increase. As the variational approach enables one to: 1. rapidly isolate the ground state following the source, thus enabling an earlier current insertion, and 2. rapidly isolate the ground state again after inserting the current enabling an earlier Euclidean time fit, the associated reduction in the error bar through this process outweighs the increased cost in constructing the matrix of cross-correlators. In our implementation, due to the construction of the complete correlation matrix of three-point functions, we not only have access to the ground state, but also to the first n − 1 excited states, where n is the dimension of our operator basis. This has been utilised in [25] to access the axial charge of nucleon excitations. In principle, if one were solely interested in the ground state properties, one could use the optimised sources generated via the two-point correlation matrix as the input for the SST inversion, providing SST propagators that couple directly with the ground state. This reduces the cost from 2n inversions down to n+1. For this calculation the cost would be reduced from 8 to 5 inversions. Further reduction in cost is demonstrated through Fig. 4. It was found that access to ground state properties can be achieved with 3 levels of smearing, provided the smearings are chosen to span the space. Therefore, we could further reduce the cost to 3 + 1 = 4 inversions per configuration, only a factor of 2 above the minimum for what is equivalent to an order of magnitude improvement in the statistics.
In Table 3 we present a comparison of our result for g R A with results by other groups on similar ensembles. The consistency between our result and those of other groups is testament to the care taken by these collaborations to minimise systematic uncertainties.
A key issue in the calculation of any three-point function is how large must one make their source-sink separation to ensure that excited state contaminations are sufficiently suppressed [16]. There is a general consensus within the community that source-sink separations 1.0 fm will suffer from Table 3: Comparison of results for g R A on ensembles with similar volumes and values of m π to our calculation. For the CLS/Mainz group we have included results for the conventional ratio method (upper) and the summation method (lower). The asterisk indicates that these results include the correction of finite-volume effects and so will tend to sit slightly higher. excited state contaminations without fine-tuning the source and sink to isolate the state. Indeed our results highlight this systematic effect when using the conventional approach. Here the source-sink separation of ∼ 1.0 fm is too small and the extracted value for g A suffers from excited state effects as illustrated in Fig. 3. The underlying issue is that there is insufficient time to isolate the ground state prior to current insertion and again isolate the ground state before annihilation. Based on our earlier arguments regarding a more suitable current insertion time, we would expect an suitable sink time would be t s = 29, increasing the source-sink separation to ∼ 1.2 fm. This result is consistent with the source-sink separations used by the other groups in Table 3.
Using the variational approach, due to rapid onset of ground state dominance through ideal interpolators, we are able to use much smaller source-sink separations. For our variational results, ground state dominance after the current insertion occurs as early as t s = 23 resulting in a temporal separation between source and sink of only 0.64 fm. Thus, by applying the variational technique to fixed sink methods, one could consider source-sink separations ∼ 0.7 fm which would result in small statistical errors.

Conclusion
In this letter we have illustrated how the variational approach can be used to eliminate excited state effects from the calculation of the nucleon axial-vector coupling constant g A . These effects act to suppress lattice simulation results for g A . The use of optimised interpolators results in rapid ground state dominance allowing for earlier insertion of the current and earlier fit windows resulting in smaller statistical uncertainty. The key advantage to this approach is that once a suitable basis has been chosen, optimised sources are constructed automatically eliminating the need to tune smearing parameters and source-sink separations.
The method is general and would be ideally suited to calculations of form factors where the variational approach could be applied separately for each choice of source-sink momentum combination. Another quantity that has so far proved elusive for lattice calculations and could benefit from our approach is the quark momentum fraction, x , which is notorious for producing lattice results that are more than 50% larger than phenomenological determinations (see [36] for a review).
Future investigations will accurately calculate g A at a variety of quark masses and connect these results to Nature via finitevolume chiral effective field theory.