Excited-state effects in nucleon structure on the lattice using hybrid interpolators

It would be very useful to find a way of reducing excited-state effects in lattice QCD calculations of nucleon structure that has a low computational cost. We explore the use of hybrid interpolators, which contain a nontrivial gluonic excitation, in a variational basis together with the standard interpolator with tuned smearing width. Using the clover discretization of the field strength tensor, a calculation using a fixed linear combination of standard and hybrid interpolators can be done using the same number of quark propagators as a standard calculation, making this a cost-effective option. We find that such an interpolator, optimized by solving a generalized eigenvalue problem, reduces excited-state contributions in two-point correlators. However, the effect in three-point correlators, which are needed for computing nucleon matrix elements, is mixed: for some matrix elements such as the tensor charge, excited-state effects are suppressed, whereas for others such as the axial charge, they are enhanced. The results illustrate that the variational method is not guaranteed to reduce the net contribution from excited states except in its asymptotic regime, and suggest that it may be important to use a large basis of interpolators capable of isolating all of the relevant low-lying states.


I. INTRODUCTION
One of the most challenging sources of systematic uncertainty faced by lattice QCD calculations of nucleon structure is excited-state contamination: the failure to isolate the ground-state nucleon from the tower of higher-energy states to which the interpolating operator can couple. Although the unwanted excited states can be exponentially suppressed by Euclidean time evolution, this is hindered by an exponentially decaying signal-to-noise ratio [1] that makes it impractical to evolve long enough in Euclidean time.
The variational method [2][3][4] provides a way of improving the interpolating operator such that the lowest-lying excited states can be systematically removed. Variational approaches have been used to study nucleon structure in Refs. [5][6][7][8][9], which used bases of interpolators with different smearing widths and different site-local spin structures; Refs. [10][11][12], which used bases with the standard interpolator evolved by different Euclidean time intervals, i.e., the generalized pencil-of-function method [13,14]; and Ref. [15], which used the distillation method to enable the use of interpolators with a variety of local structures including covariant derivatives. In these cases, the variational setup was more computationally expensive than a standard calculation because of the need for additional quark propagators with different smeared sources or (for time-evolved operators) additional source-sink separations. 1 In this paper, we present a study of a variational setup that requires the same number of quark propagators as a standard calculation, for a fixed choice of optimized interpolator. This is accomplished by supplementing the standard interpolator with hybrid ones [16] that contain a gluonic excitation. 2 In Ref. [16], it was found that the latter have the next-largest overlaps onto the ground state, after the standard interpolator, but overlap much more strongly than the standard interpolator with certain high-lying excited states. The use of hybrid interpolators presents the possibility of an improvement over the standard approach at low computational cost.
It should be stressed that this study, along with all previous ones in the nucleon sector, 3 uses only local interpolating operators, which are poor at isolating multiparticle states. In practice, the true spectrum that includes Nπ-and Nππ-like states is not identified, meaning that the calculation is not in a regime where the variational method has been proven to improve the isolation of the ground state. Therefore, the question of whether one interpolator is an improvement over another is an empirical one, to be decided by examining excited-state contributions in estimators for a variety of observables.
This paper is organized as follows. Section II discusses our lattice setup, the basis of interpolating operators, and two tuning runs. Results for two-point correlators are discussed in Sec. III A, forward matrix elements are shown in Sec. III B, and form factors are presented in Sec. III C. Our conclusions are given in Sec. IV.

II. LATTICE SETUP
As this is an exploratory calculation, we use a single lattice ensemble with a coarse lattice spacing at a heavier-thanphysical pion mass and a relatively small box size; its parameters are summarized in Table I. This has 2 þ 1 flavors of tree-level improved Wilson-clover quarks coupled to the gauge links via two levels of HEX smearing [18].
Aside from the interpolating operator, the methodology used here for computing nucleon matrix elements and form factors is unchanged from previous work such as Ref. [11]; the reader is referred to that earlier work for details. Our focus is on seeing whether excited-state contamination can be reduced, and therefore we use three relatively short source-sink separations, T, ranging from 0.70 to 1.16 fm. We use two methods for determining matrix elements: the ratio method-for which the asymptotically leading excited-state contributions decay as e −ΔET=2 , where ΔE is the energy gap to the lowest excited state-and the summation method, for which they decay as Te −ΔET .
Given a set of N interpolating operators fχ i g, one would like to find a linear combination χ ¼ v i χ i that has a reduced coupling to excited states. The standard approach is to compute a matrix of two-point correlators, and then solve a generalized eigenvalue problem (GEVP) for some choice of ðt 1 ; t 2 Þ. It has been shown [4] that by suitably increasing t 1 and t 2 to remove contributions from higher excited states in the determination of v, one can define improved estimators for the ground-state energy and matrix elements, for which the leading excited-state effect depends on the energy gap to state N þ 1 rather than the second (i.e., first excited) state. However, it is known that for light pion masses and large volumes, the number of low-lying excited states with the quantum numbers of the nucleon is large due to the presence of multiparticle (Nπ, Nππ, etc.) states [19][20][21][22][23][24][25][26]. Removing the effects of these states would require that the basis include at least one operator for each state. In addition, it has been found in meson spectroscopy calculations that nonlocal multiparticle interpolators must be included in order to correctly identify the multiparticle spectrum: see, e.g., Ref. [27]. For nucleons, the need for nonlocal operators is also supported by chiral perturbation theory, which predicts at leading order that the ratio of couplings for single nucleon and nucleon-pion states is the same for all local operators [21]. This defeats the diagonalization procedure of Eq.
(2) such that Nπ states cannot be removed. The challenge of systematically removing all contributions from the lowest-lying excited states will be left to future work. Instead, we hope to find an improved local operator that can be used in existing software with minimal modifications and with little additional computational cost. Our standard operator is where P þ ¼ ð1 þ γ 4 Þ=2 is a positive parity (nonrelativistic) projector andq is a smeared quark field. When used in a two-point or three-point correlator with a polarization matrix that includes a factor of P þ , the projector is applied to all three quarks. This allows for computational cost savings in the quark propagators used for constructing correlators: only half of the propagator solves are required [28]. Of the three possible site-local nucleon operators, χ 1 is the only one that can be constructed using positiveparity-projected quark fields (see e.g., Appendixes B and C of Ref. [29]). We also consider hybrid operators, introduced in Ref. [16], which include an insertion of the chromomagnetic field B i ¼ − 1 2 ϵ ijk F jk and are interpreted as having a nontrivial gluonic excitation. Using the clover discretization of F μν [30], no additional quark propagators are needed for constructing two-point correlators using hybrid operators. Two different nucleon operators exist that use positive-parity-projected quark fields: where These differ in the spin of the three quarks: for χ 2 it is 1 2 and for χ 3 it is 3 2 . In both cases, this is combined with the chromomagnetic field to produce an overall spin of 1 2 . Our production strategy begins with two tuning runs where only two-point correlators are computed. The first uses only χ 1 and serves to select the quark smearing parameters. The second serves for determining the coefficients v i of an optimized operator χ opt ¼ v Ã i χ i . These are followed by a production run with higher statistics in which both two-point and three-point correlators are computed. The three-point correlators are computed using the standard operator χ 1 and the linear combination χ opt , in both cases keeping the same operator at the source and the sink. For three-point correlators, using χ opt at the source and the sink requires a different sequential propagator than for χ 1 , but the total number of propagators needed in each case is the same. We chose not to compute all nine combinations of source and sink interpolators in three-point correlators because this would require nine times as many sequential propagators as a standard calculation.

A. Tuning of quark smearing
We use Wuppertal smearing [31],q ∝ ð1 þ αHÞ N q, where H is the nearest-neighbor gauge-covariant hopping matrix constructed using the same smeared links used in the fermion action. The parameter α is fixed to 3.0, and N is varied to produce different smearing widths. The smearing radius is determined by taking a color field φð⃗ xÞ with support only at the origin and then defining a density from the squared norm of the smeared field: ρð⃗ xÞ ¼ jφð⃗ xÞj 2 . Finally, we take the root-mean-squared radius: For the tuning run, we used N ∈ f20; 40; 70; 110; 160g, which correspond roughly to r=a ∈ f3; 4; 5; 6; 7g. Figure 1 shows the effective mass am eff ðtÞ ¼ log CðtÞ CðtþaÞ for each smearing width. For t ¼ 2a and 3a, we can see that the minimum lies near N ¼ 40 and 70. Based on this, we decided to use the same smearing parameter, N ¼ 45, that was used in a previous calculation [11,32,33].

B. Tuning of variational operator
The hybrid operators are constructed using a chromomagnetic field made from smeared gauge links. We take the smeared links to which the fermions couple in the action and then apply additional three-dimensional stout smearing 4 [34]: 20 steps with ρ ¼ 0.1. The traceless part of the clover definition of the field strength tensor is used.
For the second tuning run, we computed the full 3 × 3 matrix of two-point correlators. Solving the GEVP yielded the coefficients v i for χ opt ; we did this at the largest available time separations before the noise became too large. For our choice of operators, the correlator matrix is real and thus v i are also real. Normalizing such that v 1 ¼ 1, we selected χ opt ¼ χ 1 − 4.4χ 2 − 7.3χ 3 . The determination of coefficients from the subsequent full-statistics run is shown in Fig. 2 for a range of t 1 and t 2 in the GEVP. Our selection based on the tuning run is consistent with the values determined at large times in the full run.

A. Two-point correlators
The full 3 × 3 matrix of two-point correlators was computed in the full-statistics run, allowing for a more detailed analysis. We begin by determining the excited energy levels using the variational method. Solving the GEVP at ðt 1 ; t 2 Þ ¼ ð3a; 5aÞ yields an eigenvector v n for each state n. This allows us to define projected correlators C n ðtÞ ¼ v † n CðtÞv n and then compute their effective energies; these are shown in Fig. 3. The two excited energies are nearly degenerate and lie in the range from 1 to 1.5 GeV above the ground-state nucleon; this is similar to the hybrid states observed in Ref. [16]. However, we stress that this should not be considered a reliable determination of the spectrum, as many lower-lying excitations are expected and cannot be identified using our small basis of three operators.
The presence of residual excited-state effects in the effective mass (Fig. 1) suggests that more than three states are needed to describe the two-point correlators. We employ the fit model  4 We have not studied the effect of varying this.
with the energies ordered E 1 < E 2 < Á Á Á. We obtained a good fit to the range t=a ∈ ½3; 12 using N states ¼ 4, which yielded χ 2 =dof ¼ 0.89 (p ¼ 0.68). The four energy levels are shown in Fig. 3. States 3 and 4 are consistent with the effective energies from the GEVP-projected correlators for the two excited states with t=a ¼ 3. The additional energy level, state 2, sits below the two excited states identified by the GEVP. This identification of states between the GEVP and the four-state fit is also supported by Fig 4, which shows the overlap factors normalized to the ground state,  Fig. 1. Note that χ 1 and χ opt were tuned at zero momentum, whereas the full variational analysis is retuned at nonzero momentum. Z i;n =Z i;1 . It also shows the normalized overlap factors for the operator χ opt , which are determined via Z opt;n ¼ v i Z i;n . The fit results indicate that χ 1 has a significant overlap with state 4, which is eliminated in χ opt . The overlap of χ 1 with state 3 is consistent with zero, and this is preserved in χ opt despite the large overlaps of χ 2 and χ 3 with state 3. The operators show no significant difference in the relative overlaps with states 1 and 2; because of this, the GEVP was insensitive to state 2 and was unable to eliminate the coupling of χ opt to it.
Finally, we can compare effective energies produced using different operators. In addition to χ 1 and χ opt , which were tuned at zero momentum using reduced statistics, we also perform a new variational analysis based on the fullstatistics run. In Fig. 1, there is no significant difference between χ opt and the new variational operator. Both of them have smaller effective masses than χ 1 at early times, indicating a significant reduction in excited-state contributions. However, they also suffer from increased statistical uncertainty. Figure 5 shows the same comparison at nonzero momentum. In all cases, χ opt shows smaller excitedstate effects than χ 1 . At the smallest values of t, the new variational operator also shows an improvement over χ opt , and this effect grows with momentum. This is not surprising, as the new variational operator is tuned for each momentum, whereas χ opt was chosen at ⃗p ¼ 0.

B. Forward matrix elements
We have only computed two combinations of source and sink interpolators in three-point correlators: those with the same interpolator at the source and the sink, which is chosen to be either χ 1 or χ opt . For comparing these two setups, we start by considering observables computed from isovector matrix elements at zero momentum: the axial, FIG. 6. Isovector axial, tensor, and scalar charges (g A , g T , and g S ) and isovector momentum fraction hxi u−d determined using the standard interpolator χ 1 (open symbols) and the linear combination of standard and hybrid interpolators χ opt (filled symbols). Ratiomethod data are shown for source-sink separations T=a ¼ 6 (green diamonds), 8 (orange circles), and 10 (blue squares), and are plotted versus the operator insertion time τ, shifted by half the source-sink separation. Summation-method (magenta triangles) data are based on the discrete derivative of sums, ½SðT þ 2aÞ − SðTÞ=ð2aÞ; the two points correspond to T=a ¼ 6 and 8. tensor, and scalar charges (g A , g T , and g S ), and the average momentum fraction hxi u−d .
Results are shown in the "plateau plots" of Fig. 6. The behavior depends strongly on the observable. For g A , χ opt produces a much stronger dependence on the operator insertion time, τ, indicating a significant enhancement of excited-state contributions. The opposite is true for g T , where excited-state effects appear to be significantly suppressed when using χ opt . The difference between the two operators is relatively small for g S and hxi u−d . It is particularly problematic that χ opt , which appears to be an improved operator based on the two-point correlator, produces significantly enhanced excited-state effects in g A . In addition, across all observables χ opt produces consistently larger statistical uncertainties.
We have also explored the use of simultaneous four-state fits to two-point and three-point correlators. However, for each operator O the corresponding fit model has ten independent unknown matrix elements hn 0 jOjni. Since we have computed three-point correlators with only two combinations of source and sink interpolators, the fits are unable to constrain most of the operator matrix elements. This prevents the use of four-state fits to understand what causes the differences between χ 1 and χ opt for estimating matrix elements.

C. Form factors
The isovector form factors F 1 and F 2 of the electromagnetic current are shown in Fig. 7 and the form factors G A and G P of the axial current are shown in Fig. 8. The results are rather mixed. For F 2 near a 2 Q 2 ¼ 0.3 and for F 1 , χ opt produces a weaker dependence on the source-sink separation than χ 1 , indicating the suppression of excitedstate contributions. However, the opposite is true for F 2 at lower Q 2 and for G A . For G P , the two operators produce similar excited-state effects, which are very large at low Q 2 . Chiral perturbation theory predicts this large excited-state contribution [25] as a result of nucleon-pion states, which the hybrid operators are unlikely to help remove.

IV. CONCLUSIONS
The use of a variational basis comprising a standard interpolating operator and hybrid ones presents the possibility of reducing excited-state contamination in nucleon structure calculations at low computational cost. In twopoint correlators, this is borne out, as seen in Figs. 1 and 5. However, in nucleon structure observables that depend on three-point correlators there is no consistent result. The tensor charge shows significantly reduced excited-state effects, whereas the axial charge shows increased effects. 5 Other observables show little change and for form factors the result can depend on Q 2 . If one also takes into account the increased statistical uncertainty, then the setup using hybrid interpolators appears to be not worth pursuing further in its current form.
Results from the four-state fit indicate that the variational procedure succeeds at suppressing some excited states while another lower-lying state is unaffected. The presence of many relevant excited states suggests ways of understanding the results: one possibility is that for the axial charge, the variational procedure spoils a partial cancellation of contributions between two different excited states. Since operator matrix elements do not play a role in the GEVP, there is no guarantee that reducing the net contribution from excited states in the two-point correlator will do the same in three-point correlators.
Reliably improved results using the variational method can only be obtained by being in the asymptotic regime of Ref. [4]. This means that all low-lying excitations must be identified, including multiparticle states, which in practice require nonlocal interpolators. Doing so amounts to a challenging computation; however, as this work has shown, half-measures (such as using a small number of local interpolators) do not consistently reduce excited-state contamination in nucleon structure observables.