Model-independent determination of the nucleon charge radius from lattice QCD
Constantia Alexandrou, Kyriakos Hadjiyiannakou, Giannis Koutsou, Konstantin Ottnad, Marcus Petschlies
AA model-independent determination of the nucleon charge radius from latticeQCD
Constantia Alexandrou,
1, 2
Kyriakos Hadjiyiannakou,
1, 2
GiannisKoutsou, Konstantin Ottnad, ∗ and Marcus Petschlies Department of Physics, University of Cyprus, PO Box 20537, 1678 Nicosia, Cyprus Cyprus Computation-based Science and Technology Research Center, 20 C. Kavafi Street, 2121 Nicosia, Cyprus PRISMA Cluster of Excellence and Institut f¨ur Kernphysik,Johann-Joachim-Becher-Weg 45, University of Mainz, 55099 Mainz, Germany HISKP and BCTP, Rheinische Friedrich-Wilhelms Universit¨at Bonn, 53115 Bonn, Germany
Lattice QCD calculations of nucleon form factors are restricted to discrete values of the Euclideanfour-momentum transfer. Therefore, the extraction of radii typically relies on parametrizing andfitting the lattice QCD data to obtain its slope close to zero momentum transfer. We investigatea new method, which allows to compute the nucleon radius directly from existing lattice QCDdata, without assuming a functional form for the momentum dependence of the underlying formfactor. The method is illustrated for the case of the isovector mean-square charge radius of thenucleon (cid:104) r (cid:105) and the quark-connected contributions to (cid:104) r p (cid:105) and (cid:104) r n (cid:105) for the proton and neutron,respectively. Computations are performed using a single gauge ensemble with N f = 2 + 1 + 1maximally twisted mass clover-improved fermions at physical quark mass and a lattice spacing of a = 0 .
08 fm. ∗ Electronic address: [email protected] a r X i v : . [ h e p - l a t ] J un I. INTRODUCTION
The radius of the proton is a fundamental quantity for atomic, nuclear, and particle physics. In atomicphysics, it enters in the determination of the Rydberg constant, the most precisely known constant in na-ture, as well as in precision tests of quantum electrodynamics. In nuclear physics, it characterizes the sizeof the most abundant hadron in nature and in particle physics it is an input for beyond the standard modelphysics, testing lepton universality and the possible existence of new particles. Electron-proton scatteringwas traditionally used to determine the proton radius that is extracted from the slope of the electric Sachsform factor G E ( q ) extrapolated to zero momentum transfer. The proton radius value extracted from preci-sion measurements of elastic electron scattering cross sections at the Mainz Microtron found a proton radiusof r p = 0 . stat (4) sys fm [1], where the systematic errors are added quadratically. Combining electronscattering data with more accurate data from the hyperfine structure of electronic hydrogen, the value rec-ommended by CODATA was r p = 0 . r p = 0 . r p = 0 . +15 − [4] compatible with older determi-nations; cf. [5]. In a pioneering experiment in 2010, the radius of the proton was extracted from the Lambshift measured for muonic hydrogen [6]. The muonic hydrogen determination is much more accurate and thevalue extracted r p = 0 . r p = 0 . r p = 0 . r p = 0 . stat (12) sys fm [10] in agreement with the muonic hydrogen extraction, toward a possibleresolution of the puzzle. We also refer to Ref. [11] for a summary of the current state of this debate.A theoretical determination of the proton radius starting from the fundamental theory of the strong interactionrequires a nonperturbative framework. In this work, we used the lattice QCD formulation to determine theradius of the proton. The standard procedure in lattice QCD is to compute the electric form factor for finitevalues of the momentum transfer and then perform a fit to determine the slope at zero momentum transfer.However, on a finite lattice, the smallest nonzero momentum is 2 π/L where L is the spatial size of the lattice.Therefore, to reach very small momentum transfers one needs very large lattices. In this work, we use a directmethod to extract the proton radius that does not depend on fitting the form factor. We illustrate the validityof our method by analyzing one ensemble of gauge configuration generated at the physical pion mass.The structure of this paper is as follows: In Sec. II, we describe the general lattice setup, in Sec. III, we explainour new direct approach to extract the radius, in Sec. IV, we discuss our results and in Sec. V, we compare toother lattice QCD determinations and give our conclusions. II. LATTICE SETUP
Our lattice calculations are performed on a single ensemble with N f = 2 + 1 + 1 flavors of twisted mass[12] Wilson clover-improved fermions at maximal twist [13, 14]. The quark masses on this ensemble aretuned to their physical masses; however, the simulations do not account for any isospin splitting for thelight quarks. Furthermore, the ensemble has been tuned to maximal twist to implement automatic O ( a )improvement [15, 16]. The gauge ensemble is simulated with a spatial volume of ( L/a ) = 64 and a temporalextent of T /a = 128. The lattice spacing and pion mass have been determined in Ref. [13] with values of a = 0 . stat (4) sys fm and M π = 138 . stat (6) sys MeV, respectively. The value of M π is very close to itsphysical value in the isospin limit [17]. For further details on the simulation we refer to Ref. [13].Measurements are performed on 750 gauge configurations, measuring on only every fourth configuration whichcorresponds to a separation of four hybrid Monte Carlo trajectories between consequent measurements. Allstatistical errors in our analysis are computed from a binned jackknife procedure taking into account correla-tions and possible residual affects of autocorrelations. A. Electric form factor of the nucleon
The aim of this study is an extraction of the nucleon charge radius; hence we consider the electromagneticmatrix element of the nucleon (cid:104) N ( p f , s f ) | j µ | N ( p i , s i ) (cid:105) = ¯ u ( p f , s f ) (cid:20) γ µ F ( Q ) + σ µν Q ν m N F ( Q ) (cid:21) u ( p i , s i ) , (1)where on the left-hand side N ( p i , s i ) ( N ( p f , s f )) label nucleon states with initial (final) state momentum p i ( p f ) and spin s i ( s f ), and j µ is a vector current insertion which will be discussed below. On the right-handside, u ( p i , s i ), ¯ u ( p f , s f ) denote Dirac spinors, m N the nucleon mass, and we have introduced the Dirac andPauli form factors F ( Q ) and F ( Q ), which depend on the Euclidean four-momentum transfer Q = − q with (cid:126)q = (cid:126)p f − (cid:126)p i . We will work in Euclidean spacetime throughout this study and in the following discussionwe will always replace F ( Q ) and F ( Q ) by the electromagnetic Sachs form factors G E ( Q ) and G M ( Q )which are more convenient for our purposes and related to F ( Q ) and F ( Q ) via G E ( Q ) = F ( Q ) − Q m N F ( Q ) , (2) G M ( Q ) = F ( Q ) + F ( Q ) . (3)The electromagnetic vector current j emµ ( x ) can be expressed in terms of vector currents j fµ ( x ) for the individualquark flavors f through j em µ ( x ) = (cid:88) f e f j fµ ( x ) , (4)where e f denotes the electrical charge for a quark of flavor f . On the lattice, we employ the symmetrizedconserved vector current which for a quark field χ f ( x ) in the twisted mass basis is given by j fµ ( x ) = 14 (cid:2) ¯ χ f ( x + a ˆ µ ) U † µ ( x )(1 + γ µ ) χ f ( x ) − ¯ χ f ( x ) U µ ( x )(1 − γ µ ) χ f ( x + a ˆ µ )+ ¯ χ f ( x ) U † µ ( x − a ˆ µ )(1 + γ µ ) χ f ( x − a ˆ µ ) − ¯ χ f ( x − a ˆ µ ) U µ ( x − a ˆ µ )(1 − γ µ ) χ f ( x ) (cid:3) , (5)where U µ ( x ) denotes a gauge link and ˆ µ a unit vector in the µ direction. We restrict our calculation to lightvalence quarks and the relation of the light quark doublet χ ( x ) = ( χ u ( x ) , χ d ( x )) T to the physical up and downquark at maximal twist is given by the chiral rotation χ ( x ) = 1 √ (cid:0) iγ τ (cid:1) ψ ( x ) , (6)¯ χ ( x ) = ¯ ψ ( x ) 1 √ (cid:0) iγ τ (cid:1) . (7)In the SU (2) isospin symmetry limit, the matrix elements of the electromagnetic current satisfy (cid:104) p | j emµ | p (cid:105) − (cid:104) n | j emµ | n (cid:105) = (cid:104) p | j uµ − j dµ | p (cid:105) ≡ (cid:104) p | j u − dµ | p (cid:105) , (8)where the isovector current j u − dµ ( x ) has been introduced and | p (cid:105) , | n (cid:105) denote proton and neutron states, respec-tively. Unlike the electromagnetic combination, the isovector current does not give rise to quark-disconnecteddiagrams in our lattice simulations. However, in this work, we neglect any quark-disconnected contributionsas our main goal is to present and benchmark a new method to extract nucleon radii. We remark that thesecontributions have been shown to be small compared to the quark-connected contribution in Refs. [18, 19].In the twisted mass lattice regularization, the flavor symmetry is broken from SU (2) to the subgroup U (1) and the isospin symmetry-based relations, Eq. (8), are valid at nonzero lattice spacing up to lattice artifacts.In particular, the quark-disconnected contribution from the insertion of the isovector current j u − d will bediscarded from our calculation, since it vanishes in the continuum limit.The point-split vector current in Eq. (5) is the Noether current associated with the residual flavor symmetrygroup U (1) × U (1) and thus remains conserved. Therefore, no additional multiplicative renormalization ofthe matrix element is required. B. Two- and three-point functions
The lattice QCD evaluation of the nucleon matrix element in Eq. (1) requires the computation of two- andthree-point functions C ( (cid:126)p, t f − t i ) = Γ βα (cid:88) (cid:126)x f e − i(cid:126)p · ( (cid:126)x f − (cid:126)x i ) (cid:104) J N,α ( (cid:126)x f , t f ) ¯ J N,β ( (cid:126)x i , t i ) (cid:105) , (9) C µ (Γ ν , (cid:126)q, t op − t i , t f − t i ) = Γ βαν (cid:88) (cid:126)x f ,(cid:126)x op e − i(cid:126)p f · ( (cid:126)x f − (cid:126)x op ) e i(cid:126)p i · ( (cid:126)x op − (cid:126)x i ) (cid:104) J N,α ( (cid:126)x f , t f ) j u − dµ ( (cid:126)x op , t op ) ¯ J N,β ( (cid:126)x i , t i ) (cid:105) , (10)where t i , t op , and t f label initial (source), operator insertion, and final (sink) time slices, respectively. Thespin projectors Γ ν are given by Γ = (1 + γ ) for ν = 0 and Γ k = Γ iγ γ k for ν = k = 1 , ,
3. By timetranslation invariance, only time differences are of physical relevance and it is convenient to introduce thesource-sink time separation t sep = t f − t i and the shorthand t = t op − t i . Our kinematic setup is chosen suchthat the final state is produced at rest, i.e., (cid:126)p f = 0, (cid:126)p i = − (cid:126)q . Finally, the proton interpolating field is given inthe physical basis by J N,α ( x ) = (cid:15) abc (cid:0) u Ta ( x ) Cγ d b ( x ) (cid:1) u c,α ( x ) , (11)where C denotes the charge conjugation matrix and is transformed into the twisted basis in terms of χ u/d by using the chiral rotation in Eq. (7). Since nucleon structure calculations are known to be hampered bya severe signal-to-noise problem, it is crucial to increase the overlap of the interpolating operator with thedesired nucleon ground state, effectively suppressing excited states and allowing to extract a signal at smallerEuclidean time separations. To this end, we apply Gaussian smearing to the quark fields [20, 21], q a ( (cid:126)x, t ) → ˜ q a ( (cid:126)x, t ) = (cid:88) (cid:126)y (cid:2) + α G H ab ( (cid:126)x, (cid:126)y ; U ( t )) (cid:3) N G q b ( (cid:126)y, t ) , q = u, d , (12)with a choice of α G = 0 . N G = 125 smearing steps corresponding to a smearing radius of r G = 0 .
47 fmin physical units [22]. The hopping matrix H ab ( (cid:126)x, (cid:126)y ; U ( t )) is given by H ab ( (cid:126)x, (cid:126)y ; U ( t )) = (cid:88) i =1 (cid:104) U abi ( x ) δ x,y − ˆ i + U ∗ bai ( x − ˆ i ) δ x,y +ˆ i (cid:105) , (13)where a , b are color indices. Furthermore, we use APE-smeared [23] gauge links U in the construction of H with a smearing parameter of α APE = 0 . N APE = 50 smearing steps.For the computation of three-point functions, we use sequential inversions through the sink [24]. Therefore,we need to perform separate inversions for each choice of the source-sink time separation, the projector indexand—in principle—the momentum at sink. However, the latter is fixed to be zero, as mentioned before. Sincewe are interested in the nucleon charge radius, it is sufficient to consider three-point functions projected withΓ . Three-point functions are computed for five values of t sep as listed in Table I. In order to obtain comparableeffective statistics at each value of t sep , we add additional source positions for increasing values of t sep . Thenumber of source positions N s per source-sink time separation has also been included in Table I together withthe total number of measurements N meas . The source positions are randomly and independently chosen oneach gauge configuration. Two-point functions are computed with matching statistics for each value of t sep t sep /a t sep / fm N s N meas
12 0.97 4 300014 1.13 6 450016 1.29 16 1200018 1.46 48 3600020 1.62 64 48000TABLE I: Overview of the values of t sep used in the computation of three-point functions in lattice and physical units. N s is the number of source positions for each source-sink time separation on each of the 750 gauge configuration and N meas the corresponding, total number of measurements. from the forward propagators obtained in the calculation of three-point functions; hence, there is a possiblechoice of either using matching statistics for two- and three-point functions at each value of t sep or always usingthe full available statistics for the two-point functions. However, we found that fully preserving correlation byusing exactly matching statistics yields a slight advantage with respect to the resulting statistical fluctuations.Calculations are performed using an appropriately tuned multigrid algorithm [25–27] for the efficient inversionof the Dirac operator that is required for the computation of the quark-connected diagrams. C. Ratio method
Extracting the physical matrix elements in Eq. (1) requires the cancellation of unknown overlap factors inthe three-point function, which can be achieved by forming an optimized ratio involving two-point functions[28–30] R µ (Γ ν , (cid:126)q, t, t sep ) = C µ (Γ ν , (cid:126)q, t, t sep ) C ( (cid:126) , t sep ) (cid:115) C ( (cid:126)q, t sep − t ) C ( (cid:126) , t ) C ( (cid:126) , t sep ) C ( (cid:126) , t sep − t ) C ( (cid:126)q, t ) C ( (cid:126)q, t sep ) . (14)At large Euclidean time separations t and t f − t , the ground state contribution is expected to dominateasymptotically and the ratio approaches a plateau,lim t →∞ lim t f − t →∞ R µ (Γ ν , (cid:126)q, t, t sep ) = Π µ (Γ ν , (cid:126)q ) . (15)The electromagnetic Sachs form factors can be extracted from Π µ (Γ ν , (cid:126)q ) for appropriate choices of insertionand projector indices, Π (Γ , (cid:126)q ) = − C E N ( (cid:126)q ) + m N m N G E (cid:0) Q (cid:1) , (16)Π i (Γ , (cid:126)q ) = − C i m N q i G E (cid:0) Q (cid:1) , (17)Π i (Γ k , (cid:126)q ) = − C m N (cid:15) ijk q j G M (cid:0) Q (cid:1) , (18)where C = (cid:113) m N E N ( (cid:126)q )( E N ( (cid:126)q )+ m N ) . For the computation of the nucleon electric charge radius, we use the firstrelation, which gives by far the best signal for G E ( Q ) and in addition allows to obtain a result directly atzero momentum transfer.The determination from the second relation involving G E ( Q ) is impeded due to the momentum factor appear-ing on the right-hand side of the equation and similarly in the last equation for G M ( Q ). A direct method forderivatives of form factors at zero momentum has been put forward in [31], based on the algebraic definition ofthe momentum derivative and quark propagator expansion, which promotes the n-point correlation functionby one point per derivative and is thus very costly.In Ref. [32] we have explored model-independent position space methods to remedy this issue for G M ( Q )without resorting to fits. We remark that one of these methods called momentum elimination in the plateau-region is similar to the approach we will introduce in the next section to allow for a direct, model-independentcomputation of the nucleon electric charge radius, which is otherwise hindered by the discrete nature ofmomenta on the lattice. D. Summation method
In nucleon structure calculations it is notoriously difficult to reach ground state dominance due to an exponen-tial signal-to-noise problem; hence, a careful analysis of the corresponding systematics is required. While wehave lattice data available for five values of t sep as listed in Table I, it is not a priori clear that this is sufficientto control excited state effects. Therefore, we use the so-called summation method [33–35] in addition to thedirect plateau method, which allows for a stronger suppression of excited states t sep − t ex (cid:88) t = t ex R µ (Γ ν , (cid:126)q, t, t sep ) = const + ( t sep − t ex + a ) · Π µ (Γ ν , (cid:126)q ) + O ( e − ∆ t sep ) , (19)where t ex = 2 a for the conserved vector current insertion. The contribution from the next higher state witha mass gap ∆ is now suppressed by an additional factor e − ∆ t sep in contrast to the plateau method for whichthe suppression is only ∼ e − ∆ t . III. POSITION SPACE METHOD FOR (cid:104) r E (cid:105) The mean-square charge radius of the nucleon (cid:104) r E (cid:105) is defined from the electric Sachs form factor in Eq. (2)through expansion around small Q , G E ( Q ) = G E (0) (cid:2) − Q (cid:104) r E (cid:105) / O ( Q ) (cid:3) . (20)Since for the proton and the isovector combination one has G p,u − dE (0) = 1, the mean-square charge radius canbe extracted from the expansion via computing (cid:104) r E (cid:105) = − dG E ( Q ) dQ (cid:12)(cid:12)(cid:12)(cid:12) Q =0 . (21)For the neutron, the leading term in Eq. (20) is absent and the normalization factor G nE (0) = 0 is dropped inthe definition to obtain a finite result; hence, (cid:104) r E (cid:105) n can be computed in the same way from Eq. (21). Moreover,for the proton and the isovector combination, it is possible to define a root mean-square charge radius, i.e., r p,u − dE = (cid:113) (cid:104) r E (cid:105) p,u − d , (22)which is not meaningful for the neutron as its mean-square charge radius (cid:104) r E (cid:105) n is negative.Computing (cid:104) r E (cid:105) requires an evaluation of the derivative with respect to Q which is not directly possible in afinite box as only finite, discrete momenta are accessible. In the literature, several methods have been used tocircumvent this issue, e.g., dipole fits. However, any such method introduces model dependence and potentiallylarge, systematic uncertainties considering the physical values of momenta that are typically available in latticesimulations and the steepness of G E ( Q ) close to zero momentum transfer. For a recent review on commonmethods to parametrize electromagnetic form factors, we refer to Ref. [36].In this study, we aim to evaluate the pertinent derivative in Eq. (21) by an integral method for a suitablydefined form factor in the small momentum region given below in Eqs. (28) and (29). It allows for systematicprobing of the dependence on the available supporting lattice data and has a well-defined infinite volume aswell as continuum limit without any further adjustments of parametrization. -0.04-0.0200.020.040.060.080.1 0 5 10 15 20 25 30 ¯ Π u − d ( y , q m a x ) y/a q max = 1 · (2 π/L ) q max = 2 · (2 π/L ) q max = 3 · (2 π/L ) q max = 4 · (2 π/L ) q max = 5 · (2 π/L ) -0.002-0.00100.0010.0020.0030.004 0 5 10 15 20 25 30 ¯ Π n ( y , q m a x ) y/a q max = 1 · (2 π/L ) q max = 2 · (2 π/L ) q max = 3 · (2 π/L ) q max = 4 · (2 π/L ) q max = 5 · (2 π/L ) FIG. 1: Position space lattice data Π( y ) for different choices of the momentum cutoffs q max and t sep = 0 .
97 fm. Resultsare shown for the isovector combination (left panel) and the neutron (right panel). Data for different values of q max have been displaced horizontally to improve readability. We present the calculation in its simplest one-dimensional form based on (numerical) form factor data foron-axis three-momenta ( (cid:126)q ∝ (cid:126)e k ). It can be extended to incorporate arbitrary momentum directions, which,however, require taking into account anisotropy effects and analytic continuation.In the continuum we have Q = − m N ( E N ( q ) − m N ) given our momentum setup for the nucleon at sourceand sink, and we can straightforwardly replace the Q derivative by ddQ = − E N ( q ) m N ddq (23)for on shell E N = (cid:112) m N + q . In particular, we can use ddQ G E ( Q ) (cid:12)(cid:12)(cid:12)(cid:12) Q =0 = − lim q → G E ( q ) − G E (0) q . (24)Thus, on the lattice, we define the form factorΠ( q ) = − (cid:115) m N E N ( q ) E N ( q ) + m N (cid:88) σ = ± , k =1 , , Π (Γ , (cid:126)q = σ q (cid:126)e k ) , q = 2 π/L × { , , , . . . , q max } . (25)In Eq. (25), q max denotes the maximal value of on-axis momentum, for which numerical data for the ratio R µ and thus the form factor Π can be obtained given the achieved statistical precision for two- and three-point correlation functions. Apart from the dependence on the source-sink time separation t sep , checking thesaturation of the integral defining r E under variation of q max will be a major point of the systematic erroranalysis. To this end, we will also use model data to study the influence of the large- Q tail, i.e. allowing usto implement larger values of q max than the statistical precision of the available lattice data would otherwisepermit.Together with Π( q ) in Eq. (25) we obtain its counterpart in position space ¯Π( n ) from the discrete Fouriertransform for − N/ ≤ n ≤ N/ n = y/a , N = L/a and ¯Π( − y ) = ¯Π( y ). Sample data for the positionspace form factor ¯Π( y, q max ) with varying cutoff q max in the discrete Fourier transform are shown in Fig. 1.Upon inverse Fourier transform, we can express the resulting form factor as a function of the lattice momentumˆ k = 2 sin( k/
2) via Π (cid:16) ˆ k (cid:17) = N/ (cid:88) n =0 c n ¯Π( n ) cos (cid:16) n asin (cid:16) ˆ k/ (cid:17)(cid:17) (26) c n = (cid:40) n = 1 , . . . , N/ − n = 0 , N/ . Using the relation cos (cid:16) n asin(ˆ k/ (cid:17) = ( − n T n (cid:16) ˆ k/ (cid:17) (27)where T n are Chebyshev polynomials of the first kind, this givesΠ (cid:16) ˆ k (cid:17) = N/ (cid:88) n =0 c n ¯Π( n ) ( − n T n (ˆ k/
2) (28)and with respect to Eq. (24) 1ˆ k (cid:16) Π(ˆ k ) − Π(0) (cid:17) = N/ (cid:88) n =1 c n ¯Π( n ) P n (ˆ k ) (29)where P n (ˆ k ) = (cid:104) ( − n T n (ˆ k/ − (cid:105) / ˆ k are polynomials in ˆ k of degree n −
1. The difference quotient inEq. (29) is an analytic function with a well-defined infinite volume and continuum limit and gives in particular (cid:104) r E (cid:105) = lim ˆ k → k (cid:16) Π(ˆ k ) − Π(0) (cid:17) = N/ (cid:88) n =1 c n ¯Π( n ) P n (0) . (30)Since for the Chebyshev polynomials T n ( x ) = ( − n +1 (2 n ) x + O (cid:0) x (cid:1) , we also have the familiar integral /summation form (cid:104) r E (cid:105) ∝ N/ (cid:88) n =1 c n ¯Π( n ) (2 n ) , (31)or equivalently in terms of the original data Π( q ) in Eq. (25), (cid:104) r E (cid:105) ∝ q max (cid:88) q =0 w ( q ) Π( q ) (32) w ( q ) = N/ (cid:88) n =1 c n (2 n ) cos( qn ) . (33)The kernel weight w in Eq. (32) for the lattice setup used in this study and exemplary data for the integrandin Eq. (33) are shown in Fig. 2. Inspection of the integrands in the right-hand panel of Fig. 2 shows that thecontribution from the form factor at on-axis momenta beyond q/ (2 π/L ) = 5 ∼ IV. RESULTS
The lattice data for the electric Sachs form factor that are used as input in our calculation of the nucleon radiiare shown in Fig. 3 as a function of Q . We have included data from the ratio method in Eq. (15) at all fiveavailable source-sink time separations as well as from the summation in Eq. (19). Results are shown for theproton, neutron and the isovector combination. In any case, we observe that the summation method yieldsresults compatible with the ratio method for the largest value of t sep .In Fig. 4, we show our lattice data for the difference − G u − dE ( Q ) − G u − dE (0)) /Q as a function of Q inphysical units for all five values of t sep /a and the summation method. The extrapolation bands are computedby evaluating Eq. (29) for any given value ofˆ k ( Q ) = 4 sin (cid:32) (cid:114) Q m + Q (cid:33) (34) -400.0-300.0-200.0-100.00.0100.0200.0300.0400.0 0 5 10 15 20 25 30 35 w ( q ) q/ (2 π/L )kernel weight, L/a = 64 -400.0-300.0-200.0-100.00.0100.0200.0300.0400.0 0 5 10 15 20 25 30 35 w ( q ) × G E ( q ) q/ (2 π/L ) dipole1 / log FIG. 2: Left: kernel weight w ( q ) of Eq. (33); right: integrand (kernel weight times evaluated fit for G E ( q )) fordipole- (red open squares) and inverse logarithmic-type decay in the large momentum tail region ( q/ (2 π/L ) (cid:38) and multiplying the expression by ˆ k /Q leading toΠ( Q ) − Π(0) Q = − N/ (cid:88) n =1 c n ¯Π( n ) 2 sin (cid:18) n (cid:113) Q m + Q (cid:19) Q , (35)where (cid:113) Q m + Q = q . Note that the factor related to substituting ˆ k → Q is one for the nucleon radius,i.e., at Q = ˆ k = 0. However, at nonzero Q , the change of variable ˆ k → Q has to be taken into accountexplicitly to make contact with the discrete lattice data, which is achieved by the above expression. Here wehave consistently chosen a rather small value of q max = 4 · (2 π/L ) for the extrapolations in all six plots. Still, theresulting bands describe well the lattice data for Q (cid:46) . Q = 1 GeV . This is expected as thevalues of q entering the initial Fourier transform from momentum to position space are restricted to values of Q below 0 . for q max = 4 · (2 π/L ). The behavior is particularly visible for the data at smaller values of t sep which are statistically more precise. A. Study of the large- Q tail contribution In order to further study the systematics associated with the truncation of the Fourier transform for Π( q )we have explored ways of generating synthetic data for the tail of the form factor to avoid this truncationaltogether. To this end, we use several fit models that are applied to the lattice data at low values of Q . In asecond, step we use the jackknife samples for the fit parameters to generate synthetic data samples for valuesof Q corresponding to on-axis momenta (cid:126)q = ( q, , T with q > q max . This synthetic data are then used toperform the Fourier transform of Π( q ) for all q > q max in addition to the actual lattice data at q ≤ q max . Asimple choice for the proton and isovector form factor is the dipole fit model G p,u − dE ( Q ) = A (1 + BQ ) , (36)0 G u − d E ( Q ) Q / GeV t sep = 0 . t sep = 1 . t sep = 1 . t sep = 1 . t sep = 1 . G p E ( Q ) Q / GeV t sep = 0 . t sep = 1 . t sep = 1 . t sep = 1 . t sep = 1 . -0.15-0.1-0.0500.050.10.15 0 0.5 1 1.5 2 G n E ( Q ) Q / GeV t sep = 0 . t sep = 1 . t sep = 1 . t sep = 1 . t sep = 1 . FIG. 3: Lattice data for the ratio in Eq. (14) for G E ( Q ) at all five source-sink time separations for all three isospincombinations of the operator insertion. Only quark-connected contributions are included. − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 0 . t sep = 0 . − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 1 . t sep = 1 . − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 1 . t sep = 1 . − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 1 . t sep = 1 . − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 1 . t sep = 1 . − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresultsummation data not used in extrapolationsummation data used in extrapolation FIG. 4: Lattice data for − G E ( Q ) u − d − G u − dE (0)) /Q and continuous band computed from Eq. (35) using q max =4 · π/L for all five source-sink time separations and the summation method. where A , B are free parameters of the fit. For the neutron, we use a Galster-like parametrization [37, 38]instead of the dipole fit model. This model takes into account that G nE (0) = 0, G nE ( Q ) = CQ (1 + DQ ) (cid:16) Q Λ (cid:17) . (37)The free parameters of this model are denoted by C and D , while Λ = 0 .
71 GeV is treated as a constant.Since the dipole fit model approaches zero quadratically for increasing values of Q , one might anticipate thatusing synthetic data from this model yields a result not very different from the one obtained by a simpletruncation in the Fourier transform of Π( q ) at a reasonable value of q max . Therefore, we consider a 1 / log fitmodel with an unphysically slower decrease at large values of Q , G E ( Q ) = E log (1 + F Q ) (38)2with two fit parameters E and F , which allows us to further investigate the dependence on the choice of themodel. This model can be applied to all three isospin combinations, but it requires a careful choice of thelower bound of the fit range as the model is divergent for Q → G E ( Q ).Corresponding results for the isovector insertion are very similar to the ones for the proton. As expected, thebands from the dipole and Galster-like fits exhibit a much sharper dropoff in the large- Q tail of the formfactor than the corresponding bands from the 1 / log-fit. In general, fit ranges have been chosen such that Q ∈ (cid:2) , . (cid:3) for the dipole and Galster-like fits and Q ∈ (cid:2) . , . (cid:3) for the 1 / log-fit.We found this choice to yield good values of χ / d . o . f . (degrees of freedom) regardless of the source-sink timeseparation as well as for data from the summation method. G p E ( Q ) Q / GeV dipole fit1 /log -fitoriginal data t sep = 1 . -0.04-0.0200.020.040.060.080.1 0 2 4 6 8 10 G n E ( Q ) Q / GeV Galster fit1 /log -fitoriginal data t sep = 0 . FIG. 5: Different fit models of the lattice data for G E ( Q ). Left panel: resulting fit bands for the proton at t sep ≈ .
29 fmfrom a dipole- and 1 / log-fit. Right panel: resulting fit bands for the neutron at t sep ≈ .
97 fm from a Galster-like fitand a 1 / log-fit. As mentioned before, the results for the fit parameters can be used to generate synthetic data for any desiredvalue of Q in the large- Q tail of G E ( Q ). This allows us to compute the position space form factor ¯Π( y, q max )for any value of q max substituting missing lattice data in the Fourier transform by synthetic data from a givenfit model. Typical results from this procedure are shown in Fig. 6. Regardless of the model, the signal becomesincreasingly peaked at y = 0 for larger values of q max while the large- y tail becomes flatted out, which is theexpected behavior. Moreover, including data generated from the 1 / log fit model enhances this effect furtherat large values of q max compared to the dipole or Galster-like fit.Some results for the final extrapolations together with lattice and model data as a function of Q are shownin Fig. 7 for all three isospin combinations and different fit models. We find that there is hardly any differencevisible between the extrapolations obtained from using synthetic data from the dipole or Galster-like fit modelin the left panels and the 1 / log fit model in the right panels of Fig. 7. In the presented examples we haveemployed model data for either q ≥ · (2 π/L ) or q ≥ · (2 π/L ) depending on whether results have beenobtained at finite values of t sep or from the summation method, respectively. At any rate, this confirms theexpectation from the kernel weight function in Fig. 2 that the final result is indeed saturated by the first fewvalues of on-axis momentum q . Therefore, residual systematic effects related to the truncation of the Fouriertransform must be small. B. Further systematics and final values
The smallness of systematic effects due to truncation of the Fourier transform is further corroborated by adirect comparison of the results for r u − d,pE and (cid:104) r E (cid:105) n at different cutoffs q max as well as using synthetic datafrom the fit models. An overview of results for all available values of t sep together with the summation methodis shown in Fig. 8. The corresponding numerical values are listed in Table II together with the mean-squaredradii for the proton and isovector combination. In general, the results using data from the two fit models arein excellent agreement at any given value of t sep . A significant difference is only observed between the two3 -0.0500.050.10.150.2 0 5 10 15 20 25 30 ¯ Π u − d ( y , q m a x ) y/a original data only, q max = 1 · (2 π/L )original data only, q max = 2 · (2 π/L )original data only, q max = 4 · (2 π/L )dipole fit data, q max = 8 · (2 π/L )dipole fit data, q max = 16 · (2 π/L )dipole fit data, q max = 32 · (2 π/L ) -0.0500.050.10.150.2 0 5 10 15 20 25 30 ¯ Π u − d ( y , q m a x ) y/a original data only, q max = 1 · (2 π/L )original data only, q max = 2 · (2 π/L )original data only, q max = 4 · (2 π/L )1 /log -fit data, q max = 8 · (2 π/L )1 /log -fit data, q max = 16 · (2 π/L )1 /log -fit data, q max = 32 · (2 π/L ) -0.00200.0020.0040.0060.008 0 5 10 15 20 25 30 ¯ Π n ( y , q m a x ) y/a original data only, q max = 1 · (2 π/L )original data only, q max = 2 · (2 π/L )original data only, q max = 4 · (2 π/L )Galster fit data, q max = 8 · (2 π/L )Galster fit data, q max = 16 · (2 π/L )Galster fit data, q max = 32 · (2 π/L ) -0.00200.0020.0040.0060.008 0 5 10 15 20 25 30 ¯ Π n ( y , q m a x ) y/a original data only, q max = 1 · (2 π/L )original data only, q max = 2 · (2 π/L )original data only, q max = 4 · (2 π/L )1 /log -fit data, q max = 8 · (2 π/L )1 /log -fit data, q max = 16 · (2 π/L )1 /log -fit data, q max = 32 · (2 π/L ) FIG. 6: Position space data ¯Π( y, q max ) for the isovector combination (upper row) and the neutron (lower row) fordifferent choices of the momentum cutoff q max and t sep = 0 .
97 fm. Results shown for q max = 8 · (2 π/L ) (pink), q max = 16 · (2 π/L ) (light blue) and q max = 32 · (2 π/L ) (black) use synthetic data generated from the fitted parametersof the models in Eqs. (36)–(38) in addition to the lattice data for q ≤ · (2 π/L ). Data for different values of q max havebeen displaced horizontally to improve readability. t sep / fm (cid:104) r E (cid:105) u − d (cid:104) r E (cid:105) p (cid:104) r E (cid:105) n r u − dE r pE (cid:104) r E (cid:105) in physical units for the isovector combination as well as the proton and the neutronmeasured only from lattice data for all five available values of t sep and the summation method. For the isovectorcombination and the proton results for the root-mean-square radius are included as well. All results have been obtainedfor q max = 5 · (2 π/L ) using only actual lattice data. Errors are statistical only. smallest values of the cutoff q max = 4 · (2 π/L ) and q max = 5 · (2 π/L ) for the statistically most precise data at thesmallest three source-sink time separations. However, the results obtained with a cutoff of q max = 5 · (2 π/L )in the Fourier transform are compatible with the results obtained from modeling the large- Q tail of the formfactor.Regarding excited state contamination, there is a weak trend visible for the first few values of t sep . We findthat for the proton and isovector radius the results at the largest available value of t sep = 1 .
62 fm are in good4 − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 1 . t sep = 1 . − ( G u − d E ( Q ) − G u − d E ( )) / Q / f m Q / GeV extrapolationresult t sep = 1 . t sep = 1 . − ( G p E ( Q ) − G p E ( )) / Q / f m Q / GeV extrapolationresultsummation data not used in extrapolationsummation data used in extrapolationsynthetic data dipole fit − ( G p E ( Q ) − G p E ( )) / Q / f m Q / GeV extrapolationresultsummation data not used in extrapolationsummation data used in extrapolationsynthetic data 1/log fit -0.15-0.1-0.0500.05 0 2 4 6 8 10 − ( G n E ( Q ) − G n E ( )) / Q / f m Q / GeV extrapolationresult t sep = 0 . t sep = 0 . -0.15-0.1-0.0500.05 0 2 4 6 8 10 − ( G n E ( Q ) − G n E ( )) / Q / f m Q / GeV extrapolationresult t sep = 0 . t sep = 0 . FIG. 7: Lattice (filled symbols) and synthetic data (open symbols) at large Q for − G E ( Q ) u − d − G u − dE (0)) /Q together with the continuous extrapolation band computed from Eq. (35). The first, second, and last rows showresults for the isovector data at t sep = 1 .
29 fm, the proton data from the summation method and the neutron dataat t sep = 0 .
97 fm, respectively. In the left column the dipole fit and Galster-like fit models have been used in thegeneration of synthetic data, while for the right column synthetic data from a 1 / log-fit have been used. For theisovector combination and the neutron, synthetic data have been used for q ≥ · (2 π/L ) in the construction of theextrapolation, while for the more noisy proton data from the summation method the switch to model data occurs at q ≥ · (2 π/L ). agreement with the summation method. Therefore, we quote the value obtained at the largest available valueof t sep and for q max = 5 · (2 π/L ) as our final results for the isovector and proton radius, i.e., r u − dE = 0 . stat (05) a fm , (39) r pE = 0 . stat (04) a fm , (40)where the first error is statistical and the second error refers to the scale setting uncertainty. Note that all finalresults are obtained from using lattice data only and that synthetic data is only used in the cross-checkingthe impact of the neglected tail contribution. At the current level of precision, the statistical error clearly5 r u − d E / f m t sep / fm q max = 4 · π/Lq max = 5 · π/Lq max = 6 · π/L final result q max = π , log fit q max = π , dipole fit 00.20.40.60.81 summationmethod0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 r p E / f m t sep / fm q max = 4 · π/Lq max = 5 · π/Lq max = 6 · π/L final result q max = π , log fit q max = π , dipole fit -0.3-0.25-0.2-0.15-0.1-0.0500.05 summationmethod0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 h r E i n / f m t sep / fm q max = 4 · π/Lq max = 5 · π/Lq max = 6 · π/L final result q max = π , log fit q max = π , Galster fit FIG. 8: Overview of results for r u − dE , r pE and (cid:104) r E (cid:105) n in physical units. Results are shown as a function of t sep /a and forthe summation method using different values of q max as well as for two possible choices of modeling the large- Q tailof the lattice data for each isospin combination. Results shown for q max / (2 π/L ) = 4 , , q max = π use model data for q ≥ · (2 π/L ). dominates and systematics due to the scale settings are negligible. The data for the neutron combination aremore noisy and the signal for (cid:104) r E (cid:105) n is essentially lost after the third source-sink time separation. For thatreason, we quote the value obtained at t sep = 1 .
29 fm and q max = 5 · (2 π/L ) as the final result (cid:104) r E (cid:105) n = − . stat (01) a fm . (41) C. Comparison with other studies
Although electromagnetic form factors have been studied within lattice QCD since many years, it is onlyrecently that they have been extracted using simulations with physical values of the light quark masses, referredto as physical point ensembles. Also, while there are a number of studies for the isovector combination, namely,the difference between proton and neutron form factors, for the proton and neutron themselves the resultsare scarce. This is due to the complexity of evaluating accurately contributions from disconnected quarkloops. Therefore, we compare the results of our direct method with recent results extracted using simulationswith nearly physical pion masses neglecting the quark-disconnected contributions. Such contributions yield acorrection at the 15% level for the neutron electric form factors, while for the proton the correction is only atthe one percent level. We summarize below recent simulations used by various collaborations for the extractionof the nucleon electromagnetic form factors.(i) ETMC analyzed three physical point ensembles [18, 39] of twisted mass fermions: a N f = 2 + 1 + 1ensemble with m π L = 3 .
62 which has also been used in the present study, and two N f = 2 ensembleswith m π L = 2 .
98 and m π L = 3 .
97 and the same lattice spacing of a = 0 . O ( a )-improved; hence, cutoff effects are of O ( a ).6(ii) LHPC analyzed two N f = 2 + 1 ensembles simulated with two step HEX-smeared clover fermions: anensemble with m π = 149 MeV, lattice spacing a = 0 .
116 fm, and m π L = 4 .
21 [40] and one ensemble at m π = 135 MeV and m π L = 4 with a finer lattice spacing of a = 0 .
093 fm [41]. In the latter study theyemployed a momentum derivative method to extract directly the radii. They use improved currents;therefore, their cutoff effects are O ( a ).(iii) The PACS Collaboration analyzed two ensembles of N f = 2 + 1 stout-smeared clover fermions: oneensemble of m π = 146 MeV, a = 0 .
085 fm, and m π L = 6 [42], and one ensemble with m π = 135 MeV, a = 0 . m π L = 7 . O ( a ).(iv) The PNDME Collaboration [44] analyzed eleven ensembles of N f = 2 + 1 + 1 highly improved staggeredquarks simulated by MILC Collaboration. They used a mixed action setup with clover fermions in thevalence sector. The ensembles have lattice spacings a (cid:39) . , . , . , .
15 fm and the pion massesare m π (cid:39) , ,
315 MeV. Using a combined fit analysis, they performed a chiral, continuum andinfinite volume extrapolation. We limit ourselves here to their results obtained using the m π = 135 MeVensemble with a = 0 . m π L = 3 .
7. No improved currents have been used; thus, cutoffeffects are of O ( a ). q › r fi u d [fm] H C O D A T A This workETMC(a) '19 ETMC(b) '19ETMC '17 LHPC '14LHPC '17 PACS '18 PACS '19 PNDME '19 q › r fi p [fm] H C O D A T A › r fi n [fm ] P D G FIG. 9: Results from this study (red circles) for (cid:112) (cid:104) r E (cid:105) u − d (left panel), (cid:112) (cid:104) r E (cid:105) p (middle panel), and (cid:104) r E (cid:105) n (right panel)compared to results from other lattice QCD analyses. Results are shown for ETMC for the N f = 2 + 1 + 1 ensemble(green left-pointing triangles), N f = 2 and m π L (cid:39) N f = 2 and m π L (cid:39) N f = 2 + 1 ensemble with m π L = 4 .
21 [40] (brown downward-pointing triangles) and an N f = 2 + 1 ensemble with m π L = 4 [41] (square magenta); PACS using an N f = 2 + 1ensemble with m π L = 6 [42] (cyan rhombus) and an ensemble with m π L = 7 . N f = 2 + 1 + 1 mixed action and their physical point ensemble with m π L = 3 . In Fig. 9 we compare the results for the radii as extracted from this study with the aforementioned studies.For the isovector combination we find that the value computed within our direct method is compatible with7all other studies. We note that one of the three ETMC values is extracted by analyzing the same ensemblebut using a dipole fit to extract the radius instead. For the proton, we find very good agreement with mostother determinations. Only the study by PACS using a m π L = 7 . µH experiment and that from CODATA. This holdsalso true for other recent lattice results, e.g., Ref. [45] that has not been included in Fig. 9. Therefore, it isnot yet possible to make a distinction between these two values based on lattice results, even if disconnectedcontributions are included. In fact, disconnected contributions have been included in Ref. [18] and shownto yield a correction of at most 1% for the proton radius which is significantly below the current statisticalprecision. For the neutron, the relative errors on lattice data are much larger and there is overall agreementamong lattice results. The value extracted in this work is in fact very close to the PDG value, despite thatquark-disconnected contributions are not yet included. V. SUMMARY AND OUTLOOK
In this work we extract the rms charge radii avoiding a fit
Ansatz to the electric form factors of the nucleonthat may introduce a model error. The method has been shown to be insensitive toward the large- Q tailof the form factor by testing different Ans¨atze to model the large Q dependence of the form factors wherelattice QCD data are not available. In particular, even using an unphysical Ansatz that falls very slowly with Q the results are unaffected, demonstrating that lattice QCD results up to about Q ∼ determinethe radii for the currently available lattice volume and statistical precision. Comparing the results for r p − nE , r pE and r nE extracted by applying this approach to the values when the traditional approach of fitting theelectric form factors is used for the same ensemble we find consistent values. Moreover, our values agreewith the model-independent determination using the approach of algebraic derivative and quark propagatorexpansion for the direct calculation of form factor derivatives [31] employed by LHPC [41]. In all approaches,the statistical errors on the proton charge radius are still large and lattice results can currently not distinguishbetween the values extracted from muonic hydrogen and older electron scattering experiments. Besides, mostcalculations do not yet contain quark-disconnected contribution.We plan to implement this approach to study the magnetic radii and moments where an increased precision inthe lattice QCD data will be required. For the neutron, one has also to include the disconnected contributionwhich may bring the value closer to the experimental one. Similarly, disconnected diagrams will have to beincluded if aiming for (cid:46)
2% precision on the proton radius. Finally, nonzero lattice spacing and finite volumeartifacts need to be evaluated. This can only take place when ensembles at physical quark mass values aresimulated for smaller lattice spacings and larger volumes. Such a program will be possible as these ensemblesare simulated and analyzed.
ACKNOWLEDGMENTS [1] J. C. Bernauer et al. ( A1 Collaboration), Phys. Rev. Lett. , 242001 (2010), arXiv:1007.5076 [nucl-ex] .[2] P. J. Mohr, B. N. Taylor and D. B. Newell, Rev. Mod. Phys. , 1527 (2012), arXiv:1203.5425[physics.atom-ph] .[3] I. T. Lorenz, H. W. Hammer and U.-G. Meißner, Eur. Phys. J. A48 , 151 (2012), arXiv:1205.6628 [hep-ph] .[4] I. T. Lorenz, U.-G. Meißner, H. W. Hammer and Y. B. Dong, Phys. Rev.
D91 , 014023 (2015), arXiv:1411.1704[hep-ph] .[5] P. Mergell, U. G. Meißner and D. Drechsel, Nucl. Phys.
A596 , 367 (1996), arXiv:hep-ph/9506375 [hep-ph] .[6] R. Pohl et al. , Nature (London) , 213 (2010).[7] A. Antognini et al. , Science , 417 (2013).[8] A. Beyer et al. , Science , 79 (2017).[9] H. Fleurbaey, S. Galtier, S. Thomas, M. Bonnaud, L. Julien, F. Biraben, F. Nez, M. Abgrall, and J. Gu´ena, Phys.Rev. Lett. , 183001 (2018), arXiv:1801.08816 [physics.atom-ph] .[10] W. Xiong et al. , Nature (London) , 147 (2019).[11] H.-W. Hammer and U.-G. Meißner, Sci. Bull. , 257 (2020), arXiv:1912.03881 [hep-ph] .[12] R. Frezzotti, P. A. Grassi, S. Sint and P. Weisz ( Alpha
Collaboration), J. High Energy Phys. , 058 (2001), arXiv:hep-lat/0101001 [hep-lat] .[13] C. Alexandrou et al. , Phys. Rev. D98 , 054518 (2018), arXiv:1807.00495 [hep-lat] .[14] A. Abdel-Rehim et al. ( ETM
Collaboration), Phys. Rev.
D95 , 094515 (2017), arXiv:1507.05068 [hep-lat] .[15] R. Frezzotti and G. C. Rossi, J. High Energy Phys. , 007 (2004), arXiv:hep-lat/0306014 [hep-lat] .[16] R. Frezzotti and G. C. Rossi, Nucl. Phys. Proc. Suppl. , 193 (2004), arXiv:hep-lat/0311008 [hep-lat] .[17] S. Aoki et al. , Eur. Phys. J. C77 , 112 (2017), arXiv:1607.00299 [hep-lat] .[18] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, and A.V. Aviles-Casco, Phys. Rev.
D100 , 014509 (2019), arXiv:1812.10311 [hep-lat] .[19] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, and G. Koutsou, arXiv:1909.10744 [hep-lat] .[20] S. Gusken, Nucl. Phys. Proc. Suppl. , 361 (1990).[21] C. Alexandrou, S. Gusken, F. Jegerlehner, K. Schilling and R. Sommer, Nucl. Phys. B414 , 815 (1994), arXiv:hep-lat/9211042 [hep-lat] .[22] C. Alexandrou, et al. , Phys. Rev.
D101 , no.3, 034519 (2020) arXiv:1908.10706 [hep-lat] .[23] M. Albanese et al. ( APE
Collaboration), Phys. Lett.
B192 , 163 (1987).[24] D. Dolgov et al. ( LHPC, TXL
Collaboration), Phys. Rev.
D66 , 034506 (2002), arXiv:hep-lat/0201021[hep-lat] .[25] S. Bacchio, C. Alexandrou, J. Finkenrath, A. Frommer, K. Kahl, and M. Rottmann, Proc. Sci.
LATTICE2016 ,259 (2016), arXiv:1611.01034 [hep-lat] .[26] S. Bacchio, C. Alexandrou and J. Finkerath, EPJ Web Conf. , 02002 (2018), arXiv:1710.06198 [hep-lat] .[27] C. Alexandrou, S. Bacchio, J. Finkenrath, A. Frommer, K. Kahl, and M. Rottmann Phys. Rev.
D94 , 114509(2016), arXiv:1610.02370 [hep-lat] .[28] C. Alexandrou, G. Koutsou, J. W. Negele and A. Tsapalis, Phys. Rev.
D74 , 034508 (2006), arXiv:hep-lat/0605017 [hep-lat] .[29] C. Alexandrou, M. Brinet, J. Carbonell, M. Constantinou, P. A. Harraud, P. Guichon, K. Jansen, T. Korzec, andM. Papinutto, Phys. Rev.
D83 , 094502 (2011), arXiv:1102.2208 [hep-lat] .[30] C. Alexandrou, M. Constantinou, S. Dinter, V. Drach, K. Jansen, C. Kallidonis, and G. Koutsou, Phys. Rev.
D88 , 014509 (2013), arXiv:1303.5979 [hep-lat] .[31] G. M. de Divitiis, R. Petronzio and N. Tantalo, Phys. Lett.
B718 , 589 (2012), arXiv:1208.5914 [hep-lat] .[32] C. Alexandrou, M. Constantinou, G. Koutsou, K. Ottnad and M. Petschlies (
ETM
Collaboration), Phys. Rev.
D94 , 074508 (2016), arXiv:1605.07327 [hep-lat] .[33] L. Maiani, G. Martinelli, M. L. Paciello and B. Taglienti, Nucl. Phys.
B293 , 420 (1987).[34] S. J. Dong, K. F. Liu and A. G. Williams, Phys. Rev.
D58 , 074504 (1998), arXiv:hep-ph/9712483 [hep-ph] .[35] S. Capitani, M. Della Morte, G. von Hippel, B. J¨ager, A. J¨uttner, B. Knippschild, H. B. Meyer, and H. Wittig,Phys. Rev.
D86 , 074502 (2012), arXiv:1205.0180 [hep-lat] .[36] V. Punjabi, C. F. Perdrisat, M. K. Jones, E. J. Brash and C. E. Carlson, Eur. Phys. J.
A51 , 79 (2015), arXiv:1503.01452 [nucl-ex] .[37] S. Galster, H. Klein, J. Moritz, K. H. Schmidt, D. Wegener, and J. Bleckwen Nucl. Phys.
B32 , 221 (1971). [38] J. J. Kelly, Phys. Rev. C70 , 068202 (2004).[39] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou, K. Jansen, C. Kallidonis, G. Koutsou, and A. V. Aviles-Cascio, Phys. Rev.
D96 , 034503 (2017), arXiv:1706.00469 [hep-lat] .[40] J. R. Green, J. W. Negele, A. V. Pochinsky, S. N. Syritsyn, M. Engelhardt, and S. Krieg, Phys. Rev.
D90 , 074507(2014), arXiv:1404.4029 [hep-lat] .[41] N. Hasan, J. Green, S. Meinel, M. Engelhardt, S. Krieg, J. Negele, A. Pochinsky, and S. Syritsyn, Phys. Rev.
D97 , 034504 (2018), arXiv:1711.11385 [hep-lat] .[42] K.-I. Ishikawa, Y. Kuramashi, S. Sasaki, N. Tsukamoto, A. Ukawa, and T. Yamazaki(
PACS
Collaboration), Phys.Rev.
D98 , 074510 (2018), arXiv:1807.03974 [hep-lat] .[43] E. Shintani, K.-I. Ishikawa, Y. Kuramashi, S. Sasaki and T. Yamazaki, Phys. Rev.
D99 , 014510 (2019), arXiv:1811.07292 [hep-lat] .[44] Y.-C. Jang, R. Gupta, H.-W. Lin, B. Yoon and T. Bhattacharya, Phys. Rev.
D101 , 014507 (2020), arXiv:1906.07217 [hep-lat] .[45] S. Capitani, M. Della Morte, D. Djukanovic, G. von Hippel, J. Hua, B. J¨ager, B. Knippschild, H. B. Meyer, T. D.Rae, and H. Wittig, Phys. Rev.
D92 , 054511 (2015), arXiv:1504.04628 [hep-lat] .[46] P. J. Mohr, D. B. Newell and B. N. Taylor, Rev. Mod. Phys. , 035009 (2016), arXiv:1507.07956[physics.atom-ph] .[47] M. Tanabashi et al. ( Particle Data Group
Collaboration), Phys. Rev.
D98 , 030001 (2018).[48] K. Jansen and C. Urbach, Comput.Phys.Commun. , 2717 (2009), arXiv:0905.3331 [hep-lat] .[49] A. Deuzeman, S. Reker, and C. Urbach, Comput.Phys.Commun. , 1321 (2012), arXiv:1106.4177 [hep-lat]arXiv:1106.4177 [hep-lat]