Controlling Excited-State Contributions with Distillation in Lattice QCD Calculations of Nucleon Isovector Charges $g_S^{u-d}$, $g_A^{u-d}$, $g_T^{u-d}$

We investigate the application of the distillation smearing approach, and the use of the variational method with an extended basis of operators facilitated by this approach, on the calculation of the nucleon isovector charges $g_S^{u-d}$, $g_A^{u-d}$, and $g_T^{u-d}$. We find that the better sampling of the lattice enabled through the use of distillation yields a substantial reduction in the statistical uncertainty in comparison with the use of alternative smearing methods, and furthermore, appears to offer better control over the contribution of excited-states compared to use of a single, local interpolating operator. The additional benefit arising through the use of the variational method in the distillation approach is less dramatic, but nevertheless significant given that it requires no additional Dirac inversions.


I. INTRODUCTION
The past decade has seen tremendous improvement in the ability of Lattice QCD (LQCD) calculations to provide results that can confront experiment. However, lattice computations of some key quantities remain at odds with experimental determinations, including the momentum fraction carried by partons in the nucleon, and, notably, the axial-vector charge, g u−d A , of the nucleon. These discrepancies are often attributed to finitevolume effects, and to the contribution of excited states to ground-state matrix elements. The calculation of the axial-vector charge has been a particular focus within the lattice community, and a dedicated effort to resolve these discrepancies has demonstrated success [1], using a method to control excited states inspired by the Feynman-Hellman theorem. In this paper we explore a different means of controlling excited states, through the use of a novel smearing method, "distillation" [2], and the variational method for the case of the nucleon charges g u−d A , g u−d S and g u−d T .
The matrix element of the flavor non-singlet axialvector current A a µ = ψγ µ γ 5 τ a 2 ψ, where ψ is the isospin doublet of u, d quark fields and τ the Pauli-spin matrices acting in isospin space, between nucleons of momenta p and p can be expressed in terms of the axial-vector and the induced pseudoscalar and tensor form factors where q µ = p µ − p µ the momentum transfer. At zero momentum transfer, the contribution of the pseudoscalar and tensor form factors to the matrix element of equation 1 vanish, and we are left with the definition of the isovector axial charge of the nucleon g u−d A = G A (0). The axialvector charge g A quantifies the n → pe − ν e coupling, the degree of low-energy chiral symmetry breaking in QCD, and even the difference between the u and d quark contributions to the proton spin g A = ∆u − ∆d. The breadth of different phenomena dependent on a knowledge of g A highlights the need for precise agreement between theoretical and experimental determinations. Moreover, the straightforward definition of g A serves as a useful proving ground for any new lattice algorithm purporting to calculate meaningful quantities in QCD.
Precise neutron β decay experiments have measured the nucleon isovector axial charge g u−d A = 1.2724±0.0023 [3]. Yet, historically, the bulk of lattice calculations have systematically under-determined g u−d A by roughly 10-15%. This has led to intense discussion in the community on the role of finite-volume effects, the lack of chiral symmetry in most lattice formulations of QCD, the influence of heavy quark flavors, the use of local currents with large discretization effects, and excitedstate contamination might have on the calculation of g u−d A and nucleon matrix elements in general. Calculations have been performed using the Highly Improved Staggered Quark (HISQ) [4] and Domain Wall (DW) [5] fermion actions which, despite the restoration of chiral symmetry, present challenges, namely, identifying quark flavors and the numerical cost compared to Wilson-type actions, respectively. Calculations have even been performed with O a 2 improved currents with novel noise reduction techniques [6,7], and others still with 2+1+1 flavor QCD accounting for several excited-states in requisite fits [8,9]. In all such cases the calculated value of g u−d A differs from experiment by ∼ 5 − 10%. Clearly standard techniques of calculating and subsequently fitting two-and three-point correlation functions to extract g u−d A routinely come up short. Recent work by [1] employs a methodology inspired by the Feynman-Hellman theorem to control excited-state effects by summing over all current insertion times, engendering extrapolation in a single Euclidean temporal variable rather than two -agreement was found to within 1% of experiment. We remark that a few 2-flavor lattice QCD calculations employing the Wilson plaquette/clover fermion actions have been performed for g u−d A whose results are within ∼ 1% of the experimental determination [10,11]. These studies calculated g u−d A for source-sink separations greater than 1 fm. Furthermore it was found g u−d A has a strong dependence on the spatial lattice volume V 3 . Taken together, suggesting the largest contributors to uncertainty in g u−d A lie in small lattice volumes and excited-state contamination.
In this paper, we investigate a means of overcoming one of the dominant systematic uncertainties in the calculation of nucleon charges, namely that arising from our inability to isolate the ground-state nucleon from its excitations. The variational method, an algorithm to improve overlap onto a desired state, was applied in [12][13][14] to a small basis of Jacobi-smeared nucleon interpolators of the form χ (x) = abc u a (x) Cγ 5 d b (x) u c (x), where it was found variationally improved interpolators can reduce excited-state effects. We will show that an amalgam of a suitable basis of interpolating operators, the use of the variational method with that basis, and an efficient means of computing the needed correlation functions through the application of "distillation", affords a powerful and computationally efficient method of taming excited-state contributions to those charges. Furthermore, we will show that "distillation" alone, enabling a momentum projection to be performed at each time slice in the two-and three-point functions, provides a significant increase in the statistical precision whilst being an effective smearing operator with only a single interpolating operator. This allows the reliable extraction of matrix elements at earlier source-sink separations when the nucleon signal-to-noise ratio is exponentially better. As lattice calculations proceed to ever more complicated quantities, exemplified by quasi-PDFs and pseudo-PDFs, the need for tamed excited-state effects is paramount.
In this work, we abstain from calculating and presenting renormalized isovector charges of the ground-state nucleon, and from performing continuum, finite-volume, and chiral extrapolations. Throughout this work we consider forward-scattering matrix elements of nucleons at rest. Although much of our motivating comments have centered around g u−d A , we extend our technique to include the scalar g u−d S and tensor g u−d T charges of the nucleon, whose precise determination will constrain BSM searches at the TeV scale and dark matter directdetection searches. A future study will explore forwardscattering matrix elements of moving states, and generalizations to scalar, axial, and tensor nucleon form factors. This paper is organized as follows. In section II we discuss excited-state contamination and present the computational methodology employed throughout this work. In particular, we review distillation as a low-mode approximation of the more standard Jacobi and Wuppertal smearing algorithms, and discuss the variational method as a means of improving overlap onto a desired state. Section III describes the lattice ensemble, the explicit construction of the interpolating fields, and then finishes with a discussion of the strategy used to extract the nucleon charges from our data. In section IV we first feature a comparison of the nucleon effective masses obtained using a Jacobi-smeared interpolator, a single "local" distilled interpolator, and different variationally improved interpolators from distinct bases of distilled interpolators. We then proceed to present determinations of effective masses of the nucleon via fits to our data and ultimately our extracted charges. We then conclude with a discussion of our results, a cost-benefit analysis of standard smearing techniques and distillation, implications for assorted studies in nucleon structure, and directions for future research.

II. COMPUTATIONAL METHODOLOGY
We begin this section with some definitions we will use throughout this work. Isovector charges g Γ of the nucleon are measured experimentally through the neutron to proton transitions p (p, s)| O ud Γ |n (p, s) , where O ud Γ = uΓd, or via proton and neutron charge differences. Imposing isospin symmetry, one can show where O u−d Γ = uΓu − dΓd is an external current. The isovector charges of the nucleon are thus defined as where we normalize the nucleon spinors in Euclidean space according to Indeed the calculation of isovector quantities is computationally less demanding than the calculation of isoscalar and flavor-diagonal quantities, in which the cancellation of disconnected quark lines in the isospin limit does not occur.

A. Excited-State Contamination
The calculation of hadronic matrix elements within LQCD requires the construction of operators that attempt to interpolate a given hadronic state from the vacuum, with minimal overlap with neighboring states. As the precise wavefunction of any hadronic state is not known, any operator construction of a desired continuum J P C is merely a "best guess" and necessarily overlaps with other hadronic states of the same quantum numbers -most notably, excited-states and multi-particle states. This problem is compounded with the explicit breaking of Lorentz symmetry, in which continuum operators residing in different irreducible representations (irreps) of the Lorentz group can mix under subduction to the same lattice irrep thereby increasing the number of states contributing within a given lattice irrep.
Explicitly, consider a two-point correlation function projected to zero momentum where O is an interpolating operator for the state of interest. Performing a spectral decomposition, this can be expressed as where the sum is over all states, of energy E n , that can be created with the operator O. In order to extract reliably the ground-state mass M 0 ≡ E 0 , one can study the large time behavior of the two-point correlator, wherein contributions of excited states are negligible, or make a judicial choice of interpolating operator such that the overlap factors Z n = 0| O |n , for states n > 0, are greatly suppressed relative to 0| O |n = 0 , and thereby determine E 0 at small temporal separations. Given that lattice calculations of baryon properties are constrained by an exponentially increasing noise-to-signal ratio with increasing Euclidean time, the latter approach is far preferable. This issue is further compounded when considering matrix elements of an external current J , where the suppression of excited-states is needed both between the source interpolator and the current, and the current and the sink interpolator. Thus there is strong encouragement to develop operators that couple predominantly to the ground-state, and only weakly to the excited-states.

B. Smearing
Interpolating fields constructed of point-like quark and gluonic fields couple to hadronic states at all energy scales. "Smearing" is a well-established technique to increase the overlap of interpolators with the low-lying states of the spectrum (i.e. confinement scale physics), and to reduce the contribution of the high-energy modes to correlation functions, through the use of spatially extended operators of hadronic size. Specifically, we can replace the quark fields ψ( x, t) occurring in the path integral with spatially extended quark fields where S[U ]( x, y) is a gauge-covariant scalar "smearing" kernel that is functionally dependent on the underlying gauge fields U on some time slice t.
A frequently utilized smearing operator is J σ,nσ (t) = 1 + σ∇ 2 (t) nσ nσ defining the Jacobi-smearing method [15], where ∇ 2 is a gauge-covariant discretization of the Laplacian, σ is the smearing "width", and n σ represents the number of applications of the Laplacian. In the large n σ limit, thus approaching a Gaussian of width σ, characteristic of the size of a hadron.

C. Distillation & Operator Construction
Distillation [2] is a low-rank approximation to a gaugecovariant smearing kernel. Specializing to the case of the Laplacian, we begin by seeking solutions to −∇ 2 (t) ξ k (t) = λ k (t) ξ k (t). Ordering solutions by the magnitude of the eigenvalues λ k (t), the distillation operator is constructed as the outer-product of two eigenvectors on a given time slice where the color indices {a, b} are made explicit. The distillation operator is then applied to each quark or antiquark field both at the source and at the sink.
Distillation affords several advantages over Jacobi-like smearing methods. Distillation, firstly, factorizes the construction of the interpolating or current operators from quark propagation. Secondly, distillation allows operators not only at the sink location but also at the source to incorporate more elaborate spatial structure, and in particular derivatives, without additional inversions of the Dirac operator. Finally, distillation enables momentum projection to be performed both at the source and sink time slices, thus providing a more thorough sampling of a given lattice configuration. In this case, equation 3 becomes These features have been key to the precise determination of the many energy levels in the QCD spectrum needed for the study of resonances in lattice calculations. However, these features are also of advantage in studies of hadron structure, by facilitating more varied interpolator constructions, and by enabling momentum projection at each time slice in a correlation function.
Specializing to the case of baryons, we will construct our interpolating operators O i following the procedure of ref. [16,17] where D 1,2,3 are spatial operators constructed from covariant derivatives, introduced to probe the radial/angular structure of the nucleon wavefunction, and S αβγ are subduction matrices that project a state of definite continuum spin into irreps of the hypercubic lattice (explicit construction given in Section III). The building blocks of the two-point and three-point correlation functions employing these interpolating operators are where it should be noted that the inversion of the lattice Dirac operator against a smeared point-source in standard techniques, is replaced in distillation by inversion of the lattice Dirac operator against the k th eigenvector of time slice t. The principle disadvantage of the distillation method is that the number of distillation eigenvectors N needed to construct a correlation function of the same resolution is expected to scale as the lattice spatial volume V 3 . Since the evaluation of the Wick contractions for mesons and baryons scales as N 3 and N 4 , respectively, the volume-scaling cost is potentially severe.

D. Variational Analysis
The variational method is an often employed technique in lattice spectroscopy calculations that seeks to disentangle the contributions of eigenstates to the two-point correlation of two interpolating operators. The variational method begins by constructing a matrix of correlation functions where O i and O j belong to some basis B of appropriately constructed interpolators with identical quantum numbers. In practice, these interpolators transform with definite symmetries in quark flavor, derivative structure, and Dirac structure -where the Dirac structure of the operator is encoded in the subduction of a continuum operator into lattice irreps. The actual variational method begins by considering the system of generalized eigenvalue equations with n ∈ {1, · · · , dim (B)} and t > t 0 . It can be shown that at large times t this system of equations is described by the "principal correlators" λ n (t, where the trivial solution v n = 0 is avoided by imposing the normalization condition v † n' C (t 0 ) v n = δ n'n . Equation 9 is solved independently for each t > t 0 , with λ n (t 0 , t 0 ) = 1 by construction. It is possible that in the process of solving Equation 9, what appears to be the m th eigenvector on time slice t m is not the m th eigenvector on time slice t m+1 . To remedy this, we follow the procedure of [16] by associating states on different time slices using the relative similarity between their associated eigenvectors. This is accomplished by selecting some reference time slice t ref and is applied to the principal correlator to obtain the masses M n , M n , and amplitude A n . The choice for t 0 is made by attempting to reconstruct the original correlation matrix C ij (t) using the masses M n , M n and the extracted overlap factors Z n i [18]. The degree of agreement for t > t 0 dictates the choice for t 0 . Solving the variational method introduces a slight time-dependence in the overlap factors. The overlap factors are thus evaluated on a time slice t Z > t 0 that minimizes differences between original/reconstructed correlators [18]. The resulting eigenvectors v n for each t > t 0 yield the optimal linear combination of O i ∈ B to project the n th -state from the vacuum The variational method may equally be applied to correlation matrices comprised of three-point functions of different interpolating fields. This would necessarily produce a better determined O opt n , but due to the high Wick contraction cost of using distillation we elected to perform the variational method on a correlation matrix of two-point functions, thereby fixing O opt n for later use.

III. COMPUTATIONAL DETAILS
Our analysis considers a 350 configuration ensemble of 2+1 flavor QCD using the clover-Wilson fermion action, where the associated gauge links are smeared by one application of the stout smearing [19] algorithm. This smearing yields a tadpole-improved tree-level clover coefficient that is nearly identical to the corresponding non-perturbative determination. The reader is referred to [12,13] for a discussion regarding the gauge action used to generate this ensemble. Calculations were performed on a 32 3 × 64 lattice with periodic (spatial) and anti-periodic (temporal) boundary conditions and an inverse coupling of β = 6.3. In the three flavor theory, the strange quark mass was set by requiring the ratio 2M 2 K + − M 2 π + /M Ω − to assume its physical value.
This lattice ensemble was found to have a lattice spacing of a = 0.09840(4) fm via the Wilson-flow scale w 0 , and with am π ∼ 0.176803 yielding am π L ∼ 5.658 and m π = 356 MeV. We explore the efficacy of four different types of operators used to interpolate ground-state nucleons from the vacuum, with particular consideration given to distillation. In this work we only study zero-momentum nucleons, polarized along the z-axis. A future work will study nucleons with non-zero momentum.

Jacobi smeared Interpolator
Prior to construction of our quark sources comprising our selected interpolators, the background gauge links are smeared via 20 applications of stout smearing with smearing parameter: ρ ij = 0.08 and ρ µ4 = ρ 4µ = 0, where ρ µν quantifies the weight given to staple links aligned in the µν-plane when constructing the smeared links. Such gauge smearing is essential for reducing noise present in the resultant correlation functions due to source fluctuations. Before inverting the Dirac operator against the smeared sources, we apply a single interation of stout smearing with ρ = 0.125 to the gauge links -thereby avoiding potentially small eigenvalues in the inversion.
As a benchmark with which to compare to distillation, we begin with the simplest nucleon interpolator consistent with the nucleon J P C where u, d are the two flavors of (degenerate) light quarks, {a, b, c} color indices, C = γ 2 γ 4 the charge conjugation matrix, and a suppressed Dirac index. The nonrelativistic projector (1 ± γ 4 ) /2 is included in the operator construction to reduce the noise-to-signal ratio in forward (backward) propagating states. To make contact with previous work (e.g. [12,13]) N (x) is smeared with 60 hits of Jacobi smearing, with width σ = 5.0. We refer the reader to [12,13] for an extensive analysis of the effect different Jacobi smearing parameters has on the determination of nucleon isovector charges. Herein we refer to the Jacobi smeared interpolator as "Jacobi-SS". Correlation functions that employ N as the source/sink interpolator are constructed via application of appropriate projection operators where P 2pt = (1 + γ 4 ) /2 is used to project onto the forward-propagating positive-parity nucleon, and where P 3pt = P 2pt (1 + iγ 5 γ 3 ) is used for the corresponding connected insertions. The spatial sums serve to project each correlation function to zero momentum. A standard spectral decomposition demonstrates that the Dirac structure of O u−d Γ must be 1, γ 4 , γ 3 γ 5 , γ 1 γ 2 to extract the scalar, vector, axial, and tensor charges, respectively. Lastly, the sequential source method is implemented to calculate C 3pt , thereby minimizing the number of distinct inversions of the Dirac operator.

Distilled Interpolators
We use a distillation space of rank 64, from which the perambulators and solution vectors are constructed. The distillation space on each time slice is calculated only after 10 iterations of stout smearing is applied to the gauge links with smearing factor ρ ij = 0.08 and ρ µ4 = ρ 4µ = 0.
When expressed in a form exposing the permutational symmetry of the flavor (F), spatial (D) and Dirac (S) structures, our distilled interpolators take the form where P (· · · ) expresses the symmetric (S), mixedsymmetric (M), and anti-symmetric (A) character of the given structure. Explicitly our employed distilled interpolators are where the superscript J P indicates the overall spin-parity quantum numbers of the interpolator [16]. For brevity, we have expressed the interpolators in a compact spectroscopic notation 2S+1 L P J P , where S represents the Dirac spin, L the angular momentum induced by any derivatives, P the permutational symmetry of the derivatives, and J P the total angular momentum and parity of the interpolator. The first distilled interpolator we consider is the 2 S S 1 2 + interpolator, which is the closest nonrelativistic analogue of the standard nucleon interpolator given in Eq. 12.
Our first application of the variational method is to a basis of three distilled interpolators where we note 4 P M  [17] where it was found these hybrid interpolators, in addition to 2 S S 1 2 + , had predominant overlap onto the ground-state nucleon. The variational method applied to B 3 leads to a variationally improved interpolator that we define to beP 3 . The final interpolator we consider is found by expanding the basis B 3 to include distilled interpolators that probe the radial/orbital structure of the nucleon (see above) We refer to this variationally improved interpolator aŝ P 7 . The superficially redundant 2 S S 1 2 + and 2 S S 1 2 + interpolators, with the same spectroscopic notation but differing derivative constructions, correspond to different radial extents of the interpolator. As outlined in Section II D, the construction of a variationally-improved interpolator relies on careful selection of a reference time slice t 0 , and the time slice t Z at which to evaluate the overlap factors. Dividing out the ground-state time dependence, if a single exponential were to contribute to the principal correlator, a plateau of unity should be expected. Based on the applied fits of equation 10, a good determination of the ground-state within our basis of interpolators is thus indicated by a plateau in the re-scaled correlator around unity -A 0 1 and ∆m n 0 = m n − m 0 1. For the two variationallyimproved interpolators we consider, we found led to ideal reconstruction of the original correlator.

A. Matrix Element Extraction
The extraction of matrix elements on the lattice typically proceeds by calculation of some 3-point correlator, over a range of source-sink interpolator separations, and where an external current is inserted for all intermediate times. Under the presumption of no excited-state contamination (i.e. in the limits 0 τ t sep ), the 3-point correlator is then divided per ensemble average by a two point correlator with the same source/sink interpolators of the same source-sink separation. This division removes overlap factors, masses, and exponential source/sink time dependence from the matrix element signal. A plateau in this ratio, for fixed interpolator separations and varying current insertion times, should be the desired matrix element. However, as we are interested in ground-state matrix elements J 00 , the plateau necessarily contains contributions from matrix elements of all excited-states.

Correlator Behavior
Our two-point correlation function using Jacobismeared interpolators is defined by Inserting the non-relativistic projector, and performing a spectral decomposition exposes the competing contributions of all states in the spectrum: where the sum is only over eigenstates with quantum numbers of N . To quantify and control excited-state contributions to C 2pt , we elect to perform a 2-state fit of the form for each of the four distinct interpolators we consider. The factoring of the ground-state time dependence aids in stabilizing our fits and in the determination of our extracted charges, as explained later. In this manner we obtain determinations of M 0 , M 1 , |Z 0 | 2 and |Z 1 | 2 , while simultaneously quantifying the efficacy of each interpolator to separate the ground-state from its excitations viz. ∆m = M 1 − M 0 . The correlator behavior when using distilled interpolators is identical to that above, except for the addition of an overall volume factor V 3 due to the momentum projection at the source implied by eqn. 6.
A zero momentum projected three-point correlation function using Jacobi-smeared interpolators is given by where τ is the insertion time slice, restricted to the timeordering 0 < τ < t sep , and J the external current. An analogous spectral decomposition with the appropriate projector leads to Again retaining the lowest two energy eigenstates in the sum, we have where J n n = n , s | J |n, s , with n , n ∈ N. By isolating the current insertion time dependence of the threepoint correlator, it becomes clear that calculation of a three-point correlation function for a single source-sink separation t sep is insufficient to reliably extract the matrix elements J 00 and J 11 . As Z n are real and J 01 = J 10 for zero-momentum states, the above can be reorganized into We then apply a 2-state fit to the three-point correlation functions, where τ is the current insertion time and ∆m = M 1 −M 0 . We note that we retain M 0 and M 1 as parameters in our fit, rather than the difference ∆m. As with the functional form employed to fit the two-point correlation functions, we factor the ground-state time dependence from the functional form of the three-point fits. This factoring makes manifest the desired ground-state matrix element in the limits 0 τ t sep . The correlator behavior when using distilled interpolators is again identical to that above, except for the addition of an overall volume factor V 3 to Equation 20.
Determining the precise relationship between the fitted parameters to extract the ground-state matrix element J 00 , requires a more detailed look at our interpolators. Although the constructed distilled interpolators contain no free spinor indices, our use of positive-parity nucleons polarized along the z-direction can be viewed as the standard application of projectors on the nucleon interpolating field of Equation 12. At zero-momentum we again have, Proceeding with the spectral decomposition of C 2pt we have, From which it is easily shown that the constant coefficients of the 2-state fit applied our two-point correlators are of the form |a| 2 = 2 |Z 0 | 2 and |b| 2 = 2 |Z 1 | 2 . The spectral decomposition of C 3pt proceeds analogously, It then follows A = 2g 00 Γ |Z 0 | 2 and B = 2g 11 Γ |Z 1 | 2 . We can then extract the ground-state matrix element via To extract the masses, overlap factors and matrix elements from our data, we perform simultaneous fits to the two-point and three-point correlation functions according to Equations 19 and 22, respectively.

IV. RESULTS
We compute the two-point functions averaged over three source positions. The extraction for the matrix elements is obtained from three-point functions for t sep ∈ {8, 12, 16}, and τ ∈ [0, t sep − 1]. Rather than view the calculated two-point correlation functions directly, we judge the quality of the two-point correlators for each interpolator by plotting the effective mass as a function of the source-sink separation t. Calculating each two-point correlator, averaged over three different source positions, leads to the effective masses seen in Figure 1. The lack of a plateau in the effective mass of the Jacobi-SS interpolator until t ∼ 10 is indicative of excited-state contamination for source-sink separations of 1 fm. Use of the 2 S S 1 2 + distilled interpolator leads to an earlier onset of a plateau in the effective mass, with statistical uncertainty that is at least 50% smaller than that of the Jacobi-SS interpolator. This plateau is seen to begin for t ∼ 6, or ∼ 0.6 fm. It is also worth noting that the expected exponentially increasing noise in the nucleon effective mass is substantially suppressed at larger source-sink separations, when compared to the Jacobi-SS interpolator. Use of a variationally improved interpolator derived from bases of distilled interpolators, (P 3 orP 7 ), leads to a more rapid exponential decay of excited states at early Euclidean times. The effective mass induced by theP 3 interpolator exhibits a plateau that has nearly the same statistical precision as that of the 2 S S 1 2 + interpolator, while the excited-states are seen to decay more rapidly for t < 5. Expanding the basis of distilled interpolators, thê P 7 interpolator leads to an even more rapid decay of excited states for t < 5, consistent with a variationally improved interpolator receiving excited-state contributions from states higher in energy than those within the basis. As with the 2 S S 1 2 + andP 3 interpolators, the effective mass ofP 7 begins around t = 6, however the plateau is noticeably lower than the former. In general, the statistical precision of all distilled interpolators appears to be comparable, except for large Euclidean times wherein the determination of the variationally improved interpolators becomes increasingly uncertain. We attribute this increased uncertainty to stem from the variationally improved interpolators being relatively unconstrained at large source-sink separations, where elements of the correlation matrix (Eq. 8) are themselves dominated by noise.

A. Two-State Analysis
To quantify the discussion above and to guide our future simultaneous fits to the two-and three-point functions, we explore the efficacy of one-and two-state fits have in describing the calculated effective masses. All fits are restricted such that 2 ≤ t fit ≤ 16 thereby avoiding contact terms in the clover Wilson action and a collective fluctuation of the effective masses for t 16. The results of these fits are collected in Tables I and II. Although single-state fits make explicit the improvements gained by using distillation over standard Jacobismearing, notably a ∼ 3% difference in the determination of the ground-state mass, the improvements at this stage hardly seem worth the cost of constructing the required distillation basis. By including a second state in  the fits performed to the two-point correlators, the gains produced by distillation are quite encouraging. (blue) distilled interpolator shows considerable improvement over the Jacobi-SS interpolator, while applying the GEVP to a basis of distilled interpolators of order 3 (green) and of order 7 (red) show further improvement.
The inclusion of a second state in fits to the two-point correlators yields slightly smaller determinations for the ground-state mass and overlap factors when using the Ja-cobi, 2 S S Evidently distillation and the variational method lead to greater elimination of excited-state contributions to the two-point correlators, where the mass gap is ∼ 44%, ∼ 58%, and ∼ 80% larger for the 2 S S 1 2 + ,P 3 , andP 7 interpolators, respectively, when compared to the Jacobi-SS interpolator. The compounded improvements of the variational method applied to our bases of distilled interpolators is entirely consistent with [20].
In this section we seek a more reliable determination of the masses, overlap factors, and likewise nucleon matrix elements, by performing simultaneous fits to three-point and two-point correlators with interpolator O and current structure Γ. We fix the window over which the twopoint correlators are fit to be t fit 2pt ∈ [2,16], while several three-point fit windows are studied. The three-point fit windows are identified by τ buff -that is for a given t sep , To illustrate the extracted isovector charges and to quantify the degree of excited-state contamination, we plot an effective charge defined as where C 3pt Γ (t sep , τ ) is the three-point correlation function for a given source-sink separation and current insertion time, and C 2pt fit (t sep ) is the corresponding best fit applied to the two-point function of the same interpolators and source-sink separation. Errors on the effective charges are computed via a simultaneous jackknife re-sampling of the two-point fit and three-point correlator, to account for the correlations between fitted parameters. Shown together with the effective charges are our extracted values for the isovector charges using the methodology described above. Error (gray) bands of the extracted charges (black line) are determined through the ratio of fit parameters A/ |a| 2 , likewise accounting for parameter covariance. We highlight that nearly the same masses for the ground-and first-excited states are found from the 2-state fits to the two-point correlators and from the simultaneous fits to the two-and three-point correlators.
Coupling of the isovector scalar current S 3 = q τ 3 2 q to the Jacobi-SS interpolator N α produces an effective matrix element (denoted here as M eff S ) that is determined to within no less than ∼ 10% uncertainty. Although M eff S is symmetric about the midpoint τ − t sep /2 for t sep = 8, indicating equal excited-state contamination on the source/sink side of the insertion, M eff S exhibits antisymmetric behavior for t sep /a = 12 and is largely statistical noise for source-sink separations greater than 1.5 fm.
Immediately noticeable with the use of a single distilled interpolator ( 2 S S are not necessarily on the curve reflects that the fit is largely constrained by data with t sep 12. Introduction of two hybrid interpolators to construct P 3 appears to return the t sep = 8, 12 determinations of M eff S to a form structurally similar to that found for the Jacobi-SS interpolator. It is curious to note that the t sep = 16 determination of M eff S is nearly identical to the corresponding determination with 2 S S  2 q is the insertion that is most sensitive to excited-state contamination and finite volume effects. As with g u−d S , determinations of the effective matrix element (denoted here as M eff A ) for different t sep using the Jacobi-SS interpolator are plagued with poor statistics. Clearly spatial Gaussian smearing of N α is not sufficient to suppress excited-state effects.
Employing the local distilled interpolator 2 S S 1 2 + yet again leads to dramatic reduction in statistical uncertainties in determinations M eff A , and an increase in the extracted value of g u−d A by ∼ 7%. Notably broad plateaus in M eff A can likewise be seen for t sep = 8, 12 that are consistent within error, thereby reducing the minimal sourcesink separation ( 1 fm) to reliably extract g u−d A . The behavior of theP 3 interpolator is in line with that found for g u−d S -reduced agreement between t sep = 8, 12 determinations, and fits most heavily constrained by data at smaller source-sink separations. Curiously, the slight oscillation seen in the t sep = 16 2 S S 1 2 A is amplified following application of the variational method. Prior to τ − t sep < −2, the t sep = 16 effective matrix element seems to fall into agreement with the t sep = 12 determination. We speculate the abrupt decrease in the t sep = 16 effective matrix element to be caused by neglect of backward-propagating positive-parity states in the corresponding two-point functions.
TheP 7 interpolator perpetuates the double-plateau feature in the t sep = 16 determination of M eff A . Remarkably, however, application of the variational method to this enlarged basis of distilled interpolators demonstrates absolute agreement between the t sep = 8, 12 determina- The isovector tensor current T 3 µν = q i 2 [γ µ , γ ν ] τ 3 2 q is among one of the best determined quantities in the nucleon, particularly at zero virtuality defining g u−d T . This facet is supported by models that show excitations of excited-states in the nucleon are suppressed when coupling the tensor current T 3 µν [21]. As with previous charges, the 2 S S 1 2 + interpolator leads to a dramatic reduction in statistical uncertainty of the effective matrix elements (denoted here as M eff T ) when compared to Jacobi-SS. By t sep ∼ 1.18 fm, a definite plateau is present in M eff T for several insertion times τ ; this same plateau is shared with the t sep = 16 determination. The absolute agreement of the central values and error of the t sep = 12, 16 determinations suggest that when using distillation, a source-sink separation of roughly 1fm is sufficient to reliably extract g u−d T . The variationally improvedP 3 expands upon the improvements seen with the 2 S S the t sep = 12, 16 determinations depict an even broader plateau that are again entirely consistent within error. These enlarged plateaus for numerous insertion times τ follow naturally from a determination of ∆m that is greater than that determined with 2 S S 1 2 + .
Proceeding to theP 7 interpolator continues the trends seen with other distilled interpolators, wherein the statistical uncertainties are greatly reduced and the plateaus in the effective matrix elements are seen to become ever larger. Remarkably, the t sep = 16 determination usinĝ P 7 resembles that of the vector charge -a matrix element that is constant in τ within minor statistical fluctuations. We emphasize here that the effective matrix element for τ −t sep /2 = −8 is coincident with the source interpolator and thus should not be considered a reflection of g u−d T .
Recalling the functional form of the fit applied to the three-point correlator, Eq. 22, the coefficient B encodes the first excited-state matrix element and C captures the ground to first excited-state transition matrix element. As evident from Table V, the determinations of B for each interpolator are generally consistent with zero, indeed supporting the notion that excited-state contributions to g u−d T are greatly suppressed. On the other hand, the largest source of contamination appears to stem from the transition matrix element contained in C.
As the isovector vector current V 3 µ = qγ µ τ 3 2 q simply counts the number of quarks within the nucleon, and all its excitations, there is no surprise in the extraction of g u−d V for each interpolator. We have included the nucleon vector charge here as a useful sanity check, thus ensuring such an observable that is independent of state and τ is indeed recovered. The continuum requirement trivially sets the vector current renormalization to be Z V ∼ 0.856, consistent with determinations on finer lattice ensembles with lighter m π [12,13]. As for the previous charges, distillation alone affords considerably improved statistics over Jacobi-SS. That said, application of the variational method to a basis of distilled interpolators appears to produce a curious "spreading" in the determined effective vector matrix element.

C. Numerical Cost
Given the substantial benefits incurred by use of distillation and the variational method in such calculations of hadronic structure, it is worth while to pause and consider the numerical cost of doing so. We highlight that a true one-to-one comparison of Distillation to standard techniques is not entirely possible, as distillation is markedly distinct from traditional spatial smearing techniques via sampling of entire time slices.
The calculation of point-to-all propagators in standard lattice calculations proceeds by inverting the chosen discretization of the Dirac operator against a point source in coordinate, spinor, and color space This operation requires 12 distinct inversions of the Dirac operator, one for each spinor and color index {α, a}. In our case of the nucleon, with degenerate u and d-quarks, this captures the propagation of the nucleon from some source point to all other points on the lattice. Implementing the sequential-source method, as we did for the Jacobi-SS interpolator, reduces the number of required inversions by combining the computed point-to-all propagator with the sink interpolator and using this object as a source for further inversions (deemed sequential propagators). As the u-quarks can be combined into one sequential source and the d-quarks another, the computation of a three-point function in the sequential-source framework requires two additional inversions against a color and spinor point source. As the sequential propagators must be recomputed for each new source-sink separation t sep , we arrive at where we note that this is now proportional to N 2 ops . In the case of distillation, the inversion of the Dirac operator against a point source is replaced with inversion against an eigenvector on some time slice t: As the eigenvectors are determined by solution of −∇ 2 (t) ξ (k) (t) = λ (k) (t) ξ (k) (t) given some gauge covariant discretization of ∇ 2 , the calculation of the k th solution vector requires 3 inversions of the Dirac operator, one for each (suppressed) color index. In practice, for a given t sep , we calculate the solution vectors forward (backward) from the source (sink), where the solution vectors from the source are used in the two-point and three-point calculations. Thus for a single distilled interpolator the total number of inversions becomes where N eigs is the dimension of the distillation space.
With N eigs = 64 in our case, this cost at first seems excessive. However, once these solution vectors have been computed, any number of interpolating fields, variationally-optimized or otherwise, can be correlated without additional cost. We note that we have not taken into account the cost of the Wick comparisons when using distillation in this study, and a future work will include the stability of our extracted matrix elements as the rank of the distillation space is reduced.

D. Conclusions
We have investigated the use of distillation, and an extended basis of interpolators, in the calculation of the scalar, axial and tensor isovector charges of the nucleon, and made comparisons with a calculation on the same ensemble using the conventional Jacobi-smearing method of a single smearing radius. We find that distillation affords a considerable improvement in the statistical quality of the data when compared with calculations using Jacobi smearing. We attribute this improvement to be a consequence of momentum projection performed at both source and sink in the case of a two-point function, and at source, sink and current in the case of a three-point function. More surprisingly, even the use of a single, local distilled interpolating operator results in a suppression of the contribution of excited-states relative to that of the ground-state in both two-point and three-point functions.
For our variational analysis, we employed a basis of operators comprising the non-relativistic operators that can be constructed with up to two covariant derivatives, together with so-called "hybrid" operators where the gluons play a manifest role. Whilst the improvement was not as dramatic as that between a single Jacobi-smeared and a single distillation-smeared operator, the use of the variational method and the extended basis provided more consistency and fidelity in the matrix elements for different source-sink separations. Furthermore, the extended basis can be introduced without further propagator calculations, in contrast to the case of Jacobi smearing where the use of the variational method demands a considerable increase in the number of propagators to be computed.
The next step in our program is to extend our investigation from matrix elements between nucleons at rest to those in motion, and from forward matrix elements to off-forward matrix elements. The former is key to the efficient application of the quasi-PDF, pseudo-PDF and current-current correlator methods to the calculation of parton distribution functions in the nucleon, whilst the latter is important for form factors at high momenta, and off-forward matrix elements such as generalized parton distributions. Nonetheless, the work presented here clearly demonstrates the efficacy of distillation as a means both of decreasing the statistical un-certainty, and of reducing excited-state contributions, in calculations of nucleon properties.

Acknowledgments
Calculations were performed using the Chroma [22] software suite on the computing clusters at Jefferson Lab. We are grateful to Jo Dudek and Stefan Meinel for the use of their fitting codes, and to Robert Edwards for invaluable discussions on the feasibility of these calculations. CE extends thanks to Balint Joó for invaluable discussions on building software on the varied machine architectures on the Jefferson Lab clusters. This material is based upon work supported by the U.S. Department