Proton and neutron electromagnetic form factors from lattice QCD

TL;DR

This paper presents lattice QCD calculations of proton and neutron electromagnetic form factors at physical quark masses, including disconnected contributions, to improve understanding of nucleon structure.

Contribution

It provides the first high-precision computation of disconnected contributions to nucleon form factors at the physical point using twisted mass fermions.

Findings

Disconnected contributions are non-negligible, up to 15% for neutron electric charge radius.

Form factors, radii, and magnetic moments are determined with multiple sink-source separations.

Results enhance understanding of nucleon electromagnetic structure from first principles.

Abstract

The electromagnetic form factors of the proton and the neutron are computed within lattice QCD using simulations with quarks masses fixed to their physical values. Both connected and disconnected contributions are computed. We analyze two new ensembles of $N_f = 2$ and $N_f = 2 + 1 + 1$ twisted mass clover-improved fermions and determine the proton and neutron form factors, the electric and magnetic radii, and the magnetic moments. We use several values of the sink-source time separation in the range of 1.0 fm to 1.6 fm to ensure ground state identification. Disconnected contributions are calculated to an unprecedented accuracy at the physical point. Although they constitute a small correction, they are non-negligible and contribute up to 15% for the case of the neutron electric charge radius.

Tables

Table 1. Table 1: Simulation parameters for the N f = 2 + 1 + 1 subscript 𝑁 𝑓 2 1 1 N_{f}{=}2{+}1{+}1   Alexandrou  et al. ( 2018a ) and N f = 2 subscript 𝑁 𝑓 2 N_{f}{=}2   Abdel-Rehim  et al. ( 2017 ) ensembles used in this work. When two errors are given, the first error is statistical and the second is systematic. The systematic in the lattice spacing of the N f = 2 subscript 𝑁 𝑓 2 N_{f}=2 ensembles is obtained as described in the text, while the systematic in the pion mass provided in physical units is propagated from the lattice spacing.

ensemble	$c_{\mathrm{SW}}$		$\beta$		$N_f$	$a$ [fm]	Vol.		$a{m}_{\pi }$	$m_{\pi }L$	$a{m}_N$	$m_N/m_{\pi }$	$m_{\pi }$ [GeV]	$L$ [fm]
cB211.072.64	1.	69	1.	778	2+1+1	0.0801(4)	${64}^3$$\times$	$128$	0.05658(6)	3.62	0.3813(19)	6.74(3)	0.1393(7)	5.12(3)
cA2.09.64	1.	57551	2.	1	2	0.0938(3)(1)	${48}^3$$\times$	$96$	0.06208(2)	2.98	0.4436(11)	7.15(2)	0.1306(4)(2)	4.50(1)
cA2.09.48	1.	57551	2.	1	2	0.0938(3)(1)	${64}^3$$\times$	$128$	0.06193(7)	3.97	0.4421(25)	7.14(4)	0.1303(4)(2)	6.00(2)

Table 2. Table 2: Statistics for the evaluation of the connected three-point functions for the N f = 2 + 1 + 1 subscript 𝑁 𝑓 2 1 1 N_{f}{=}2{+}1{+}1 cB211.072.64 ensemble. Columns from left to right are the sink-source time separations, the number of configurations analyzed, the number of source positions per configuration chosen randomly and the total number of measurements for each time separation.

$t_s/a$	$N_{\mathrm{cnfs}}$	$N_{\mathrm{srcs}}$	$N_{\mathrm{meas}}$
12	750	4	3000
14	750	6	4500
16	750	16	12000
18	750	48	36000
20	750	64	48000

Table 3. Table 3: Details on the set-up for the evaluation of the light disconnected diagrams. N cnfs subscript 𝑁 cnfs N_{\rm cnfs} is the number of configurations analyzed, N def subscript 𝑁 def N_{\rm def} is the number of low modes we deflate, N r subscript 𝑁 𝑟 N_{r} the number of noise vectors, and N Had subscript 𝑁 Had N_{\rm Had} the number of Hadamard vectors. N sc subscript 𝑁 sc N_{\rm sc} corresponds to spin-color dilution and N inv subscript 𝑁 inv N_{\rm inv} /conf is the total number of inversions per configuration. N srcs subscript 𝑁 srcs N_{\rm srcs} is the number of randomly distributed smeared point sources per configuration used to obtain the nucleon two-point functions and N meas subscript 𝑁 meas N_{\rm meas} the total number of measurements.

	Loops					Two-point
ensemble	$N_{\mathrm{cnfs}}$	$N_{\mathrm{def}}$	$N_r$	$N_{\mathrm{Had}}$	$N_{\mathrm{sc}}$	$N_{\mathrm{inv}}$/conf	$N_{\mathrm{srcs}}$	$N_{\mathrm{meas}}$
cB211.072.64	750	200	1	512	12	6144	200	150000
cA2.09.48	2120	-	2250	-	-	2250	100	212000

Table 4. Table 4: Statistics for the evaluation of the connected three-point functions for the cA2.09.64 and cA2.09.48 ensembles. For the latter, for t s / a = 16 , 18 subscript 𝑡 𝑠 𝑎 16 18 t_{s}/a=16,18 the connected three-point functions have been computed only for the unpolarized projector. The notation is as in Table  2 .

$t_s/a$	$N_{\mathrm{cnfs}}$	$N_{\mathrm{srcs}}$	$N_{\mathrm{meas}}$
cA2.09.64: $N_f={\mathrm{2\hspace{0.33em}\hspace{0.33em}64}}^3\times 128$ ensemble
12	333	16	5328
14	515	16	8240
16	1040	16	16640
cA2.09.48: $N_f={\mathrm{2\hspace{0.33em}\hspace{0.33em}48}}^3\times 96$ ensemble
10,12,14	578	16	9248
16	530	88	46640
18	725	88	63800

Table 5. Table 5: Our results for the electromagnetic radii and the magnetic moment using the N f = 2 + 1 + 1 subscript 𝑁 𝑓 2 1 1 N_{f}{=}2{+}1{+}1 ensemble for the isovector combination ( p − n ) 𝑝 𝑛 (p{-}n) , isoscalar ( p + n ) 𝑝 𝑛 (p{+}n) , the proton and neutron. The first error is statistical, the second is a systematic due to the fit Ansätz, and the third a systematic due to excited states, derived as explained in the text.

	$\sqrt{⟨r_E^2⟩}$ [fm]	$\sqrt{⟨r_M^2⟩}$ [fm]	$\mu$
$p-n$	0.796(19)(12)(12)	0.712(27)(87)(5)	3.97(15)(2)(5)
$p+n$	0.691(9)(7)(14)	0.695(36)(80)(13)	0.89(4)(3)(3)
$p$	0.742(13)(9)(14)	0.710(26)(80)(6)	2.43(9)(1)(3)
$n$	$⟨r_E^2⟩$ [fm2]	0.716(29)(44)(24)	-1.54(6)(2)(3)
	-0.074(16)(16)(8)

Table 6. Table 6: Results for the electromagnetic form factors using the cB211.072.64 ensemble for the isovector combination G E , M p − G E , M n superscript subscript 𝐺 𝐸 𝑀 𝑝 superscript subscript 𝐺 𝐸 𝑀 𝑛 G_{E,M}^{p}-G_{E,M}^{n} , for the proton G E , M p superscript subscript 𝐺 𝐸 𝑀 𝑝 G_{E,M}^{p} and neutron G E , M n superscript subscript 𝐺 𝐸 𝑀 𝑛 G_{E,M}^{n} , including the disconnected contributions for the latter two form factors.

$Q^2\left[Ge{V}^2\right]$	$G_E^{p-n}\left(Q^2\right)$	$G_E^p\left(Q^2\right)$	$G_E^n\left(Q^2\right)$	$G_M^{u-d}\left(Q^2\right)$	$G_M^p\left(Q^2\right)$	$G_M^n\left(Q^2\right)$
0.000	0.997(3)	0.998(2)	0.001(1)	NA	NA	NA
0.057	0.858(10)	0.874(6)	0.016(5)	3.516(101)	2.156(62)	-1.361(43)
0.113	0.752(11)	0.775(8)	0.023(5)	3.105(78)	1.903(48)	-1.202(33)
0.167	0.662(14)	0.694(9)	0.032(7)	2.801(80)	1.719(47)	-1.082(35)
0.219	0.601(17)	0.631(10)	0.030(8)	2.583(82)	1.580(50)	-1.003(35)
0.270	0.534(14)	0.575(9)	0.040(7)	2.430(62)	1.485(38)	-0.945(26)
0.320	0.482(16)	0.529(11)	0.046(7)	2.224(70)	1.367(43)	-0.857(29)
0.417	0.405(23)	0.450(16)	0.045(10)	1.943(78)	1.200(48)	-0.743(32)
0.464	0.374(21)	0.420(15)	0.047(9)	1.789(71)	1.104(44)	-0.684(29)
0.510	0.363(25)	0.404(18)	0.041(11)	1.655(69)	1.012(42)	-0.644(31)
0.554	0.345(27)	0.385(21)	0.040(11)	1.610(84)	0.994(50)	-0.616(37)
0.598	0.310(43)	0.351(33)	0.041(16)	1.472(127)	0.923(79)	-0.549(51)
0.642	0.291(33)	0.336(27)	0.045(13)	1.495(109)	0.926(67)	-0.570(44)
0.684	0.263(34)	0.308(27)	0.046(13)	1.386(112)	0.866(71)	-0.521(44)
0.767	0.239(60)	0.283(48)	0.043(27)	0.893(149)	0.532(90)	-0.361(65)
0.807	0.250(39)	0.278(31)	0.028(15)	0.942(101)	0.583(62)	-0.361(43)
0.847	0.213(43)	0.249(35)	0.035(16)	1.006(117)	0.605(70)	-0.402(51)
0.886	0.158(61)	0.203(51)	0.048(22)	1.028(211)	0.662(135)	-0.367(80)
0.925	0.150(58)	0.205(50)	0.056(24)	0.891(190)	0.539(115)	-0.353(79)
0.963	0.172(48)	0.200(41)	0.029(18)	0.765(153)	0.461(94)	-0.302(61)

Equations

$$⟨N(p^{\prime },s^{\prime })|j_{\mu }|N(p,s)⟩=\sqrt{\frac{m_N^2}{E_N(\vec{p}{}^{\prime })E_N(\vec{p})}}\times$$

$${\bar{u}}_N(p^{\prime },s^{\prime })\left[{\gamma }_{\mu }{F}_1(q^2)+\frac{i{\sigma }_{\mu \nu }{q}^{\nu }}{2m_N}{F}_2(q^2)\right]u_N(p,s).$$

$$j_{\mu }=\sum _f{e}_f{j}_{\mu }^f=\sum _f{e}_f{\bar{q}}_f{\gamma }_{\mu }{q}_f,$$

$$j_{\mu }^f(x)=\frac{1}{4}$$

$$G_E(q^2)=F_1(q^2)+\frac{q^2}{4m_N^2}{F}_2(q^2),$$

$$G_M(q^2)=F_1(q^2)+F_2(q^2).$$

$$\langle r^{2}\rangle=\frac{6}{G(0)}\frac{\partial G(q^{2})}{\partial q^{2}}\Big{|}_{q^{2}=0}\,.$$

$$J_N(\vec{x},t)={ϵ}^{abc}{u}^a(x)[u^{b⊺}(x)\mathcal{C}{\gamma }_5{d}^c(x)],$$

$$C({\Gamma }_0,\vec{p};t_s,t_0)=\sum _{{\vec{x}}_s}\mathrm{Tr}\left[{\Gamma }_0⟨J_N(t_s,{\vec{x}}_s){\bar{J}}_N(t_0,{\vec{x}}_0)⟩\right]e^{-i({\vec{x}}_s-{\vec{x}}_0)\cdot \vec{p}},$$

$$C_{\mu }({\Gamma }_{\nu },\vec{q},\vec{p}{}^{\prime };t_s,t_{\mathrm{ins}},t_0)=\sum _{{\vec{x}}_{\mathrm{ins}},{\vec{x}}_s}{e}^{i({\vec{x}}_{\mathrm{ins}}-{\vec{x}}_0)\cdot \vec{q}}{e}^{-i({\vec{x}}_s-{\vec{x}}_0)\cdot \vec{p}{}^{\prime }}\times$$

$$\mathrm{Tr}\left[{\Gamma }_{\nu }⟨J_N(t_s,{\vec{x}}_s)j_{\mu }(t_{\mathrm{ins}},{\vec{x}}_{\mathrm{ins}}){\bar{J}}_N(t_0,{\vec{x}}_0)⟩\right].$$

$$q_{\mathrm{smear}}^a(t,\mathbf{x})$$

$$F$$

$$H(\mathbf{x},\mathbf{y};U(t))$$

$$R_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p};t_s,t_{\mathrm{ins}})=\frac{C_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p};t_s,t_{\mathrm{ins}})}{C({\Gamma }_0,\vec{p}{}^{\prime };t_s)}\times$$

$$\sqrt{\frac{C({\Gamma }_0,\vec{p};t_s-t_{\mathrm{ins}})C({\Gamma }_0,\vec{p}{}^{\prime };t_{\mathrm{ins}})C({\Gamma }_0,\vec{p}{}^{\prime };t_s)}{C({\Gamma }_0,\vec{p}{}^{\prime };t_s-t_{\mathrm{ins}})C({\Gamma }_0,\vec{p};t_{\mathrm{ins}})C({\Gamma }_0,\vec{p};t_s)}},$$

$$R_{\mu }({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p};t_s;t_{\mathrm{ins}})\underset{t_{\mathrm{ins}}\gg 1}{\overset{t_s-t_{\mathrm{ins}}\gg 1}{\to }}{\Pi }_{\mu }({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p}).$$

$$L(t_{\mathrm{ins}},\vec{q})=\sum _{{\vec{x}}_{\mathrm{ins}}}\mathrm{Tr}\left[D^{-1}(x_{\mathrm{ins}};x_{\mathrm{ins}})\mathcal{G}\right]e^{i\vec{q}\cdot {\vec{x}}_{\mathrm{ins}}},$$

$$\frac{1}{N_r}\sum _r|{\xi }_r⟩⟨{\xi }_r|=1+\mathcal{O}\left(\frac{1}{\sqrt{N_r}}\right),$$

$$C(\vec{p},t_s)=c_0(\vec{p})e^{-E_0(\vec{p})t_s}+c_1(\vec{p})e^{-E_1(\vec{p})t_s},$$

$$C_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p},t_s,t_{\mathrm{ins}})=$$

$$A_{00}^{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p})e^{-E_0(\vec{p}{}^{\prime })(t_s-t_{\mathrm{ins}})-E_0(\vec{p})t_{\mathrm{ins}}}$$

$$+A_{01}^{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p})e^{-E_0(\vec{p}{}^{\prime })(t_s-t_{\mathrm{ins}})-E_1(\vec{p})t_{\mathrm{ins}}}$$

$$+A_{10}^{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p})e^{-E_1(\vec{p}{}^{\prime })(t_s-t_{\mathrm{ins}})-E_0(\vec{p})t_{\mathrm{ins}}}$$

$$+A_{11}^{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p})e^{-E_1(\vec{p}{}^{\prime })(t_s-t_{\mathrm{ins}})-E_1(\vec{p})t_{\mathrm{ins}}}.$$

$${\Pi }_{\mu }({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p})=\frac{A_{00}^{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p})}{\sqrt{c_0(\vec{p}{}^{\prime })c_0(\vec{p})}}.$$

$$R_{\mu }^{sum}({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p};t_s)$$

$$c+{\Pi }_{\mu }({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p})\times t_s+\cdots$$

$$d{R}_{\mu }^{sum}({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p};t_s)=$$

$$\frac{R_{\mu }^{sum}({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p};t_s+d{t}_s)-R_{\mu }^{sum}({\Gamma }_{\nu };\vec{p}{}^{\prime },\vec{p};t_s)}{d{t}_s}$$

$$G(Q^2)=\frac{G(0)}{(1+\frac{Q^2}{M^2}{)}^2},$$

$$⟨r^2⟩=\frac{12}{M^2}.$$

$$G(Q^2)=\frac{Q^{2}A}{4m_N^2+Q^{2}B}\frac{1}{(1+\frac{Q^2}{0.71{\mathrm{GeV}}^2}{)}^2},$$

$$⟨r^2⟩=-\frac{3A}{2m_N^2}.$$

$$G(Q^2)=\sum _{k=0}^{k_{\max }}{a}_k{z}^k,$$

$$z=\frac{\sqrt{t_{\mathrm{cut}}+Q^2}-\sqrt{t_{\mathrm{cut}}}}{\sqrt{t_{\mathrm{cut}}+Q^2}+\sqrt{t_{\mathrm{cut}}}}$$

Full text

Nucleon structure, Nucleon electromagnetic form factors, Disconnected, Lattice QCD

††preprint: DESY 18-033

I Introduction

Nucleons, being composite particles, have a non-trivial internal structure that can be probed by measuring their electromagnetic form factors. These fundamental quantities have been extensively studied both theoretically and experimentally. However, open issues still persist and there are on-going experimental efforts to determine them at higher precision and over a wider range of momentum transfers and to describe them theoretically. The proton electric form factor is extracted to high precision from electron proton scattering Bernauer et al. (2014). Its slope at vanishing momentum transfer squared yields the proton charge root mean square (r.m.s) radius. Prior to 2010, the charge r.m.s. radius of the proton was considered a well-determined quantity (see Ref. Punjabi et al. (2015) for a recent review). A pioneering experiment using Lamb shifts in muonic hydrogen surprisingly found a value smaller by five standard deviations Pohl et al. (2010), triggering the so-called proton radius puzzle. The origin of this discrepancy is not yet understood, and potential systematic uncertainties related to the analysis methodologies in the two types of experiments have not been excluded. Another quantity of interest is the neutron electric form factor Golak et al. (2001), which is accessed indirectly experimentally through electron-deuteron or electron-helium scattering and therefore remains poorly-known. It is of substantial importance to compute these fundamental quantities from first principles using lattice QCD, which provides an ideal formulation for such an investigation and with simulations at physical values of the QCD parameters.

Within this work, we compute the proton and neutron electromagnetic form factors including light quark disconnected contributions. We use an ensemble of twisted mass fermions with two degenerate light quarks, a strange and a charm quark ($N_{f}{=}2{+}1{+}1$) with masses fixed to their physical value (referred to hereafter as physical point). A clover term is added to the action to suppress isospin breaking effects that come quadratically with the lattice spacing. Details on the simulation can be found in Ref. Alexandrou et al. (2018a). We will refer to this ensemble as cB211.072.64. In addition, we present an analysis of a twisted mass ensemble of two degenerate light quarks with masses fixed to their physical values ($N_{f}=2$) to assess finite volume artifacts by comparing to previous results obtained using an $N_{f}=2$ ensemble with a smaller volume and same pion mass and lattice spacing Alexandrou et al. (2018b, 2017a). Comparison between $N_{f}=2$ and $N_{f}=2+1+1$ also sheds light on any possible unquenching effect of the strange and charm quarks. The momentum dependence of the form factors is fitted using two Ansätze, namely either a dipole or the Galster-like parameterization Galster et al. (1971) and the model independent z-expansion Hill and Paz (2010). The fits allow for the extraction of the magnetic moment and the electric and magnetic r.m.s radii of the proton and neutron and provide a measure of the systematics due to the choice of the fit method.

A crucial component of our analysis is the use of hierarchical probing Stathopoulos et al. (2013) combined with deflation of the lower lying eigenvalues Gambhir et al. (2017) that enables us to calculate the light quark disconnected contributions to the form factors at an unprecedented accuracy at the physical point. This allows us to obtain the proton and neutron form factors at the physical point without neglecting disconnected contributions.

The remainder of this paper is organized as follows: In Section II, we describe the nucleon matrix elements required to extract the electromagnetic form factors and in Section III we provide details on the lattice QCD techniques employed for the computation of the connected and disconnected diagrams. In Section IV, we discuss the analysis of the data paying particular attention to the identification of the ground state matrix element. In Section V we include an assessment of finite volume and unquenching effects using results from the analysis of the two $N_{f}{=}2$ ensembles. In Section VI, we fit the isovector and isoscalar form factors to extract the magnetic moments and radii. We compare to other lattice QCD studies using simulations close to physical pion masses in Section VII Capitani et al. (2015); Green et al. (2014); Hasan et al. (2018); Ishikawa et al. (2018). Our final results for the proton and neutron electromagnetic form factors are given in Section VIII. Finally, in Section IX, we summarize our findings and conclude. For completeness, we summarize in Appendix A the decomposition of the nucleon matrix elements in terms of the form factors and in Appendix B we provide a table with the numerical results for the electric and magnetic form factors as a function of the momentum transfer squared.

II Electromagnetic form factors

The nucleon matrix element of the electromagnetic current is parameterized in terms of the Dirac ($F_{1}$) and Pauli ($F_{2}$) form factors given in Minkowski space by,

[TABLE]

$N(p,s)$ is the nucleon state with initial (final) momentum $p$ ($p^{\prime}$) and spin $s$ ($s^{\prime}$), with energy $E_{N}(\vec{p})$ ($E_{N}(\vec{p}\,^{\prime})$) and mass $m_{N}$. $q^{2}{\equiv}q_{\mu}q^{\mu}$ is the momentum transfer squared $q_{\mu}{=}(p_{\mu}^{\prime}-p_{\mu})$ and $u_{N}$ is the nucleon spinor. The local vector current $j_{\mu}$ is given by,

[TABLE]

where $q_{f}$ is the quark field of flavor $f$ and $e_{f}$ its electric charge, and the summation runs over all the quark flavors. Instead of the local vector current, we instead use the symmetrized lattice conserved vector current given by

[TABLE]

which, unlike the local vector current, does not need renormalization. The electric and magnetic Sachs form factors $G_{E}(q^{2})$ and $G_{M}(q^{2})$ are alternative Lorentz invariant quantities and are expressed in terms of $F_{1}(q^{2})$ and $F_{2}(q^{2})$ via the relations,

[TABLE]

[TABLE]

In the isospin limit, where the up and down quarks are degenerate, we consider the isovector combination $\langle p|j_{\mu}^{u}-j_{\mu}^{d}|p\rangle$ that gives the difference between the proton and neutron form factors and the isoscalar combination $\langle p|j_{\mu}^{u}+j_{\mu}^{d}|p\rangle/3$ for the sum of the proton and neutron form factors. The electric form factor at zero momentum yields the nucleon charge, i.e. $G_{E}^{p}(0){=}1$ and $G_{E}^{n}(0){=}0$ which, when using the lattice conserved current, holds by symmetry, even prior to gauge averaging. The magnetic form factor at $q^{2}{=}0$ gives the magnetic moment, while the radii can be extracted from the slope of the electric and magnetic form factors as $q^{2}\rightarrow 0$, namely:

[TABLE]

III Calculation on the Lattice

III.1 Nucleon matrix element

Extraction of nucleon matrix elements within the lattice QCD formulation requires the evaluation of two- and three-point correlation functions in Euclidean space. We thus give all quantities in Euclidean space from here on. We use the standard nucleon interpolating field

[TABLE]

where $u$ and $d$ are up- and down-quark spinors and $\mathcal{C}{=}\gamma_{0}\gamma_{2}$ is the charge conjugation matrix. The two-point function in momentum space is given by

[TABLE]

and the three-point function is given by

[TABLE]

The initial position and time, $x_{0}$, is referred to as the source, the position and time of the current $j_{\mu}$ couples to a quark is denoted by $x_{\rm ins}$ and referred to as the insertion and the final position, $x_{s}$, as the sink. $\Gamma_{\nu}$ is a projector acting on spin indices, with $\Gamma_{0}{=}\frac{1}{2}(1{+}\gamma_{0})$ yielding the unpolarized and $\Gamma_{k}{=}\Gamma_{0}i\gamma_{5}\gamma_{k}$ the polarized matrix elements. Inserting complete sets of states in Eq. (8), one obtains the nucleon matrix element as well as additional matrix elements of higher energy states with the quantum numbers of the nucleon multiplied by overlap terms and time dependent exponentials. For large enough time separations, the excited state contributions are suppressed compared to the nucleon ground state and one can then extract the desired matrix element. In order to increase the overlap with the nucleon state and decrease overlap with excited states we use Gaussian smeared quark fields Alexandrou et al. (1994); Gusken (1990) for the construction of the interpolating fields:

[TABLE]

In addition, we apply APE-smearing to the gauge fields $U_{\mu}$ entering the hopping matrix $H$.

The Gaussian smearing parameters are optimized using the nucleon two-point function. We set $\alpha=0.2$ and $n=125$ Alexandrou et al. (2008). The values are: $\alpha=4.0$ and $n=50$, $70$ and $90$ for $\beta=3.9$, $4.05$ and $4.2$ respectively. For the APE smearing Albanese et al. (1987) we use 50 iteration steps and $\alpha_{APE}{=}0.5$.

An optimized ratio Alexandrou et al. (2013, 2011a, 2006) of the three-point function over a combination of two-point functions is used to cancel time dependent exponentials and overlaps, given by

[TABLE]

where $t_{s}$ and $t_{\rm ins}$ are taken to be relative to the source $t_{0}$ for simplicity. In the limit of large time separations, $(t_{s}{-}t_{\rm ins})\gg 1$ and $t_{\rm ins}\gg 1$, the lowest state dominates and the ratio becomes time independent

[TABLE]

$G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ are extracted from linear combinations of $\Pi_{\mu}(\Gamma_{\nu};\vec{p}\,^{\prime},\vec{p})$ as expressed in Appendix A, with $Q^{2}{\equiv}-q^{2}$ the Euclidean momentum transfer squared.

Contracting the quark fields in Eq. (8) gives rise to two types of diagrams depicted in Fig. 1, namely the so-called connected and disconnected contributions.

In the case of the connected diagram, the insertion operator couples to a valence quark and an all-to-all propagator arises between sink and insertion. We use sequential inversions through the sink that require keeping the sink-source time separation $t_{s}$, the projector, and the sink momentum $\vec{p}^{\prime}$ fixed. We perform additional sets of inversions to compute the three-point function for several values of $t_{s}$, for both the unpolarized and polarized projectors. We set $\vec{p}^{\prime}{=}\vec{0}$. We use an appropriately tuned multigrid algorithm Bacchio et al. (2018, 2016); Alexandrou et al. (2016) for the efficient inversion of the Dirac operator entering in the computation of the connected diagram. The disconnected diagram involves the disconnected quark loop correlated with the nucleon two-point correlator. The disconnected quark loop is given by

[TABLE]

where $D^{-1}(x_{\rm ins};x_{\rm ins})$ is the quark propagator that starts and ends at the same point $x_{\rm ins}$ and $\mathcal{G}$ is an appropriately chosen $\gamma$-structure. For the local vector current, which we use for the disconnected diagram, $\mathcal{G}=\gamma_{\mu}$. A direct computation of quark loops would need inversions from all spatial points on the lattice, making the evaluation unfeasible for our lattice size. We therefore employ stochastic techniques to estimate it combined with dilution schemes Wilcox (1999) that take into account the sparsity of the Dirac operator and its decay properties. Namely, in this work, we employ the hierarchical probing technique Stathopoulos et al. (2013), which provides a partitioning scheme that eliminates contributions from neighboring points in the trace of Eq. 12 up to a certain coloring distance $2^{k}$. Using Hadamard vectors as the basis vectors for the partitioning, one needs $2^{d*(k-1)+1}$ vectors, where $d{=}4$ for a 4-dimensional partitioning. Note that the computational resources required are proportional to the number of Hadamard vectors, and therefore in $d$=4 dimensions increase 16-fold each time the probing distance $2^{k}$ doubles. Contributions entering from points beyond the probing distance are expected to be suppressed due to the exponential decay of the quark propagator and are treated with standard noise vectors which suppress all off-diagonal contributions by $1/\sqrt{N_{r}}$, i.e.

[TABLE]

where $N_{r}$ is the size of the stochastic ensemble. Hierarchical probing has been employed with great success in previous studies Green et al. (2015, 2017) for an ensemble with a pion mass of 317 MeV. For simulations at the physical point, it is expected that a larger probing distance is required since the light quark propagator decays more slowly at smaller quark masses. We avoid the need of increasing the distance by combining hierarchical probing with deflation of the low modes Gambhir et al. (2017). Namely, we construct the low mode contribution to the light quark loops by computing exactly the 200 smallest eigenvalues and corresponding eigenvectors of the squared Dirac operator and combine them with the contribution from the remaining modes, which are estimated using hierarchical probing. Additionally, we employ the one-end trick McNeile and Michael (2006), used in our previous studies Alexandrou et al. (2014, 2017b, 2017c, 2017d) and fully dilute in spin and color.

III.2 Gauge Ensembles and Statistics

For the extraction of the electromagnetic form factors we analyze one $N_{f}{=}2{+}1{+}1$ Alexandrou et al. (2018a) and one $N_{f}=2$ ensemble. For both ensembles the quark masses are tuned to their physical values. The fermion action is the twisted mass fermion action with a clover term. Automatic ${\cal O}(a)$ improvement is achieved by tuning to maximal twist Frezzotti and Rossi (2004a, b). The $N_{f}=2+1+1$ cB211.072.64 ensemble is simulated using a lattice of size $64^{3}{\times}128$ with $Lm_{\pi}{=}3.62$ Alexandrou et al. (2018a), where $L$ is the spatial extent of the lattice. We determine the nucleon mass by fitting the effective mass in the large-time limit where the ground state dominates. The final value is chosen within a fit range where the value extracted is within half a standard deviation from the one determined by including in the fit the first excited state (two-state fit). The ratio of the nucleon to pion mass is $m_{N}/m_{\pi}=6.74(3)$ compared to the physical ratio of 6.8. Therefore, we use directly the average proton and neutron mass of 0.9389 GeV to set the scale. We find $a=0.0801(4)$ fm. For the pion mass we find $m_{\pi}=0.1393(7)$ MeV consistent with the average physical pion mass. Our current value of the lattice spacing is an update compared to the one given in Ref. Alexandrou et al. (2018a) using higher statistics where the values are consistent.

To assess finite volume effects, we use two $N_{f}{=}2$ ensembles, which only differ in their volume, namely one has $Lm_{\pi}=2.98$ and the other $Lm_{\pi}{=}3.97$. We will refer to them as the cA2.09.48 and cA2.09.64 ensembles, respectively. We note that since the pion mass is not exactly at the physical value we interpolate to the physical pion mass using one-loop chiral perturbation theory. We include a systematic error on the extracted lattice spacing, determined as the difference in the mean value obtained using one-loop chiral perturbation theory and heavy baryon chiral perturbation theory. This systematic error on the lattice spacing appears for the two $N_{f}{=}2$ ensembles, while it is absent in the case of the cB211.072.64 ensemble. Results on the form factors for the $N_{f}=2$ ensemble with $Lm_{\pi}{=}2.98$ are from Ref. Alexandrou et al. (2017a) while results for the other two ensembles are reported here for the first time. The simulation parameters of all three ensembles considered in this work are tabulated in Table 1.

For the analysis of the cB211.072.64 ensemble we use 750 configurations separated by 4 trajectories. For the connected contributions we evaluate the three-point function for five sink-source time separations in the range 0.96 fm to 1.60 fm increasing the number of source positions per configuration as we increase the time separation so as to keep the statistical error approximately constant. In Table 2 we give the statistics used in the calculation of the connected three-point functions.

For the evaluation of the disconnected contributions we use $N_{\rm srcs}=200$ source positions to generate the nucleon two-point functions that are correlated with the quark loop to produce the disconnected contribution to the three-point function. We find that the volume is sufficiently large so that the data extracted from this large number of randomly distributed source positions on the same configuration is statistically independent. Nevertheless, we average over all source positions for each configuration and take the averaged correlation function as one statistic in our jackknife error analysis. As mentioned in the previous section, for the evaluation of the light quark loops we use the first 200 low modes of the squared Dirac operator to reconstruct exactly part of the loop. The contribution from the high modes is estimated stochastically using one noise vector per configuration combining hierarchical probing, one-end trick and spin-color dilution. For the hierarchical probing we use distance eight coloring resulting in 512 Hadamard vectors, which when combined with spin-color dilution leads to 6144 inversions per configuration. We note that the next coloring distance would demand 8192 Hadamard vectors, resulting in 98304 inversions per configurations after combining with spin-color dilution, making such a computation more than an order of magnitude more expensive.

For the computation of the disconnected contributions for $N_{f}=2$ cA2.09.48 ensemble computed previously we used only the one-end trick and 2250 noise vectors for the calculation. Two-point functions were computed for 100 source positions. More detail can be found in Ref. Alexandrou et al. (2018b). in Table 3 we summarize the parameters for the computation of the disconnected three-point functions.

The cA2.09.64 ensemble is used to check for finite volume effects, comparing the connected contributions to those of the cA2.09.48 ensemble. For the latter, the setup is reported in Ref. Alexandrou et al. (2017a) and summarized in Table 1. For the larger lattice size ensemble, we analyze three sink-source time separations in the range of 1.1 fm to 1.5 fm. We fix the number of source positions per configuration to 16 and we use more configurations for the larger time separations to control statistical error. In Table 4 we summarize the statistics for both $N_{f}=2$ ensembles.

III.3 Excited states contamination

Assessment of excited state effects is imperative for the proper extraction of the desired nucleon matrix element. However, ensuring ground state dominance is a delicate process due to the exponentially increasing statistical noise with increasing sink-source separation. We use four methods to study the effect of excited states and identify the final results based on a critical comparison among these methods. Only by employing these different methods can one reach a reliable assessment of excited state contributions and extract the nucleon matrix element of interest. The methods employed are as follows:

Plateau method: In this method we use the ratio in Eq. (10) and identify a time-independent window (plateau) as we increase $t_{s}$. The converged plateau value then yields the desired matrix element.

Two-state method: Within this method we fit the two- and three-point functions considering contributions up to the first excited state using the expressions

[TABLE]

[TABLE]

In Eqs. (14) - (15) $E_{0}(\vec{p})$ and $E_{1}(\vec{p})$ are the energies of the ground and first excited states with total momentum $\vec{p}$, respectively. The ground state corresponds to a single particle state and therefore one can use the continuum dispersion relation, $E_{0}(\vec{p})=\sqrt{\vec{p}^{2}+m_{N}^{2}}$, with $\vec{p}=\frac{2\pi}{L}\vec{n}$ with $\vec{n}$ a lattice vector with components $n_{i}\in(-\frac{L}{2a},\frac{L}{2a}]$. The continuum dispersion relation is satisfied for all $Q^{2}$ values considered in this work. The first excited state, on the other hand, can be a two-particle state. We, thus, fit simultaneously the two-point functions with momenta $\vec{p}$ and $\vec{p}^{\prime}$ and the three-point function involving in total eleven parameters. Note that for non-zero momentum transfer, $A_{01}^{\mu}(\Gamma_{\nu},\vec{p}\,^{\prime},\vec{p})\neq A_{10}^{\mu}(\Gamma_{\nu},\vec{p}\,^{\prime},\vec{p})$. This allows us to extract the matrix element given by

[TABLE]

Summation method: Summing over $t_{\rm ins}$ in the ratio of Eq. (10) yields a geometric sum  Maiani et al. (1987); Capitani et al. (2012) from which we obtain,

[TABLE]

where the ground state contribution, $\Pi_{\mu}(\Gamma_{\nu};\vec{p}\,^{\prime},\vec{p})$, is extracted from the slope of a linear fit with respect to $t_{s}$. The sink-source time separation $t_{s}$ considered in the fit should be large enough to suppress higher order contributions.

Derivative Summation method: Instead of performing a linear fit in Eq. (17) to extract the matrix element, one can take finite differences to the summed ratio Savage et al. (2017) as follows

[TABLE]

and fit to a constant to extract the desired matrix element.

IV Analysis of lattice results

IV.1 Isovector and connected isoscalar form factors

The isovector combination gives the difference between the proton and neutron form factors, and in this case, only the connected diagram contributes since disconnected contributions cancel, up to cut-off effects of ${\cal O}(a^{2})$. For the connected diagram we use a frame where the nucleon final momentum is zero, thus $\vec{q}{=}{-}\vec{p}$. In Figs. 2 and 3 we show the ratios defined in Eq. (10) as a function of the sink-source time separation, and for three values of the momentum transfer squared, that is for $Q^{2}{=}0.057$ GeV2, $Q^{2}{=}0.219$ GeV2 and $Q^{2}{=}0.554$ GeV2. In a frame where the final momentum of the nucleon $\vec{p}^{\prime}$ is zero, the expressions in Appendix A given by Eqs. (39) and (40) reduce to Eqs. (42,43), and (44), giving separately the electric and magnetic form factors. We note that for the electric form factor Eq. (42) leads to much more precise results compared to Eq. (43) and therefore we use only Eq. (42). In the case of the ratio determining $G_{E}^{u-d}(Q^{2})$, as $t_{s}$ increases, the plateau value decreases with larger deviations as $Q^{2}$ increases. This shows that at larger $Q^{2}$ values contamination due to excited states is more severe. In the case of $G_{M}(Q^{2})$, excited states are suppressed and only a small variation with $t_{s}$ is observed.

We further investigate effects due to excited states by employing the summation and two-state fits. In Fig. 4 we show linear fits to the summed ratio for three different values of $Q^{2}$. The slope gives the nucleon matrix element. All three momenta follow well the linear behavior, within the statistical error, indicating that contributions from higher order terms are suppressed. In the right panel of Fig. 4, we demonstrate the plateaus for the derivative summation method, fitting to a constant to extract the matrix element of the ground state. Within statistical accuracy all three momenta are indeed flat, and are thus described well by a constant.

In Fig. 5 we show the results extracted using two-state fits for both electric and magnetic form factors. The data correspond to the ratio of Eq. (10) and the curves are obtained by fitting simultaneously the three- and two-point functions to Eqs. (15) and (14). The gray horizontal band shows the nucleon matrix element value and error extracted from the two-state fit as in Eq. (16). For the electric form factor, the ratio shows a trend towards lower values as we increase the sink-source separation, with $t_{s}/a{=}20$ becoming compatible with the value extracted from the two-state fit. In the case of the magnetic form factor, the value extracted from the two-state fit is compatible with the ratio for all time separations considered confirming the weak dependence of the matrix element on the sink-source time separation observed in the plateau method.

In Fig. 6 we show the extracted values for the matrix element yielding the isovector electromagnetic form factors. We compare the plateau, summation, derivative summation and two-state fit methods. For the plateau method we show the value extracted from the constant fit for all sink-source separations available. For the other cases we vary the lower fit range, keeping the upper range fixed to $t_{s}/a{=}20$. We seek for the earliest agreement between the plateau method and the other three cases. As already pointed out, the isovector electric form factor shows more severe excited state effects for large $Q^{2}$-values and we therefore take the largest time separation for the plateau method to fulfill our criterion for agreement with the other methods. For the isovector magnetic form factor, although excited state effects are mild, we still observe a shift to larger values for the smallest $Q^{2}$ and, therefore, we conservatively use the largest time separation available also in this case. An additional observation is that summation and derivative summation methods produce compatible results with similar accuracy, as can be seen in Fig. 6, and thus from now on we will restrict to showing results only from the summation method.

In Fig. 7, we present our results for $G_{E}^{u-d}(Q^{2})$ and $G_{M}^{u-d}(Q^{2})$ as a function of the momentum transfer squared $Q^{2}$. We limit the plot to up $Q^{2}{=}0.5$ GeV2 to make visible the values extracted using the plateau at the four largest separations, the summation, and the two-state fit approaches. For the summation and two-state fit we show the values indicated by the filled symbols in Fig. 6. As can be seen, for the electric form factor effects of excited states are small for small values of $Q^{2}$ but become more severe for higher $Q^{2}$ values, with the extracted value decreasing with increasing time separation in line with the observation made in Fig. 2. For $G_{M}^{u-d}(Q^{2})$, excited state effects are small for larger $Q^{2}$ values whereas for smaller $Q^{2}$ values there is a systematic increase in the values of the form factor with the time separation.

For the extraction of the connected isoscalar form factors we follow a similar analysis procedure as in the isovector case. In Fig. 8 we present the connected contribution to isoscalar electric and magnetic form factors comparing the plateau, summation and two-state fit methods. Excited states have a smaller effect on the isoscalar form factors being detectable only for the magnetic at small values of $Q^{2}$ where the two-state fit yields systematically larger values. Given that there is agreement between the plateau values for the largest time separation and the two-state fit we will use the plateau value as the final result for the form factors. The deviation from the values determined from the two-states fits are then taken as an estimate of the systematic error due to excited states. Since we will be using the plateau values in the case of disconnected since two-state fits are not stable in that case we do the same for the connected for consistency.

IV.2 Disconnected contributions

A major component of this work is the evaluation of the disconnected contributions shown diagrammatically in Fig. 1 that enter in the evaluation of the isoscalar as well as in the proton and neutron form factors.

The disconnected quark loops are computed using the formalism described in Section III.2 with the statistics summarized in Table 3. As already discussed, the hierarchical probing method, combined with deflation of the low eigenmodes, provides an accurate determination of the diagonal of the quark propagator entering in the evaluation of the quark loops. It is thus preferable to use the local vector current for the evaluation of the disconnected contributions since the conserved current includes non-diagonal terms. We therefore need the renormalization function $Z_{V}$, which is determined non-perturbatively, in the RI*′*-MOM scheme, employing momentum sources. We perform a perturbative subtraction of $\mathcal{O}(g^{2}a^{\infty})$-terms, as described in Refs. Alexandrou et al. (2017e, 2011b), which subtracts the leading cut-off effects leaving only a weak dependence on the renormalization scale $(a\mu)^{2}$, as shown in Fig. 9. We find a value of $Z_{V}{=}0.728(1)$ where the error is statistical. Alternatively, $Z_{V}$ can be determined at $Q^{2}=0$, by taking the ratio of $G^{u-d}_{E}(0)$ computed with the local current to $G^{u-d}_{E}(0)$ computed using the lattice conserved current. This ratio yields a value of 0.715(3). Although this is 2% smaller than $Z_{V}$ as determined from the vertex function, the difference between them is still an order of magnitude smaller as compared to the statistical errors for the disconnected contributions. In what follows we use $Z_{V}{=}0.728(1)$ to renormalize the matrix elements computed using the local current, since this determination has taken into account higher order cut-off effects as compared to the one determined from the ratio. We note that $Z_{V}$ only enters in the disconnected three-point function. A more detailed description of the renormalization procedure including other renormalization functions will be provided in a future publication.

Disconnected quark loops are evaluated for every time-slice allowing us to compute the three-point function for every combination of $t_{s}$ and $t_{\rm ins}$. As in the case of the connected, we are seeking for a reasonable window in $t_{s}$ to extract the nucleon matrix elements, where excited states are sufficiently suppressed and noise is not prohibitively large. In contrast to the connected diagram, where we have results only for the case $\vec{p}\,^{\prime}{=}\vec{0}$, for the disconnected diagrams we have all sink momenta at no additional cost. We analyze, besides $\vec{p}\,^{\prime}{=}\vec{0}$, the matrix element for the six final momenta with $\vec{p}\,^{\prime}{=}\pm\frac{2\pi}{L}\hat{n}$, with $\hat{n}=\hat{x}$, $\hat{y}$, or $\hat{z}$, i.e. the unit vector in one of the three spatial directions. Given that the statistical errors in the case of the disconnected diagrams are larger as compared to the connected diagrams, we restrict ourselves in using the plateau method for different values of $t_{s}$ in order to check for ground state dominance. This is because the two-state fits are problematic given the larger errors of the disconnected diagrams. Whenever they work they yield large errors and are consistent with the plateau extraction.

In Fig. 10 we present our results for the disconnected contributions to $G_{E}^{u+d}(Q^{2})$ and $G_{M}^{u+d}(Q^{2})$ up to $Q^{2}{=}1$ GeV2 for three time separations, $t_{s}=0.96,1.12$ and 1.28 fm. We can achieve a relative statistical error that is less than 20$\%$ for up to $t_{s}{=}14a{\sim}1.12$ fm, which is unprecedented given that we are using a physical pion mass ensemble. As we increase the time separation from $t_{s}/a{=}12$ to $t_{s}/a{=}14$ we observe, for both $G_{E}^{u+d}(Q^{2})$ and $G_{M}^{u+d}(Q^{2})$, that there is a trend for larger values, while the results extracted for $t_{s}/a{=}16$ are in a very good agreement with those extracted for $t_{a}=14a$ for most $Q^{2}$ values, albeit with larger errors. We, therefore, take as our final result for the disconnected contribution the value extracted using $t_{s}/a{=}14$ for both $G_{E}^{u+d}(Q^{2})$ and $G_{M}^{u+d}(Q^{2})$. We use the difference between the central value of the results at $t_{s}=14a$ and $t_{s}=16a$ as an estimate of the systematic error from excited state effects when we quote quantities that include disconnected contributions.

V Assessment of lattice artifacts

We collectively discuss here lattice artifacts that may lead to systematic errors. Since we use simulations with physical values of the light quark masses no chiral extrapolation is needed eliminating one of the biggest uncertainty present in past lattice QCD computations of these quantities.

• •

Disconnected contributions: The main novelty of this work is the accurate computation of the light quarks disconnected contributions using simulations with quark masses tuned to their physical values. This enables us, for the first time, to eliminate this systematic uncertainty in the determination of the proton and neutron form factors at the physical point. Strange quarks loop contribution is not included in this study, but we know from previous studies Alexandrou et al. (2018b); Sufian et al. (2017a); Green et al. (2015) that it is much smaller compared to the statistical error of the connected contribution.

• •

Quenching effects: The analysis of the $N_{f}=2+1+1$ ensemble and two $N_{f}=2$ ensembles allows us to check for unquenching effects due to the strange and charm quarks. In Figs. 11 and 12 we compare results using the $N_{f}{=}2$ cA2.09.48 ensemble to the $N_{f}{=}2{+}1{+}1$ cB211.072.64 ensemble. The results are extracted from the plateau method with time separation $t_{s}=1.3$ fm for both isovector and isoscalar electric and magnetic form factors. We observe consistent results between the two ensembles. Therefore, to the accuracy of our data, no quenching effects due the strange and charm quarks can be detected. This corroborates our previous study where we found consistent results when comparing an $N_{f}=2$ and an $N_{f}=2+1+1$ ensemble at a pion mass of about 370 MeV Alexandrou et al. (2011a).

• •

Isolation of ground state matrix element: An analysis of excited state contributions is carried out by performing the calculation using several time separations $t_{s}$. For the target $N_{f}=2+1+1$ ensemble we use five values of $t_{s}$ tabulated in Table 2. We probe convergence to the ground state matrix element by demanding that the matrix element extracted using the plateau and two-state fits are consistent, as explained in detail in Section III.3. The value of $t_{s}\sim 1.6$ fm is the largest utilized in this study and to our knowledge in any other study at the physical point. We increase the statistics as $t_{s}$ increases to keep the errors under control so that a meaningful analysis can be performed to isolate the ground state matrix element. For all of our results we observe agreement between the values extracted using the plateau and two-state fits. Despite this agreement, residual contamination can still lead to a systematic error within our current statistical. We give an estimate of such a systematic error by comparing the values obtained with the plateau and two-state fits.

• •

Finite volume effects: For the assessment of finite volume effects we compare the two $N_{f}{=}2$ physical point ensembles cA2.09.48 and cA2.09.64 that yield respectively $m_{\pi}L\simeq 3$ Alexandrou et al. (2017a) and $m_{\pi}L\simeq 4$. The lattice spacing and pion mass are the same for these two ensembles. We also use the same time separation $t_{s}$ for each observable when comparing between the two ensembles. The isovector electric and magnetic form factors extracted using the plateau method are shown in Fig. 13. The results fall on the same curve indicating no significant finite volume effects between the two volumes of $m_{\pi}L\simeq 3$ and $m_{\pi}L\simeq 4$. The same behavior is observed for the isoscalar form factors shown in Fig. 14.

We would like to stress once more that our statement of detecting no volume effects can only be made within the current accuracy and some residual volume effects can still lead to a systematic effect. One complication as the volume increases is the contamination due to higher excited states, since the number of multi-hadron states allowed increases Bar (2019). Such multi-hadron states are not expected to affect results at larger pion masses but are expected to be more severe at the physical pion mass. Such effects can be modeled within chiral perturbation theory for the axial form factors Bar (2019). For the electromagnetic form factors these effects are not known but an interplay between volume and excited state effects may account for the deviations observed between lattice QCD data and experimental values.

• •

Finite lattice spacing, a: Since in this work we are using the twisted mass formulation at maximal twist our results are automatically ${\cal O}(a)$ improved without any need to improve the current. This is different from clover fermions where the current must be improved in order to eliminate oder $a$ contributions. Therefore, our results only have corrections of ${\cal O}(a^{2})$. Continuum extrapolation cannot be performed given that we have analyzed only one $N_{f}=2+1+1$ ensemble. The two $N_{f}=2$ ensembles analyzed have the same lattice spacing and so again finite lattice spacings affects cannot be assessed. Previous studies done using ensembles with pion mass spanning about 460 MeV to 260 MeV and three values of the lattice spacings have indeed demonstrated that the ${\cal O}(a^{2})$ correction is negligible Alexandrou et al. (2011a). We thus do not expect large systematic cut-off effects on our results. However, an analysis of cut-off effects will need to be carried out in the future when additional ensembles are available.

In summary, there maybe a slow convergence as a function of the volume in conjunction with residual excited state effects. This may explain the few $\sigma$ discrepancy observed between lattice QCD results and the experimental values. In particular, we note that the electric form factors has increasing excited state effects for larger values of $Q^{2}$, whereas for $G_{M}$ these effects are bigger at small $Q^{2}$. As we will see these are the ranges of momenta where we see discrepancies with the experimental values.

VI $Q^{2}$-dependence of the isovector and isoscalar form factors

In this section we discuss the $Q^{2}$-dependence of the form factors using standard parameterizations as described in the next section.

VI.1 Parameterizations of the $Q^{2}$-dependence

Assuming vector meson pole dominance for $Q^{2}<0$, one expects that for small $Q^{2}>0$ the behavior will be dominated by the poles in the time-like region. One would then expect a dipole form given by Perdrisat et al. (2007)

[TABLE]

where $M$ is the mass of the vector meson that parameterizes the $Q^{2}$ dependence. The value of the form factor at zero momentum transfer gives the electric charge in the case of the electric form factor and the magnetic moment in the case of the magnetic form factor. Combining Eq. (19) and Eq. (5) one can relate $M$ to the mean square radius as

[TABLE]

The neutron electric form factor and disconnected contributions to the electric form factors are zero for $Q^{2}{=}0$ and we treat them as special cases, fitting them using the Galster-like parameterization Galster et al. (1971); Alberico et al. (2009), given by

[TABLE]

with $A$ and $B$ fit parameters. In this case the radius is given by

[TABLE]

Another fit form, which has been applied recently to experimental data of both electromagnetic and axial form factors, is the model independent z-expansion Hill and Paz (2010). In this case, the form factor is expanded in a series given by,

[TABLE]

where

[TABLE]

and $t_{\rm cut}$ is the time-like cut of the form factor. We take $t_{\rm cut}{=}4m_{\pi}^{2}$ for the isovector combination $G^{u}-G^{d}$ and $t_{\rm cut}{=}9m_{\pi}^{2}$ for the isoscalar combination $G^{u}+G^{d}$ Hill and Paz (2010). For convergence of the truncated series of Eq. (23), the coefficients $a_{k}$ should be bounded in size and convergence should be demonstrated by increasing $k_{\rm max}$. The interested reader is referred to Ref. Alexandrou et al. (2018b) for details about our procedure. The mean square radius is given by

[TABLE]

while the value of the form factor at zero momentum transfer is $G(0){=}a_{0}$.

VI.2 Fits to lattice QCD results

We consider first the isovector form factors where only the connected diagram contributes. In Figs. 15 and 16 we show fits using the dipole form, comparing between results from the plateau method at $t_{s}/a=20$ and from two-state fits for $G_{E}^{u-d}(Q^{2})$ and $G_{M}^{u-d}(Q^{2})$, respectively. As can be seen, fits using the plateau and two-state methods are fully consistent and do not show any significant systematic effect on the determination of the $Q^{2}$-dependence of the form factors, indicating that excited states are sufficiently suppressed. Since results are in agreement, from now on we will use the plateau method at $t_{s}/a=20$ to extract final results on the form factors, r.m.s radii, and magnetic moment. We will use the results extracted from the two-state fits to estimate the systematic error due to excited states.

In Figs. 17 and 18, we show fits to $G_{E}^{u-d}(Q^{2})$ and $G_{M}^{u-d}(Q^{2})$, respectively using the dipole form and the z-expansion of Eq. (23) and compare to experiment. For the z-expansion, we check convergence by increasing $k_{\rm max}$. The resulting magnetic moment and r.m.s radii are shown in Fig. 19, where we observe convergence for $k_{\rm max}{=}4$. For $G_{E}^{u-d}(Q^{2})$, we see from Fig. 17 that the slope of the lattice QCD data is less as compared to the experimental values. Therefore, although both dipole form and z-expansion describe very well our data, shown in separate panels for clarity, they lie consistently above the experimental values. A study using a larger volume with a careful examination of excited state effects is planned to understand the origin of this remaining discrepancy. In extracting the r.m.s radius, we see from Fig. 19 that results obtained from using the dipole fit and z-expansion are compatible, and yield

[TABLE]

where the central value and the statistical error are taken from the dipole fit, the second error is a systematic computed as the difference in the mean values between dipole and z-expansion for $k_{\rm max}{=}4$ and the third error is the systematic error due to excited states obtained from the difference when fitting the form factor extracted from the plateau and from the two-state fit method. Subsequent quantities given in the paper will have statistical and systematic errors quoted using the same convention as in Eq. 26.

For $G_{M}^{u-d}$, shown in Fig. 18, we observe that our results are in agreement with the experimental values for $Q^{2}>0.2$ GeV2, whereas for small $Q^{2}$ they tend to be lower. A possible explanation for this discrepancy is that effects from the pion cloud, expected to be prominent for small momenta Sato et al. (2003), are suppressed in our calculation due to our finite volume. The fact that we have seen no volume effects when we increase the volume from $Lm_{\pi}{\simeq}3$ to $Lm_{\pi}{\simeq}4$ for our two $N_{f}=2$ ensembles may indicate that pion cloud suppression may not be detectable for these volume sizes requiring larger volumes to unfold. Indeed, preliminary results by PACS using a physical point ensemble with $Lm_{\pi}{\simeq}7.4$ Shintani et al. (2018) finds a higher value that may point to a finite volume effect. This would need further investigation to confirm.

The isovector magnetic moment and mean square magnetic radii are shown in Fig. 19. As can be seen, the mean values extracted for $\mu^{u-d}$ using the dipole and z-expansion are compatible, while for $\langle r_{M}^{2}\rangle^{u-d}$ the z-expansion produces a slightly higher mean value, which, however, is consistent within errors. Quoting the values from the dipole fit, we find

[TABLE]

Here we have included a fourth systematic error computed as the difference in the values of $\mu^{u-d}$ and $\langle r_{M}^{2}\rangle^{u-d}$ when fitting $G_{M}^{u-d}(Q^{2})$ including and excluding the lowest $Q^{2}$ value from the fit. The error is asymmetric, since the expectation is that pion cloud effects will increase the value of the magnetic form factor. It is also small compared to the systematic error due to excited states. In what follows we will not include this fourth systematic error.

Before presenting fits to the total isoscalar form factors we discuss separately the $Q^{2}$-dependence of the disconnected contributions. In Fig. 20 we show the disconnected contribution to the isoscalar electric form factor $G_{E}^{u+d}$, accompanied by fits to the Galster-like parameterization and z-expansion. We note that in the case of the z-expansion we take $a_{0}{=}0$, since $G_{E}^{u+d}(0){=}0$ for the disconnected contribution. Both parameterizations describe well our results with the z-expansion yielding a larger error for the larger $Q^{2}$ values.

The disconnected contribution to $G_{M}^{u+d}(Q^{2})$ is shown in Fig. 21. We find that both dipole and z-expansion are in good agreement. In particular, they yield compatible values at zero momentum transfer. Like in the case of the disconnected contribution to $G_{E}^{u+d}(Q^{2})$, for large $Q^{2}$ the dipole fit has a smaller error band as compared to the z-expansion. The values extracted from fitting the disconnected contributions alone are

[TABLE]

where we have not normalized with the value of the form factor at zero momentum transfer, i.e. the radii are extracted from \langle r^{2}\rangle=-6\frac{\partial G(Q^{2})}{\partial Q^{2}}\Big{|}_{Q^{2}=0}\, rather than from Eq. 5.

In Fig. 22 we show the isoscalar form factors when including and excluding disconnected contributions. Although the effect is small for both $G_{E}^{u+d}(Q^{2})$ and $G_{M}^{u+d}(Q^{2})$ there is a systematic shift affecting the parameters of the fits. This comparison shows that disconnected contributions although small are important to include and that their omission would result in an uncontrolled systematic error comparable to the statistical uncertainty. Such systematics need to be under control for precision results required for distinguishing e.g. the two experimental determinations of the charge radius of the proton.

In Figs. 23 and 24 we show the fits of the total isoscalar electric and magnetic form factors using the dipole form and z-expansion. Both fits describe well the data with the dipole fit being more precise at larger $Q^{2}$, a behavior also observed for the isovector form factors. For intermediate $Q^{2}$ values our results are systematically higher compared to experiment, which is then reflected in the fit bands. Since for low $Q^{2}$ there is agreement, the extracted value for the isoscalar magnetic moment agrees with the experimental value. On the other hand, the slope of our lattice data is not as steep as in the experimental results, which leads to a smaller value for the corresponding radii.

In the top panel of Fig. 25 we show the isoscalar electric square radius. As can be seen, the z-expansion fit yields values that are within errors for $k_{\rm max}>1$ but with twice larger errors than the dipole. In Fig. 25 we also show results for the magnetic moment and the magnetic radius where convergence of the z-expansion is observed already for $k_{\rm max}{=}2$. In general, there is agreement between the results extracted from the dipole and z-expansion. In what follows we will quote the values determined from the dipole fits and quote as a systematic error the difference between the mean values of the dipole and the z-expansion fits. We find

[TABLE]

Note that from our definition of the isoscalar combination, the proton plus neutron magnetic moment is obtained by: $\mu^{p}+\mu^{n}=\mu^{u+d}/3$.

VII Comparison with other studies

Before we discuss our final results for the proton and neutron form factors we compare with results by other groups using different lattice QCD ensembles and discretization schemes. These mainly exist for the isovector electromagnetic form factors allowing us to qualitatively assess lattice artifacts. This is useful since most groups use a single ensemble and thus infinite volume and continuum extrapolations are lacking. We summarize the lattice QCD discretized actions used by different groups for the computation of the electromagnetic form factors, restricting ourselves only to published works and results that were obtained using simulations with pion mass less than 170 MeV:

• •

LHPC analyzed one ensemble of $N_{f}{=}2{+}1$ with two levels of HEX-smeared clover fermions with $m_{\pi}{=}149$ MeV, lattice spacing $a{=}0.116$ fm and $Lm_{\pi}{=}4.21$ at three sink-source time separations from 0.93 fm to 1.39 fm Green et al. (2014). They give as their final results the ones extracted using the summation method, which leads to larger statistical errors. Additionally, they analyzed an $N_{f}{=}2{+}1$ ensemble with two levels of HEX-smeared clover fermions with $m_{\pi}=135$ MeV, lattice spacing $a{=}0.093$ fm and $Lm_{\pi}{=}4$ Hasan et al. (2018). They analyzed three lattice separations from 0.93 fm to 1.5 fm and they have extracted results using the summation method. A momentum derivative method has been used to extract the magnetic moment and the electric radius directly from the correlation functions avoiding a fitting procedure.

• •

The PACS collaboration analyzed an ensemble of $N_{f}{=}2{+}1$ stout-smeared clover fermions with $m_{\pi}{=}146$ MeV, $a\simeq 0.085$ fm and a spatial extend of 8.1 fm or $Lm_{\pi}\simeq 6$ allowing access to relatively small momenta Ishikawa et al. (2018); Shintani et al. (2018). PACS has computed three-point functions for one sink-source time separation of 1.27 fm and they used the plateau method to identify the ground state matrix element.

• •

The $\chi$QCD collaboration Sufian et al. (2017b) computed only the disconnected contributions to the nucleon electromagnetic form factors using a hybrid action of overlap valence quarks and $N_{f}{=}2{+}1$ domain wall sea quarks produced by RBC/UKQCD. Their analysis includes an ensemble with pion mass $m_{\pi}=139$ MeV, $a{=}0.1141(2)$ fm and $Lm_{\pi}{=}3.86$. They computed nucleon two-point functions stochastically using $Z_{3}$-noise grid sources and disconnected quark loops with $Z_{4}$-noise grids applying even-odd and time dilution as well as low-mode average.

• •

Our results obtained using the three ensembles of Table 1 simulated by the Extended Twisted Mass Collaboration (ETMC). These include the two analyses of this work, namely the $N_{f}$=2+1+1 ensemble with $m_{\pi}$=139 MeV, $a$=0.0801(4)(3) fm and $Lm_{\pi}{\simeq}3.6$ and the $N_{f}{=}2$ cA2.09.64 ensemble with $m_{\pi}{=}130(1)$ MeV, $a{=}0.0938(3)(1)$ fm and $Lm_{\pi}{\simeq}4$ as well as our results from Refs. Alexandrou et al. (2018b, 2017a), which were obtained using the $N_{f}{=}2$ cA2.09.48 ensemble with $Lm_{\pi}{\simeq}3$ and same pion mass and lattice spacing.

In Fig. 26 we show a comparison of lattice QCD results for $G_{E}^{u-d}(Q^{2})$ up to $Q^{2}{=}0.5$ GeV2 from the analyses mentioned above. As can be seen, ETMC and PACS results are in good agreement but systematically higher than the experimental values. LHPC results were obtained using the summation method and in general have larger statistical errors making them compatible with both our results and the experimental values.

In Fig. 26, we also show the corresponding results for $G_{M}^{u-d}(Q^{2})$. The ETMC results of this work are the most precise and in good agreement with those obtained from other studies. We note the very good agreement of lattice QCD results and experiment for $Q^{2}{>}0.2$ GeV2. As pointed out, the underestimation of lattice QCD results compared to experimental values at smaller $Q^{2}$ may indicate that a larger spatial volume is required to develop fully the pion contributions. Although our study using two ensembles of $N_{f}=2$ showed no detectable volume effects when we increase the spatial extent from 4.5 fm to 6 fm (or equivalently from $Lm_{\pi}\sim 3$ to $Lm_{\pi}\sim 4$) the volume dependence could be weak and require a larger volume to manifest itself. The new PACS results may indicate such a trend Shintani et al. (2018). A conclusion that we can, however, draw from these lattice QCD studies is that there is agreement among them for both the electric and magnetic form factors. Given the different discretization schemes employed, this agreement indicates that cut-off effects are smaller than the statistical errors.

In Fig. 27 we show a comparison of the disconnected contributions to $G_{E}^{u+d}$ and $G_{M}^{u+d}$ using results obtained from our $N_{f}=2+1+1$ and $N_{f}=2$ twisted mass ensembles and from the hybrid action as analyzed by the $\chi$QCD collaboration Sufian et al. (2017b). We would like to stress the accuracy of the results of the current work using the $N_{f}$=2+1+1 twisted mass ensemble. In our previous evaluation of the disconnected contributions for the $N_{f}=2$ twisted mass ensemble we used 2120 configurations with 100 source positions for the computation of the two-point functions and 2250 stochastic vectors for the disconnected loops Alexandrou et al. (2018b). This is approximately the same number of inversions (and thus cost) as for the $N_{f}$=2+1+1 ensemble (see Table 3), which demonstrates the effectiveness of the hierarchical probing method employed in the current analysis of the $N_{f}$=2+1+1 ensemble.

The proton and neutron form factors can be extracted from the isovector and isoscalar form factors discussed in Section VI, using the linear combinations

[TABLE]

In Fig. 28, we show lattice QCD results for the proton electromagnetic form factors. To extract these, one needs both the isovector and isoscalar combinations. The latter includes disconnected contributions, which have only been computed by ETMC for ensembles with physical pion masses. We still provide a comparison with the lattice results by LHPC which however do not include these disconnected contributions. We use filled symbols to indicate lattice QCD results that include disconnected contributions. For both the proton electric and magnetic form factors LHPC results are in agreement with ours, with the LHPC results exhibiting larger errors due to the usage of the summation method. The accurate ETMC results are higher than the experimental values for $G_{E}^{p}(Q^{2})$, while for $G_{M}^{p}(Q^{2})$ they are in agreement except for the two lowest $Q^{2}$ values. Unfortunately, LHPC results carry large errors and in general are compatible both with our values and the experimental ones prohibiting any definite conclusions as to the nature of the discrepancy with the experimental values. As discussed volume and residual excited state effects may lead to a slow convergence of the lattice data that can account for the discrepancies with the experimental values.

Results for the neutron electromagnetic form factors are only provided by the ETMC for pion masses below 170 MeV. They are compared to the experimental values in Fig. 29. We observe that results for the electric form factor extracted from the cB211.072.64 ensemble that includes disconnected contributions are in agreement with the experimental values. This is also true for the cA2.09.48 ensemble that includes disconnected contributions although they carry larger errors. For the cA2.09.64 ensemble, where disconnected contributions have not been included, underestimate the electric neutron form factor. This clearly indicates the significance of including disconnected contributions, especially for this quantity, an observation consistent with the conclusion reached also in Ref. Sufian et al. (2017b). For the magnetic form factor, results using the cB211.072.64 twisted mass ensemble with disconnected contributions are closer to experiment as to compared to the $N_{f}=2$ ensembles, but there is still a discrepancy with the experiment for small $Q^{2}$ values that needs to be further investigated.

In Fig. 30, we compare the lattice QCD values of the isovector r.m.s radii $\sqrt{\langle r_{E}^{2}\rangle^{u-d}}$, and $\sqrt{\langle r_{M}^{2}\rangle^{u-d}}$ finding agreement among them. As expected by the less steep fall-off of the electric isovector form factor, lattice QCD results are systematically lower than the experimental values. We note that the ETMC results have errors that are already the same as the difference between the two experimental determinations showing that the statistical accuracy required can be achieved. A high-statistics dedicated study to better assess the remaining systematics can thus yield valuable insights on the r.m.s. charge radius from a first principles calculation. In the case of $\langle r_{M}^{2}\rangle^{u-d}$ the errors are larger and lattice QCD results are both in good agreement among them and compatible with the PDG value Patrignani et al. (2016).

In Fig. 31 we show the corresponding quantities for the proton. Only the ETMC results include disconnected contributions, which, although small, have a systematic effect. We observe a similar behavior as for the isovector case, namely smaller values for the electric and magnetic r.m.s radii. LHPC results extracted using the summation method have larger errors and are thus compatible with both the muonic and electron scattering determinations of the r.m.s. radii. For the neutron radii we have only results from ETMC and LHPC. They are displayed in Fig. 32. ETMC results on the electric r.m.s. radius are determined at high accuracy and include all contributions. Although they are still smaller in magnitude than the experimental values, the discrepancy is within one standard deviation. We note that including disconnected contributions brings better agreement in particular in the case of $\langle r_{E}^{2}\rangle^{n}$.

VIII Proton and neutron electromagnetic form factors

Having compared with other groups and with the $N_{f}$=2 results from ETMC, we collect here our final results on the proton and neutron form factors using the $N_{f}$=2+1+1 ensemble, which has the most accurate results at the physical point. In Fig. 33 we show our results for the proton electric and magnetic form factors compared to experimental data. As expected from the behavior observed for the isovector and isoscalar electric form factors, the proton electric form factor is consistently higher than the experimental results. The proton magnetic form factor agrees with the experiment for all $Q^{2}$ except the lowest two. This may be due to finite volume or residual excited state effects as discussed in Section V.

In Fig. 34 we show our results for the neutron form factors. The determination of $G_{E}^{n}(Q^{2})$ directly from lattice QCD is very promising: we find good agreement with the experimental values but more importantly, at low $Q^{2}$, the errors from lattice QCD are smaller by up to a factor of four in some cases, allowing for a more precise description of its $Q^{2}$ dependence. The lattice QCD determination yields also accurate results for $G_{M}^{n}(Q^{2})$ that are in agreement with experiment for $Q^{2}{>}0.2$ GeV2. At small $Q^{2}$ we observe the same discrepancy as that observed for the isovector case. Such an underestimation have been seen also for the induced pseudo-scalar form factor where leading order chiral perturbation theory can show that is due to multi-hadron state contributions with pions. Whether this is the explanation also for the neutron magnetic form factor remains an open question.

Our results for the proton radii and magnetic moment, as extracted from the dipole fit, are

[TABLE]

The corresponding quantities for the neutron using the Galster-like parameterization for the electric and the dipole form for the magnetic are

[TABLE]

As already explained, the first error is statistical, the second is an estimate of the systematic due to the Ansätz chosen for the fit and the third an estimate of excited state effects. We note here that disconnected contributions to $\langle r^{2}_{E}\rangle^{n}$ are non-negligible. If we were to neglect them we would obtain $\langle r^{2}_{E}\rangle^{n,conn.}{=}{-}0.063(15)$ fm2, namely more than a 15% shift in the mean value, i.e. comparable to the other quoted systematic errors.

IX Summary and conclusions

The nucleon electromagnetic Sachs form factors are computed using an $N_{f}{=}2{+}1{+}1$ ensemble of maximally twisted mass fermions with quark masses tuned to their physical values as well as an ensemble of $N_{f}{=}2$ twisted mass fermions simulated at a pion mass of 130 MeV. Comparing results calculated using $N_{f}$=2 and $N_{f}$=2+1+1 twisted mass ensembles leads to the conclusion that no quenching effects are detected within the accuracy of the results that is within a couple of a percentage.

A main novelty of this work is the computation to an unprecedented accuracy of the disconnected light quark contributions, allowing us to extract the individual proton and neutron electromagnetic form factors. This is accomplished by using state-of-the-art techniques that combine hierarchical probing and deflation of the lowest eigen-modes and a large number of randomly distributed smeared point sources in order to suppress gauge noise. In particular, we find that disconnected contributions to the neutron electric form factor are non-negligible and need to be taken into account to bring agreement with the experimental values.

Excited states are thoroughly investigated using five sink-source time separations in the range of [0.96-1.60] fm allowing the identification of the ground state to good precision and the determination of a systematic error due to the excited states by comparing results from the plateau method with the two-state fit method. The summation method is used as a confirmation of the results extracted from the plateau and two-states fits.

Our values for the electric and magnetic r.m.s. radii as well as the magnetic moments for the isovector, isoscalar, proton and neutron are collected in Table 5.

The results are extracted using the dipole ansätz or the Galster-like parameterization and a systematic error on the parameterization used is extracted by comparing with the model independent z-expansion. Our result for the proton electric r.m.s radius is underestimated due to the slower decay of $G_{E}^{p}(Q^{2})$. Similarly there is an underestimation of the magnetic moments for the proton and neutron. A most plausible explanation for these remaining discrepancies may come from a combination of residual volume and multi-hadron contributions. Finite volume effects are investigated in this work by comparing two $N_{f}{=}2$ twisted mass ensembles with pion mass of 130 MeV with the same lattice spacing but $Lm_{\pi}{\simeq}3$ and $Lm_{\pi}{\simeq}4$. Although we observe consistent results between these two volumes, we cannot exclude finite volume effects that may affect the magnetic form factor for small $Q^{2}$ values as well as the electric form factor. A slow convergence of the results as a function of the volume in combination with effects of multi-hadron states maybe difficult to detect. A theoretical investigation within chiral perturbation theory can shed light on multi-hadronic contributions. Furthermore, a study on a larger volume will also help to probe adequately volume effects. Thus, further studies are required to be able to take the infinite volume limit and make definite conclusions on the small $Q^{2}$ behavior of the magnetic form factor and on the slope of the electric form factor. Finite lattice spacing effects, although are expected to be small, need to also be investigated. Before this program is completed one cannot make final statements on the two experimental results for the proton charge radius. The ETM collaboration is generating further ensembles in order to enable the investigation of these issues that will require large computational resources.

Acknowledgements.

We would like to thank all members of ETMC for a very constructive and enjoyable collaboration. M.C. acknowledges financial support by the U.S. National Science Foundation under Grant No. PHY-1714407. This project has received funding from the Horizon 2020 research and innovation program of the European Commission under the Marie Skłodowska-Curie grant agreement No 642069. S.B. is supported by this program as well as from the project COMPLEMENTARY/0916/0015 funded by the Cyprus Research Promotion Foundation. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding the project pr74yo by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). Results were obtained using Piz Daint at Centro Svizzero di Calcolo Scientifico (CSCS), via the project with id s702. We thank the staff of CSCS for access to the computational resources and for their constant support. This work also used computational resources from Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number TG-PHY170022. This work used computational resources from the John von Neumann-Institute for Computing on the Jureca system at the research center in Jülich, under the project with id ECY00.

Appendix A Expressions relating nucleon vector matrix elements to electromagnetic form factors

In this Appendix we give a summary of the expressions relating the Sachs form factors $G_{E}{\equiv}G_{E}(Q^{2})$ and $G_{M}{\equiv}G_{M}(Q^{2})$ to the ratio of three-point and two-point functions. The expressions are given for a general frame with initial (final) momentum ${\vec{p}}$ (${\vec{p}}^{\prime}$) and initial (final) energy $E$ ($E^{\prime}$). All expressions are given in Euclidean space.

[TABLE]

[TABLE]

where $C$ is a kinematic factor given by

[TABLE]

In the case where $\vec{p}\,^{\prime}=\vec{0}$ the expressions simplify as follows

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Appendix B Numerical results for the electromagnetic form factors