Lattice QCD calculation of the pion charge radius using a model-independent method

Xu Feng, 2, 3, 4, ∗ Yang Fu, and Lu-Chang Jin 6, † School of Physics, Peking University, Beijing 100871, China Collaborative Innovation Center of Quantum Matter, Beijing 100871, China Center for High Energy Physics, Peking University, Beijing 100871, China State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, China Department of Physics, University of Connecticut, Storrs, CT 06269, USA RIKEN-BNL Research Center, Brookhaven National Laboratory, Building 510, Upton, NY 11973 (Dated: November 12, 2019)

Introduction. -In particle physics, hadron is a bound state of quarks and gluons, which are held together by the strong interaction force. Different from a point-like particle, hadron has a rich internal structure. One intrinsic property of a hadron is its charge radius, which corresponds to the spatial extent of the distribution of the hadron's charge. The accurate determination of the charge radius not only leaves us useful information on the size and the structure of the hadron, but also provides crucial precision tests of the Standard Model at low energy. It is of special importance in resolving the proton radius puzzle [1], where two recent experiments report results which agree with the previous ones obtained by spectroscopy of muonic hydrogen [2,3] and represent a decisive step towards solving the puzzle for a decade.
In the theoretical study, the charge radius of the hadron is essentially a non-perturbative quantity. It is highly appealing to have a reliable calculation of this quantity with robust uncertainty estimate using lattice QCD. The traditional approach for determining the charge radius on the lattice involves the extrapolation of the expression (F (q 2 ) − 1) q 2 to zero momentum transfer q 2 = 0, where F (q 2 ) is the vector form factor. The choices of the fit ansatz and fitting window would inevitably bring systematic uncertainties from modeling the momentum dependence of F (q 2 ). To reduce such uncertainties, twisted boundary conditions [4,5] and momentum derivatives of quark propagators [6,7] are proposed and used.
In this work, we use an approach to directly determine the charge radius without the momentum extrapolations. The method is easy to implement on the lattice calculation with no requirement on twisted boundary conditions or sequential-source propagators containing momentum derivatives. As an example, we apply the method to the calculation of pion's charge radius. This quantity has been determined by various groups using the traditional method [8,9], twisted boundary conditions [10][11][12][13][14][15][16] as well as the momentum extrapolations from the timelike region [17][18][19]. In our study we find that the statistical errors can be reduced to 2.1% -4.6% at the physical point. At m π ≈ 340 MeV, we obtain a sub-percent statistical uncertainty 0.8%, which is about 6 -13 times smaller than that from the previous calculations at the similar pion masses [9][10][11][12][13][14].
Upon finishing this work, we note that a similar idea has been proposed by C. C. Chang et. al. in Ref. [20] to calculate the proton charge radius. The earlier work along this direction can be traced back to mid 90's to calculate the slope of the Isgur-Wise function at zerorecoil [21]. [22] When utilizing this method, we find that it cannot be used directly in the calculation of the pion charge radius as it suffers from significant finite-volume effects. We therefore develop the techniques to solve the problems, which are described in the following context.
Charge radius in the continuum theory. -We start with a Euclidean hadronic function in the infinite volume where π( ⃗ 0)⟩ is a pion initial state carrying zero spatial momentum. J µ is an electromagnetic vector current. φ is an interpolating operator, which can annihilate a pion state. It can be chosen as e.g. a pseduoscalar operator uγ 5 d or an axial vector currentūγ µ γ 5 d. In this study we use the hadronic function At large time t, H(x) is saturated by the single pion state where the symbol ≐ denotes the omission of the excited states. The decay constant f π ≈ 130 MeV is from PCAC relation ⟨0 A 4 (0) π(⃗ p)⟩ = Ef π and E = m 2 π + ⃗ p 2 the pion's energy. The pion form factor F π (q 2 ) can be extracted from the matrix element ⟨π(⃗ p) J 4 (0) π( ⃗ 0)⟩ = (E + m π )F π (q 2 ), with q 2 = (E − m π ) 2 − ⃗ p 2 . In the Taylor expansion c 0 = 1 is required by the charge conservation and c 1 is related to the mean-square charge radius via c 1 = m 2 π 6 ⟨r 2 π ⟩. The spatial Fourier transform of Eq. (2) yields The derivative ofH(t, ⃗ p) at ⃗ p 2 = 0 leads to withH π (t, ⃗ 0) = f π m π e −mπt . Combining Eq. (5) and (6), one can determine c 1 using H(x) as input through Charge radius on the lattice. -In a realistic lattice QCD calculation with a lattice size of ∼ 5 fm, the finite volume truncation effects are very large as at the edge of box the integrand in Eq. (5) scales as (7) are too sloppy to be used in a precision calculation. On the lattice with a size L and a lattice spacing a, the hadronic function H (L) (x) is approximated by where Γ indicates a set of discrete momenta with ⃗ x ∈ L 3 running through x i = −L 2, −L 2+a, ⋯, L 2−a for i = 1, 2, 3.
Considering the lattice discretization, we propose to use the lattice dispersion relationÊ 2 =m 2 + ∑ ip 2 i with aÊ = 2 sinh(aE 2), am = 2 sinh(am π 2) and ap i = 2 sin(ap i 2). The notation ∑ i indicates the summation over all spatial directions. We further adopt the latticemodified relations with ap 0 = 2 sinh(ap 0 2) and aẼ = sinh(aE). The square of momentum transferq 2 is given byq 2 = (Ê−m) 2 −∑ ip 2 i . As a next step, we construct a ratio with the coefficients β (L) Here we have used the relations in Eq. (10). The value of c 1 can be approximated by . Note that when a → 0 and L → ∞, all the coefficients β (L) n for n ≥ 2 vanish as in Eq. (7). We can consider the contamination from c n≥2 terms as the systematic effects, which are well under control by using the fine lattice spacings and large volumes. Therefore Eq. (13) provides a direct way to calculate the pion charge radius using the lattice quantity H (L) (x) as input.
Error reduction. -The hadronic function H (L) (x) is exponentially suppressed at large ⃗ x and thus the lattice data near the boundary of the box mainly contribute to the noise rather than the signal. To reduce the statistical error, we introduce an integral range ξL with ξ ≤ √ 3 2 (For ξ ≤ 1 2 the range has a spherical shape.) and define which are related to c n through To remove the systematic contamination from the c 2 term, we use both D (L,ξ) 1 (t) and D (L,ξ) 2 (t) to construct the ratio R (L,ξ) (t) where the parameters f 1 , f 2 and h are chosen to remove the c 0 and c 2 terms. Namely, we impose three conditions Under these conditions R (L,ξ) (t) is given by Although R (L,ξ) (t) still receives the contamination from c n≥3 terms, we expect these effects are negligibly small. In the vector meson dominance model, c n is given by where m ρ is the rho meson mass. For n ≥ 3, c n is estimated to be less than 0.1% of c 1 .
Correlator construction. -We use four gauge ensembles at the physical pion mass together with an additional one at m π ≈ 340 MeV, generated by the RBC and UKQCD Collaborations using domain wall fermion [23,24]. The ensemble parameters are shown in Table I. We calculate the correlation function ⟨A 4 (x)J 4 (0)φ † π (−t π )⟩ using wall-source pion interpolating operators φ † π , which have a good overlap with the π ground state. We find the ground-state saturation for t π ≳ 1 fm. In practise the values of t π are chosen conservatively as shown in Table I.  Table I. Ensembles used in this work. For each ensemble we list the pion mass mπ, the spatial and temporal extents, L and T , the inverse of lattice spacing a −1 , the number of configurations used, N conf , the number of point-source light-quark propagator generated for each configuration, Nr, and the time separation, tπ, used for the π ground-state saturation.
We produce wall-source light-quark propagators on all time slices and point-source ones at N r random spacetime locations {x 0 }. The values of N r are shown in Table I. For each configuration we perform 4 N r measurements of the correlator and obtain an average of where t 0 and t are the time component of x 0 and x, respectively.
The hadronic function H (L) (x) can be obtained from C(x; t π ) through with the factor N π defined as N π = 1 2mπ ⟨π( ⃗ 0) φ † π 0⟩ and Z V A the renormalization factor which converts the local vector/axial-vector current to the conserved one. Note that the overall factor N −1 π Z V Z A e mπtπ cancels out when building the ratio R (L,ξ) (t).
Numerical analysis. -The results of R (L,ξ) (t) as a function of t are shown in Fig. 1 for each ensemble. By using ξL = 1.5 fm, we find that the statistical uncertainties of R (L,ξ) (t) are reduced by a factor of 1.3 -1.8 comparing to the results using ξ = √ 3 2 . We expect that the error reduction can be much more significant in the calculation of the nucleon charge radius, where the signalto-noise ratio decreases as e ( 3 2 mπ−m N ) x at large x, with m N the nucleon's mass. At large t, we perform a correlated fit of R (L,ξ) (t) to a constant and determine c 1 . The corresponding results for ⟨r 2 π ⟩ are listed in Table II. To make a comparison, we also calculate the charge radius using the tradition method. We perform the discrete spatial Fourier transform and calculateH (L) (t, ⃗ p) 0.476 (18) Table II. Charge radii ⟨r 2 π ⟩ from the new method by fitting R (L,ξ) (t) to a constant and from the traditional method by using the momentum extrapolation of Fπ(q 2 ). In the last row, the PDG value of ⟨r 2 π ⟩ = 0.434(5) [25] is listed for a comparison. ) as a function ofq 2 m 2 π are shown in the right panel. using where N R = ∑R ∈O h 1 and O h is the full cubic group for all lattice ratotations and reflectionsR. We then construct the ratio with ⃗ p = 2π L ⃗ n for ⃗ n = (0, 0, 1), (0, 1, 1), (1, 1, 1) and (0, 0, 2). Here we use the ensemble 24D-340 with smallest statitical uncertainty as an example and show the t dependence of M (L) (t, ⃗ p) in the left panel of Fig. 2 as well as theq 2 dependence of F π (q 2 ) in the right panel. We perform a correlated fit of the lattice data to a polynomial function The fitting results are shown in Table II. These results are consistent with the ones from the new method, but the errors are 1.5 -1.9 times larger. Therefore, we use ⟨r 2 π ⟩ from the new method in the following analysis. Systematic effects. -To examine the finite-volume effects, we use the ensembles, 24D and 32D, which have the same pion mass and lattice spacing but different lattice sizes, L = 4.7 and 6.2 fm. For these two ensembles, the results for ⟨r 2 π ⟩ are very consistent, suggesting that the finite-volume effects are mild. This is not surprising since the coefficients β (L,ξ) k,n in Eq. (17) are introduced to treat the finite-volume effects properly.
We have four ensembles nearly at the physical pion mass. The remaining systematic effects from the unphysical pion mass are small and can be corrected by using the information of the fifth ensemble, 24D-340, at m π ≈ 340 MeV. We adopt the chiral extrapolation formula [26,27] (4πF0) 2 + const with the constant term including the possible lattice artifacts. We fix F 0 = 87 MeV, a value estimated by using F π = 92.2(1) MeV and F π F 0 = 1.062(7) [27], and use the ensembles 24D, 32D and 24D-340 with the same lattice spacing to study the pion mass dependence. By extrapolating to the physical point, ⟨r 2 π ⟩ for ensembles 24D, 32D, 32D-fine and 48I are shifted by 0.2%, 0.3%, 0.5%, −0.1%, respectively. These changes are very small compared the statistical errors.
The largest systematic uncertainties in our study arise from the lattice discretization effects. The values of ⟨r 2 π ⟩ for 24D and 32D are 13% larger than that for 32D-fine, suggesting a large lattice artifacts. Unfortunately, the result from 48I cannot be used in the continuum extrapolation together with 24D, 32D and 32D-fine ones, as the 48I ensemble is simulated with Iwasaki gauge action, while the other three use Iwasaki+DSDR action. Considering the fact that 48I has the finest lattice spacing, we quote its value of ⟨r 2 π ⟩ as the final result and attribute to it a ∼ 3% discretization error by an order counting O((aΛ QCD ) 2 ) with Λ QCD = 300 MeV Conclusion.
-We have used a model-independent method to calculate the hadron's charge radius using lattice QCD. It has three peculiar features.
1. Simplicity: The method does not require the additional quark propagator inversion on twisted boundary conditions or sequential-source propagators with momentum derivatives. It does not require the modeling of the momentum dependence of the form factor. The charge radius can be simply extracted from R (L,ξ) (t) at large time separation to avoid the excited-state contamination.
2. Flexibility: In the whole calculation, it only requires the generation of the wall-source and point-source propagators. These propagators can be used to calculate other correlation functions in the future projects. Besides, the hadronic function ⟨0 A µ (x)J ν (0) π( ⃗ 0)⟩ constructed in this study can be used for other relevant physics processes, such as the radiative corrections to the pion's decay.
3. Precision: The statistical uncertainties of ⟨r 2 π ⟩ from the new method are about 1.5 -1.9 smaller times than that from the traditional method. We expect the method is more efficient in the nucleon sector where the hadronic function near the boundary of box contribute significant noise. Besides for the reduction of the statistical uncertainty, the model dependence from the choices of the fit ansatz is also avoided by using the new method.
In this study, we find that the largest source of the uncertainty is from the lattice discretization. This can be controlled by using gauge configurations with finer lattice spacings and performing the continnum extrapolations. With the developments of supercomputers, technologies as well as the new ideas and methods, we can foresee that in the near future lattice QCD calculations can provide the determinations of ⟨r 2 π ⟩, which has the similar precision as the current PDG value or even surpasses it. These developments also shed the light on precise determinations of the proton charge radius from first-principle theory that can distinguish between the conflicting experimental values.