Excited-state effects in nucleon structure on the lattice using hybrid interpolators
Jeremy R. Green Michael Engelhardt Nesreen Hasan Stefan Krieg Stefan Meinel John W. Negele Andrew V. Pochinsky Sergey N. Syritsyn
hep-lat




  

jeremy.green@desy.de
  NIC, Deutsches Elektronen-Synchroton, 15738 Zeuthen, Germany


  Department of Physics, New Mexico State University, Las Cruces, NM 88003-8001, USA


  Bergische Universität Wuppertal, 42119 Wuppertal, Germany
  IAS, Jülich Supercomputing Centre, Forschungszentrum Jülich, 52425 Jülich, Germany


  Bergische Universität Wuppertal, 42119 Wuppertal, Germany
  IAS, Jülich Supercomputing Centre, Forschungszentrum Jülich, 52425 Jülich, Germany


  Department of Physics, University of Arizona, Tucson, AZ 85721, USA
  RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA


  Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA


  Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

 
  RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA
   Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA



  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.


Excited-state effects in nucleon structure on the lattice using hybrid interpolators
    Sergey N. Syritsyn
    December 30, 2023
====================================================================================





§ 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 <cit.> that
makes it impractical to evolve long enough in Euclidean time.

The variational method <cit.> 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. <cit.>, which used bases of
interpolators with different smearing widths and different site-local
spin structures, and Ref. <cit.>, 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[The comparison is
  more difficult when using the distillation method, which uses
  timeslice sources rather than the point sources used in standard
  calculations.].

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 <cit.> that contain a gluonic
excitation[The resulting basis of interpolators is similar to
  the one called ℬ_3 in Ref. <cit.>.]. In
Ref. <cit.>, it was found that the latter have the
next-largest overlaps onto the ground state, after the standard
interpolator. 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,
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 <ref> discusses
our lattice setup, the basis of interpolating operators, and two
tuning runs. Results for two-point correlators are discussed in
Section <ref>, forward matrix elements are shown in
Section <ref> and form factors are presented in
Section <ref>. Our conclusions are given in
Section <ref>.



§ LATTICE SETUP



As this is an exploratory calculation, we use a single lattice
ensemble with a coarse lattice spacing at a heavier-than-physical pion
mass and a relatively small box size; its parameters are summarized in
Table <ref>. This has 2+1 flavors of tree-level
improved Wilson-clover quarks coupled to the gauge links via two
levels of HEX smearing <cit.>.



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. <cit.>; 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^-Δ E T/2, where
Δ E is the energy gap to the lowest excited state — and the
summation method, for which they decay as Te^-Δ E T.

Given a set of N interpolating operators {χ_i}, 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,

    C_ij(t) = ⟨χ_i^(t) χ_j^†(0) ⟩,

and then solve a generalized eigenvalue problem (GEVP)

    C(t_2)v = λ C(t_1)v

for some choice of (t_1,t_2). It has been
shown <cit.> 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 <cit.>. 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. <cit.>. 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 <cit.>. This defeats the diagonalization
procedure of Eq. (<ref>) 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

    χ_1 = ϵ_abc(ũ_a^T Cγ_5P_+ d̃_b) ũ_c,

where P_+=(1+γ_4)/2 is a positive parity (nonrelativistic)
projector and q̃ 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 <cit.>. Of the three
possible site-local nucleon operators, χ_1 is the only one that
can be constructed using positive-parity-projected quark fields (see
e.g. Appendices B and C of Ref. <cit.>).

We also consider hybrid operators, introduced in
Ref. <cit.>, which include an insertion of the
chromomagnetic field B_i=-1/2ϵ_ijkF_jk and are
interpreted as having a nontrivial gluonic excitation. Using the
clover discretization of F_μν <cit.>, 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:

    χ_2    = ϵ_abc[ (B_i ũ)_a^T Cγ_j P_+ d̃_b] γ_iγ_j ũ_c
     - ϵ_abc[ ũ_a^T Cγ_j P_+ (B_i d̃)_b] γ_iγ_j ũ_c,
    χ_3    = ϵ_abc[ (B_i ũ)_a^T Cγ_j P_+ d̃_b] P_ijũ_c
     - ϵ_abc[ ũ_a^T Cγ_j P_+ (B_i d̃)_b] P_ijũ_c,

where P_ij=δ_ij-1/3γ_iγ_j. 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.



 §.§ Tuning of quark smearing


We use Wuppertal smearing <cit.>,
q̃∝ (1+α H)^Nq, 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⃗)=|φ̃(x⃗)|^2. Finally, we take the
root-mean-squared radius:

    r^2 = ∑_x⃗ |x⃗|^2ρ(x⃗)/∑_x⃗ρ(x⃗).

For the tuning run, we used N∈{20,40,70,110,160}, which
correspond roughly to r/a∈{3,4,5,6,7}.



Figure <ref> shows the effective mass
am_eff(t)=logC(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 <cit.>.



 §.§ 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[We have not studied the
  effect of varying this.] <cit.>: 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. <ref> 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.



§ RESULTS




 §.§ 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. <ref>. 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. <cit.>. However, we stress
that this should not be considered a reliable determination of the
spectrum, as many lower-lying excitations are expected and can not be
identified using our small basis of three operators.



The presence of residual excited-state effects in the effective mass
(Fig. <ref>) suggests that more than three states are needed
to describe the two-point correlators. We employ the fit model

    C_ij(t) = ∑_n=1^N_states Z_i,n Z_j,n e^-E_n t,

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. <ref>. 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 <ref>, which shows the overlap
factors normalized to the ground state, 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 full-statistics run. In
Fig. <ref>, 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 <ref> shows the same comparison at
nonzero momentum. In all cases, χ_opt shows smaller
excited-state 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.



 §.§ 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, tensor, and scalar charges
(g_A, g_T, and g_S), and the average momentum fraction
⟨ x⟩_u-d.



Results are shown in the “plateau plots” of
Fig. <ref>. 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
⟨ x⟩_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
𝒪 the corresponding fit model has ten independent unknown
matrix elements ⟨ n'|𝒪|n⟩. 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.



 §.§ Form factors







The isovector form factors F_1 and F_2 of the electromagnetic
current are shown in Fig. <ref> and the form factors G_A
and G_P of the axial current are shown in Fig. <ref>. The
results are rather mixed. For F_2 near a^2Q^2=0.3 and for F_1,
χ_opt produces a weaker dependence on the source-sink
separation than χ_1, indicating the suppression of excited-state
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 <cit.> as a result of nucleon-pion states,
which the hybrid operators are unlikely to help remove.



§ 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 two-point correlators, this is borne out, as
seen in Figs. <ref> and <ref>. 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[The opposite effect was seen using a different basis
  of interpolators in Ref. <cit.>, where the variational
  analyses produced reduced excited-state effects for g_A and
  enhanced excited-state effects for g_T, compared with the
  interpolator S_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. <cit.>. 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.



We thank the Budapest-Marseille-Wuppertal collaboration for making
their configurations available to us. Calculations for this project
were done using the Qlua software suite <cit.>, and made use of
the QOPQDP adaptive multigrid solver <cit.>.

This research used resources on the supercomputer JUWELS <cit.>
at Jülich Supercomputing Centre (JSC); we acknowledge computing time
granted by the John von Neumann Institute for Computing (NIC).

SM is supported by the U.S. Department of Energy (DOE), Office of
Science, Office of High Energy Physics under Award Number
DE-SC0009913. SM and SS are also supported by the RIKEN BNL Research
Center under its joint tenure track fellowships with the University of
Arizona and Stony Brook University, respectively. ME, JN, and AP are
supported in part by the Office of Nuclear Physics of the U.S. DOE
under grants DE-FG02-96ER40965, DE-SC-0011090, and DE-SC0018121,
respectively. SK and NH received support from Deutsche
Forschungsgemeinschaft grant SFB-TRR 55.



