Excited states in nucleon structure calculations

Excited state contributions represent a formidable challenge for hadron structure calculations in lattice QCD. For physical systems that exhibit an exponential signal-to-noise problem they often hinder the extraction of ground state matrix elements, introducing a major source of systematic error in lattice calculations of such quantities. The development of methods to treat the contribution of excited states and the current status of related lattice studies are reviewed with focus on nucleon structure calculations that are notoriously affected by excited state contamination.


Introduction
The lattice formulation of quantum chromodynamics (QCD) is the only known ab initio method to study the properties of hadrons from first principles. Through the development of more powerful computers and algorithms over the last decades it has become possible to perform precise lattice calculations for a broad range of applications. However, in order to make contact with experimental results from a controlled extrapolations of lattice results to the physical point, it is crucial to not only improve statistical precision but also to achieve control over systematics. A lot of progress has been made regarding chiral, continuum and finite size extrapolations which are now part of many current lattice studies, as simulations at physical quark mass and including multiple values of the lattice spacing and volumes have become feasible in the last few years. Still, in the presence of an exponential signal-to-noise problem, which is typical for observables involving baryons, a fourth systematic effect often remains a challenge, i.e. residual excited state contamination. The reason for this is that the signal is lost in noise at Euclidean time separations that are still insufficient to gain ground state a e-mail: kottnad@kph.uni-mainz.de (corresponding author) dominance. Therefore, isolation of the ground state cannot be achieved in a direct way without taking further measures.
Most notably affected by this kind of systematic effect are lattice QCD calculations of nucleon structure, that cover a rich variety of observables. Among the most basic such observables are nucleon charges that can be computed from forward matrix elements at zero momentum transfer. An important example is the nucleon axial charge g u−d A = 1.2724 (23) [1] which is experimentally measured in neutron β-decay and often serves as a benchmark observable for nucleon structure calculations in lattice QCD. On the other hand, there is less experimental information available on scalar and tensor charges which may give beyond the Standard Model (BSM) contribution to the nucleon β-decay [2]. Therefore, control over excited states is all the more important in lattice QCD calculations of these observables which may provide crucial input on searches for BSM physics as they are relevant to dark matter searches [3]. Furthermore, the tensor charge plays a role in BSM searches for C Pviolation as it controls the contribution of quark electric dipole moments to the neutron electric dipole moment [4]. Another example for an observable of great phenomenological interest at zero momentum transfer is the average quark momentum fraction which contributes to the nucleon spin decomposition [5]. At non-vanishing momentum transfer various form factors are studied, e.g. electromagnetic and axial form factors and the corresponding radii that parametrize the slope of these form factors at vanishing momentum. In particular, the proton radius has received a lot of attention in the last decade due to the so-called proton radius puzzle, see e.g. refs. [6][7][8][9][10][11][12], although this has likely been resolved by now [13] and the most recent measurement from electron-proton scattering [14] indeed agrees with the very precise results obtained from measurements of the Lamb shift measured for muonic hydrogen [8,12]. Axial form factors in turn are experimentally less well-known [15][16][17] but may provide critical input for future experiments related to neutrino physics [18,19], which makes them attractive observables to be studied in lattice QCD.
The aim of this article is to review the state-of-the-art of methods that are applied to tame residual excited state contamination in modern nucleon structure calculations and to perform a critical assessment of the efficacy of these methods. Therefore, further systematics related to chiral, continuum and finite size extrapolations that affect these calculations are not within the scope of this review. From a physics point of view, a subset of nucleon structure observables will be considered for which excited state contamination is known to be an important, or even the dominant systematic. The selection comprises observables that can be computed with sufficiently good statistical precision and for which dedicated studies of excited state systematics and related methods can be found in the literature.
The review is organized as follows: In the next section basic methods for nucleon structure calculations in lattice QCD are summarized and observables relevant to this review are defined. In Sect. 3 theoretical aspects of excited states and their expected effects on such calculations are discussed. The following sections are dedicated to the review of various methods that have been employed in the treatment of excited state contamination and their respective applications in recent lattice studies. More specifically, methods that aim at additional suppression of excited states by summation over the operator insertion are reviewed in Sect. 4, while Sects. 5 and 6 deal with the various application of multi-state fits and the variational approach, respectively. In the end, a summary on the advantages and shortcomings of the individual methods and a brief outlook on possible future developments is given.

Nucleon structure calculations in lattice QCD
The calculation of nucleon properties in lattice QCD is based on the numerical evaluation of n-point functions in discretized Euclidean spacetime. For the simple example of a two-point function with source at the origin and sink at (x, t), the corresponding spectral decomposition contains all possible hadronic states |k, p with integer label k and energies E k ( p) that share the same continuum quantum numbers compatible with the choice of the interpolating operator χ(x, t). This typically includes excitations of the ground state as well as matching multi-particle states. The reason for this is that the wave functions of the hadronic states |k are unknown and thus it is impossible to directly construct interpolating operators that couple only to a desired state k while the remaining overlap factors 0| χ |k for k = k vanish exactly. Moreover, rotational symmetry is broken at finite values of the lattice spacing which may lead to additional mixing between operators that would otherwise fall into different irreducible representation in the continuum limit. Therefore, the extraction of the ground state energies and overlap factors requires to compute C 2pt ( p, t) at large Euclidean time separation t, such that all higher terms are exponentially suppressed This motivates the definition of the effective energy which is called effective mass for p = 0. For large values of t it asymptotically approaches a plateau from which the energy can be extracted. The parameter τ is commonly set to τ/a = 1. In practice, the plateau can be identified once the residual slope has become negligible compared to the statistical error of the individual data points.

Nucleon matrix elements and ratio method
While hadron masses can be readily obtained from two-point functions, the study of the structure of hadrons from lattice QCD relies on the computation of matrix elements. With respect to excited states a particularly relevant application are nucleon matrix elements (NMEs) where N ( p i , s i ), N ( p f , s f ) denote nucleon states with initial (final) state momentum p i ( p f ) and spin s i (s f ). On the right hand side of u( p i , s i ),ū( p f , s f ) are the corresponding Dirac spinors, while W X μ 1 ...μ n (Q 2 ) contains form factors and kinematic factors. The form factor decomposition depends on the choice of the operator insertion that we restrict to a bilinear operator of the form O X μ 1 ...μ n =q(x)Γ X μ 1 ...μ n q(x).
where Γ X μ 1 ...μ n collects Dirac matrices and possibly derivatives, and q(x),q(x) denote quark fields. Actual examples of operator insertions relevant to nucleon structure calculations will be discussed in Sect. 2.3. In practice, most lattice calculations are performed assuming exact isospin symmetry in the light quark sector. In this case, the insertion of the isovector combination leads to the cancellation of quark-disconnected diagrams that are inherently more noisy then quark-connected contributions. For the isoscalar counterpart O X,u+d μ 1 ...μ n this cancellation does not occur, which further increases the signal-to-noise problem. Moreover, for studies of observables that probe the contributions of strange and charm quarks to nucleon structure such as e.g. strange and charm electromagnetic form factors, quark-disconnected diagrams give the only contribution apart from possible mixing with quark-connected contributions under renormalization.
The lattice determination of NMEs as defined in Eq. (5) generally requires the computation of spin-projected twopoint functions and three-point functions where are suitable spin projectors and the operator is inserted at . A common choice for the nucleon interpolating field is given by e.g.
where C is the charge conjugation matrix and u, d denote up-and down quark fields that are usually smeared. Assuming that the source time is zero, the two-point function in momentum space for t > 0 reads where the overlap of the interpolating operator χ with the k-th state A k ( p) = Ω| χ |k, p of momentum p has been introduced in the spectral decomposition in the second line. Similarly, the three-point function can be expressed as where the source-sink separation t sep = t f − t i and the shorthand t ins = t O − t i for the insertion time have been defined and again t i = 0 has been assumed. The extraction of the ground state matrix element 0, p f O X μ 1 ...μ n 0, p i requires cancellation of the unknown overlap factors A k ( p f ) and A l ( p i ). This can be achieved by taking an appropriate ratio of two-and three-point functions. A common choice is the ratio [20] which has been shown to be particularly beneficial with respect to statistical errors in Ref. [21]. For asymptotically large Euclidean time separations ground state dominance is achieved giving access to the matrix element from a plateau of the effective form factor, similar to the extraction of the ground state energy from Eq. (4). For this reason the ratio method for matrix elements is also sometimes referred to as plateau method.

Lattice techniques
The computation of NMEs from the ratio in Eq. (13) requires three-point functions at several source-sink separations t sep to allow for a controlled study of excited state effects. The de facto standard approach is to use sequential inversion through a fixed sink, which gives access to all values of the insertion time t ins for any given value of t sep . However, in this setup each source-sink separation needs additional sequential inversions. Besides, the number of inversion per sourcesink separation increases if non-local operator insertions are used, e.g. for point-split currents or the derivative operators in Eqs. (28)- (30). Another possibility is to perform sequential inversions through a fixed operator insertion at t ins = const [22]. This approach is complementary to the previous one as it gives access to all values of t sep without the need for additional inversions while allowing for also freely choosing momentum projection and operator at sink without incurring extra computational cost. The latter makes this an attractive choice for certain setups used for employing a dedicated variational analysis that will be discuss in Sect. 6. However, additional inversions are needed for each value of t ins as well as every operator insertion which usually outweighs the aforementioned advantages in terms of computational cost.
A third method that effectively eliminates t ins is obtained by summation over the operator insertion. This approach can be considered a particular case of the summation method [23] and will be discussed in Sect. 4.2. In this setup it is possible to change the operator at sink without the need for additional inversion, but again not the operator insertion. It has only been employed in a few studies in recent years [24][25][26].
Most calculations of connected contributions use pointlike sources in combination with smearing of the quark fields q(x) to increase the ground state overlap, hence improving suppression of excited state contamination. This setup gives also full flexibility with respect to momentum projection at the source unlike e.g. timeslice sources that would require new inversion for each momentum. A common choice for the smearing of quark fields is Gaussian smearing [27,28] q(x) →q(x) = (1 + κ G Δ) N smear (15) together with spatial APE smearing [29] of the gauge fields that enter through the three-dimensional Laplacian Δ. In general, the smearing parameter κ G and the number of smearing steps N smear need to be tuned to give optimal results. Besides, smearing can be used as a comparably cheap way to define additional interpolating operators which will be discussed in the context of variational techniques in Sect. 6. For the fixed sink and fixed insertion method the same point-to-all forward propagators can be used for two-and three-point functions.
Since control over excited states depends on large sourcesink separations, methods to reduce the computational cost are of particular importance for nucleon structure calculations as this requires high statistics. Commonly used approaches for the sequential method are the truncated solver method [30] or all-mode averaging (AMA) [31], that allow to increase statistics at reduced cost compared to performing only exact inversions. Further methods that have been investigated include replacing the conventional sequential propagators by a stochastic estimator [32][33][34] which facilitates the reuse of the propagator at the expense of introducing stochastic noise, the coherent sink method [35,36] that allows to invert simultaneously on multiple, temporally separated sequential sources, and finally distillation [37] which is suitable to systematically build a basis of interpolating operators and that has been employed in a recent study in Ref. [38].
Furthermore, a lot of progress has been made in the last few years regarding the inclusion of disconnected diagrams which has become feasible through the development and application of variance reduction techniques such as hierarchical probing [39], low-mode deflation [40], the one-end trick [41], the hopping parameter expansion [42,43] and combinations thereof [44,45]. This lead to a substantial reduction in computational cost compared to e.g. using naive stochastic all-to-all estimators. However, the relative statistical error of quark-disconnected contributions remains inherently larger than the one of purely quark-connected diagrams even if gauge noise is reached.

Observables
There exists a large variety of nucleon structure related observables that are computed in lattice QCD. However, for the purpose of this review mainly a subset of observables will be considered for which uncontrolled excited state contributions have been shown to produce a bias at the current level of precision. As a first example, the insertion of a local vector current for a generic quark flavor or the corresponding point-split current gives rise to the electromagnetic form factors Usually, mean square radii r 2 are defined from the expansion of a form factor in the Euclidean four-momentum transfer Q 2 around zero For example the electric and magnetic mean square charge radii of the nucleon are given by where F 1 (Q 2 ), F 2 (Q 2 ) enter through the Sachs form factors  [46]) with M π ≈ 350 MeV, a = 0.06426 fm and T × L 3 = 128a × (48a) 3 as part of a study published in Ref. [47] for a fixed number of measurements on each configuration independent of the value of t sep Note that for the neutron G n E (0) = F n 1 (Q 2 ) = 0 and the factor 1/F n 1 (0) is dropped from Eq. (19) to obtain a finite definition. Since radii are not explicitly momentum dependent they are more readily compared to other lattice calculations than the full momentum dependence of the form factors themselves. Moreover, radii can be directly compared to experimental results as well.
Another set of observables for which excited state effects are particularly relevant are nucleon charges defined at zero momentum transfer. However, for the vector current the electric charge G E (0) is conserved and the computation of the magnetic moment μ M = G M (0) is hindered by a momentum prefactor in the form factor decomposition which vanishes at Q 2 = 0. The latter introduces additional systematic effects as the extraction of μ M requires either an extrapolation in Q 2 or position space methods, see e.g. Refs. [48,49]. The situation is different for NMEs from other local bilinear operators, i.e. the axial, scalar and tensor currents for which the associated charges g A,S,T are not conserved. While radii require knowledge of the momentum dependence of the form factor, the axial, scalar and tensor charges are readily computed at zero-momentum transfer (up to renormalization) from a simple ratio of a three-and a single twopoint function for equal initial and final state momentum p ≡ p i = p f . For example, the isovector nucleon axial charge g u−d A that is defined as for p f = p i ≡ p where P( p, s f ) and N ( p, s i ) refer to the final proton and initial neutron state, is obtained from Some example data for the ratio as a function of t sep for the three charges are shown in the first three panels of Fig. 1.
Since radii are less straightforward to compute, it is the experimentally well-known isovector axial charge g u−d A that is considered a benchmark quantity for lattice QCD, although the statistical quality for the signal of a vector current insertion is superior to any of the operators in Eqs. (22)- (24). Still, the axial charge has rather good signal quality compared to e.g. the scalar charge shown in the upper right panel of Fig. 1, which makes statistical errors of order O(1%) achievable in modern lattice simulations. Beyond charges also the axial and (induced) pseudoscalar form factors G A (Q 2 ) and G P (Q 2 ) have very recently received attention with focus on excited states from both, the theory [50][51][52] and the lattice side [36,53,54]. This will be discussed in more detail in Sects. 3.2 and 5.2.
In addition to nucleon charges, the related form factors at non-zero Q 2 and radii, NMEs have been studied for operators involving derivatives, e.g. the twist-2 operators where {...} and [...] denote symmetrization with subtraction of the trace and antisymmetrization, respectively, and the symmetric derivative is defined as . Matrix elements of these operators occur for example in the computation of the average quark momentum fraction of the nucleon x q , helicity and transversity moments x Δq and x δq or for the generalized parton distribution functions contributing in the computation of the spin decomposition of the proton, see e.g. Refs. [55,56]. The latter involves also a gluonic operator insertion. In principle, any such NME may exhibit sizable excited state effects, as will be discussed in the next section. A rather well-studied example is the average quark momentum fraction of the nucleon for which example data is shown in the lower right panel of Fig. 1 for the isovector combination. However, for many of these observables a systematic treatment of excited states is more difficult due to the signal quality which especially becomes an issue if quarkdisconnected contributions are involved, see e.g. Refs. [57][58][59]. It is for this reason that methods to treat excited states are currently of particular importance compared for simple observables like g u−d A,S,T , x u−d or in precision calculations of the nucleon radius because the achievable statistical error for these quantities does not allow to conceal the systematic effect of residual excited state contamination in final results.

Excited states
While excited states are present in essentially any lattice calculation, they become especially an issue in calculations of nucleon structure. In the following the reasons for this are discussed in more detail and some theoretical expectations from chiral perturbation theory (χ PT) and modeling excited states are reviewed.

The signal-to-noise problem
The primary cause why excited states are a persistent issue in nucleon structure calculation is an exponential signal-tonoise problem that prevents computation of n-point functions at large Euclidean time separations as required for e.g. obtaining a plateau in the effective mass in Eq. (4) or for the ratio method in Eq. (14). For the nucleon this problem has been formulated a long time ago [60,61] from a field theoretical computation of the variance σ 2 stat (t) of the nucleon two-point function, showing that the dominating contribution at large values of t is a three pion state, i.e.
For the signal-to-noise behavior of a nucleon two-point function this implies unlike the pion, for which a constant signal-to-noise ratio is expected. An example is shown in Fig. 2 for the effective mass of the nucleon and the pion computed on an ensemble at physical quark mass. Clearly, the signal for the nucleon is lost between 1.5 fm and 2 fm, while for the pion at rest the ratio M π (t)/ΔM π (t) remains constant as expected.
Regarding the behavior of excited states and the signal-tonoise problem in the ratio method, it is instructive to consider the two-state truncation for the zero-momentum case (initial and final state at rest) of the ratio in Eq. (25) where momentum arguments as well as Dirac indices have been suppressed and Δ = E 1 − E 0 denotes the energy gap between the first excited state and the ground state. Overlap factors and matrix elements from the spectral decomposition of the two-and three-point functions in Eqs. (11,12) have been collected in the A kl factors. In particular, the first term A 00 is proportional to the ground state matrix element, e.g. the nucleon charges g S,A,T if the generic operator O X μ 1 ,...,μ n in Eq. (25) is chosen appropriately and A 01 = A 10 holds for the zero-momentum case. Using the midpoint of the ratio as an estimate at any given source-sink separation, the leading excited state contamination σ esc is therefore expected to scale as Note that at non-zero momentum transfer the behavior of the ratio midpoint is similar, however, the time dependence in Eq. (33) is no longer symmetric if initial and final state differ by momentum. Considering the forward scattering case and neglecting the N π interaction one may assume at least for a qualitative analysis that the leading gap is close to 2M π , i.e. the scaling is approximately ∼ exp(−M π t sep ). In Fig. 3 this behavior is plotted for several pion masses as a function of t sep and it is obvious that the issue of residual excited states for the ratio method becomes more severe towards the physical quark masses. For example, at a source-sink separation of 1.5 fm that is typically reached in NME lattice calculations, σ esc is roughly one order of magnitude larger at physical light quark mass compared to the situation at M π = 400 MeV. Assuming that the matrix element is of O(1) this implies corrections of up to ∼ 40% for t sep = 1.5 fm at physical pion mass. Therefore, it is doubtful if the values of t sep that can be reached in lattice calculations are sufficiently large to ensure control over excited states in the ratio method particularly at small light quark mass. Furthermore, excited-state effects can be strongly operator-dependent which is also observed empirically in lattice calculations; see e.g. Fig. 1. Consequently, it cannot be concluded from observing ground state dominance for an observable in the ratio method that similar values of t sep will be sufficient for a different observable.  On the other hand, the statistical error of the effective form factor behaves as i.e. it exponentially grows with t sep . Even at moderately heavy quark mass the effect is sizable as shown in Fig. 4 for a selection of zero-momentum NMEs measured with fixed statistics at several values of t sep . In fact, going from t sep = 1 fm to t sep = 1.5 fm the signal-to-noise ratio becomes roughly four times smaller. This means that any desired reduction in residual excited state contamination by increasing t sep competes with the statistical precision of the final result. Moreover, lowering the light quark mass increases both σ stat and σ esc individually. If one demands that the statistical error σ stat and the systematic uncertainty σ esc scale in the same way, i.e. σ esc = const · σ stat ≡ σ , the required statistics N for a target uncertainty σ can be inferred in the asymptotic regime from Eqs. (34,35) leading to [62] In the review in Ref. [62] it has been further pointed out that at physical quark mass and again assuming Δ = 2M π the exponent is approximately −13 which must be compared to the naive factor of −2 obtained when neglecting the effects of excited states. It is in this sense that the presence of excited states strongly enhances the issue of the signal-tonoise problem in lattice calculations of NMEs, particularly when approaching physical quark mass.

Multi-particle states and theory predictions
A further complication arises from the particular excited state pattern expected in nucleon structure calculations, i.e. a dense spectrum of multi-particle states in addition to resonances, which limits the efficacy of a variational treatment or fits. In finite volume only discrete momenta p = 2π L n with n = (n 1 , n 2 , n 3 ) T and integer n i are allowed and any combination of a single nucleon with an arbitrary number of pions that is compatible with symmetries and momentum conservation will occur in the spectrum. Increasing the (spatial) box size L while keeping the remaining physics parameter fixed, will reduce the size of a momentum unit 2π/L, thus leading to a denser spectrum. Neglecting the interactions between the nucleon and the pions, an approximation of the multiparticle spectrum can be obtained that is shown in the left panel of Fig. 5 for physical pion mass. Up until M π L = 4 the lowest lying excited state is a N ππ state, which justifies the previously made assumption of Δ = 2M π for the lowest gap. Note that at heavier pion mass this remains true for even larger values of M π L. Moreover, corrections to the lower, non-interacting N π levels are rather small as can be seen in the right panel of Fig. 5. The results shown in this figure have been obtained in Ref. [63] using the χ PT in infinite volume together with the Lellouch-Lüscher formalism. For higher N π states E 1400 MeV deviations become more apparent, although the general pattern is not affected.
At heavier-than-physical quark mass the spectrum of multi-particle states is thinned out while keeping M π L fixed. Therefore, systematics related to excited states are again expected to become more severe at lighter quark masses. This is demonstrated in Fig. 6, where the ratio M N (t)/M N of the effective mass and the fitted ground state value is shown for four ensembles with quark masses covering a range from ∼ 135 MeV up to 350 MeV. At small values of t a clear hierarchy is observed before the signal becomes to noisy, i.e. a plateau is approached more rapidly at heavier quark masses.
In the last few years there have been several studies [50][51][52][64][65][66][67][68][69] carried out using χ PT to predict effects of excited states in nucleon structure calculations. While recent lattice calculations are performed near or at physical quark mass which already removes a major systematic effect, other restrictions with respect to the applicability of χ PT results to lattice data remain. This concerns e.g. the size of the available source-sink separations in case of three-point functions.
Regarding the use of smeared interpolating operators it has been pointed out in Ref. [66] that they are mapped onto the point-like nucleon field in the effective theory, provided that the smearing radius is small compared to the Compton wavelength of the pion. The resulting effective operators containing the pion-nucleon coupling in the second term then only differ by the value of a low energy constant (LEC) α for different smearings. Furthermore, at leading order this LEC is canceled in ratios, hence χ PT predictions for excited state corrections are independent of the actual choice of smearing at leading order. At the very least, these studies provide qualitative insight into the behavior of excited state contamination, but in more recent work χ PT predictions have also been used to systematically remove excited state contamination from lattice data; see e.g. Ref. [50]. Excited states in nucleon two-point functions have been studied in Refs. [64][65][66]69]. For the nucleon effective mass it was found that the excited state correction due to N π contributions are expected to be below 2% at t ≥ 0.5 fm and to become a sub-percent effect for t ≥ 1 fm. This is roughly consistent with empirical findings in lattice studies e.g. considering the behavior of the (total) relative excited state contribution in Fig. 6 as a function of t. In the most recent Ref. [69] the study has been extended to three-particle (N ππ) states which where found to contribute at most at the permille level and thus considered to be negligible for all practical purposes in the foreseeable future.
For three-point functions and the resulting matrix elements the situation is more complicated. In the past, the main focus has been on N π contributions in the three-point function with an axial vector insertion relevant for g u−d A [64,65,67,68]. The predicted effect on g u−d A is an overestimation of at least several percent at typical values of t sep 1.5 fm that are accessible in lattice simulations. An effect of similar size has been predicted for g u−d T , while for g u−d S the overestimation was determined to be ∼ 50% larger [67]. In Ref.
[68] the first moments of parton distribution functions, x u−d , x Δu−Δd and x δu−δd have been analyzed in addition to charges and it has been found that the contribution of the N π state lead to an overestimation of 5-10% at t sep = 2 fm in all of these observables. Furthermore, it has been pointed out that the source-sink separations that are currently accessible in lattice simulations are too small for a direct application of χ PT and hence these predictions are of limited applicability. In particular, resonances that are not included in χ PT might give non-negligible contributions. This may explain why the predicted effect for g u−d A is opposite to what is actually observed in lattice simulations, i.e. excited states lead to smaller values for g u−d A . An example is shown in the first panel of Fig. 1 where the effective form factor for g u−d A appears to approach the asymptotic value from below.
In Ref. [63] the issue of excited states in lattice calculations of g A has been revisited using the Lellouch-Lüscher formalism and experimental results for the N π scattering phases. While the inclusion of information on the spectrum has been found to only give a small correction to the χ PT prediction for the excited state contamination in g A , it has been argued that plausible model assumptions for the overlap factor and particularly postulating a sign change in the infinite volume axial-vector transition amplitude can qualitatively reproduce the observed behavior of lattice data which is shown in Fig. 7. The mechanism behind this indeed turned out to be a cancellation between the (positive) contribution of lower lying states and the (negative) contribution of higher excited states to the overall excited state contamination.
Recently, χ PT studies of N π state contamination have been extended beyond zero-momentum for the axial form factor G A (Q 2 ) and the related, induced pseudoscalar form factorG P (Q 2 ) in Ref. [50]. For G A (Q 2 ) it has been found that the leading order χ PT prediction for the N π contribution remains at the 5% level at t sep = 2 fm and is basically independent of Q 2 . This is in stark contrast to the findings forG P (Q 2 ) for which an underestimation of 10 to 40% due to N π states has been predicted that is strongly dependent on Q 2 . This is in good agreement with numerical findings of lattice calculations, unlike the prediction of an overestimation Fig. 7 Prediction of the excited state contamination in the ratio method for g A from Ref. [63]. The two upper curves represent the leading order χPT prediction, but taking into account interacting values of the Lellouch-Lüscher factors. The remaining two pairs of curves correspond to certain variations of model assumptions for the overlap factor and axial-vector transition amplitude. The lower pair of curves qualitatively reproduces the behavior observed in lattice calculations and the two values of M π L roughly cover the range found in modern lattice calculations. The figure has been originally published in Ref. [63] and is reproduced under the Creative Commons CC-BY license in case of the axial form factor. This very different behavior of the two axial form factors has been attributed to how N π states contribute in the two observables: For G A (Q 2 ) the entire tower of N π state contributes, while forG P (Q 2 ) the N π excited state contamination consists of only a single N π state associated with a pion with its spatial momentum determined by the momentum transfer Q 2 . Consequently, the prediction forG P (Q 2 ) remains valid for much smaller source-sink separations than for G A (Q 2 ) and it has been demonstrated that this can used to reliably correct lattice data at source-sink separations well below 1.5 fm. An example for this correction of lattice data from Ref. [70] is shown in Fig. 8 together with experimental data and predictions from the pion pole dominance model.
In the same study in Ref. [50] N π states have also been found to be a likely cause for the strong excited state contamination observed in the effective form factor obtained from an operator insertion of the time component of the axial vector current, see e.g. Refs. [53,71]. While the excellent agreement of lattice data with the χ PT prediction demonstrated in Fig. 9 may be accidental to some extent, it can be concluded that the N π contribution to the excited state contamination is large and dominant in this case and source-sink separations even beyond 3 fm might be required to sufficiently suppress them. Note that this is also the reason why lattice computations of the axial charge and G A (Q 2 ) commonly employ three-point functions from spatial components of the axial vector insertion, which exhibit much milder excited state contamination.  [70] with and without correction for N π states. The lattice data is generally well described by the pion pole dominance (ppd) model before and after correction for the N π states and is found to be in agreement with experimental results after the correction. The figure has been originally published in Ref. [50] and is reproduced under the Creative Commons Attribution 4.0 International license  Furthermore, in Ref. [51] an analysis of the N π -state contamination has been performed for the pseudoscalar form factor G P (Q 2 ). Similar to the induced form factorG P (Q 2 ) a large contamination of −20% up to 50% depending on Q 2 has been predicted at t sep = 2 fm. Moreover, it has been argued that a N π -state contamination is the most relevant source for the large violations of the generalized Goldberger-Treiman relation that have been observed in lattice simulations [70,[72][73][74][75]. This relation establishes a connection between the form factors G A (Q 2 ) andG P (Q 2 ) associated with the partially conserved axial-vector current (PCAC) and the form factor G P (Q 2 ) of the pseudoscalar density. A quantitative measure of its violation is given by the deviation from unity for the following expression [53,72] wherem denotes the average bare PCAC quark mass. A second, closely related ratio that can be directly tested from lattice data is the pion-pole dominance (PPD) hypothesis relating G A (Q 2 ) andG P (Q 2 ). An example for the violation of these relations by lattice results are the lower sets of data points in Fig. 10 that have been obtained using conventional fit ansätze to two-and three-point functions in two recent studies in Refs. [36,54]. Again, the theoretical findings in Ref. [51] were corroborated by lattice data, even at source-sink separations as small as t sep = 1.3 fm, although no definite conclusion has been drawn in this study. This has motivated further lattice investigations of the issue in Refs. [36,53,54]. In Ref. [53] a projection method has been introduced to remove excited state contamination in nucleon form factor calculations, however, this only lead to improvement for r PCAC but not for r PPD . Subsequently, it has been argued in Ref. [52] using chiral perturbation theory that this projection method in fact enhances the N π contamination in G P (Q 2 ) and the improvement observed for the PCAC relation in Eq. (38) is caused by the enhanced N π contribution in G P (Q 2 ) compensating the underestimation ofG P (Q 2 ). However, no such cancellation takes place for r PPD which does not depend on G P (Q 2 ). Finally, in Ref. [54] it has been demonstrated that the deviations from unity for Eqs. (38) and (39) observed in previous lattice studies can indeed be attributed to a low-lying excited state which is missed in commonly used fits, e.g. fit strategy S 2pt in the left panel of Fig. 10. This result and the correspondingly modified fit procedure S A 4 determining the energy of this state from a fit to the three-point function will be discussed in some more detail in Sect. 5.2 in the context of multi-state fits. A different ansatz explicitly modeling the excited state contamination in presence of the pion pole has been introduced in Ref. [36] and has been found to lead to very similar results. Results from the modified analysis strategies in both studies are visible in Fig. 10, i.e. the upper sets of data points in both panels, which are compatible with one.

Summed operator insertions
The ratio method introduced in Sect. 2.1 depends on two Euclidean time separations to become large such that excited states are sufficiently suppressed. As discussed in the previous section this can hardly be achieved in the presence of the nucleon signal to noise problem. While it may be possible to explicitly remove specific excited state contamination for certain matrix elements like e.g. N π states in axial and pseudoscalar form factors, there is no general way how this can be achieved using additional, theoretical input only. Therefore, methods are needed that improve the suppression of excited states using data from the available source-sink separations and that are equally applicable to a broad class of nucleon structure observables. A simple approach that satisfies this requirement is the summation method that has been originally published in Ref. [23]. It is based on summing the operator insertion over t ins , leaving t sep as the only time dependence.

Summation method
In its commonly used version [76,77] the summation is performed at the level of the ratio in Eq. (13) running only over timeslices between source and sink, i.e. t ins ∈ t ex , t sep − t ex with some additional freedom of leaving out further timeslices parameterized by t ex The leading excited state contamination related to the gap Δ is still present but more strongly suppressed compared to the excited state contamination at the midpoint of the effective form factor in the ratio method that scales with O(e −Δt sep /2 ). It should be noted, that only at zero-momentum transfer the gap Δ is actually the same as in Eq. (33). However, the resulting suppression of the leading excited state contamination is always of O(e −Δt sep ), where Δ is the smallest gap to the ground state, see also Refs. [78,79]. Figure 11 shows an example for the improved suppression of excited states compared to the ratio method for x u−d . While the midpoint of the effective form factor approaches the result from the summation method, the errors overlap only at the largest source-sink separation of t sep = 1.54 fm and the central value from the ratio method is still larger. Key advantages of the summation method are that it is trivial to implement, does not introduce model dependence or additional parameters (apart from the choice of t ex ) and that it is often possible to reuse data generated for other methods e.g. the ratio method or multi-state fits. Therefore, it has seen very Fig. 11 Comparison of results for x u−d from the summation method and from the midpoint (t ins = t sep /2) in the ratio method. Example data shown for the same ensemble as in Fig. 1 widespread use in lattice calculations, often for crosschecking results obtained from other methods, but also for obtaining final results. On the other hand, the summation method tends to produce larger statistical errors in the matrix element of interest compared to the ratio method; cf. Fig. 11. Besides, the excited state suppression in the summation method still depends on the values of t sep used. In actual lattice calculations the results is often dominated by the smallest sourcesink separations that enter the linear fit because the effective statistics rapidly deteriorates with increasing t sep if the number of measurements is kept constant independent of t sep . Since a similar issue arises with other methods, it may be desirable to spend extra computational effort to keep effective statistics constant, i.e. scale the number of measurements at different source-sink separations so as to achieve (approximately) constant statistical error. This has been done in some recent lattice studies, see e.g. Refs. [44,80]. Note that even for the ratio method it is not sufficient to check convergence in the presence of exponential error growth by inspecting the behavior of the observable as a function of t sep because the signal is typically lost before the asymptotic value is reached.

Feynmann-Hellmann inspired approach
Another realization of the summation over the operator insertion that appeared in the original paper for the summation method [23] has been revisited and used a few times in the last years [24,25,48,81,82]. Most recently it has been employed in a dedicated calculation of the axial charge [26,83], leading to a precise estimate of g u−d A = 1.271 (13) in agreement with the experimental result. The method is inspired by the Feynman-Hellmann theorem that relates the energy shift due to a perturbation in the action for a local bilinear operator O(x) with Dirac structure Γ and some parameter λ to the matrix element of a state |k Following the presentation in Refs. [24,83] the spectral decomposition of the nucleon two-point function in Eq. (11) at takes the following form where we have again assumed the source to be at t i = 0 and |Ω λ denotes the modified vacuum in the presence of the perturbation. The effective energy in Eq. (4) where the leading excited state behavior from the smallest energy gap Δ has been indicated by the second term in the last line. The correlation function ∂ λ C 2pt λ ( p, t) = ∂ ∂λ C 2pt λ ( p, t) can be straightforwardly computed w.r.t. the original vacuum at λ = 0 by replacing one of the propagators in a standard two-point function by a so-called Feynman-Hellmann propagator [24] S(y, x) = z=(t z ,z) i.e. a sequential propagator that is summed over the insertion time. Note that unlike in the previously discussed version of the summation method, here the summation over the operator insertion time is performed over the entire lattice, which is also the result originally derived in Ref. [23]. An advantage of this approach over the sequential method is that it requires only the computation of two-point functions which depend on the source-sink separation t = t sep but not on the insertion time. Therefore, it does not require new inversions for each value of t sep , and the computation of the usual sequential propagators for three-point functions is treated for the computation of the sequential propagator in Eq. (45). This mitigates some of the signal-to-noise problem when performing a dedicated calculation for a single observable. However, it allows only to compute results for a single operator insertion and a single momentum transfer at a time, which changes the cost comparison in favor of the standard method when computing multiple observables or especially for computing form factors at non-zero Q 2 .
In Fig. 12 a comparison of results for g u−d A is shown for the Feynmann-Hellmann inspired approach (left panel) used in Ref. [26] and the sequential method (right panel) from an analysis in Ref. [86] on a common ensemble, but using a different action in the valence sector. In both cases the final result has been obtained from a two-state fit to the data; see also Sect. 5. While the results are compatible, the statistical error in the left panel is smaller by roughly a factor of two. In Refs. [26,83] it has been argued that a key advantage of the Feynmann-Hellmann inspired approach is that the leading excited state contamination is not just of O e −Δt sep as for the commonly used version of the summation method but that it is further suppressed by the difference between contributions separated by τ as given in Eq. (44). This is used to justify the inclusion of data at much smaller values of t sep in these fits than for the sequential method, which causes the final statistical errors to be smaller. In fact, the fit in the left panel of Fig. 12 includes values of t sep 0.3 fm, while for the sequential method in the right panel only data for t sep /a 1.0 fm has been used; see also the discussion in the 2019 FLAG review in Ref. [87]. Note that in principle it is possible to perform a similar fit for the summation method by including the leading excited state contamination in Eq. 40, which would allow to include of data at smaller values of t sep as well; see the discussion in Sect. 5.1. However, the analysis and the resulting error estimate for g u−d A published in ref. [26] has been subject to criticism in Ref. [88] by the PNDME collaboration (re-)analyzing an extended set of ensembles and in Ref. [62] because the distribution of fit qualities fails the Kolmogorov-Smirnov test, unlike the updated results from the sequential method published in Ref. [88]. Note that there has also been disagreement regarding the central value of g u−d A obtained in Ref. [26] and the result from the sequential method in Refs. [86,88], which has been attributed to the chiral and continuum extrapolation as discussed in Ref. [88] leading to the claim that an overall error of ∼ 5% for g u−d A is realistic. A more detailed assessment of this claim is beyond the scope of this review.

Multi-state fits
A widespread approach to deal with excited states in nucleon structure calculations are fit models that explicitly include the effects of some of the excited states in the determination of the ground state NME. Since fits can be applied to the same data as the ratio and the summation method they are The vertical gray bands mark data excluded from the fit and the blue, green and red bands highlight the values of t sep /a = 10, 12, 14 that have been used for the data from the sequential method shown in the right panel. The curves in the right panel correspond to a standard twostate fit ansatz; for further details see the original analysis published in Ref. [86]. In both plots the horizontal line indicates the result for the ground state matrix element. The figure is reproduced with kind permission of the authors of Ref. [83] straightforward to implement and are currently used in most nucleon structure calculations. Furthermore, fits are quite flexible with respect to the actual fit ansatz, the choice of fit ranges and priors as well as e.g. the choice of simultaneously fitted observables. Therefore, basically every major lattice collaboration performing nucleon structure calculations has adopted their own favored method in the last years. In this section several approaches will be discussed that have been used in modern lattice calculations.

Two-state truncation
Fit models for excited state contamination in NME calculations are based on truncations of the spectral decomposition of the individual two-and three-point functions in Eqs. (11) and (12). A popular choice are the respective twostate truncations that are given by and where the superscripts i and f indicate the dependence of the overlap factors and energies on the initial and final state momenta p i and p f , e.g.
In principle, all the overlap factors, matrix elements and energies in this expression are free parameters that need to be determined by the fit. Correlated fits are performed to the lattice data by minimizing where C denotes the covariance matrix and χ = Y − f (X 0 , . . . , X n ) the difference between the vector of lattice data Y entering the fit and the fit model f (X 1 , . . . , X n ) depending on a set of fit parameters X 1 , . . . , X n that are to be determined by the fit. Although it is possible to perform a direct, simultaneous fit to the model given by Eqs. (46) and (47) in its most general form, see e.g. the nucleon form factor calculations in Refs. [73,89], in practice it can be difficult to obtain reliable results in this way considering the finite statistical precision of lattice data. The reason for this are the fairly large number of parameters even in case of the two-state truncation and possibly the size of the covariance matrix. This is particularly true at non-zero momentum transfer because the number of independent fit parameter is larger than for the zero-momentum case while the signal quality deteriorates with increasing Q 2 . Therefore, additional assumptions and knowledge about parameters are often used to stabilize fits. A common way to increase stability of the fit are assumptions on the energy gaps between the first excited state and the ground state. A rather basic ansatz has been explored several years ago by the Mainz group for a calculation of isovector electromagnetic form factors in Ref. [79]. Instead of fitting Right panel: results for including all terms in Eq. (33). The data for both figures has been generated in the study in Ref. [47] two-and three-point functions separately, the corresponding two-state truncations have been used to parametrize the leading excited state contribution to the effective form factors obtained from the ratio Eq. (13) as In this expression the energy gaps have been fixed to M π and 2M π assuming that multi-particle states are responsible for the leading excited state contamination while neglecting nucleon-pion interactions. The choice of M π for the first gap is motivated by the initial state nucleon having momentum p i allowing for a corresponding N π state with a moving nucleon and a pion at rest. On the other hand the final state nucleon is produced at rest, i.e. p f = 0 which implies that the lowest state for the second gap is a N ππ state with two pions in an S-wave, motivating the choice of 2M π for this gap. While this approach helps to stabilize the fit to lattice data, it comes at the price of introducing additional model dependence.
In a more recent study of nucleon charges and moments of twist-2 operators in Ref. [47] similar fits to ratio data have been carried out, which amounts to fitting the expression in Eq. (33). Although the number of fit parameters is reduced for the analysis of data at Q 2 = 0, this alone turned out insufficient to yield stable results when fitting the data for a single observable. The same problem has been observed in another calculation of nucleon charges in Ref. [90] that used single observable fits to the two-state truncated ratio among other methods. In order to remedy this issue, simultaneous fits to all six observables, i.e. g u−d A,S,T , x u−d , x Δu−Δd and x δu−δd have been introduced in Ref. [47]. This exploits correlations in the data which was found to further stabilize the fits and improve statistical precision without the need resorting to fur-ther assumptions or priors. However, one must keep in mind that the covariance matrix entering Eq. (48) can become a limiting factor for such simultaneous fits, as it may be poorly estimated for the available number of independent measurements if too many observables (or timelices) are included. Still, for the two-state truncation fitting several observables was sufficient to track the convergence of the gap as a function of the lower value of the fit range t fit , t sep /2 on plateau data symmetric around t ins = t sep /a. Results from this procedure are shown in Fig. 13 on an ensemble with a pion mass of M π = 200 MeV. In the left panel the last term in Eq. (33) has been neglected and clear convergence towards the expected gap of Δ ≈ 2M π is observed for M π t fit 0.4. Unlike approaches that use information from e.g. fitting the two-point function to fix the gap in the fit of the ratio or two-and three-point functions, this gives additional confidence that the fit model is applied in a regime for which ground state dominance is reached and higher states beyond the leading contamination become sufficiently suppressed.
The right panel of Fig. 13 shows results for Δ from a fit including the excited-to-excited state contribution in the last term of Eq. (33). This leads to significantly larger errors on Δ and to larger values of Δ itself although at least the trend towards 2M π remains visible 1 . The primary reason for this is again related to the signal-to-noise problem, i.e. the additional parameter A 11 in the fit is rather poorly constrained by the lattice data due to the exponential error growth with t sep and the limited number of different t sep values that can be included in the fit. Similar to what has been discussed before for the linear fits in the summation method that tend to be dominated by data at the smallest values of t sep , the situation could be improved by increasing the number of measurement for every step in t sep such that the effective statistics remains At any rate, going beyond zero-momentum transfer with this method without additional assumptions would remain difficult due to the increasing number of parameters and less precise data.
A different approach that employs a two-state truncation on the summation method [91] has been explored in Refs. [90] to complement the usual ratio and summation methods as well as the aforementioned (single observable) two-state fits to the ratio. The method corresponds to an extension of the usual summation method by supplementing the linear behavior of the summed ratio in t sep as given in Eq. (40) with the first excited state correction. This leads to the following fit form where the operator-dependent coefficients c X,i μ 1 ...μ n and the operator-independent gap Δ are left as free parameters of the fit. Exemplary results for the unrenormalized, axial charge as a function of the minimal source-sink separation entering the fit are shown in Fig. 14. Regarding the central value there is no dependence observed on the minimal value of t sep within the rapidly growing error. This kind of fit is similar to the ones performed in the determination of g u−d A using the Feynman-Hellmann inspired approach as discussed in Sect. 4.2, i.e. from a purely technical point of view the two calculations mainly differ in how the summation over the operator insertion has been implemented for the data entering the fit. However, the gauge ensembles used and most likely the resulting computational cost are very different, which makes a conclusive comparison impossible, even if one were to ignore the fact that in the sequential method multiple observables and momenta are computed simultaneously. It is interesting to note though that in both studies similarly small statistical errors are obtained for including data starting from t sep 0.3 fm. Still, the preferred value quoted in Ref. [90] from this method has been chosen such that only data with t sep 0.58 fm have been included due to concerns regarding the statistical quality of the fit including data at smaller values of t sep . In a comprehensive study of results obtained from different methods, this lead to a less favorable statistical signal quality for this method compared to other approaches. Since the gap Δ in Eq. (50) is independent of the observable it might be worth investigating if a simultaneous fit across multiple observables helps to improve stability and statistical precision of the fit as it was found for fits to the ratio in Ref. [47].
Instead of explicitly fixing an energy gap using model assumptions, or leaving it entirely as a free parameter as in the previously discussed approaches used in calculations of NMEs at zero-momentum transfer, it is possible to use information on energies from a separate analysis of the nucleon two-point function. This arguably leads to milder model dependence than e.g. the explicit choice Δ = 2M π made in Ref. [79] while still enabling two-state fits of NMEs at non-zero momentum transfer. A very recent account of this kind of approach is found in an analysis of the decomposition of the proton spin and momentum fraction by the Extended Twisted Mass Collaboration (ETMC) in Ref. [56]. In their favored approach the final state is always produced at rest p f = 0 and the fit is split up in several steps, starting with the two-point function at zero-momentum to extract the nucleon mass E 0 (0) = M N . The mass is used to fix all other ground state energies E 0 ( p) from the dispersion relation E 0 ( p) = M 2 N + p 2 removing the parameter from the fits of the two-state function at non-zero momentum. Since M N can be computed with rather good precision the values from the dispersion relation at higher momenta are much more precise than actual lattice data at the same momenta, thus one may assume that this helps to stabilize the fit, which is likely the motivation for this choice. Subsequently, fits to the two-point function at non-zero momentum are carried out to determine the energy gap Δ( p) = E 1 ( p) − E 0 ( p) for all required values of p, as well as the (ratio of) overlaps. Finally, the remaining four parameters (i.e. the matrix elements) are extracted from a simultaneous fit to the ratio data for several values of t sep . The fit ansatz is given by plugging Eqs. (46) and (47) in Eq. (13) using the already known val-  (7) MeV, a = 0.0801(4) fm and T × L 3 = 128a × (64a) 3 . The figure has been originally published in Ref. [56] and is reproduced under the Creative Commons Attribution 4.0 International license ues of the remaining fit parameters to remove them from the final fit. Note that in another variation of this method used in e.g. a calculation of electromagnetic form factors in Ref. [44] the final fit as been applied directly to the three-point function instead of the ratio. This is in principle equivalent, but the ratio may be preferred as it is usually shown in figures instead of the actual three-point functions.
Some results from the procedure in Ref. [56] for the quark-connected contribution to the isoscalar average quark momentum fraction are shown in Fig. 15 including a comparison with the summation method and the resulting values of χ 2 cor /d.o.f.. While the fit ranges in t ins have been fixed to t ins ∈ t sep /2 − 5a, t sep /2 + 5a , the lowest source-sink separation entering the fit t low s has been varied and for the two-state fit a slight upwards trend is observed. The fit qualities improve with increasing value of t low s , indicating that at small t low s higher excited states likely are non-negligible. A similar, but slightly downwards pointing trend is observed for the summation method which yields higher values than the two-state fit at t low s 1.4 fm. For the final choice of t low s = 0.96 fm in the two-state fit there appears to be some tension between the predicted dependence on the source-sink separation and data from the ratio method. Besides, the (correlated) errors of the summation method and the two-state fit overlap only at the largest value of t low s , which gives a hint that results from the two methods may only converge at even larger values of t sep . Since neither the fit range in t ins has been varied nor the resulting energy gaps are given in Ref. [56] it is difficult to judge if the fit is performed in the regime where contamination due to higher states are in fact already negligible, or a larger result as favored by the summation method and indicated by the trend in both methods would be more trustworthy. At any rate, the example highlights the importance of carefully crosschecking results from multi-state fits with other methods, apart from only demanding a good fit quality.
Ideally, also the energy gaps should be compared against theoretical expectations. In particular, the gap determined by the fit should converge towards the lowest gap in the spectrum for sufficiently large Euclidean time separations of the data in the fit. While corresponding results for the fitted gaps are rarely included in the literature, results from an older calculation of isovector nucleon charges by the PNDME collaboration in Ref. [86] give a strong hint that convergence for the energy gap between ground and first excited state is difficult to achieve for this kind of two-state fits that use information on the energies from the two-point function. In this study a similar fit ansatz for the two-state truncations of twoand three-point functions has been used, albeit only at zeromomentum transfer. Results for the energies of the first two states from this approach are shown in Fig. 16. Clearly, the results for the energy E 1 of the first excited state lie systematically above the non-interacting N ππ energies. This is in line with what has been reported in Ref. [90] for the leading gap obtained from a two-state fit to the two-point function on two ensembles at physical quark mass. Note that depending on the box size a N π state might have even lower energy, however, for typical lattices the energies are usually close to each other. The authors of Ref. [86] also performed some variations with respect to the fit ranges, i.e. the fit labeled "case 4" starts at larger values of t and t ins for the two-and three-point functions, respectively. While this tends to increase the statistical errors as expected, there is only a minor shift (if any) towards smaller values of E 1 observed. Therefore, it must Fig. 16 Energies from a simultaneous two-state fit of two-and threepoint functions from a calculation of g u−d A,S,T in Ref. [86] as a function of (valence) M π . The filled and open symbols refer to the ground state energy E 0 and the energy of the first excited state E 1 , respectively. Results are shown for two different fits "case 1" and "case 4" which differ by the number of timeslices included when fitting two-and threepoint functions as defined in table V of Ref. [86]. Whenever two different results have been quoted on the same ensemble in Ref. [86], i.e. from using only high precision solves and the full AMA results, only the result with full statistics has been included be considered doubtful whether such two-state fits correctly describe the contribution from the excited states that need to be subtracted to identify the ground state.

Including additional states
All of the fits discussed in the previous subsection share as a common feature that they only take into account the contribution of the lowest excited state as they are based on the two-state truncation in Eqs. (46) and (47). In principle, this can be extended by including further states in the truncation. This is the approach favored by the PNDME collaboration which they have used in several recent NME calculations, e.g. for isovector nucleon charges in Refs. [88,92], moments of of twist-2 operator insertion in Ref. [93], as well as electromagnetic and axial form factors in Refs. [54,72,94].
As outlined in Refs. [72,88] the two-point function is generally fitted to the four-state truncation In order to stabilize these fits beyond the two-state truncation a (sequential) empirical Bayesian analysis with Gaussian priors [95,96] is carried out for the masses E k≥1 and amplitudes A k≥1 as described in Ref. [72]. Although the results for E k>1 and amplitudes A k>1 remain sensitive on both the priors and the lower bound t min in the fit, the four-state truncation is found to describe the data sufficiently well as can be seen  Fig. 17 for four ensembles covering several values of the pion mass and two lattice spacings. However, as discussed in Ref. [88] the resulting energy gaps E 1 (0) − E 0 (0) at zeromomentum transfer are again found to be mostly incompatible with the theoretical expectation of E 1 − E 0 ≈ 2M π as it has been the case for the two-state truncation. This observation corroborates the conclusion that fits to a single two-point function are unable to fully resolve the excited state spectrum. Note that this behavior differs from what has been found in the previously mentioned simultaneous, direct fits to the ratio introduced in Ref. [47] for which convergence towards 2M π is in fact observed. Therefore, while using masses and amplitudes obtained from the four-state truncation of the two point function in the fits of the three-point function can yield an effective description of the data, this procedure should not be expected to allow for a systematic elimination of the contribution from the lowest excited states. For the truncation of the three-point function several ansätze have been used in recent calculations by PNDME that are referred to as 2 * -, 2-and 3 * -fits. In their naming scheme, the 2 * -fits include only the ground state contribution ∼ 0, p f O X μ 1 ...μ n 0, p i and the 0 → 1 transition matrix elements, while the 2-fit includes the full two-state truncation. Furthermore, the 3 * fits incorporate all terms involving the 1 → 2 and 0 → 2 transition matrix elements. In a recent calculation of electromagnetic form factors and radii in Ref. [94] final results are mostly obtained from the 3 * fits, however, in other cases these fits turn out unstable and the results from the two-state fit are preferred, see e.g. Ref. [88].
In fact, any matrix elements beyond the ones in the 2 * -fits are found to be poorly constrained by these fits.
Including further states in the truncations of two-and three-point functions does still not guarantee the efficacy of multi-state fits for removal of excited states, which remains dependent on the observable, the available statistical precision and a careful choice of fit ranges and priors. In particular, for observables requiring non-zero momentum transfer like r u−d E,M and μ u−d N , the excited state contamination is Q 2dependent For example, the result for the nucleon electric radius r u−d E = 0.769(27)(30) fm in Ref. [94] was found to be 16% smaller than the experimental value. This may be explained by the fact that the excited state correction in G u−d E (Q 2 ) increases with Q 2 and its convergence is from above. If the multi-state fit does not sufficiently remove these effects for data at larger Q 2 that enter the extraction of the radius, this may lead to an underestimation of the radius. On the other hand for μ u−d M the excited state contamination is large at small values of Q 2 and convergence is from below, which likely explains the low result of μ u−d N = 3.939 (86)(138) if the multi-state fits fails to fully remove the contamination.
Anyhow, the underestimation of electromagnetic radii and the magnetic moment is a common feature of several recent lattice calculations, see e.g. Refs. [44,80,89,94], and it is difficult to judge if deviations from experimental results are dominated by residual excited state contamination alone. For example, for the electric form factor also finite size corrections are expected to be Q 2 -dependent and can thus be difficult to disentangle from excited state effects. Two studies by the PACS collaboration in Refs. [70,97] give a hint that finite size corrections may indeed play a role, as their results obtained on two fairly large physical volumes (M π L ≈ 6 and M π L ≈ 7.2) are in better agreement with experimental results, albeit with larger errors. Moreover, final results for r E,M and μ M may exhibit further dependence on the method used to extract them from lattice data, see e.g. Ref. [49], the available lattice momenta and the chiral, continuum and finite size extrapolations.
Another interesting application of this kind of fit ansatz is found in the recent study of axial form factors in Ref. [54] that has already been mentioned in Sect. 3.2. In this case the conventional multi-state fit strategy that uses information on the energy gaps from the two-point function, labeled by S 2pt in the notation of Ref. [54], fails badly to achieve ground state dominance. The reason for this is that the four-state fit applied to the two-point function does not resolve N π states which are responsible for a large excited state contamination in the pseudoscalar form factor as discussed at the end of Sect. 3.2. In Ref. [54] an alternative fit strategy labeled S A 4 has been adopted that uses information on this additional, low-lying energy gaps through a fit to the three-point function from an insertion of the temporal component of the Fig. 18 Energy gaps in units of M π determined from the two fit strategies S 2pt and S A4 in Ref. [54]. The dotted lines represent estimates for noninteracting N π states with back-to-back momentum (red line) and with the nucleon at rest (blue line). The figure has been originally published in Ref. [54] and is reproduced under the Creative Commons Attribution 4.0 International license axial vector current. The result for the relevant energy gaps obtained from both procedures as a function of the momentum transfer in integer units are shown in Fig. 18 together with the noninteracting estimates for N π states with back-to-back momentum and with the nucleon at rest. The noninteracting estimates agree rather well with the behavior of the energy gaps extracted using strategy S A 4 . Using the fit strategy S A 4 the lattice data for r PCAC and r PPD defined in Eqs. (38) (39), respectively, are found to be compatible with unity as shown in the left panel of Fig. 10.
A conceptually different approach that aims at modeling N π state contributions in a many-state fit has been investigated in Ref. [90] for the extraction of g u−d A,S,T . It relies on using the noninteracting energies for the first N π states to fix the energy gaps Δ n = E n − M N for states with relative momentum p = 2π n/L up to some cutoff value |n| 2 ≤ |n max | 2 ≡ n 2 max . This approximation may be justified by the observation in Ref. [63] that the deviation between interacting and noninteracting energy levels is small compared to Δ n itself. The fit ansatz used in ref. [90] reads assuming that all ground-to-excited state transitions enter with the same coefficient and that off-diagonal transitions between different excited states are negligible, while excited state contamination in the two-point function is considered important. In Fig. 19 results from this procedure are shown for the unrenormalized axial charge. For small values of the lower bound τ 0 in the insertion time a strong dependence on the choice of n max is observed as expected in the presence of a tower of N π states contributing to the signal. With increasing values of τ 0 the values appear to converge across different n max and the result is found to be consistent with other methods albeit with larger errors. However, in Ref. [90] it has been concluded that the strong dependence on n 2 max in modeling the excited states at small τ 0 likely make the approach unreliable.

Variational techniques
A conceptually different approach to tame excited state contamination is the variational method [98,99], which has become a standard tool in spectroscopy calculations as it allows to systematically remove the lowest excited state contributions [100]. The method is based on the computation of a matrix of correlations functions for a basis of N operators χ = (χ 1 (t), . . . , χ N (t)) T with suitable quantum numbers, and solving the generalized eigenvalue problem (GEVP) for t > t 0 and k ∈ [0, . . . , N − 1]. Energy levels are obtained at large t from the principal correlators λ k (t, t 0 ) ∼ e −E k ( p)t , while the eigenvectors v k (t, t 0 ) carry information on matrix elements. Note that solving the GEVP generally results in an t 0 )) on each timeslice and performing the state assignment going from timeslice t to t + 1 can be a non-trivial task particularly in the presence of an exponentially deteriorating signal-to-noise ratio; for a discussion of methods to sort the states see Ref. [101]. The most important feature of the variational approach is that ground state energies and matrix elements are improved with respect to the leading excited state contamination which now depends on the gap E N ( p)− E 0 ( p) to the N th state in the spectrum [100] instead of the smallest gap in the spectral decomposition of a single two-point function E 1 ( p) − E 0 ( p). Therefore, contamination from excited states is more strongly suppressed which can be systematically improved by adding more (independent) interpolating operators to the basis χ . Beyond spectroscopy the variational approach can be applied in calculations of hadronic matrix elements [81]. In this case the generalized eigenvectors that diagonalize the two-point correlation function matrix are used to construct optimized interpolating operators for the kth state, that are used in the computation of optimized three-point functions. For nucleon structure calculations an issue arises again from the dense spectrum of multiparticle states, particularly towards physical quark masses and in large volumes. In principle, each state to be removed in the variational approach requires an additional operator, which may require a very large basis. Furthermore, it is known from spectroscopy calculations [102] using distillation [37] that multiparticle operators must be included to systematically resolve multiparticle states, which for NME calculations is also supported by chiral perturbation theory [66]; see the discussion in Sect. 3.2. However, these are numerically expensive to implement for the computation of nucleon three-point functions. In fact, there has so far only been a single, published study [38] using distillation in the context of NME calculations 2 , which will be discussed in some more detail in Sect. 6.3. Although nonlocal operators have been used in this calculation, multiparticle operator have still not been included.

Fig. 20
Effective form factor of the isovector axial charge for three different smearing levels and a variationally improved interpolator. Results from smeared interpolators are denoted by sm32, sm64 and sm128 corresponding to N smear = 32, 63, 128. For the variationally optimized interpolator t 0 and Δt refer to the sink times t = t 0 and t = t 0 +Δ used in the construction of the optimized two-point correlator, see Ref. [105] for technical details. The figure has been originally published in Ref. [105] and is reproduced under the Creative Commons Attribution 4.0 International license

Smeared interpolators
Since building a basis that systematically accounts for the lowest-lying excited states has not been attempted for NME calculations due to the computational cost induced by the need for multiparticle operators, all existing NME studies using the variational method have aimed at constructing an improved interpolator from a basis of computationally affordable interpolators. A straightforward way is to use differently smeared operators, which has been explored in Refs. [104][105][106]. However, it remains an open question if this approach is actually beneficial compared to the sequential method using a single interpolator with properly tuned smearing in combination with other methods for controlling excited states, e.g. two-state fits or the summation method. This concerns the resulting suppression of excited states as well as the computational cost required to achieve a given target precision. Figure 20 illustrates the first part of this issue for g u−d A : the choice of smearing steps N smear has a large impact on the resulting excited state suppression for the single interpolator approach and careful tuning of the smearing is required for a meaningful comparison, because too small a value of N smear leads to enhanced excited state contamination. Although the excited state contamination appears to be slightly more reduced for the variational method in this example from Ref. [105], it is close to the one from the single interpolator approach for the largest value of N smear and a further increase of N smear might have resulted in a compatible value.
While in Refs. [105,107] some indication has been presented that the variational method is more robust with respect to excited state contamination than the two-state fit and summation method, this has been challenged by the study in Ref. [106]. In this study it has been concluded that the efficacy of the variational method and a two-state fit applied to the single interpolator approach for reducing excited state contamination is similar. In particular, it has been stressed that the comparison to the two-state method in Refs. [105,107] did not use the optimal smearing size and that the resulting statistical errors might be artificially large due to the choice of the values of t sep in the fit. Furthermore, it has been pointed out in Ref. [106] that the computational cost depends on the setup and that lighter quark masses and larger temporal lattice sizes T work in favor of a (properly tuned) two-state fit approach. At any rate, it appears fair to conclude that no clear advantage of the variational method using a basis of smeared interpolators over other commonly used methods has been found in existing studies.

Generalized pencil of function (GPOF)
Another, very cost-efficient approach to obtain a variational basis is the so-called generalized pencil of function method, that has been applied for NME calculations in Refs. [108][109][110][111]. It offers an alternative way to analyze the data that are computed using the sequential method for several sourcesink separations, as required for other methods such as twostate fits or the summation method. In the GPOF method a set of linearly independent interpolating operators is created by time-shifting an existing operator χ(t) and used to build a (N + 1) × (N + 1) matrix of two-point functions where any momentum dependence has been suppressed in the notation. Typically, for NME calculations only a single, original operator is used, e.g. the standard interpolator in Eq. (10). For the case of a nucleon three-point function a non-symmetric but otherwise similar correlation function matrix is built by applying the same procedure to its respective source and sink operators where all further indices and parameters besides the timedependence have been suppressed. The matrix of eigenvectors obtained from the GEVP applied to C 2pt (t) is then used to diagonalize the three-point function matrix The principal correlator Λ 0 (t ins , t sep ) for the ground state can be treated in the same way as the standard three-point functions in a subsequent analysis, e.g. the ratio method. Anyhow, NME calculations are just a particular application of the GPOF method and it can in fact be applied to any problem in hadron spectroscopy and structure calculations whenever a variational basis is desired. Further details on GPOF based approaches and the relation between several families of methods, i.e. the GEVP and GPOF, the Prony method [112][113][114][115] and the Gardner method [116] have recently been discussed in Ref. [101], where also a new combination of GEVP and GPOF has been proposed. In practice, the signal-to-noise problem and availability of different values of t sep limits the application of the GPOF in NME calculations to a single additional interpolator and thus a 2 × 2 correlation function matrix. Figure 21 shows corresponding results from a 2 × 2 GPOF with δt = 2a for the nucleon effective mass on an ensemble at physical quark mass. The excited state contamination is found to be drastically reduced in the GPOF principal correlator at small values of t/a. However, the point errors for the GPOF ground state are larger than from the single. As pointed out in Ref. [110] this behavior is expected, because for a system with exactly two-states with energies E 0,1 the two-and three-point ground state principal correlators read and respectively, implying that the statistical uncertainties are dominated by the correlators at the largest value of t sep . This strongly limits the statistical precision for the three-point case due to the signal-to-noise problem, at least if effective statistics are not kept constant for increasing values of t sep . Still, for the nucleon mass this effect is compensated by the much earlier onset of the plateau at ∼ 0.4 fm. Furthermore it has been pointed out in Ref. [110], that the GPOF might not be efficient for contamination from transition matrix elements. This is because they are more strongly suppressed in the two-point function than the corresponding contributions to the three-point functions by transition to the ground state. However, within the statistical precision of that study it has been concluded that GPOF and summation method lead to comparable results. This has also been confirmed by a more recent study in Ref. [111] of moments of twist-2 operator insertions, which are notorious for large excited state contamination. In the context of NME calculations the main advantage of the GPOF is clearly its simple implementation compared to e.g. multi-state fits and that it can usually be applied on existing data without adding computational cost. However, despite its simplicity the method has not seen widespread use in recent NME calculations. Still it might be worthwhile exploring the GPOF approach for NME calculations further, particularly if the required data are readily available from other methods such that it can e.g. serve as a crosscheck.

Hybrid interpolators
The search for an affordable basis of interpolators has recently been extended to so-called hybrid interpolators that include an insertion of a chromomagnetic field B i = − 1 2 i jk F jk and that can be implemented at similar computational cost compared as the standard approach. In Ref. [117] the following basis of operators has been studied Fig. 22 Comparison of results for effective M N from a low-statistics smearing scan using the interpolator χ 1 defined in Eq. (63) (open, gray circles), the variationally improved operator χ opt from a second tuning run and results for χ 1 (filled, green circles) and χ opt (blue diamonds) from full statistics. The figure has been originally published in Ref. [117] and is reproduced under the Creative Commons Attribution 4.0 International license whereũ,d are smeared quark fields, P i j = δ i j − γ i γ j /3, and χ 1 is the standard interpolator in Eq. (10) up to an additional (positive) parity projector P + = (1 + γ 0 )/2. Two tuning runs have been performed to optimize the smearing for χ 1 and construct a variationally improved interpolator χ opt for the ground state as introduced in Eq. (56). Some results for the nucleon effective mass from Ref. [117] are shown in Fig. 22. The variational basis is found to lead to a reduction of excited state contamination compared to the standard approach, albeit less than the reduction that is usually observed in the previously discussed GPOF approach. Furthermore, these results once more demonstrate the impor-tance of tuning the smearing for optimal results and a meaningful comparison to other methods. However, for the actual NME calculation the results found in Ref. [117] were mixed, as can be seen from the results for g u−d A,T in Fig. 23. While the excited state contamination is significantly reduced for g u−d T using the variationally improved interpolator χ opt , it is actually increased for g u−d A . The latter result might be explained by a partial (accidental) cancellation of excited states present in the single interpolator approach, which is weakened or prevented in the variational approach. Besides, the statistical errors for χ opt are always larger than for χ 1 . These findings corroborate the conclusion that a small variational basis without multi-particle operators does not allow for a systematic treatment of excited states in nucleon structure calculations.
Similar results have been found previously in Ref. [38]. In addition to a standard, Jacobi smeared [118] interpolating field and a single distilled operator 2 S S 1 2 + that resembles the nucleon interpolator in Eq. (63), two bases of distilled operators B 3 , B 7 have been constructed using covariant derivatives to obtain variationally improved interpolators denoted byP 3 andP 7 , respectively. Basis B 3 includes 2 S S 1 2 + and two hybrid interpolators which have been found to have large ground state overlap in Ref. [119], while basis B 7 expands B 3 by four additional operators that probe the radial structure of the nucleon. Note that B 3 is similar to the basis used in Ref. [117] to obtain the variationally improved interpolator.  to what is observed in Ref. [117]. A second observation is that further increasing the basis indeed improves the situation as one should expect. Again, for g u−d T (not shown here) the excited state contamination has been found to be reduced by including hybrid interpolators when compared to the single interpolator approach. It should be further noted that the standard, smeared interpolator has been found to give results qualitatively consistent with using 2 S S 1 2 + as expected, albeit with much larger errors. However, it is not clear if the smearing has been tuned in Ref. [117]. Besides, it is not possible to directly compare computational cost in a meaningful way to distillation. Still, the findings of the two studies in Refs. [38,117] using hybrid interpolators qualitatively agree in that a small variational basis without multiparticle operators is insufficient to reliably deal with excited state contamination in NME calculations across multiple observables.

Parity-expanded variational analysis (PEVA)
Yet another way to build a variational basis has been developed and applied in a series of papers [120][121][122] for studies of baryons at non-zero momentum. Originally introduced in Ref. [120] at the level of nucleon two-point functions, the method aims at resolving parity mixing between nucleon states at non-zero momentum, hence its name "parityexpanded variational analysis". Unlike other approaches that attempt to find a generally applicable variational basis to construct improved operators, the PEVA approach is designed to specifically deal with potential contamination caused by this mixing with states of opposite parity.
The basic idea of this method is to start from a basis of conventional interpolators denoted by {χ i } and construct an extended basis using a helicity projector At zero-momentum the two operators transform as eigenstates of parity, i.e. χ i,± p → +χ i,± p and χ i,± p → −χ i,± p , respectively, while at non-zero momentum they do not have definite parity. The authors of Refs. [120][121][122] always consider an initial basis made up two interpolators with four different levels of Gaussian smearing to generate the smeared quark fieldsũ,d, effectively resulting in an 8×8 basis for the conventional variational approach and a 16 × 16 basis for the PEVA. In Ref. [121] the method has been applied to the calculation of electromagnetic form factors while the most recent study using the PEVA approach in ref. [122] extends this further to (elastic) form factors of the first two parityodd excitations of the nucleon and its lowest-lying parityeven excitation. Here we shall focus on some of the results from Ref. [121]. In this study the PEVA method has been first applied for G E (Q 2 ) for which no significant difference to results from the conventional approach has been observed. However, for the magnetic form factor PEVA was found to lead to additional reduction of the excited state contamination. This is shown in Fig. 25 where the ratio G Conv.  Results are shown separately for up-and down-quark contributions and have been obtained on an ensemble with M π ≈ 156 MeV, a = 0.0933 (13) fm and T × L 3 = 64a × (32a) 3 . The figure has been originally published in Ref. [121] and is reproduced under the Creative Commons Attribution 4.0 International license Fig. 26 Comparison of results for the isovector magnetic moment from the conventional approach (open symbols) and the PEVA method (filled symbols) across several ensembles. The data have been finite-volume corrected using the corrections in Ref. [123]. The values at physical quark mass have been obtained from a chiral extrapolation. The figure has been originally published in Ref. [121] and is reproduced under the Creative Commons Attribution 4.0 International license ensemble with the lightest quark mass in this study, which can be inferred from Fig. 26 showing results for the isovector magnetic moment as a function of M 2 π . Note that the data in this figure have been corrected for finite volume effects using chiral perturbation theory [123] which is found to give a correction of ∼ 10% for the isovector magnetic moment.
Although it is not obvious why one would expect contamination by opposite parity excited states for the magnetic form factor particularly at lighter pion masses, the authors of Ref. [121] give some possible explanation for this fea-ture from an investigation of the negative-parity spectrum, i.e. they infer an increasing role of multiparticle states in the negative-parity spectrum at lighter pion mass. They further speculate that this could lead to a change in the coupling to the localized operators they use and thus affect the opposite-parity contamination in the ground state matrix element. However, if this explanation is correct, it is not a priori clear that the PEVA method would still yield a significant correction at light pion masses once multiparticle operators are actually included and it might in fact be more relevant for studies of parity-odd states as in the exploratory study in Ref. [122].

Summary and outlook
In the last few years studies of nucleon structure at physical quark mass have become feasible, thus systematic effects due to the chiral extrapolation can be considered well under control. Furthermore, discretization effects are empirically found to be rather small and continuum extrapolations are now part of state-of-the-art calculations as well. The later also applies to some extent to finite size extrapolations, although in this case the picture may be less clear. On the other hand, excited states remain arguably as the most important source of systematic uncertainty in current lattice simulations. One reason for this is that they are strongly observable-dependent and can be large, which in practice makes it difficult to come up with a generally applicable approach to treat them. Moreover, excited state contamination is intimately related to the signalto-noise problem, which is why there is usually a trade-off between statistical precision and accuracy related to this systematic. This is problematic as it can easily lead to an underestimation of the overall error.
In this review an overview of approaches to mitigate excited state contamination in contemporary lattice QCD calculations of nucleon structure has been presented. These methods can be roughly divided in three categories, i.e. summed operator insertions, multi-state fits and applications of the variational method. While direct, quantitative comparisons of different methods taking into account the computational cost are rarely found in the literature and often difficult or even impossible to perform, each of these methods has its individual advantages and shortcomings that can be summarized as follows: 1. Methods based on the summation over the operator insertion aim at increased suppression of excited states. At least in the commonly used form this approach is simple to implement, rather robust and less affected by human bias than e.g. multi-state fits, because there are not many tunable parameters involved. However, statistical errors are typically larger than for (naive) fits which might explain why the summation method has often been merely used as a crosscheck for other methods. Moreover, there is no a priori guarantee that the additional suppression of excited states is indeed sufficient for a given observable and target precision.
2. Multi-state fits have become rather sophisticated over the last years and are by now predominantly used to obtain final results. They aim at explicitly correcting lattice data for the leading excited state contamination. Their main advantages are flexibility with respect to the choice of fit ansatz and -in principle -the possibility of tracking convergence of the resulting energy gaps. In practice, most implementations use information from the two-point functions to determine energy gaps and make heavy use of priors if more than two states are included, making them more prone to human bias. Moreover, fits naively relying on information on energies from two-point function should be considered unreliable as the fitted energy gaps notoriously fail to converge to the theoretically expected lowest-lying state and it has been shown that these fits may indeed miss the lowest-lying (multiparticle) states which can lead to large residual excited state contamination as observed for e.g. axial form factors.
3. The variational method is potentially the most powerful approach as it allows to systematically extract states and suppress excited state contamination in a given matrix element. However, in current implementations the efficacy of the method remains limited and strongly observabledependent as only fairly small bases of operators have been used and multiparticle operators have not been included at all. The last point causes an issue similar to the one for multi-state fits using information on energy gaps from two-point functions, i.e. a small basis particularly without multiparticle may simply miss certain low-lying excited states. Still, such states can yield significant residual contamination depending on the matrix element.
Since no approach is clearly favorable across multiple observables, it remains crucial to crosscheck whatever method is used to obtain final results and perform a careful assessment of the residual excited state contamination. In particular, it is insufficient to claim agreement between two methods if one of them has much larger errors while the result with the smaller statistical error is quoted as the final estimate. Ideally, one should study correlated differences and assign a systematic error due to residual excited state contamination. If multi-state fits are used on a sufficient number of ensembles, non-parametric criteria may be used to further test the plausibility of results. Additional care is required when using methods that are likely to be affected by human bias due to a large number of free parameters or that make model assumptions, such as using the spectrum from the twopoint functions to determine energy gaps of the three-point function.
One way to systematically improve over existing analyses of excited states would be the inclusion of multiparticle operators in a variational analysis. This should yield a more reliable resolution of the spectrum and thus lead to a more comprehensive and observable independent removal of excited state contamination provided that a large enough basis is used. However, such an extension would still be be demanding from a computational as well as technical point of view. Another possibility for future improvement would be a solution or mitigation of the exponential signal-to-noise problem itself. This might be achieved by using multilevel methods as introduced in Refs. [124,125] which have recently been further explored in the context of the hadronic vacuum polarization contribution to g − 2 [126].