Nucleon Electromagnetic Form Factors in the Continuum Limit from 2+1+1-flavor Lattice QCD

TL;DR

This study uses high-statistics lattice QCD simulations to accurately determine nucleon electromagnetic form factors, extrapolating to the continuum limit and analyzing their $Q^2$ dependence with various methods, providing precise radius and magnetic moment estimates.

Contribution

First comprehensive lattice QCD analysis of nucleon electromagnetic form factors across multiple lattice spacings and quark masses with continuum extrapolation.

Findings

Form factors consistent with experimental data within uncertainties.

$z$-expansion and dipole form yield compatible radii and magnetic moments.

Data plotted versus $Q^2/M_N^2$ reduces discretization errors.

Abstract

Results are presented for the nucleon isovector electromagnetic form factors using 11 ensembles generated by the MILC collaboration using the 2+1+1-flavors HISQ action. They span 4 lattice spacings $a \sim$ 0.06, 0.09, 0.12 and 0.15~fm and 3 values of $M_\pi \sim 135, 225$ and 315 MeV. High-statistics estimates are used to perform a simultaneous extrapolation in the lattice spacing, lattice volume and light-quark masses. The $Q^2$ dependence over the range 0.05-1.4 ${\rm GeV}^2$ is investigated using both the $z$-expansion and the dipole form. Final $z$-expansion estimates for the isovector r.m.s. radius are $r_E = 0.769(27)(30)$ fm $r_M = 0.671(48)(76)$ fm and $\mu^{p-n} = 3.939(86)(138)$ Bohr magneton. The first error is the combined uncertainty from the leading-order analysis, and the second is an estimate of the additional uncertainty due to using the leading order chiral-continuum-finite-volume fits. The dipole estimates, $r_E = 0.765(11)(8)$ fm, $r_M = 0.704(21)(29)$ fm and $\mu^{p-n} = 3.975(84)(125)$, are consistent with those from the $z$-expansion but with smaller errors. Our analysis highlights three points. First, all data from the eleven ensembles and existing lattice data on, or close to, physical mass ensembles from other collaborations collapses more clearly onto a single curve when plotted versus $Q^2/M_N^2$ as compared to $Q^2$ with the scale set by quantities other than $M_N$. The difference between these two analyses is indicative of discretization errors, some of which presumably cancel when the data are plotted versus $Q^2/M_N^2$. Second, the size of the remaining deviation of this common curve from the Kelly curve is small and can be accounted for by statistical and possible systematic uncertainties. Third, to improve lattice estimates, high statistics data for $Q^2 < 0.1$ ${\rm GeV}^2$ are needed.

Tables

Table 1. Table 1: The values of the space-like four-momentum squared, Q 2 superscript 𝑄 2 Q^{2} , transferred to the ground state nucleon, in units of GeV 2 . The data for the thirteen calculations defined in Table  12 are labeled by the 3-momentum vector n → → 𝑛 \vec{n} . The ground state energy is obtained from a 4-state fit.

$\overrightarrow{n}$	$a15m310$	$a12m310$	$a12m220L$	$a12m220$	$a12m220S$	$a09m310$	$a09m220$
$(1,0,0)$	0.2519(5)	0.1765(5)	0.0670(1)	0.1047(4)	0.1747(15)	0.1834(3)	0.0861(2)
$(1,1,0)$	0.4831(14)	0.3415(13)	0.1318(2)	0.2060(15)	0.3386(35)	0.3558(12)	0.1685(4)
$(1,1,1)$	0.7034(25)	0.4982(24)	0.1947(4)	0.3012(20)	0.4905(61)	0.5198(43)	0.2479(8)
$(2,0,0)$	0.9111(60)	0.6459(35)	0.2565(8)	0.3909(25)	0.6358(87)	0.6735(44)	0.3244(14)
$(2,1,0)$	1.1020(67)	0.7871(42)	0.3159(10)	0.4824(37)	0.774(10)	0.8186(79)	0.3983(18)
$(2,1,1)$	1.2971(91)	0.9202(52)	0.3740(13)	0.5678(47)	0.910(13)	0.9610(127)	0.4703(23)
$(2,2,0)$	1.6372(215)	1.178(9)	0.4872(21)	0.7321(81)	1.178(23)	1.1974(92)	0.6077(37)
$(2,2,1)$	1.8026(222)	1.293(10)	0.5413(25)	0.8077(103)	1.307(25)	1.3229(131)	0.6743(44)
$(3,0,0)$	1.7896(289)	1.315(19)	0.5412(28)	0.8064(118)	1.238(33)	1.3248(168)	0.6713(46)
$(3,1,0)$	1.9171(314)	1.435(18)	0.5950(32)	0.8845(124)	1.358(36)	1.4210(144)	0.7357(51)
$\overrightarrow{n}$	$a09m130W$	$a06m310$	$a06m310W$	$a06m220$	$a06m220W$	$a06m135$
$(1,0,0)$	0.0492(2)	0.1888(13)	0.1899(6)	0.1101(3)	0.1093(3)	0.0513(2)
$(1,1,0)$	0.0974(5)	0.3648(33)	0.3653(15)	0.2159(11)	0.2132(9)	0.1014(6)
$(1,1,1)$	0.1450(9)	0.5322(70)	0.5277(29)	0.3175(24)	0.3130(19)	0.1510(12)
$(2,0,0)$	0.1913(15)	0.6828(99)	0.6895(48)	0.4142(46)	0.4120(55)	0.1975(15)
$(2,1,0)$	0.2373(18)	0.8457(118)	0.8402(65)	0.5087(57)	0.5045(61)	0.2459(22)
$(2,1,1)$	0.2824(23)					0.2941(32)
$(2,2,0)$	0.3704(33)					0.3866(47)
$(2,2,1)$	0.4108(41)					0.4323(51)
$(3,0,0)$	0.4067(48)					0.4259(60)
$(3,1,0)$	0.4490(50)					0.4703(65)

Table 2. Table 2: Results for the bare G E ​ ( Q 2 ) subscript 𝐺 𝐸 superscript 𝑄 2 G_{E}(Q^{2}) extracted from Re ( V 4 ) Re subscript 𝑉 4 \mathop{\rm Re}(V_{4}) are listed for the thirteen calculations defined in Table  12 . The results are obtained using 4-state fits to the 2-point functions and 3 ∗ superscript 3 ∗ 3^{\ast} -state fits to the 3-point functions (2-state fits for the a ​ 12 ​ m ​ 220 ​ S 𝑎 12 𝑚 220 𝑆 a12m220S ensemble) as described in the text. The value, G E ​ ( 0 ) = 1 / Z V subscript 𝐺 𝐸 0 1 subscript 𝑍 𝑉 G_{E}(0)=1/Z_{V} , given in the first row provides one estimate of the renormalization constant for the vector current. The momentum transfer Q 2 superscript 𝑄 2 Q^{2} , in units of GeV 2 , associated with each n → → 𝑛 \vec{n} is given in Table  1 .

$\overrightarrow{n}$	$a15m310$	$a12m310$	$a12m220L$	$a12m220$	$a12m220S$	$a09m310$	$a09m220$
$(0,0,0)$	1.069(4)	1.061(8)	1.067(4)	1.071(9)	1.081(18)	1.045(3)	1.049(4)
$(1,0,0)$	0.650(4)	0.728(8)	0.908(12)	0.840(11)	0.706(17)	0.735(4)	0.859(5)
$(1,1,0)$	0.440(4)	0.536(9)	0.789(12)	0.666(23)	0.513(17)	0.549(6)	0.718(7)
$(1,1,1)$	0.321(4)	0.407(10)	0.694(11)	0.553(16)	0.402(20)	0.423(12)	0.614(8)
$(2,0,0)$	0.261(8)	0.332(11)	0.618(12)	0.469(15)	0.324(20)	0.348(7)	0.538(9)
$(2,1,0)$	0.212(5)	0.279(8)	0.553(10)	0.396(17)	0.279(18)	0.285(9)	0.472(8)
$(2,1,1)$	0.167(6)	0.239(9)	0.499(9)	0.349(15)	0.234(17)	0.240(9)	0.417(8)
$(2,2,0)$	0.140(15)	0.176(15)	0.413(8)	0.280(19)	0.156(24)	0.186(4)	0.338(8)
$(2,2,1)$	0.114(12)	0.161(12)	0.380(7)	0.260(16)	0.155(21)	0.162(4)	0.307(7)
$(3,0,0)$	0.110(30)	0.157(35)	0.387(8)	0.203(32)	0.148(42)	0.177(8)	0.315(9)
$(3,1,0)$	0.088(20)	0.155(24)	0.357(7)	0.200(23)	0.154(27)	0.153(5)	0.290(7)
$\overrightarrow{n}$	$a09m130W$	$a06m310$	$a06m310W$	$a06m220$	$a06m220W$	$a06m135$
$(0,0,0)$	1.052(6)	1.043(6)	1.035(11)	1.050(7)	1.039(9)	1.042(10)
$(1,0,0)$	0.937(6)	0.700(16)	0.711(9)	0.822(8)	0.811(9)	0.919(10)
$(1,1,0)$	0.836(6)	0.502(21)	0.521(8)	0.670(10)	0.654(11)	0.814(13)
$(1,1,1)$	0.756(6)	0.373(24)	0.395(9)	0.552(14)	0.536(14)	0.716(18)
$(2,0,0)$	0.680(8)	0.306(24)	0.318(13)	0.465(17)	0.440(26)	0.664(15)
$(2,1,0)$	0.624(8)	0.232(23)	0.260(11)	0.398(17)	0.384(20)	0.588(18)
$(2,1,1)$	0.571(8)					0.528(20)
$(2,2,0)$	0.497(9)					0.433(21)
$(2,2,1)$	0.455(9)					0.399(19)
$(3,0,0)$	0.439(15)					0.422(21)
$(3,1,0)$	0.418(12)					0.380(20)

Table 3. Table 3: Results for the bare G E ​ ( Q 2 ) subscript 𝐺 𝐸 superscript 𝑄 2 G_{E}(Q^{2}) extracted from Im ( V i ) Im subscript 𝑉 𝑖 \mathop{\rm Im}(V_{i}) are listed for the thirteen calculations defined in Table  12 . The rest is the same as in Table  2 .

$\overrightarrow{n}$	$a15m310$	$a12m310$	$a12m220L$	$a12m220$	$a12m220S$	$a09m310$	$a09m220$
$(1,0,0)$	0.610(15)	0.818(50)	0.871(30)	0.761(61)	0.774(59)	0.699(17)	0.814(30)
$(1,1,0)$	0.435(10)	0.592(32)	0.786(27)	0.634(69)	0.564(38)	0.536(16)	0.694(26)
$(1,1,1)$	0.336(10)	0.448(27)	0.711(28)	0.575(42)	0.439(36)	0.413(29)	0.619(25)
$(2,0,0)$	0.262(21)	0.410(32)	0.654(28)	0.527(42)	0.432(39)	0.358(14)	0.551(27)
$(2,1,0)$	0.205(13)	0.340(24)	0.585(24)	0.429(36)	0.316(32)	0.296(17)	0.479(22)
$(2,1,1)$	0.158(14)	0.284(24)	0.543(23)	0.334(37)	0.249(37)	0.255(14)	0.422(21)
$(2,2,0)$	0.135(10)	0.193(27)	0.464(22)	0.284(37)	0.181(33)	0.206(8)	0.362(18)
$(2,2,1)$	0.103(25)	0.166(30)	0.429(21)	0.250(43)	0.161(45)	0.182(8)	0.308(19)
$(3,0,0)$	0.139(43)	0.188(66)	0.445(25)	0.249(59)	0.202(59)	0.195(15)	0.342(17)
$(3,1,0)$	0.083(58)	0.178(76)	0.405(25)	0.150(74)	0.220(92)	0.172(17)	0.330(25)
$\overrightarrow{n}$	$a09m130W$	$a06m310$	$a06m310W$	$a06m220$	$a06m220W$	$a06m135$
$(1,0,0)$	0.871(43)	0.733(59)	0.641(44)	0.718(46)	0.778(65)	0.793(65)
$(1,1,0)$	0.791(34)	0.515(49)	0.500(29)	0.600(40)	0.660(49)	0.688(57)
$(1,1,1)$	0.710(31)	0.367(64)	0.400(30)	0.534(43)	0.553(51)	0.626(64)
$(2,0,0)$	0.679(31)	0.390(48)	0.270(42)	0.444(53)	0.413(88)	0.643(49)
$(2,1,0)$	0.625(27)	0.219(34)	0.229(37)	0.420(42)	0.366(67)	0.529(52)
$(2,1,1)$	0.582(27)					0.446(55)
$(2,2,0)$	0.494(26)					0.399(54)
$(2,2,1)$	0.461(26)					0.354(54)
$(3,0,0)$	0.467(38)					0.367(60)
$(3,1,0)$	0.435(33)					0.311(58)

Table 4. Table 4: Results for the bare magnetic form factor G M ​ ( Q 2 ) subscript 𝐺 𝑀 superscript 𝑄 2 {G}_{M}(Q^{2}) for the thirteen calculations defined in Table  12 . Values of G M ​ ( 0 ) subscript 𝐺 𝑀 0 G_{M}(0) are obtained by a linear extrapolation of the data for G M ​ ( Q 2 ) / ( G E ​ ( Q 2 ) × Z V ) subscript 𝐺 𝑀 superscript 𝑄 2 subscript 𝐺 𝐸 superscript 𝑄 2 subscript 𝑍 𝑉 G_{M}(Q^{2})/(G_{E}(Q^{2})\times Z_{V}) to Q 2 = 0 superscript 𝑄 2 0 Q^{2}=0 as discussed in the text. The rest is the same as in Table  2 .

$\overrightarrow{n}$	$a12m310$	$a12m310$	$a12m220L$	$a12m220$	$a12m220S$	$a09m310$	$a09m220$
$(0,0,0)$	4.596(61)	4.553(107)	4.538(107)	4.465(144)	4.597(217)	4.324(32)	4.505(76)
$(1,0,0)$	2.968(29)	3.318(51)	4.018(61)	3.657(90)	3.139(85)	3.207(19)	3.749(51)
$(1,1,0)$	2.160(32)	2.597(43)	3.557(41)	3.082(76)	2.352(75)	2.513(24)	3.249(40)
$(1,1,1)$	1.665(26)	2.092(46)	3.172(32)	2.660(68)	1.940(82)	2.041(42)	2.863(36)
$(2,0,0)$	1.255(45)	1.728(59)	2.874(31)	2.251(75)	1.566(98)	1.684(32)	2.468(46)
$(2,1,0)$	1.155(28)	1.532(41)	2.615(30)	2.009(66)	1.389(74)	1.471(35)	2.225(38)
$(2,1,1)$	0.959(31)	1.363(31)	2.407(29)	1.818(67)	1.257(68)	1.283(41)	2.034(36)
$(2,2,0)$	0.817(50)	1.106(54)	2.032(37)	1.643(76)	1.164(92)	1.074(18)	1.698(35)
$(2,2,1)$	0.760(47)	0.969(55)	1.913(32)	1.448(75)	0.963(82)	0.982(19)	1.554(36)
$(3,0,0)$	0.717(102)	1.241(128)	1.883(41)	1.422(88)	0.851(136)	0.947(35)	1.585(44)
$(3,1,0)$	0.734(33)	0.911(104)	1.771(41)	1.424(89)	0.917(111)	0.913(49)	1.467(36)
$\overrightarrow{n}$	$a09m130W$	$a06m310$	$a06m310W$	$a06m220$	$a06m220W$	$a06m135$
$(0,0,0)$	4.297(82)	4.163(168)	4.303(134)	4.138(102)	4.293(142)	4.229(123)
$(1,0,0)$	3.956(67)	3.083(73)	3.181(70)	3.405(65)	3.505(99)	3.824(105)
$(1,1,0)$	3.547(50)	2.440(59)	2.491(57)	2.865(55)	2.873(83)	3.413(90)
$(1,1,1)$	3.281(47)	1.984(73)	2.020(60)	2.493(59)	2.411(89)	3.051(97)
$(2,0,0)$	2.992(48)	1.591(83)	1.655(69)	2.214(68)	2.196(117)	2.838(89)
$(2,1,0)$	2.820(39)	1.341(83)	1.459(58)	1.862(66)	1.837(104)	2.612(84)
$(2,1,1)$	2.616(37)					2.340(97)
$(2,2,0)$	2.286(39)					2.064(96)
$(2,2,1)$	2.156(39)					1.850(97)
$(3,0,0)$	2.158(59)					1.849(110)
$(3,1,0)$	2.041(46)					1.764(96)

Table 5. Table 5: The isovector electric mean-square charge radius ⟨ r E 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝐸 2 \mathop{\langle r_{E}^{2}\rangle}\nolimits in units of fm 2 from the seven fits (dipole, z 3 superscript 𝑧 3 z^{3} , z 3 + 4 superscript 𝑧 3 4 z^{3+4} , z 4 superscript 𝑧 4 z^{4} , z 4 + 4 superscript 𝑧 4 4 z^{4+4} , z 5 superscript 𝑧 5 z^{5} and z 5 + 4 superscript 𝑧 5 4 z^{5+4} ) to the isovector form factor G E ​ ( Q 2 ) subscript 𝐺 𝐸 superscript 𝑄 2 G_{E}(Q^{2}) . The bottom half of the table gives the results of the four CCFV fits discussed in the text, with and without the leading finite volume term c 4 E superscript subscript 𝑐 4 𝐸 c_{4}^{E} and for each Q 2 superscript 𝑄 2 Q^{2} fit.

Ensemble	Dipole	$z^3$	$z^4$	$z^5$	$z^{3+4}$	$z^{4+4}$	$z^{5+4}$
$a15m310$	0.535(6)	0.523(10)	0.519(10)	0.519(10)	0.492(17)	0.531(8)	0.514(12)
$a12m310$	0.561(17)	0.542(23)	0.542(23)	0.542(23)	0.513(31)	0.573(22)	0.528(26)
$a12m220L$	0.575(10)	0.562(32)	0.575(33)	0.575(33)	0.574(45)	0.588(32)	0.562(35)
$a12m220$	0.596(23)	0.557(40)	0.558(40)	0.558(40)	0.546(59)	0.572(35)	0.548(44)
$a12m220S$	0.609(30)	0.686(69)	0.686(67)	0.686(67)	0.688(93)	0.690(57)	0.682(74)
$a09m310$	0.487(6)	0.485(8)	0.485(8)	0.485(8)	0.480(12)	0.494(10)	0.480(10)
$a09m220$	0.580(14)	0.575(26)	0.574(28)	0.573(28)	0.566(31)	0.576(29)	0.568(28)
$a09m130W$	0.587(15)	0.503(51)	0.507(67)	0.506(67)	0.577(85)	0.421(115)	0.546(72)
$a06m310$	0.548(34)	0.537(32)	0.537(32)	0.537(32)	0.533(41)	0.542(36)	0.534(33)
$a06m310W$	0.532(14)	0.502(21)	0.502(21)	0.502(21)	0.483(36)	0.522(20)	0.493(25)
$a06m220$	0.538(22)	0.560(40)	0.561(40)	0.561(40)	0.604(57)	0.543(35)	0.567(44)
$a06m220W$	0.565(22)	0.546(35)	0.546(35)	0.546(35)	0.548(54)	0.551(32)	0.542(39)
$a06m135$	0.599(25)	0.529(64)	0.545(83)	0.545(82)	0.735(135)	0.398(159)	0.649(101)
13-pt $c_4^E\ne 0$	0.592(17)	0.570(39)	0.596(40)	0.595(40)	0.658(55)	0.604(39)	0.601(43)
13-pt $c_4^E=0$	0.581(13)	0.565(30)	0.597(34)	0.597(34)	0.674(47)	0.609(36)	0.618(37)
11-pt $c_4^E\ne 0$	0.586(17)	0.564(39)	0.591(41)	0.590(41)	0.653(56)	0.604(41)	0.597(44)
11-pt $c_4^E=0$	0.572(14)	0.554(32)	0.588(36)	0.587(36)	0.665(49)	0.606(39)	0.609(39)
10-pt $c_4^E\ne 0$	0.587(36)	0.552(69)	0.567(70)	0.567(70)	0.632(93)	0.530(68)	0.592(76)
10-pt $c_4^E=0$	0.558(24)	0.540(48)	0.570(53)	0.570(53)	0.662(74)	0.554(58)	0.618(59)
10∗-pt $c_4^E\ne 0$	0.595(18)	0.587(40)	0.609(42)	0.608(42)	0.668(57)	0.607(41)	0.609(45)
10∗-pt $c_4^E=0$	0.571(14)	0.546(32)	0.577(36)	0.577(36)	0.657(50)	0.587(40)	0.600(40)

Table 6. Table 6: Isovector magnetic charge radius ⟨ r M 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝑀 2 \mathop{\langle r_{M}^{2}\rangle}\nolimits in units of fm 2 from the seven fits to the isovector form factor G M ​ ( Q 2 ) subscript 𝐺 𝑀 superscript 𝑄 2 G_{M}(Q^{2}) . The derived value for G M ​ ( 0 ) subscript 𝐺 𝑀 0 G_{M}(0) is included in the fits as discussed in the text. The rest is the same as in Table  5 .

Ensemble	Dipole	$z^3$	$z^4$	$z^5$	$z^{3+4}$	$z^{4+4}$	$z^{5+4}$
$a15m310$	0.437(11)	0.466(15)	0.464(15)	0.464(15)	0.451(23)	0.472(13)	0.461(17)
$a12m310$	0.414(15)	0.457(30)	0.457(30)	0.457(30)	0.438(41)	0.484(28)	0.447(34)
$a12m220L$	0.456(16)	0.475(41)	0.475(41)	0.475(41)	0.485(58)	0.472(39)	0.473(44)
$a12m220$	0.442(25)	0.419(51)	0.425(50)	0.424(50)	0.406(74)	0.451(46)	0.408(56)
$a12m220S$	0.454(33)	0.597(89)	0.599(87)	0.599(87)	0.566(118)	0.617(77)	0.588(95)
$a09m310$	0.410(8)	0.409(11)	0.409(11)	0.409(11)	0.400(17)	0.423(12)	0.402(14)
$a09m220$	0.469(14)	0.489(30)	0.492(31)	0.492(31)	0.498(38)	0.499(32)	0.488(31)
$a09m130W$	0.478(18)	0.437(56)	0.386(84)	0.384(83)	0.478(132)	0.092(173)	0.485(85)
$a06m310$	0.409(31)	0.408(41)	0.408(41)	0.408(41)	0.436(55)	0.396(46)	0.413(43)
$a06m310W$	0.407(21)	0.406(32)	0.406(32)	0.405(32)	0.388(52)	0.416(33)	0.398(36)
$a06m220$	0.425(24)	0.485(49)	0.484(49)	0.484(49)	0.546(78)	0.438(44)	0.502(55)
$a06m220W$	0.474(35)	0.473(53)	0.473(53)	0.473(53)	0.477(88)	0.465(48)	0.473(59)
$a06m135$	0.519(34)	0.456(76)	0.427(97)	0.427(97)	0.634(214)	0.115(195)	0.567(133)
13-pt $c_4^M\ne 0$	0.497(29)	0.454(61)	0.462(64)	0.461(64)	0.560(90)	0.443(65)	0.496(68)
13-pt $c_4^M=0$	0.482(19)	0.449(41)	0.458(49)	0.457(49)	0.575(77)	0.441(59)	0.520(54)
11-pt $c_4^M\ne 0$	0.495(29)	0.445(62)	0.450(65)	0.449(65)	0.548(92)	0.434(67)	0.486(69)
11-pt $c_4^M=0$	0.480(20)	0.436(43)	0.443(51)	0.441(51)	0.558(80)	0.432(63)	0.505(56)
10-pt $c_4^M\ne 0$	0.537(44)	0.471(100)	0.467(102)	0.466(102)	0.617(141)	0.363(102)	0.539(110)
10-pt $c_4^M=0$	0.497(28)	0.444(62)	0.447(71)	0.446(71)	0.610(112)	0.378(87)	0.549(80)
10∗-pt $c_4^M\ne 0$	0.503(32)	0.481(65)	0.484(67)	0.484(67)	0.563(94)	0.449(67)	0.505(71)
10∗-pt $c_4^M=0$	0.480(20)	0.431(43)	0.433(52)	0.432(51)	0.550(80)	0.407(64)	0.498(57)

Table 7. Table 7: Isovector magnetic moment of the nucleon, μ p − n ≡ μ p − μ n superscript 𝜇 𝑝 𝑛 subscript 𝜇 𝑝 subscript 𝜇 𝑛 \mu^{p-n}\equiv\mu_{p}-\mu_{n} , in units of the Bohr magneton from the seven fits to the isovector form factor G M ​ ( Q 2 ) subscript 𝐺 𝑀 superscript 𝑄 2 G_{M}(Q^{2}) and including the derived value for G M ​ ( 0 ) subscript 𝐺 𝑀 0 G_{M}(0) . The rest is the same as in Table  5 .

Ensemble	Dipole	$z^3$	$z^4$	$z^5$	$z^{3+4}$	$z^{4+4}$	$z^{5+4}$
$a15m310$	4.280(57)	4.295(57)	4.295(57)	4.295(57)	4.296(57)	4.296(57)	4.295(57)
$a12m310$	4.205(85)	4.303(94)	4.303(94)	4.303(94)	4.297(94)	4.310(94)	4.301(94)
$a12m220L$	4.215(65)	4.253(84)	4.253(84)	4.253(84)	4.257(85)	4.252(83)	4.253(84)
$a12m220$	4.103(125)	4.143(130)	4.143(130)	4.143(130)	4.134(130)	4.135(130)	4.141(130)
$a12m220S$	4.005(155)	4.256(182)	4.255(182)	4.255(182)	4.262(182)	4.262(182)	4.257(182)
$a09m310$	4.141(30)	4.141(31)	4.141(31)	4.141(31)	4.141(31)	4.140(31)	4.141(31)
$a09m220$	4.260(60)	4.292(66)	4.292(66)	4.292(66)	4.291(66)	4.292(67)	4.292(66)
$a09m130W$	4.086(71)	4.088(75)	4.087(75)	4.087(75)	4.095(74)	4.084(75)	4.089(75)
$a06m310$	4.044(149)	3.985(159)	3.985(159)	3.985(159)	3.989(157)	3.986(160)	3.985(159)
$a06m310W$	4.145(124)	4.163(126)	4.163(126)	4.163(126)	4.163(125)	4.161(126)	4.163(125)
$a06m220$	3.938(93)	3.941(94)	3.941(94)	3.941(94)	3.940(94)	3.941(94)	3.941(94)
$a06m220W$	4.119(131)	4.113(132)	4.113(132)	4.113(132)	4.113(132)	4.112(132)	4.113(132)
$a06m135$	4.100(115)	4.078(113)	4.077(113)	4.077(113)	4.078(112)	4.079(113)	4.076(113)
13-pt $c_4^{\mu }\ne 0$	3.962(79)	3.930(81)	3.929(81)	3.929(81)	3.932(81)	3.927(81)	3.930(81)
13-pt $c_4^{\mu }=0$	3.950(79)	3.918(80)	3.917(80)	3.917(80)	3.920(80)	3.915(80)	3.917(80)
11-pt $c_4^{\mu }\ne 0$	3.975(84)	3.940(86)	3.939(86)	3.939(86)	3.942(86)	3.937(86)	3.939(86)
11-pt $c_4^{\mu }=0$	3.968(84)	3.933(86)	3.932(86)	3.932(86)	3.935(85)	3.930(86)	3.933(86)
10-pt $c_4^{\mu }\ne 0$	4.167(151)	3.982(164)	3.982(164)	3.982(164)	3.987(164)	3.976(164)	3.982(164)
10-pt $c_4^{\mu }=0$	3.999(122)	3.879(129)	3.877(129)	3.877(129)	3.884(129)	3.874(129)	3.879(129)
10∗-pt $c_4^{\mu }\ne 0$	3.970(85)	3.942(86)	3.941(86)	3.941(86)	3.945(86)	3.940(86)	3.942(86)
10∗-pt $c_4^{\mu }=0$	3.964(84)	3.933(86)	3.932(86)	3.932(86)	3.936(85)	3.930(86)	3.933(86)

Table 8. Table 8: Values of the parameters, c i E , M , μ superscript subscript 𝑐 𝑖 𝐸 𝑀 𝜇 c_{i}^{E,M,\mu} , defined in Eqs. ( 25 ), ( 26 ) and ( 27 ) for the 11-point fit used to obtain ⟨ r E 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝐸 2 \mathop{\langle r_{E}^{2}\rangle}\nolimits , ⟨ r M 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝑀 2 \mathop{\langle r_{M}^{2}\rangle}\nolimits and μ 𝜇 \mu in the continuum limit from the dipole and z 4 superscript 𝑧 4 z^{4} data.

$⟨r_E^2⟩$	$c_1^E$	$c_2^E$	$c_3^E$	$c_4^E$
	(fm2)	(fm)	(fm2)	(fm2)
dipole	0.32(3)	0.62(12)	-0.08(1)	0.33(26)
$z^4$	0.31(5)	0.50(23)	-0.08(2)	0.11(69)
$⟨r_M^2⟩$	$c_1^M$	$c_2^M$	$c_3^M$	$c_4^M$
	(fm2)	(fm)	(fm2 GeV)	(fm2 GeV)
dipole	0.31(3)	0.32(18)	0.024(6)	-0.14(19)
$z^4$	0.28(6)	0.77(32)	0.023(15)	-0.09(48)
$\mu$	$c_1^{\mu }$	$c_2^{\mu }$	$c_3^{\mu }$	$c_4^{\mu }$
		(fm-1)	(GeV-1)	(GeV-1)
dipole	3.93(11)	2.28(89)	0.33(40)	-44(39)
$z^4$	3.91(11)	3.10(98)	0.22(42)	-51(45)

Table 9. Table 9: Our final results for ⟨ r E 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝐸 2 \mathop{\langle r_{E}^{2}\rangle}\nolimits , ⟨ r M 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝑀 2 \mathop{\langle r_{M}^{2}\rangle}\nolimits and μ 𝜇 \mu from the 11-point CCFV fit to the dipole and the z 4 superscript 𝑧 4 z^{4} analysis data. The determination of the second error in these two estimates is explained in the text. The combined analysis, defined by Eq. ( 28 ) with the z 4 superscript 𝑧 4 z^{4} truncation for the Q 2 superscript 𝑄 2 Q^{2} behavior, has a single overall error. The bottom half of the table gives results for the Dirac, ⟨ r 1 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 1 2 \mathop{\langle r_{1}^{2}\rangle}\nolimits , and Pauli, ⟨ r 2 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 2 2 \mathop{\langle r_{2}^{2}\rangle}\nolimits , radii obtained using Eq. ( 12 ).

	$⟨r_E^2⟩$	$\sqrt{⟨r_E^2⟩}$	$⟨r_M^2⟩$	$\sqrt{⟨r_M^2⟩}$	$\mu$
	(fm2)	(fm)	(fm2)	(fm)	(Bohr Magneton)
dipole fit	0.586(17)(13)	0.765(11)(8)	0.495(29)(41)	0.704(21)(29)	3.975(84)(125)
$z^4$ fit	0.591(41)(46)	0.769(27)(30)	0.450(65)(102)	0.671(48)(76)	3.939(86)(138)
Combined fit	0.564(114)	0.751(76)	0.459(189)	0.678(140)	3.922(83)
	$⟨r_1^2⟩$	$\sqrt{⟨r_1^2⟩}$	$⟨r_2^2⟩$	$\sqrt{⟨r_2^2⟩}$
	(fm2)	(fm)	(fm2)	(fm)
dipole fit	0.389(18)(15)	0.623(15)(12)	0.531(44)(63)	0.729(30)(43)
$z^4$ fit	0.396(42)(49)	0.629(33)(37)	0.469(90)(141)	0.685(66)(103)
Combined fit	0.370(115)	0.609(94)	0.490(258)	0.700(184)

Table 10. Table 10: Lattice parameters of calculations that have presented results for G E ​ ( Q 2 ) subscript 𝐺 𝐸 superscript 𝑄 2 G_{E}(Q^{2}) and G M ​ ( Q 2 ) subscript 𝐺 𝑀 superscript 𝑄 2 G_{M}(Q^{2}) at or near the physical pion mass.

Ensemble ID	$a$ (fm)	$M_{\pi }$ (MeV)	$L^3\times T$	$M_{\pi }^{\mathrm{val}}L$	$\tau /a$	$N_{\text{conf}}$	$N_{\mathrm{meas}}$	Action
$a09m130$ (this work)	0.0871(6)	138(1)	${64}^3\times 96$	3.90	$\{8,10,12,14,16\}$	1290	165,120	clover-on-2+1+1-HISQ
$a06m135$ (this work)	0.0570(1)	136(2)	${96}^3\times 192$	3.7	$\{16,18,20,22\}$	675	43,200	clover-on-2+1+1-HISQ
LHPC’17 Hasan et al. (2018)	0.093	135	${64}^3\times 64$	4.08	$\{10,13,16\}$	442	56,576	2+1-clover
ETMC’18 Alexandrou et al. (2018)	0.0809(4)	138(1)	${64}^3\times 128$	3.62	$\{12,14,16,18,20\}$	750	3K$-$48K	2+1+1-Twisted Mass
ETMC’17 Alexandrou et al. (2017b)	0.0938(3)	130(1)	${48}^3\times 96$	2.98	$\{10,12,14,16,18\}$	578$-$725	9K$-$64K	2-Twisted Mass
ETMC’18 Alexandrou et al. (2018)	0.0938(3)	130(2)	${64}^3\times 128$	3.97	$\{12,14,16\}$	333$-$1040	5K$-$17K	2-Twisted Mass
PACS’18 Ishikawa et al. (2018)	0.0846(7)	146	${96}^3\times 96$	6.01	$\left\{15\right\}$	200	12,800	2+1-clover
PACS’18A Shintani et al. (2019)	0.0846(7)	135	${128}^3\times 128$	7.41	$\{10,12,14,16\}$	20	2.5K$-$10K	2+1-clover

Table 11. Table 11: Results for r E subscript 𝑟 𝐸 \mathop{r_{E}}\nolimits , r M subscript 𝑟 𝑀 \mathop{r_{M}}\nolimits , μ 𝜇 \mu and the nucleon mass from published calculations at or near the physical pion mass. The quantity used to set the lattice scale is given in the third column, with r 0 2 ​ F ​ ( r 0 ) superscript subscript 𝑟 0 2 𝐹 subscript 𝑟 0 r_{0}^{2}F(r_{0}) and r 1 subscript 𝑟 1 r_{1} extracted from the heavy quark potential  Sommer ( 2014 ) . ETMC’18  Alexandrou  et al. ( 2018 ) results are derived from a single fit in Q 2 superscript 𝑄 2 Q^{2} to the combined 2- and 2+1+1-flavor data, i.e., neglecting the dependence on the number of flavors N f subscript 𝑁 𝑓 N_{f} and the difference in the lattice spacing a 𝑎 a . The LHPC’17  Hasan  et al. ( 2018 ) results are from a single ensemble and taken from their analysis using the summation method to control ESC. The calculation of the scale used in LHPC’17 is given in Ref.  Hasan  et al. ( 2019 ) , and that by the PACS collaboration in Ref.  Ishikawa  et al. ( 2016 ) .

Ensemble ID	$M_N$ (MeV)	$a$ from	$Q^2$ Fit	$r_E$ (fm)	$r_M$ (fm)	$\mu$
$a09m130W$	953(4)	$r_1$	$z^4$	0.769(27)(30)	0.671(48)(76)	3.94(9)(14)
$a06m135$	951(10)	$r_1$	$z^4$	0.765(11)(8)	0.704(21)(29)	3.98(8)(13)
LHPC’17 Hasan et al. (2018)	912(8)	$M_{\Omega }$	$z^5$	0.887(49)		4.75(15)
ETMC’18 Alexandrou et al. (2018)	929(6)	$r_0^{2}F\left(r_0\right)=1.65$	dipole	0.802(19)(12)(1)	0.714(26)(88)(16)(${}_0{}^{}{}_{}{}^1$)	3.96(14)(3)(7)(${}_0{}^{}{}_{}{}^1$)
ETMC’17 Alexandrou et al. (2017b)	941(2)	$r_0^{2}F\left(r_0\right)=1.65$	dipole	0.808(30)(19)	0.732(36)(45)	4.02(21)(28)
PACS’18 Ishikawa et al. (2018)	958(10)	$M_{\Omega }$	$z^8|z^7$	0.915(99)	1.437(409)	4.81(79)
PACS’18A Shintani et al. (2019)	942(11)	$M_{\Omega }$	dipole	0.875(15)(28)	0.805(32)(274)	4.417(138)(317)

Table 12. Table 12: Parameters, including the Goldstone pion mass M π sea superscript subscript 𝑀 𝜋 sea M_{\pi}^{\rm sea} , of the eleven 2+1+1- flavor HISQ ensembles generated by the MILC collaboration and analyzed in this study are quoted from Ref.  Bazavov  et al. ( 2013 ) . All fits are made versus M π val superscript subscript 𝑀 𝜋 val M_{\pi}^{\rm val} and finite-size effects are analyzed in terms of M π val ​ L superscript subscript 𝑀 𝜋 val 𝐿 M_{\pi}^{\rm val}L . Estimates of M π val superscript subscript 𝑀 𝜋 val M_{\pi}^{\rm val} , the clover-on-HISQ pion mass, are the same as given in Ref.  Bhattacharya  et al. ( 2015 ) and the error is governed mainly by the uncertainty in the lattice scale. In the last four columns, we give, for each ensemble, the values of the source-sink separation t sep subscript 𝑡 sep t_{\rm sep} used in the calculation of the three-point functions, the number of configurations analyzed, and the number of measurements made using the high precision (HP) and the low precision (LP) truncation of the inversion of the clover operator. The smearing size used in the calculation of the quark propagator is given in Table  13 .

Ensemble ID	$a$ (fm)	$M_{\pi }^{\mathrm{sea}}$ (MeV)	$M_{\pi }^{\mathrm{val}}$ (MeV)	$L^3\times T$	$M_{\pi }^{\mathrm{val}}L$	$\tau /a$	$N_{\text{conf}}$	$N_{\mathrm{meas}}^{\mathrm{HP}}$	$N_{\mathrm{meas}}^{\mathrm{LP}}$
$a15m310$	0.1510(20)	306.9(5)	320.6(4.3)	${16}^3\times 48$	3.93	$\{5,6,7,8,9\}$	1917	7668	122,688
$a12m310$	0.1207(11)	305.3(4)	310.2(2.8)	${24}^3\times 64$	4.55	$\{8,10,12\}$	1013	8104	64,832
$a12m220S$	0.1202(12)	218.1(4)	225.0(2.3)	${24}^3\times 64$	3.29	$\{8,10,12\}$	946	3784	60,544
$a12m220$	0.1184(10)	216.9(2)	227.9(1.9)	${32}^3\times 64$	4.38	$\{8,10,12\}$	744	2976	47,616
$a12m220L$	0.1189(09)	217.0(2)	227.6(1.7)	${40}^3\times 64$	5.49	$\{8,10,12,14\}$	1000	4000	128,000
$a09m310$	0.0888(08)	312.7(6)	313.0(2.8)	${32}^3\times 96$	4.51	$\{10,12,14,16\}$	2264	9056	114,896
$a09m220$	0.0872(07)	220.3(2)	225.9(1.8)	${48}^3\times 96$	4.79	$\{10,12,14,16\}$	964	7712	123,392
$a09m130W$	0.0871(06)	128.2(1)	138.1(1.0)	${64}^3\times 96$	3.90	$\{8,10,12,14,16\}$	1290	5160	165,120
$a06m310$	0.0582(04)	319.3(5)	319.6(2.2)	${48}^3\times 144$	4.52	$\{16,20,22,24\}$	1000	8000	64,000
$a06m31W$						$\{18,20,22,24\}$	500	2000	64,000
$a06m220$	0.0578(04)	229.2(4)	235.2(1.7)	${64}^3\times 144$	4.41	$\{16,20,22,24\}$	650	2600	41,600
$a06m22W$						$\{18,20,22,24\}$	649	2600	41,600
$a06m135$	0.0570(01)	135.5(2)	135.6(1.4)	${96}^3\times 192$	3.7	$\{16,18,20,22\}$	675	2700	43,200

Table 13. Table 13: The parameters used in the calculation of the clover propagators. The hopping parameter for the light quarks, κ l subscript 𝜅 𝑙 \kappa_{l} , in the clover action is given by 2 ​ κ l = 1 / ( m l + 4 ) 2 subscript 𝜅 𝑙 1 subscript 𝑚 𝑙 4 2\kappa_{l}=1/(m_{l}+4) . m l subscript 𝑚 𝑙 m_{l} is tuned to achieve M π val ≈ M π sea superscript subscript 𝑀 𝜋 val superscript subscript 𝑀 𝜋 sea M_{\pi}^{\rm val}\approx M_{\pi}^{\text{sea}} . The parameters used to construct Gaussian smeared sources  Güsken  et al. ( 1989 ) , { σ , N KG } 𝜎 subscript 𝑁 KG \{\sigma,N_{\text{KG}}\} , are given in the fourth column where N KG subscript 𝑁 KG N_{\text{KG}} is the number of applications of the Klein-Gordon operator and the width of the smearing is controlled by the coefficient σ 𝜎 \sigma , both in Chroma convention  Edwards and Joo ( 2005 ) . The resulting root-mean-square radius of the smearing, defined as ∫ r 2 ​ S † ​ S ​ 𝑑 r / ∫ S † ​ S ​ 𝑑 r superscript 𝑟 2 superscript 𝑆 † 𝑆 differential-d 𝑟 superscript 𝑆 † 𝑆 differential-d 𝑟 \sqrt{\int r^{2}\sqrt{S^{\dagger}S}dr/\int\sqrt{S^{\dagger}S}dr} , is given in the last column.

ID	$m_l$	$c_{\text{SW}}$	Smearing	RMS smearing
			Parameters	radius
$a15m310$	$-0.0893$	1.05094	{4.2, 36}	4.69
$a12m310$	$-0.0695$	1.05094	{5.5, 70}	5.96
$a12m220S$	$-0.075$	1.05091	{5.5, 70}	5.98
$a12m220$	$-0.075$	1.05091	{5.5, 70}	5.96
$a12m220L$	$-0.075$	1.05091	{5.5, 70}	5.96
$a09m310$	$-0.05138$	1.04243	{7.0,100}	7.48
$a09m220$	$-0.0554$	1.04239	{7.0,100}	7.48
$a09m130W$	$-0.058$	1.04239	{7.0,100}	7.50
$a06m310$	$-0.0398$	1.03493	{6.5, 70}	7.22
$a06m310W$	$-0.0398$	1.03493	{12, 250}	12.19
$a06m220$	$-0.04222$	1.03493	{5.5, 70}	6.22
$a06m220W$	$-0.04222$	1.03493	{11, 230}	11.24
$a06m135$	$-0.044$	1.03493	{9.0,150}	9.56

Table 14. Table 14: Nucleon ground and excited-state masses in GeV GeV \mathop{\rm GeV}\nolimits extracted from a 4-state fit.

ID	$M_0$	$M_1$	$M_2$	$M_3$
a15m310	1.0848(28)	2.038(62)	2.40(8)	2.88(8)
a12m310	1.0888(44)	1.576(89)	2.55(16)	3.18(16)
a12m220L	1.0165(35)	1.691(175)	2.76(26)	3.44(26)
a12m220	1.0133(52)	1.634(116)	2.75(27)	3.42(27)
a12m220S	0.9915(86)	1.499(78)	3.00(22)	3.66(22)
a09m310	1.1001(31)	2.065(128)	3.61(32)	4.78(33)
a09m220	1.0172(45)	1.718(91)	2.56(16)	3.44(17)
a09m130W	0.9532(39)	1.761(88)	2.98(15)	3.81(15)
a06m310	1.1014(102)	1.646(105)	2.79(16)	3.73(22)
a06m310W	1.1109(61)	2.054(145)	3.00(24)	3.99(27)
a06m220	1.0365(65)	1.874(73)	3.05(11)	3.96(19)
a06m220W	1.0345(72)	1.816(144)	2.70(24)	3.69(31)
a06m135	0.9512(100)	1.734(89)	3.01(13)	4.01(17)

Table 15. Table 15: Values for the parameters of the two CCFV fits to the nucleon mass defined in Eq. ( 29 ).

Fit	$c_0$	$c_1$	$c_2$	$c_3$	$c_4$	$c_5$	$M_N$	${\chi }^2$/dof
	$[GeV]$	$\left[GeV{fm}^{-1}\right]$	$\left[GeV{fm}^{-2}\right]$	$\left[{GeV}^{-1}\right]$	$\left[{GeV}^{-2}\right]$	$\left[{GeV}^{-1}\right]$	$[GeV]$	[$p$-value]
1	0.931(22)	-0.273(391)	0.34(1.96)	2.722(487)	-2.2(1.3)	-10.9(5.1)	0.9755(202)	0.45 [0.87]
2		0.430(80)	-3.14(53)	2.875(48)	-2.6(1.3)	-5.6(4.2)	0.939	0.81 [0.59]

Equations

$$\begin{gathered} ⟨N({\vec{p}}_f)|V_{\mu }^{\mathrm{em}}(\vec{q})|N({\vec{p}}_i)⟩= \\ {\overline{u}}_N({\vec{p}}_f)\left(F_1(Q^2){\gamma }_{\mu }+{\sigma }_{\mu \nu }{q}_{\nu }\frac{F_2(Q^2)}{2M_N}\right)u_N({\vec{p}}_i), \end{gathered}$$

$$\sum _s{u}_N(\vec{p},s){\bar{u}}_N(\vec{p},s)=\frac{E(\vec{p}){\gamma }_4-i\vec{\gamma }\cdot \vec{p}+M}{2E(\vec{p})}.$$

$$V_{\mu }^{\mathrm{em}}=\frac{2}{3}\overline{u}{\gamma }_{\mu }u-\frac{1}{3}\overline{d}{\gamma }_{\mu }d.$$

$$\begin{gathered} ⟨p({\vec{p}}_f)|\overline{u}{\gamma }_{\mu }u-\overline{d}{\gamma }_{\mu }d|p({\vec{p}}_i)⟩= \\ ⟨p({\vec{p}}_f)|V_{\mu }^{\mathrm{em}}(\vec{q})|p({\vec{p}}_i)⟩-⟨n({\vec{p}}_f)|V_{\mu }^{\mathrm{em}}(\vec{q})|n({\vec{p}}_i)⟩. \end{gathered}$$

$$G_E(Q^2)$$

$$G_M(Q^2)$$

$$g_V=G_E{|}_{Q^2=0}=F_1{|}_{Q^2=0}$$

$${\mu }^p-{\mu }^n=G_M{|}_{Q^2=0}=(F_1+F_2){|}_{Q^2=0}=1+{\kappa }_p-{\kappa }_n.$$

$${\kappa }_p$$

$${\kappa }_n$$

$$⟨r_{E,M}^2⟩=-6\frac{d}{d{Q}^2}{\left(\frac{G_{E,M}(Q^2)}{G_{E,M}(0)}\right)|}_{Q^2=0}.$$

$$⟨r_{1,2}^2⟩=-6\frac{d}{d{Q}^2}{\left(\frac{F_{1,2}(Q^2)}{F_{1,2}(0)}\right)|}_{Q^2=0}.$$

$$⟨r_1^2⟩$$

$$\kappa ⟨r_2^2⟩$$

$$r_E^p$$

$$r_M^p=0.776(38)\mathrm{fm}\mathrm{electrons}.$$

$$G_i(Q^2)=\frac{G_i(0)}{(1+Q^2/{\mathcal{M}}_i^2{)}^2}⟹⟨r_i^2⟩=\frac{12}{{\mathcal{M}}_i^2},$$

$$\frac{G_{E,M}(Q^2)}{G_E(0)}=\sum _{k=0}^{\infty }{a}_{k}z(Q^2{)}^k,$$

$$z=\frac{\sqrt{t_{\text{cut}}+Q^2}-\sqrt{t_{\text{cut}}+\overline{t_0}}}{\sqrt{t_{\text{cut}}+Q^2}+\sqrt{t_{\text{cut}}+\overline{t_0}}},$$

$$\sum _{k=n}^{k_{\max }}k(k-1)\dots (k-n+1)a_k=0,n=0,1,2,3.$$

$$C_{\Gamma }^{(3\text{pt})}$$

$$|{\mathcal{A}}_0^{\prime }||{\mathcal{A}}_0|⟨0^{\prime }|{\mathcal{O}}_{\Gamma }|0⟩e^{-E_{0}t-M_0(\tau -t)}+$$

$$|{\mathcal{A}}_0^{\prime }||{\mathcal{A}}_1|⟨0^{\prime }|{\mathcal{O}}_{\Gamma }|1⟩e^{-E_{0}t-M_1(\tau -t)}+$$

$$|{\mathcal{A}}_1^{\prime }||{\mathcal{A}}_0|⟨1^{\prime }|{\mathcal{O}}_{\Gamma }|0⟩e^{-E_{1}t-M_0(\tau -t)}+$$

$$|{\mathcal{A}}_1^{\prime }||{\mathcal{A}}_1|⟨1^{\prime }|{\mathcal{O}}_{\Gamma }|1⟩e^{-E_{1}t-M_1(\tau -t)}+$$

$$|{\mathcal{A}}_0^{\prime }||{\mathcal{A}}_2|⟨0^{\prime }|{\mathcal{O}}_{\Gamma }|2⟩e^{-E_{0}t-M_2(\tau -t)}+$$

$$|{\mathcal{A}}_2^{\prime }||{\mathcal{A}}_0|⟨2^{\prime }|{\mathcal{O}}_{\Gamma }|0⟩e^{-E_{2}t-M_0(\tau -t)}+$$

$$|{\mathcal{A}}_1^{\prime }||{\mathcal{A}}_2|⟨1^{\prime }|{\mathcal{O}}_{\Gamma }|2⟩e^{-E_{1}t-M_2(\tau -t)}+$$

$$|{\mathcal{A}}_2^{\prime }||{\mathcal{A}}_1|⟨2^{\prime }|{\mathcal{O}}_{\Gamma }|1⟩e^{-E_{2}t-M_1(\tau -t)}+$$

$$|{\mathcal{A}}_2^{\prime }||{\mathcal{A}}_2|⟨2^{\prime }|{\mathcal{O}}_{\Gamma }|2⟩e^{-E_{2}t-M_2(\tau -t)},$$

$$\chi (x)={ϵ}^{abc}\left[{q_1^a}^T(x)C{\gamma }_5\frac{(1\pm {\gamma }_4)}{2}{q}_2^b(x)\right]q_1^c(x)$$

$${\mathcal{R}}_{\mu }$$

$${\left[\frac{C^{(2\text{pt})}(t,{\mathbf{p}}^{\prime })C^{(2\text{pt})}(\tau ,{\mathbf{p}}^{\prime })C^{(2\text{pt})}(\tau -t,\mathbf{p})}{C^{(2\text{pt})}(t,\mathbf{p})C^{(2\text{pt})}(\tau ,\mathbf{p})C^{(2\text{pt})}(\tau -t,{\mathbf{p}}^{\prime })}\right]}^{1/2},$$

$$\sqrt{2E_p(M_N+E_p)}\mathrm{Re}({\mathcal{R}}_i)=$$

$$\sqrt{2E_p(M_N+E_p)}\mathrm{Im}({\mathcal{R}}_i)=$$

$$\sqrt{2E_p(M_N+E_p)}\mathrm{Re}({\mathcal{R}}_4)=$$

$$⟨r_E^2⟩(a,M_{\pi },L)=$$

Full text

nucleon form factors, lattice QCD, charge radii

††preprint: LA-UR-19-25275††preprint: MSUHEP-19-006

I Introduction

Experiments studying electron scattering off protons and neutrons have a long history of providing an understanding of the structure of nucleons Perdrisat et al. (2007); Punjabi et al. (2005). Quantitative understanding of the distribution of charge is described by the electric and magnetic form factors, $G_{E}$ and $G_{M}$ Bernauer et al. (2014). Quantities of phenomenological interest obtained from the slope of the form factors at space-like four momentum transfer squared $Q^{2}=0$ are the electric and magnetic charge radii of the nucleons. At present there is a 6$\sigma$ discrepancy between the electric charge radius of the proton obtained from electronic energy levels combined with electron scattering data Mohr et al. (2016) versus that from the Lamb shift in muonic hydrogen $E_{\mu p}(2S-2P)$ Antognini et al. (2013, 2016). A second issue that needs resolution is the behavior of the ratio $G_{E}/G_{M}$ at $Q^{2}>1$ GeV2 Ye et al. (2018), and whether this ratio crosses zero at about $8$ GeV2 as indicated by experiments at JLab Gayou et al. (2002); Electro-Magnetic Form Factor collaboration, Jefferson Lab() (EMFF). In this work, we focus on determining the electromagnetic form factors in the range $0.05\lesssim Q^{2}\lesssim 1$ GeV2 and extracting the charge radii from them.

The electric and magnetic form factors, $G_{E}$ and $G_{M}$, of the nucleon can be calculated directly from large scale simulations of lattice QCD. In recent years, advances in algorithms and computing power have allowed the community to push the calculations towards physical masses for the light $u$ and $d$ quarks, and on lattice spacings that are small enough that discretization effects are expected to be at the few percent level Bhattacharya et al. (2014); Green et al. (2014); Capitani et al. (2015); Alexandrou et al. (2017a). In this paper we present results from thirteen calculations on eleven ensembles that cover a range of lattice spacings (0.06\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }a\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }0.15 fm), pion masses (135\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }M_{\pi}\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }320 MeV) and lattice volumes (3.3\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }M_{\pi}L\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }5.5). These ensembles were generated using $2+1+1$-flavors of highly improved staggered quarks (HISQ) Follana et al. (2007) by the MILC collaboration Bazavov et al. (2013). This suite of calculations allows us to understand and assess various sources of systematic errors. The analysis is carried out using both the dipole ansatz and the $z$-expansion, which give consistent estimates for the isovector mean-square charge radii $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and the magnetic moment $\mu^{p-n}$.

Our final results for the isovector mean-square charge radii $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ (also for Dirac, $\mathop{\langle r_{1}^{2}\rangle}\nolimits$, and Pauli, $\mathop{\langle r_{2}^{2}\rangle}\nolimits$, radii) defined in Eqs. (10), (11) and (12), and for the magnetic moment $\mu$ are given in Tab. 9. We also present a comparison with other lattice data obtained close to the physical pion mass and with the Kelly parameterization of the experimental data Kelly (2004) in Fig. 22. Our estimates for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu^{p-n}$ are about 17%, 19%, and 16% smaller than the phenomenological values given in Eq. (30) and the precise experimental value in Eq. (9). Throughout this paper, we have paid attention to the size of possible statistical and systematic errors, and find that a linear combination of these is large enough to explain the deviations.

We analyze the world data for $G_{E}$ and $G_{M}$ in Sec. VII and find that data from all 13 of our calculations and those from other collaborations done at or near the physical pion mass fall roughly onto a single curve when plotted versus versus $Q^{2}$ or $Q^{2}/M_{N}^{2}$. However, there is a noticeable shift between the two curves when compared to the Kelly fit. The difference between the two ways of analyzing the data is a discretization artifact: specifically, it is a consequence of the difference in values of the lattice scale obtained from different observables. The size of the difference again indicates that the present underestimate of $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ should not be considered significant. Our overall conclusion is that to significantly reduce the systematics and improve the precision with which these observables can be extracted will require high statistics data at smaller values of the lattice spacing and of $Q^{2}<0.1$ GeV2.

We stress that the long-term goal of lattice QCD is to directly predict the form factors and not to reproduce the Kelly curve, a parameterization of the experimental data. Throughout this paper, we use the Kelly curve to provide a reference point for comparison, and for discussing systematics and trends in the lattice data. We do not show an error band on the Kelly curve as it is negligible on the scale of the errors in the lattice data.

This paper is organized as follows. In Sec. II, we review the theory, computational approach and the status of the experimental and phenomenological results. In Sec. III, we describe the salient features of the calculation. The fits used to isolate excited-state contamination (ESC) and extract the form factors are described in Sec. IV. Fits to quantify the $Q^{2}$ behavior of the $(p-n)$ form factors are discussed in Sec. V, and the extraction of our final results for the isovector mean-square charge radii, $\langle r_{E}^{2}\rangle$ and $\langle r_{M}^{2}\rangle$, and the anomalous magnetic moment $\mu^{p-n}$ are presented in Sec. VI. Comparisons with form factors extracted from experiments and with previous lattice QCD calculations are made in Sec. VII. We end with conclusions in Sec. VIII. Some further details of the calculations are given in four Appendices: lattice parameters in Appendix A, analysis of nucleon mass in Appendix B, ESC in Appendix C, and a review of the experimental data for the form factors in Appendix D.

II Electromagnetic Form Factors of the Nucleon

The Dirac, $F_{1}$, and Pauli, $F_{2}$, form factors are extracted from the matrix elements of the electromagnetic current within the nucleon state $N$ through the relation

[TABLE]

where $\vec{q}=\vec{p}_{f}-\vec{p}_{i}$ is the momentum transfer. The discrete lattice momenta are given by $2\pi{\bf n}/La$ with the entries of the vector ${\bf n}\equiv(n_{1},n_{2},n_{3})$ taking on integer values, $n_{i}\in\{0,L\}$. The spacing between the momenta is controlled by the spatial lattice size, $La$. The normalization used for the nucleon spinors in Euclidean space is

[TABLE]

and in Eq. (1), the electromagnetic current is

[TABLE]

In the isospin symmetric limit, the difference of its matrix elements between a proton and a neutron state are related to the isovector form factors of the proton by the relation

[TABLE]

The quantity we calculate on the lattice is the left hand side of Eq. (4), i.e., the isovector form factors of the proton. Throughout this paper, the term isovector form factors of the proton and the $(p-n)$ form factors refer to the same quantities as defined in Eq. (4). These will henceforth be analyzed in terms of the space-like 4-momentum squared, $Q^{2}={\vec{p}}^{2}-(E-m)^{2}=-q^{2}$.

Another common set of definitions of the electromagnetic form factors, widely used in the analysis of experimental data, are the Sachs electric, $G_{E}$, and magnetic, $G_{M}$, form factors that are related to the Dirac and Pauli form factors as

[TABLE]

From these, the vector charge is given by

[TABLE]

and the difference between the magnetic moment of the proton and the neutron by

[TABLE]

The anomalous magnetic moments of the proton and the neutron, in units of the Bohr magneton, are known very precisely Patrignani et al. (2016):

[TABLE]

The electric and magnetic size of the nucleon are defined as the slope of the form factors with respect to $Q^{2}$ at $Q^{2}=0$ Miller (2019):

[TABLE]

The form factors $G_{E,M}$ are normalized by their values at $Q^{2}=0$: $G_{E}(Q^{2}=0)\equiv g_{V}$ and $G_{M}(Q^{2}=0)/g_{V}\equiv\mu$. This definition makes them independent of the renormalization constant, $Z_{V}$, of the lattice vector current, and improves the signal because some of the systematics cancel in the ratios. Therefore, in this work, we will use Eq. (10) when calculating $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$. Note that $Z_{V}g_{V}=1$ as the electric charge is conserved. A second independent estimate of $Z_{V}$, obtained using nonperturbative lattice calculations in the RI-sMOM scheme, is given in Ref. Gupta et al. (2018), where the difference between the two estimates was shown to be $\lesssim 3\%$.

One similarly defines the isovector Dirac and Pauli mean-square radii as

[TABLE]

These are related to $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu\equiv 1+\kappa$ as

[TABLE]

Our analysis of the lattice data is carried out in terms of $G_{E}$ and $G_{M}$. Results for $\mathop{\langle r_{1}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{2}^{2}\rangle}\nolimits$ are also given in Table 9 in Sec. VI, where we extract $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$.

The electric root-mean-square charge radius $\mathop{r_{E}}\nolimits\equiv\sqrt{\langle r_{E}^{2}\rangle}$ of the proton has been measured in three ways: (i) laser spectroscopy of the Lamb shift in muonic hydrogen Antognini et al. (2013, 2016); Krauth et al. (2017), (ii) continuous-wave laser spectroscopy of hydrogen Pohl et al. (2017), and (iii) elastic scattering of electrons off protons Sick and Trautmann (2014, 2017). Results using electrons, i.e., the latter two ways, are included in the CODATA-2014 world average Tanabashi et al. (2018); Mohr et al. (2016):

[TABLE]

and the third result is from muonic hydrogen. The large difference between the CODATA-2014 and muonic-hydrogen values was termed the “proton radius puzzle”. The new CODATA-2018 value CODATA-2018 resolves the puzzle in favor of the muonic-hydrogen result. The magnetic radius of the proton extracted from experiments using electrons is Tanabashi et al. (2018); Mohr et al. (2016)

[TABLE]

Values for the isovector charge radii, extracted from the experimental data and used to compare lattice data against, are given in Eq. (30) in Appendix D.

To reduce the uncertainty in results from electron scattering experiments, which have been done down to $Q^{2}\approx 0.004$ GeV2, new experiments to constrain the low $Q^{2}$ behavior have been initiated Weber et al. (2016); Gasparian (2017). Similarly, for lattice QCD calculations to help resolve the puzzle, we need to calculate the form factors to $Q^{2}\approx 0.004$ GeV2 to extract $\mathop{r_{E}}\nolimits$ with better than 1% accuracy.

A challenge to the direct extraction of $\langle r_{i}^{2}\rangle$ from the lattice data is that the value of the smallest momenta, $2\pi/La$, is large in typical lattice simulations. In our calculations, it is \raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle>}}{{\sim}} }220 MeV, and the range of $Q^{2}$ values, given in Table 1, are between 2–10 $M_{\pi}^{2}$. It is, therefore, traditional to fit the data for the $G_{i}$ to an ansatz, and then use the fit to evaluate the derivative given in Eq. (10). Both, using an ansatz and estimating its parameters from fits to data with \sqrt{Q^{2}}\ \raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle>}}{{\sim}} }200 MeV introduces systematic uncertainties when evaluating the derivative at $Q^{2}=0$. We estimate the dependence of $\langle r_{i}^{2}\rangle$ on the choice of the ansatz by comparing results for each ensemble obtained using two different fits, the dipole model and the $z$-expansion.

Two alternate approaches are, one, to calculate the form factors at fixed $Q^{2}$ and extrapolate these to the continuum limit first and then fit the $Q^{2}$ behavior. Unfortunately, the values of $Q^{2}$ are different on each ensemble. Second, combine the dipole or the $z$-expansion parameterization of the $Q^{2}$ behavior with the chiral-continuum-finite volume (CCFV) ansatz for one overall fit. This combined fit is discussed in Sec. VI.2. The central analysis presented here consists of first fitting the data versus $Q^{2}$ using the dipole model and the $z$-expansion to extract $\langle r_{i}^{2}\rangle$ and $\mu$ on each ensemble and then get the physical results from a CCFV fit in $a$, $M_{\pi}$ and $M_{\pi}L$ that addresses the associated systematics.

It is important to note that both the electron scattering experiments and lattice QCD calculations suffer from paucity of data close to $Q^{2}=0$ that impacts the extraction of the charge radii. However, there is a large range, $0.004\lesssim Q^{2}\lesssim 1$ GeV2 over which accurate experimental data exist. Thus, more than just extracting the charge radii, our goal is to directly compare the lattice and the experimental data over this range of $Q^{2}$ as discussed in Sec. V.

An ansatz that is commonly used to fit the experimental data is the dipole. It arises if one assumes an exponentially falling charge distribution. The resulting form factor is characterized by a single parameter, the mass ${\cal M}$,

[TABLE]

and normalized to $F_{1}=G_{E}=g_{V}$ at $Q^{2}=0$. It goes as $Q^{-4}$ in the $Q^{2}\to\infty$ limit in accord with perturbation theory Lepage and Brodsky (1980).

The second ansatz is a model-independent parameterization called the $z$-expansion Hill and Paz (2010); Bhattacharya et al. (2011):

[TABLE]

where the $a_{k}$ are fit parameters and $z$ is defined as

[TABLE]

with $t_{\text{cut}}=4M_{\pi}^{2}$ denoting the nearest singularity in $G_{E,M}(Q^{2})$. In terms of $z$, the domain of analyticity of $G_{E}(Q^{2})$ is mapped into the unit circle with the branch cut at $Q^{2}=-4M_{\pi}^{2}$ Bhattacharya et al. (2011). We analyzed the data with $\bar{t_{0}}=0$ and $\bar{t_{0}}^{\rm mid}=\{0.12,\ 0.20,\ 0.40\}$ GeV2 for the $M_{\pi}\approx\{135,\ 220,\ 315\}$ MeV ensembles. By choosing the value of the constant $\bar{t_{0}}$ to lie in the middle of the range of $Q^{2}$ at which we have data, one reduces $z_{\rm max}$. By reducing the value of $z_{\rm max}$ we hope to improve the stability of the estimates, with improvement judged by comparing result from different truncations of the series. In practice, for our data set, we find that the quality of the fits and the results are insensitive to the choice of $\bar{t_{0}}$. The final results for the charge radii and magnetic moment are obtained from fits using $\bar{t_{0}}^{\rm mid}$.

The values of $Q^{2}$ for the thirteen calculations are given in Table 1. Note that in four cases the number of nonzero values are only five. The data for $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ versus $z$ with $\bar{t_{0}}^{\rm mid}$ are shown in Fig. 8. As discussed in Sec. V.2, we restrict our fits to $Q^{2}\leq 1$ GeV2 because the reliability of some of the higher $Q^{2}$ data is questionable.

To implement the perturbative behavior $G_{i}(Q^{2})\to Q^{-4}$ as $Q^{2}\to\infty$ Lepage and Brodsky (1980) in the $z$-expansion requires $Q^{n}G_{i}(Q^{2})\to 0$ for $n=0,1,2,3$. These constraints can be incorporated into the z-expansion as four sum rules Lee et al. (2015)

[TABLE]

For $n=0$ it reduces to $\sum_{k=0}^{k_{\rm max}}a_{k}=0$. A priori, using these sum rules ensures that the $a_{k}$ are not only bounded but must also decrease at large $k$ Lee et al. (2015).

A key issue in the $z$-expansion analysis is the value of $k^{\rm max}$ required to obtain results with a certain precision. The analysis of the experimental data carried out in Appendix D shows that results stabilize for $k_{\rm max}\approx 4$ with and without sum rules. For the lattice data, the choice has to take into account the number of values of $Q^{2}$ at which data have been generated to not over-parameterize the fit. For our data and fits without priors, the $a_{k}$ fluctuate and the higher order coefficients ($k\geq 4$) are ill determined due to the over-parameterization of the fits. To avoid the resulting large fluctuations in $a_{k}$, we put a bound on them as suggested in Lee et al. (2015). For $G_{E}$ and $G_{M}/5$, we constrain $|a_{k}|\lesssim 5.0$ for all $k$ by using Gaussian priors with central value zero and width five. With this constraint, results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ do not change significantly for $k_{\rm max}\leq 3$ and stabilize for $k_{\rm max}\geq 4$ as shown in Fig. 9. The convergence of estimates from fits with sum rules is slower and occurs for $k_{\rm max}\geq 7$ as also shown in Fig. 9. We, therefore, use the fits with $k_{\rm max}=4$ and without sum rules for our final results as they converge faster. Results with sum rules, which converge for $k_{\rm max}\geq 7$, are used only as consistency checks. Since $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ are best extracted from data at small $Q^{2}$, the sum rule constraints imposed to guarantee the large $Q^{2}$ behavior are not essential for their determination.

Overall, the fits to $G_{E}(Q^{2})$ are more stable than those to $G_{M}(Q^{2})$. The main reason is the extra data point at $G_{E}(Q^{2}=0)$ which pins down the sign of the slope of $G_{E}(Q^{2})$ at small $Q^{2}$. Using a value for $G_{M}(0)$, derived from the ratio $G_{M}(Q^{2})/G_{E}(Q^{2})$ as discussed in Sec. IV.2, greatly improved the stability of fits to $G_{M}(Q^{2})$.

III Lattice Methodology

The parameters of the thirteen calculations done on eleven HISQ ensembles are the same as used in Ref. Gupta et al. (2018) for the calculation of isovector charges. To keep the paper self-contained, the lattice parameters of the calculations and the number of measurements made are summarized in Table 12 in the appendix A. The parameters used to generate the Wilson-clover quark propagators using the multigrid algorithm Babich et al. (2010) are also given in Table. 13. We remind the reader that two ensembles, $a06m310$ and $a06m220$, have been analyzed twice with different smearing parameters giving a total of 13 calculations. Also, compared to Refs. Bhattacharya et al. (2016); Gupta et al. (2017), six ensembles ($a12m220S$, $a12m220$, $a12m220L$, $a09m310$, $a09m220$ and $a09m130W$) have been simulated afresh with randomly chosen source points on each configuration to increase their statistical independence, and data at a larger number of momenta have been accumulated.

To increase the statistics cost-effectively, we used the truncated solver with bias correction method Bali et al. (2010); Blum et al. (2013). We also used the coherent source method to construct sequential propagators from the sink time slice, at which a zero-momentum nucleon state is inserted Bratt et al. (2010); Yoon et al. (2016).

The details of our strategy for the calculations and the analysis have been published in earlier works Bhattacharya et al. (2016); Gupta et al. (2017, 2018). Here we provide a brief summary of the points relevant to the calculation of the electric and magnetic form factors:

• •

All errors are determined using a single elimination jackknife method over configurations, i.e., we first construct the bias corrected average for each configuration and then carry out the fits to the two- and three-point functions within the same jackknife procedure over these configuration averages.

• •

To control excited-state contamination, we use the same toolkit as in Ref. Gupta et al. (2018). The 2-point functions are fit keeping four states in the spectral decomposition. The amplitudes and the masses obtained from these fits are input into the analysis of three-point functions. The results for the masses are given in Table 14 in Appendix B.

• •

On each ensemble, we calculate the three-point functions at multiple values of source-sink separation $\tau$. These values of $\tau$, given in Table 12, are the same as in Ref. Gupta et al. (2018).

• •

The insertion of the vector current at definite momenta $\bm{p}$ is carried out on each time slice $t$ between the source and the sink, and for each value of $\tau$. These data for the three-point functions, $C_{\Gamma}^{(3\text{pt})}(t;\tau;\bm{p}^{\prime},\bm{p})$, at a large number of values of $t$ and $\tau$ are fit using three states in the spectral decomposition:

[TABLE]

where the source point is translated to $t=0$, the operator is inserted at time $t$, and the nucleon state is annihilated at the sink time slice $\tau$, which numerically is also the source-sink separation. In this relation, the numbers refer to the state $|n\rangle$, a state with superscript ′ denotes that it could have nonzero momentum $\bm{p}^{\prime}$, and the momentum $\bm{p}$ at the sink is fixed to zero.

• •

With our data, the term $\langle 2^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 2\rangle$ could not be resolved. So, in all the fits we set the contribution of the term with $\langle 2^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 2\rangle$ equal to zero, and call these $3^{\ast}$-state fits.

• •

In the case of $a12m220S$ data, the $\bm{p}^{\prime}=0$ data are analyzed using $3^{\ast}$-state fits, while the $\bm{p}^{\prime}\neq 0$ data are fit using two states because the $3^{\ast}$-state fits for $Q^{2}\neq 0$ are unstable. Having stated this caveat, we will, for brevity, use the label $3^{\ast}$-state to describe the excited-state fits to all data, even those for this ensemble.

• •

The values of $Q^{2}$ at which the form factors are calculated are collected in Table 1. These are obtained using the nucleon ground-state energy $E_{p}$ extracted using 4-state fits to the 2-point functions.

• •

To extract the desired matrix element $\langle 0^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 0\rangle$ using Eq. (19), the masses $M_{i}$, energies $E_{i}$, and the amplitudes $\lvert\mathcal{A}_{i}\rvert$ and $\lvert\mathcal{A}_{i}^{\prime}\rvert$ are taken from the fit to the two-point function within one overall jackknife procedure. This assumes that the ordering of the coupling to the excited states is the same as in two-point functions. To improve the signal, the amplitude ${\cal A}_{0}^{\prime}$ with which the nucleon interpolating operator at the source time slice couples to the ground state $|0^{\prime}\rangle$ with energy $E_{0}$ and momentum $\bm{p}^{\prime}$ should be large while the coupling to excited states should be small. We find that for the smearing parameters given in Table. 13, the signal in all the ten momentum channels analyzed is good.

• •

Off diagonal terms with nonzero momentum transfer such as $\lvert\mathcal{A}_{i}^{\prime}\rvert\lvert\mathcal{A}_{j}\rvert\langle i^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert j\rangle$ are related to $\lvert\mathcal{A}_{j}^{\prime}\rvert\lvert\mathcal{A}_{i}\rvert\langle j^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert i\rangle$ by a combination of Lorentz boost, parity and hermitian transformation provided the tower of states and the coupling to them are the same on either side of the operator. In our calculation, the nucleon operator used is

[TABLE]

with color indices $\{a,b,c\}$, charge conjugation matrix $C=\gamma_{0}\gamma_{2}$, and $q_{1}$ and $q_{2}$ denoting the two different flavors of light Dirac quarks. The quark propagator is smeared both at the source and the sink using a gauge invariant Gaussian smearing procedure Güsken et al. (1989) described in Appendix A. The nonrelativistic projection $(1\pm\gamma_{4})/2$, inserted to improve the signal Bhattacharya et al. (2015, 2016); Yoon et al. (2017), as well as the smearing of the quark fields, breaks Lorentz covariance. Also, the sink is explicitly constructed to have $\vec{p}=0$. We, therefore, treat all such pair of matrix elements as independent free parameters in the fits.

• •

The data for 3-point functions at nonzero momentum transfer are not symmetric about the midpoint, $\tau/2$, between the source and the sink. Nevertheless, in the simultaneous 3-state fit to the data with multiple source-sink separations $\tau$ and intermediate times $t$, we skip the same $\mathop{t_{\rm skip}}\nolimits$ points adjacent to the source and the sink for every $\tau$ to remove points with the largest ESC. Two considerations motivated this choice: (i) the time slice of the onset of the plateau in the nucleon effective mass plot is roughly independent of the momentum as shown in Refs. Gupta et al. (2017, 2018), and (ii) because we choose the values of $\mathop{t_{\rm skip}}\nolimits$ to be as small as possible based on the stability of the covariance matrix used in the fits. The values of $\mathop{t_{\rm skip}}\nolimits$ used here are the same as in Ref. Gupta et al. (2018).

• •

The vector current in the continuum theory is conserved, however the local vector current used in our lattice calculations is not. The renormalization constant $Z_{V}$ for this current has been determined in two ways: (i) nonperturbatively in the RI-sMOM scheme and then converted to $\overline{MS}$ using perturbation theory and (ii) measured directly from the matrix element of $V_{4}$ at $Q^{2}=0$, i.e., $1/g_{V}$. The two sets of values are compared in Ref. Gupta et al. (2018) and differ by up to 3%. This size of difference is not unreasonable in our clover-on-HISQ formulation which has discretization effects starting at $O(\alpha_{s}a)$. Here, we implement method (ii) by forming ratios $G_{i}(Q^{2})/G_{E}(0)$, in which some of the systematics cancel. The discretization errors in $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ are addressed by the continuum extrapolation, a part of the CCFV fit.

The key input, other than statistical precision of the 3-point data, that impacts the stability of the n-state fits to control ESC and obtain the ground state matrix elements is the energy of the first excited state since the terms with $\langle 1^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 0\rangle$ and $\langle 0^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 1\rangle$ give the dominant contribution. Once the ground-state matrix elements have been determined, the procedure for obtaining the form factors from them is described in the next section.

IV Extracting form factors from matrix elements

The following ratios, ${\cal R}_{\mu}$, of the three-point to the two-point correlation functions,

[TABLE]

give the desired ground state matrix elements (ME) $\langle 0^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 0\rangle$, introduced in Eq. 19, in the limits $t\to\infty$ and $(\tau-t)\to\infty$. In the calculation of the nucleon three-point functions, we use the spin projection operator ${\cal P}_{3}=(1+\gamma_{4})(1+i\gamma_{5}\gamma_{3})/2$. With this ${\cal P}_{3}$, and the vector current defined in Eqs. (3) and (4) with Euclidean $\gamma_{\mu}$, the following quantities have a signal and give either the electric or the magnetic form factors:

[TABLE]

Note that, in practice, these ratios are used only to plot the data. Our results are obtained by making n-state fits to the correlation functions.

Exploiting the cubic symmetry under spatial rotations, we construct two averages over equivalent 3-point correlators before doing fits to get the ground state matrix elements: over $\mathop{\rm Re}({C}_{1})$ and $\mathop{\rm Re}({C}_{2})$ for $G_{M}(Q^{2})$ and over $\mathop{\rm Im}({C}_{1})$, $\mathop{\rm Im}({C}_{2})$ and $\mathop{\rm Im}({C}_{3})$ for $G_{E}(Q^{2})$. We label these form factors as $G_{M}^{V_{i}}$ and $G_{E}^{V_{i}}$. Together with $G_{E}^{V_{4}}$ extracted from Eq. (24), they constitute the three form factors analyzed. Their extraction is straightforward as each of the three is given by a distinct three-point function. It is important to note that the discretization artifacts and the excited-state contaminations in each can be very different.

The data for the ratio defined in Eq. (21) and the results of $3^{\ast}$ fits to the three 3-point correlators are illustrated in Figs. 24–30 and Figs. 33–34. The ideal expected behavior of all 3-point functions with large $t$ and $\tau-t$, is a flat region near $\tau/2$ that becomes independent of $\tau$. Our data show that this is not manifest even at $\tau\approx 1.4$ fm. We, therefore, use $3^{\ast}$-state fits to data at the various values of $t$ and $\tau$ to obtain estimates of the ground state matrix elements. Results for the three sets of form factors, $G_{E}^{V_{4}}$, $G_{E}^{V_{i}}$ and $G_{M}^{V_{i}}$, extracted from these matrix elements using Eqs. (22), (23) and (24) are given in Tables 2, 3 and 4 for the thirteen calculations.

IV.1 Extraction of $G_{E}(Q^{2})$

The pattern of the ESC in the extraction of $G_{E}^{V_{i}}$ versus $G_{E}^{V_{4}}$ can be, and is found to be, very different as shown in Figs. 24 and 25 for the two physical mass ensembles. The data for $G_{E}^{V_{4}}$ show a clear monotonic but slow convergence from above, and a flattish region near the middle. The estimates of the $\tau\to\infty$ values given by the $3^{\ast}$ fits are found to be stable under variations in $\mathop{t_{\rm skip}}\nolimits$ and the values of $\tau$ included in the fits.

The data for $G_{E}^{V_{i}}$ show much larger ESC and the ME $\langle 0^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 1\rangle$ and $\langle 1^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 0\rangle$ are an order of magnitude larger for $n^{2}=1$ as compared to those from $G_{E}^{V_{4}}$. The resulting pattern versus $t$ is essentially linear for each $\tau$. As $\tau$ is increased, this “line” rotates towards becoming flat, but the rotation is slow. The pivot point is approximately the point of intersection of the various $\tau$ lines and converges to the ground state estimate as $t$ and $(\tau-t)\to\infty$.

The difference in the shape of the ESC between $G_{E}^{V_{i}}$ and $G_{E}^{V_{4}}$ can be explained by the behavior of the transition matrix elements under parity transformation and hermitian conjugation. The imaginary parts of the matrix elements of $V_{i}$ at nonzero momentum pick up a negative sign under the combined transformations. As a result, for example, the term $\lvert\mathcal{A}_{0}^{\prime}\rvert\lvert\mathcal{A}_{1}\rvert\langle 0^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 1\rangle e^{-E_{0}t-M_{1}(\tau-t)}$ has opposite sign to that of its partner $\lvert\mathcal{A}_{1}^{\prime}\rvert\lvert\mathcal{A}_{0}\rvert\langle 1^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 0\rangle e^{-E_{1}t-M_{0}(\tau-t)}$. Thus, each such pair of terms give a “sinh”-like correction, that makes the data looks like a straight line at an angle to the extracted ground-state result. On the other hand, the matrix elements in the related pairs of terms from the real parts of $V_{i}$ and $V_{4}$ have the same sign, and therefore exhibit a “cosh”-like correction. Even in this case, the magnitudes of the two ME in such pairs of terms are not the same. Therefore, in fits to the three-point data using Eq. (19), we leave all the matrix elements as free parameters. In fact, in practice, it is the product of the amplitudes and the ME, such as $\lvert\mathcal{A}_{0}^{\prime}\rvert\lvert\mathcal{A}_{1}\rvert\langle 0^{\prime}\lvert\mathcal{O}_{\Gamma}\rvert 1\rangle$, that are free parameters in the fits. In these cases, only the energies are free parameters and these are taken from the two-point functions.

It is also evident from Figs. 24 and 25 that the ESC in $G_{E}^{V_{i}}$ is the largest at the smallest nonzero momentum, i.e., the “angle” the data make with the horizontal line is the largest. On the other hand, the ESC in $G_{E}^{V_{4}}$ increases with momentum. By comparing the data in the two figures, we also conclude that the ESC increases with decreasing $a$ for both $G_{E}^{V_{i}}$ and $G_{E}^{V_{4}}$.

A consequence of this difference in the ESC behavior is that the errors in $G_{E}^{V_{i}}$ are 3–10 times larger than in $G_{E}^{V_{4}}$ (see data in Tables 2 and 3). Also, since one cannot extract a value for $G_{E}(Q^{2}=0)$ using the operators $V_{i}$ due to kinematic constraint, the fits to $G_{E}^{V_{i}}$ versus $Q^{2}$, discussed in Sec. V.2, are less stable because they are not anchored at $Q^{2}=0$. As a result, the extraction of the electric charge radius from the $G_{E}^{V_{i}}$ data has much larger errors. Because of these two reasons, it has been common to analyze only $G_{E}^{V_{4}}(Q^{2})$. With our high-statistics data, we are able to compare the ESC, the efficacy of the $3^{\ast}$ fits, and the discretization errors between $G_{E}^{V_{4}}$ and $G_{E}^{V_{i}}$.

A comparison of results for $G_{E}^{V_{i}}$ and $G_{E}^{V_{4}}$ is presented in Fig. 1 for the thirteen calculations. As stated above, the errors in $G_{E}^{V_{i}}$ are much larger than those in $G_{E}^{V_{4}}$, however, there are two additional noteworthy patterns. First, the data for $G_{E}^{V_{i}}$ for $Q^{2}\lesssim 0.2$ GeV2 on the $a12m220$, $a09m220$, $a06m220$ and the two physical mass ensembles $a09m130W$ and $a06m135$, have the largest errors and mostly lie below those from $G_{E}^{V_{4}}$. On the other hand, the data for $Q^{2}\gtrsim 0.2$ GeV2 overlap in most cases. Our conclusion, based on these data, is that for $Q^{2}\gtrsim 0.2$ GeV2 the two measurements can be considered to have the same mean but with different variance.

The pattern of data at $Q^{2}\lesssim 0.2$ GeV2 is puzzling and we do not have an explanation for the larger errors or the systematic differences. In particular, we cannot discern whether they are due to residual ESC, statistical fluctuations and/or different discretization errors. In summary, while our high-statistics data have allowed us to quantify the larger errors and fluctuations in $G_{E}^{V_{i}}$, we do not have a resolution for the difference. Operationally, using a weighted average of the nonzero $Q^{2}$ data from $G_{E}^{V_{i}}$ and $G_{E}^{V_{4}}$, i.e., assuming that the differences are statistical fluctuations, gives results that are essentially identical to those from $G_{E}^{V_{4}}$. We, therefore, analyze only the data from $G_{E}^{V_{4}}$ in the rest of this paper. To establish full control over all systematics, future calculations should demonstrate consistency between $G_{E}^{V_{i}}$ and $G_{E}^{V_{4}}$.

IV.2 Extraction of $G_{M}(Q^{2})$

Examples of the size and shape of the ESC in the extraction of $G_{M}^{V_{i}}$ are shown in Figs. 28 and 29. For small momentum transfer, the convergence is monotonic from below as shown in Fig. 28 for $n^{2}=1$. The ESC is observed to grow with decreasing $a$ and $M_{\pi}$.

The pattern of convergence changes with $Q^{2}$: for small $n^{2}$ it is from below but by about $n^{2}=6$, it has changed to from above in most cases as illustrated in Fig. 29. As a result, removing ESC increases the value of $G_{M}(Q^{2})$ at small momentum transfers and decreases it at larger momenta. Consequently, if ESC is not removed, both the magnetic charge radius and the magnetic moment extracted are underestimated.

The results of the $3^{\ast}$-fits to the data for the bare form factor $G_{M}(Q^{2})$ are summarized in Table 4. A key shortcoming of the analysis of the lattice $G_{M}(Q^{2})$ is the lack of data at $Q^{2}=0$. To overcome this, we note that the ratio $G_{M}/G_{E}$, shown in Fig. 12, is, within errors, linear in $Q^{2}$ for $Q^{2}\lesssim 0.6$ GeV2. We, therefore make a linear fit to the ratio of the form factor data, $G_{M}(Q^{2})/G_{E}(Q^{2})$, with momenta up to $\vec{n}=(2,1,1)$ to obtain an estimate for the renormalized $G_{M}(Q^{2}=0)$. The corresponding unrenormalized values, which we call derived $G_{M}(Q^{2}=0)$, are also given in Table 4. These values are indistinguishable from those obtained from taking a ratio of the two correlators and then making a linear fit versus $Q^{2}$ to these data. Including these values of $G_{M}(Q^{2}=0)$ improved the stability of the $z$-expansion fits. Note that, the extrapolation of $G_{M}/G_{E}$, inclusion of the extrapolated value of $G_{M}(0)$, and the fit to $G_{M}$ are done within a single jackknife loop, therefore, the statistical errors are accounted for correctly.

To estimate the importance of using the derived point $G_{M}(Q^{2}=0)$, which anchors the fits to data, especially on ensembles with largish values of the minimum $Q^{2}$, we performed the following test. We fit the nonzero $Q^{2}$ data for $G_{E}^{V_{4}}$ to extract the value and the slope at $Q^{2}=0$ for each ensemble. Comparing the value for $g_{V}$ from this fit with the data given in Table 2, we find the magnitude of the difference for the dipole and $z^{4}$ fit is between 0.01–0.04 for the 13 calculations. The difference in the slope, $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, compared to the data in Table 5 is up to 9% for the dipole fit and up to 20% for the $z^{4}$ fit. Based on this test, it is not unreasonable that an uncertainty of similar size can be present in the extraction of $\mu$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$. Thus, to get high precision results without resorting to a derived value for $G_{M}(0)$ or without using priors, requires having data at smaller values of $Q^{2}$.

IV.3 Dependence of $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ on the lattice parameters

In Figs. 2–6, we explore the dependence of the renormalized form factors $G_{E}^{V_{4}}(Q^{2})/g_{V}$ and $G_{M}^{V_{i}}(Q^{2})/g_{V}$, which we henceforth label $G_{E}(Q^{2})/g_{V}$ and $G_{M}(Q^{2})/g_{V}$ for brevity, as a function of the pion mass, lattice spacing, lattice volume and the smearing size. The significant features are:

• •

The dependence of $G_{E}(Q^{2})/g_{V}$ on the pion mass, keeping the lattice spacing roughly constant, is shown in Fig. 2. The data show a steeper fall off as the quark mass is lowered. The behavior of $G_{M}(Q^{2})/g_{V}$ is similar as shown in Fig. 4.

• •

The data for $G_{E}(Q^{2})/g_{V}$ do not show any significant dependence on the lattice spacing $a$ for fixed pion mass as shown in Fig. 3. A similar insensitivity to change in $a$ is exhibited by $G_{M}(Q^{2})/g_{V}$ as shown in Fig. 5.

Estimates for $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ from $z$-expansion fits without including our derived value for $G_{M}(0)$ are, in many cases unstable even for the $z^{3}$ or $z^{3+4}$ fits, i.e., estimates for $\mathop{r_{M}}\nolimits$ become negative. We conclude that the fits in these cases are over-parameterized. Including the derived value of $G_{M}(0)$ and imposing the constraint on $a_{k}$ discussed in Sec. II greatly improved the $z$-expansion fits. On the other hand, the dipole fits give consistent estimates with or without using a value for $G_{M}(0)$. Our final results for both types of $Q^{2}$ fits are obtained including the $G_{M}(0)$ points.

Lastly, the comparison of the lattice data with the Kelly fit to the experimental data is shown in Figs. 3 and 5. Both $G_{E}(Q^{2})/g_{V}$ and $G_{M}(Q^{2})/g_{V}$ move towards the Kelly curve as $M_{\pi}$ and $a$ are reduced. However, $G_{E}(Q^{2})/g_{V}$ from the two physical mass ensembles still shows significant deviations from the Kelly fit. The data for $G_{M}(Q^{2})/g_{V}$ show a different curvature from the Kelly curve and points with $Q^{2}\lesssim 0.2$ GeV from the physical mass ensembles move below the Kelly curve. This change in behavior in $G_{M}(Q^{2})/g_{V}$ results in an underestimate of both $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and the magnetic moment $\mu^{p-n}$ as discussed in Sec. VI.

IV.3.1 Dependence on lattice size

Simulations on large lattices are not only important for reducing finite volume effects but also provide the simplest solution to obtaining data at smaller $Q^{2}$ for fixed $a$ and $M_{\pi}$. To demonstrate the improvement possible, we compare data from the $a12m220S$, $a12m220$ and $a12m220L$ ensembles in Fig. 30 in Appendix C. As the data move to smaller $Q^{2}$ with increasing $L$, the statistical quality of the signal also improves for a fixed number of measurements.

In Fig. 31, we show $G_{E}$ and $G_{M}$ versus $Q^{2}$ for these three ensembles. The data on the two larger volumes, $a12m220$ ($M_{\pi}L=4.38$) and $a12m220L$ ($M_{\pi}L=5.49$), overlap for both $G_{E}(Q^{2})/g_{V}$ and $G_{M}(Q^{2})/g_{V}$, indicating that finite volume effects are small for $M_{\pi}L\gtrsim 4.4$. On the smaller volume $a12m220S$ ($M_{\pi}L=3.29$), $G_{M}(Q^{2})/g_{V}$ falls off faster with $Q^{2}$.

In Fig. 32, we compare the results of three fits to $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ given in Tables 2 and 4 versus $Q^{2}$ for these three ensembles. For the $z$-expansion fits, the results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ from the two larger volumes are consistent within $1\sigma$, while those on $a12m220S$ differ. We find no significant difference in the dipole fits. These comparisons indicate that finite volume corrections are smaller than the statistical errors on the two larger volumes corresponding to $M_{\pi}L\gtrsim 4.4$. For this reason, we carry out CCFV fits including (11-point fit) and discarding the $a12m220S$ point ($10^{\ast}$-point fit). Operationally, the fits are insensitive to the $12m220S$ point due to the larger errors in it. Nevertheless, our final results, presented in Sec. VI, are from the 11-point fit.

The bottom line is that increasing $L$ for fixed $M_{\pi}$ and $a$ improves the analysis in a number of ways because the values of $Q^{2}$ for a given ${\vec{n}}$ decrease. First, the statistical errors for a fixed number of measurements decrease. The reduction in errors roughly compensates for the increase in cost of each measurement due to a larger volume. Second, with the decrease in $Q^{2}$, the ESC in $G_{E}^{V_{4}}$ becomes smaller, while that in $G_{M}$ becomes easier to control using n-state fits. Lastly, the extraction of $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ improves since the fit parameters are determined from data with values of $Q^{2}$ closer to zero.

IV.3.2 Dependence on smearing size

In Figs. 33 and 34 in Appendix C, we compare the ESC in $G_{E}^{V_{4}}$ and $G_{M}^{V_{i}}$ for two different smearing sizes using data from the $a06m310$ and $a06m220$ ensembles. The data show that the ESC is smaller with the larger smearing size.

The results of the dipole, $z^{4}$ and $z^{5+4}$ fits to $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ versus $Q^{2}$ for these two ensembles are shown in Fig. 35. Results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ are consistent within $1\sigma$ for the two smearings. The data in Fig. 6, however, show that estimates of $\mu$ can differ by about 5% between the two calculations with different smearing size. This level of difference can be explained by a combination of statistical and possible systematic uncertainties.

IV.3.3 Dependence on lattice scale setting

The two places the lattice scale enters our calculation is in converting $Q^{2}a^{2}$ to physical units and in the CCFV fits. In Table 12, we give the values of $a$ for the HISQ ensembles obtained by the MILC collaboration using the Sommer scale $r_{1}$ Bazavov et al. (2013); Sommer (2014). In Table 14 in Appendix B, we give the value of $M_{N}$ obtained on each ensemble using these values of $a$ and fit them using the leading order CCFV fit defined in Eq. (29). The result in the continuum limit is $M_{N}=976(20)$ MeV. The deviation of about 4% from the experimental value indicates a systematic uncertainty of 2–6% in the scale obtained from $r_{1}$ versus $M_{N}$, the latter analyzed using the leading order CCFV fit. The question then is, how does this difference impact the analysis of the form factors and the extraction of $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$?

The lattice data plotted in Figs. 2–5 show that the dependence of the form factors on $M_{\pi}^{2}$ and $a$ is small. To explore the dependence further, we remove the use of $a$ taken from the analysis of the Sommer scale $r_{1}$ on the HISQ ensembles by plotting the data versus $Q^{2}/M_{N}^{2}$ in Fig. 7 (bottom) where the lattice values of $M_{N}$ are used to construct the dimensionless ratio $Q^{2}/M_{N}^{2}$ for the lattice data and $M_{N}=939$ MeV for the Kelly curve. The relative movement between the data and the Kelly curve, when plotted versus $Q^{2}/M_{N}^{2}$ as compared to $Q^{2}$, brings the data closer together onto a single curve as can be seen by comparing the top and bottom set of panels. For the physical mass ensembles, the size of the relative movement of data depends only on the discretization errors, i.e., the value of $M_{N}$ at that value of $a$, assuming finite volume corrections are negligible. Presuming a cancellation of some of the systematics when the data are plotted versus $Q^{2}/M_{N}^{2}$, this comparison indicates that the observed larger deviation from the Kelly curve, when the data are plotted versus $Q^{2}$, can be explained partly as a systematic effect due to discretization errors, i.e., variations in the lattice scale set using different observables. This systematic is avoided if data at a given $Q^{2}$ are first extrapolated to $M_{\pi}=135$ MeV and $a=0$ and then compared with the Kelly curve. An attempt at doing this is described in Sec. VI.2.

It is important to note that $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ extracted for each ensemble are unchanged whether one calculates them using $Q^{2}$ or $Q^{2}/M_{N}^{2}$ as the independent variable in Eq. (10). The result would be different if the product $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits$ is calculated on each ensemble, and extrapolated to the continuum limit first, and the result divided by the experimental value for $M_{N}$. We discuss this analysis in Sec. VI.1.

Having made clear that part of the noticeable spread in the behavior of the form factors shown in Fig. 7 can be accounted for, in a large part, as due to discretization errors, the question is–what is the more robust way of analyzing the data? Should we use the scale set using $r_{1}$ or work with dimensionless variables in units of $M_{N}$? While our analysis has exposed this systematic, our conclusion is that a larger data set, or the use of a lattice action with much smaller discretization errors or a better determined extrapolation ansatz are needed to significantly reduce such systematics. Having highlighted the size of this systematic uncertainty, most of the analysis presented below is carried out versus $Q^{2}$. We provide comparison with results plotted versus $Q^{2}/M_{N}^{2}$ at appropriate places, and analyze data for $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits$ in Sec. VI.1.

V Characterizing the $Q^{2}$ behavior of the form factors

In order to extract the charge radii defined in Eq. (10) and the magnetic moment in Eq. (8), we need to parameterize the form factors versus $Q^{2}$. The two fits we explore are the dipole and the $z$-expansion truncated at some power $k$ as discussed in Sec. II. Since the dipole ansatz is the solution to an exponentially falling charge distribution (thus a model) and the $z$-expansion involves a truncation plus a constraint on the size of the coefficients $a_{k}$, it behooves us to first test these ansätze on the high-precision experimental data as discussed next.

V.1 Experimental data for the form factors and their $Q^{2}$ behavior

Electromagnetic form factors of nucleons are extracted from differential cross-sections measured in the scattering of electrons off nuclei. The process of going from measurements of the differential cross-sections to nucleon form factors is nontrivial and involves modeling Perdrisat et al. (2007); Lee et al. (2015); Yan et al. (2018). As already stated, we have two reasons to analyze the experimental data: to compare them against the lattice data over the range $0<Q^{2}\lesssim 0.8$ GeV2, and to test the efficacy of the dipole and $z$-expansion fit ansätze. For these purposes, we have collected together compiled experimental data for the proton and the neutron in Appendix D (see Figs. 36 and 37). From these, we have determined the Kelly parameterization for the isovector combinations, $G_{E}^{p}-G_{E}^{n}$ and $G_{M}^{p}-G_{M}^{n}$. Henceforth, for brevity, we will continue to use $G_{E}(Q^{2})/g_{V}$ and $G_{M}(Q^{2})/g_{V}$ to represent the $(p-n)$ combinations when comparing the lattice and the experimental data.

Next, we test the fit ansätze on the experimental data. The results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ for the proton are shown in Fig. 36. Based on the $\chi^{2}/{\rm DOF}$, the dipole fit works surprisingly well for $G_{E}(Q^{2})$, and the deviation from the data is less than a percent over the range $0<Q^{2}\leq 1$ GeV2. This difference is far less than the precision of our lattice data. For $G_{M}(Q^{2})$, the deviation is larger (up to 6%) and the $\chi^{2}/DOF$ of the fit is poor. In the the $z$-expansion fits with constraints, results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ stabilize for $k\geq 5$ as shown in Fig. 38.

Based on this analysis, and as noted in Appendix D, one should not expect a match between our lattice and the experimental data to better than about 5% or be able to resolve differences between the dipole and the $z$-expansion fits at or below this level. These comparisons provide a framework for our lattice analyses using the $z$-expansion: extract $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ from $k=4$ to avoid over-parameterization for some of the ensembles.

In the final estimates, we have assigned an additional systematic uncertainty to account for the fact that the CCFV fits have been made using just the leading order corrections. This is discussed further in Sec. VI.

V.2 Analysis of the lattice QCD data for the form factors

A comparison of the form factors $G_{E}(Q^{2})/g_{V}$ and $G_{M}(Q^{2})/g_{V}$ from all thirteen simulations with the Kelly parameterization of the experimental $(p-n)$ data is shown in Fig. 7. The data for $G_{E}(Q^{2})/g_{V}$ lie above the Kelly curve with those from the two physical mass ensembles being the closest as shown in Fig. 3, whereas the data for $G_{M}(Q^{2})/g_{V}$ lies about the Kelly curve for $Q^{2}\gtrsim 0.2$ GeV2 and then falls below it for smaller $Q^{2}$ as highlighted in Fig. 5. In both cases, these deviations from the Kelly curve impact the slope at $Q^{2}=0$, i.e., both $\mathop{r_{E}}\nolimits$ and $\mathop{r_{M}}\nolimits$ come out smaller than the phenomenological estimates. More importantly, the very precisely measured magnetic moment, $G_{M}(0)=\mu_{p}-\mu_{n}$, is underestimated by about 16%. As remarked above in Sec. IV.1 and Sec. IV.2, removing the ESC using the $3^{\ast}$-fits increases the value of all three, nevertheless, the final results presented in Sec. VI are smaller than the experimental values. Furthermore, deviations of the lattice form factors from the Kelly curve are apparent over a range of $Q^{2}$.

As discussed in Sec. IV.3.3, part of the difference between the Kelly curve and the data is due to the mismatch in the scale set by $r_{1}$ and $M_{N}$. This is highlighted in Fig. 7 where data are plotted versus $Q^{2}$ (top panels), evaluated using the lattice scale set by $r_{1}$, and versus the dimensionless variable $Q^{2}/M_{N}^{2}$ (bottom panels). Note that this change of variable does not impact the results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu^{p-n}$ on each individual ensemble and thus their extrapolated values.

The data for $G_{E}$ and $G_{M}/(\mu=4.7058)$ versus $z$ are shown in Fig. 8. Our overall strategy for extracting $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu^{p-n}$ is the following: We first determine by eye the largest value of $Q^{2}$ up to which the data are smooth in $z$. Next, since we are interested in the value and slope of the fits at $Q^{2}=0$, we restricted the data to $Q^{2}\leq 1$ GeV2, except for the $a15m310$ ($Q^{2}\leq 1.4$ GeV2) and $a12m220$ ($Q^{2}\leq 0.8$ GeV2) ensembles. The allowed range [math]–$Q^{2}|_{\rm max}$, where $Q^{2}|_{\rm max}$ is the largest value allowed by the cuts defined above, is marked by the two vertical red lines in Fig. 8. With these cuts, the points at all $Q^{2}$ are retained for most of the ensembles. Only the high $Q^{2}$ data for $a15m310$, $a12m310$, $a12m220S$ and $a12m220$ ensembles are removed. These show a break in the smooth behavior in $z$ as is clear from Fig. 8. Going back to the ESC analysis, the reliability of these points is questionable since the data have large errors and the ESC fits were poor.

The results from the $z$-expansion fits are stable for $k\geq 4$ as shown in Fig. 9. Results from fits including the sum rules are similar, except that stability is reached only for $k\geq 7$. Estimates from fits with and without sum rules are consistent, however the errors are larger with the sum rules. The values and $\chi^{2}/{\rm DOF}$ of the dipole fits have been stable under increase in statistics for all thirteen calculations. On the other hand, the results from the $z$-expansion fits required high statistics to exhibit convergence with the order of the truncation.

The results from seven fit ansatz are collected together in Tables 5, 6 and 7. Overall, the seven estimates are consistent within errors. Since $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu^{p-n}$ should be extracted from the small $Q^{2}$ behavior, our final results are from the $z^{4}$ fits. Estimates with sum rules are used only as consistency checks. The dipole, $z^{4}$ and $z^{5+4}$ fits and results are shown in Figs. 10 and 11 for ten ensembles.

The data for the ratio $\mu^{p-n}\times G_{E}(Q^{2})/G_{M}(Q^{2})$ are shown in Fig. 12. Experimental data indicate that this ratio for the proton is estimated to cross zero around $Q^{2}=8$ GeV2 Punjabi and Perdrisat (2014). Our data for the isovector combination do show a negative slope over the region $Q^{2}\lesssim 0.6$ GeV2, nevertheless, data at larger $Q^{2}$ are needed to determine if and where the ratio crosses zero. As discussed in Sec. IV.2, we have used this “linear” behavior at small $Q^{2}$ to estimate $G_{M}(0)$ from the ratio, including which helped stabilize the fits to $G_{M}(Q^{2})$. In the next section, we discuss the continuum-chiral-finite-volume (CCFV) fits used to get the physical estimates.

VI Results for $\langle r_{E}^{2}\rangle$, $\langle r_{M}^{2}\rangle$ and $\mu$

To obtain results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ in the limits $a\to 0$, $M_{\pi}\to 135$ MeV and $M_{\pi}L\to\infty$, we make a simultaneous (CCFV) fit in these three variable to the data given in Tables 5, 6 and 7. Given the spread in the lattice parameters of the 11 ensembles analyzed, we include the leading order correction term in each of the three variables, i.e., fits with four free parameters, $c_{i}^{E,M,\mu}$. The fit ansatz for the electric mean-square charge radius used is

[TABLE]

where the mass scale $\lambda$ is chosen to be $M_{\rho}=775\,\mathop{\rm MeV}\nolimits$ and the form of the chiral and FV corrections are taken from Refs. Bernard et al. (1998); Kubis and Meissner (2001); Gockeler et al. (2005). For the magnetic mean-square charge radius, we use

[TABLE]

where the leading dependence on $M_{\pi}$ is taken from Ref. Bernard et al. (1998); Kubis and Meissner (2001). Lastly, the ansatz used for the magnetic moment is

[TABLE]

where the forms of the chiral and finite volume correction terms are taken from Ref. Kubis and Meissner (2001); Beane (2004). We express all masses in units of GeV and the lattice spacing in fm.

In all three CCFV fit ansatz, Eqs. (25)–(27), heavy baryon chiral perturbation theory ($\chi$PT) has been used only to determine the form of the leading order chiral correction. For example, for $\mu$, $\chi$PT predicts the slope, $c_{3}^{\mu}$, of the linear dependence on $M_{\pi}$ as $M_{N}g_{A}^{2}/(4\pi F_{\pi}^{2})$ with $F_{\pi}=92.2$ MeV Beane and Savage (2003), however, we leave $c_{3}^{\mu}$ a free parameter. For $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$, we do not have data at enough values of $M_{\pi}$ to test the contribution of the different terms in the $\chi$PT prediction Kubis and Meissner (2001) as discussed later in this section. To avoid over parameterization of the fit we, therefore, include only the nonanalytical term in Eqs. (25) and (26). Our focus is on obtaining estimates at $M_{\pi}=135$ MeV, and this is achieved by relying on the data from the two physical mass ensembles to anchor the chiral part of the fit.

In Tables 5, 6 and 7, we also give the results of the CCFV fits for the following four combinations of the thirteen data points:

• •

13-point fit. All the thirteen calculations as considered to be independent, even though the a06m310 and a06m220 ensembles have been analyzed twice with different smearing sizes.

• •

11-point fit. We use the average of the two values for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ on the $a06m310$ and $a06m220$ ensembles as these have been analyzed twice. In this averaging, we assume maximum correlation between data

• •

10-point fit. We remove the coarsest ensemble, $a15m310$, from the eleven data points defined above.

• •

$10^{\ast}$-point fit. We remove the smallest volume ensemble, $a12m220S$, from the eleven data points defined above.

For each of these fits, we give results with (labeled extrap $c_{4}^{X}\neq 0$) and without (labeled extrap $c_{4}^{X}=0$) the finite volume correction. The values of the coefficients are given in Table 8. In the limit $a\to 0$ and $M_{\pi}L\to\infty$, only the terms proportional to $c_{1}^{X}$ and $c_{3}^{X}$ contribute. From these fits, we observe the following:

• •

Of our estimate $\mathop{\langle r_{E}^{2}\rangle}\nolimits\approx 0.59$ fm2, roughly half comes from $c_{1}^{E}$ and the other half from $c_{3}^{E}$. Compared to the experimental value $\mathop{\langle r_{E}^{2}\rangle}\nolimits\approx 0.86$ fm2 (see Eq. (30)), about $0.27$ fm2 is missing.

• •

Of $\mathop{\langle r_{M}^{2}\rangle}\nolimits\approx 0.46$ fm2, roughly 60% comes from $c_{1}^{M}$ and the rest from $c_{3}^{M}$. Compared to the experimental value $\mathop{\langle r_{M}^{2}\rangle}\nolimits\approx 0.85$ fm2 (see Eq. (30)), about $0.39$ fm2 is missing.

• •

There is a significant dependence of $\mu$ on the lattice spacing $a$. As a result, we get a low value, $\mu\approx 4$ Bohr magneton, in the continuum limit.

• •

The coefficient of the finite volume term is poorly determined, which is reflected in the larger error estimates with $c_{4}^{X}\neq 0$. In all cases, the two types of results overlap. To be conservative, we quote all final results including the finite volume term.

In Figs. 13, 14 and 15, we show the CCFV fits versus $a$, $M_{\pi}$ and $M_{\pi}L$ for three analyses: dipole, $z^{4}$ and $z^{5+4}$. In addition, we show fits versus a single variable $a$ or $M_{\pi}^{2}$ (gray bands). When the pink and gray bands are close or overlap, it means that the dominant sensitivity of the CCFV fit is with respect to the single variable of the gray band.

For $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, we also show the fit using the $\chi$PT expression given in Ref. Kubis and Meissner (2001) as a solid red line in Fig. 13. The variation with $M_{\pi}^{2}$ in the dipole, $z^{4}$ and the $z^{5+4}$ data is small over the range $135<M_{\pi}<350$ MeV, and consistent with the prediction of $\chi$PT. The singular behavior is expected to dominate for $M_{\pi}<135$ MeV. As shown in Fig. 16, over the range $135<M_{\pi}<350$ MeV, the decrease in the “log” part is partially compensated for by the increase in the “analytical” contribution as $M_{\pi}\to 0$.

The shape of the CCFV fit bands are similar for the dipole, $z^{4}$ and the $z^{5+4}$, except that the $z$-expansion data and the fits have larger errors. A visual overview of all 13 individual results for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and of the four CCFV fits is presented in Fig. 17. The variation with $a$ and $M_{\pi}$ in the thirteen individual calculations is small and somewhat smaller in the dipole than in the various $z$-expansion estimates.

For $\mathop{\langle r_{M}^{2}\rangle}\nolimits$, the variation with $a$, $M_{\pi}$ and $M_{\pi}L$ for each of the three cases, the dipole, $z^{4}$ and the $z^{5+4}$, is small as shown in Figs. 14 and 18. Presumably, the expected $1/M_{\pi}$ chiral behavior (see Eq. (26)) sets in at $M_{\pi}<135$ MeV. Again, the CCFV fits for the three cases shown in Fig. 14 are similar.

The largest variation in $\mu$ is versus $a$ as shown in Fig. 15. Again, the fit bands are similar for the dipole, $z^{4}$ and the $z^{5+4}$ data. The positive slope versus $a$ lowers the continuum limit result with respect to the experimental value $\mu|_{\rm expt}=4.7058$. The size of the difference between lattice data and experimental results suggests that discretization and other systematic errors in $G_{M}(Q^{2})$ are underestimated. From the summary of the results presented in Fig. 19, it is clear that the largest uncertainty is in the smallest volume, $a12m220S$, and the two physical mass, a09m130W and $a06m135$, points.

It is instructive to compare our data to the predictions of chiral perturbation theory shown in detail Fig. 16. The chiral expansion for the isovector $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mu\mathop{\langle r_{M}^{2}\rangle}\nolimits$, given in Ref. Kubis and Meissner (2001), are shown as the sum of three terms: those independent of $M_{\pi}$ (labeled constant), proportional to $\ln M_{\pi}^{2}/M_{N}^{2}$ (labeled log) and the remaining terms proportional to powers of $M_{\pi}^{2}$ (labeled analytical). In making these plots, the low energy constants (LEC) used are $c_{4}=3.4$ GeV*-1*, $c_{6}=4.77$, $d_{6}^{r}=0.74$, $e_{74}^{r}=1.65$, $g_{A}=1.276$ and $F_{\pi}=92$ MeV.

For $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, the “constant” term is negative ($\approx-0.49$ fm2). In the range $M_{\pi}=$ 135-350 MeV, $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ is approximately constant at $0.9$fm2: in this interval, the growth in the $\ln M_{\pi}^{2}/M_{N}^{2}$ term compensates for the decrease in the “analytical” terms. Below $M_{\pi}\approx 135$, the log term drives the rise in the sum. As shown in Fig. 13, lattice data are significantly smaller in magnitude but show a similar small variation between $M_{\pi}=135$–$350$ MeV. Because of this small variation, and having data at only three $M_{\pi}$ values, even including an additional analytical term proportional to just $M_{\pi}^{2}$ in our CCFV fits would over-parameterize the fit. Furthermore, over this range, a simple $M_{\pi}^{2}$ term would equally well mimic the sum of the log and the analytical terms. For this reason, we have included only one of the possible $M_{\pi}$ dependent terms, the log, in our CCFV fits.

In the case of $\mu\mathop{\langle r_{M}^{2}\rangle}\nolimits$ shown in Fig. 16 (right), the “log” and “analytical” terms are small and the “log” shows little variation. The dominant contribution comes from the $1/M_{\pi}$ and “constant” terms. Since the $1/M_{\pi}$ term provides the largest variation with $M_{\pi}$, we have only included it in the CCFV fit defined in Eq. (26).

In short, even though the $\chi$PT based expressions used for both $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$, given in Eqs. (25) and (26), use only the leading chiral correction term from the expressions in Ref. Kubis and Meissner (2001), the variation in our data at three values of $M_{\pi}$ between 135–315 MeV is small, and including more terms would result in over-parameterization. In this situation, having lattice data at $M_{\pi}\approx 135$ MeV, is crucial for controlling the uncertainty in the chiral fit to the lattice data. Note that the errors we quote in the CCFV fit results are comparable to those in the two physical mass points.

As is evident from the data in Tables 5, 6 and 7, and shown in Fig. 9, the $z$-expansion results without sum rules converge for $k\geq 4$. Since $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ should ideally be extracted from the small $Q^{2}$ behavior, our final results are obtained as follows: We take the $z^{4}$ result for the central value and the first error in it represents the analysis uncertainty, i.e., including the ESC, $Q^{2}$ and CCFV fits. We also quote a second systematic uncertainty to account for having used just the leading order CCFV fits. This is taken to be the largest of the following:

• •

The difference between the two values on the physical mass ensembles, $a09m130W$ and $a06m135$. The second error estimate for $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ is given by this difference.

• •

The difference between the value at $a06m135$ and the continuum value given by the CCFV fit. This gives the second error estimate for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and for $\mu^{p-n}$, which show the largest variation versus $a$.

• •

For the $z$-expansion, we also considered the difference between the $z^{3}$ ($z^{5}$) and $z^{4}$ values. These turn out to be smaller than the estimates from the previous two cases.

The final results, obtained by applying this prescription to the 11-point CCFV fit values summarized in Tables 5, 6 and 7, are given in Table 9. For completeness, we also give the results for the Dirac and Pauli radii derived from these $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ using Eq. (12) in Table 9.

The central values for $\mathop{r_{E}}\nolimits^{p-n}$, $\mathop{r_{M}}\nolimits^{p-n}$ and $\mu^{p-n}$ are about 17%, 19% and 16% smaller than the phenomenological values given in Eq. (30) and the precise experimental value in Eq. (9). Estimates from the dipole and $z$-expansion fits, given in Table 9, are consistent, however, the errors in $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ from the $z$-expansion fits are much larger, and about half the difference from the experimental/phenomenological values. The errors in the dipole fits are small compared to the difference between the lattice and the phenomenological/experimental estimates. As discussed in a number of places above, differences between the lattice and phenomenological estimates can be accounted for if a linear combination of the statistical and the various systematic errors is taken.

The encouraging results from our analysis are: (i) the data for both $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ is seen to converge towards the Kelly parameterization as $a$ and $M_{\pi}$ are decreased; (ii) the stability of the $z$-expansion fits improves with statistical precision, however constraints on the coefficients $a_{k}$ are still needed; (iii) while it is hard to test the nonanalytical chiral behavior predicted in Eqs. (25) and (26) with data at only three values of $M_{\pi}^{2}$, having data at the two physical mass ensembles anchors the CCFV fit and provides control over the uncertainty in values obtained from the fits.

A weakness of the lattice analysis is that $G_{M}(0)$ cannot be calculated directly due to kinematic constraints. We have motivated the use of a derived value of $G_{M}(0)$ to stabilize fits to $G_{M}(Q^{2})$. Looking ahead, the most significant improvement needed for extracting all three quantities with higher precision is generating data at smaller values of $Q^{2}$. This, unfortunately, requires ensembles with larger spatial volumes and/or new approaches such as a lattice formulation of the Dirac action with twisted boundary conditions de Divitiis et al. (2004); [ALPHA 01] R. Frezzotti et al. (2001). Both options are beyond the scope of this work as they require new simulations.

Two variants of the analysis presented above are described briefly next.

VI.1 Analysis of $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits$

The dimensionless quantities $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits$ are plotted in Fig. 20 versus $M_{\pi}^{2}$. For comparison, the phenomenological values for the isovector mean-square charge radii given in Eq. (30), imply $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits\approx 19.5$ and $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits\approx 17.3$. A priori, if some of the systematics cancel in the product, then one would get smaller variation with $a$ and $M_{\pi}$. The data for $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits$ in Fig. 20 show that there is significant variation with $M_{\pi}$. Lacking a well motivated fit ansatz, a reasonable option is to take the average of the values from the two physical mass ensembles. These estimates are again low: $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits|_{\rm dipole}=13.75(32)$, $M_{N}^{2}\mathop{\langle r_{E}^{2}\rangle}\nolimits|_{z^{4}}=12.16(1.22)$, $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits|_{\rm dipole}=11.34(38)$, and $M_{N}^{2}\mathop{\langle r_{M}^{2}\rangle}\nolimits|_{z^{4}}=9.40(1.50)$. Multiplying the results given in Table 9 by $M_{N}^{2}=22.7$ fm*-2* gives similar values. Given that the errors are also similar and because these estimates neglect possible $a$ dependence, we do not find this variant of the analysis as providing an obvious improvement.

VI.2 Combined $Q^{2}-$CCFV Fit

We also carried out a combined $Q^{2}-$CCFV fit to the $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ data from the 13 calculations using a product of the $z$-expansion for the $Q^{2}$ behavior and the functional forms for the CCFV ansatz given in Eqs. (25) and (26):

[TABLE]

Here $\bm{\eta}$ represents the vector of variables in the CCFV fit, and each coefficient $d_{k}$ of the $z$-expansion has a CCFV expansion of the form given in Eq. (25) or in Eq. (26). For example, for the four term CCFV ansatz given in Eqs. (25), $\bm{\eta}=(1,a,\log(M_{\pi}^{2}/\lambda^{2}),\log(M_{\pi}^{2}/\lambda^{2})e^{-M_{\pi}L})$, the combined fit has twenty parameters for the $z^{4}$ analysis. In performing these fits, we used Gaussian priors with mean 0 and width 5, in their appropriate units, for all the parameters. The resulting central values of the parameters were within this range.

The central values of the results with and without the finite volume term are consistent, however, the errors with the finite volume correction term included are about a factor of two larger. In Fig. 21, we show, for the $z^{4}$ case, the combined fits neglecting the finite volume correction term. The results of these combined fits are summarized in Table 9, and found to be consistent with those obtained by doing the $z^{4}$ and CCFV fits separately (labeled $z^{4}$-fit).

VII Comparison with Previous Work

There have been a number of lattice QCD calculations of electric and magnetic isovector form factors of the nucleon. Recent ones include LHPC’14 Green et al. (2014), Mainz’15 Capitani et al. (2015), LHPC’15 Green et al. (2015), ETMC’17 Alexandrou et al. (2017b), LHPC’17 Hasan et al. (2018), PACS’18 Ishikawa et al. (2018), PACS’18A Shintani et al. (2019) and ETMC’18 Alexandrou et al. (2018). In this work, we restrict the comparison to calculations that have presented results at or near the physical pion mass. Their lattice parameters are given in Table 10 and the data for $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ are plotted in Fig. 22. We focus on comparing the data for $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$ as these are the primary quantities calculated. Since the calculations have been done with different lattice actions, the data, even at the physical pion mass, are only expected to agree in the continuum limit. We find that, in fact, they agree remarkably well, much better than our analyses of various systematics would indicate.

All the data included in the comparison are shown in the upper two panels in Fig. 22. From the plot versus $Q^{2}$, we draw the following overall conclusions:

• •

The $G_{E}(Q^{2})$ data approach the Kelly curve from above, while $G_{M}(Q^{2})$ from below for $Q^{2}<0.2$ GeV2.

• •

No significant dependence on the number of flavors or the lattice spacing $a$ is manifest.

• •

The PACS’18A data at $Q^{2}<0.1$ GeV2, obtained using a large volume, show a qualitatively different behavior and lie closer to the Kelly curve. In this range of $Q^{2}$, the data are almost linear and highly correlated. They give a larger slope in both $G_{M}(Q^{2})$ and $G_{E}(Q^{2})$ and thus larger $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$.

• •

We can also compare data at $Q^{2}\approx 0.05$ and 0.1 GeV2 from our $a09m130W$ and $a06m135$ ensembles, from ETMC’18 Alexandrou et al. (2018), and the low error LHPC’17 Hasan et al. (2018) points with that from PACS’18A. Given the size of the statistical and systematic errors in individual data points, it is not clear if the observed small differences are significant at these two $Q^{2}$ values.

Two points become clear on plotting the data versus $Q^{2}/M_{N}^{2}$, as shown in the bottom two panels of Fig. 22. First, the collapse of all data into a single curve over the whole range $Q^{2}/M_{N}^{2}\lesssim 0.8$ GeV2 becomes even more pronounced. Second, the deviation of this common curve from the Kelly curve is smaller. Thus, not only do all our data from the 13 calculations fall on a common curve when plotted versus $Q^{2}/M_{N}^{2}$, as shown in Fig. 7, but so do data from four other collaborations using different lattice actions and volumes. A priori, a common curve would suggest that all the systematics cancel, and the apparent differences between the various calculations when the data are plotted versus $Q^{2}$ was largely a consequence of how the lattice scale is set.

Deviations from the Kelly curve are, however, significant in $G_{E}(Q^{2})$ for $Q^{2}>0.1$ GeV2. Data for $G_{M}(Q^{2})$ data undershoot for $Q^{2}<0.2$ GeV2 and are consistent with the Kelly curve above it. These differences are a 2–3$\sigma$ effect, and comparable to the size of the shift when the data are plotted versus $Q^{2}/M_{N}^{2}$ or $Q^{2}$. While an understanding of how the different systematics contribute, and whether there is one that dominates requires future more detailed calculations, we remind the reader that during the course of our analyses, we have pointed out systematics, for example due to ESC and the deteriorating signal in both the 2- and 3-point correlation functions at large $\vec{q}^{2}$, could give rise to uncertainties of this size.

We have already shown that the pattern of ESC in our data is sensitive to the value of $Q^{2}$. In particular, as discussed in Sec. IV, the ESC in correlators from which $G_{E}$ is extracted increases with momentum and the convergence is from above. On the other hand, it is large at small $\vec{q}^{2}$ in correlators from which we get $G_{M}$ and the convergence is from below. Thus, possible residual ESC could account for the observed deviation from the Kelly curve.

For $G_{E}$, there is a clear benefit to performing calculations at small $Q^{2}$. As illustrated in Figs. 24 and 25 for the physical mass ensembles, the ESC in $G_{E}(Q^{2})$ is still small for $\vec{n}^{2}=2$ corresponding to $Q^{2}\approx 0.1$ GeV2. (Note that for $\vec{q}^{2}=0$, the ESC is essentially zero as the vector charge is conserved and the local current has no $O(a)$ correction in forward matrix elements.) There is, however, an increase in the ESC with decreasing $a$ as shown in the bottom panels in Fig. 26. On the other hand, the ESC in $G_{M}(Q^{2})$ is large at small $Q^{2}$ as shown in Fig. 28, and the resulting larger errors in $G_{M}(Q^{2})$ reflect that uncertainty. In contrast, the PACS’18A calculation indicates that the ESC is removed in both form factors by using a tuned simple exponentially falling smearing of sources for generating quark propagators compared to the excited-state pattern that results from using a gauge invariant Gaussian smearing used in our and the other four calculations summarized in Table 10. Clearly, the efficacy of the exponential source used by PACS’18A to remove essentially all ESC needs to be validated.

The collapse of the data into a single curve indicates that finite volume corrections are already small for $M_{\pi}L\geq 4$, and the main advantage of the large volume used in the PACS’18A Shintani et al. (2019) study is it gives data at low $Q^{2}$. These data for $Q^{2}<0.1$ GeV2 represent a qualitative change in the behavior of both $G_{M}(Q^{2})$ and $G_{E}(Q^{2})$ which leads to larger values for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$. Note that since the PACS’18A estimate $M_{N}=942(11)$ MeV is close to $M_{N}^{\rm phy}=939$ MeV, their data do not move with respect to the Kelly curve when plotted versus $Q^{2}/M_{N}^{2}$ or $Q^{2}$. The authors attribute the much smaller errors, compared to the much higher statistics PACS’18 Ishikawa et al. (2018) calculation, to the use of the all-mode-averaging method and to a better tuned smearing ansatz (exponential) for the quark sources used to calculate the quark propagators. Since the advantage of simulations on large volume lattices to get data at low $Q^{2}$, and thus reliable estimates for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$, is obvious, it is important to validate the relatively low statistics PACS’18A calculation.

VIII Conclusions

We have presented calculations of the isovector electric and magnetic form factors, $G_{E}^{p-n}$ and $G_{M}^{p-n}$, using thirteen calculations on eleven ensembles of 2+1+1-flavors of HISQ Follana et al. (2007) fermions generated by the MILC collaboration Bazavov et al. (2013). These ensembles are at four lattice spacings, $a\approx 0.06$, $0.09$, $0.12$ and $0.15$ fm, three values of pion masses, $M_{\pi}\approx 135$, $220$ and 310 MeV, and the lattice size covers the range 3.3\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }M_{\pi}L\,\raisebox{-3.01385pt}{\stackrel{{\scriptstyle\textstyle<}}{{\sim}} }5.5. Each of these ensembles have been analyzed using $O(10^{5})$ measurements using the truncated solver method with bias correction. Using these high-statistics data we demonstrate control over excited-state contamination and perform a simultaneous fit in lattice spacing $a$, pion mass $M_{\pi}$ and lattice size $M_{\pi}L$ to get results at the physical point that can be compared with experimental values.

Our work constitutes three improvements:

• •

The much higher statistics allowed us to understand and control ESC better by keeping three states in the spectral decomposition of the 3-point correlation functions.

• •

Calculations at multiple values of $a$ and $M_{\pi}$ show that the variations in the data versus these two parameters is small for $Q^{2}\gtrsim 0.1$ GeV2 as illustrated in Figs. 2, 3, 4, 5 and in Figs. 7.

• •

We have presented first results with a CCFV fit to control the lattice artifacts due to discretization, chiral and finite volume effects. The data for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and the CCFV fits in Figs. 13 and 14 show the variation versus $M_{\pi}$ is small and consistent with the predictions of chiral perturbation theory Kubis and Meissner (2001) as shown in Fig. 16. In the $\chi$PT prediction, the nonanalytical term in $M_{\pi}$, included in Eqs. (25) and (26), becomes significant only for $M_{\pi}<135$ MeV, whereas over the range $350>M_{\pi}>135$ MeV, its growth is compensated for by the decrease in the analytical corrections. With such competing contributions in $M_{\pi}$, the data on the two physical pion mass ensembles at $a\approx 0.09$ and $0.06$ fm play a significant role in controlling the uncertainty. The CCFV fit in Fig. 15 shows a significant $a$ dependence in $\mu$ that leads to an underestimate by $\sim 16\%$.

Our final results for the mean-square charge radii, $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ (or equivalently the Dirac, $\mathop{\langle r_{1}^{2}\rangle}\nolimits$, and Pauli, $\mathop{\langle r_{2}^{2}\rangle}\nolimits$, radii derived from them), and the magnetic moment $\mu$ are given in Table 9. Using the dipole ansatz and the $z$-expansion to fit the $Q^{2}$ dependence give consistent results, however, the combined errors in the latter approach are about 2–3 times larger. The central values for $\mathop{\langle r_{E}^{2}\rangle}\nolimits$, $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ and $\mu$ are, about 17%, 19% and 16%, respectively, smaller than the phenomenological values given in Eq. (30) and the precise experimental value in Eq. (9). The trend in the data for the form factors, however, is towards the experimental values as the $Q^{2}$, lattice spacing and the light quark mass are decreased.

With higher precision data, the major improvement observed has been in the $z$-expansion estimates. Including constraints on the fit parameters, $|a_{k}|\lesssim 5$, the $z$-expansion fits for different truncations became more consistent. Based on an analysis of the experimental data with the same fit ansätze and evaluation of various systematics in the lattice calculations, the extraction of charge radii and magnetic moment could have $O(10\%)$ errors due to the modeling of the $Q^{2}$ behavior. Errors of similar size could also be due to statistics and ESC fits. Keeping in mind these estimates of the magnitude of possible systematics, the total uncertainty in estimates given in Tab. 9, especially for the dipole fit, are likely underestimated. Consequently, we do not consider the current deviations from the experimental values significant.

The magnitude of the systematic associated with what variable is used to set the lattice scale is exposed by plotting the data versus $Q^{2}/M_{N}^{2}$. As shown in Fig. 7 (bottom), our data from the 13 calculations fall on a common curve when plotted versus $Q^{2}/M_{N}^{2}$. In Fig. 22, we further show that both $G_{E}$ and $G_{M}$ from all lattice calculations done close to the physical pion mass also collapse onto this curve. The shift in the data when $G_{E}$ and $G_{M}$ are plotted versus $Q^{2}/M_{N}^{2}$ as compared to $Q^{2}$ is a discretization effect, i.e., the scale obtained from $M_{N}$ is different from that obtained by the various collaborations using the quantities shown in Table 11. This is remarkable considering that the number of quark flavors, lattice size and lattice spacing are different in the various calculations. Also, the deviation of the combined lattice data from the Kelly curve is significantly reduced.

Given our demonstration in Sec. VII that a large part of the difference between data obtained by various collaborations is an artifact of scale setting, the major advantage of the PACS’18A Shintani et al. (2019) calculation is that the large volume provides data at $Q^{2}<0.1$ GeV2; the agreement between data from different collaborations presented in Sec. VII indicates that finite volume corrections are already small for $M_{\pi}L\approx 4$. The PACS’18A data show no movement with respect to the Kelly curve because the estimate of the nucleon mass is consistent with the physical value. This may be because the lattice scale is set using the Omega baryon mass, $M_{\Omega}$, which is likely correlated with $M_{N}$, rather than indicating that discretization errors are already small at $a\approx 0.09$ fm. It is important to validate their data at $Q^{2}<0.1$ GeV2 in future calculations and confirm that the resulting estimates of $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ are consistent with the experimental values.

To conclude, our analysis highlights three points. First, all lattice data are remarkably consistent and the form factors show little dependence on the number of flavors, lattice spacing, quark mass or the lattice volume, at least for data with $M_{\pi}\lesssim 300$ MeV and $M_{\pi}L\gtrsim 4$. Second, the size of the remaining deviations in $G_{E}$ and $G_{M}$ between the lattice data and the Kelly curve are consistent with the various quantifiable systematics such as excited-state contamination and deteriorating signal at large $\vec{q}$. Third, current results provide confidence that there are no hidden systematics that afflict the calculations of form factors on the lattice.

With the lattice methodology in place, improved estimates for form factors will be obtained in future high-statistics calculations that provide data at $Q^{2}<0.1$ GeV2 and use nucleon interpolating operators that have smaller excited-state contamination.

Appendix A Lattice parameters

In this Appendix, we summarize, in Table 12, the parameters of the eleven ensembles used in the calculation. Two ensembles, $a06m310$ and $a06m220$ have been analyzed twice with different smearing sizes as listed in Table 13, where we give the parameters used in the generation of the clover propagators. These two tables are essentially the same as in Ref. Gupta et al. (2018), and have been reproduced here to keep the discussion self-contained.

Appendix B Nucleon Mass

The masses of the nucleon ground and three excited states given by our 4-state fit are summarized in Table 14. The ground state masses are found to be stable under changes in the number of states kept in the spectral decomposition of the two-point function and the Euclidean time interval used in the fits since the data exhibit a reasonable plateau in the effective mass plot for all the ensembles. On the other hand, the excited-state energies are sensitive to the details of the fits. The main reason is the small number, 6–10, of points at short times that are available to determine the six excited-state parameters before the ground state dominates the two-point function. This is particularly true of the $a15m310$, $a09m310$ and $a06m310W$ ensembles. Overall, the estimates for the excited-state masses are larger than values expected based on phenomenological arguments. For example, the first excited-state mass for the “Roper’, and the $N\pi$ and the $N\pi\pi$ multiparticle states for our physical mass ensembles, should all be between 1.3–1.7 GeV for our lattice parameters. Having a reliable estimate of the first excited-state energy is the key variable in the 3∗-fits to control ESC.

In the fits to the nucleon two-point function, we find a strong correlation between the excited-state energies and the amplitudes. This poses a challenge: what priors to choose in the 3- and 4-state fits, especially when there are near flat directions in the parameter space. We chose priors with a large width and aimed for stable first excited-state energy and amplitude. These are the most important input for the analysis of the ESC as the $0\leftrightarrow 1$ transition matrix elements are found to be the dominant artifact. Note that the priors are used only to stabilize the fits and the errors are given by the jackknife procedure.

The two physical mass ensembles give estimates for $M_{N}$ that are about 13 MeV larger than the physical value $M_{N}^{\rm phy}=939$ MeV. To investigate the dependence of $M_{N}$ on the lattice spacing, pion mass and lattice size, we have carried out two fits:

[TABLE]

where the second relation enforces $M_{N}^{\text{phys}}=939$ MeV at $M_{\pi}=135$ MeV. The values of $a$ for the HISQ ensembles used to convert the lattice data to GeV are taken from Ref. Bazavov et al. (2013) and given in Table 12. The fits to the ground-state nucleon mass $M_{0}$, given in Table 14, using Eqs. (29) are shown in Fig. 23, and the values of the fit parameters are given in Table 15. Note that the $\chi$PT predicted value for the coefficient $c_{4}=3g_{A}^{2}/(32\pi F_{\pi}^{2})=-5.716$ using $g_{A}=1.276$ and $F_{\pi}=92.2\,\mathop{\rm MeV}\nolimits$, whereas both fits give smaller values.

These fits imply that the first, unconstrained, CCFV fit to $M_{N}$ requires higher order correction terms, but with just three values of the pion mass and four values of $a$, most of the fits parameters $c_{i}$ are already poorly determined as shown in Table 15. While the results for $M_{0}$ from the two physical mass ensembles are about $13$ MeV larger than the physical value, the unconstrained fit gives an even larger value $M_{N}=976(20)$ MeV. It is clear that better control over systematics via calculations on a larger number of ensembles is needed in future calculations. The impact of the resulting mismatch between the scales set using $r_{1}$ calculated on the HISQ ensembles and from $M_{N}$ calculated using the Wilson-clover fermions, on the form factors is shown in Fig. 7 and discussed in Sec. V.2.

Appendix C ESC in the extraction of the form factors

In this Appendix, we show the data and the $3^{\ast}$-state fits used to control the ESC in the extraction of the electric and magnetic form factors. There are three sets of figures:

• •

The comparison of the ESC on the various ensembles and at different values of the momenta are shown in Figs. 24–29.

• •

The improvement in the quality of the signal with increase in the lattice size $L$ is shown in Fig. 30 using data from for the three ensembles $a12m220L$, $a12m220$ and $a12m220S$. A study of finite size effects in the form factors using these three ensembles are examined in Fig. 31, and in the extraction of $\mathop{\langle r_{E}^{2}\rangle}\nolimits$ and $\mathop{\langle r_{M}^{2}\rangle}\nolimits$ using the dipole, $z^{4}$ and $z^{5+4}$ fits in Fig. 32.

• •

A comparison of the ESC with two different smearing sizes is shown in Figs. 33 and 34 for the $a06m310$ and $a06m220$ ensembles, respectively.

Appendix D Summary of Experimental form factors

In this Appendix, we collect in one place the experimental data for the form factors for the proton and the neutron. In Fig. 36, we show the data for $G_{E}^{p}(Q^{2})$ and $G_{M}^{p}(Q^{2})$ compiled by Douglas Higinbotham Yan et al. (2018); Alarcón et al. (2019); Higinbotham and McClellan (2019) from the cross sections provided in the Lee-Arlington-Hill supplemental material Lee et al. (2015), who rebinned the original data obtained by the A1 Collaboration using the MAMI beam at Mainz Bernauer et al. (2010, 2014). The neutron data, $G_{E}^{n}(Q^{2})$, are collected from Refs. Gentile and Crawford (2011); Schiavilla and Sick (2001); Sulkosky et al. (2017), and $G_{M}^{n}(Q^{2})$ from Refs. Anderson et al. (2007); Anklin et al. (1994, 1998); Bartel et al. (1973, 1972); Bermuth et al. (2003); Bruins et al. (1995); Budnitz et al. (1968); Eden et al. (1994); Esaulov et al. (1987); Gao et al. (1994); Golak et al. (2001); Glazier et al. (2005); Hanson et al. (1973); Herberg et al. (1999); Kubon et al. (2002); Madey et al. (2003); Markowitz et al. (1993); Meyerhoff et al. (1994); Ostrick et al. (1999); Passchier et al. (1999); Plaster et al. (2006); Rohe et al. (1999); Stein et al. (1966); Warren et al. (2004); Zhu et al. (2001). These are shown in Fig. 37. From these data, we evaluate the isovector form factors $G_{E}^{p-n}(Q^{2})$ and $G_{M}^{p-n}(Q^{2})$ to which our lattice data are compared.

The construction of the isovector form factors is done as follows: we first fit the four sets of experimental data for $Q^{2}\lesssim 1$ GeV2 using the Kelly parameterization as shown in Figs. 36 and 37. Using the resulting Kelly fits, we then construct the isovector combinations, $G_{E}^{p}-G_{E}^{n}$ and $G_{M}^{p}-G_{M}^{n}$. This parameterization is used throughout the paper to compare the lattice data against. From this procedure we get

[TABLE]

whereas, using the parameter values given in the original Kelly fit Kelly (2004) gives

[TABLE]

Lattice results for the isovector combination of the radii should be compared to the values given in Eq. (30). Note that our less sophisticated analysis, which is also used to analyze the lattice data, gives larger errors.

The results of the dipole fits to the proton data shown in Fig. 36 give $r_{E}^{p}\sim 0.833$ and $r_{M}^{p}\sim 0.795$, which are roughly consistent with the careful analysis of electron-experiment results Yan et al. (2018) given in Eqs. (13) and (14) and the Kelly fits shown in the bottom row of Fig. 36. Overall, the dipole ansatz does a good job of fitting the experimental $G_{E}(Q^{2})$ data, and the deviation is less than 1% for $Q^{2}<1$ GeV2. The dipole fit to $G_{M}(Q^{2})$ is less good as shown by the large $\chi^{2}/DOF$.

The convergence of the $z$-expansion fits, with constraints on the $a_{k}$, versus $k$ is shown in Fig. 38. Estimates with $k\geq 5$ are stable for all three quantities.The results from $z$-expansion fits, also shown in Fig. 36, are marginally larger than those from the dipole and differ by a few percent from those in Eqs. (13) and (14). The errors in the $z$-expansion estimates are larger, especially with the inclusion of the sum rules. The overall lesson from this exercise is that an uncertainty of $O(5\%)$ could be present in our analysis using either the dipole or the $z$-expansion fits.

Acknowledgements.

We thank the MILC Collaboration for providing the 2+1+1-flavor HISQ lattices used in our calculations. The calculations used the Chroma software suite Edwards and Joo (2005). R. Gupta thanks D. Higinbotham for discussions and for providing the experimental data on the form factors. Simulations were carried out on computer facilities of (i) the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231; and, (ii) the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725; (iii) the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy, and (iv) Institutional Computing at Los Alamos National Laboratory. T. Bhattacharya and R. Gupta were partly supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics under Contract No. DE-AC52-06NA25396. T. Bhattacharya, R. Gupta, J.-C. Jang and B.Yoon were partly supported by the LANL LDRD program. The work of H.-W. Lin is supported by the US National Science Foundation under grant PHY 1653405 “CAREER: Constraining Parton Distribution Functions for New-Physics Searches”.