Pion form factor and charge radius from Lattice QCD at physical point
Xiang Gao, Nikhil Karthik, Swagato Mukherjee, Peter Petreczky, Sergey Syritsyn, Yong Zhao
PPion form factor and charge radius from Lattice QCD at physical point
Xiang Gao,
1, 2, ∗ Nikhil Karthik,
3, 4
Swagato Mukherjee, Peter Petreczky, Sergey Syritsyn,
5, 6 and Yong Zhao
7, 1 Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA Physics Department, Tsinghua University, Beijing 100084, China Department of Physics, College of William & Mary, Williamsburg, VA 23185, USA Thomas Jefferson National Accelerator Facility, Newport News, VA 23606, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA RIKEN-BNL Research Center, Brookhaven National Lab, Upton, NY, 11973, USA Physics Division, Argonne National Laboratory, Lemont, IL 60439, USA (Dated: February 12, 2021)We present our results on the electromagnetic form factor of pion over a wide range of Q usinglattice QCD simulations with Wilson-clover valence quarks and HISQ sea quarks. We study theform factor at the physical point with a lattice spacing a = 0 .
076 fm. To study the lattice spacingand 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 withthe experimental results. Through fits to the form factor, we estimate the charge radius of pion forphysical pion mass to be (cid:104) r π (cid:105) = 0 . . I. INTRODUCTION
Pion is one of the most prominent strongly interact-ing particle next to the nucleon since it is a Goldstoneboson of QCD. For this reason it is important to studythe pion internal structure and find out if there is a con-nection between its internal structure and its Goldstoneboson nature. This issue is particularly relevant for un-derstanding the origin of mass generation in QCD, seee.g. discussions in Refs. [1, 2].The knowledge of internal structure of the pion is muchmore limited than that of the nucleon. On the par-tonic level the parton distribution function (PDF) of thepion can been studied through the global analysis of theDrell-Yan production in pion-nucleon collisions and intagged deep inelastic scattering (DIS), for recent analy-ses see Refs. [3, 4]. Recently there have been many ef-forts in lattice QCD to study the pion PDF [5–10] whichuse the quasi-PDF in Large Momentum Effective The-ory [11, 12], the pseudo-PDF [13, 14] and current-currentcorrelator [15–17] (also referred to as good lattice crosssection) approaches, see Refs. [18–21] for recent reviews.Lattice calculations of the lowest moments of pion PDF[22–27] are also available and can be used as additionalconstraints in the global analysis.Form factor, defined as (cid:104) P | J µ | P (cid:105) = ( P + P ) µ F π ( Q ) , (1)with J µ being the electro-magnetic current and Q = − ( P − P ) , provide a different insight into pion struc-ture, namely the charge distribution. It can be in prin-ciple measured in electron pion scattering. Generalizedparton distribution (GPD) combine the information con-tained in PDF and form factors and provide a three-dimensional image of a hadron. In the case of the nucleon ∗ [email protected] the study of GPD is the subject of large experimentaland theory efforts (see e.g. Ref. [28] for a recent review).The experimental study of the pion GPD is far more chal-lenging and will be only possible at Electron-Ion Collider(EIC), if at all. Fortunately GPDs can be studied on thelattice using LaMET, including pion GPDs [29–32].Experimentally, the pion form factor was measuredby scattering of pions off atomic electrons in Fermilab[33, 34] and CERN [35, 36]. This allowed to determinethe pion form factor for momentum transfer Q up to0 .
253 GeV [33–36]. For larger Q one has to deter-mine the pion form factor from the electro-productionof charged pions off nucleons. The corresponding exper-iments have been performed in Cornell [37–39] DESY[40, 41], and Jlab [42–46]. This determination, however,is model dependent. Recent determination of the pionform factor up to Q of 2 .
45 GeV is carried out by the F π collaboration using data both from DESY and JLab[46]. Experiments at the future EIC facility will allow toprobe even higher Q up to 30 GeV and possibly see thepartonic structure in a exclusive elastic process and makecontact with asymptotic large Q perturbative behavior[47].Lattice QCD calculations allow to obtain the pion formfactor from first principles, i.e. without any model de-pendence, up to relatively large Q . Therefore, they willprovide an important cross-check for the experimentaldeterminations. The first lattice calculations of the pionform factor date back to late 80s and have been per-formed in the quenched approximation [48, 49]. Morerecently lattice calculations of the pion form factor havebeen performed with two flavors ( N f = 2) of dynamicalquarks [50–54], with physical strange quark and two lightquark flavors ( N f = 2+1) [55–61], as well as with dynam-ical charm quark, strange quark and two flavors of thelight quarks with nearly physical masses ( N f = 2 + 1 + 1)[62]. Most of the lattice studies focused on the small Q behavior of the pion form factor and the extraction ofthe pion charge radius. The pion charge radius is very a r X i v : . [ h e p - l a t ] F e b sensitive to the quark mass. Chiral perturbation the-ory predicts a logarithmic divergence of the pion chargeradius when the quark mass goes to zero [63]. There-fore, one has to work at the physical quark mass or havecalculations performed in an appropriate range of quarkmasses to perform chiral extrapolations. Furthermore,studies have been performed for lattice spacing a > . Q . Therefore, we perform calculationsfor small lattice spacings, namely a = 0 . .
06 fmwith valence pion mass of about 300 MeV. Furthermore,to study quark mass effect we also perform calculationsat the physical pion mass, though at somewhat largerlattice spacing, a = 0 .
076 fm. Unlike previous studieswe also perform calculations for highly boosted pion inorder to extend them in the future for the pion GPD.
II. LATTICE SETUP
In this study we use Wilson-Clover action with hy-percubic (HYP) [64] smearing on 2+1 flavor staggered L t × L s lattices generated with highly improved staggeredquark (HISQ) action by HotQCD collaboration [65, 66].For the clover coefficient we use the tree-level tadpoleimproved value c sw = u − / , with u being the HYPsmeared plaquette expectation value. This setup is thesame as the one used by us to study the valence partondistribution of the pion [9, 10]. As in Refs. [9, 10] weuse two lattice spacings a = 0 .
04 fm and a = 0 .
06 fmand the valence pion mass of 300 MeV. The lightest pionmass for these gauge configurations is m seaπ = 160 MeVand the lattice spacings was fixed through r scale [65]using the value r = 0 . .
076 fmand valence pion mass of 140 MeV using the gauge con-figurations that correspond to the lightest pion mass of m seaπ = 140 MeV [66]. The lattice spacings was set bythe kaon decay constant, f K [66]. The lattice ensemblesused in this study and the corresponding parameters aresummarized in Table I. Because of the use of the HISQaction the taste splitting in the pion sector is small forlattice spacings a ≤ .
076 fm. For a = 0 .
076 the rootmean square pion mass is only 15% higher than the light-est pion mass, while the heaviest pion mass is only 25%above the lightest pion mass [66]. In what follows for a = 0 .
076 fm ensemble will will not make a difference be-tween the sea and the valence pion mass and refer to thisensemble as m π = 140 MeV ensemble or the ensemblewith physical pion mass. The effects of partial quench-ing will persist at finite lattice spacings but will go awayin the continuum limit.To obtain the form factor we calculate the pion two-point and three point functions. We consider two pointfunctions defined as C ss (cid:48) ( t ; P z ) = (cid:68) π s ( P , t ) π † s (cid:48) ( P , (cid:69) , (2) 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 Gaus-sian sources in Coulomb gauge with boost along the z -direction k z = 2 π/ ( aL s ) · (0 , , j z ). The radius of theGaussian sources r G is also given in Table I. The three-point function is defined as C ( P f , P i , τ, t s ) = (cid:68) π S ( P f , t s ) O γ t ( τ ) π † S ( P i , (cid:69) , (3)with O γ t ( τ ) = (cid:88) x e − i ( P f − P i ) x (cid:20) u ( x ) γ t u ( x ) − d ( x ) γ t d ( x ) (cid:21) , x = ( x , τ )(4)being the iso-vector component of the electric charge op-erator. Note that the iso-singlet component of the elec-tric charge vanishes between the pion states. The initialmomentum in the above expression is P i = 2 π/ ( aL s ) · (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 asthe corresponding boost parameter j z are summarizedin Table I. We calculated the three point functions forthree values of the source-sink separations, t s for the twocoarser lattices. For the finest lattice we used four source-sink separations. The source-sink separations used in ourstudy are also listed in Table I.The calculations of the two and three point functionswere performed on GPU, with the QUDA multi-grid al-gorithm used for the Wilson-Dirac operator inversions toget the quark propagators. We used many sources to-gether with All Mode Averaging (AMA) technique [68]to increase the statistics. The stopping criterion for AMAwas set to be 10 − and 10 − for the exact and sloppyinversions respectively. Since the signal is deteriorat-ing with increasing momenta we use different number ofsources and number of gauge configurations for differentmomenta. The number of gauge configurations and num-ber of sources used in the analysis are given in the lasttwo columns of Table I for each value of the momenta.For the study of the form-factor it is useful to use theBreit frame, where | P i | = | P f | . The use of the Breitframe is essential when studying the GPD within LaMET[29–32]. Therefore we also calculated the pion form factorusing the Breit frame. The parameters of this set-up aresummarized in Table II. III. TWO-POINT FUNCTION ANALYSIS
Since the source-sink separation used in this study arenot very large it is important to quantify the contribu-tions of the excited states when extracting pion matrixelements. This in turn requires a detailed study of thepion two point functions. For a = 0 .
04 fm and 0 .
06 fm
Ensemble: m valπ (GeV) c sw t s /a r G fm n z n i ( i = x, y ) j z a = 0 .
076 fm, m seaπ = 0 .
14 GeV, 0.14 1.0372 6, 8, 10 0.342 [0,3] ± ± , × [4,7] ± ± , a = 0 .
06 fm, m seaπ = 0 .
16 GeV, 0.3 1.0336 8, 10, 12 0.312 [0,1] ± ± , × [2,3] ± ± , ± ± , ± , a = 0 .
04 fm, m seaπ = 0 .
16 GeV 0.3 1.02868 9,12, 0.208 [0,1] ± , × ± , ± , m valπ t s /a n pz n pi n qi a, L t × L s (GeV) i = x, y i = x, ya = 0 .
06 fm, m seaπ = 0 .
16, 0.3 8, 2 ± ∓ , × a = 0 .
04 fm, m seaπ = 0 .
16, 0.3 9,12, 2 ± ∓ , × P i ⊥ = 2 πn pi / ( L s a ), and final state P f = P i + q have the same energy. lattices and m valπ = 300 MeV the pion two point func-tions have been studied for different momenta along the z -direction in Refs. [9, 10]. Furthermore, this analy-sis was very recently extended to include momenta alsoalong the x and y -directions for a = 0 .
04 fm [69]. Weextended this analysis to a = 0 .
076 fm and the physicalpion mass.The pion two point function in Eq. (2) has the spectraldecomposition: C ss (cid:48) ( t ) = N state − (cid:88) n =0 A sn A s (cid:48) ∗ n ( e − E n t + e − E n ( aL t − t ) ) , (5)where E n +1 >E n , with E being the energy of the pionground state. A n is the overlap factor (cid:104) Ω | π s | n (cid:105) of thestate n and the state created by operator π s from thevacuum state | Ω (cid:105) . Thanks to the Gaussian smearing, theexcited state contribution is suppressed. So we truncatethe Eq. (5) up to N state = 3 and then fit the data in arange of t between [ t min , a L t /2]. The one-state fit re-sults of the smeared-smeared (SS) two point function areshown in left panels of Fig. 1 for three different momenta.As one can see, the ground state energies, E reach aplateau when t min (cid:38) a . And at P = 0, the plateau isaround m π = 140 MeV, which is the physical pion mass.The lines in the plots are computed from the dispersionrelation E ( P ) = (cid:112) P + m π , and show good consistencywith the plateaus. Thus for the determination of thenext energy level, we can fix the ground state energy E to be from the dispersion relation, and perform a 2-state fit. Interestingly, as shown in right panels of Fig. 1,we can also observe plateaus for E when t min > E ( P ) = (cid:112) P + m π (cid:48) wellwithin the errors. This could imply a single-particle ex-cited state, namely the first radial excitation of the pion π (1300) [69]. Since the first excited state energy, E doesnot reach a plateau for t min < t/a < t values. To performthe three state fit, we fix E to the dispersion relationand put a prior to E using the best estimates from SSand smeared-point (SP) correlators [10] together with theerrors from the two-state fit. This way we get the thirdexcited state energy, E , which does not depend on t min within the statistical errors. However, the value of E isvery large, about 3 GeV. This implies that E does notactually belong to a single state but rather to a tower ofmany higher excited states. The situation is similar forother two 300 MeV ensembles [10].Now we understand that a two-state spectral modelcan describe our two-point functions well when t min (cid:38) a ,while three-state can describe t min (cid:38) a . This will be im-portant to keep in mind when analyzing the three pointfunction and pion matrix elements in the next section.To summarize this section, in Fig. 2 we show the dis-persion relation obtained from the above analysis. Wealso extended the analysis for a = 0 .
06 fm [10] by in- t min / a E ( G e V ) SS, P = (0, 0, 0) E , N state = 1 t min / a E ( G e V ) SS, P = (0, 0, 0) E , N state = 2 E , N state = 3 E , N state = 32.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 t min / a E ( G e V ) SS, P = (0, 0, 0.76) GeV E , N state = 1 2 3 4 5 6 7 8 t min / a E ( G e V ) SS, P = (0, 0, 0.51) GeV E , N state = 2 E , N state = 3 E , N state = 32 4 6 8 10 12 14 16 t min / a E ( G e V ) SS, P = (0.51, 0, 0.76) GeV E , N state = 1 2 3 4 5 6 7 8 t min / a E ( G e V ) SS, P = (0.51, 0, 0.76) GeV E , N state = 2 E , N state = 3 E , N state = 3 FIG. 1. E from 1-state fit (left) and E , E from constrained 2-state and 3-state fit (right) for three different momentum areshowed as a function of t min . The lines are computed from the dispersion relation E ( P ) = (cid:112) P + E ( P = 0) , with E ( P = 0)to be 0.14 GeV for E and 1.3 GeV for E . As can be observed, the E and E reach a plateau for large enough t min . cluding additional momenta with non-zero componentsalong the x and y -directions. The corresponding resultsare also shown in Fig. 2. We clearly see the effect ofthe quark masses. For the physical pion mass the firstexcited state has the correct mass of π (1300), while forthe larger quark mass ( a = 0 .
06 fm) the excited state isabout 200 MeV higher.
IV. EXTRACTION OF BARE MATRIXELEMENT OF PION GROUND STATE
To obtain the bare pion form factor we consider thefollowing standard ratio of the three point and two point pion functions [70, 71] R fi ( τ, t s ) ≡ (cid:113) P f P i P f + P i C ( P f , P i , τ, t s ) C ( t s , P i ) × (cid:18) C ( t s − τ, P f ) C ( τ, P i ) C ( t s , P i ) C ( t s − τ, P i ) C ( τ, P f ) C ( t s , P f ) (cid:19) / . (6)This ratio gives the bare pion form factor in the limit t s → ∞ : h B ( P f , P i ) = lim t s →∞ R fi ( τ, t s ).As explained in Sec. II, we calculated the three-pointfunctions with P i along with z direction, and multiplemomentum transfer q = P f − P i for each P i . Thus thereis no difference for q with same magnitude of the trans-verse momentum transfer. In other words, there shouldbe transverse symmetry for the three-point function data.We find that indeed our numerical results for R fi ( τ, t s )with same | n qx | and | n qy | are consistent within the error. | P |[GeV] E ( P ) [ G e V ] a = 0.076 fm E E E | P |[GeV] E ( P ) [ G e V ] a = 0.06 fm E E E FIG. 2. Dispersion relation determined by the plateau of Fig. 1 for the physical pion mass ensemble (left) and a = 0.06 fmensemble (right). The lines are dispersion relation calculated by E ( P ) = (cid:112) P + E ( P = 0) . Therefore, we average the three-point functions data withsame magnitude of the transverse momentum transfer inthe following analysis.Since the temporal extent of our lattices is not largeit is important to consider thermal contamination, alsocalled wrap-around effects, caused by the periodic bound-ary condition [10]. To remove the wrap-around ef-fects in the two point function we replaced C ( t ) by C ( t ) − A e − E ( aL t − t ) using the best estimate of A and E from the 2-point function analysis. To under-stand wrap-around effects in the three point function weconsider the spectral decomposition of C in Eq. (6) (cid:104) π S ( P f , t s ) O γ t ( τ ) π † S ( P i , (cid:105) = (cid:88) m,n,k (cid:104) m | π S | n (cid:105)(cid:104) n | O γ t | k (cid:105)(cid:104) k | π † S | m (cid:105)× e − τE k e − ( t s − τ ) E n e − ( aL t − t s ) E m , (7)where m, n, k = Ω , , , . . . , with 0 being the pionground state. In general, terms with non-zero E m will be highly suppressed by e − ( aL t − t s ) E m (we assume E Ω = 0). Therefore, in most studies such terms are ne-glected. However for the P = 0 case e − ( aL t − t s ) E m ( P =0) = e − aL t m π is not very small. We have e − aL t m π ∼ e − ( aL t − t s ) E m are smaller than0.003 and can be neglected. Therefore, for a = 0 .
04 fmand 0 .
076 fm calculations we only consider non-zero mo-menta and limit the sum over index m in Eq. (7) toinclude only the vacuum state in what follows. We need,however, to consider the wrap-around effects when deal-ing with the renormalization, as discussed in the nextsection. In this work, we use multi-state fit to extract thebare matrix elements of the ground state (cid:104) P f | O γ t | P i (cid:105) ≡(cid:104) P f | O γ t | P i (cid:105) by inserting the spectral decompositionof the two point function in Eq. (5) and the three pointfunction in Eq. (7) with m = Ω, and the sum over n trun-cated to N state terms. Furthermore, we take the best esti-mate of A n and E n from the two-point function analysis.and put them into Eq. (6). In the following we will referto this method as Fit( N state , n sk ), in which N state is thenumber of states in the corresponding 2-point functionanalysis and n sk labels how many τ we skipped on thetwo sides of t s . We consider N state = 2 and N state = 3that imply four and nine fit parameters, respectively.In Fig. 3, we show the examples of ratio R fi ( τ, t s ) aswell as the 2-state and 3-state fit results. As one cansee, for large momentum, the 2-state and 3-state fit re-sults are consistent with each other because of the largestatistical errors, while this is not the case for smaller mo-mentum, where the data are more precise. The 3-statefit can describe the ratio data better than the 2-state fit.Thus for the following analysis, we will take the 3-statefit results as the central value and use the correspondingstatistical errors. However, even when using the 3-statefit there is no guarantee that we are free from excitedstate contamination. Therefore, we take the differencebetween the 2-state fit and the 3-state fit results as thesystematic errors in the following analysis. V. THE PION FORM FACTORS
To obtain the form factor from the bare form factordetermined in the previous section it needs to be multi-plied by the vector current renormalization factor, Z V . ( t s /2) a R f i ( , t s ) n p = (0, 0, 1) n q = (0, 0, 0) Fit(2, 3)Fit(3, 2)Fit(2, 3)Fit(3, 2) t s =6a t s =8a t s =10a ( t s /2) a R f i ( , t s ) n p = (0, 0, 3) n q = (0, 0, 0) Fit(2, 3)Fit(3, 2)Fit(2, 3)Fit(3, 2) ts6ts8ts10 ( t s /2) a R f i ( , t s ) n p = (0, 0, 1) n q = (2, 0, 0) Fit(2, 3)Fit(3, 2)Fit(2, 3)Fit(3, 2) t s =6a t s =8a t s =10a ( t s /2) a R f i ( , t s ) n p = (0, 0, 3) n q = (2, 0, 0) Fit(2, 3)Fit(3, 2)Fit(2, 3)Fit(3, 2) t s =6a t s =8a t s =10a ( t s /2) a R f i ( , t s ) n p = (0, 0, 1) n q = (2, 2, 0) Fit(2, 3)Fit(3, 2)Fit(2, 3)Fit(3, 2) t s =6a t s =8a t s =10a ( t s /2) a R f i ( , t s ) n p = (0, 0, 3) n q = (2, 2, 0) Fit(2, 3)Fit(3, 2)Fit(2, 3)Fit(3, 2) t s =6a t s =8a t s =10a FIG. 3. R fi ( τ, t s ) for n p i = (0 , ,
1) (left) and (0 , ,
3) (right) for n q = (0 , , , (2 , , , (2 , ,
0) of physical ensemble areshown. The curves are the central value of multi-state fit Fit(2,3) (dashed) and Fit(3,2) (solid), and the bands are the estimatedbare matrix elements.
The simplest way to obtain this is to calculate the for-ward matrix element h B ( P i , P i ) = (cid:104) P i | O | P i (cid:105) = Z − V .But one needs to keep in mind the wrap-around effectdiscussed in the previous section. The other issue is cut-off dependence of h B ( P i , P i ) at large values of P i . InFig. 4, we show h B ( P i , P i ) for a = 0 .
076 fm as func-tion of P i . In absence of discretization effects h B ( P i , P i )should be independent of P i since after renormalizationit gives the charge of the pion. In other words Z V shouldnot depend on the momentum of the external state. Fol-lowing Ref. [10] we model the discretization effects usingthe form h B ( P i , P i ) = h B ( P i = 0 , P i = 0) + r ( aP iz ) .As one can see from Fig. 4 this form describes the dataquite well, except for P i = 0. The anomalously largevalue of h B ( P i , P i ) at P i = 0 is due to the wrap-aroundeffects as discussed in the previous section. This meansthat h B ( P i , P i ) is contaminated by a small contributionproportional to e − aL t m π mentioned in the previous sec- tion. This contribution is also proportional to matrixelements containing two or more pion states with the ap-propriate qunatum numbers. Constraining such matrixelements is difficult in practice. However, under somephysically well motivated assumption it is possible to es-timate the corresponding contributions and remove themfrom h B ( P i , P i ) [10]. Therefore, we follow the procedureexplained in Appendix. A of Ref. [10] to remove this con-tribution from the matrix element. The corrected resultfor h B ( P i = 0 , P i = 0) is shown as the blue point inFig. 4 and is not very different from the result obtainedby the fit. Thus we understand the discretization effectsin the forward matrix element h B ( P i , P i ). We also cal-culated Z V for a = 0 .
076 fm using RI-MOM scheme andobtained Z V = 0 . h B ( P i = 0 , P i = 0) shown in Fig. 4 within errors.From Fig. 4 we also see that the discretization errorsare smaller than 1% for P iz < | P i |(GeV) h B ( P i , P i ) h B ( P i , P i )periodicity corrected FIG. 4. The forward matrix elements h B ( P i , P i ). The P iz de-pendence can be described by h B ( P i , P i ) = h iiB ( P i = 0 , P i =0) + r ( aP iz ) shown as the line.
2% for P iz < . P iz will be similar for off-forward matrixelement it is convenient to obtain the renormalized pionform factor by simply dividing h B ( P f , P i ) by h B ( P i , P i ).Then we have F π ( Q = 0) = 1 by construction and thediscretization errors for large P iz are removed. We stillmay have discretization errors proportional to ( aQ ) . As-suming that these discretization errors are similar to the( aP iz ) discretization errors we can neglect them. This isbecause other sources of errors for the form factors aresignificantly larger for the considered Q range as we willsee below. We comment further on the cutoff dependencein the form-factor in Appendix A.In Fig. 5, we show the renormalized pion form factorsobtained for the m π = 140 MeV ensemble and comparedto the experimental data from CERN [36], as well as theresults from F π collaboration [46]. We see good agree-ment between the lattice results and the experimentaldata within the estimated error bars at low Q . It is ex-pected that for low Q the pion form factors can be welldescribed by a simple monopole Ansatz motivated by theVector Meson Dominance (VMD) model [72] F π ( Q ) = 11 + Q /M . (8)The monopole mass should be close to the ρ meson mass.Therefore, in Fig. 5 we show the inverse of the pion fromfactor, 1 /F π ( Q ) as function of Q . We see clearly that inthe studied range of Q the inverse form factor can be welldescribed by a linear function, as expected from mono-ple form, at least up to Q = 0 . Q also agrees with the pion form factor ob-tained by F π collaboration [46], possibly indicating thatthe monopole form may work in an extended range of Q .At very low Q the the pion form factor can be char- Q (GeV ) F ( Q ) CERN F collaborationa = 0.076 fm Q (GeV ) / F ( Q ) CERNa = 0.076 fm
FIG. 5. Pion form factors (upper panel) and the inverse formfactors (lower panel) derived from the a = 0.076 fm ( m π = 140MeV) ensemble (blue points), compared with the experimentdata from CERN (red points) [36] and F pi collaboration (yel-low points) [46]. The bands are the fit result of our a = 0.076fm data, in which the filled band is from z -expansion fit andthe dashed band is from monopole fit. acterized in terms of the pion charge radius r π = − dF π ( Q ) dQ | Q =0 . (9)As mentioned in the introduction the pion charge radiusis very sensitive to the quark mass, and is clearly seenin the lattice calculations. In fact, it appears to be chal-lenging to obtain the correct pion charge radius from thelattice results [50–62]. The lattice calculations at the un-physical quark masses lead to smaller pion charge radiusthan the experimental results. If the monopole form (8)could describe the pion form factor for all Q the pioncharge radius would be related to the monopole mass as r π = √ M . (10)Therefore, it is convenient to represent the lattice resultsin terms of the effective charge radius defined as [50] r eff ( Q ) = 6(1 /F π ( Q ) − Q . (11)In Fig. 6 we show the effective radius for a = 0 . m valπ = 300 MeV. We see from the figure that r eff isconstant as function of Q for all three lattice spacings.For the smallest lattice spacing, a = 0 .
04 fm the resultson the effective radius are Q -independent for Q as highas 1 . . This is consistent with earlier findings [50].We also clearly see the quark mass dependence of r eff .The effective radius is smaller for the heavier pion massas expected. Comparing the results at a = 0 .
06 fm and a = 0 .
04 fm we see no clear lattice spacing dependenceof r eff . Therefore, we conclude that for a = 0 .
06 fm thediscretization errors for the pion form factor are smallerthan the estimated lattice errors in the range of Q stud-ied by us. Finally, for the two finer lattices we also showthe results from the calculations using Breit frame, whichagree with the non-Breit frame results.While the monopole Ansatz seems to describe the pionform factor very well, it is useful to consider an alterna-tive parameterizations of the pion form factor. An alter-native way to fit the form factors is the model indepen-dent method called the z -expansion [73]. Here the formfactor is written as F π ( Q ) = k max (cid:88) k =0 a k z k z ( t, t cut , t ) = √ t cut − t − √ t cut − t √ t cut − t + √ t cut − t (12)where t = − Q , a k are the fit parameters with con-strain condition F π ( Q = 0) = 1, and t cut = 4 m π isthe two-pion production threshold. Furthermore, t ischosen to be the optimal value t opt0 ( Q ) = t cut (1 − (cid:112) Q /t cut ) to minimize the maximum value of | z | ,with Q the maximum Q used for the fit. We useAIC model selection rules to determine k max , which are1 for a = 0.06 fm, and 2 for a = 0.04, 0.076 fm data andfor the Q under consideration. The z expansion resultsare also shown in Fig. 5 and appear to agree well withthe monopole fit, but for larger Q it has larger errors,c.f. Fig. 5. We also show the fits with the z -expansionin Fig. 6. From this figure we see that this fit workswell also for the valence pion mass of 300 MeV and nat-urally reproduces the Q independence of the effectiveradii. To better understand the quark mass dependenceof the pion form factor as well to facilitate the compar-ison with the experimental results in Fig. 7 we show allthe results for the pion form factor in terms of the ef-fective radius r eff ( Q ). We see that the effective radiusobtained for the physical pion mass is clearly larger thanthe one obtained for m valπ = 300 MeV and is much closerto the CERN data. Furthermore, the fits of r eff for m valπ = 300 MeV for the two lattice spacing agree withinerrors. While the individual lattice data and the CERNdata appear to agree within errors we also see from thefigure that there is a tendency for the CERN data tolie higher than the lattice data. This leads so a slightdifference in the pion charge radius as discussed below.The z -expansion provides an model independent wayto obtain the pion charge radius. In Table III we show Q (GeV ) r e ff ( Q ) [ f m ] a = 0.076 fm r M r Z n z =1 n z =2 n z =3 n z =4 Q (GeV ) r e ff ( Q ) [ f m ] a = 0.06 fm r M r Z n z =0 n z =1 n z =2 n z =3 Breit Q (GeV ) r e ff ( Q ) [ f m ] a = 0.04 fm r M r Z n z =1 n z =2 n z =3Breit FIG. 6. The effective radius as function of Q . The smallererror bars are the statistical errors, while the larger error barsalso include the systematic errors. We show results for a =0 .
076 fm (top panel), a = 0 .
06 fm (middle panel) and a = 0 . z -expansion fit results as a function of Q , whilethe green band is a constant from monopole fit combining datafor different momenta. the charge pion radius for the three lattice spacings usedin our study obtained from the monopole fit and from the z -expansion fit. The statistical error are often smaller forthe monopole fit, but this fit has larger systematic errorscompared to the fit based on z -expansion. Within theestimated errors the two fit forms give consistent results.Thus, the model uncertainty in our determination of thepion charge radius is small. As expected the calculationsfor the heavier quark mass give smaller pion charge ra- Q / m r e ff ( Q ) [ f m ] CERNa=0.076 fm a=0.06 fma=0.04 fm
FIG. 7. The comparison of effective radius between CERNand our lattice data as a function of Q /m π . The bands arethe z expansion fit results of lattice data (blue, green andorange). Data n z (cid:104) r M (cid:105) [fm ] (cid:104) r Z (cid:105) [fm ]a=0.076fm [1,3] 0.402(6)(23) 0.412(10)(19)a=0.06fm [0,3] 0.339(4)(18) 0.329(3)(14)a=0.04fm [1,3] 0.313(5)(27) 0.319(9)(11)TABLE III. The charge radius computed from monopole fit( (cid:104) r M (cid:105) ) and z -expansion fit ( (cid:104) r Z (cid:105) ). The first error is statistical,while the second error is systematic. dius. As our final estimate of the pion charge radius forphysical point we take the result from the z -fit: (cid:104) r π (cid:105) = 0 . , (13)where we added the statistical and systematic errors inquadrature. This result is one sigma lower than the pioncharge radius quoted by Particle Data Group (PDG), (cid:104) r π (cid:105) PDG = 0 . [74], but agrees well with theHPQCD determination that uses HISQ action both inthe sea and the valence sector in 2+1+1 flavor QCD, (cid:104) r π (cid:105) = 0 . [62]. The most precise latticedetermination of the pion charge radius in 2+1 flavorQCD using overlap action in the valence sector and do-main wall action in sea has (cid:104) r π (cid:105) = 0 . [61]. The 2+1 flavor domain wall calculation gives (cid:104) r π (cid:105) = 0 . [60]. Finally, the other 2+1flavor lattice determinations of the pion charge radiushave significantly larger errors [58, 59]. Since both ourand HPQCD’s calculations give lower value of the pioncharge radius compared to the PDG value, while calcu-lations with chiral quarks agree very well with with thePDG result we may speculate whether the uncontrolledeffects of partial quenching are responsible for this. Cal-culations at smaller lattice spacing will be needed to clar-ify this issue. VI. CONCLUSIONS
In this paper we studied the pion form factor in 2+1flavor lattice QCD using three lattices spacings a =0 . a = 0 .
06 and a = 0 .
04 fm. The calculations onthe coarsest lattice have been performed with the physi-cal value of the quark masses, while for the finer two lat-tices the valence pion mass was 300 MeV. We have foundthat the pion form factor is very sensitive to the quarkmass, as expected. We showed that lattice discretizationeffects are quite small for lattice spacings smaller than0 .
06 fm. For the physical quark masses our lattice re-sults on the pion form factor appear to agree with theexperimental determinations. Unlike other lattice stud-ies we also considered highly boosted pions in the initialstate using momentum boosted Gaussian sources. In ad-dition we performed calculations also in the Breit frame.We demonstrated that the calculations of the pion formfactor performed at different momenta of the pion as wellas in the Breit frame give consistent results. This is veryimportant for extending the calculations to pion GPDs.An important outcome of our analysis is that themonopole Ansatz can describe the pion form factor inlarge range of Q , up to Q = 1 . . In the fu-ture it will be important to extend the calculations toeven higher momentum transfer given the experimentalefforts in Jlab and EIC. To do this we should use boostedsources that also depend on the value of Q . At presentthe momentum boost was optimized only according tothe pion momentum in the initial state.From the low Q dependence of the pion form fac-tor we determined the pion charge radius, which is onesigma lower that the experimental result. We speculated,whether this is due to the effect of partial quenching. Tofully resolve this issue calculations at smaller lattice spac-ing with the physical value of the pion masses are needed. AKNOWLEDGEMENTS
This material is based upon work supported by: (i)The U.S. Department of Energy, Office of Science, Of-fice of Nuclear Physics through the Contract Nos. DE-SC0012704 and DE-AC02-06CH11357; (ii) The U.S. De-partment of Energy, Office of Science, Office of NuclearPhysics and Office of Advanced Scientific ComputingResearch within the framework of Scientific Discoverythrough Advance Computing (ScIDAC) award Comput-ing the Properties of Matter with Leadership ComputingResources; (iii) X.G. is partially supported by the NSFCGrant Number 11890712. (iv) N.K. is supported by Jef-ferson Science Associates, LLC under U.S. DOE ContractNo. DE- AC05-06OR23177 and in part by U.S. DOEgrant No. DE-FG02-04ER41302. (v) S.S. is supported bythe National Science Foundation under CAREER AwardPHY- 1847893 and by the RHIC Physics Fellow Programof the RIKEN BNL Research Center. (vi) This researchused awards of computer time provided by the INCITE0 Q (GeV ) r e ff ( Q ) [ f m ] a = 0.076 fm r M r Z n z =1 n z =2 n z =3 n z =4 FIG. 8. Similar plot as Fig. 6 for a = 0 .
076 fm ensemble butusing constant Z − V for renormalization is shown. and ALCC programs at Oak Ridge Leadership Comput-ing Facility, a DOE Office of Science User Facility op- erated under Contract No. DE-AC05-00OR22725. (vii)Computations for this work were carried out in part onfacilities of the USQCD Collaboration, which are fundedby the Office of Science of the U.S. Department of En-ergy. Appendix A: Discretization errors
As is shown in Fig. 4, there are (cid:46)
2% discretizationeffects of Z − V ( P i ) = h B ( P i , P i ), and we choose to ob-tain the renormalized pion form factor by simply dividingthe h B ( P f , P i ) by h B ( P i , P i ). To estimate the impactof the discretization errors to the charge radius, insteadwe can divide the h B ( P f , P i ) by a constant Z − V such as Z − V (0.25 GeV) of a = 0.076 fm ensemble. The effectiveradius for a = 0.076 fm ensemble is shown in Fig. 8, and inthis case we estimate the charge radius from monopole fitand z -expansion fit as 0.406(6)(25) fm and 0.418(11)(23)fm , which shift (cid:46)
2% but are consistent with the esti-mates in Table III. [1] Z.-F. Cui, M. Ding, F. Gao, K. Raya, D. Binosi,L. Chang, C. D. Roberts, J. Rodr´ıguez-Quintero,and S. M. Schmidt, Eur. Phys. J. A , 5 (2021),arXiv:2006.14075 [hep-ph].[2] C. D. Roberts and S. M. Schmidt, Eur. Phys. J. ST ,3319 (2020), arXiv:2006.08782 [hep-ph].[3] P. C. Barry, N. Sato, W. Melnitchouk, and C.-R. Ji,Phys. Rev. Lett. , 152001 (2018), arXiv:1804.01965[hep-ph].[4] I. Novikov et al. , Phys. Rev. D , 014040 (2020),arXiv:2002.02902 [hep-ph].[5] J.-H. Zhang, J.-W. Chen, L. Jin, H.-W. Lin, A. Sch¨afer,and Y. Zhao, Phys. Rev. D , 034505 (2019),arXiv:1804.01483 [hep-lat].[6] R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu,and D. G. Richards, Phys. Rev. D , 074507 (2019),arXiv:1901.03921 [hep-lat].[7] B. Jo´o, J. Karpie, K. Orginos, A. V. Radyushkin, D. G.Richards, R. S. Sufian, and S. Zafeiropoulos, Phys. Rev.D , 114512 (2019), arXiv:1909.08517 [hep-lat].[8] R. S. Sufian, C. Egerer, J. Karpie, R. G. Edwards, B. Jo´o,Y.-Q. Ma, K. Orginos, J.-W. Qiu, and D. G. Richards,Phys. Rev. D , 054508 (2020), arXiv:2001.04960 [hep-lat].[9] T. Izubuchi, L. Jin, C. Kallidonis, N. Karthik, S. Mukher-jee, P. Petreczky, C. Shugert, and S. Syritsyn, Phys. Rev.D , 034516 (2019), arXiv:1905.06349 [hep-lat].[10] X. Gao, L. Jin, C. Kallidonis, N. Karthik, S. Mukherjee,P. Petreczky, C. Shugert, S. Syritsyn, and Y. Zhao, Phys.Rev. D , 094513 (2020), arXiv:2007.06590 [hep-lat].[11] X. Ji, Phys. Rev. Lett. , 262002 (2013),arXiv:1305.1539 [hep-ph].[12] X. Ji, Sci. China Phys. Mech. Astron. , 1407 (2014),arXiv:1404.6680 [hep-ph].[13] A. V. Radyushkin, Phys. Rev. D , 034025 (2017),arXiv:1705.01488 [hep-ph]. [14] K. Orginos, A. Radyushkin, J. Karpie, andS. Zafeiropoulos, Phys. Rev. D , 094503 (2017),arXiv:1706.05373 [hep-ph].[15] V. Braun and D. M¨uller, Eur. Phys. J. C , 349 (2008),arXiv:0709.1348 [hep-ph].[16] Y.-Q. Ma and J.-W. Qiu, Phys. Rev. D , 074021(2018), arXiv:1404.6860 [hep-ph].[17] Y.-Q. Ma and J.-W. Qiu, Phys. Rev. Lett. , 022003(2018), arXiv:1709.03018 [hep-ph].[18] K. Cichy and M. Constantinou, Adv. High Energy Phys. , 3036904 (2019), arXiv:1811.07248 [hep-lat].[19] Y. Zhao, Int. J. Mod. Phys. A , 1830033 (2019),arXiv:1812.07192 [hep-ph].[20] A. V. Radyushkin, Int. J. Mod. Phys. A , 2030002(2020), arXiv:1912.04244 [hep-ph].[21] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao,(2020), arXiv:2004.03543 [hep-ph].[22] C. Best, M. Gockeler, R. Horsley, E.-M. Ilgenfritz,H. Perlt, P. E. L. Rakow, A. Schafer, G. Schierholz,A. Schiller, and S. Schramm, Phys. Rev. D , 2743(1997), arXiv:hep-lat/9703014.[23] M. Guagnelli, K. Jansen, F. Palombi, R. Petronzio,A. Shindler, and I. Wetzorke (Zeuthen-Rome (ZeRo)),Eur. Phys. J. C , 69 (2005), arXiv:hep-lat/0405027.[24] S. Capitani, K. Jansen, M. Papinutto, A. Shindler, C. Ur-bach, and I. Wetzorke, Phys. Lett. B , 520 (2006),arXiv:hep-lat/0511013.[25] A. Abdel-Rehim et al. , Phys. Rev. D , 114513(2015), [Erratum: Phys.Rev.D 93, 039904 (2016)],arXiv:1507.04936 [hep-lat].[26] M. Oehm, C. Alexandrou, M. Constantinou, K. Jansen,G. Koutsou, B. Kostrzewa, F. Steffens, C. Urbach,and S. Zafeiropoulos, Phys. Rev. D , 014508 (2019),arXiv:1810.09743 [hep-lat].[27] C. Alexandrou, S. Bacchio, I. Cloet, M. Constantinou,K. Hadjiyiannakou, G. Koutsou, and C. Lauer, Phys. Rev. D , 014508 (2021), arXiv:2010.03495 [hep-lat].[28] J. Dudek et al. , Eur. Phys. J. A , 187 (2012),arXiv:1208.1244 [hep-ex].[29] Y.-S. Liu, W. Wang, J. Xu, Q.-A. Zhang, J.-H. Zhang,S. Zhao, and Y. Zhao, Phys. Rev. D , 034006 (2019),arXiv:1902.00307 [hep-ph].[30] J.-W. Chen, H.-W. Lin, and J.-H. Zhang, Nucl. Phys. B , 114940 (2020), arXiv:1904.12376 [hep-lat].[31] H.-W. Lin, (2020), arXiv:2008.12474 [hep-ph].[32] C. Alexandrou, K. Cichy, M. Constantinou, K. Had-jiyiannakou, K. Jansen, A. Scapellato, and F. Steffens,Phys. Rev. Lett. , 262001 (2020), arXiv:2008.10573[hep-lat].[33] E. B. Dally et al. , Phys. Rev. D , 1718 (1981).[34] E. B. Dally et al. , Phys. Rev. Lett. , 375 (1982).[35] S. R. Amendolia et al. , Phys. Lett. B , 116 (1984).[36] S. Amendolia et al. (NA7), Nucl. Phys. B , 168(1986).[37] C. J. Bebek, C. N. Brown, M. Herzlinger, S. D. Holmes,C. A. Lichtenstein, F. M. Pipkin, S. Raither, and L. K.Sisterson, Phys. Rev. D , 25 (1976).[38] C. J. Bebek et al. , Phys. Rev. Lett. , 1326 (1976).[39] C. J. Bebek et al. , Phys. Rev. D , 1693 (1978).[40] H. Ackermann, T. Azemoon, W. Gabriel, H. D. Mertiens,H. D. Reich, G. Specht, F. Janata, and D. Schmidt, Nucl.Phys. B , 294 (1978).[41] P. Brauel, T. Canzler, D. Cords, R. Felst, G. Grind-hammer, M. Helm, W. D. Kollmann, H. Krehbiel, andM. Schadlich, Z. Phys. C , 101 (1979).[42] J. Volmer et al. (Jefferson Lab F(pi)), Phys. Rev. Lett. , 1713 (2001), arXiv:nucl-ex/0010009.[43] V. Tadevosyan et al. (Jefferson Lab F(pi)), Phys. Rev. C , 055205 (2007), arXiv:nucl-ex/0607007.[44] T. Horn et al. (Jefferson Lab F(pi)-2), Phys. Rev. Lett. , 192001 (2006), arXiv:nucl-ex/0607005.[45] H. P. Blok et al. (Jefferson Lab), Phys. Rev. C , 045202(2008), arXiv:0809.3161 [nucl-ex].[46] G. Huber et al. (Jefferson Lab), Phys. Rev. C , 045203(2008), arXiv:0809.3052 [nucl-ex].[47] G. P. Lepage and S. J. Brodsky, Phys. Lett. B , 359(1979).[48] G. Martinelli and C. T. Sachrajda, Nucl. Phys. B ,865 (1988).[49] T. Draper, R. M. Woloshyn, W. Wilcox, and K.-F. Liu,Nucl. Phys. B , 319 (1989).[50] D. Br¨ommel et al. (QCDSF/UKQCD), Eur. Phys. J. C , 335 (2007), arXiv:hep-lat/0608021.[51] R. Frezzotti, V. Lubicz, and S. Simula (ETM), Phys.Rev. D , 074506 (2009), arXiv:0812.4042 [hep-lat].[52] S. Aoki et al. (JLQCD, TWQCD), Phys. Rev. D ,034508 (2009), arXiv:0905.2465 [hep-lat]. [53] B. B. Brandt, A. J¨uttner, and H. Wittig, JHEP , 034(2013), arXiv:1306.2916 [hep-lat].[54] C. Alexandrou et al. (ETM), Phys. Rev. D , 014508(2018), arXiv:1710.10401 [hep-lat].[55] F. D. R. Bonnet, R. G. Edwards, G. T. Fleming,R. Lewis, and D. G. Richards (Lattice Hadron Physics),Phys. Rev. D , 054506 (2005), arXiv:hep-lat/0411028.[56] P. A. Boyle, J. M. Flynn, A. Juttner, C. Kelly, H. P.de Lima, C. M. Maynard, C. T. Sachrajda, and J. M.Zanotti, JHEP , 112 (2008), arXiv:0804.3971 [hep-lat].[57] O. H. Nguyen, K.-I. Ishikawa, A. Ukawa, and N. Ukita,JHEP , 122 (2011), arXiv:1102.3652 [hep-lat].[58] H. Fukaya, S. Aoki, S. Hashimoto, T. Kaneko, H. Mat-sufuru, and J. Noaki, Phys. Rev. D , 034506 (2014),arXiv:1405.4077 [hep-lat].[59] S. Aoki, G. Cossu, X. Feng, S. Hashimoto, T. Kaneko,J. Noaki, and T. Onogi (JLQCD), Phys. Rev. D ,034504 (2016), arXiv:1510.06470 [hep-lat].[60] X. Feng, Y. Fu, and L.-C. Jin, Phys. Rev. D , 051502(2020), arXiv:1911.04064 [hep-lat].[61] G. Wang, J. Liang, T. Draper, K.-F. Liu, and Y.-B.Yang (chiQCD), (2020), arXiv:2006.05431 [hep-ph].[62] J. Koponen, F. Bursa, C. Davies, R. Dowdall,and G. Lepage, Phys. Rev. D , 054503 (2016),arXiv:1511.07382 [hep-lat].[63] J. Bijnens, G. Colangelo, and P. Talavera, JHEP , 014(1998), arXiv:hep-ph/9805389.[64] A. Hasenfratz and F. Knechtli, Phys. Rev. D , 034504(2001), arXiv:hep-lat/0103029.[65] A. Bazavov et al. (HotQCD), Phys. Rev. D , 094503(2014), arXiv:1407.6387 [hep-lat].[66] A. Bazavov et al. , Phys. Rev. D , 094510 (2019),arXiv:1908.09552 [hep-lat].[67] A. Bazavov et al. (MILC), PoS LATTICE2010 , 074(2010), arXiv:1012.0868 [hep-lat].[68] E. Shintani, R. Arthur, T. Blum, T. Izubuchi, C. Jung,and C. Lehner, Phys. Rev. D , 114511 (2015),arXiv:1402.0244 [hep-lat].[69] X. Gao, N. Karthik, S. Mukherjee, P. Petreczky, S. Syrit-syn, and Y. Zhao, (2021), arXiv:2101.11632 [hep-lat].[70] S. Capitani et al. , Nucl. Phys. B Proc. Suppl. , 294(1999), arXiv:hep-lat/9809172.[71] W. Wilcox, T. Draper, and K.-F. Liu, Phys. Rev. D ,1109 (1992), arXiv:hep-lat/9205015.[72] H. B. O’Connell, B. C. Pearce, A. W. Thomas, andA. G. Williams, Phys. Lett. B , 14 (1995), arXiv:hep-ph/9503332.[73] G. Lee, J. R. Arrington, and R. J. Hill, Phys. Rev. D , 013013 (2015), arXiv:1505.01489 [hep-ph].[74] P. A. Zyla et al. (Particle Data Group), PTEP2020