Pion form factor and charge radius from Lattice QCD at physical point

We present our results on the electromagnetic form factor of pion over a wide range of $Q^2$ using lattice QCD simulations with Wilson-clover valence quarks and HISQ sea quarks. We study the form factor at the physical point with a lattice spacing $a=0.076$ fm. To study the lattice spacing and quark mass effects, we also present results for 300 MeV pion at two different lattice spacings $a=0.04$ and 0.06 fm. The lattice calculations at the physical quark mass appear to agree with the experimental results. Through fits to the form factor, we estimate the charge radius of pion for physical pion mass to be $\langle r_{\pi}^2 \rangle=0.42(2)~{\rm fm}^2$.


I. INTRODUCTION
Pion is one of the most prominent strongly-interacting particles next to the nucleon since it is a Goldstone boson of QCD. For this reason, it is important to study the pion internal structure and find out if there is a connection between its internal structure and its Goldstone boson nature. This issue is particularly relevant for understanding the origin of mass generation in QCD, see e.g. discussions in Refs. [1,2].
Knowledge of internal structure of the pion is much more limited than that of the nucleon. On the partonic level, the parton distribution function (PDF) of the pion has been studied through the global analysis of the Drell-Yan production in pion-nucleon collisions and in tagged deep inelastic scattering (DIS), for recent analyses see Refs. [3,4]. Recently, there have been many efforts in lattice QCD to study the pion PDF [5][6][7][8][9][10], which have used the quasi-PDF in Large Momentum Effective Theory [11,12], the pseudo-PDF [13,14] and current-current correlator [15][16][17] (also referred to as a "good lattice cross-section") approaches, see Refs. [18][19][20][21] for recent reviews. Lattice calculations of the lowest moments of pion PDF [22][23][24][25][26][27] are also available and can be used as additional constraints in the global analysis.
Form factor, defined as with J µ being the electromagnetic current and Q 2 = −(P 2 − P 1 ) 2 , provide a different insight into pion structure, namely the charge distribution. It can be, in principle, measured in electron-pion scattering. Generalized parton distribution (GPD) combine the information contained in PDF and form factors and provide a three- * xgao@bnl.gov dimensional image of a hadron. In the case of the nucleon, the study of the GPDs is the subject of large experimental and theory efforts (see e.g. Ref. [28] for a recent review). Experimental study of the pion GPD is far more challenging and will be only possible at Electron-Ion Collider (EIC), if at all. Fortunately, GPDs can be studied on the lattice using LaMET, including pion GPDs [29][30][31][32]. Experimentally, the pion form factor was measured by scattering of pions off atomic electrons in Fermilab [33,34] and CERN [35,36]. This allowed determination of the pion form factor for momentum transfer Q 2 up to 0.253 GeV 2 [33][34][35][36]. For larger Q 2 , one has to determine the pion form factor from the electroproduction of charged pions off nucleons. The corresponding experiments have been performed in Cornell [37][38][39] DESY [40,41], and Jlab [42][43][44][45][46]. These determinations, however, were model-dependent. The recent determination of the pion form factor up to Q 2 of 2.45 GeV 2 is carried out by the F π collaboration using data both from DESY and JLab [46]. Experiments at the future EIC facility will allow us to probe even higher Q 2 up to 30 GeV 2 and possibly see the partonic structure in an exclusive elastic process and make contact with asymptotic large-Q 2 perturbative behavior [47]. In the timelike region, the pion form factor can be determined by analyzing e + e − → π + π − process [48] (see also references therein). This analysis also constrains the form factor in the spacelike region.
Lattice QCD calculations allow one to obtain the pion form factor from first principles, i.e. without any model dependence, up to relatively large Q 2 . Therefore, they will provide an important cross-check for the experimental determinations. The first lattice calculations of the pion form factor date back to late 80s and were performed in the quenched approximation [49,50]. More recently, lattice calculations of the pion form factor have been performed with two flavors (N f = 2) of dynamical quarks [51][52][53][54][55], with physical-mass strange-and two light-quark flavors (N f = 2 + 1) [56][57][58][59][60][61][62], as well as with a dynamical charm quark, a strange quark and two flavors of the light quarks with nearly-physical masses (N f = 2 + 1 + 1) [63]. Most of the lattice studies focused on the small Q 2 behavior of the pion form factor and the extraction of the pion charge radius. The pion charge radius is very sensitive to the quark mass. Chiral perturbation theory predicts a logarithmic divergence of the pion charge radius when the quark mass goes to zero [64]. Therefore, one has to work at the physical quark mass or have calculations performed in an appropriate range of quark masses to perform chiral extrapolations. Furthermore, studies have been performed for lattice spacing a > 0.09 fm. Constrained by the analyticity and unitarity, the charge radius is correlated with the phase of form factors in the timelike region. It is proposed in Ref. [65] that high-precision determinations of the pion form factor and the charge radius have potential to shed light on the discrepancy of hadronic vacuum polarization (HVP) derived from e + +e − → hadron cross-sections and lattice calculations [66].
The aim of this paper is to study the pion form factor in a wide range of Q 2 . Therefore, we perform calculations for small lattice spacings, namely a = 0.04fm and 0.06 fm, with valence pion mass of about 300 MeV. Furthermore, to study quark-mass effect, we also perform calculations at the physical pion mass, though at somewhat larger lattice spacing, a = 0.076 fm. Unlike previous studies, we also perform calculations for highly boosted pion in order to extend them in the future to the pion GPD.

II. LATTICE SETUP
In this study, we use Wilson-Clover action with hypercubic (HYP) [67] link smearing on (2+1)-flavor L t × L 3 s lattice ensembles generated by HotQCD collaboration [68,69] with highly-improved staggered quark (HISQ) sea action. For the clover coefficient we use the tree-level tadpole improved value c sw = u −3/4 0 , with u 0 being the HYP-smeared plaquette expectation value. This setup is the same as the one used by us to study the valence parton distribution of the pion [9,10]. As in Refs. [9,10], we use two lattice spacings a = 0.04 fm and a = 0.06 fm and the valence pion mass of 300 MeV. The lightest pion mass for these gauge configurations is m sea π = 160 MeV and the lattice spacings were fixed with the r 1 scale [68] using the value r 1 = 0.3106(18) fm [70]. In addition, we performed calculations at a lattices spacings of 0.076 fm and valence pion mass of 140 MeV using gauge configurations that correspond to the lightest pion mass of m sea π = 140 MeV [69]. The lattice spacing was set by the kaon decay constant, f K [69]. The lattice ensembles used in this study and the corresponding parameters are summarized in Table I. Due to the HISQ action, the taste splitting in the pion sector is small for lattice spacings a ≤ 0.076 fm. For a = 0.076 the root mean square pion mass is only 15% higher than the lightest pion mass, while the heaviest pion mass is only 25% above the lightest pion mass [69]. In what follows for a = 0.076 fm ensemble, will will not make a difference between the sea and the valence pion mass and refer to this ensemble as m π = 140 MeV ensemble or the ensemble with physical pion mass. The effects of partial quenching will persist at finite lattice spacings but will go away in the continuum limit.
To obtain the form factor we calculate the pion twopoint and three-point functions. We consider two-point functions defined as where π s (P, t) are either smeared or point sources, s = S, P , with spatial momentum P = 2π aL s · (n x , n y , n z ).
As in the previous studies [9,10], we used boosted Gaussian sources in Coulomb gauge with boost along the zdirection k z = 2π/(aL s ) · (0, 0, j z ). The radius of the Gaussian sources r G is also given in Table I. The threepoint function is defined as (4) being the isovector component of the electric charge operator. Note that the isosinglet component of the electric charge vanishes between the pion states. The initial momentum in the above expression is P i = 2π/(aL s ) · (0, 0, n z ), while the final momentum is P f = P = P i + q. The values of the momenta used in this study as well as the corresponding boost parameter j z are summarized in Table I. We calculated the three-point functions for three values of the source-sink separations, t s for the two coarser lattices. For the finest lattice we used four sourcesink separations. The source-sink separations used in our study are also listed in Table I.
The calculations of the two-and three-point functions were performed on GPUs with the QUDA multi-grid algorithm [71] used for the Wilson-Dirac operator inversions to get the quark propagators. We used multiple sources per configuration together with All Mode Averaging (AMA) technique [72] to increase the statistics. The stopping criterion for AMA was set to be 10 −10 and 10 −4 for the exact and sloppy inversions, respectively. Since the signal is deteriorating with increasing momenta, we use different number of sources and number of gauge configurations for different momenta. The number of gauge configurations and number of sources used in the analysis are given in the last two columns of Table I for each value of the momenta.

Ensemble:
m val π (GeV) csw rG fm ts/a nz ni (i = x, y) jz #cfgs (#ex,#sl) a = 0.076 fm, m sea π = 0.14 GeV, 0.14 1.0372 0. 59   For the study of the form-factor, it is convenient to use the Breit frame, where |P i | = |P f |. Using the Breit frame is essential when studying the GPD within LaMET [29][30][31][32], therefore we also calculated the pion form factor using the Breit frame. The parameters of this set-up are summarized in Table II.

III. TWO-POINT FUNCTION ANALYSIS
Since the source-sink separation values used in this study are not very large, it is important to quantify the contributions of the excited states when extracting pion matrix elements. This in turn requires a detailed study of the pion two-point functions. For a = 0.04 fm and 0.06 fm lattices and m val π = 300 MeV, the pion two-point functions have been studied for different momenta along the z-direction in Refs. [9,10]. Furthermore, this analysis was very recently extended to include momenta also along the x and y-directions for a = 0.04 fm [73]. We have extended this analysis to a = 0.076 fm and the physical pion mass.
The pion two-point function in Eq.
(2) has the following spectral decomposition: where E n+1 >E n , with E 0 being the energy of the pion ground state. A n is the overlap factor Ω|π s |n of the state n and the state created by operator π s from the vacuum state |Ω . Thanks to the Gaussian smearing, the excited state contribution is suppressed. So we truncate the Eq. (5) up to N state = 3 and then fit the data in a range of t ∈ [t min , aL t /2]. In the left panels of Fig. 1, we show the extracted E 0 for three different momenta. As one can see, the ground-state energies, E 0 reach a plateau when t min 10a, 5a and 2a for 1-state, 2-state and 3-state fits, respectively. The horizontal lines in the plots are computed from the dispersion relation Here the value of m π was obtained by considering the pion masses from the fits with t min ∈ [10a, 20a], and then fitting these results to a constant. The fit to a constant has there is no statistically significant t min dependence of the pion mass. The ground-state energies for different momenta agree with the horizontal lines for sufficiently large t min , i.e. follow the dispersion relation. Thus for the determination of the next energy level, we can fix the groundstate energy E 0 to be from the dispersion relation, and perform a 3-state fit. Interestingly, as shown in right panels of Fig. 1, we can also observe plateaus for E 1 when t min >5a. The energy of the first excited state also follows the dispersion relation E 1 (P) = P 2 + m 2 π with m π = 1.3 GeV. This could imply that the first excited state is single particle state, namely the first radial excitation of the pion π(1300) [73]. We cannot rule out, however, the possibility that it is a multi-pion states within the large errors. Since the first excited state energy, E 1 does not reach a plateau for t min <5a, we conclude that for t/a < 5 the contribution of higher excited states in the two-point function is significant. Therefore, we need to consider 3-state fits for these t values. To perform a 3-state fit, we fix E 0 to the dispersion relation and put a prior to E 1 using the best estimates from SS and smeared-point (SP) correlators [10] together with the errors from the 2state fit. This way we get the third excited state energy, E 2 , which does not depend on t min within the statistical errors. However, the value of E 2 is very large, about 3 GeV. This implies that E 2 does not actually refer to a single state but rather to a tower of many higher excited states. The situation is similar for other two 300 MeV ensembles [10]. Now we understand that a 2-state spectral model can describe our two-point functions well when t min 5a, while 3-state can describe t min 2a. This will be important to keep in mind when analyzing the three-point function and pion matrix elements in the next section.
(1, 32) 64 × 48 3 10 a = 0.04 fm, m sea π = 0.16, 0.3 9,12, 2 ±1 ∓2 120 (1, 32) 64 × 64 3 15,18 TABLE II. Two sets of measurements in the Breit frame on the two heavy-pion ensembles are shown. Using the notation similar to Table I, the initial pion state with transverse momentum P i ⊥ = 2πn p i /(Lsa), has the same energy as the final state with momentum P f = P i + q.   The lines are computed from the dispersion relation E(P) = P 2 + E(P = 0) 2 , with E(P = 0) to be 0.14 GeV for E0 and 1.3 GeV for E1. As can be observed, the E0 and E1 reach a plateau for large enough tmin.
To summarize this section, in Fig. 2 we show the dispersion relation obtained from the above analysis. We also extended the analysis for a = 0.06 fm [10] by including additional momenta with non-zero components along the x and y-directions. The corresponding results are also shown in Gaussian sources. These sources have poor overlap with the scattering states.

IV. EXTRACTION OF BARE MATRIX ELEMENTS OF PION GROUND STATE
To obtain the bare pion form factor we consider the following standard ratio of the three-point and two-point pion correlation functions [74,75] This ratio gives the bare pion form factor in As explained in Sec. II, we calculated the three-point functions with P i along theẑ direction, and multiple values of momentum transfer q = P f − P i for each P i . Thus there is no difference for q with same magnitude of the transverse momentum transfer. In other words, there should be transverse symmetry for the three-point function data. We find that indeed our numerical results for R f i (τ, t s ) with same |n q x | and |n q y | are consistent within the error. Therefore, we average the three-point functions data with same magnitude of the transverse momentum transfer in the following analysis.
Since the temporal extent of our lattices is not large, it is important to consider thermal state contaminations, also called wrap-around effects caused by the periodic boundary condition in time [10]. To remove the wraparound effects in the two-point function we replaced C 2pt (t) by C 2pt (t) − A 0 e −E0(aLt−t) using the best estimate of A 0 and E 0 from the two-point function analysis. To understand wrap-around effects in the three-point function we consider the spectral decomposition of C 3pt in Eq. (6) where m, n, k = Ω, 0, 1, . . . , with 0 being the pion ground state. In general, terms with non-zero E m will be highly suppressed by e −(aLt−ts)Em (we assume E Ω = 0). Therefore, in most studies such terms are neglected. However for the P = 0 case e −(aLt−ts)Em(P =0) = e −aLtmπ is not very small. We have e −aLtmπ ∼ 0.03, 0.003, 0.02 for a = 0.076, 0.06 and 0.04 fm lattices, respectively. On the other hand, for non-zero momenta the terms proportional to e −(aLt−ts)Em are smaller than 0.003 and can be neglected. Therefore, for a = 0.04 fm and 0.076 fm calculations we only consider non-zero momenta and limit the sum over index m in Eq. (7) to include only the vacuum state in what follows. We need, however, to consider the wrap-around effects when dealing with the renormalization, as discussed in the next section.
In this work, we use multi-state fit to extract the bare matrix elements of the ground state P f |O γt |P i ≡ 0P f |O γt |P i 0 by inserting the spectral decomposition of the two-point function in Eq. (5) and the three-point function in Eq. (7) with m = Ω, and the sum over n truncated to N state terms. Furthermore, we take the best estimate of A n and E n from the two-point function analysis. and put them into Eq. (6). In the following, we will refer to this method as Fit(N state , n sk ), in which N state is the number of states in the corresponding two-point function analysis and n sk labels how many τ points are skipped on the two sides of t s . We consider N state = 2 and N state = 3 that have four and nine fit parameters, respectively.
We perform multi-state fit using bootstrap method with time separations t s = 6a, 8a, 10a. The data with t s = 20a and n p i = (0, 0, 1) are used only to cross-check our analysis. Since the ratio defined in Eq. (6) is a derived quantity not defined on a single gauge configuration we used uncorrelated fits. The statistical correlation between the different data points is taken into account through the bootstrap procedure. In Fig. 3, we show the examples of ratio R f i (τ, t s ) as well as the 2-state and 3-state fit results. As one can see, for large momentum with large statistical errors, the reconstructed curves go through the data points well, and the 2-state and 3-state fit results are consistent with each other. However, this is not the case for smaller momentum, where the data are more precise. The 3-state fit is required to describe the ratio data with χ 2 /dof < 1, while the 2-state fit result in χ 2 /dof 1. Thus for the following analysis, we will take the 3-state fit results as the central value and use the corresponding statistical errors. However, even when using the 3-state fit there is no guarantee that we are free from excited state contamination. Therefore, we take the difference between the 2-state fit and the 3-state fit results as the systematic errors in the following analysis. It can be also observed that the data points of t s = 20a show plateau around t s /2 within the errors, and are also consistent with the 3-state fit results, which support our estimate of bare matrix elements. In App. B, we discuss   4. The forward matrix elements hB(P i , P i ). The P i z dependence can be described by hB(P i , P i ) = h ii B (P i = 0, P i = 0) + r(aP i z ) 2 shown as the line.
the plateau fit results using t s = 20a data.

V. THE PION FORM FACTORS
To obtain the form factor from the bare form factor determined in the previous section it needs to be multiplied by the vector current renormalization factor, Z V . The simplest way to obtain this is to calculate the forward matrix element h B (P i , P i ) = 0P i |O|P i 0 = Z −1 V . However, one needs to keep in mind the wrap-around effect discussed in the previous section. The other issue is cutoff dependence of h B (P i , P i ) at large values of P i . In Fig. 4, we show h B (P i , P i ) for a = 0.076 fm as a function of P i . In absence of discretization effects, h B (P i , P i ) should be independent of P i since after renormalization it gives the charge of the pion. In other words, Z V should not depend on the momentum of the external state. Following Ref. [10], we model the discretization effects using the form h B (P i , P i ) = h B (P i = 0, P i = 0) + r(aP i z ) 2 . As one can see from Fig. 4 this form describes the data quite well, except for P i = 0. The anomalously large value of h B (P i , P i ) at P i = 0 is due to the wrap-around effects as discussed in the previous section. This means that h B (P i , P i ) is contaminated by a small contribution proportional to e −aLtmπ mentioned in the previous section. This contribution is also proportional to matrix elements containing two or more pion states with the appropriate quantum numbers. Constraining such matrix elements is difficult in practice. However, under some physically well-motivated assumptions it is possible to estimate the corresponding contributions and remove them from h B (P i , P i ) [10]. Therefore, we follow the procedure explained in Appendix A of Ref. [10] to remove this contribution from the matrix element. The corrected result for h B (P i = 0, P i = 0) is shown as the blue point in Fig. 4 and is not very different from the result obtained by the fit. Thus we understand the discretization effects in the forward matrix element h B (P i , P i ). We also calculated Z V for a = 0.076 fm using RI-MOM scheme and MeV) ensemble (blue points), compared with the experiment data from CERN (red points) [36] and Fπ collaboration (green points) [46]. The purle bands are the dispersive analysis results of experimental data from Ref. [48], which also included form factors in time-like region. Our fit results of a = 0.076 fm data are shown as the blue bands, in which the filled band is from z-expansion fit and the dashed band is from monopole fit. The errors in this plot have included the systematic errors.
obtained Z V = 0.946 (12) which agrees with the results on h B (P i = 0, P i = 0) shown in Fig. 4 within errors. From Fig. 4 we also see that the discretization errors are smaller than 1% for P i z < 1 GeV , and are less than 2% for P i z < 1.6 GeV. Since the discretization effects as functions of P i z will be similar for off-forward matrix element it is convenient to obtain the renormalized pion form factor by simply dividing h B (P f , P i ) by h B (P i , P i ). Then we have F π (Q 2 = 0) = 1 by construction and the discretization errors for large P i z are removed. We still may have discretization errors proportional to (aQ) 2 . Assuming that these discretization errors are similar to the (aP i z ) 2 discretization errors we can neglect them. This is because other sources of errors for the form factors are significantly larger for the considered Q 2 range as we will see below. We comment further on the cutoff dependence in the form-factor in App. A.
In Fig. 5, we show the renormalized pion form factors obtained for the m π = 140 MeV ensemble and compared to the experimental data from CERN [36], as well as the results from F π collaboration [46]. The purple bands are the dispersive analysis results of experimental data from Ref. [48], which also included form factors in timelike region. We see good agreement between the lattice results and the experimental data within the estimated error bars at low Q 2 . It is expected that at low Q 2 , the pion form factors can be described well by a simple monopole Ansatz motivated by the Vector Meson Dominance (VMD) model [76] F π (Q 2 ) = The monopole mass M should be close to the ρ meson mass. Therefore, in Fig. 5 we show the inverse of the pion from factor, 1/F π (Q 2 ), as a function of Q 2 . We see that in the studied range of Q 2 the inverse form factor can be roughly described by a linear function up to Q 2 = 0.4 GeV within the errors, as expected from monople form. The monopole fit of the lattice data (dashed band in Fig. 5) extended to higher Q 2 also agrees with the pion form factor obtained by F π collaboration [46], possibly indicating that the monopole form may work in an extended range of Q 2 within the current precision. At very low Q 2 , the pion form factor can be characterized in terms of the pion charge radius As mentioned in the introduction, the pion charge radius is very sensitive to the quark mass, and it is clearly seen in the lattice calculations. In fact, it appears to be challenging to obtain the correct pion charge radius from the lattice results [51][52][53][54][55][56][57][58][59][60][61][62][63]. The lattice calculations at the unphysical quark masses lead to smaller pion charge radius than the experimental results. If the monopole form (8) could describe the pion form factor for all Q 2 the pion charge radius would be related to the monopole mass as It is convenient to represent the form factors in terms of the effective charge radius defined as [51] In Fig. 6 we show the effective radius for a = 0.076 fm ensemble as well as for the two finer ensembles with m val π = 300 MeV. We see from the figure that r 2 ef f is roughly constant as a function of Q 2 for all three lattice spacings. For the smallest lattice spacing, a = 0.04 fm the results on the effective radius are Q 2 -independent for Q 2 as high as 1.4 GeV 2 . This is consistent with earlier findings [51]. We also clearly see the quark mass dependence of r 2 ef f . The effective radius is smaller for the heavier pion mass as expected. Comparing the results at a = 0.06 fm and a = 0.04 fm we see no clear lattice spacing dependence of r 2 ef f . Therefore, we conclude that for a = 0.06 fm the discretization errors for the pion form factor are smaller than the estimated lattice errors in the range of Q 2 studied by us. Finally, for the two finer lattices we also show the results from the calculations using Breit frame, which agree with the non-Breit frame results.
While the monopole Ansatz seems to describe the pion form factor well and was used to obtain the pion charge radius in the past (see e.g. Ref. [51]) there is no strong theoretical reason why it should describe the pion form factor. Therefore, one has to consider an alternative and more flexible parameterization of the pion form factor. An alternative way to fit the form factors is the model independent method called the z-expansion [77]. Here the form factor is written as where t = −Q 2 , a k are the fit parameters with constrain condition F π (Q 2 = 0) = 1, and t cut = 4m 2 π is the two-pion production threshold. Furthermore, t 0 is chosen to be the 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 used for the fit. In the timelike region near the two pion threshold, the leading singularity of form factor should be proportional to (4m 2 π −t) 3/2 due to the P-wave nature of the π−π scattering [48,78,79], which leads to the additional constraint kmax k=1 (−1) k ka k = 0. We use AIC model selection rules to determine k max , which are 2 for a = 0.06 fm, and 3 for a = 0.04, 0.076 fm data and for the Q 2 under consideration. The z expansion results are also shown in Fig. 5 and appear to overlap with the monopole fit, but for larger Q 2 it has larger errors. We also show the fits with the z-expansion in Fig. 6. From this figure we see that this fit works well also for the valence pion mass of 300 MeV and naturally reproduces little Q 2 dependence of the effective radii. To better understand the quark mass dependence of the pion form factor as well to facilitate the comparison with the experimental results, in Fig. 7 we show all the results for the pion form factor in terms of the effective radius r ef f (Q 2 ). We see that the effective radius obtained for the physical pion mass is clearly larger than the one obtained for m val π = 300 MeV and is much closer to the CERN data. Furthermore, the fits of r ef f for m val π = 300 MeV for the two lattice spacings agree within errors. While the individual lattice data and the CERN data appear to agree within errors we also see from the figure that there is a tendency for the CERN data to lie higher than the lattice data. This leads to a slight difference in the pion charge radius as discussed below.
The pion charge radius can be derived from zexpansion fit results using Eq. (9), which are summarized in Table III for the three lattice spacings used in this work. We also discuss the radius obtained from the monopole fit for comparison in App. C. As expected the calculations for the heavier quark mass give smaller pion charge radius. Since the z-expansion provides a model independent way to obtain the pion charge radius, for our final estimate of the pion charge radius at the physical point we take the result from the z-expansion fit: where we added the statistical and systematic errors (defined by the difference between the results from 2state and 3-state fit of matrix elements) in quadrature. This result is consistent the pion charge radius quoted by Particle Data Group (PDG), r 2 π PDG = 0.434(5) fm 2 [80], which is averaged from determination from t-channel πe→πe scattering data [34,36,81] and schannel e + e − →π + π − data sets [48,82]. The HPQCD determination that uses HISQ action both in the sea and the valence sectors of (2 + 1 + 1)-flavor QCD is r 2 π = 0.403(18)(6) fm 2 [63]. The most precise lattice determination of the pion charge radius in 2+1 flavor QCD using overlap action in the valence sector and domain wall action in the sea sector has r 2 π = 0.436(5)(12) fm 2 [62]. The 2+1 flavor domain wall calculation gives r 2 π = 0.434(20)(13) fm 2 [61]. Finally, the other 2+1 flavor lattice determinations of the pion charge radius have significantly larger errors [59,60]. We summarize the comparison in Fig. 8.

VI. CONCLUSIONS
In this paper we studied the pion form factor in 2+1 flavor lattice QCD using three lattices spacings a = 0.076, a = 0.06 and a = 0.04 fm. The calculations on the coarsest lattice have been performed with the physical value of the quark masses, while for the finer two lattices the valence pion mass was 300 MeV. We have found that the pion form factor is very sensitive to the quark mass, as expected. We showed that lattice discretization effects are quite small for lattice spacings smaller than 0.06 fm. For the physical quark masses our lattice results on the pion form factor appear to agree with the experimental determinations. Unlike other lattice studies we also considered highly boosted pions in the initial state using momentum boosted Gaussian sources. In addition we performed calculations also in the Breit frame. We demonstrated that the calculations of the pion form factor performed at different momenta of the pion as well as in the Breit frame give consistent results. This is very important for extending the calculations to pion GPDs.
An important outcome of our analysis is that the monopole Ansatz can describe the pion form factor in large range of Q 2 , up to Q 2 = 1.4 GeV 2 . In the future it will be important to extend the calculations to even higher momentum transfer given the experimental efforts in Jlab and EIC. To do this we should use boosted sources that also depend on the value of Q 2 . At present the momentum boost was optimized only according to the pion momentum in the initial state.
From the low Q 2 dependence of the pion form factor we determined the pion charge radius, which is one sigma lower that the experimental result. We speculated, whether this is due to the effect of partial quenching. To fully resolve this issue calculations at smaller lattice spacing with the physical value of the pion masses are needed.  As is shown in Fig. 4, there are 2% discretization effects of Z −1 V (P i ) = h B (P i , P i ). We chose to divide h B (P f , P i ) by h B (P i , P i ) so that the renormalized pion form factors could reduce such effects. To estimate the impact of the discretization errors to the form factors as well as pion charge radius, instead we can renormalize the bare form factors h B (P f , P i ) by a constant Z −1 V such as Z −1 V (0.25 GeV) of a = 0.076 fm ensemble. The effective radius for a = 0.076 fm ensemble is shown in Fig. 9, and in this case we estimate the charge radius from monopole fit and z-expansion fit as 0.406(6)(25) fm 2 and 0.427(10)(22) fm 2 , which shift 2% but are consistent with the estimates in Table III. It has been observed in Sec. IV that the ratio R f i (τ, t s ) of t s = 20a shows plateau around t s /2 which is also consistent with the results from Fit(3,2) method, implying that the smallness of excited-state contribution in this region. Therefore it is reasonable to perform a one-state fit, namely plateau fit, to extract the bare matrix elements. We denote this method by Plateau(τ min , τ max ) which fit R f i (τ, t s = 20a) of τ ∈ [τ min , τ max ] to a constant.
The fit results from Plateau(τ min , τ max ) are shown in Fig. 10 as the blue bands where the multi-state fit results are also shown for comparison. Clearly, the plateau fit shows good agreement with 3-state fit results. In Fig. 11, we show the distribution of difference between plateau fit and multi-state fit using bootstrap samples. In the main text, we have taken the difference between 2-state and 3-state fit as the systematic errors of excited-state contamination. It can be seen that such an estimate is larger than the difference between plateau fit and 3-state fit which should give a sufficiently conservative total error.
We also determined the pion form factor from the plateau fits for t s = 20 The corresponding results in terms of the effective radius are shown in Fig. 12. Once again, consistent results between Plateau(τ min , τ max ) and Fit (3,2) can be observed.

Appendix C: Model dependence of radius extraction
In this work, we used z-expansion Ansatz to obtain the charge radius from the pion form factors shown in Table III. For comparison, in Table IV we also show the radius obtained from monopole fit whose statistical error are often smaller, but this fit has larger systematic errors compared to z-expansion. Both fits produce good χ 2 /df . For the a = 0.076 fm ensemble, for example we get, χ 2 /df = 0.56 for monopole fit, and χ 2 /df = 0.51 for z-expansion fit. Within the estimated errors the two fit forms give consistent results but only marginal. In Fig. 13, we show the effective radius (c.f. Eq. (11)) calculated from the z-expansion fit (blue band) as well as monopole fit (red band). Clearly the z-expansion fit is more flexible so that the effective radius is a function of Q 2 rather than a constant. At Q 2 = 0 where the charge radius is defined, the result from z-expansion fit ( r 2 Z ) is higher than monopole fit ( r 2 M ). We show the distribution of r 2 M nst3 − r 2 Z nst3 from bootstrap samples in Fig. 14, where the N-state fit is denoted by nstN. The central value of this distribution is 0.02 fm 2 .