Opposite-Parity Contaminations in Lattice Nucleon Form Factors

The recently-introduced Parity Expanded Variational Analysis (PEVA) technique allows for the isolation of baryon eigenstates at finite momentum free from opposite-parity contamination. In this paper, we establish the formalism for computing form factors of spin-1/2 states using PEVA. Selecting the vector current, we compare the electromagnetic form factors of the ground state nucleon extracted via this technique to a conventional parity-projection approach. Our results show a statistically significant discrepancy between the PEVA and conventional analyses. This indicates that existing calculations of matrix elements of ground state baryons at finite momentum can be affected by systematic errors of ~20% at physical quark masses. The formalism introduced here provides an effective approach to removing these systematic errors.


I. INTRODUCTION
Lattice QCD investigations of baryon matrix elements form a rich and varied field of study. In such investigations, it is necessary to isolate the state of interest from the tower of energy levels observed on the lattice. For ground states, this can be achieved through simple Euclidean time evolution, but it is common to use more advanced techniques such as multi-exponential fits, the summation method, or variational analysis in order to extract the signal of interest from earlier time slices and avoid the noise present in the tails of correlation functions. Regardless of the specific technique used, when working with spin-1 /2 baryons it is typical to perform a simple zero-momentum parity projection to the parity sector of interest, significantly reducing the number of possible contaminating states before even beginning the analysis. This technique works perfectly for at-rest baryons, completely removing opposite-parity contaminations, but for matrix elements where at least one of the initial or final state is boosted to non-trivial momentum this is no longer the case. Since eigenstates with nonzero momentum are not eigenstates of parity, the parity sectors are no longer well defined and a naïve parity projection allows opposite-parity states to re-enter the correlation functions.
In Ref. [1], we introduced the Parity Expanded Variational Analysis (PEVA) technique to address this issue, applying it to the extraction of the nucleon spectrum from two point correlation functions. In this paper, we extend the formalism to the computation of matrix elements from three-point correlation functions, and apply it to the specific example of calculating the Sachs electric and magnetic form factors G E (Q 2 ) and G M (Q 2 ) of the ground-state proton and neutron. The Sachs form factors describe the response of a baryon to the vector current. At low Q 2 , these form factors give information about the large-scale electromagnetic structure of the state, such as its charge radius and magnetic moment; at high Q 2 they give information about the short-distance internal structure of the state. These form factors can be determined experimentally to high accuracy. Computing them in ab-initio lattice QCD provides an important confrontation of theory with experiment.
In addition, computing these form factors on the lattice gives us important insight into the underlying physics. For example, on the lattice it is possible to separately compute the contributions to the form factors from connected diagrams (as studied in this paper) and disconnected diagrams, giving insight into the role sea quarks play in the structure of the proton and the neutron. One can also alter the electric charges of the quarks, readily illustrating the quark-flavour structure of the nucleon.
We probe the values of these form factors by creating an incoming nucleon on the lattice, having it interact with a vector current with some momentum transfer q, and then annihilating the outgoing nucleon with a fixed momentum p ′ . By momentum conservation, the incoming state must have momentum p = p ′ − q. Due to the way we include the vector current on the lattice, we only consider a small number of fixed momentum transfers q. By varying the three-momenta of the outgoing state and hence the incoming state, we gain access to the form factors at a range of (1) In particular, these boosts provide access to values close to Q 2 = 0, well below |q min | 2 = (2π/L s ) 2 , without requiring the use of twisted boundary conditions. By accessing a range of values, we gain insight into the Q 2 dependence of the form factors, and can make a comparison with various models and experiment. By studying the low-Q 2 dependence of the electric form factor, we can make an ab-initio determination of the charge radius of the proton. In addition, we observe that when considering the contributions from each quark flavour independently, G E (Q 2 ) and G M (Q 2 ) have a similar Q 2 dependence in the range considered. Hence, we can access the magnetic dipole moments of the proton and neutron by taking ratios of the quark-sector form factors.

II. PARITY EXPANDED VARIATIONAL ANALYSIS
The process of extracting a baryonic excited-state spectrum via the PEVA technique is presented in full in Ref [1]. We present here a brief summary of this process to introduce the notation and concepts necessary to describe the extension to the computation of baryonic matrix elements.
We begin with a basis of n conventional spin-1 /2 operators {χ i (x)} that couple to the states of interest. Adopting the Pauli representation of the gamma matrices, we introduce the PEVA projector Γ ±p ≡ 1 4 I + γ 4 I ∓ iγ 5 γ k p k , and construct a set of basis operators We note that we use a Euclidean metric δ µν , and hence there is no need to distinguish between contravariant and covariant indices.
We then seek an optimised set of operators φ ±p α (x) that each couple strongly to a single energy eigenstate α. These optimised operators are constructed as linear combinations of the basis operators where the sum is over both the primed and unprimed operators. The coefficients v α i (p) and u α i (p) can be found as the left and right generalised eigenvectors [2, 3] of G(p ; t + ∆t) and G(p ; t), where the correlation matrix with i and j ranging over both the primed and unprimed operators.
The choice of sign in the PEVA projector has no effect on the values of these two-point correlation functions. This is because the basis operators occur twice in the correlation function expression, for both creation and annihilation. The diagonal blocks of the Dirac matrix pick up zero or two minus signs. In the off-diagonal blocks, the sign changes for the primed operator in Eq. (2b) and the coupling to the parity mixing via ∓iγ 5 γ k p k cancel out. As a result, the coefficients for constructing φ +p α (x) and φ −p α (x) are identical, up to a choice of overall sign of the eigenvector. We choose the eigenvector sign to ensure the operators match at zero momentum, where the choice of ±p has no effect on the physics.
Using the optimised operators, we can construct the eigenstate-projected two-point correlation function

III. BARYON MATRIX ELEMENTS
A. General matrix elements In this section, we establish the formalism to extend the PEVA technique to the computation of baryon form factors. To perform the extension, we consider threepoint correlation functions, inspecting their Dirac structure to extract the signal of interest. We then take ratios with two point functions to remove the time dependence and cancel out dependence on the interpolator couplings. The calculations are performed in the most general kinematics that can be realised.
Due to a lattice Ward identity associated with the conserved current, the three-point correlation functions for the electric form factor normalised to unit charge must approach the two-point correlation functions exactly on a configuration-by-configuration basis as Q 2 − → 0. As a result, the two-and three-point correlation functions are highly correlated at low Q 2 . The ratios we take facilitate the cancellation of statistical fluctuations, significantly reducing the statistical uncertainties in our extracted form factors.
By performing a parity-expanded variational analysis as described in Sec. II, we construct optimised operators φ ±p α (x) that couple to each state α. We can use these operators to calculate the three point correlation functions where J (x) is some current operator, which is inserted with a momentum transfer q = p ′ − p. The consideration of G 3 − (J ; p ′ , p ; t 2 , t 1 ; α) (where the sink operator uses the opposite PEVA projector sign convention to the source operator) is required to optimise the extraction of the form factors for general kinematics. We note that it is sufficient to consider this change of projector for the sink operator alone, leaving the source operator as φ +p α (0) in all cases considered.
By inserting the complete set of states on either side of the current, and noting the use of Euclidean time and fixed boundary conditions (or negligible backward-running state contributions), we can rewrite this three point correlation function as Note, the formalism presented here assumes perfect state isolation such that each optimised operator couples only to a single state. We see that the time dependence of this three point correlator is entirely contained within exponentials of the energy, and the remaining structure depends on both the overlap of the optimised operator with its corresponding state and the matrix element for the current operator, α ; p ′ ; s ′ | J (0) |α ; p ; s . As we will see below, this Euclidean time dependence, along with the scalar factors relating to the coupling between the optimised operator and the state α can be removed by taking appropriate ratios with the two point correlation functions In this paper, we investigate the electromagnetic properties of the proton and neutron by choosing the current operator J (x) to be the vector current. In particular, we use the O(a)-improved [4] conserved vector current used in Ref. [5], where r is the Wilson parameter, and This current is derived from the standard conserved vector current for the Wilson action This choice of current operator gives the matrix element where 2 is the squared fourmomentum with the conventional sign, and the invariant scalar functions F 1 (Q 2 ) and F 2 (Q 2 ) are respectively the Dirac and Pauli form factors.
Hence, using Eqs. (9) and (14) we can rewrite the correlator as Using the spin sum the three-point function is To extract our desired signal from this spinor structure, we can take the spinor trace with some spin-structure projector Γ S . This trace is then called the spinor-projected three-point correlation function If we consider the function where the prime on F ′ ± (Γ S , J ) denotes the presence of the PEVA projectors, then we can express Eq. (18) as These spinor-projected correlation functions have a nontrivial time dependence, which can be removed by constructing the ratio [6] where Γ 4 = (I + γ 4 )/2 and Γ k = (I + γ 4 )(i γ 5 γ k )/2 form the basis for the spin projectors we use, and r µ and s µ are coefficients selected to determine the form factors. In addition, as the momentum transfer q − → 0, charge conservation requires that the temporal component of the three point correlator for the conserved vector current becomes exactly proportional to the two point correlator on each gauge field configuration, that is Because of this, taking the ratio in Eq. (21) facilitates the cancellation of statistical fluctuations in the two-and three-point correlators, providing results with small statistical uncertainties, at least in the case of G E (Q 2 ). We can then define the reduced ratio, Taking this reduced ratio and substituting in the expressions for the projected correlation functions, we obtain By investigating the r µ and s σ dependence of this ratio, we find that the clearest signals are given by where s is chosen such that p · s = 0 = p ′ · s , r is equal to q × s , and the sign ± in Eq. (25) is chosen such that 1 ± p · p ′ is maximised. This choice maximises the signal in the lattice determination of the correlation function ratios.
We can then find the Sachs electric and magnetic form factors, through appropriate linear combinations of R T ± and R S ∓ . A similar procedure can be applied to extract the relevant form factors from any current.
We have shown how the PEVA technique can be applied to the calculation of baryon form factors for arbitrary kinematics. Doing so ensures that these form factors are free from opposite parity contaminations.

IV. SACHS ELECTRIC FORM FACTOR
We now apply this technique to calculate the Sachs electric form factors of the proton and the neutron. This gives us insight into the distribution of charge within these states.
The results in this paper are calculated on the PACS-CS (2+1)-flavour full-QCD ensembles [7], made available through the ILDG [8]. These ensembles use a 32 3 ×64 lattice, and employ a renormalisation-group improved Iwasaki gauge action with β = 1.90 and non-perturbatively O(a)-improved Wilson quarks, with C SW = 1.715. We use five ensembles, with stated pion masses from m π = 702 MeV to 156 MeV [7], and set the scale using the Sommer parameter with r 0 = 0.4921(64) fm [7]. More details of the individual ensembles are presented in Table I, including the squared pion masses in the Sommer scale. When fitting correlators, the χ 2 /dof is calculated with the full covariance matrix, and the χ 2 values of all fits are consistent with an appropriate χ 2 distribution.
At the lighter pion masses considered, the relativistic components of the baryon spinor will be enhanced. As a result, the parity-mixing at finite momentum will be increasingly problematic. However, at lighter pion masses, the gauge noise is more significant, and can occlude the parity-mixing effects if the statistics are insufficient.
For the variational analyses in this paper, we begin with an eight-interpolator basis is formed from the conventional spin-1 ⁄2 nucleon interpolators with 16, 35, 100, or 200 sweeps [9] of gauge-invariant Gaussian smearing [10] with a smearing fraction of α = 0.7, applied at the quark source and sinks in creating the propagators. For the PEVA analyses, this basis is expanded to sixteen operators as described in Section II. Before performing the Gaussian smearing, the gauge links to be used are smoothed by applying four sweeps of three-dimensional isotropic stout-link smearing [11] with ρ = 0.1.
To extract the form factors, we fix the source at time slice N t /4 = 16, and (utilising the sequential source technique [12]) invert through the current, fixing the current insertion at time slice 21. We choose time slice 21 by inspecting the projected two-point correlation functions associated with each state and observing that excited-state contaminations are suppressed by time slice 21. This is evaluated by fitting the effective mass in this region to a single state ansatz verifying that the full covariance χ 2 /dof is satisfactory. We then extract the form factors as outlined in Section III for every possible sink time and once again look for a plateau consistent with a singlestate ansatz.
Performing the sequential source technique through the current requires us to choose our current operators and momentum transfers at inversion time. However, this allows us to vary the sink momentum, and by extension the source momentum, as well as varying the form of the interpolation functions at the sink. This gives us access to several states, as well as a range of values of In particular, values of Q 2 well below that encountered in the frames with p, or p ′ = (0, 0, 0) are accessed via kinematics such as p = (1, 0, 0), p ′ = (2, 0, 0). The main alternative approach to accessing Q 2 values in this region is to use twisted boundary conditions to change the momentum discretisation, allowing the valence quarks to take different momentum values from the sea quarks.
In our approach, all momenta attained by the valence quarks are present in the sea, and thus we avoid the complexities of partial quenching effects inherent in the twisted boundary approach. Table II summarises the kinematics considered herein. To begin, we inspect the Euclidean time dependence of G E (Q 2 ), extracted as outlined in Section III. We consider independently the connected contributions to G E (Q 2 ) from single valence quarks of unit charge. The two flavours considered are the doubly represented quark flavour, or the up quark in the proton (u p ); and the singly represented quark flavour, or the down quark in the proton (d p ).
In the case of perfect optimised operators, there should be no Euclidean time dependence, and the extracted form factors should be perfectly constant (up to statistical fluctuations) after the current insertion. However, in practice a finite operator basis is insufficient to perfectly isolate each state, leading to residual excited-state contaminations. These show up as enhanced or suppressed form factors at early Euclidean times. In light of this, care must be taken to select a Euclidean time region in which these excited-state contaminations are suppressed and the single state ansatz is satisfied.
In Figs. 1 and 2 we plot both PEVA and conventional extractions of G E (Q 2 ) with respect to Euclidean sink time at the lightest pion mass of m π = 156 MeV and the lowest-momentum kinematics of p = (0, 0, 0) and p ′ = (1, 0, 0). We see that starting from time slice 22, Table I. Details of the gauge field ensembles used in this analysis. For each ensemble we list both the pion mass given in Ref. [7], with the lattice spacing set by hadronic inputs, and our determination of the the squared pion mass with the lattice spacing listed in the table, which is set by the Sommer parameter with r0 = 0.4921(64) fm [7].  Table II. Different kinematics used in our analysis to access a range of Q 2 values. The Q 2 value listed is for the ground-state nucleon at the middle pion mass of mπ = 411 MeV. The statistical error listed for Q 2 comes from both the determination of the mass of the state and the conversion to physical units.  which is immediately after the source, both extractions of G E (Q 2 ) are quite flat across all time slices considered. However, the errors on G E (Q 2 ) are sufficiently small to identify a small Euclidean time dependence at early time slices. We find that this dependence is suppressed by time slice 24 and are able to find a clear and clean plateau from 24-27 for both extractions. For both quark flavours considered, there is no significant difference in the fit ranges, extracted values or errors between the two extractions.
The conventional thought is that the opposite-parity contaminations are small. Because they are from heavier states, these contaminations are suppressed by Euclidean time evolution. We will see that this is not the case for the magnetic form factor, where we will present direct evidence for important opposite parity contamination. For the electric form factor, agreement between the PEVA and conventional analyses can be maintained provided the opposite-parity contaminations contribute to the form factor in a manner similar to that of the parity sector under examination. As a result, despite their continued presence, the opposite-parity contaminations do not significantly perturb the value of the electric form factor.
At the heavier pion masses, the statistical noise in the form factor extractions decreases, and the plateau region   shifts somewhat. However for all five masses, the qualitative behaviour described above remains true, save for the following anomalies. At m π = 570 MeV, the plateaus from PEVA start two time slices earlier than those from the conventional analysis. For example, Fig. 3 shows the plateaus for the doubly represented quark flavour. This is potentially a signal of opposite-parity contaminations entering into the analysis. However, there is no statistically significant difference in the fit values from the two methods and the different plateaus do not show up at   ) dp (PEVA) dp (Conv.) Figure 5. Contributions of dp to the ground state GE(Q 2 ). Pion mass and kinematics are as in Fig. 4 above. The conventions used in this plot are the same as in Fig. 1. The PEVA fit starts from time slice 23, but the conventional analysis starts from time slice 24. any of the other masses considered, so it is inconclusive. We can also consider changing the momenta of the initial and final states, both by changing the momentum transfer, and by boosting both the initial and final states without changing the three-momentum transfer.
If we do this for the m π = 156 MeV ensemble, where we previously found consistent plateaus between PEVA and a conventional analysis, we find some discrepancies. In Figs. 4 and 5, we boost the initial state momentum to p = (−1, 0, 0) and increase the momentum transfer to q = (2, 0, 0), leading to a significant increase in Q 2 . In this case, we find that the PEVA plateaus start one time The y-axis intercept is fixed to one, as we are using an improved conserved vector current and the quarks are taken with unit charge. The errors on these fits are small enough that the shaded bands are barely distinguishable from the central lines. The fits correspond to a charge radius of 0.684 (19) fm for the doubly represented quark (up) and 0.659(21) fm for the singly represented quark (dp).
slice earlier than the conventional plateaus. They have consistent plateau values, but due to the earlier onset of the PEVA plateaus the statistical error is reduced. These results suggest that there are contaminations in the extraction of G E (Q 2 ) with the conventional analysis at this mass. However the differences are not consistent across all higher-momentum kinematics, and are not enough to categorically ascribe these problems to opposite-parity contamination. We do note that when there is a difference in the onset of the plateaus, the PEVA plateau always starts earlier.
For the other four masses, almost all kinematics have identical plateaus in G E (Q 2 ) from both analyses, save for m π = 570 MeV, which once again has consistently earlier plateaus for PEVA than the conventional analysis. It is unclear why m π = 570 MeV has inconsistent plateaus at a range of kinematics when the other three heavy masses don't. However, it is clear that whatever opposite-parity contaminations are occurring, they are not affecting the G E (Q 2 ) values extracted, at least within our current statistical uncertainties.
Across all five masses, we consistently find that at higher momenta there is more statistical noise in the extraction of G E (Q 2 ).
In Fig. 6, we take the plateau values from each of the kinematics listed in Table II at m π = 156 MeV and plot their Q 2 dependence. We exclude any kinematics for which we are unable to find a clear plateau, or the variational analysis produces a negative generalised ei- genvalue (as negative eigenvalues indicate issues with the variational analysis, and can cause problems with state identification). We see the contributions from both quark flavours are very similar and each agrees well with a dipole ansatz with G 0 fixed to one, as we are working with single quarks of unit charge. These fits correspond to an RMS charge radius of r 2 1 /2 = √ 12/Λ = 0.684 (19) fm for the doubly represented quark flavour and 0.659(21) fm for the singly represented quark flavour. That these values are smaller than the physical expressions can be ascribed to the finite volume of the lattice [13]. For brevity, we omit similar plots for the other four masses.
In order to compute the form factors of the proton, G E p (Q 2 ), and neutron, G E n (Q 2 ), we need to take the correct linear combinations of the contributions from the doubly and singly represented quark flavours to reintroduce the multiplicity of the doubly represented quark and the physical charges of the up and down quarks.
In Fig. 7 and Table III, we present the electric form factors obtained by these combinations for the lightest pion mass considered here.
In this work, we only consider connected contributions to the nucleon form factors. There is no a priori reason that the disconnected loops could not be included in a PEVA calculation. They were simply omitted from the analysis presented here for computational efficiency. The disconnected loop contributions to the proton and neutron should be approximately the same (exactly the same in our lattice calculations, as we are in the isospin symmetric limit). Hence, if we take the isovector combination G E p (Q 2 ) − G E n (Q 2 ), the disconnected loop con- Table III. GE(Q 2 ) at mπ = 156 MeV for all acceptable kinematics. We present results for the ground-state proton and neutron, as well as isovector combination GE p(Q 2 ) − GE n(Q 2 ) (which is insensitive to disconnected loop corrections).  (17) 0.011 (11) 0.350 (25) tributions will cancel. The form factor values for this combination are also presented in Table III. The form factor for the neutrally charged neutron is close to zero for all masses considered, as expected. Similar to the linear combinations taken for the form factors, we can combine the squared charge radii from the individual quark sectors with the appropriate multiplicities and charge factors to obtain the squared charge radius of the neutron. For all five pion masses, we find a small negative value. For example, at m π = 156 MeV, the neutron's squared charge radius is −0.022 (9) fm. This is qualitatively consistent with the negative squared charge radii observed in experiment. A more quantitative discussion of this effect requires knowledge of the disconnected loop contributions, which are not considered in this work.
The form factor of the proton matches well with a dipole fit with G 0 fixed to one (the charge of the proton). As expected, the charge radii extracted from these dipole fits approach the experimentally measured proton charge radius from below as the pion mass is reduced towards the physical point.
As discussed in Refs. [14,15], the exact physical value of the proton radius has been a puzzle for the last seven years, since precision laser spectroscopy of muonic hydrogen yielded a proton radius of 0.840 87(39) fm [16] in 2010. This value is 4.6 %, or 5.6σ lower than the codata 2014 world average of 0.8751(61) fm [17], from a combination of laser spectroscopy of electronic hydrogen and deuterium, and elastic electron scattering. Recent precision results from new laser spectroscopy of electronic hydrogen provide a proton radius of 0.8335(95) fm [18], which agrees with the muonic hydrogen radius. This suggests that the discrepancy is likely due to systematic errors in the existing results for electronic hydrogen and elastic electron scattering.
Returning to our results, in Fig. 8, we plot the values of dipole fits to the isovector combination on various Q 2 slices as a function of the squared pion mass. We see that the pion-mass dependence is quite smooth, suggesting that the structure of the state is fairly consistent at all five masses considered here. It has a clear trend downwards as the mass is reduced, corresponding to the increasing charge radius. This effect is in accord with the expectations of finite-volume chiral perturbation theory [13].
For all pion masses and kinematics considered in this paper, in the specific case of the electric form factor, there is no conclusive evidence of opposite parity contaminations. Both the PEVA and conventional variational analysis show clear and clean plateaus in G E (Q 2 ) with good excited state control. This supports previous work demonstrating the utility of variational analysis techniques in calculating baryon matrix elements [19,20]. By using such techniques we are able to cleanly isolate precise values for the Sachs electric form factor of the ground-state proton and neutron.

V. SACHS MAGNETIC FORM FACTOR
Moving on to G M (Q 2 ), we once again begin with the lightest pion mass and the lowest momenta. In Fig. 9,  up (Conv.) dp (Conv.) Figure 10. The contributions to the magnetic form factor from single quarks of unit charge for the ground-state nucleon at mπ = 570 MeV for the lowest-momentum kinematics, providing Q 2 = 0.1444(44) GeV 2 . The plateau regions for both analyses are consistent, starting from time slice 23 for all four fits, but the value of the conventional plateau for the singly represented quark (dp) has a magnitude slightly lower than the PEVA plateau.
we see that while the signal is noisier than G E (Q 2 ), the excited-state contaminations present at early Euclidean times are less significant, and for both the PEVA and conventional analyses we are able to find a plateau from time slice 23 to 25. We are cautious in fitting noisy data and restrict fit regimes to avoid large fluctuations. Fig. 10 illustrates a similar plot for m π = 570 MeV. Here the extended plateau is from time slice 23 to 27 and is more representative of the three heaviest pin masses considered. Contrary to the electric case, there is a statistically significant difference in the values of the plateaus from the PEVA and conventional analysis for the singly represented quark flavour. If we take the correlated ratio of G M (Q 2 ) from the conventional analysis to G M (Q 2 ) from the PEVA analysis, we get a value of 0.66 (9). This ratio is clearly less than 1, indicating that the magnitude of the form factor is being significantly underestimated in the conventional analysis. This suggests that despite finding a plateau, the conventional analysis is being affected by opposite-parity contaminations that are introducing a systematic error of approximately 35 %. This is the strongest effect we see across all kinematics considered. However, the conventional plateaus for other kinematics still show a statistically significant deviation from the PEVA plateaus despite having the same fit regions and acceptable χ 2 values. In Fig. 11, we plot the correlated ratio discussed above for the kinematics that give the least noisy extractions of G M (Q 2 ) with acceptable plateaus and positive generalised eigenvalues. We see that for the singly represented quark sector, all such kinematics give a ratio that sits significantly below unity, indicating a statistically significant underestimation of the plateau value by the conventional analysis. While the strength of the effect appears to depend on the specific kinematics, all these ratios indicate a systematic underestimation of the form factor contributions. For the doubly represented quark sector, the effect is smaller and more easily lost in noise, but there are still multiple kinematics where the plateaus are significantly underestimated.
As the states become less relativistic at larger quark masses, we see a reduction in the amount of parity mixing that occurs, and consequentially in the size of the systematic errors, particularly at the heaviest two masses. However, we still observe statistically significant deviations of the ratio below unity. For the heaviest two masses of 570 MeV and 702 MeV, we see a systematic underestimation of the singly represented quark contributions by 5-10 % and at the remaining masses of 411 MeV and 296 MeV, we see a 10-15 % underestimation.
These results provide strong evidence for oppositeparity contaminations in conventional extractions. These contaminations have a clear effect on the extracted magnetic form factor at all five pion masses, on the order of 10 % for the doubly represented quark sector (u p ) and 20 % for the singly represented quark sector (d p ) at the lighter masses. Moving forward, use of the PEVA technique will be critical in precision calculations of G M (Q 2 ) for the ground-state nucleon, for which such systematic up dp Figure 11. Ratios of conventional plateaus to PEVA plateaus for ground state GM (Q 2 ) at mπ = 156 MeV. For clarity, ratios with large statistical errors have been excluded from the plot. If the plateaus were consistent, the points should be distributed about 1.0. For the doubly represented quark flavour (up), some kinematics match this expectation, but others sit more than one standard deviation low, suggesting that the conventional variational analysis may be contaminated by opposite-parity states. The singly represented quark flavour (dp) shows an even larger effect, strongly supporting the presence of opposite-parity contaminations. These results indicate systematic errors of 10-20 % for dp, and perhaps more for some specific kinematics.
errors are unacceptable. We now proceed to examine the extracted form factors. In light of the opposite-parity contaminations present in the conventional extractions, we focus only on the PEVA results. Fig. 12 shows the Q 2 dependence of the contribution to G M (Q 2 ) from each quark flavour at m π = 156 MeV. We see good agreement with a dipole ansatz, with magnetic radii of 0.514(30) fm for the doubly represented quark flavour and 0.85(11) fm for the singly represented quark flavour.
Similar to the electric form factor case described in Section IV, we take linear combinations of the contributions from the doubly and singly represented quark flavours to obtain the magnetic form factors of the proton (G M p (Q 2 )) and neutron G M n (Q 2 ). In addition, we can take the isovector combination (G M p (Q 2 ) − G M n (Q 2 )) to cancel out disconnected loop contributions.
We plot these combinations for the lightest pion mass in Fig. 13, and present the values in Table IV. At all five masses, the magnetic form factors of both the proton and the neutron agree well with a dipole fit. The magnetic radius obtained from each of these fits is close to the electric charge radius of the proton extracted from G E (Q 2 ) at the same pion mass.
In Fig. 14, we plot the quark-mass dependence of the dipole fits to the isovector combination G M p (Q 2 ) − G M n (Q 2 ). We can once again see a quark-mass depend-  Figure 13. GM (Q 2 ) for the ground-state proton and neutron at mπ = 156 MeV. The shaded region corresponds to a dipole fit to the form factor, with a magnetic radius of 0.551(29) fm for the proton and 0.618(31) fm for the neutron. We also include the isovector combination (GM p(Q 2 ) − GM n(Q 2 )), which is insensitive to disconnected loop corrections. ence, with the different Q 2 slices fanning out at lower masses, due to the simultaneous increase in the size of the state and the increasing value of G M (0), as expected in chiral perturbation theory [21,22].
In this section, we demonstrated the importance of the PEVA technique in controlling systematic errors arising from opposite-parity contaminations in extractions of the magnetic form factor for the ground-state nucleon. Due to these opposite-parity contaminations, the conventional analysis underestimates the contribution to the magnetic form factor from the singly represented quark sector by ∼ Table IV. GM (Q 2 ) at mπ = 156 MeV for all acceptable kinematics. We present results for the ground-state proton and neutron, as well as isovector combination GM p(Q 2 ) − GM n(Q 2 ) (which should be free from disconnected loop corrections).  Figure 14. Quark-mass dependence of dipole fits to GM (Q 2 ) for the ground-state isovector combination (GM p(Q 2 ) − GM n(Q 2 )). As in Fig. 8, the marker shapes represent different Q 2 slices and the dashed line corresponds to the physical pion mass. 20 % at light pion masses, whereas the PEVA technique removes the contaminations and gives improved results.

VI. MAGNETIC DIPOLE MOMENT
Returning to the individual quark flavour contributions, and noting that the observed Q 2 dependence of G E (Q 2 ) and G M (Q 2 ) is very similar, we hypothesise that G M (Q 2 ) and G E (Q 2 ) have the same scaling in Q 2 over the range considered here. If this hypothesis is valid, then the ratio of G M (Q 2 ) to G E (Q 2 ) should be independent of Q 2 . Since we are working with an improved conserved vector current, and single quarks of unit charge, G E (0) = 1 exactly, and G M (0) is the contribution of the quark flavour to the magnetic moment (up to scaling by the physical charge). Hence, the ratio is expected to be constant in Q 2 , and equal to the contribution to the magnetic moment from the given quark flavour.
Experimental results show that at high Q 2 , µ G E (Q 2 )/G M (Q 2 ) diverges significantly from unity [23], so our hypothesis must break down at sufficiently high Q 2 . However, over the low Q 2 range we consider here, these experimental results show that µ G E (Q 2 )/G M (Q 2 ) is close to one, and hence within this range G M (Q 2 )/G E (Q 2 ) approximates the magnetic moment.
For all five pion masses, we find that µ Eff (Q 2 ) is indeed approximately constant across the Q 2 range considered. For example, Fig. 15 shows the Q 2 dependence of µ Eff (Q 2 ) at the lightest pion mass. The remaining masses show very similar Q 2 dependence. By taking a constant fit across all kinematics we obtain a robust estimate for the contributions to the magnetic moment of the nucleon from single quarks of unit charge. In the graphs shown here, the statistical errors on this fit are small, and the shaded band showing these errors is almost indistinguishable from the solid line indicating the central value of the fit. Fig. 16 shows the pion mass dependence of these fits. These individual quark-flavour contributions show a smooth pion-mass dependence with an enhancement of the magnetic moments at low pion mass consistent with chiral perturbation theory [21,22,24]. We can take the linear combinations discussed in Section IV to obtain the magnetic moments of the groundstate proton and neutron. The quark mass dependence of these combinations is illustrated in Fig. 17, as well as in Table V. We also present the equivalent magnetic moment extractions from the conventional analysis. At the heavier pion masses, the conventional analysis slightly underestimates the magnetic moment values. At the lighter pion masses, this discrepancy grows rapidly, reaching approximately 10 % at the lightest pion mass considered here.
The magnetic moments of the proton and neutron extracted by PEVA have a similar quark mass dependence to the individual quark-flavour contributions and are close to the experimental values of 2.792 847 350 8(85) µ N for the proton, and −1.913 042 73(45) µ N for the neutron [17]. The small discrepancy between our results and the physical values is due to a combination of disconnected loop contributions which are not included in our calculation, and finite-volume effects. To avoid the disconnected loop corrections, we compare the isovector combination µ p − µ n to the equivalent combination of the experimentally determined magnetic moments. Doing this, we find that we underestimate the experimental value by around 10 %. This remaining discrepancy can be attributed to finite volume corrections.
To address these finite volume corrections, we consider the chiral effective field theory study presented in Ref. [24]. By adjusting our magnetic moment extrac-  Figure 17. Pion-mass dependence of the extracted magnetic moment for the ground-state proton and neutron. To cancel out any disconnected loop contributions, we plot the isovector combination µp −µn. As the physical point is approached, the trend in this combination approaches but doesn't quite reach the physical value of 4.70 µN [17], represented by a black star.
tions at each pion mass by the difference between their chiral contributions for magnetic moments in the infinitevolume and in a 3 fm box, we produce Fig. 18. We see that this correction brings our results very well in line with the experimental value for the isovector magnetic moment.

VII. CONCLUSION
In this paper, we extended the parity-expanded variational analysis (PEVA) technique to the calculation of elastic form factors, and applied it to calculating the Sachs electric and magnetic form factors of the groundstate proton and neutron. This required inspection of the Dirac structure of the three point correlation function and careful selection of appropriate spinor projectors to extract the desired form factors with maximised signal.
Nucleon structure is a vibrant and rich field of study, and there have been investigations of the Sachs electric and magnetic form factors of the ground state nucleon spanning decades. In this paper we focused specifically on the application of the PEVA technique [1] to form factor calculations and on the systematic errors intro- Table V. Magnetic moments at mπ = 156 MeV for all five pion masses. We present results for the ground-state proton and neutron, for both the PEVA analysis and the conventional parity projected analysis. We see that the conventional analysis consistently underestimates the magnetic moments, with the largest effect at smaller pion masses, where it reaches approximately 10 %.   Figure 18. Pion-mass dependence of the extracted magnetic moment for the isovector combination with finite-volume corrections from Ref. [24]. As the physical point is approached, the trend in this combination approaches physical value of 4.70 µN [17], represented by a black star.
duced by untreated opposite-parity contaminations in conventional analyses. We demonstrated the efficacy of variational analysis techniques in general, and PEVA specifically, at controlling excited-state contaminations in the electric form factor. Both the PEVA and conventional variational analysis show clear and clean plateaus, supporting previous work demonstrating the utility of variational analysis in calculating baryon matrix elements [19,20].
In the particular case of the magnetic form factor, we found evidence that the conventional analysis was contaminated by opposite-parity states. For the kinematics considered here, we observe ∼ 20 % underestimation of the magnitude of the contributions to the magnetic form factor from the singly represented quark flavour at the lighter pion masses. This indicates that existing calculations that do not take into account opposite-parity contaminations may be affected by systematic errors on the order of 20 % at physical quark masses. As such, the PEVA technique is critical for precision measurements of the nucleon form factors.
By utilising the PEVA technique and boosted-frame techniques, we are able to successfully extract the Sachs form factors of the ground-state nucleon at a range of Q 2 values. These extractions allow us to investigate the Q 2 dependence of these form factors. By taking ratios of the form factors, we are also able to extract the magnetic moments of both the ground-state proton and the ground-state neutron.
This paper has established the groundwork for applying the PEVA technique to calculating baryonic matrix elements. The applications for future research are broad. The techniques presented here could be applied directly to the examination of other nucleon observables or excited state observables. An extension is underway for the calculation of nucleon transition form factors, and approaches to applying PEVA to spin-3 /2 states. A particularly interesting development in the field that could benefit greatly from the application of the PEVA technique is the computation of non-forward matrix elements from two point correlation functions via the Feynman-Hellmann theorem [25].