$I=0$ $\pi\pi$ $s$-wave scattering length from lattice QCD

We deliver lattice results for the $I=0$ $\pi\pi$ elastic $s$-wave scattering length calculated with the MILC $N_f=3$ flavors of the Asqtad-improved staggered fermions. The scattering phase shifts are determined by L\"uscher's formula from the energy-eigenvalues of $\pi\pi$ systems at one center of mass frame and four moving frames using the moving wall source technique. Our measurements are good enough to resolve the scattering length $a$ and effective range $r$, moreover, it allows us to roughly estimate the shape parameter $P$. Using our lattice results, the scattering length $a$ and effective range $r$ at the physical point are extrapolated by chiral perturbation theory. Our results are reasonably consistent with the Roy equation determinations and the newer experimental data. Numerical computations are carried out with two MILC fine ($a\approx0.09$~fm, $L^3 \times T = 40^3\times 96$) and one MILC superfine ($a\approx0.06$~fm, $L^3 \times T = 48^3\times 144$) lattice ensembles at three pion masses of $m_\pi\sim247~{\rm MeV}$, $249~{\rm MeV}$, and $314~{\rm MeV}$, respectively.

It is well-known that the I = 0 ππ channel harbors the lowest resonance: the scalar σ or f 0 (500) meson, which have been recently reconfirmed with dispersive analyses and new experimental data [17]. Hence, it is highly desirable to investigate the isospin-0 ππ interaction properties directly from lattice QCD.
However, only a few of lattice studies in the isospin-0 channel are reported so far. Kuramashi et al. carried out a pioneering work for the isospin-0 channel with moving wall source technique [18,19]. A first attempt to extract the sigma from lattice was made in Refs. [20,21]. With vacuum diagram, we attempted to crudely calculate the isospin-0 ππ scattering, and estimated the value of ℓ I=0 ππ , which is a LEC appearing in the NLO χPT expression of I = 0 ππ scattering length [23][24][25]. Liu et al. study the isospin-0 ππ s-wave scattering length from twisted mass lattice QCD, and evaluate ℓ I=0 ππ [26]. A chiral extrapolation of the Hadron Spectrum results is performed in the isoscalar sector [22]. Recently, R. A. Briceno investigated the isoscalar ππ scattering and σ meson resonance [27]. Guo et al. studied the elastic phase-shifts for ππ scattering in the scalar and isoscalar channel [28].
It is very impressive that, as it is demonstrated for the isospin-2 ππ scattering [10], the near threshold behavior of the inverse partial wave amplitude can be exploited to determine the scattering length, effective range and shape parameters, which can be written in terms of a series of threshold parameters which satisfy low-energy theorems imposed by chiral symmetry and encrypt the chief momentum dependence of partial wave amplitude [5,10]. Thus, one can predict the scattering parameters at the physical point by relevant NLO χPT expressions [10]. Lattice predictions of scattering parameters are observed to be consistent with Roy equation determinations [10]. In this work, we will extend this technique to the isospin-0 ππ scattering. As it is shown later, this is not a trivial development since it exhibits some novel features.
Additionally, this approach is only valid in the elastic region, so the pion mass values should be small enough to be below threshold where σ meson approaches stable [26]. Although the strict value of this threshold is not clear, the one-loop inverse amplitude method [29] indicates that the threshold m π < 400 MeV should be secure (see more details in Refs. [30][31][32]). In this work, our lattice studies are calculated at pion masses: 247 MeV, 249 MeV and 314 MeV, respectively. Since pion mass values are all below this threshold, the influence of the σ meson can be reasonably ignored [26]. Therefore, we will only use ππ correlator to calculate the I = 0 ππ elastic s-wave scattering lengths, as it is done in Ref. [26].
According to the rule of thumb discussions [48,49], if we use the fine gauge configurations, employ lattice ensembles with relatively large spatial dimensions L, and sum the correlators over all time slices (96 or 144), the signals are anticipated to be significantly improved [49]. Consequently, the signals of vacuum diagram are found to be remarkably improved as compared with our previous works in Refs. [23,24]. It allows us to not only measure the scattering length, but also explore the effective range. The chiral extrapolations of the scattering length m π a I=0 0 is performed using NLO χPT. Extrapolated to the physical value of m π /f π , our final results give rise to m π a I=0 0 = 0.217(9)(5), ℓ I=0 ππ = 45.6(7.6)(3.8), which are in reasonable agreement with the recent experimental and theoretical determinations as well as the lattice calculations.
Most of all, after the chiral extrapolations of the effective range m π r to the physical point, we can probe its lattice result as m π r = −6.07 (44) (36), which is also in fair accordance with the Roy-equation determination [5,6].
This paper is organized as follows. The Lüscher's finite volume method, lattice setup and the computation of the finite volume spectrum of ππ system are discussed in Sec. II. The lattice results are given in Sec. III, together with relevant fits, which are employed to determine the threshold parameters and the effective range parameters. In Sec. IV, a summary of the relevant χPT formulas at NLO and the chiral extrapolation of latticemeasured data are provided. A brief summary and some discussions are shown in Sec. V.

A. Center of mass frame
In the center-of-mass frame, the energy levels of two free pions are provided by where p = 2π L n, and n ∈ Z 3 . The lowest energy E for n = 0 (e.g., n = (0, 0, 1)) is not within the elastic region 2m π < E < 4m π , or beyond the t-channel cut, which starts at k 2 = m 2 π [10]. We should remark at this point that the finite-volume methods are only valid for elastic scattering, consequently, we are only interested in n = 0 for the current study.
Due to the interaction between two pions, the energy levels of ππ system are shifted by the hadronic interaction from E to E, where the dimensionless scattering momentum q ∈ R. The most important irreducible representation is A + 1 . It is the Lüscher formula that relates the energy E to the s-wave ππ scattering phase δ [38,39], where the zeta function is formally defined by The zeta function Z 00 (s; q 2 ) can be efficiently evaluated by the method described in Ref. [50]. We notice an equivalent Lüscher formula has been recently developed in Ref. [52], where the influence of the d-wave mixing with s-wave for the boosted case is found to be small. The lowest energy levels in the center-of-mass frame are below the threshold (namely, k 2 < 0) due to the attractive interaction in the isospin-0 ππ scattering, as it is already noticed in Refs. [18,19,23,24,26]. It should be worthwhile to stress that the scattering phase shift δ(k) in the continuum is solely denoted for k 2 > 0 [38]. As for the case of k 2 < 0, it is usual to usher in a phase σ(k), which is associated with δ(k) via the analytic continuation of tan σ(k) = −i tan δ(k) [38]. Consequently, in the rest of the analysis, it is convenient to always adopt the notation k cot δ(k), as it is already done in Ref. [18,19,23,24,26].

B. Moving frame
Using a moving frame with non-zero total momentum P = (2π/L)d, d ∈ Z 3 , the energy levels of two free pions where p 1 , p 2 denote the three-momenta of the pions, which obey the periodic boundary condition, p 1 = 2π L n 1 , p 2 = 2π L n 2 , n 1 , n 2 ∈ Z 3 , and total momentum P meets P = p 1 + p 2 [40].
In the presence of the interaction between two pions, the energy E CM is where the dimensionless momentum q ∈ R, k = |p|, and p are quantized to the values p = 2π L r, r ∈ P d , and the set P d is where γ −1 is the inverse Lorentz transformation operating in the direction of the center-of-mass velocity v, where p and p ⊥ are the ingredients of p parallel and perpendicular to v, respectively. Using the Lorentz transformation, the energy E CM is connected to the E MF through The first moving frame (MF1) are taken with d = e 3 . We implemented the second moving frame (MF2) with d = e 2 + e 3 . In order to acquire more eigenenergies, we considered the third moving frame (MF3) with d = e 1 + e 2 + e 3 . The fourth moving frame (MF4) with d = 2e 3 is also taken into account. For these evaluations, the scattering phase shifts can be acquired from the total energy of the two-particle system enclosed in a cubic torus with total momentum P = 2π L d by the generalized Lüscher relation [37][38][39][40] where the γ-factor is denoted by γ = E MF /E CM . The most important irreducible representation is A + 1 , which is relevant for two-particle s-wave scattering states in infinite volume.
The evaluation procedure of zeta functions Z d 00 (1; q 2 ) is discoursed in Appendix A of Ref. [50]. As a consistency check, we exploit the formula (5) described in Ref. [10] to calculate the zeta functions Z d 00 (1; q 2 ) as well, both methods are found to arrive at the same results. Moreover, we have verified that two approaches are mathematically equivalent for the choice of Λ = 1 in Ref. [10].

C. ππ correlator function
In this section, we follow original derivations and notations in Refs. [18,19,34] to study ππ scattering of two Nambu-Goldstone pions in Asqtad-improved staggered dynamical fermion formalism. We build the isospin-0 ππ state with following interpolating operator [18,19] with the interpolating pion operators denoted by Note that the operator O I=0 ππ (p, t) belongs to A + 1 . The ππ four-point correlation function with the momentum p can be expressed as, where we usually choose t 1 = 0, t 2 = 1, t 3 = t, and t 4 = t + 1 to prevent the intricate color Fierz rearrangement of the quark lines [18,19] .
The wall sources are represented by short bars, and open circles indicate the wall sinks for local pion operators. 1 In this operator, the momentum p is associated with time t + 1, and there is another similar operator where the momentum p is associated with time t. After averaging over the two operators, we build an operator, which is symmetric by swapping the momentum 0 and p, and preserve the parity symmetry. Please consult the technical issue of constructing a parity conserving operator in Ref. [51]. Note that this parity violation in Eq. (6) vanishes when t + 1 is replaced by t, consequently, it tends to be a small effect. It should be incorporated into more sophisticated lattice study in the future.
In the present study, we employ moving wall source technique to evaluate four quark-line diagrams [18,19]. In our previous studies [24,25], we present a detailed procedure to express these diagrams in the center-of-mass frame [24] with the light quark propagator G [18,19], and the relevant expressions in the moving frame are also provided in Ref. [25]. To write them in the generic frame (i.e., the momenta p = [0, 0, 0], [0, 0, 1], [0, 1, 1], [1,1,1] and [0, 0, 2]), we used an up quark source with 1, and an anti-up quark source with e ip·x (except for V , where we use 1) on each site for two pion creation operator, respectively, then relevant expressions can be written as where we note the fact that for the momenta The combinations of light quark propagator G that we apply for ππ four-point functions are shown in Fig. 1. Note that the vacuum diagram is not accompanied by a vacuum subtraction for non-zero momenta p [25]. It should be worthwhile to mention that the treatment of vacuum part is actually similar to that of the disconnected piece for sigma operator in our previous study [53], where Professor Carleton DeTar has especially designed a FFT algorithm to calculate the disconnected part [53]. In a same manner, the vacuum diagram can be calculated with a FFT algorithm, which is courteously dedicated to Appendix A. In practice, we found that it indeed saves computer resource significantly.
The ππ correlation functions can be expressed in terms of four diagrams [18,19,34], where the staggered-flavor factor N f = 4 is introduced due to the number of tastes natural to the Kogut-Susskind formulation [34]. Using the quadruple root of the fermion determinant, the four-fold degeneracy of the staggered sea quarks is believed to be neatly removed [34]. We should remember that, the contributions of non-Nambu-Goldstone pions in the intermediate states is exponentially suppressed for large t [18,19,34]. Hence, we think that ππ interpolator does not greatly couple to other ππ tastes, and neglect this systematic errors.

D. Lattice Calculation
We employed the MILC gauge configurations with three Asqtad-improved staggered sea quarks [54,55]. The simulation parameters are summarized in Table I. By MILC convention, lattice ensembles are referred to as "fine" for a ≈ 0.09 fm, and "super-fine" for a ≈ 0.06 fm. For easy notation, it is handy to adopt (am l , am s ) to categorize lattice ensembles. We should bear in mind that MILC gauge configurations are generated using the staggered formulation of lattice fermions [56] with the fourth root of fermion determinant [33]. All the gauge configurations were gauge fixed to Coulomb gauge before calculating light quark propagators.
Although it is expensive, the moving wall source technique [18,19] is believed to be able to calculate the relevant correlators with high quality. We extend this method to two-particle system with nonzero momenta to examine the κ, σ, and K ⋆ (892) meson decays [25,57,58], and meson-meson scattering [23,24,47]. From these works, we found that moving-wall source technique can calculate four-point correlators with desirable quality.
To compute ππ correlators, the conjugate gradient method is used to get the matrix element of the inverse Dirac fermion matrix. We compute the correlators on all the T time slices, and explicitly combine theses results, namely, the correlator C ππ (t) is measured by I. Simulation parameters of the MILC lattice ensembles. Lattice dimensions are described in lattice units with spatial (L) and temporal (T ) size. The gauge coupling β is shown in Column 3. The fourth block give bare masses of the light and strange quark masses in terms of am l and ams, respectively. Column 5 gives pion masses in MeV. The lattice spatial dimension (L) in fm and in units of the finite-volume pion mass are given in Column 6 and 7 respectively. We also list the mass ratio mπ/fπ. The values of the calculated ππ correlators for each of the lattice ensembles are shown in Column 9, and the last Column gives the number of gauge configurations used in this work.
Ensemble After averaging the propagators over all the T values, the statistics are found to be remarkably improved 2 .
According to the analytical arguments in Refs. [19,48] and the semi-empirical discussions in Ref. [49], the noiseto-signal ratio of ππ correlator is improved approximately where L is the lattice spatial dimension, and N slice is the number of the time slices calculated the light propagators for a given lattice ensemble. In the present study, we exploit the MILC lattice ensembles with the relatively large L (40 or 48), and sum ππ correlators over all the time slices (96 or 144), consequently, it is natural that the signals of ππ correlators should be significantly improved. Admittedly, the most efficient way to improve the relevant noise-to-signal ratio is to use the finer gauge configurations or equivalent anisotropic gauge configurations [27].
We compute two-point pion correlators with the zero and none-zero momenta (0 and p) as well, where π is pion point-source operator and W π is pion wall-source operator [33,35]. To simplify notations, the summation over lattice space point in sink is not written out. It should be worthwhile to stress that the summations over all the time slices for π propagators guarantee the extraction of pion mass with high precision (see the semi-empirical discussions in Ref. [49]). Discarding the contributions from the excited states, the pion mass m π and energy E π (p) can be robustly ex-tracted at large t from the two-point pion correlators (10) and (11), respectively [55], where the ellipses indicate the oscillating parity partners, and A π (0) and A π (p) are two overlapping amplitudes, which are used to evaluate the wrap-around contributions for ππ correlators subsequently [59][60][61].
Using the method described in MILC's work [62,63], we can calculate the pion decay constants f π for three MILC ensembles. Our fitted values of the pion decay constants are found to be in well agreement with the same quantities which are professionally computed on the same lattice ensembles by the MILC collaboration [55,62,63]. The MILCs determinations on pion decay constants are directly quoted in this work, and are listed in Table I with m π /f π .

E. Extraction of energies
The energy E ππ of ππ system can be extracted from ππ four-point function as [64,65] for a large t to suppress the excited states, the terms alternating in sign is a symbol of a staggered scheme, and the ellipsis implies the contributions from the excited states which are reduced exponentially. In practice, the pollution due to the "wraparound" effects [59][60][61] should be considered. Actually, we will withdraw the "wraparound" pollution before fitting with this formula, as practiced in our former works [23][24][25].

A. Lattice Phase Shift
In order for a more intuitive presentation of our lattice results, we evaluate the ratios [18,19], where C π (0, t) and C π (1, t+1) are pion correlators with a given momentum [18,19]. We should remark at this point that our visualizations of various diagrams for the ππ correlator is somewhat analogous to those of the Hadron Spectrum Collaboration in Ref. [27], where the time dependence of ππ correlator is weighted by e E0t with E 0 the energy of the lightest state. As a matter of fact, both visualized methods indicate the corresponding energy shifts of ππ system.
In Fig. 2, we display the various contributions to an example ππ correlator for a MILC ensemble (0.0031, 0.031) at P = [0, 0, 0] as the functions of t, which are illustrated as the individual ratios R X (t), X = D, C, R, and V . The vacuum deduction, which is time independent, is carried out with the technique described in Ref. [24]. It is worth stressing once again that the vacuum diagram is not accompanied by a vacuum subtraction for the cases with the non-zero momenta [25].
The ratio values of the direct amplitude R D are quite close to oneness [18,19,[23][24][25][26]. In fact, the I = 0 ππ scattering is perturbative at low momentum and at small light-quark masses, as mandated by χPT. Consequently, in a finite volume, two-pion energies deviate slightly from the noninteracting energies; namely, the sum of the pion masses (or the boosted pion masses for moving frames). Moreover, in lattice practice, the ratios of these energy shifts to total energies are found normally to be about 5% [18,19,[23][24][25][26].
The systematically oscillating behavior of rectangular amplitude is clearly observed, which is a typical feature of the staggered formulation of lattice fermions [64]. In this work, we make use of DeTar's strategy [65] (namely, Eq. (14)) to remove this oscillating contribution. For the vacuum amplitude of this MILC lattice ensemble, we nicely obtain a good signal up to t ∼ 26, beyond which signals are quickly lost (see Fig. 2). We notice that all the four diagrams are measured with good signals, and make a remarkable progress as compared with our previous works [23 -25].
As already explained in previous studies [23][24][25], a persuasive way to process our lattice data is the resort to the "effective energy" plot, which is a variant of the effective mass plot, and very similar to the "effective scattering length" plot [66]. In practice, the isospin-0 ππ four-point functions were fit by altering the minimum fitting distances D min , and putting the maximum distance D max either at T /2 or where the fractional statistical errors exceed about 20% for two sequential time slices [33]. The example "effective energy" plots for the MILC ensemble (0.0031, 0.031) as the functions of D min are illustrated in Fig. 3. For P = [0, 0, 0], the plateau is clearly observed from D min = 8 ∼ 16. We also notice that this plateau is distorted starting from t = 17, since the relevant vacuum contribution becomes noisy after that, as compared with the other contributions, please see Fig. 2 for details.
The energies aE ππ of ππ system are secured from the "effective energy" plots, and usually chosen by looking for a combination of a "plateau" in the effective energy plots as the function of D min , a good confidence level, and D min large enough to suppress the excited states [59,66]. It should be worthwhile to stress once again that the"wraparound" contamination [59][60][61] are removed before fitting by subtracting the latticecalculated wraparound contribution from the relevant correlators. In the present study, our measured quantities from two-point functions are sufficiently precise to allow us to subtract the wraparound contributions, and the unwanted finite-T effects are anticipated to be neatly removed. The details of how to calculate the "wraparound" contamination in center-of-mass frame and the moving frames are provided in our previous work [49].
The fitted values of the energies aE ππ of ππ system, fit range and fit quality (χ 2 /dof) are given in Table II. The fit qualities χ 2 /dof are turned out to be reasonable. Now it is straightforward to substitute these fitted energies aE ππ into Lüscher formula (1) or (5) to secure the relevant values of the k cot δ/m π . All of these values are also summarized in Table II, where the statistical errors of k 2 are calculated from both the statistical errors of the energies aE ππ and pion masses am π . Note that the TABLE II. Summaries of the lattice results for the fitted energies Eππ of the I = 0 ππ system. The third block shows the fitted energies Eππ in lattice units. Column four gives the fitting range, and Column five indicates the number of degrees of freedom (dof) for the fit. The six block is the center-of-mass scattering momentum k 2 in terms of m 2 π , and Column seven gives the values of k cot δ/mπ, which are calculated by Lüscher formula (1)

B. The effective range approximation parameters
As it is shown latter, the effective range approximation is an expansion of the real part of the inverse partial wave scattering amplitude [Re t(k)] −1 in powers of k, the magnitude of the center-of-mass three-momentum of each pion [10], 3 namely, where m π a, m π r, and P are called as the scattering length, effective range, and shape parameter, respectively. Note that m π a and m π r are cited in units of m −1 π , and an alterative way is employed in the present study to denote the dimensionless quantity shape parameter P , as compared with that of the isospin-2 ππ scattering in Ref. [10], where P is scaled with (m π r) 3 . Moreover, there is no minus sign in the first term in Eq. (16) since the scattering length m π a is positive value for the isospin-0 ππ scattering [2,5,10]. For simple notations, ma ≡ ma I=0 0 , similar for m π r and P . As it is pointed out in Ref. [10], the afore-mentioned effective range approximation is believed to be convergent for the energies below the t-channel cut, which starts at k 2 = m 2 π . a Four data are within the region k 2 /m 2 π < 0.5 for this ensemble.
In the present study, we just have five lattice data for each lattice ensemble at disposal. This is mainly due to our currently limited computational resources. Of course, the lack of lattice ensembles with the different spatial extent L for a given pion mass is an another important reason. The lattice-determined values of k cot δ/m π are summarized in Table II, we also exhibit these values in Figs. 4, 5, 6 for the (0.0031, 0.031), (0.0031, 0.0031) and (0.0036, 0.018) ensembles, respectively. Fortunately, the lattice-determined values of k cot δ/m π are turned out to be all within the t-channel cut k 2 = m 2 π . Moreover, we observe that the values of k cot δ/m π are not roughly linear in k 2 during the region k 2 /m 2 π < 1.0, which reflects the fact that the shape parameter P indeed has a impact on the curvature. Actually, according to our analytical discussions in Sec. IV, the second term and third term in Eq. (16) both contribute significantly for the values of k cot δ/m π , 4 Consequently, to be safe, we will include both terms in the analyses of the values of k cot δ/m π throughout the remaining analyses.
So far, the influences of the higher order terms (NNLO, etc) from χPT is quite limited, which indicates the contributions from O(k 6 ) term or higher are not clear [26]. 4 According to the discussions in Sec. IV, we can estimate that the ratio of the effective range mπ r to the shape parameter P at the physical point is −2.27 (39), which partially confirmed the assumption [26] that the contribution of O(k 4 ) term is not big than that of O(k 2 ) term at least within t-channel cut k 2 = m 2 π . Additionally, using just five lattice data for each ensemble, we are not able to effectively consider such high order contribution.
In the region k 2 /m 2 π < 0.5, we only have three or four data at disposal. As is done for the isospin-2 ππ scattering in Ref. [10], the scattering length m π a and effective range m π r are fit (Fit A) with Eq. (16), where the shape parameter P and the other higher orders terms are fixed to be zero. The fitted values of m π a and m π r are given in the third Column of Table III. In the region k 2 /m 2 π < 1, lattice calculations indicates that the curvatures have the quadratic (and higher) dependence on k 2 . Therefore, three leading effective range expansion parameters in Eq. (16)  It is explicit from Table III that the fit parameters of Fit B are reasonable consistent with those of Fit A within the statistical uncertainties. In what follows, we will make use of χPT to predict the scattering parameters at physical pion mass, with the fit parameters listed in Table III as input. We note that our lattice-measured values of k cot δ/m π for each lattice ensemble have relatively large errors, which give rise to the large statistical errors for the extracted quantities (m π a, m π r and P ). On the same time, it leads to the relatively low χ 2 /dof values. The straightforward way to move the χ 2 /dof into a more reasonable range is to use more gauge configu- rations for each lattice ensemble [19,48,49], which is certainly beyond the scope of this work since it requires the huge of the computer resource. We note that some other fitting strategies to improve the statistical errors are discussed in Refs. [68,69], Usually, the truncation of the effective range r is considered to be an important source of systematic error, and in Ref. [24], we only use one point (center of mass P = [0, 0, 0]) to evaluate the scattering length (m π a) 000 , which is also calculated using only one data from center of mass frame in the pioneering works [18,19]. To make this difference of these results more intuitive, we also listed the relative error, R a = m π a − (m π a) 000 m π a × 100% , in Table III for three lattice ensembles. From Table III, we find that the relative errors for three lattice ensembles are not small, which means that this kind of the systematic error can not be ignored. We view this as one of the most important results in the present study.
Considering the relatively strong interaction in the isospin-0 ππ channel, Liu et al. realized that the contribution of O(k 2 ) and higher order terms in the effective range expansion is very important for the reliable lattice calculation of scattering length [26]. Since they have one energy level for each pion mass, the values of the effective range r are directly determined from NLO χPT [26]. Hence, our investigation in this work exactly follows Liu et al.'s work [26], and makes an improvement towards a calculation of I = 0 ππ scattering length by directly estimating the effective range r from lattice simulation.
In addition, we found that the shape parameter P should be also included for the successful fit of the low momentum lattice values of k cot δ/m π during the region k 2 /m 2 π < 1, which are illustrated in Figs. 4, 5, 6, respectively. This will be clearly interpreted by the chiral The same of Fig. 4, but for the (0.0036, 0.018) ensemble.
perturbation prediction at NLO in Sec. IV.
Admittedly, our fitted values of the shape parameter P listed in Table III contain the rather large statistical errors, as a result, it is not convincing to use these data to perform the chiral extrapolation to the physical pion mass. This kinds of work should be waiting for the more robust lattice data in the future. Of course, the most efficient way to improve the statistical errors of the shape parameter P is working on the lattice ensembles with the different size L for a given pion mass, as is done for the isospin-2 ππ scattering in Ref. [10], since it can produce more reliable data in the region k 2 /m 2 π < 1.0. Furthermore, according to the analytical arguments in Refs. [19,48] and rule of thumb discussions in Ref. [49], most straightforward way to improve the noise-to-signal ratios of the relevant correlators is to employ very fine gauge configurations or anisotropic gauge configurations [27]. In addition, if we use lattice ensembles with relatively large L, and sum ππ correlators over all time slices, the signals are anticipated to be significantly improved [49]. With this aim in mind, our ongoing lattice investigation on the isospin-0 ππ scattering will be carried out with MILC ultrafine gauge configuration (a ≈ 0.45 fm and L 3 × T = 64 3 × 192).
About five years ago, when we sincerely request Professor Eulogio Oset to give us some comments for our previous work in Ref. [24], he suggests us to consider the effective range parameter r in the calculation of the I = 0 ππ scattering length 5 . This work is specially to answer his inquiries, and we here admire his sharp insight of the physical essence.

IV. CHIRAL EXTRAPOLATIONS
The low-energy theorems imposed by chiral symmetry indicates that each scattering parameter can be related to the relevant LEC which appears at NLO in χPT [10]. In this section, we will exploit NPLQCD Collaboration's technique in Sec. V of Ref. [10] and follow the original derivations and notations in Refs. [2,5,10] to provide the NLO χPT expressions for the effective range expansion and threshold parameters in the isospin-0 ππ scattering. For easy notation, the Mandelstam variables (s, u, t) are expressed in units of the physical pion mass squared m 2 π .

A. Threshold parameters in χPT
In the elastic region (4 ≤ s ≤ 16), unitarity indicates that the isospin-0 ππ s-wave scattering amplitude t(s) ≡ t I=0 ℓ=0 (s) can be described by real phase shift δ [2,5] According to the discussions in Appendix B, the I = 0 ππ s-wave NLO scattering amplitude can then be expressed in terms of only three independent low-energy constants: C 1 , C 2 , and C 3 [2,5]: where, k = |k| is the magnitude of the center-of-mass three-momentum of each pion, i.e., s = 4(1 + k 2 /m 2 π ). and the constants C i (i = 1, 2, 3) can be described in terms of the renormalized, quark mass independent lowenergy constants ℓ r i (µ = f π,phy ) [2, 5] It should be worthwhile to stress that three constants: 5/(384π 3 ), 247/(576π 3 ), and 35/(48π 3 ) in Eq. (18) will be neatly cancelled out in Taylor expansion of last two terms in Eq. (18) for small k 2 . Note that the renormalization scale µ is fixed to be the physical pion decay constant f π,phy . To express the relevant formulae as a function of m π /f π , we change f π,phy with f π , which only results in the corrections at NNLO [26]. The near threshold (to be specific, k 2 → 0) behavior for the real part of the partial wave amplitude t(k) can be expressed as a power-series expansion in the centerof-mass three-momentum [2,5] Re t(k) = m π a + k 2 b + k 4 c + O(k 6 ), (20) where the threshold parameters a and b are referred to as the scattering length and slope parameter, respectively. Of course, the threshold parameter c is also regarded as an another slope parameter. Note that there is no minus sign before m π a in Eq. (20), which is consistent with the definition of the effective range approximation in Eq. (16). We should remark at this point that three independent constants C 1 , C 2 and C 3 are directly related to the scattering length a, slope parameter b, and slope parameter c, respectively. This is a little bit different from those for the isospin-2 ππ scattering in Ref. [10], where three independent constants C 1 , C 2 and C 4 are directly related to the scattering length a, effective range r, and shape parameter P , respectively.
It is easy to show that the inverse form of real part of scattering amplitude t(k) in Eq. (17) can be written as an elegant form, (21) Plugging Eq. (16) into Eq. (21), and comparing them with the near threshold behavior of the partial wave amplitude t(k) denoted in Eq. (20), the effective range m π r and the shape parameter P can be nicely described only in terms of three threshold parameters [5]: Equations (22) and (23) can be used inversely to secure the slope parameters b and c from the lattice-measured effective range expansion parameters [10].
To simplify the notation, it is convenient to follow the original notation in Ref. [10] to denote z ≡ m 2 π /f 2 π . After expanding the NLO partial wave scattering amplitude denoted in Eq. (18) in powers of k 2 , it is straightforward to obtain the NLO χPT expressions for the threshold parameters: It is interesting and important to notice that no any contributions are accepted from LO χPT for the slope parameter c, as it is already noticed that for the isospin-2 ππ scattering in Ref. [10]. In our previous work [24], we directly make use of the explicit results in Appendix C of Ref. [5] to get the χPT formula for the scattering length m π a, we find it is nice to consistent with the relevant result in Eq. (24), which was recently used to extrapolate the lattice-measured data to the physical pion mass by Liu et al. in Ref. [26]. Note that the constant C 1 is related to the constant ℓ I=0 ππ denoted in Ref. [24] by As it demonstrated in Appendix B, we also nicely reproduce the relevant results in Appendix C of Ref. [5] for slope parameter b. Note that there is no explicit formula for the slope parameter c in Ref. [5]. From Eq. (22), we note that the effective range m π r contains three components. It is easy to show that the second term indeed dominates the effective range m π r for the small z-values, as it is observed in Ref. [26], on the other hand, for the large z-values, the third term absolutely dominates the effective range m π r, and the first term approaches zero. To be safe, we include all the three parts in the present study.
Plugging these formulas for the threshold parameters [namely, the equations (24), (25) and (26)] into the equations (22) and (23), after some strenuous algebraic operations, it is now straightforward to achieve the NLO χPT descriptions for the effective range approximation parameters: where the constants C 4 , C 5 and C 6 are solely denoted in terms of the constants C 1 , C 2 , and C 3 via Note that we should keep last two terms in NLO χPT expressions for m π r in Eq. (28), which naturally come from the third part of m π r denoted in Eq. (22). In the same manner, the last two terms in NLO χPT expressions for m 2 π ar in Eq. (29) are required to hold as well. On the other hand, the relevant two terms for the isospin-2 case can be reasonably overlooked [10]. Similarly, the last six terms in NLO χPT expressions for P in Eq. (30) should be reserved.
We should remark at this point that Eqs. (28), (29), and (30) are valid for the range of interest in the present study, and more terms should be added into these equations for large z-values. Moreover, the factor −9/7 in Eq. (29) indicates the relevant LO χPT prediction.
With the estimated values of C 1 , C 2 , and C 3 in Eq. (37), and using equations (22) and (23), the ratio of the effective range m π r to the shape parameter P at the physical pion mass can be estimated as −2.27 (39). This indicates the second term and third term in Eq. (16) both contribute significantly for the lattice-measured values of k cot δ/m π . Meanwhile, it partially confirmed the assumption in Ref. [26] that the contribution of the O(k 4 ) term is not big than that of O(k 2 ) term at least within the t-channel cut k 2 = m 2 π . Admittedly, at NLO χPT level, the contributions from O(k 6 ) term or higher are not clear for us as well [26], and it definitely needs the knowledge of the higher order terms (NNLO, etc) from χPT, which is certainly beyond the scope of this work.
It should be worthwhile to stress that, for the isospin-2 ππ scattering, the m π r always has positive value [10]. However, for the isospin-0 case, it always holds negative value. To verify the valid of the Eq. (28), we straightforwardly use Eq. (22) to measure the corresponding results, where the Eqs. (24) and (25) are used, and the values of C 1 , C 2 , and C 3 in Eq. (37) are plugged in. The relevant results are displayed in top panel of Fig. 7. We also show the relevant results using Eq. (28), which are displayed in bottom panel of Fig. 7. One thing greatly comforting us is that both methods give similar results, which indicates that Eq. (28) is valid at least for the range of interest in this work. Analogously, we can discuss the validity of Eq. (29). In what follows, we can use Eq. (28) and Eq. (29). Similarly, to check the applicable scope of Eq. (30), we directly make use of Eq. (23) to measure the corresponding results. It is very inspiring that both methods also result in the somewhat similar results, which indicates that Eq. (30) is also valid at least for the range of interest in the present study. Admittedly, Eq. (30) is not simple style. For the large z-values (for example, z > 6), it deviates seriously from these using Eq. (23). In the rest of the analysis, we can use Eq. (30).
It is worth mentioning that the NLO expressions given in Eqs. (28), (29), and (30) are much more complicate than those of the isospin-2 [10]. Actually, it partially reflects the fact that the LO χPT can nicely reproduce the I = 2 s-wave ππ scattering length with just a 0.5% deviation as compared with the relevant experimental and theoretical results [1,8,9]. Nevertheless, in the isospin-0 channel, the agreement between LO χPT and corresponding experiments and theoretical results is departed by about 30% [1,8,9].

B. Chiral extrapolation of threshold parameters
In this work, lattice calculations are performed at pion masses: 247 MeV, 249 MeV and 314 MeV, respectively, according to the previous discussions, the influence of the σ meson can be reasonably ignored. In principle, we can exploit all of our data to carry out the relevant chiral extrapolation. Nevertheless, as it is pointed out in Ref. [26], the pion mass values should be small enough to make the NLO chiral expansion valid. For this purpose, we carry out the chiral extrapolation only using the data with two lower pion masses (247 MeV and 249 MeV).
Using NLO χPT expressions for the scattering length m π a in Eq. (24), the value of the constant C 1 can be obtained by fitting with lattice-measured m π a from Fit B provided in Table III. We can then translate the constant C 1 into the familiar ℓ I=0 ππ by Eq. (27). Moreover, the scattering length m π a at physical pion mass can be predicted using NLO χPT expressions denoted in Eq. (24). The chiral extrapolation of the scattering length m π a is shown in Fig. 8, and the extrapolated value at physical point is indicated by black circle on the physical line. The relevant results are given in Table IV as Fit-1.
On the same time, the corresponding fitting results with all the tree pion masses are also provided in Table IV as Fit-2, which agree with Fit-1 within statistical errors. The differences of two fits are considered as the estimated systematic uncertainties [26]. This leads to our ultimate results for scattering length m π a I=0 0 , ℓ I=0 ππ and C N LO where the superscript in the constant C 1 indicates that it is estimated at NLO χPT. Our lattice result of m π a I=0 0 is in reasonable agreement with the newer experimental and theoretical determinations as well as lattice calculations. In Table V, we compare this result, together with ℓ I=0 ππ , to these relevant results accessible in the literature.  With NLO χPT expressions for the effective range m π r in Eq. (28) and m 2 π ar in Eq. (29), the C 4 and C 5 values can be obtained by fitting with lattice-determined m π r and m 2 π ar from Fit B listed in Table III, where the C 1 is fixed to the cental value of C N LO 1 in Eq. (38), since we just have two lattice data at disposal. Moreover, the effective range m π r and m 2 π ar at physical pion mass can be predicted by NLO χPT. The relevant results are provided in Table IV as Fit-1.
According to Eq. (28), the effective range m π r is divergent in the chiral limit, this can partially explain the extrapolated value of m π r to the physical point usually has relative large statistical uncertainty, as compared with those of m π a I=0 0 and m 2 π ar. In practice, the statistical error of m π r is roughly estimated by the error of C 4 , and it is not dependent on the pion mass (or z-value).
The technique of the "moving" wall source introduced in Refs. [18,19] is exploited to calculate four diagrams classified for the I = 0 ππ scattering in Refs. [18,19]. Consequently, the signals of vacuum diagram are remarkably improved as compared with our previous studies [24,25], which enables us to not only measure the scattering length, but also explore the effective range. The chiral extrapolations of m π a I=0 0 and m π r are performed by NLO χPT. Extrapolated to the physical value of m π /f π , our final outcomes yield m π a I=0 0 = 0.217 (9)(5), m π r = −6.07 (44)(36), which are in reasonable agreement with the newer experimental and theoretical determinations as well as the lattice calculations.
In the interpolation of a fit to the lattice-measured values of k cot δ/m π during the region k 2 /m 2 π < 1.0, we include three leading parameters in the effective range expansion, which implement Liu et al.'s strategies [26] directly from lattice QCD. In particular, we confirmed that the effective range r and shape parameter P should be included for the successful fit [26].
For each lattice ensemble, we just calculate five points, simply due to the lack of computational resources, thus, robust extraction of shape parameter P definitely need more lattice data for each lattice ensemble. Admittedly, the most efficient way to improve the statistical errors of P is working on lattice ensembles with different size L for a given pion mass, as is done for the isospin-2 ππ scattering in Ref. [10].
So far, the influences of the higher order terms from χPT is quite limited, and we also cannot rule out that such contributions are significant [26]. It will be very interesting to systematically study the NNLO χPT expressions for leading three terms in the effective range expansion, and investigate ultimate numerical results due to its modifications from the NLO χPT expressions, we reserve this challenging work in the future.
The σ meson is clearly presented in low energy [47,72], and it is necessary to map out "avoided level crossings" between σ resonances and isospin-0 ππ states to secure the reliable scattering length as investigated in the πK scattering [47,72]. Luckily, according to Liu et al.'s discussions [26], the contaminations from σ meson for three lattice ensembles with small pion masses are negligible. Of course, sigma meson should be incorporated into more sophisticated lattice computation in the future.

ACKNOWLEDGMENTS
This work is partially supported by both the National Magnetic Confinement Fusion Program of China (Grant No. 2013GB109000) and the Fundamental Research Funds for the Central Universities (Grant No. 2010SCU23002). This project was partially supported by the High Performance Computing Center of Sichuan University, China. We would like to deliver our sincere gratitude to MILC Collaboration for permitting us to use the MILC gauge configurations and MILC codes. We deeply thank Carleton DeTar for his instructing Ziwen the necessary theoretical knowledge and computational skills for this work. We especially thank Michael Doering for his enlightening and constructive comments. We cordially express our boundless gratitude for the support of Professor Hou Qing and Professor He Yan. We would like to express our appreciation to Professor He Yan, Professor Huang Ling, Professor Wang Jun, and Professor Fu Zhe, and other warm-hearted people for donating us enough removable hard drives to store the propagators for the present work. We also express gratitude to the Institute of Nuclear Science and Technology, Sichuan University, from which the computer resources and electricity costs were furnished, and numerical calculations were partially carried out at both PowerLeader Clusters and AMAX, CENTOS, HP, and ThinkServer workstations. In this Appendix, we follow the original notations in Ref. [53] to transplant the FFT algorithm for the disconnected piece for sigma operator [53] into the vacuum diagram of ππ operator.
According to Eq. (8), the first part of the vacuum diagram can be expressed with quark propagators G as C V ππ (p, t 4 , t 3 , t 2 , t 1 ) = Re x2,x3 then Eq. (A1) can be rewritten as C V ππ (p, t 4 , t 3 , t 2 , t 1 )=Re x2,x3 e ip·(x2−x3) σ(x 2 , t 2 )σ(x 3 , t 3 ) , If we define the Fourier transform: then Eq. (A1) can be recast as where t = t 3 −t 1 , and we sum over all the time slice to improve the statistics. Note that t 2 = t 1 + 1 and t 4 = t 3 + 1. Of course, in the center-of-mass frame, the vacuum diagram is still accompanied by a vacuum subtraction [25].
For the elastic ππ scattering, the Mandelstam variables are written in units of physical pion mass squared m 2 π as where θ is the scattering angle, and k = |k| is the magnitude of the center-of-mass three-momentum of each pion. Note that the Legendre polynomials P 0 (cos θ) = 1 and P 2 (cos θ) = 3/2 cos 2 θ − 1.
To calculate the amplitudes on the isospin-0 ππ scattering, one expands the combinations with I = 0 in the s-channel as T 0 (s, t) = 3A(s, t, u) + A(t, u, s) + A(u, s, t).
For later simple notation, we denote Using the power representations of the loop integrals (B.1) in Ref [6], we can writē where we consider that the values of u and t are small due to the small k 2 value. Note that J(u) +J(t) = 1 96π 2 u + t + (u + t) 2 10 Now it is ready to show that Hence, it is obvious to show that It is obvious that we nicely reproduce the relevant results in Appendix C of Ref. [5] for the scattering length a I ℓ and the slope parameter b I ℓ at NLO, as expected. Note that we also give the explicit formula for the slope parameter c 0 0 in Eq. (B9) using the CGL language in Ref. [5]. In this appendix, to double-check our relevant results for the s-wave, we also list the results for d-wave, since the scattering amplitude in Eq. (B4) includes both information of them.