Proton and neutron electromagnetic form factors from lattice QCD

The electromagnetic form factors of the proton and the neutron are computed within lattice QCD using simulations of maximally twisted mass fermions with masses fixed to their physical values. We analyse two new ensembles of $N_f=2$ and $N_f=2+1+1$ and 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. The light quark disconnected contributions are evaluated to very high accuracy enabling us to include them in the determination of the proton and neutron form factors, the electric and magnetic radii, and the magnetic moments. We find that disconnected contributions are non-negligible and contribute up to 15% for the case of the neutron electric charge radius.


I. INTRODUCTION
Nucleons, being composite particles, have a non-trivial internal structure that can be probed by measuring their electromagnetic form factors. Despite having been extensively studied both theoretically and experimentally, there is still interest in obtaining these quantities at higher precision and over a wider range of momentum transfers. The proton electric form factor is extracted to high precision from electron proton scattering [1]. 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. [2] for a recent review). A pioneering experiment using Lamb shifts in muonic hydrogen surprisingly found a value smaller by five standard deviations [3], 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 [4], which is accessed indirectly experimentally through electrondeuteron or electron-helium scattering and therefore remains poorly-known. It is of substantial importance to compute these fundamental quantities from first principles, and lattice QCD, with simulations at physical values of the QCD parameters, provides an ideal formulation for such an investigation.
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. [5]. 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 [6,7]. Comparison between N f = 2 and N f = 2 + 1 + 1 also sheds light on any possible unquenching effects 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 [8] and the model independent z-expansion [9]. 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 choice of the fit method.
A crucial component of our analysis is the use of hierarchical probing [10] combined with deflation of the lower lying eigenvalues [11] 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. We include an assessment of finite volume and unquenching effects using results from two N f =2 ensembles. In Sec-tion V, 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 VI [12][13][14][15]. Our final results for the proton and neutron electromagnetic form factors are given in Section VII. Finally, in Section VIII, we summarize our findings and conclude. For completeness, we include in Appendix A the decomposition of the nucleon matrix elements in terms of the form factors and in Appendix B a table with the numerical results for the electric and magnetic form factors as a function of the square of the momentum transfer.

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, N (p, s) is the nucleon state with initial (final) momentum p (p ) and spin s (s ), with energy E N ( p) (E N ( p )) and mass m N . q 2 ≡q µ q µ is the momentum transfer squared q µ =(p µ − p µ ), u N is the nucleon spinor and j µ is the vector current. The local vector current is given by, where q f and e f is the quark field of flavor f and its electric charge respectively and the summation runs over all the quark flavors. We use the symmetrized lattice conserved vector current given by 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, G M (q 2 ) = F 1 (q 2 ) + F 2 (q 2 ) .
In the isospin limit, where the up and down quarks are degenerate, we consider the isovector combination p|j u µ − j d µ |p that gives the difference between the proton and neutron form factors and the isoscalar combination p|j u µ + j d µ |p /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 p E (0)=1 and G n E (0)=0 which, when using the lattice conserved current, holds by symmetry, even prior 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 → 0, namely:

III. CALCULATION ON THE LATTICE
A. 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 where u and d are up-and down-quark spinors and C=γ 0 γ 2 is the charge conjugation matrix. The two-point function in momentum space is given by and the three-point function is given by The initial position and time (x 0 ) is referred to as the source, the position and time of the current j µ is (x ins ) as the insertion and the final position (x s ) as the sink. Γ ν is a projector acting on spin indices, with Γ 0 = 1 2 (1+γ 0 ) yielding the unpolarized and Γ k =Γ 0 iγ 5 γ 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. To increase the overlap of the interpolating field with the nucleon ground state, so that excited states are suppressed at earlier time separations, we use Gaussian smeared point sources [16,17] with 125 steps and parameter α G =0.2 and we smear the gauge links which enter the Gaussian smearing operator with APE [18] smearing, using 50 steps and α AP E =0.5. An optimized ratio [19][20][21] of the three-point function over a combination of two-point functions is used to cancel time dependent exponentials and overlaps, given by where t s and t ins are taken to be relative to the source t 0 for simplicity. In the limit of large time separations, (t s −t ins ) 1 and t ins 1, the lowest state dominates and the ratio becomes time independent G E (Q 2 ) and G M (Q 2 ) are extracted from linear combinations of Π µ (Γ ν ; p , p) as expressed in Appendix A, with Q 2 ≡ − 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 p 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, setting p = 0. We use an appropriately tuned multigrid algorithm [22][23][24] 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 where D −1 (x ins ; x ins ) is the quark propagator which starts and ends at the same point x ins and G is a general γ-structure, which for this work is the vector current. A direct computation of quark loops would need inversions from all-to-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 [25] that take into account the sparsity of the Dirac operator and its decay properties. Namely, in this work, we employ the hierarchical probing technique [10], which provides a partitioning scheme that eliminates contributions from neighboring points in the trace of Eq. 11 up to a certain 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 increases by one. 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/ √ N r , i.e.
where N r is the size of the stochastic ensemble. Hierarchical probing has been employed with great success in previous studies [26,27] 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 [11], 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 with the contribution from the remaining modes which are estimated using hierarchical probing. Additionally, we employ the one-end trick [28], used in our previous studies [29][30][31][32] and fully dilute in spin and color.

B. Gauge Ensembles and Statistics
For the extraction of the electromagnetic form factors we analyze one N f =2+1+1 ensemble [5] 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 O(a) improvement is achieved by tuning to maximal twist [33,34]. The N f = 2 + 1 + 1 ensemble is simulated using a lattice of 64 3 ×128, lattice spacing a = 0.0809(2)(4) fm, pion mass m π = 138(1) MeV, and Lm π =3.62 [5]. To assess finite volume effects, we use two N f =2 ensembles, which only differ in their volume, namely one has Lm π = 2.98 and the other Lm π =3.97 while both have pion mass 130(1) MeV and a=0.0938(3)(1) fm. Results for the former ensemble are from Refs. [5,7] while results for the latter are reported here for the first time. The simulation parameters are tabulated in Table I. For the analysis of the N f =2+1+1 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.97 fm to 1.62 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 II we give the statistics used in the calculation of the connected three-point functions. For the evaluation of the disconnected contributions we use N 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 men-tioned 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 spincolor dilution. For the hierarchical probing we use distance eight coloring resulting in 512 Hadamard vectors, which 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. Table III summarizes the parameters for the computation of the disconnected three-point functions. The N f =2, 64 3 × 128 ensemble is used to check for finite volume effects, comparing the connected contributions to those of the N f =2, 48 3 × 96 ensemble. For the latter, the setup is reported in Ref. [7] and summarized in Table I. For the former, larger lattice, 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 analyze more configurations for higher separations to control statistical error. In Table IV we summarize the statistics. Assessment of excited state effects is imperative for the proper extraction of the desirable nucleon matrix element. However, ensuring ground state dominance is a delicate process due to the exponentially increasing statistical noise with the 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. (9) 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 twoand three-point functions considering contributions up to the first excited state using the expressions In Eqs. (13) -(14) E 0 ( p) and E 1 ( p) are the energies of the ground and first excited state with total momentum p, respectively. Fitting simultaneously the two-and three-point functions involves twelve parameters. Note that for non-zero momentum transfer, A µ 01 (Γ ν , p , p) = A µ 10 (Γ ν , p , p). The desired nucleon matrix element is obtained from Summation method: Summing over t ins in the ratio of Eq. (9) yields a geometric sum [36,37] from which we obtain, where the ground state contribution, Π µ (Γ ν ; p , 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 excited states.
Derivative Summation method: Instead of performing a linear fit in Eq. (16) to extract the matrix element, one can take finite differences to the summed ratio [38] as follows and fit to a constant to extract the desired matrix element.

D. Momentum dependence fits of form factors
Electromagnetic form factors are described by the nucleon matrix element of the vector current, and 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 [39] 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. (18) and Eq. (5) one can relate M to the mean square radius as The neutron electric form factor and disconnected contributions to the electric form factors are zero for Q 2 =0 and we treat them as a special cases, fitting them using the Galster-like parameterization [8,40], given by with A, B fit parameters and m N the nucleon mass. In this case the radius is given by Another fit form, which has been applied recently to experimental data of both electromagnetic and axial form factors, is the model independent z-expansion [9]. In this case, the form factor is expanded in a series given by, where and t cut is the time-like cut of the form factor. We take t cut =4m 2 π for the isovector combination G u − G d and t cut =9m 2 π for the isoscalar combination G u + G d [9]. For convergence of the truncated series of Eq. (22) the coefficients a k should be bounded in size and convergence should be demonstrated by increasing k max . The interested reader is referred to Ref. [6] for details about our procedure. The mean square radius is given by while the value of the form factor at zero momentum transfer is G(0)=a 0 .

A. 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. For the connected diagram we use a frame where the nucleon final momentum is zero, thus q=− p.
In Figs. 2 and 3 we show the ratios defined in Eq. (9) as a function of the sink-source time separation, and for three values of the momentum transfer squared, that is for Q 2 =0.056 GeV 2 , Q 2 =0.216 GeV 2 and Q 2 =0.547 GeV 2 . In a frame where the nucleon final momentum is zero, the expression given in Eq. (A1) and Eq. (A2) reduce to Eqs. (A4,A5), and (A6) in Appendix A, giving separately the electric and magnetic form factors. In the case of G u−d E (Q 2 ) the plateau value decreases as t s increases with the effect becoming larger as Q 2 increases showing 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 is observed. we show linear fits to the summed ratio for three different values of Q 2 . The slope of the lines 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 flat within errors, 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. (9) and the curves are ob-tained by fitting simultaneously the three-and two-point functions to Eqs. (14) and (13). The gray horizontal band shows the nucleon matrix element value and error extracted from the two-state fit as in Eq. (15). 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.  16) as a function of the sink-source time separation for three momenta, namely n 2 q =1 (blue circles), n 2 q =4 (green squares) and n 2 q =10 (red triangles) corresponding to Q 2 =0.056, 0.216, and 0.547 GeV 2 for the isovector electric (top) and isovector magnetic (bottom) form factors. The bands are fits to a linear form. Right panel: the derivative of the summed ratio as in Eq. (17), using the same notation as that of the left panel. The bands are fits to a constant.
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 2values 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 u−d E (Q 2 ) and G u−d M (Q 2 ) as a function of the momentum transfer squared Q 2 . We show up to Q 2 =0.5 GeV 2 to facilitate comparing the values extracted using the plateau at the four largest separations, the summation, and the twostate fit approaches. The summation and two-state fit are obtained using the fit ranges indicated by the open symbols in Fig. 6. As can be seen, the results follow the observation made in Fig. 2, namely that 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. For G u−d M (Q 2 ), excited state effects are small for larger Q 2 values whereas for smaller Q 2 values there is an 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 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. In what follows, we will consider the plateau values for the largest time separation as the final results on the form factors. We will use the deviation from the values determined from the two-states fits as an estimate of the systematic error due to excited states.

B. Assessment of finite volume and unquenching effects
For the assessment of finite volume effects we compare two N f =2 physical point ensembles with m π L 3 [7] and m π L 4. The lattice spacing and pion mass are the same for the these two ensembles. The isovector electric and magnetic form factors are shown in Fig. 9, where results are extracted using the plateau method. The results fall, within errors, on the same curve indicating no significant finite volume effects between the two volumes available to us, namely m π L 3 and m π L 4. The same behavior is observed for the isoscalar form factors shown in Fig. 10.
To check for unquenching effects, we compare the N f =2 ensemble with m π L 3 to the N f =2+1+1 ensemble. From this comparison, we see that both isovector, shown in Fig. 11, and isoscalar, shown in Fig. 12, show no sensitivity to quenching effects of the strange and charm quarks, at least to the accuracy of our data.

C. 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 B with the statistics summarized in Table III. 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 nonperturbatively, in the RI -MOM scheme, employing momentum sources. We perform a perturbative subtraction of O(g 2 a ∞ )-terms, as described in Refs. [41,42], which subtracts the leading cut-off effects leaving only a weak dependence on the renormalization scale (aµ) 2 , as shown in Fig. 13. We find a value of Z V =0.728(1) where the error is statistical. Alternatively, Z V can be determined at 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 is 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 timeslice allowing us to compute the three-point function for every combination of t s and t 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 p = 0, for the disconnected diagrams we have all sink momenta at no additional cost. We analyze, besides p = 0, the matrix element for the six final momenta with p = ± 2π Ln , witĥ n =x,ŷ, orẑ, 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. Two-state fits yield too large errors and are in general consistent with the plateau extraction.
In Fig. 14 we present our results for the disconnected contributions to G u+d E (Q 2 ) and G u+d M (Q 2 ) up to Q 2 =1 GeV 2 for three time separations in the range 0.97 fm to 1.29 fm. As can be seen we can achieve a relative statistical error that is less than 20% for up to t s =14a∼1.15 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 u+d E (Q 2 ) and G u+d M (Q 2 ), that the absolute value tends to increase. 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 much larger errors. We therefore take as our final result for the disconnected contribution the value extracted using t s /a=14 for both G u+d E (Q 2 ) and G u+d M (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 ef- , ts/a=18 (squares) and ts/a=20 (circles), compared to the summation method (diamonds) and using two-state fits (stars). Results from different methods are slightly shifted to the right for clarity.
fects when we quote quantities that include disconnected contributions.
In Fig. 15 we show the disconnected contributions to G u+d E and G u+d M comparing between results obtained using the N f = 2 + 1 + 1 and N f = 2 with Lm π ∼ 3 twisted mass ensembles. As can be seen, the values are consistent but the N f =2+1+1 results have up to four times smaller statistical errors. In evaluation of the disconnected contributions of the N f = 2 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 [6]. This is approximately the same number of inversions (and thus cost) as for the N f =2+1+1 ensemble (see Table III), which demonstrates the effectiveness of the hierarchical probing method employed in the current analysis of the N f =2+1+1 ensemble [6].

V. Q 2 -DEPENDENCE OF THE ISOVECTOR AND ISOSCALAR FORM FACTORS
In this section we fit the Q 2 -dependence of form factors to the forms introduced in Section III D. We consider first the isovector form factors where only the connected diagram contributes. In Figs. 16 and 17 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 u−d E (Q 2 ) and G u−d M (Q 2 ), respectively. As can be seen, the values extracted from the plateau and twostate fits are consistent and do not show any 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 twostate fits to estimate the systematic error due to excited states.
In Figs G u−d M (Q 2 ), respectively using the dipole form and the z-expansion of Eq. (22) and comparing to experiment. For the z-expansion, we check convergence by increasing k max and observing the resulting magnetic moment and r.m.s radii, shown in Fig. 20, where we observe convergence for k max =4. In the case of G u−d E (Q 2 ), we see from Fig. 18 that, while our results are compatible with the experimental results for Q 2 < 0.1 GeV 2 , in the range 0.1-0.6 GeV 2 they are consistently higher. Therefore, although both dipole form and z-expansion describe very well our data, shown in separate panels for clarity, they lie about one to two standard deviations above the experimental values. In extracting the r.m.s radius, we see from Fig. 20 that results obtained from using the dipole fit and z-expansion are compatible, and yield 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 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. 25.
For G u−d M , shown in Fig. 19, we observe that our results are in agreement with the experimental values for Q 2 > 0.15 GeV 2 , 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 [57], 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 π 3 to Lm π 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 π 7.4 [58] finds a higher value that may (green triangles) and N f =2+1+1 (red circles). Results are extracted using the plateau method for sink-source separation ts 1.6 fm for the N f =2+1+1 and for the electric for the N f =2 while for the magnetic N f =2 the separation is ts 1.3 fm. 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. 20. As can be seen, the mean value extracted for µ u−d using the dipole and z-expansion is the same, while for r 2 M u−d the z-expansion produces a slightly higher value, which however is consistent within errors. Using the dipole fit we find Here we have included a fourth systematic error computed as the difference in the values of µ u−d and r 2 M u−d when fitting G u−d M (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  [7] and N f =2+1+1 wit the same notation as in Fig. 11.
what follows we will not include this fourth systematic error.
Before presenting fits to the total isoscalar form factors we discuss the Q 2 -dependence of the disconnected contributions. In Fig. 21 we show the disconnected contribution to the isoscalar electric form factor G u+d E , 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 u+d E (0)=0 for the disconnected contribution. Although both parameterizations describe well our results, we observe different behavior in the errors, namely the z-expansion has a larger error band for larger Q 2 values.
The disconnected contribution to G u+d M (Q 2 ) is shown in Fig. 22. 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 u+d E (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 The isovector electric form factor as a function of Q 2 (circles). We show fits to our results using a dipole form (top) and using the z-expansion (bottom) for kmax=4. Black crosses are experimental results taken from the A1 collaboration [1] for the proton and from Refs. [4,[43][44][45][46][47][48][49][50][51][52][53][54][55][56] for the neutron.
where we have not normalized with the value of the form factor at zero momentum transfer, i.e. the radii are  rather than from Eq. 5. In Fig. 23 we compare results between including and excluding disconnected contributions. Although the effect is small for both G u+d E (Q 2 ) and G u+d M (Q 2 ) there is a shift in the fits affecting the parameters of the fits. This comparison shows that disconnected contributions are thus non-negligible 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. 24 and 25 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 by one to two standard deviations 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. 26 we show the isoscalar electric square radius. As can be seen, the z-expansion fit yields values that are within errors for k max > 1 but with two times larger errors than the dipole. In Fig. 26 we also show results for the magnetic moment and the magnetic radius where convergence of the z-expansion is observed already for k max =2. For the magnetic moment,  24. Isoscalar electric form factor (circles) as a function of Q 2 . We combine the connected contribution from the plateau for ts/a=20 with the disconnected contribution for ts/a=14. The remaining notation is as in Fig. 18. Results for the isoscalar charge square radius, magnetic moment, and magnetic square radius. The notation is as in Fig. 20.
rors while for the magnetic radius the z-expansion tends to be slightly higher but with bigger errors. In what follows we will quote the values determined from the dipole fits. They are r 2 E u+d = 0.697(9)(10)(0) fm (31) µ u+d = 2.64(12)(13)(6) and (32) Note that from our definition of the isoscalar combination, the proton plus neutron magnetic moment is obtained by: µ p + µ n = µ u+d /3.

VI. 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 work 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 π =149 MeV, lattice spacing a=0.116 fm and Lm π =4.21 at three sink-source separations from 0.93 fm to 1.39 fm [13]. 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 π = 135 MeV, lattice spacing a=0.093 fm and Lm π =4 [14]. 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 π =146 MeV, a 0.085 fm and a spatial extend of 8.1 fm or Lm π 6 allowing access to relatively small momenta [15]. 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. Recently a followup study from PACS appeared but this work has not yet been published [58].
In the top panel of Fig. 27 we show lattice QCD results from the aforementioned analyses for G u−d E (Q 2 ) up to Q 2 =0.5 GeV 2 . As can be seen, ETMC and PACS results, which use similar techniques to identify the nucleon matrix element are in good agreement but systematically higher than the experimental values. LHPC results were obtained using the summation method and thus have larger statistical errors making them compatible with the experimental values. In the lower panel of Fig. 27, we show the corresponding results for G u−d M (Q 2 ). The recent ETMC results 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 GeV 2 . 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 π ∼ 3 to Lm π ∼ 4) the volume dependence could be weak and require a larger volume to manifest itself. The yet unpublished PACS results may indicate such a trend. A conclusion that we can, however, draw from these lattice QCD studies is that the fact that there is agreement among them for both the electric and magnetic form factors when different discretization schemes are employed , the N f =2 twisted mass ensemble with mπL 3 from Ref. [7] (green triangles), LHPC using N f =2+1 stout-smeared clover fermions from Ref. [13] (left purple triangles) and Ref. [14] (right yellow triangles) and from PACS using N f =2+1 stoutsmeared clover fermions from Ref. [15] (cyan rhombuses).
indicates that cut-off effects are smaller than the statistical errors.
The proton and neutron form factors can be extracted from the isovector and isoscalar form factors discussed in Section V, using the linear combinations 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 include disconnected contributions, which have only been computed by ETMC at close to physical pion masses. We use filled symbols to indicate lattice results that include disconnected contributions. For both the electric and magnetic form factors our results and those from LHPC are within error bars, with the LHPC results exhibiting larger errors at higher Q 2 possibly due to the use of the summation method. The accurate ETMC results are higher than the experimental values for G P E 0Q 2 ) and lower for G p M (Q 2 ) for the two lowest Q 2 values. Although LHPC results do not show such discrepancies, unfortunately, they carry large errors and thus do not allow for any conclusions to be drawn on the possible reasons for the discrepancy with experiment.
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 N f =2+1+1 ensemble that include disconnected contributions are in agreement with the experimental values. This is also true for the N F = 2 ensemble with Lm π ∼ 3 that includes the disconnected contributions although they have larger errors. N f =2 results, for which 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. For the magnetic form factor, results using the N f =2+1+1 twisted mass ensemble with disconnected contributions are closer to experiment compared to N f = 2, but there is still a discrepancy with the experiment for small Q 2 that needs to be further investigated.
In Fig. 30 [4,[43][44][45][46][47][48][49][50][51][52][53][54][55][56] for the electric form factor and from Refs. [59][60][61][62][63][64] for the magnetic form factor. 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 by increasing statistics. 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 r 2 M u−d the errors are larger and lattice QCD results are both in good agreement among them and compatible with the PDG value [66].
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 when errors are small. 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 M n and µn. The notation is as in Fig. 31.

VII. PROTON AND NEUTRON ELECTROMAGNETIC FORM FACTORS
Having compared with other groups and with 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 one, which may be due to finite volume effects as discussed in Section V.
In Fig. 34 we show our results for the neutron form factors. The determination of G n E (Q 2 ) directly from lattice QCD is remarkable; there is first of all good agreement with experiment 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 determination of G n M (Q 2 ) is also very accurate yielding agreement with experiment for Q 2 >0.2 GeV 2 while at small Q 2 we observe the same discrepancy as that observed for the isovector case.
Our results for the proton radii and magnetic moment, as extracted from the dipole fit, are as follows The corresponding quantities for the neutron using the Galster-like parameterization for the electric and the 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 r 2 E n are non-negligible. If we were to neglect them we would obtain r 2 E n,conn. =−0.064(15) fm 2 , namely more than a 15% shift in the mean value, i.e. comparable to the other quoted systematic errors.

VIII. 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. The 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  [4,[43][44][45][46][47][48][49][50][51][52][53][54][55][56] for the case of the electric form factor and from Refs. [59][60][61][62][63][64] for the case of the magnetic form factor. The fits to our results use Eq. (20) for the electric form factor and Eq. (18) for the magnetic form factors.
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 separation in the range of [0.97-1.62] fm allowing the identification of the ground state to good precision and the determination of a reliable 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.
Finite volume effects are investigated by comparing two N f =2 twisted mass ensembles with pion mass of 130 MeV with the same lattice spacing but Lm π 3 and Lm π 4. We observe consistent results between these two volumes, but we cannot exclude finite volume effects that may affect the magnetic form factor for small Q 2 values. Further studies are required to take the infinite volume limit and make definite conclusions on the small Q 2 behavior of the magnetic form factor. Comparing results calculated using N f =2 and N f =2+1+1 twisted mass ensembles we observed no quenching effects. 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 V. The results   TABLE V. Our results for the electromagnetic radii and the magnetic moment using the N f =2+1+1 ensemble for the isovector combination (p−n), isoscalar (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.  (16)(12)(11) 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 p E (Q 2 ) and further studies are required in order for lattice QCD results to reach an accuracy that can distinguish between the experimental results for the proton charge radius. Π µ (Γ 0 , p , p) = −iC G E 2m(4m 2 + Q 2 ) (p µ + p µ ) m (E + E + m) − p ρ p ρ + C G M 4m 2 (4m 2 + Q 2 ) δ µ0 4m 4 + m 2 Q 2 + 4m 2 p ρ p ρ + Q 2 p ρ p ρ +2im ε µkρσ p ρ p σ 2m + E + E + Q 2 2m −2m ε µ0ρσ p ρ p σ (p k + p k ) + 2m ε µk0ρ p σ p σ (p ρ − p ρ ) , where C is a kinematic factor given by In the case where p = 0 the expressions simplify as follows and C = 2m 2 E(E + m) . (A7)