Lattice Calculation of Pion Form Factor with Overlap Fermions

We present a precise calculation of the pion form factor using overlap fermions on seven ensembles of 2+1-flavor domain-wall configurations with pion masses varying from 139 to 340 MeV. Taking advantage of the fast Fourier transform and other techniques to access many combinations of source and sink momenta, we find the pion mean square charge radius to be $\langle {r_\pi^2} \rangle= 0.430(5)(13)\ {\rm{fm^2}}$, which agrees well with the experimental result, and includes the systematic uncertainties from chiral extrapolation, lattice spacing and finite-volume dependence. We also find that $\langle {r_\pi^2} \rangle$ depends on both the valence and sea quark masses strongly and predict the pion form factor up to $Q^2 = 1.0 \ {\rm{GeV^2}}$ which agrees with experiments very well.

Since lattice QCD is an ab initio calculation and the experimental determination of r 2 π from the πe scatter- * genwang27@uky.edu ing is very precise, it provides a stringent test for lattice QCD calculations to demonstrate complete control over the statistical and systematic errors in estimates of the relevant pionic matrix element in order to enhance confidence in their reliability to calculate other hadronic matrix elements where further technical complications occur. Over the years, the pion form factor has been calculated with the quenched approximation [11,12], and for the N f = 2 [13][14][15][16][17], N f = 2 + 1 [18][19][20][21][22][23] and N f = 2 + 1 + 1 [24] cases.
In this work, we use valence overlap fermions to calculate the pion form factor on seven ensembles of domain-wall fermion configurations with different sea pion masses, including three at the physical pion mass, four lattice spacings and different volumes to control the systematic errors. Due to the multi-mass algorithm available for overlap fermions, we can effectively calculate several valence quark masses on each ensemble [25][26][27][28] and also O(100) combinations of the initial and final pion momenta with little overhead with the use of the fast Fourier transform (FFT) algorithm [29] in the three-point function contraction. This allows us to study both the sea and the valence quark mass dependence of r 2 π in terms of partially quenched chiral perturbative theory, besides giving an accurate result at the physical pion mass. This work is based on Ref. [30] with more statistics on the ensembles at the physical pion mass.
The paper is organized as follows: In Section II, we present the numerical details of this calculation and a brief description of the FFT on stochastic-sandwich arXiv:2006.05431v4 [hep-ph] 31 Aug 2021 method. Fits and extrapolations are discussed in Section III with results compared with other studies. A brief summary is given in Sec. IV.

II. NUMERICAL DETAILS
We use overlap fermions on seven ensembles of HYP smeared 2+1-flavor domain-wall fermion configurations with Iwasaki gauge action (labeled with I) [31,32] and Iwasaki with Dislocation Suppressing Determinant Ratio (DSDR) gauge action (labeled with ID) [33,34] as listed in Table I. The effective quark propagator of the massive overlap fermions is the inverse of the operator (D c + m) [35,36], where D c is chiral, i.e., {D c , γ 5 } = 0 [37]. It can be expressed in terms of the overlap Dirac operator D ov as D c = ρD ov /(1 − D ov /2), with ρ = −(1/(2κ) − 4) and κ = 0.2. A multi-mass inverter is used to calculate the propagators with 2 to 6 valence pion masses varying from the unitary point to ∼ 390 MeV. On 24I, 32I and 24IDc (c stands for coarse lattice spacing), Gaussian smearing [38] is applied with root mean square (RMS) radii 0.49 fm, 0.49 fm and 0.53 fm, respectively, for both source and sink. On 48I, 32ID and 32IDh (h for heavier pion mass), box-smearing [39,40] with box half sizes 0.57 fm, 1.0 fm and 1.0 fm, respectively, is applied as an economical substitute for Gaussian smearing.
To extract pionic matrix elements, the three-point where χ π + ( x, t) =d( x, t)γ 5 u( x, t) is the interpolating field of the pion with u and d the up and down quark spinors, S(y|x) is the quark propagator from x to y, z ≡ {τ, z}, x f ≡ {t f , x f }, p i and p f are the initial and final momenta of the pion, respectively, q = p f − p i is the momentum transfer, and G is the smeared Z 3 -noise grid source [41]. The disconnected insertions in Eq.(3) vanish in the ensemble average [12]. In practice, S(G|z) in Eq. (3) is calculated using γ 5 hermiticity, i.e., S(G|z) = γ 5 S † (z|G)γ 5 , and S(z|x f ) is usually obtained in the sequential source method with γ 5 S(x f |G) as the source [42,43]. The calculation of the sequential propagators would need to be repeated for different p f and different quark mass m, so that the cost would be very high when dozens of momenta and multiple quark masses are calculated. Instead, we use the stochastic-sandwich method [44,45], but without lowmode substitution (LMS) for S(x f |G) since it is not efficient for pseudoscalar mesons [25]. However, the separation of sink position x f and current position z in splitting the low and high modes for the propagator S(z|x f ) between the current and sink can facilitate FFT along with LMS which is still useful here. More specifically, the propagator from the sink at x f to the current at z, S(z|x f ), can be split into the exact low-mode part based on the low lying overlap eigenvalues λ i and eigenvectors v i of the ith eigenmode of D c , plus the noise-source estimate S H noi of the high-mode part, where λ c is the highest eigenvalue in LMS and is much larger than the quark mass m with the typical number of eigenmodes n v ∼ 400 on 24I and 32I, and n v ∼ 1800 on 32ID, 32IDh, 24IDc, 32IDc and 48I; and S H noi (z, η j ) is the noise-estimated propagator for the high modes with the low-mode deflated Z 3 noise η j (x f ) [44,45]. Sink smearing is applied on all the sink spatial points x f of noise η j (x f ) and eigenvectors v † i (x f ). Thus C 3pt can be decomposed into factorized forms within the sums of the eigenmodes for the low modes and the n f number of noises η j for the high modes, where which are calculated by using FFTs on the spatial points z and x f for each of G L i , F L i , G H j and F H j to obtain any q and p f with the computational complexity O(V logV ), with V the lattice spatial volume. Compared with the stochastic-sandwich method for a fixed p f which also includes the summation over the spatial points z and x f , eigenvectors v i and noises η j , the additional cost factor of using FFTs, namely O(logV ), is only of order ∼ 7 for our largest 48I lattice. This allows us to calculate any combination of q and p f without much additional cost compared to the traditional stochastic-sandwich method; this is of order ∼ 10 times less expensive if we calculate more than seven different sink momenta p f and average over different directions. In practice, since larger p i or p f will lead to worse signals, we choose three cases so that for a given Q 2 we use as small p f and p i as possible: (1) p i = 0 with q = p f or p f = 0 with q = − p i which probes small Q 2 ; (2) p f = − p i with q = 2 p f which probes reasonably high Q 2 ; (3) For a given q, we calculate q/2 and choose lattice momenta p f and − p i which are close to q/2. This can also probe high Q 2 which fills Q 2 between the previous two cases.
We use the lattice dispersion relationÊ 2 =m 2 + ip 2 i with aÊ = 2sinh(aE/2), am = 2sinh(am/2) and ap i = 2sin(ap i /2) to define Q 2 for all ensembles so that there is a well defined physical limit. This is also used in Ref. [23] to calculate the pion charge radius. More details about checking the dispersion relation are included in Appendix B.

A. Three-point function fit
The source-sink separations t f used in this work with different ensembles are collected in Table II. The largest t f is ∼ 2.0 fm on the coarsest lattice 24IDc and the smallest one is ∼ 0.7 fm on the finest lattice 32I.
With the use of Wick contractions and gauge invariance, the three-point function (3pt) with two coherent sources (we have put a source at each of t = 0 and t = T /2 for most ensembles to increase statistics) has contributions from the three diagrams shown in Fig. 1.
which includes the first excited-state contamination, where Z p is the spectral weight and E and E 1 are the ground state and first excited state energies, respectively.
i and E 1 f are constrained by the joint fit with the corresponding two-point function (2pt). Z V is the finite normalization constant for the local vector current and is determined from the forward matrix element as Z V ≡ 2E π(p)|V4|π(p) . C 1 , C 2 and C 3 are free parameters for the excited-state contamination. The diagram 1.(2) contributes in which we have ignored the excited-state contamination from the source at T /2 since such terms are suppressed by e −E 1 i T /2 which is of order ∼ 10 −8 with E 1 i ≈ 1.3 GeV estimated with the experimental value of the first excited state of the pion, and the diagram 1. (13) in which this term corresponds to the creation of a hadron state with operator V 4 =qγ 4 q at time slice τ with momentum q as h(q)| V 4 |0 , an annihilation of a pion state at time slice T /2 with momentum p i as 0| χ † π + |π − (p i ) and an unknown matrix element π − (p i )| χ π + |h(q) . The excited-state contamination from E 1 i is ignored for the same reason as in the previous discussion and the excitedstate contamination from E 1 h is ignored under current statistics.
In order to test the functional form of C 3pt,(3) (τ, t f , p i , p f ), we construct 3pt with one source at time slice T /2 = 32 and sink time t f at 20, 21, 22 with in which E eff i is evaluated by a simultaneous change of τ and t f to single out E i from the exponential (1) (2) The correlation function is a rising exponential which confirms that E h > 0 in Eq. (13). The plots in the middle and right panels show the corresponding effective masses E eff h and E eff i , respectively, obtained with Eq. (14).

Thus the final functional form is
The associated 2pt is fitted with with A 1 being a free parameter for the excited-state contributions and the exponential terms with T /2 account for contributions from the source at T /2. An example of fitted energies is shown in Fig. 3. It can be seen that the first excited state energy E 1 is close to the experimental value 1.3 GeV and it has been used to constrain that of the 3pt by the joint fit of 2pt and 3pt to extract f ππ (Q 2 ). For the special | p i | = | p f | case, one can simply calculate the ratio of 3pts, and obtain the pion form factor by the following parametrization of the ratio R 1 , where the terms with B 1 and B 2 are the contributions from the excited-state contamination, and ∆E = E 1 ( p i ) − E( p i ) is the energy difference between the pion energy E( p i ) and that of the first excited state E 1 ( p i ). These energies are also constrained by the joint fit with the corresponding 2pt. Since the excited-state contamination of the forward matrix element in the denominator is known to be small and the contribution from C 4 term in Eq. (15) is suppressed by e −E( pi)T /2 with p i = 0 for both the denominator and numerator, we have dropped them in the parametrization of the ratio and our fits can describe the data with χ 2 /d.o.f. ∼ 1. Fig. 4 shows a sample plot for 32ID with the unitary pion mass of 174 MeV at Q 2 = 0.146 GeV 2 . In view of the fact that the data points are symmetric about τ = t f /2, within uncertainty, it reassures that the sink smearing implemented under the FFT contraction has the same overlap with the pion state as that of the source smearing. In order to test the fit function of 3pt in Eq. (15), a comparison of the fit of the one-source result with the source at t = 0 and that of the two-source result with a source at each of t = 0 and 32 in the same inversion is shown in Fig. 5. For illustrative purpose, the data points on the left and middle panels are shown with ratio R 2 , in which Z p and E are determined from the fit of 2pt and Z V from 3pt at zero momentum transfer. The data for the top panels use n i = 4 with {n s , n s , n s , n t } = {2, 2, 2, 1} and the data for the lower panels use n i = 2 with {n s , n s , n s , n t } = {2, 2, 2, 2} so that their statistics are matched. The case with one source and p i = 0 and q = p f is shown in the top left panel and the gray band is close to the data points due to small excited-state contamination. The similar case with two coherent sources is shown in the lower left panel and the gray band is far away from the rising data points due to the additional C 4 term with fitted E h = 820(110) MeV, which is consistent with the result of Fig. 2. The two results agree with each other within uncertainty which confirms our fit formula, but a comparison of the statistical errors reveals that the factor of two lower cost from using two coherent sources versus one source produces no net benefit for this case of p i = 0. Since the contribution from the C 4 term is suppressed significantly for 3pts with p i = 0, the data points and results from one source and two coherent sources agree with each other very well for the cases with p f = 0, q = − p i and p f = − p i , q = 2 p f which are shown in the middle, and right panels, respectively. In these kinematical cases, however, the statistical errors are the same for one source versus two coherent sources and thus the full factor of two lower cost (in computation and storage) for the latter is fully realized.
Thus for the general momentum setup which corresponds to the exchange of initial and final momenta. Fig. 6 shows an example plot on 32ID. The data points are fitted well ( (15) and the fit results are shown with colored bands. The data points for C 3pt (τ, t f , p, 0) are lower and closer to the gray band since the C 4 term has a negative contribution with a suppression factor e −E( p)T /2 compared to the case of C 3pt (τ, t f , 0, p) in which the C 4 term has a positive and large contribution with only a suppression factor e −E( 0)T /2 .

B. z-expansion fit
To obtain f ππ (Q 2 ), we have done a model-independent z-expansion [46] fit using the following equation with k max ≥ 3. where corresponds to the two-pion production threshold, with m π,mix the mass of the mixed valence and sea pseudoscalar meson calculated in Ref [47,48] on each ensemble directly with one valence domain wall propagator and one valence overlap propagator for each valence quark mass; and t 0 is chosen to be its "optimal" value t opt 0 (Q 2 max ) = t cut (1 − 1 + Q 2 max /t cut ) to minimize the maximum value of |z|, with Q 2 max the maximum Q 2 under consideration.
In order to remove the model dependence of the zexpansion fit, we need to take k max to be large enough A non-linear least squares fit of this analytical function with z-expansion fit at k max = 10 gives |a k /a 0 | max < 1.03, in which we used t cut = 4m 2 π,phys , t opt 0 (Q 2 max ) = t cut (1− 1 + Q 2 max /t cut ) and Q 2 max = 1.0 GeV 2 . Also, by investigating the z-expansion fits with k max = 3 without priors of our data, we find |a k /a 0 | max < 3.0. Thus we propose the use of the conservative choice of Gaussian prior [46] with |a k /a 0 | max = 5 (use "a k /a 0 = 0(5)" as a Gaussian prior for all a k , k > 1, in the fits) for the pion form factor. The z-expansion fitted pion form factors up to Q 2 ∼ 1.0 GeV 2 for the seven lattices with the same valence and sea pion mass are shown in Fig. 7 with Another way to reach higher k max and control the model dependence of fits is to use the fact that at the Q 2 → ∞ limit f ππ (Q 2 ) falls as 1/Q 2 up to logarithms [49,50]. Thus we have Q k f ππ (Q 2 ) → 0 for k = 0, 1 and follow the same argument in [46], which implies with z = 1 corresponding to the Q 2 → ∞ limit. These equations lead to the two sum rules for pion form factors We have explored this alternative results shown in Fig. 10.

C. Chiral extrapolation of pion radius
With the z-expansion fit of the form factor using Eq. (19), the charge radius of pion can be obtained through for all the valence masses of each lattice. Fig. 8 shows the results on 32ID and 32IDh as a function of valence pion mass m 2 π,v and mixed pion mass m 2 π,mix in the left panel and right panel, respectively. We see that there is a strong dependence on the valence pion masses from the data points on these two ensembles. Also, the disagreement in the left panel evinces a strong dependence on the sea pion mass. In contrast, the right panel shows an agreement of 32ID results and 32IDh results at similar m 2 π,mix which guides us to use m 2 π,mix , as proposed by Partially Quenched Chiral Perturbation Theory [51], as a basic variable for the chiral extrapolations.
The r 2 π on different lattices with different valence pion masses are plotted in Fig. 9. The following fit form as a function of m 2 π,mix is used which includes an essential divergent log term from the SU (2) NLO ChPT [10,51], where π,phys )) shows the pion mass dependence of the pion decay constant from partially quenched NLO SU (2) ChPT [52] with M 2 π ≡ m 2 π,mix , F andl 6 are free parameters for fitting, m π,phys = 139.57 MeV is the physical pion mass, L is the spatial size of the lattice, the b I/ID 2 terms reflect the lattice spacing dependence for the two sets of ensembles with different gauge actions (Iwasaki and Iwasaki plus DSDR), and the b 3 term accounts for the finite-volume effect [17,[53][54][55]. Instead of fitting both the low-energy constantsl 4 and F as free parameters, which leads to unstable fits, we usel 4 = 4.40 (28), as given by its FLAG average [56], as a prior and treat F as free parameter.
The results of the fits are shown in Fig. 9. The colored bands show our prediction based on the global fit of r 2 π with χ 2 /d.o.f. = 0.85; the inner gray band shows our prediction for the unitary case of equal pion mass in the valence and the sea in the continuum and infinite volume limits and the outer band includes the systematic uncertainties from excited-state contamination, zexpansion fit, chiral extrapolation, lattice spacing, and finite-volume dependence. Since the kaon mass only varies a little in the current pion mass range, we do not include the kaon log term in the fit. The discretization errors across the Iwasaki gauge ensembles are small, while those across the Iwasaki plus DSDR gauge ensembles are obvious; this is consistent with what was found in the previous work with the DWF valence quark on similar RBC ensembles [23]. The fit gives F π = 96.2(4.3) MeV, which is consistent with 92.2(1) MeV, andl 6 = 17.1(1.4), which is also consistent with the FLAG average [56] valuē l 6 = 15.1(1.2). The systematic uncertainties considered are listed as follows: • Fit results for the radius from different z-expansion fits using Eq. (24) are shown in Fig. 10. Since b I 2 and b 3 have no statistical significance, we use only three free parameters F ,l 6 and b ID 2 in these fits and treat the low-energy constantl 4 = 4.40 (28) appearing in F π as a prior. All the fits have good   Table III. The maximum difference between the result shown in black in Fig. 10 and those of the other fitted cases is treated as the systematic uncertainty from the z-expansion fit.
• The systematic uncertainty from the excited-state contamination is estimated by changing the fit ranges of 2pt and 3pt on 32ID with pion mass 174 MeV at the smallest momentum transfer which results in f ππ (Q 2 = 0.051 GeV 2 ) = 0.9158 (14)(13); the second error corresponds to the systematic uncertainty from excited-state contamination. This case is chosen because of its good signal to noise ratio which has the most control of the final result at close to the physical pion mass, and the smallest momentum transfer is chosen due to its largest influence on the radius. In order to estimate the systematic uncertainty of the radius from the form factor at only one small momentum transfer, we solve the VMD model in Eq. (20), with m as a free parameter. The predicted radius is r 2 π = 6.0/m 2 = 0.4190(74)(68) fm 2 . The second error 0.0068 fm 2 , which propagates from the systematic uncertainty of the form factor, is treated as the systematic uncertainty from the change of fit ranges for the extrapolated charge radius.
• We added a linear dependence term between the charge radius of the pion and the pion mass squared as b 4 M 2 π to Eq. (24) proposed by SU (2) NNLO ChPT [9] and repeated the fit with four free parameters r 2 π phys , b 1 , b ID 2 and b 4 . The coefficient b 4 is consistent with zero and the prediction changes by 0.0017 fm 2 which is treated as a chiral extrapolation systematic uncertainty.
Another source of the chiral extrapolation systematic uncertainty is the lack of a kaon log term in Eq. (24). On 24I, the valence pion masses ranging from 323 MeV to 391 MeV give a range of kaon mass from 532 MeV to 554 MeV. Thus we estimate the maximum kaon mass for the pion mass range in consideration to be M K,max = 554 MeV. With the use of SU (3) NLO ChPT [10], the systematic uncertainty from the kaon log term can be given by • We repeated the fit with four free parameters F , l 6 , b ID 2 and b I 2 which includes the discretization error from the Iwasaki gauge action and the prediction changes by 0.0052 fm 2 . With this fit, we get a difference between the fit predictions in the continuum limit with those from the smallest lattice spacing (32I) to be 0.0017 fm 2 . We combined these two as the systematic uncertainty of finite lattice spacing.
• We repeated the fit with four free parameters F ,l 6 , b ID 2 and b 4 which includes the finite-volume term and the prediction changes by 0.00019 fm 2 . With the inclusion of the finite-volume term, the difference of the predictions for 24IDc (which has the smallest m π L) and 32IDc is 0.005 fm 2 . We combined these two as the systematic uncertainty of finite-volume effects.
Thus, the final result of the mean square charge radius of the pion at the physical pion mass in the physical limit reads r 2 π = 0.4298(45) stat (66) z-exp (68) fit-range (37) χ (55) a (50) V = 0.4298(45)(126) fm 2 , with statistical error (stat) and systematic uncertainty from z-expansion fit (z-exp), fit-range dependence (fit-range), chiral extrapolation (χ), finite lattice spacing (a), and finite-volume (V). The total uncertainties at heavier pion masses are estimated from the scale of the total/statistical ratio at the physical pion mass.
ChPT expansion [9,10], in which F andl 6 are free parameters for fitting, c 1 and c 2 account for possible NNLO effects, c Since the inverse of f ππ (Q 2 ) is mainly dominated by the NLO contributions considering the vector dominace of the pion form factor, fitting the inverse helps avoid the need of too many low-energy constants from NNLO corrections [17]. The fit result is shown in Fig. 12 with the central values and correlations of the fit parameters are listed in Table IV. This fit (with χ 2 /d.o.f. = 1.0) gives r 2 π = 0.433(6) fm 2 , F π = 92.1(6.2) MeV andl 6 = 16.0(1.9), which are consistent with the above analysis. Our extrapolated result at the physical pion mass and continuum and infinite volume limits for the curve f ππ (Q 2 ) including the systematic uncertainties from excited-state contamination, NNLO corrections, chiral extrapolation, lattice spacing and finite-volume dependence, is shown and compared with experiments in Fig. 13; it goes through basically all the experimental data points up to Q 2 = 1.0 GeV 2 . Also, our results are consistent with the previous experimental analysis [5] and phenomenological prediction [57]. The following systematic uncertainties are included in the analysis: • With a variation of the fit ranges of 2pt and 3pt on 32I with pion mass 312 MeV we got the form factor at large momentum transfer f ππ (Q 2 = 0.865 GeV 2 ) = 0.4347(87)(98). Along with previous analysis on 32ID at small momentum transfer f ππ (Q 2 = 0.051 GeV 2 ) = 0.9158(14)(13), we estimate the systematic uncertainty from the excitedstate contamination to be equal to the statistical uncertainty of the fitted pion form factors for all Q 2 < 1.0 GeV 2 .
• Since the c 1 and c 2 terms are just an estimation of the possible NNLO effects, we estimate the NNLO systematic uncertainty by setting c 1 and c 2 in Eq. (26) to be zero and treat the changes as the systematic uncertainty from NNLO corrections.
• The systematic uncertainty from the lack of a kaon log term proposed by SU (3) NLO ChPT is calculated with which is the difference between using M K,max and m K,phys in the ChPT formula. This is treated as the systematic uncertainty from chiral extrapolation.
• We use the difference between the fit predictions in the continuum limit with those from the smallest lattice spacing (32I) as the systematic uncertainty of finite lattice spacing.
• The systematic uncertainty from finite-volume effects is estimated by the difference between the fit predictions for 24IDc with m π L ∼ 3.33 and 32IDc m π L ∼ 4.45 with both ensembles at the physical pion mass.  [58][59][60][61][62]. The inner gray band is the statistical error and the outer band includes the systematic uncertainties.

IV. SUMMARY
We have presented a calculation of the pion form factor using overlap fermions with a range of valence pion masses on seven RBC/UKQCD domain-wall ensembles including two which have the physical pion mass. The lattice results for r 2 π in the continuum and infinite volume limits are compiled in Fig. 11 together with that of experiment. Our globally fitted pion mean square charge radius is r 2 π = 0.430(5)(13) fm 2 , which includes systematic errors from chiral extrapolation, finite lattice spacing, finite volume, and others; it agrees with experimental value of r 2 π = 0.434(5) fm 2 within one sigma. We find that r 2 π has a strong dependence on both the valence and sea pion masses. More precisely, it depends majorly on the mass of the pion with one valence quark and one sea quark. A good fit of the chiral log term confirms that the pion radius diverges in the chiral limit. We also give the extrapolated form factor f ππ (Q 2 ), and the result agrees well with the experimental data points (up to Q 2 = 1.0 GeV 2 ).
Thus this work shows that the hadron form factor and the corresponding radius can be studied accurately and efficiently by combining LMS with the multi-mass algorithm of overlap fermions and FFT on the stochastic-sandwich method. This raises the expectation of an efficacious investigation of the form factor of the nucleon and its pion-mass dependence with relatively small overhead on multiple quark masses and momentum transfers. Note that for an accurate prediction of the charge radius and form factor with 1% overall uncertainty, calculations at smaller lattice spacing and larger source-sink separation will be essential, together with the QED and isospin breaking corrections.
We have chosen a set of evenly separated configurations for measurement from the full Monte Carlo evolutions available for each ensemble. The separations are 40, 32, 10, 8, 8, 10, 32 for 24I, 32I, 48I, 24IDc, 32IDc, 32ID, 32IDh, respectively. For the error analysis, we treat measurements from different configurations as independent and average the measurements over each configuration before analysis.
In the left panels of Fig. 16, we have plotted the integrated autocorrelation time on 24IDc and 48I, defined as, The error of C(∆) is estimated by jackknife re-sampling of the average on t and the error of the integrated autocorrelation time is estimated with simple error propagation. The plot shows the three-point functions with current position τ = t f /2 for the valence pion mass 137 MeV and 148 MeV at the smallest separation C 3pt (t f /2, t f , p i = 0, p f = 0) on 24IDc and 48I, respectively. The integrated autocorrelation times are less than 1 within uncertainty for both ensembles which confirms the independence of the measurements on each of the configurations. In the right panels of Fig. 16, we plot the central values and errors of C 3pt (t f /2, t f , p i = 0, p f = 0) on 24IDc and 48I as a function of binning size t bin . The statistical errors change very little as the bin size is increased for both ensembles which again confirms the measurements' independence. In Fig. 14, we compare the pion energy E( p) obtained from the fitting of 2pts to the lattice dispersion relation where aÊ = 2sinh(aE/2), am = 2sinh(am/2) and ap i = 2sin(ap i /2). As can be seen, the dispersion relation is well satisfied under the 1% level for the momenta considered in this paper. The gray band is the z-expansion fit result. The percentage differences between the data and corresponding fit results are also plotted with red, magenta and lime colors with the scale on the right. As a reference, the inverse of signal-to-noise ratio of the fit result is displayed as percentage with the coral region.
In Fig. 15, we have plotted the pion form factors on 48I with different p i and p f cases marked with different colors. The values for different cases overlap with each other quite well at similar Q 2 at the 1% level. This confirms that the combination of p i and p f considered in this paper are consistent with each other which will lead to a well defined physical limit.
Appendix C: Normalization of the local vector current for overlap fermions The left panel of Fig. 17 shows the determination of the normalization constant Z V on 32ID by fitting the inverse of the forward matrix element as 2E π(p)|V4|π(p) with p = 0. The data points from different source-sink separations overlap well with each other under the 0.1% level, so we have done a simple linear fit with χ 2 /d.o.f. ∼ 0.7.
As we are using overlap fermions which have exact chiral symmetry on the lattice, the axial normalization (finite renormalization) constant is equal to the local vector current normalization constant, as confirmed in [63]. The axial normalization constant on 32ID was calculated in [45] from the Ward identity: Z A = 2mq 0|P |π mπ 0|A4|π with P and A 4 the pseudo-scalar quark bilinear operator and the temporal component of the axial-vector operator, respectively. As shown in the right panel of Fig. 17, the axial normalization constant agrees well with the local vector current normalization constant used in this paper very well at the massless limit.