Lattice determination of $I= 0$ and 2 $\pi\pi$ scattering phase shifts with a physical pion mass

Phase shifts for $s$-wave $\pi\pi$ scattering in both the $I=0$ and $I=2$ channels are determined from a lattice QCD calculation performed on 741 gauge configurations obeying G-parity boundary conditions with a physical pion mass and lattice size of $32^3\times 64$. These results support our recent study of direct CP violation in $K\to\pi\pi$ decay \cite{Abbott:2020hxn}, improving our earlier 2015 calculation \cite{Bai:2015nea}. The phase shifts are determined for both stationary and moving $\pi\pi$ systems, at three ($I=0$) and four ($I=2$) different total momenta. We implement several $\pi\pi$ interpolating operators including a scalar bilinear"$\sigma$"operator and paired single-pion bilinear operators with the constituent pions carrying various relative momenta. Several techniques, including correlated fitting and a bootstrap determination of p-values have been used to refine the results and a comparison with the generalized eigenvalue problem (GEVP) method is given. A detailed systematic error analysis is performed which allows phase shift results to be presented at a fixed energy.


I. INTRODUCTION
The scattering of two pions is one of the simplest hadronic processes in QCD. Since the only meson involved, the pion, is the lightest hadron in the standard model and originates from the vacuum breaking of almost exact SU (2) L × SU (2) R chiral symmetry, the behavior of this process at low energy can be well described by chiral perturbation theory (ChPT) [3]. Although this scattering is not directly measurable by experiment, information can be inferred indirectly from K → ππeν e (Ke4) decay [4] and πN → ππN scattering [5].
Unfortunately from these experiments alone, it is difficult to obtain accurate ππ scattering phase shifts for a broad range of ππ energies, because, for example, of the limited energy available in Ke4 decays.
In addition to its intrinsic interest, ππ scattering is also an important ingredient in our recent calculation of the two-pion decay of the kaon [1] and the parameter ε , a highly sensitive measure of direct CP violation, which is a key component in the search for an explanation of the dominance of matter over antimatter in the present Universe. By comparing this lattice QCD result with the experimentally measured value we can gain a better understanding of CP violation in the standard model, with possible insights into the physics beyond it. While this result for ε is reported in Reference [1], essential components of this calculation are presented in two companion papers, an extensive study of G-parity boundary conditions [6] and the discussion of ππ scattering presented in this paper.
Because of the non-perturbative nature of QCD interactions, lattice QCD provides a unique, first-principles method with controlled systematic errors to determine the properties of low-energy QCD. With lattice QCD and the finite-volume Lüscher technique [7], we can calculate the ππ scattering phase shifts within the energy region from 2m π to approximately 4m π 1 while the interaction energies near the 2m π threshold can be used to determine the scattering lengths [10]. Such calculations complement the existing determinations of the scattering lengths obtained using chiral perturbation theory [11] and the dispersive calculations [12][13][14][15] of the energy dependence of the phase shift based on the Roy equations [16] and experimental input. For a recent review discussing experimental and theoretical results see also Ref. [17]. The dispersive technique might also be applied to extrapolate the lattice determination of the energy dependence of the phase shifts above the 4m π threshold. 1 For a recent example of an alternative method to extract scattering amplitudes from Euclidean correlators at energies above 4m π see Refs. [8] and [9].
Lattice QCD calculations of the ππ scattering phase shifts have been performed by a number of groups, but with pion masses heavier than the physical one so that a chiral extrapolation is required to obtain physical results, see e.g. Refs. [18][19][20][21][22][23][24][25][26]. (The exceptions to this statement are our previous physical pion mass calculations of the I = 2 and I = 0 s-wave phase shifts at a single energy close to the kaon mass, which enter the calculation of ∆I = 3/2 [27,28] and ∆I = 1/2 K → ππ decays [1,2] and a calculation of the I = 2 scattering length [29] at the physical pion mass.) This extrapolation is most likely valid for a scattering length calculation, but may become less trustworthy at higher energies where the accuracy of chiral perturbation theory becomes less certain. In this paper, we report the first lattice QCD calculation of both I = 0 and I = 2 s-wave ππ phase shifts performed over a range of two-pion energies with a physical pion mass so that a chiral extrapolation is no longer necessary.
As explained above, a central motivation for this study of ππ scattering is its importance for the calculation of the two-pion decay of the kaon [1] where it enters in three different ways. 1) The lattice calculation of the K → ππ decay matrix elements is performed in a finite volume while the matrix elements of interest are defined in infinite volume. The Lellouch-Lüscher (LL) factor which corrects for this difference is determined by the ππ interaction or, to be more specific, the derivative of the ππ scattering phase shift with respect to energy. 2) In order to determine the K → ππ decay matrix elements we need to know the amplitude with which the two-pion interpolating operators create the normalized finite-volume ππ states. The determination of these amplitudes is made difficult by excited state contamination as discussed in Section VII. 3). We need to know the finite-volume ππ state energies and the ground state energy should be close to but will not be exactly the same as the kaon mass. As is the case for the K → ππ calculation, this work is performed on a lattice with G-parity boundary conditions (GPBC) [6]. This choice is different from the periodic boundary conditions used in most lattice QCD calculations and we will discuss the advantages and drawbacks of this choice in Section II.
Our first calculation of I = 0 ππ scattering and K → ππ decay with physical kinematics [2] was published in 2015 and used the same 32 3 × 64 lattice volume, Möbius domain wall fermions and G-parity boundary conditions as the current calculation. In the earlier calculation we used a single I = 0 ππ interpolating operator to compute the scattering phase shift at an energy near the kaon mass using 216 configurations. The resulting I = 0 s-wave phase shift, δ 0 = 23.8(4.9)(1.2) • at center of mass energy E ππ = 498 (11) MeV, was significantly lower than the dispersive prediction of approximately 39 • at the kaon mass. (Here and later in this paper when two errors are given, the first is statistical and the second systematic.) Following our 2015 calculation and in light of the discrepancy between our results and the dispersive prediction, we devoted considerable effort to increasing our statistical precision.
We found that with 1400 configurations and the same four-quark ππ interpolating operator, a single-state fit to our data continued to be accurate but gave an even lower phase shift of 19.1(2.5) • , increasing the disagreement with the dispersive result [30]. As was the case with the original 216 configurations, performing a two-state fit to the two-point Green's function obtained from this single operator gave a ground state ππ energy and resulting ππ phase shift consistent with what was found from the single state fit. In addition to increasing the statistics we also experimented with adding a second, scalar bilinear ππ interpolating operator that we refer to as the σ operator and describe in more detail in Section III A.
Performing a two-state fit to the 2 × 2 matrix of two-point Green's functions obtained by including this operator revealed the presence of a previously unrecognized, nearby excited state, leading to a substantially smaller ground state energy and larger I = 0 phase shift [30].
Our present calculation builds upon this initial effort with increased statistics, additional ππ interpolating operators and a more accurate measure of the quality of the agreement between our data and our theoretical fitting formula. We also extend our calculation beyond a single ππ energy by computing ππ two-point correlation functions with the two pions carrying several values of the total momentum, allowing for an exploration of the scattering phase shifts at center-of-mass energies between approximately 2m π up to the kaon mass. Recognizing the importance of multiple operators we further increase the number of independent interpolating operators for both the stationary and moving frame calculations of the I = 0 and 2 phase shifts. Here we present results from 741 configurations using three operators in the moving frame calculation and three (for I = 0) and two (for I = 2) operators in the stationary frame calculation. With these additional operators we obtain a significant improvement in statistical precision. We are also better able to demonstrate control over the contamination from neglected excited states and to more reliably estimate their effects. We also applied a second approach to the analysis of our multi-operator, multi-state data, the generalized eigenvalue problem (GEVP) method. This new method gave results consistent with those of our traditional fitting approach with similar statistical errors.
The moving frame calculation allows us to directly calculate the LL factor from our lattice QCD data, the results of which are presented in Section VI D and utilized in Ref [1]. As 3) of Ref. [13]. This value was obtained in Ref. [13] using M π = 139.6 MeV. However, elsewhere in the present paper, we treat the neutral pion mass of M π = 135 MeV as the "physical" pion mass, following our previous papers on K → ππ decay and the conventions of the RBC and UKQCD collaborations. We therefore make a correction for the difference between the pion mass used in our lattice calculation, m π = 142.3(0.7) MeV, and this physical 135 MeV pion mass. We will discuss how we deal with the differences between these several different pion masses in greater detail in Sec.VI. Also note that here we have corrected our result to transfer the uncertainty in the energy at which the phase shift is determined onto the phase shift itself as described in Section VII and have correspondingly evaluated the dispersive prediction at this energy rather than at the kaon mass.
In the next section we describe the properties of the ensemble of gauge configurations that are used for both the calculation of ππ scattering presented here and our companion calculation of K → ππ decay and ε [2], together with a brief discussion about the G-parity boundary conditions adopted in these calculations. In Section III we present the operators we used, the matrix of two-point Green's function we measured and the statistical methods we used. In Sections IV and V we present in detail our fitting procedures and results for single pion and ππ energies and two-point function amplitudes, together with a brief comparison with another data analysis method, the generalized eigenvalue problem. With these results, in Section VI we describe how we obtain the phase shift results at various center-of-mass energies using a generalized form of Lüscher's formula. In Section VII we explain our new approach to determining the systematic uncertainties, estimate the largest systematic errors and present the resulting error budget. Finally in Section VIII we present our conclusions.
Its use enables us to work with a large, (4.6 fm) 3 spatial volume and therefore have good control over finite-volume systematic errors, without a dramatic increase in computational cost, albeit at the cost of increased discretization errors. We use G-parity boundary conditions (GPBC) in three spatial directions in order to obtain physical kinematics for the K → ππ decay.
The lattice parameters are equal to those of the 32ID ensemble documented in Refs. [31,36], with the addition of GPBC and a slightly lower pion mass of 142 MeV versus the 172 MeV used previously. This enables us to take advantage of existing results such as the value of the lattice spacing and also to compute the non-perturbative renormalization factors for the K → ππ matrix elements, in an environment free of the complexities associated with GPBC.
Below we first discuss the generation of these ensembles and measured properties including the plaquette, chiral condensate and autocorrelation times. We then discuss GPBC and how they differ from the usual periodic boundary condition (PBC). A detailed discussion on GPBC can be found in Ref. [6].

A. Ensemble Generation
The ensemble used for our 2015 calculation comprised 864 gauge configurations (after thermalization), with measurements performed on every fourth configuration giving 216 in total. Following our publication an error was discovered [37] in the generation of the random numbers used to set the conjugate momentum at the start of each Monte Carlo trajectory that introduced small correlations between widely separated lattice sites. While the effects were found to be two-to-three orders of magnitude smaller than our statistical errors, we nevertheless do not include these configurations in our new calculation.
In order to rapidly improve the statistical precision of our calculation we generated configurations via 7 independent Markov chains, each originating from widely separated configurations in our original thermalized ensemble. To compensate for any residual effects of the random number error we discarded the first 20 configurations of each stream, which is approximately 5 times the integrated autocorrelation time (see below). These configurations were generated using the hybrid Monte Carlo technique for which the Hamiltonian can be decomposed as (1) H = T + S G + S lQ (0.0001, 1) + S hQ (0.045, 1) + S DSDR , where T is the kinetic term and "lQ" and "hQ" denote the light and heavy quarks, respectively. The fermion actions comprise ratios of determinants Here the fractional power is required by the fact that for G-parity boundary conditions, the determinant of the squared matrix represents the contribution of four quark flavors [6], hence a square-root is required for the two light flavors and a fourth-root for the strange quark in order to perform a 2+1 flavor simulation. (Note that for periodic boundary conditions, the squared-matrix determinant represents the action of two quark flavors, hence for a 2+1 flavor simulation only the square root of the strangequark determinant is typically required.) The fractional power is achieved using the rational hybrid Monte Carlo (RHMC) algorithm. The light-quark action is further decomposed into two pieces using the Hasenbusch mass splitting technique [38] as follows: (3) S lQ (0.0001, 1) → S lQ (0.0001, 0.007) + S lQ (0.007, 1) .
We use an integration scheme comprising four levels of nested Omelyan integrator with Omelyan parameter λ = 0.22, using the layout detailed in Table I. Over 2200 configurations were generated in these 7 streams.  The top-most integrator (i = 1) is integrated over 16 steps of step-size 1/16 for a trajectory length of 1 MD time unit. Each update step of the action-integrator S i of size τ is divided into n S i equal-sized steps, with n S i given in the third column. For a detailed discussion of the nested integrator technique we refer the reader to Appendix A of Ref. [36] The use of RHMC for the light quark determinant introduces a significant cost overhead, primarily because the various mixed-precision techniques that have been developed to improve the efficiency of the standard CG algorithm are not generally applicable to the underlying multi-shift CG algorithm, which requires all the starting vectors to lie within the same Krylov space thus precluding the use of restarted methods. In addition we found that tighter stopping conditions on the inversion than are typical when applied to heavy quarks were required to ensure good acceptance and that more poles (20 in this case versus the O(10) required for a typical heavy quark calculation) were required to span the measured eigenvalue range. With some effort we were able to achieve a 70% performance increase, as measured on the IBM BlueGene/Q machines upon which a majority of our ensemble generation was performed, by combining a "reliable update" step with a subsequent loop over each pole with a conventional mixed-precision restarted CG [39].
A more significant improvement in the configuration generation was obtained by implementing the "exact one-flavor algorithm" (EOFA) [40,41] formulated by the TWQCD collaboration, in which a Hermitian positive-definite action for a single species of domain wall fermion is derived. The use of the EOFA allows us to circumvent the use of RHMC in the light quark sector, opening the door for various optimizations. We determined [42] that with a suitable preconditioning, the EOFA can be reformulated in a way that is not only more efficient but also allows for the re-use of the majority of our existing high-performance code for regular domain wall fermions. Coupled with algorithm and integrator tuning we achieved a 4.2× reduction in the time to generate a gauge configuration on the same hardware [42]. Utilizing this algorithm we extended 3 of our 7 streams by a total of nearly 3000 additional gauge configurations.
For this calculation we have measured on a subset of 741 configurations with consecutive measurements separated by 4 molecular dynamics time units (MDTU), which amounts to roughly 60% of the available configurations given this measurement separation.  density. The ensembles are given in the same order as in the plot shown in Figure 1. The chiral condensate and pseudoscalar density were measured only on the first three ensembles. The last three ensembles, separated by a horizontal line, were generated using the EOFA technique. The error bars were obtained using the jackknife technique applied to data that had been binned over blocks of size 12 MDTU.

B. Ensemble Properties
In Figure 1 we plot the evolution of the chiral condensate and pseudoscalar density as well as the plaquette. We remind the reader that the chains were each generated from already-thermalized and well-separated configurations of our original ensemble and hence we expect to observe no thermalization effects in these plots. The expectation values of these quantities for each of the ensembles are listed in Table II. In Figure 2 we plot the integrated autocorrelation time where C(∆) is the autocorrelation function with Here N is the number of samples v i of some quantity with meanv and standard deviation σ v .
The standard errors on τ int shown in the figure are estimated using a bootstrap resampling procedure applied to the quantities C(i, ∆) prior to performing the average in Eq. 5, and described in Refs. [31,43].
Because this is a non-standard application of resampling we provide here a detailed description and justification of the method. Since the functions C(i, ∆) from which the integrated autocorrelation time τ int (∆ cut ) is computed, depend on samples taken at different points in the Markov chain of gauge configurations, we do not attempt to obtain our error estimate by selecting samples from this Markov chain. Instead for each i we view the set of ∆ max values, { C(i, ∆)} 1≤∆≤∆max as a stochastic sample and deduce the statistical errors in our result for τ int (∆ cut ) from the fluctuations among these samples. This is analogous to the usual treatment of a two-operator correlation function M k (t) computed on a configuration k for a range of time separations t between the two operators. Here ∆ max is the largest value of ∆ that we include in our analysis. In this analysis we use ∆ max = 40. Thus, the number Consider the quantities C(i, ∆) and C(j, ∆) for j ≥ i, which contain measurements from the pairs of configurations (i, i+∆) and (j, j+∆), respectively. Due to the autocorrelations in the underlying data these values are correlated, with the corresponding correlation function peaking when j coincides with i or with i + ∆, and falling off exponentially in the time separation j − i away from these points. The secondary peak occurring when j is close to i + ∆ is smaller than the primary peak at j = i because only one of the two configurations involved in each pair coincide. Furthermore it is straightforward to show that the secondary peak vanishes entirely for ∆ sufficiently large that the configurations i and i + ∆ become effectively independent.
Thus, while the correlations between C(i, ∆) with different values of i have an unusual form, a standard binning procedure: is sufficient to generate values C(α, ∆) that, for a large enough bin size B, are statistically independent in α and that are therefore amenable to bootstrap (or jackknife) resampling.
The full procedure is then: 1. Truncate the collection of data to be analyzed to C(i, ∆) for 1 ≤ i ≤ N − ∆ max and 2. Bin this collection of data according to Eq. (7) producing N bin samples { C(α)} 1≤α≤N bin where each sample C(α) represents the ∆ max quantities { C(α, ∆)} 1≤∆≤∆max with each comprising N bin elements where b is the bootstrap ensemble index. For this analysis we use N boot = 500.

Compute the autocorrelation function under the bootstrap,
5. Compute the integrated autocorrelation function, The standard deviation of the bootstrap distribution of τ int,b (∆ cut ) over b provides an estimate of the standard error on τ int (∆ cut ). The appropriate bin size can be found, as usual, by increasing B until the error estimates stabilize; for the present analysis a bin size B = 15 was found to be sufficient. From Figure 2 we estimate an integrated autocorrelation time of τ int ∼3-4 MDTU, which is close to the separation of 4 MDTU between our measurements.
We therefore expect minimal autocorrelation effects on our measurements, but to be certain of our error estimates we will account for any residual effects using the non-overlapping block bootstrap procedure, as we will detail in Section III.

C. G-parity boundary conditions
The most significant difference between our calculation and those of other groups is the boundary conditions: in this work we use G-parity boundary conditions in all three spatial directions. G-parity is a symmetry of the QCD Lagrangian under charge conjugation coupled with a 180 degree isospin rotation about the y-axis. Applied as a spatial boundary condition on the quark fields this transforms a quark flavor doublet (u,d) into (−Cd T , Cū T ) as it passes through the boundary, where C is the 4 × 4 charge conjugation matrix.
As with all boundary condition variants the introduction of GPBC modifies the finitevolume spectrum. As described below, we take advantage of this change to improve the accuracy of our K → ππ calculation. With GPBC applied to the up and down quarks and after introducing a fictional doublet partner s to the strange quark to which GPBC are also applied (with the additional species suitably weighted out of the path integral, cf. Ref. [1]), a kaon state can be identified which satisfies periodic boundary condition while the pion states must satisfy anti-periodic boundary conditions (APBC). Recall that all three pions are odd under G-parity. This means that we can introduce a kaon ground state which is at rest, while the ππ ground state will be composed of two moving pions, with momenta close to ±π/L in each direction (the deviations from ±π/L being due to the ππ interactions we seek to measure). We can then tune the lattice parameters so that the initial kaon and the final ππ ground state have the same energy. This makes the K → ππ calculation much easier since we can focus on the dominant ground state contribution to the K → ππ matrix element. In contrast, on a lattice with PBC the ππ ground state will be composed with two nearly stationary pions, and we must tune the lattice spacing so that the energy of an excited ππ state matches the kaon mass. The matrix element of interest in this case will be a subdominant contribution to the three-point Green's function and obtaining a precise result becomes much more challenging. For more details on performing lattice simulations with G-parity boundary conditions including further discussion of the lattice symmetries and the treatment of the strange quark, we refer the reader to Ref. [6]. For the remainder of this subsection we will focus specifically on how these boundary conditions affect the measurement of the two-pion system.
Including GPBC introduces some significant differences from a calculation with PBC. Three significant differences might be identified. First, the ππ states that can be studied with these two types of boundary condition will be different. When non-interacting pions satisfy APBC in all three directions their allowed momenta become (2n 1 + 1, 2n 2 + 1, 2n 3 + 1) π L , where n i are integers. These are different from those on a volume with PBC, where the allowed momenta are (2n 1 , 2n 2 , 2n 3 ) π L . However, if we take advantage of moving frames, there is still a correspondence between the states that we can construct on a PBC volume and those present for a volume obeying GPBC. For example, if we wish to work with a ππ state comprising two pions at rest, for a volume with PBC we can do the calculation in the stationary frame, where the two component pion operators are constructed with zero momentum. However, for a GPBC volume the calculation can be performed in a moving frame where both pions have the same momentum, e.g. (π/L, π/L, π/L). With this choice, in the center-of-mass frame these two pions are at rest.
Since a moving frame calculation relies on a distorted volume which doesn't have cubic symmetry, there will be lower angular momentum partial waves whose phase shifts will enter the quantization condition that determines the s-wave phase shift, e.g. d-waves. In the stationary frame the lowest partial waves that enter beyond the s-wave are those with with l = 4. Fortunately in this work the interaction energies involved in our moving frame calculations are relatively small (around the kaon mass), and those higher partial waves that enter will have a negligible effect on the s-wave phase shift.
A second troublesome aspect of G-parity is the breaking of cubic symmetry at the quark level even for a lattice with cubic symmetry. As discussed in Ref. [6] there is a sign convention that can be chosen when G-parity is imposed in one direction that can be changed by changing the relative sign of the up and down quark fields. However, the choice of this sign in the remaining two directions is not conventional and breaks cubic symmetry by identifying one of the four diagonals connecting two corners of the cubic lattice and passing through its center. For a cubic volume in a stationary frame, the symmetry group is broken down from O h to D 3d [44].
There are two effects of this breaking of O h symmetry that we need to consider: First, its effect on the two-pion eigenstates of the QCD transfer matrix and second its effect on the rotational properties of the quark-level interpolating operators used to create those pions. Because of confinement the relevant degrees of freedom affected by the G-parity boundary conditions are the pions. Since G-parity boundary conditions are translationally invariant, for the O h -breaking properties of the quarks which make up the pion to have an effect, a single isolated quark must propagate across the lattice and through the boundary, a phenomenon that should be highly suppressed by effects of quark confinement. While this argument suggests that the two-pion eigenstates of the transfer matrix should fall into representations of the O h group, it is possible that the four-quark interpolating operators used to create these states will couple to more than one irreducible O h representation and care must be taken when constructing translationally covariant operators to suppress the creation of finite-volume states belonging to unwanted representations of the cubic symmetry group O h . This will be discussed when we write out the explicit form of these operators in Sec. III and the remaining cubic-symmetry breaking effects are discussed in Sec. VII.
Finally around-the-world effects in a moving frame will be different in a volume with GPBC compared to one with PBC. When we are performing a moving frame calculation in a volume with GPBC with one of the three smallest allowed total momenta (those with P tot = (±2, 0, 0), (±2, ±2, 0) or (±2, ±2, ±2) in units of π/L), the first-order around-theworld contribution will come from a single pion propagating from one ππ interpolating operator to the second (leg A) and a second single pion propagating from the second, through the time boundary to the first (leg B). This behavior is shown schematically as part of a later more detailed discussion in Figure 7.
For GPBC the momentum injected by each ππ interpolating operator can change the direction but not the magnitude of the momentum carried by the pion as it moves from leg A to leg B. Thus, for GPBC this around-the-world pion can carry the same energy on each leg and so that its contribution behaves as a constant when the time separation between the two operators is changed. We refer to this case where the pions in both legs carry momenta of minimum magnitude as the "first-order" around-the-world effect. The case in which the pion propagating in one of the legs carries momentum greater than the minimum is termed "second-order". Both cases are considered when performing the fits described in Section V B 1. In contrast, for the three smallest non-zero total momenta in a calculation with periodic boundary conditions all of the around-the-world terms will be time-dependent since the pions in the two legs will have different energies.

III. OVERVIEW OF THE MEASUREMENTS
In this section we describe the details of the interpolating operators used in this calculation, the two-point functions that we study and the specific contractions that are evaluated.
In the final subsection we outline the fitting methods employed and the methods used to determine a statistical error and assign a p-value to those fits.

A. Interpolating Operators
Here we discuss the structure of the interpolating operators used in this work. There are two different types of two-pion interpolating operators. The first type are denoted as "ππ(. . .)" operators and are constructed as the product of two single-pion interpolating operators and for which the parentheses and the quantity contained within are used both to specify the pion momenta and to distinguish these labels from the general set of ππ interpolating operators which can produce two pions when acting on the vacuum, the set in which all of our operators reside. The second type has the form of a quark-bilinear scalar sigma operator which shares the same quantum number as I = 0 ππ state. We start by constructing the single pion and sigma interpolating operators with momentum P = p + q, where p and q are the momenta of the individual quarks: where, using the notation of Ref. [6], Ψ and Ψ are the quark and anti-quark isospin doublets defined as: As explained in Ref. [6] the 2×2 flavor projection matrix (1+e inqπ σ 2 ) ensures that the quark In earlier studies this smearing was found to give a two-fold reduction in statistical errors [45].
With the operators constructed above, we use the all-to-all (A2A) propagator technique [46] to perform the measurements. The A2A technique divides the quark propagator into an exact low mode contribution which we can calculate using the Lanczos algorithm and a high mode contribution which can be accessed using stochastic approximation. In our calculation, we choose the number of low mode eigenvectors to be 900. For the high mode contribution, we perform spin, color, flavor and time dilution (i.e. we perform a separate inversion for each of the 24 colors, spins and flavors for each time slice). We use the same spatial field of random numbers for these 24 inversions but a different such field for each time slice [45]. We choose the number of random hits to be 1 (i.e. we use only a single random field on each time slice) since increasing it does not reduce the uncertainty [45].
We will work with two groups of pion operators. The first is labeled as π(111) with 8 different operators. These operators create pions carrying momenta (±π/L, ±π/L, ±π/L).
The second group is labeled as π(311) and contains 24 different operators. For this group one of the momentum components is replaced by ±3π/L.
We then combine two of these single-pion interpolating operators to construct ππ( p, q) operators with momenta P = p + q, where now p and q are the momenta of the individual pions: where α and β are isospin indices. As suggested by this equation, when we construct the ππ( p, q) operators, we separate the two single-pion operators in the time direction by 4 units.
This suppresses the statistical error from the disconnected diagrams (the V diagram below) by a factor of two in the I = 0 channel [47]. For consistency, when we construct the I = 2 ππ(. . .) operators, we also separate the two pion operators by 4 units in the time direction.

Momentum decomposition
The cubic symmetry breaking mentioned in Section II manifests as differences in the overlap factors between interpolating operators and finite-volume states whose momenta are related by cubic rotations (the energies themselves are not affected). In order to obtain ππ interpolating operators that respect the cubic symmetry and that can therefore be related to the continuum s-wave states, it is vital that we control this symmetry breaking. In Ref. [6] it was demonstrated that the cubic symmetry breaking in the pion states can be heavily suppressed by averaging over pairs of pion interpolating operators of the same total momenta but with different assignments of quark momenta. We apply this technique for the present work and extend it to include the sigma operator. The two quark and anti-quark momentum pairs for each pion momentum are listed in Appendix A. In Section VII we carefully analyze our data in order to account for any residual cubic symmetry breaking effects as a systematic error.
In evaluating the Wick contractions it is often convenient to utilize the γ 5 -hermiticity of the quark propagator G: where the dagger ( †) indicates the hermitian conjugate of the matrix in its spin, color and flavor indices, in order to exchange the source and sink for a particular quark propagator.
It is worth mentioning here that γ 5 -hermiticity is not an exact symmetry between the A2A approximations to the quark propagators used here because of the asymmetric treatment of the source and sink in the A2A approach. A further implication of our use of γ 5 -hermiticity to combine related contractions arises from the effective exchange of theq and q operators appearing in a meson field when γ 5 -hermiticity is used on both the propagator leavingq and that arriving at q. By symmetrizing over the assignments of momenta to theq and q factors in each meson field, we insure that this use of γ 5 -hermiticity does not result in a different amplitude. This determines the final pion interpolating operator we use: for each pion momentum we average over a total of four quark and anti-quark momentum assignments.
For the sigma operator, since it satisfies PBC and has zero momentum we average over the eight different quark momentum assignments that are listed in Appendix A to suppress cubic symmetry breaking. (Note: this symmetrical treatment of the quark and anti-quark components of the meson field implies that the contractions presented in Refs. [45] and [1] for the case of a local pion interpolating operator can be unambiguously extended to the case of a non-local meson field.)

Total momentum
We perform both a stationary-frame calculation where the total two-pion momentum is zero and moving-frame calculations for which the total momentum is non-zero. In the stationary-frame calculation we include the scalar σ operator for the I = 0 channel and for both isospin channels two classes of bilinear pair "ππ(. . .)" operators: One class has both pions in the group π(111) but with opposite momenta, which we label ππ(111, 111). The second class is made up of pions in the group π(311), again with opposite momenta and are labeled ππ(311, 311).
Total momentum Symmetry group Angular momentum Representation In the moving frame calculation we can also construct a ππ(111, 311) operator for which the constituent pion operators belong to the two different groups described above. For the present work we did not collect data using a sigma operator with non-zero momentum; however the analysis presented in the following sections suggests the inclusion of this operator may be beneficial in future work. In summary, we therefore have three different classes of operators in the moving-frame calculation for each isospin channel, as well as in the stationary frame I = 0 channel and only two classes of operators in the stationary-frame, I = 2 calculation. (Note: our notation distinguishing the two-pion interpolating operators does not specify the total momentum that they carry which must be determined from the context.)

Angular momentum
After identifying numerous ππ operators with different total momenta, the next step is to project those ππ operators onto angular momentum eigenstates. In this work we are interested in only the s-wave phase shift but we will also use d-wave states to estimate  Table III.
The second step is to construct an operator in the representation Γ by combining the operators in one of the classes described above using the characters of Γ. The detailed procedure is as follows: HereT [ p ] means we apply symmetry operationT on momentum p. We sum over all elementsT of the finite-volume symmetry group G, P is the total momentum, and χ(T ) is the character of each group elementT in the representation Γ. We choose p so that all the ππ operators appearing in the sum belong to the i th class. After projection, for each total momentum P , instead of having three or two classes of ππ operators, we will only have three or two ππ operators, each transforming under a specific representation of G and constructed from the operators within that class. Henceforth we will use the labels ππ(111, 111), ππ(111, 311), ππ(311, 311) to refer to those projected operators rather than to the classes from which they were constructed.
In the moving-frame calculations reported here, due to the limited number of classes (two) of single pion operators, we are only able to focus on the three sets of non-zero total momenta P with the smallest individual components: P = (±2, 0, 0) π L , (±2, ±2, 0) π L and (±2, ±2, ±2) π L , so that the number of different classes of ππ operators we construct on the lattice is more than one (three in this work).

B. Matrix of two-point correlation functions
We begin a discussion of the correlation functions using a single operator constrained to a single timeslice (recall our ππ(. . .) operators have the pion bilinears on separate timeslices).
For isospin I the two-point ππ correlation function is determined by the Euclidean Green's function where . . . indicates the expectation value from a Euclidean-space Feynman path integral, performed in a finite spatial volume of side L and time extent T , obeying periodic boundary conditions for the gauge field but anti-periodic boundary conditions for the fermion fields in the time direction and G-parity boundary conditions in the three spatial directions.
Here and in our two earlier papers [1,6] the hermitian conjugate which appears on the lefthand operator in Green's functions such as shown in Eq. (16) requires some explanation. For the case that the operator involves Euclidean fields evaluated at a single time, the hermitian conjugate represents a combination of path integral field variables which corresponds to the Hermitian conjugate of the indicated operator in the time-independent Schrödinger picture which is subsequently transformed to the time-dependent Heisenberg picture operator whose expectation values are described by a Euclidean path integral. For the case that the operator O I ππ is itself the product of two such operators evaluated at different times, each operator is to be interpreted in this fashion. In this case the two operators appearing in this pair are always symmetrized to insure that the resulting two-point functions are positive as this notation suggests in spite of the fact that their order is not exchanged by this prescription.
By inserting two complete sets of intermediate states, we can rewrite this two-point function as C I (t snk , t src ) = π|O I † ππ |π π|O I ππ |π e −tE π,in e −(T −t)Eπ,out in the limit where both t ≡ t snk −t src and T −t are large so that we can neglect the contribution from excited intermediate states. Notice the first term describes the "around-the-world effect", which is exponentially suppressed in T . Here E π,in and E π,out are the energies of the pions propagating from the source along the positive and negative time directions, respectively. These two energies should be the same in a stationary frame calculation but they may be different for a moving frame. The second and third terms, which can be combined scattering. This term is the largest source of statistical error because it is time-independent and therefore results in a decreasing signal-to-noise ratio as we increase the time separation t to suppress excited state contamination.
In practice, due to the rapid reduction in signal-to-noise ratio and the finite temporal extent of the lattice it is necessary to include data in the region where t or T − t is not very large. By including data from smaller time separations our results will be affected by contamination from excited-states. One way to suppress these errors is to expand the sum over intermediate states in Eq. (17) to include not only the ground state but also one or more excited states and then to fit using this more complicated expression. However, even if we only include one more state, performing such a multi-state fit may be difficult using a single interpolating operator since we are attempting to determine an increasing number of parameters purely from the time dependence of data with a rapidly falling single-to-noise.
While increasing statistics will ultimately allow the various states to be isolated, a far more powerful technique is to introduce additional interpolating operators which all share the same quantum numbers and therefore project onto the same set of states, albeit with different coefficients. While naively equivalent to increasing statistics, the additional operators actually introduce a wealth of new information that helps constrain the fit. This additional information can also be exploited more directly using the GEVP technique (described in more detail in Section V) whereby the energies of N states can be obtained from Green's functions comprising N operators using only three timeslices. A simpler method which allows for the detection of the presence of excited states using data from only a single timeslice by looking at the "normalized determinant" will be discussed in Section V.
In order to perform a stable fit where both ground and excited states are included, we introduce additional interpolating operators which all share the same quantum numbers so that the number of operators can be larger than or equal to the number of states included in the fit. Thus, we consider the matrix of two-point correlation functions: where indices i and j distinguish the operators. We can then expand Eq. (18) to include excited-state contributions: Now the excited state contamination error has been reduced since the lightest state that we neglect is the one with energy E m+1 , which is higher than E 1 , the energy of the first state that we neglected in Eq. (17). For simplicity in the discussion above we have identified a single time t snk/src that is associated with each two-pion operator. However, these operators are constructed from two, single-pion operators evaluated at the times t snk/src +4 and t snk/src as shown in Eq. (15). In the remainder of this paper we will use the variable t = t snk −t src −4 to describe the separation between the two operators which indicates a minimum distance of propagation needed to connect the two, two-pion operators.
Assuming that the fit is able to reliably obtain the parameters then clearly the larger number of states that are included in the fit, the smaller the resulting excited state contamination. However, given the added computational cost and resulting fit complexity, we should be careful to include only operators which help to distinguish the relevant excited states. An important criterion, discussed later, is the degree to which the operators introduced overlap with the state being studied or a common set of excited states.

C. Contraction diagrams
We are interested in the scattering process for specific isospin channels. The I = 0 and I = 2 ππ state with I z = 2 can be constructed from π + , π − , π 0 states as below: The matrix of two-point correlation functions for the ππ and σ operators can be obtained from a linear combination of eight different diagrams, labeled as C, D, R, V , C σππ , V σππ , C σσ and V σσ , each corresponding to a particular Wick contraction that is identified in Fig. 3.
Their definition in terms of quark propagator is given in Appendix C. They can be combined to obtain the two-point correlation functions as follows: If we were to perform the contractions for each of the different total momenta by substituting Eq. (15) into Eqs. (22) and (23), the number of different contractions to be evaluated for each gauge configuration would be 7848, which is unnecessarily large. The technique which we employ to reduce the number of momentum combinations takes advantage of three kinds of symmetry in ππ scattering: parity symmetry, which corresponds to changing each momentum from p to − p , axis permutation symmetry, which permutes the three coordinate axes and an "auxiliary-diagram" symmetry, which relies on the combination of γ 5 hermiticity and the "around-the-world" contraction to show that two diagrams whose source and sink momenta satisfy a special relation are identical (For more details about the "auxiliary-diagram" symmetry, we refer readers to Ref. [48]). Using a subset of the gauge configurations in this study, we have found that excluding all but one of the momentum combinations that are related by these three symmetries does not increase the statistical error for the measured ππ energy. This strategy substantially reduces the number of momentum combinations from 7848 to 1037 [48].

D. Estimating statistical errors and goodness of fit
In this paper we use multi-state correlated fits to determine the energies of each state and the overlap amplitudes between the different states and operators. The fitting procedure is flexible, e.g. we can perform a fit where the number of operators and states are different and we can perform a "frozen fit" where some of the parameters are held fixed during the fit, which is useful in the excited-state error analysis. An important benefit of our fitting procedure is our ability to calculate a p-value, which is a measure of how well our data matches with our theoretical expectation for the time dependence of the two-point function being analyzed.
However, the determination of statistical errors and the calculation of a p-value are not straightforward. Not only are we performing a correlated fit where the covariance matrix is itself determined by the data and therefore has its own, often substantial uncertainties, but there are autocorrelations between configurations, since the sampling interval between neighboring configurations used in our analysis is comparable to or smaller than the autocorrelation time which separates truly independent samples. While our number of samples, 741, is relatively large compared to many lattice calculations, if we group these samples into bins of two or four and thereby reduce the autocorrelations between these binned samples, the resulting decrease in the effective number of samples loses significant information about the fluctuations which is required for adequate control of the covariance matrix upon which our correlated fits are based.
Fortunately, we have developed methods to solve both of these issues. These methods are based on a combination of the jackknife and the non-overlapping blocked-bootstrap resampling techniques [49]. The bootstrap technique uses uncorrelated, non-overlapping blocks of data for its samples and gives statistical errors unaffected by the autocorrelation between our 741 samples. However, the inner jackknife resampling introduced to calculate the covariance matrix for each outer bootstrap sample is applied to the unbinned data obtained as a union of all of the blocks in a given jackknife sample. In this paper the block size is chosen to be 8 to suppress the effects of autocorrelation. Finally the distribution of bootstrap means about the mean for the entire sample, determines the proper χ 2 distribution that can be used to correctly determine the p-value for the fit. (Recall that the usual standard χ 2 distribution is not accurate when χ 2 is determined using an uncertain covariance matrix in the presence of autocorrelations.) More details of this method can be found in Ref. [49].

IV. SINGLE PION ENERGIES AND MASS
In order to determine the pion energy and mass, we calculate a two-point function using the neutral pion operator: for all possible values of t src and t snk and then we average over t src while keeping t = t snk −t src fixed. We have in total 32 different pion momenta, 8 from the π(111) group of operators and the other 24 from the π(311) group. Up to the effects of the cubic symmetry breaking induced by the boundary conditions, which are heavily suppressed by the procedure discussed in Section III and the residual effects shown to be negligible in Section VII, the two point functions within each group are related by cubic rotations hence we average the two-point functions within each group. This leaves us with two correlation functions, C i π (t), where i ∈ {(111), (311)} represents the momentum of the pion without specifying its direction.
We then perform correlated fits of each correlation function to the form using various fit ranges, all of which share the same upper limit t max = 29. Here A i π is related to the normalization of the operator O i π while E i π is the energy of a moving pion state with momentum (1, 1, 1) π L or (3, 1, 1) π L . The fitted results for E i π plotted as a function of t min are shown in Figure 4. From both plots we can see a clear plateau starting from t min = 14. The result for E π is insensitive to t max so we make the same choice of t max = 29 as was made in Ref. [2]. For those reasons we choose the fit range to be 14 − 29 and the fit results for that choice are listed in Table IV. The good p-values for both fits suggest that our data is well described by this single-state model.
Knowledge of the mass of the pion is required for the determination of the ππ phase shifts via the Lüscher procedure. Unfortunately, with GPBC we are unable to measure this mass directly and must instead infer it from the energy of a moving state with a suitable choice of dispersion relation. In Table IV we give the results of applying the continuum dispersion relation to the (111) and (311) moving pion energies, which are labeled as m π,CD .
We can see that the resulting masses are inconsistent, which we interpret as the result of discretization effects on the dispersion relation. We also calculate the pion mass using the dispersion relation obeyed by a free particle on our discrete lattice where the pion mass is identified as the energy of a pion with zero-momentum. The results are listed in Table IV as m π,LD and are consistent between the two momenta.
The large discrepancy between the two pion masses calculated using different dispersion relations suggests that when we calculate the pion mass using the larger-momenta π(311) operators the result has not only a statistical error that is 3 times larger than that from the π(111) operators, but also a large systematic error. For the remainder of this paper, we will use m π,CD = 142.3(0.7) MeV calculated from the π(111) operators using the continuum dispersion relation as the pion mass. This 142.3(0.7) MeV value differs from the physical pion mass of 135 MeV by 7 MeV. This introduces an "unphysical pion mass" error into our results which will be discussed in Sections VI and VII. We will neglect the discretization error that remains in our determination of the pion mass since the 1 MeV discrepancy between the m π,CD and m π,LD in Table IV is small compared to the 7 MeV "unphysical pion mass" error identified above. cases. Here E π is shown in lattice units with t max fixed to be 29.
We have converted to units of MeV by using the inverse lattice spacing for this ensemble, 1/a = 1.3784(68) GeV, where the error on a also has been propagated into the errors on the energies given here.

V. FINITE-VOLUME ππ ENERGIES
In this section we describe our multi-state, multi-operator fitting strategies and the resulting fit parameters for both the stationary frame and the moving frame calculations and for both the I = 0 and I = 2 channels. Since these four situations are different, we will discuss them separately. At the end of this section we briefly discuss results obtained from another data analysis technique, the GEVP. This both provides alternative results for these quantities and an opportunity to compare these two methods. Because the primary focus of this paper is on the properties of the ππ ground state, this discussion of the GEVP method is limited to the ground state energies which it determines.
A. Stationary frame In the stationary I = 2 channel, we have two classes of operators, ππ(111, 111) and ππ(311, 311). We project them onto the trivial A 1 representation of the cubic symmetry group, which is the approximate symmetry group of a finite-volume lattice. (A discussion of possible cubic symmetry breaking effects resulting from our G-parity boundary conditions will be presented in Section VII.) This projection results in two different ππ operators, . We then calculate the matrix of two-point functions constructed from these two operators by measuring where ∆ = 4 is the time-separation between two pion fields used to construct each ππ operator. We average over all values of t src while fixing t and then average the data at t with that at t = T − t − 2∆ to improve the statistics. (The individual single-pion operators at the times t snk/src + ∆ and t snk/src that make up each two-pion operator are constructed to be identical so when taking this second average we are combining equivalent physical quantities.) We then try two different fitting strategies: 1) Fit the single two-point function C aa (t) assuming a single intermediate state and an around-the-world constant using the form where A describes the normalization of the operator, E ππ is the energy of the finite-volume ππ ground state and B is the around-the-world constant. Thus, a total of three fit parameters are required. We neglect all data related to the second operator O b so this is a one-operator, one-state fit.
2) Fit the upper triangular component of the 2 × 2 matrix of two-point functions C ab using two intermediate states and three different around-the-world constants using the form where A ix is the overlap between the i th operator and the x th state; E x is the energy of the  For each case, we perform correlated fits with various choices for t min and set t max = 25.
We plot the resulting ground state energy as a function of t min in the left panel of Figure 5.
As we can see from the plot, the introduction of the second operator does not noticeably improve the fit result, as the ground state energies given by both fitting strategies are statistically consistent for all t min and the statistical errors are also consistent. As we increase t min , the ground state energy first decreases, which suggests a non-negligible excited state contamination for small t min and then reaches a plateau for t min ≈ 10. We adopt the 2operator, 2-state fit with the fitting range of 10 − 25 for our final result. In Table V we list the p-values and the final parameters obtained from that approach. We observe an excellent p-value indicating a strong consistency between the data and our model.
The fact that B aa is 60σ resolved from zero suggests the importance of including these around-the-world constants in our fits. This conclusion can also be reached by performing a similar fit in which the only change is that these constants are excluded. These fits give p-values that are consistent with zero, suggesting that these constants are required.
We also observe that the matrix of overlap amplitudes A ix is nearly diagonal, where the where the second term represents the vacuum subtraction which removes the disconnected piece in Eq. (17), since it does not contribute to ππ scattering. We then average over all t src while fixing t = t snk − t src and average the data at t with that at We then explore three different fitting strategies: 1) Fit C aa (t) using a single state and the equation where A and E have the same physical meaning as in the stationary I = 2 fit. This is a one-operator, one-state fit and we have only two fit parameters in total. In contrast with   in the case of total momentum (0, 0, 0) that is discussed in this subsection. The next three columns from the right show the parameters obtained from three-operator, three-state fits for three non-zero values of the total momentum. the stationary I = 2 fit, here we neglect the around-the-world constant since an estimate of the size of the dominant contribution resulting from a single pion propagating through the temporal boundary gives a value which is approximately ten times smaller than the statistical error on these noisier I = 0 channel data. Note, if fit as a free parameter, the result for this around-the-world constant is consistent with zero and gives a ground-state energy consistent with the result obtained when this constant is excluded, but with a statistical error that is 50% larger.
2) Fit the upper triangular components of the 2 × 2 submatrix spanned by O a and one of the other two operators using two states and the equation where N = 2, A ix is the overlap amplitude between the i th operator and the x th state, E x is the energy of the x th finite-volume state and (i, j) takes values from either {a, b} or {a, c}.
Thus, this is a six-parameter fit. An analysis similar to that mentioned in 1) above shows that the three around-the-world constants should be excluded.
3) Fit the upper triangular component of the entire 3 × 3 matrix of two-point functions using two or three states and the fitting form given in Eq. (32) where N = 2 or 3 is the number of states we include in the fit. We neglect the around-the-world constants for the same reasons as above, resulting in 12 (N=3) or 8 (N=2) fit parameters in total.
For each fitting strategy, we perform correlated fits with various values of t min and set t max = 15. We do not extend t max to 25 as we did for the I = 2 channel since the data for t > 15 have larger statistical errors than in the I = 2 case, so including them will not benefit our fit. However, adding more fit points will destabilize the correlation matrix inversion procedure because of its increased dimension. We also risk introducing data for which the neglected around-the-world contribution may be a dominant component of the large-time data that has been introduced. This behavior is suggested because although the around-the-world constants remain statistically consistent with zero the p-value does fall as t max is increased. Note that we do not observe any corresponding statistically significant effects on the amplitudes and energies as t max is increased suggesting that our fits remain robust even in the presence of around-the-world contributions. A similar issue is encountered for the moving frame I = 0 fits and is discussed in greater detail in Section V B 2.
We plot the ground state energies from these fits as a function of t min in the right panel of Fig. 5. For the three-operator case, we perform the three-state fit for t min ≤ 4 while for t min ≥ 5 we use the two-state fit as we observed that the three-state fits with t min ≥ 5 were unstable and did not converge for many bootstrap samples, indicating that the third state can no longer be reliably resolved in the data. As we increase the number of operators, the ground state energy at fixed t min becomes significantly lower and the plateau region becomes more clear and begins earlier. We conclude that in contrast with the I = 2 channel, the introduction of the two extra operators, especially the σ interpolating operator, substantially reduces not only the statistical error but also the systematic error resulting from excited state contamination.
Since the plateau region for the three-operator fit starts at t = 6, we choose the threeoperator, two-state fit with a fitting range of 6 − 15 to determine our final results. In the right-hand column of Table VI we list the p-value and final parameters for that fit. We can see that especially for the ππ(111, 111) (a) and σ (c) the overlap amplitudes between a given operator and the two states are of comparable size, which explains the effectiveness of the multiple operators that we included. This large overlap factors between operators and states is consistent with the fact that the phase shift and hence the ππ interaction strength, is considerably larger than in the I = 2 case. Hence the exchange of momentum between the two pions required for the mixing between states is enhanced. For I = 2 the two ππ operators assign momenta with different magnitudes to the pions and would each couple to a different ππ energy eigenstate if the pions were non-interacting.
The fact that the overlap between operator O a and the first excited state is about half of the overlap of that operator with the ground state also provides a strong indication that there is likely to be non-negligible excited-state contamination in a single-operator, singlestate fit. This explains the substantial discrepancy between the phase shift at an energy near the kaon mass that we published in Ref. [2] and both the results presented here and those from the earlier dispersive prediction [13]. This can also be seen in the right panel of an energy that is consistent with our previously published value but which is substantially larger than the ground-state revealed by the introduction of the additional operators.  and ππ(311, 311). We project them onto the trivial representation of the little group of the cubic symmetry group which leaves the total momentum unchanged. These little groups   The upper, middle and lower panels are for total momenta (2, 2, 2) π L , (2, 2, 0) π L and (2, 0, 0) π L , respectively. Our final results were obtained from three-operator, three-state fits. Reading from top to bottom the values of t min for these results were for I = 2 t min = 10, 12, 11 and for I = 0 t min = 6, 8 and 7.
For each choice of P tot , this gives us three different operators, O a = ππ A 1 (111, 111), (111, 311) and O c = ππ A 1 (311, 311). We calculate the matrix of two-point functions constructed from these three operators, C ij (t snk , t src ) and combine the various values of t src and t snk in the same way as was done for the stationary I = 2 calculation, except for an extra step where for each value of | P tot |, we also average over all of the possible total momentum directions. This leaves us with three correlation matrices, one for each | P tot |. We then try three different fitting strategies for each | P tot |: 1) Fit C aa alone with a single state and an around-the-world constant, as we did in the stationary I = 2 calculation.
2) Fit the upper triangular component of the 2 × 2 submatrix spanned by O a and one of the other two operators using two states and three different around-the-world constants using the equation For each value of | P tot | and fitting strategy, we perform correlated fits with t max = 25, vary the value of t min and plot the fitted ground state energy as a function of t min in Figure 6, as in the stationary I = 2 calculation.
Similar to the stationary I = 2 calculation, for all three values of | P tot |, the introduction of the two extra operators has little impact on the ground state energy. As we increase t min , the ground state energy first decreases, suggesting a non-negligible excited state error for small t min and then reaches the plateau region. This plateau starts at t min = 11 for P tot = (±2, 0, 0) π L , t min = 12 for P tot = (±2, ±2, 0) π L and t min = 10 for P tot = (±2, ±2, ±2) π L . We choose the three-operator, three-state fit with tmax=25 and tmin fixed to the start of the plateau region identified above.
In Table V we list the p-value and the final parameters for each choice of P tot . With the chosen fit ranges we observe excellent p-values for all values of the total momentum. The fact that for each P tot , B aa is 100σ resolved from zero suggests the importance of including these around-the-world constants in the fitting. The overlap matrices are all nearly diagonal as in the stationary I = 2 calculation so that each operator is dominated by a different one of the three states. Thus, as was the case for the stationary frame I = 2 calculation, this explains why the introduction of these two additional operators does not improve the determination of the ground state energy.
It is also worth mentioning that the constant terms we include in the fit only describe the lowest-order around-the-world (ATW) effect mentioned in Sec. II C, where both the pions on leg A (direct propagation between the two single-pion operators) and leg B (propagation through the temporal boundary) carry a minimum momenta with components ±π/L. Here we refer to segments of an around-the-world propagation path identified in Fig. 7. In contrast to the stationary case, the higher-order ATW terms in the moving frame need not be described by a constant term in the Green's function. For example, one of the pions on leg A or leg B could be replaced by a pion one of whose components has the larger ±3π/L value.
This possibility still conserves momentum and will show an exponential time dependence.
When compared with the first-order ATW effect, this second-order ATW effect is exponentially suppressed by the energy difference between a pion with three ±π/L momentum components and a pion with one component increased to ±3π/L. However, in our calculation, due to the time separation ∆ between the two single-pion operators that make up our ππ operator, this second-order effect can be enhanced in some cases. For example, we can look at the Green's function constructed from two O b operators. We define the state that propagates between the two temporally-separated pion operators in our ππ operator as the "internal state". Notice we have two internal states here, since we have two ππ operators.
For the first-order case, the two internal states cannot both be the vacuum while conserving momentum, but for the second-order effect they can. This is illustrated in Figure 7. Thus, in this example the second-order effect is enhanced at least by a factor of e E∆ = 4.2, where E is the lowest energy of the internal state which we approximate by the I = 0 ππ energy.
In order to investigate the size of the higher-order ATW terms we perform a fit to the I = 2 data. It can be easily shown that third-and higher-order ATW effects are always exponentially suppressed when compared with the first-order and second-order effects. This means we can perform a fit which includes some extra parameters which represent the second-order ATW effect and neglect third-and higher-order effects. Here we fit the matrix of correlation functions with the following fit function: Compared with Eq. (33), the extra term with coefficient D ij describes the second order ATW effect. Here E π 0 and E π 1 are the energies of moving pions with momenta (1, 1, 1)π/L and (3, 1, 1)π/L, respectively. Their values can be obtained from Table IV. The fitting results for the ground state ππ energy and the sample-by-sample difference between the results with and without the second order ATW effect are shown in Table VII. Since the difference is negligible and statistically consistent with zero, we conclude that we need not include the second-or higher-order ATW effects in our fits.  strategies that are similar to the three used for the moving I = 2 case, except that we exclude the around-the-world constants from all the fits. The effect of these constants will be discussed below. We then perform correlated fits with t max = 15, vary t min and plot the ground state energy as a function of t min in Figure 6. Another interesting feature is seen in the errors of the fitted parameters when we perform a single-operator, single-state fit using only the ππ(111, 111) operator. Consider how the sizes of either the relative error of the amplitude, or the absolute error of the ground state energy change as we decrease the total momentum from (2, 2, 2) 2π L to (0, 0, 0) π L , when the fit range is fixed (e.g., 6 − 15). The pattern is that these errors increase as the total momentum decreases, as can be seen in Table VIII! This behavior conflicts with the expectation that these errors would be approximately the same based on the Lepage argument [50]. For our kinematics, the non-zero total momentum is created by reversing some of the momentum Notice that according to Eq. (15), the ππ operator with definite momentum P will contain 1, 2, 4 and 8 terms for the four cases above with total momentum containing three, two, one or zero non-zero components, respectively. Since the number of two-point function contractions that we must evaluate grows like the product of the numbers of terms in its constituent operators, there will be 1, 4, 16 and 64 different contractions needed for each type with total momentum (2, 2, 2) π L , (2, 2, 0) π L , (2, 0, 0) π L and (0, 0, 0) π L , respectively. All these terms contribute V-type diagrams and hence to the error of the Green's function, with each of approximately the same size.
However, not all of these terms contribute D-type diagrams (and hence to the central value of the Green's function), because of a mismatch between the momenta. This can be understood by looking at the non-interacting limit in which a D-type diagram will only be non-zero when the two independent single-pion Green's functions are non-zero. This happens only when each of the two pions in the source ππ operator has the opposite momentum to that carried by one of the pions in the sink ππ operator, which we call momentum matching.
Counting these terms gives the numbers of momentum-matched D-type diagrams for these four possible total momenta: 2, 4, 8 and 16, respectively. After dropping a common factor of two in this counting of D-type diagrams, these estimates suggest that the proportions between the relative errors for these four total momenta become  Table VI changes as we decrease the total momentum. A similar analysis can be applied to the I = 2 channel, which suggests that these relative errors should be approximately the same, independent of the total momentum, which is consistent with what is shown in Table V.
Similar to the moving I = 2 channel, as we increase t min , the ground state energy first decreases, which suggests a non-negligible excited-state contamination for small t min . The ground state energy then reaches a plateau region which starts with t min = 7 for P tot = (±2, 0, 0) π L , t min = 8 for P tot = (±2, ±2, 0) π L and t min = 6 for P tot = (±2, ±2, ±2) π L . For our final result we choose the three-operator, three-state fit with t min equal to the beginning of the plateau region identified above. The choice of t max is more subtle and will be discussed together with the effect of the neglected around-the-world constants below. In Table VI we list for each P tot the fit range, fit procedure, as well as the resulting p-value and final parameters. We also observe in this table a trend towards smaller overlap factor between the operators and states, i.e. a more diagonal amplitude matrix A ix , as we increase the total momentum and thus decrease the center-of-mass energy. This is again consistent with our understanding of the relation between this overlap factor and the strength of the ππ interaction, which also decreases as the center-of-mass energy is decreased. This also explains why the additional operators appear to have the largest impact on the ground-state energy for the moving frame, P tot = (2, 0, 0)π/L fits in Figure 6.
Next we discuss our treatment of the around-the-world constants in the fit. There are two potential sources of systematic error in our results that must be treated carefully: the excited state contamination and the around-the-world contributions. The first error is expected to be much more significant and is discussed in Section VII. To leading exponential order in the time extent of the lattice volume, the around-the-world contributions are time-independent constants even in this moving frame calculation because of our G-parity boundary conditions. We observe that fitting with these around-the-world constants as free parameters results in good p-values for t min ≥ 6 but gives results for the constants that are either statistically consistent with 0 (P tot = (2, 2, 2) and (2, 0, 0)), or which have an unphysical, negative sign (P tot = (2, 2, 0)). For all three cases either their errors when the constants are unresolved or the non-zero fitted values when these constants can be resolved are ten times larger than the expected size, that of the I=2 around-the-world constant.
For the case of P tot = (2, 2, 2) or (2, 0, 0), we can neglect these constants in the fit since there is no statistical inconsistency between the fitted energy with and without these constants and we expect that the effects of the true around-the-world constants will be approximately ten times smaller than these sub-statistical effects. Note that excluding theses constants from the fit gives us an improvement in the statistical error of the ground state energy by a factor of 1.2−1.5. For the second P tot = (2, 2, 0) case, the most likely explanation is that the constants are acting as "nuisance parameters" that help to partially account for the excited state contamination but do not reflect true around-the-world behavior. Rather than leaving the constants as free parameters and using an unphysical model to describe our data we choose to fix the constants to zero and to account for the systematic, excited-state contamination errors separately.
The model with zero around-the-world contributions should be a good description of the data in the window [t min , t max ] for which t min is large enough that excited state effects are small and t max small enough that the contribution of these constants is small relative to the size of the data. For t max = 15 we observe very poor p-values even for large t min ≤ 10. Reducing t max from 15 to 10 we observe a significant improvement in the goodnessof-fit, finding acceptable p-values for t min ≥ 6. This behavior is consistent with the effects of around-the-world contributions, although the excited state contributions may also play a role. Note however that, despite the dramatic improvement in p-value observed when reducing t max from 15 to 10, we observe consistency in the ground-state fit parameters and no loss of precision, suggesting that the around-the-world systematic error is negligible and that the fits are under good control. In Figure 6 we use t max = 15 to show that reasonable behavior is seen when the around-the-world constants are omitted even for this large value of t max . Further evidence that supports the argument that for P tot = (2, 2, 0) these constants are "nuisance parameters" can be found by including them in the fit, fixing t max and increasing t min . The resulting around-the-world constants monotonically decrease with increasing t min , which suggests that they likely result from excited state contamination, which is expected to decrease as t min is increased, rather than representing the effects of single-pion aroundthe-world propagation. For uniformity, in Table VI we choose to list the results for the three smallest total momenta with the same value of t max = 15. This results in the small p-value of 0.016 for the (2, 2, 0) case. However, had we used t max = 12 we would have obtained equivalent results with a p-value of 0.205.

C. Normalized determinant
It is important to emphasize that the introduction of these additional ππ(. . .) operators (in all cases) and σ operators (in the stationary I = 0 case) offers something more than a simple statistical improvement but gives new information about the underlying energy eigenstates.
The two-point Green's functions C xy for x = a and/or y = a typically have larger statistical uncertainties than C aa at the same t, suggesting that including these additional operators may lead to only a small reduction in the statistical errors of the fitting parameters. However, in some cases including these operators significantly improves the statistical error of the ground state energy, (e.g. the stationary I = 0 case shown in Figure 5). We also observe in several cases a significant reduction in the energy of the apparent plateau as well as an earlier onset of the plateau region, suggesting that the extra operators are dramatically improving our ability to resolve nearby excited states which may be very difficult to distinguish from the ground state when we have only a single operator, even with large statistics.
Some insight into how this improvement comes about can be gained by considering the "normalized determinant" of the N × N matrix of Green's functions, N (t), defined as where C(t) is the matrix of Green's functions. We normalize the determinant using the product of the diagonal elements of the matrix so that this quantity does not depend on the scale of the interpolating operators. In fact, it can be shown that 0 ≤ |N (t)|≤ 1.
If the number of intermediate states that contribute to C(t), N , is smaller than N then N (t) = 0 (since the N , N -component vectors constructed from the matrix elements of the N operators between these N states and the vacuum, which determine the N × N matrix C ij do not span the entire N dimensional space on which C ij acts). Thus, if at a given time t we find N (t) = 0, then we can be certain that at least N distinct states are contributing to C(t). However, if we find N (t) = 0 we cannot conclude that there are fewer than N states in the Hilbert space that contribute to the correlation matrix C(t). Thus, we cannot use N (t) to tell us if a sufficient number of operators has been used to distinguish all of the states that contribute to C(t). When N (t) ≈ 1, it suggests that these operators create states from the vacuum which are orthogonal to each other.
In Figure 8 we plot N as a function of t for both the 2 × 2 matrix of stationary I = 0 two-point functions comprising ππ(111, 111) and σ operators and the 3×3 matrix of Green's functions constructed from all three operators. For the two-dimensional matrix case we find at t = t min = 6, N (t) = 0.31 (7) giving unambiguous proof that more than one state must be present, while for the three-dimensional matrix case, N (t) is relatively suppressed and takes value of 0.14(13) at t = 5 and consistent with zero at t ≥ 6. The observation that the third state can no longer be distinguished from the noise for t ≥ 5 explains why we were unable to perform reliable 3-operator, 3-state fits to the I = 0 stationary two-point functions with t min ≥ 5 earlier in this section. This is closely related to the discussion of the size of the excited state systematic error, as will be explained in Section VII. We emphasize that the determinant is computed from two-point function measurements at a single time separation E 0 = 547.5(6), E 1 = 725(5) and E 2 = 948 (13). While the differences between these energies are sufficient to easily see the time dependence of shown in Figure 8, they are insufficient to be resolved in a single-operator fit, even with the statistical precision achieved with 741 configurations.

D. Comparison of multi-operator multi-state fits with the GEVP method
Multi-parameter fitting is a straightforward method to analyze the correlation functions between pairs of interpolating operators to determine the energies of finite-volume states which these operators create and the overlap amplitudes between these operators and states.
A second approach to analyze such data is the generalized eigenvalue problem (GEVP) approach [51,52]. The GEVP can be viewed as a generalization of the concept of effective mass, from single-operator to multiple-operator Green's functions. In principle, this approach has good control over the systematic error resulting from the excited states that are not included in the analysis. Following the notation of Ref. [52], the N -dimension GEVP can be defined as where C(t) is the N-dimensional matrix of two-point functions, v n , 1 ≤ n ≤ N are the eigenvectors and λ n (t, t 0 ) are the corresponding generalized eigenvalues. (In this section only, we follow the conventions of Ref. [52], and construct the correlation function C ij (t) from the product O i (t)O j (0) † .) In the limit where the lattice temporal extent, T , is large, the energy of the n th state is related to λ n by E eff n (t, t 0 ) = log(λ n (t, t 0 )) − log(λ n (t + 1, t 0 )) .
The GEVP approach can also be used to construct an operator A † n which creates the normalized lattice energy eigenstate: It has been shown that in the region where t 0 > t/2, the systematic error in the energy of the i th state resulting from states omitted from the analysis is constrained by [52] If T is not sufficiently large, we need to consider two complications to the GEVP procedures described above. The first is around-the-world propagation, which introduces timeindependent constants into the correlation functions for both isospin channels for each of our four values of total momentum. One way to eliminate this effect is to introduce a "subtracted matrix of two-point functions" D(t) defined as D(t) = C(t) − C(t + δt) and use this D matrix in the GEVP calculation [20]. Notice that this step will not affect the formula for the energy, but a modification is needed for the operator A n which is now given by The second complication comes from backward propagating states. One way to accommodate this effect is to modify the relation between the eigenvalue and the corresponding energy. For more detail, see Ref. [53]. We will not use this method here but instead work with a smaller time range where the effect of backward propagating states can be neglected.
In this paper we will compare only the ground state energies obtained in our fitting and GEVP analyses. Specifically we compare the energies obtained from our multi-operator simultaneous fits to the GEVP effective energy defined in Eq. (37). Our comparison of the fitting and GEVP approaches begins by comparing the fitting results with those obtained from the GEVP at a fixed time t chosen to be the same as the value of t min used in the fit while the value of the GEVP quantity t 0 is chosen as t 2 . We find that the GEVP energy is statistically consistent with the fit result but its statistical error is about five times larger. Actually the result of this direct comparison should not be surprising, since much information is lost when the GEVP method is applied to a single (t, t 0 ) pair.
An improvement to the GEVP method proposed in Ref. [20] is to fit the set of generalized eigenvalues λ 0 (t, t 0 ), as a function of t with t 0 fixed and to include in that fit possible correction terms from omitted excited states. This addresses the statistical noise problem identified above by including more of the correlation function data in the GEVP analysis.
We adopt a simple version of this fitting approach and perform a correlated fit to the GEVP eigenvalues of the form fitting all of the data from t = t 0 + 1 to a largest value t max . As shown in Ref. [52], this functional form for λ 0 (t, t 0 ) will contain errors from neglected higher energy states which are bounded by By choosing the smallest value of t used in this GEVP fit to be one time unit above t 0 we are minimizing the statistical error in our result for E 0 for a given choice of t 0 . We then treat t 0 in the same spirit as t min in our previous multi-parameter fitting of the matrix of Green's functions (which we will call the "usual fit" in the remainder of this section). Thus, we vary t 0 searching for a plateau region for sufficient large t 0 and then adopt as the result of this GEVP fitting that value obtained for E 0 from the smallest value of t 0 within that plateau region.
While this procedure is similar to what is used in our usual fit, we have not carried out the detailed discussion of the residual systematic errors coming from excited state contamination that is attempted in Section VII for our usual fit. Thus, we are unable to say if choosing a larger value of t than t 0 + 1 would have resulted in a sufficiently reduced systematic error to give a reduction in the total error, overcoming the increase in statistical error that would result from increasing t − t 0 above one.
In addition to examining the dependence of the GEVP result for E 0 on t 0 , we must also make sure that our choice of t max is appropriate: if t max is too large, then neglecting the backward propagating state will introduce an error; if t max is too small, then we will have a small number of input data points which makes our fit less reliable. We will choose t max to be no larger than that used in our usual fit, so that we can use those earlier results to determine how we should treat the around-the-world effects. Thus, based on the results obtained from our usual fit, we neglect the around the world effects for the I = 0 channel and work directly with the correlation matrix C(t) of two-point functions. For I = 2 these effects were found to be important so in that case we analyze the subtracted matrix D(t) defined above, using δt = 1. For the I = 2 channel, where we perform this matrix subtraction process, t max is taken to be 20, smaller than the t max = 25 in the usual fit, due to the increased noise resulting from the construction of the subtracted matrix of two-point functions D(t). The value of t max in the I = 0 channel is chosen to be 15, the largest used in the usual fit. We can then plot the ground state energy from this GEVP fitting as a function of t 0 and look for the beginning of the plateau region and also the p-value in order to determine t 0 .
Plots that include both the GEVP fit and the usual fit results are shown in Fig. 9. Notice the x-axis represents t min for the usual fit, and t 0 for the GEVP method. This choice for plotting was made to best align both the central values and statistical errors from the two methods. For all but the stationary I = 0 case, both the usual fit and the GEVP fit are P tot,I sub t 0 fitting range E 0 (GEVP) E 0 (usual fit) p-value     (I = 0). The total momenta from the top down are (2, 2, 2) π L , (2, 2, 0) π L , (2, 0, 0) π L and 0. Here the x-axis represents t min for the usual fit, and t 0 for the GEVP fit. In the legend of the lower right panel, 2D and 3D indicate a 2 × 2 and 3 × 3 GEVP matrix.
can see that for the I = 0 channel the GEVP fit results are approximately 2σ larger than the results from the usual fit if we only include the statistical error. However, the two are consistent if the systematic errors arising from excited state contamination in the usual fit are included. These excited-state error estimates were obtained by independent methods as described in Section VII and are also listed in Table XVI and match surprisingly well the differences between the results obtained from our usual and GEVP fits.
For the I = 2 channel the two results are consistent with each other with comparable statistical errors which are much smaller than in the I = 0 case. Also notice that in all cases, the GEVP method gives statistical errors that are no larger than the usual fit method, which suggests that the GEVP fitting method is a useful tool for analyzing the matrix of Green's functions for the scattering considered here. We also expect that the usual fits may be less successful than the GEVP method when we increase the number of operators in the fit due to instabilities that will likely result from the larger number of fit parameters. As used here, the GEVP fit is far simpler, being a correlated fit to a single one-parameter function of the time. The multi-operator fits will further suffer from the quadratic increase in the number of elements in the covariance matrix, the inversion of which may become unstable once it becomes too large. In the GEVP case, the size of the covariance matrix depends only on the size of the fit window.
Nevertheless, for the problem at hand, fitting the matrix of two-point functions is more direct than fitting the GEVP eigenvalues and, as we will show in Section VII, allows considerable flexibility in estimating the size of the systematic error arising from omitted excited states. In the traditional GEVP method [52] the number of operators must be larger than or equal to the number of states, while in the multi-parameter fitting approach more states than operators can be easily accommodated. Currently, this is a crucial step in estimating the excited state error as will be seen in Section VII and in identifying a plateau in the stationary I = 0 case as can be seen from the lower right panel of Fig. 9. However, for the case studied here we have not attempted to make similar estimates of the systematic errors in our GEVP results and effective methods for estimating such errors may well be possible there.

VI. DETERMINATION OF THE PHASE SHIFT
In this section we discuss in detail how we determine the ππ phase shifts from the finitevolume ππ energies. We begin with Lüscher's formula generalized to the case of anti-periodic boundary conditions 2 with a general total momentum. Next we discuss the strategy of working with energy differences to reduce discretization errors, especially for the moving frame calculations. We then calculate the ππ scattering phase shifts at various center-ofmass energies for both isospin channels using this technique. We also describe our method for specifying the energy at which these phase shifts have been determined in order to reduce the effects of the slightly unphysical pion mass used in our lattice calculation. Finally we calculate the Lellouch-Lüscher factor [54] that is needed to interpret the finite-volume K → ππ calculation [1].

A. Lüscher's quantization condition for non-zero total momentum and anti-periodic boundary conditions
Euclidean-space lattice QCD calculations determine finite-volume ππ energies from which the infinite-volume scattering phase shifts can be obtained using an approach developed by Lüscher [7]. While initially derived for the case of a stationary frame and periodic boundary conditions, this approach was later generalized to non-zero total momentum [55][56][57] and later this moving frame result was generalized to the case with general twisted boundary conditions [58]. In this section, we will write down this result for our GPBC lattice where pions satisfy the anti-periodic boundary conditions. In particular, the s-wave phase shift, δ(s) can be determined from the relation: δ(s) + φ d,γ (s) = nπ where n is an integer and typically allows for more than one solution to this quantization condition, resulting in a series of energy eigenstates in a single volume. The function φ d,γ (s) is defined by Here s is the square of the invariant mass of the two-pion system, γ is the Lorentz factor which boosts the laboratory frame to the CM frame and q is related to the magnitude 2 In this section we will focus on the case of anti-periodic boundary conditions obeyed by the pions, which result from the G-parity boundary conditions obeyed by the quarks.
of the momentum k carried by either pion in their center-of-mass frame. Each of these quantities can be determined from the finite-volume ππ energy, E ππ obtained from the lattice calculation: The vector of integers, d, is related to the total momentum P tot by and Z d,γ 00 (ŝ, q 2 ) is the generalized Lüscher's zeta function, which is defined as where the conventional argument s of the zeta function is replaced here byŝ to remove the possible confusion with the square of the center-of-mass energy and the set N d,γ is defined as Here the vector represents the effect of the boundary conditions. If the particle satisfies periodic boundary conditions in the i th direction then i = 0, while with anti-periodic boundary conditions we have i = 1. The quantityγ −1 is a linear transformation on 3vectors defined asγ using the notation of Ref. [59].
As defined in Eq. (46) the generalized zeta function diverges atŝ = 1 and needs to be expressed differently to be evaluated atŝ = 1. We use a simple generalization of a formula given in Ref. [59] to do this. Combining all of these formulae we can obtain the ππ scattering phase shift at the energy √ s from the finite-volume energy eigenvalue E ππ determined from our lattice calculation. Note that in obtaining Eq (43) we are implicitly neglecting the contributions to the scattering of partial waves with l ≥ 1. In the stationary frame, assuming the ππ operators are constructed in the trivial representation of the cubic group, cubic symmetry prevents states with 1 ≤ l < 4 from contributing [7]. The interaction strength in the l ≥ 4 channels is known to be small and these interactions can be safely neglected [7].
However in the moving frame the relativistic length contraction naturally breaks the cubic symmetry down to a smaller group, the trivial representation of which also allows for contributions from d-wave (l = 2) interactions. Previous calculations [30] have shown that the phase shifts in the l = 2 channel are small around the kaon mass, which will be further suppressed in the moving frame calculation where √ s is smaller than m K so we can therefore continue to assume s-wave dominance. As described above, the G-parity boundary conditions also break the cubic symmetry but the effects can be suppressed with a careful choice of ππ operator. Any systematic errors arising from this source are discussed further in Section VII A.

B. Calculation technique
As shown above, the ππ scattering phase shift is related to the energy of a finite-volume ππ state, or more specifically to the "pion momentum" k carried by either pion outside the range of the strong force. However, on a discretized lattice with anti-periodic boundary conditions (i.e. a case where the single-pion ground-state has non-zero momentum), the determination of k from the measured ππ energy must be performed carefully. If the ππ interaction is relatively weak then k, which is a measure of that interaction, will be close to its free field value and we must take precautions that the potentially small difference between k and its free field value is nevertheless large when compared with the discretization errors associated with the spatial momenta of the pions in our calculation. However, as can be seen in Eq. (44), k is determined from differences of larger quantities and care must be taken to insure that the quantities being subtracted have, to the degree possible, common finite lattice spacing errors so that these errors will largely cancel in the difference. Specifically the quantities being subtracted should be chosen so that their difference will vanish in the limit that the ππ interactions vanish, even when computed at finite lattice spacing.
This cancellation of finite lattice spacing errors can be accomplished by working with two related quantities determined from our calculation: ∆E = E ππ − 2E π , which measures the ππ interaction strength, and E π , the lowest energy of a moving pion. Using ∆E for example, the effects of the finite lattice spacing upon the pion dispersion relation that enter both E ππ and 2E π will largely cancel, leaving only the subtler effects of the discretization upon the two-pion interaction itself. Even for the case of non-zero total momentum P tot , we will exploit our choice of anti-periodic boundary conditions in all three directions and use for E π the ground-state, single-pion energy. Each of the three non-zero total momenta that we study can be formed from two pions carrying the minimum allowed momenta p = (±1, ±1, ±1)π/L so that 2E π will be the minimum energy of two interacting pions in the limit in which that interaction vanishes.
Thus, the quantities ∆E and E π will be computed on the lattice, and systematic errors estimated to account for the residual effects of the finite lattice spacing. The results are finite-volume predictions for ∆E and E π in the continuum limit, albeit with an unphysical pion mass, and in Section VII we will estimate and propagate the systematic discretization errors on these quantities. (While it would be better to determine ∆E and E π by performing calculations at multiple lattice spacings and taking the continuum limit, this is at present beyond our available resources.) Adopting this strategy to account for the discretization effects, we can then apply the generalization of Lüscher's finite-volume quantization condition without ambiguity using the continuum dispersion relation where we have continued to assume anti-periodic boundary conditions in three directions and note that the last two terms on the right-hand side of Eq. (49) are exactly known.
The second line in Eq. (49) demonstrates the purpose of this rearrangement: our result for k is proportional to the small quantity ∆E (a measure of the ππ interaction strength) plus other kinematic quantities that are determined without finite lattice spacing error. This guarantees that at finite lattice spacing the phase shift determined in this way from the quantization condition will vanish when ∆E → 0 so that the fractional finite lattice spacing errors expected in ∆E can be directly propagated to determine the corresponding error in δ(s).
We thereby obtain a value for the phase shift at an unphysical pion mass, from which a prediction for the physical phase shift can be obtained by assigning suitable systematic errors for discretization effects and the unphysical pion mass as will be discussed in Section VII.

C. Phase shift results (statistical error only)
In this section we tabulate our results for the ππ scattering phase shifts including their statistical errors, computed according to the method described above. For each result we must specify the energy at which the phase shift takes the quoted value and we choose to assign the appropriate √ s in such a way as to minimize the error introduced by the unphysical pion mass, 143 MeV, at which our calculation is performed, a value 6% larger than the 135 MeV which we adopt as the physical pion mass in this paper. For values of √ s on the order of the kaon mass, the error associated with this unphysical 8 MeV shift in the pion mass is small and is estimated using chiral perturbation theory in Section VII.
However, since the I = 0 and 2 phase shifts vanish when √ s = 2m π , this unphysical pion mass error can become large as √ s approaches 2m π,unphy > 2m π,phy . This effect can be easily eliminated if we view our computed phase shifts as functions of the pion momentum k rather than √ s, since a calculation of the phase shifts will give results which vanish at k = 0 independent of the pion mass.
Thus, for each computed value of the phase shift we use the measured lattice ππ energy and lattice pion mass to obtain the relative momentum k, and then when presenting our results for the phase shifts assign an energy determined by combining this momentum in the continuum limit with the physical pion mass by applying the dispersion relation s = 4(k 2 + m 2 π,phy ) .
The effect of the unphysical pion mass on the actual strength of the interaction (i.e. upon the phase shift itself) is small and is treated as a systematic error that we estimate using chiral perturbation theory, as discussed in Sec. VII.
In Table X we list the phase shifts calculated from 4(3) different momenta of the center of mass for the I = 2(0) channel and calculate the corresponding √ s using Eq. (50). Here we only include the statistical error. The full error budget will be discussed in Section VII.
We do not provide a result for the phase shift of the I = 0 channel in the case where P tot = (2, 2, 2) π L . At this lowest value of √ s the attractive interaction between two pions results in a center-of-mass ππ energy that lies below 2m π .
Instead, we calculate the scattering length for both isospin channels with the moving frame where P tot = (2, 2, 2) π L . For the I = 2 channel, we start with the following expansion of the phase shift as a function of the relative momentum k: where a 0 is the scattering length and r eff is the effective range. Since k is very small where P tot = (2, 2, 2) π L , we can keep only the leading term: The result is listed in Table X. For the I = 0 channel, we generalize Eq. (1.3) in Ref. [10] to a moving frame, and keep only the leading term while neglecting all terms with higher power of a 0 /L.
This approximation is reasonable since the contribution from these higher-order terms is of O(10) smaller than the statistical error and the systematic error discussed later. The result is also shown in Table X. We point out that the scattering length for I = 0 channel given in that Table is inconsistent with the experimental value at 3σ level if only the statistical error is considered. This inconsistency mostly comes from the excited state systematic error, as will be discussed in Sec. VII.

D. Lellouch-Lüscher factor
In our companion calculation of the I = 0 K → ππ matrix elements [1], an important ingredient is the Lellouch-Lüscher factor [54], which removes both the difference in normalization between states defined in finite and infinite volume and the leading power-law finite-volume corrections to the finite-volume matrix element. This factor is defined as: where δ I is the isospin I, s-wave ππ phase shift an φ d,γ is defined in Eq. (43). This formula should be evaluated at E ππ = m K .
The moving frame calculation enables us to determine the phase shifts at various energies, which allows us to perform an ab initio measurement of ∂δ 0 ∂k using a finite-difference approximation. We now focus on the I = 0 case since our calculation has been tuned to give s close to m 2 K for the I = 0, ππ ground state. We approximate the factor F using  the I = 2(0) channel and the corresponding √ s, together with the I = 2(0) scattering length calculated from moving frame calculation with total momentum (2, 2, 2) π L . Here the statistical error of each phase shift is obtained not by simply propagating the statistical error of E ππ , but a more elaborate method discussed in Sec. VII G which removes the uncertainty of the energy at which we quote the phase shift. two different methods. In the first we subtract the values of δ 0 at P tot = (0, 0, 0) π L and P tot = (2, 0, 0) π L , which gives ∂δ 0 ∂k = 0.372(153) (54) and in the second method we replace the second total momentum with P tot = (2, 2, 0) π L , and obtain ∂δ 0 ∂k = 0.205 (63) .
Both results are consistent with which is calculated from the dispersive analysis [1,13]. Note we have not attempted to account for systematic effects arising from the finite-difference approximation or other effects here. Nevertheless we find good agreement between our lattice results and the dispersive prediction, albeit with large statistical errors. These results are also presented in Ref. [1] where the dispersive result was used for the final analysis. Note that these values differ slightly (within errors) due to different choices of fit range and the finite-difference approxi-mation being applied there to the phase shift is a function of energy rather than a function of k.

VII. SYSTEMATIC ERROR ANALYSIS
There are several sources of systematic error which affect our results: the breaking of cubic symmetry by our G-parity boundary conditions, the non-zero lattice spacing of our single gauge ensemble, the unphysical value of our pion mass and contamination of our multi-operator, multi-state fits due to the presence of additional excited states, not included in our fit. In this section, we describe our procedure for estimating the size of these errors.
The full error budget for the phase shifts we obtain is given at the end of this section and a comparison is made with the dispersive predictions [13].

A. Cubic symmetry breaking
One distinguishing feature of our calculation is our choice of boundary conditions: we use G-parity instead of the standard periodic boundary conditions commonly used in other ππ scattering calculations. As discussed in Sec. II and in Ref. [6], G-parity boundary terms in the quark action break the usual cubic symmetry of our lattice action and cubic volume.
We will distinguish two possible effects of this breaking of cubic symmetry by the boundary conditions: the effects on the finite-volume eigenstates of the transfer matrix and the limitations on the symmetry properties of interpolating operators constructed from the quark fields.
Since the physical states in our finite volume are pions which obey cubically anti-periodic boundary conditions, we expect that the effects of this quark-level cubic asymmetry will be suppressed exponentially in the linear size of our spatial volume. Local phenomena will not be affected by these boundary terms but only phenomena which span the entire volume.
This consideration should apply to the size of the corrections to the standard ππ finitevolume quantization condition, reducing these G-parity cubic symmetry breaking effects to the size of other finite-volume corrections.
Of greater concern is our inability to confidently use cubic symmetry when interpreting the rotational quantum numbers of the states produced by our interpolating operators. The G-parity breaking of cubic symmetry limits the selection of quark momenta that can be introduced when constructing interpolating operators resulting in operators which contain a mixture of representations of the cubic group. The only solution to this problem which we have found is an empirical one: we must carefully construct pion interpolating operators to reduce the mixing of different cubic symmetry representations below the level that we are able to observe.
As described in Ref. [6], a numerical investigation on single-pion correlation functions has been performed on a smaller lattice, which suggests that if we construct these pion interpolating operators with a single choice of quark momentum assignment chosen from the set of allowed quark momenta (for example choice 1 of Appendix A of the present paper), then we observe a clear cubic symmetry breaking effect in the overall normalization of the corresponding two-point functions. In that appendix we also introduce a second choice with the same total momentum but with different assignments of quark momentum. We observe that if we construct our pion interpolating operators by averaging the two momentum choices, then the resulting cubic symmetry breaking becomes sub-statistical. Since this cubic symmetry breaking is purely due to the boundary condition, it will be further suppressed by the larger volume used in the current study, and is therefore negligible in this work. While the normalization of the two single pion operators carrying momenta which are related by cubic symmetry show small differences, the pion energies are always the same providing evidence for the assertion in the preceding paragraph that the spectrum of the transfer matrix shows only exponentially small cubic asymmetry.
We can also calculate the size of the cubic symmetry breaking in our ππ interpolating operators directly by studying the overlap between interpolating operators belonging to different representations of the cubic group. We will focus on the stationary frame since in the moving frame calculation the s and d-waves are coupled to each other even if we have exact cubic symmetry. If we have exact cubic symmetry, we can project all three groups of ππ interpolating operators (they are ππ(111, 111), ππ(311, 311) and σ) onto the A 1 and T 2 representations, which primarily map onto the l = 0 and l = 2 representations of the continuum rotation group, respectively, and for each group of operators these two representations will be orthogonal. However, if the symmetry group is reduced to the D 3d group, the T 2 representation of O h group will no longer be irreducible and will contain the can then serve as a measure of the cubic symmetry breaking.
We start by considering the two projections of the ππ(111, 111) operators and define their overlap as the average where ∆ = 4, as introduced earlier in Eq. (27). Here and below we follow a convention similar to that introduced in Section V in which the labels a, b and c correspond to ππ(111, 111), ππ(311, 311) and σ respectively. If the cubic symmetry breaking effects are negligible, then (t) will be consistent with 0.
Since we are interested in the size of these cubic symmetry breaking effects relative to the correlation functions from which we obtain our results, we will present the normalized where C I,T 2 a,a (t) and C I,A 1 a,a (t) are defined by: where R = A 1 or T 2 . The ratio R I,T 2 ,A The results are shown in Fig. 10, where the left panel shows the normalized overlap amplitude for the I = 0 channel, and the right panel shows that for the I = 2 channel. Here and in the later graphs shown in Figure 11 we choose the time ranges to best present our results. We exclude large times because the statistical errors become very large and would require a highly compressed scale to display. However, in each case sufficiently large times are shown that the signal from the states which we study should be an important contributor to the correlation function, so the small size of R I,T 2 ,A 1 aa (t) for those later times implies at most a fractional percent contamination of our results from cubic symmetry breaking. Of The results are most easily interpreted if we again examine the normalized ratio The results are shown in Fig. 11, where the upper panel shows the overlap between appendix, we can achieve accurate cubic symmetry at the meson level, despite the symmetry breaking at the quark level. We therefore do not assign any systematic error arising from cubic symmetry breaking.

B. Finite lattice spacing
Fundamental to the connection between the scattering phase shifts and the two-particle finite-volume energies is the recognition that it is the interaction between the particles, described by a non-zero scattering phase shift, that causes the two-particle energy in finite volume to be shifted away from the simple spectrum of non-interacting particles in a box.
When adopting formulae to determine the scattering phase shifts from the two-particle finite-volume energies in Section VI we were careful to preserve this connection for non-zero lattice spacing.
Specifically Eq. (49) determines the relative center of mass momentum k between the two pions in the finite-volume ππ ground state that enters Lüscher's quantization condition as a function of the energy difference ∆E between the finite volume ππ energy and that of two non-interacting pions. This difference would vanish in the absence of interactions, even at non-zero lattice spacing. We then assign a relative systematic error to this measured energy difference that is of the same size as is found for other similar quantities computed on this ensemble for which a continuum limit has been evaluated. Thus, we use where c is chosen from the finite lattice spacing errors reported in Ref. [31]. In detail we use the average ChPTFV value of the magnitudes of c ID f , c ID f (K) , c ID w 0 ,a and c ID √ t 0 ,a , given in Table  XVII in that paper which are the four coefficients which describe the a 2 finite lattice spacing errors for these four different physical quantities computed with the same lattice action and gauge coupling as used here. This gives a relative error of 1.6% for ∆E which we round up to 2%.
Note this relative error is usually a small quantity, and therefore a small absolute error when compared with E π and E ππ since we are calculating the ππ scattering phase shift at relatively low energies (near the kaon mass). The error determined from Eq. (62) is then propagated in the standard way to obtain the O(a 2 ) error for the scattering phase shift which we use for each of our four values of total momentum.

C. Finite volume
Finite volume affects the energy of ππ states in two ways. The first effect results in the quantized finite-volume energies, is described by the Lüscher quantization condition and can be viewed for large L as a power law effect. The second effect falls exponentially with the system size and is caused by the interaction radius being a finite fraction of the system size or, equivalently, the effect of off-shell singularities when the Poisson summation formula is used to estimate finite-volume effects. This second effect is usually much smaller than the first and is the source of the systematic error considered here. This exponentially suppressed correction for the I = 2 channel for periodic boundary conditions can be formulated as [60]: where for the case of near-zero relative momentum. According to Fig. 2 from Ref [60], this correction introduces an approximate 1% relative error in the scattering length for a volume with periodic boundary conditions but the same size and physical parameters as the volume with G-parity boundary conditions studied here.
For our G-parity boundary condition lattice, since the pion satisfies anti-periodic boundary conditions, the formula for ∆ F V has to be modified as follows: This leads to a relative error of approximately 0.6%. We round this number up to 1% and adopt it as an estimate of the finite volume effects for our more general case which includes the I = 0 channel, non-zero ππ relative momentum and non-zero total momentum.

D. Unphysical kinematics
The pion mass which we measured on this ensemble is 142.3 MeV, which is 5% larger than our choice for the physical pion mass (135 MeV) at which we wish to determine the scattering phase shifts. We deal with this pion-mass mismatch in two steps. In the first step we shift the ππ energy at which we quote the phase shift, as has been discussed in Sec. VI, by expressing the phase shift as a function of the two-pion relative momentum k in the center-of-mass system and then identifying this value of k with a ππ energy using the physical pion mass. We view this correction, which will be large for energies near the ππ threshold, as the most important effect of this pion-mass mismatch. In the second step we account for the remaining effects of this pion mass mismatch as a systematic error in our result for the phase shift. We estimate the remaining pion-mass-mismatch error by using ChPT to calculate the difference between the scattering phase shift evaluated at these two different pion masses but at the same value for k. The NLO ChPT prediction for the scattering amplitude of both the I=0 and I=2 channels for small relative momenta are listed in Appendix B, and the predicted phase shift difference as a function of √ s is plotted in Fig. 12.
There is a remaining uncertainty in this approach that must be resolved. The ChPT calculation is only valid for small k, a condition not valid for our stationary calculation, which results in the rapid rise of the phase shift difference in the I = 0 channel when √ s > 380 MeV. Similar behavior is seen also in the I = 2 channel, although the breakdown appears to occur more slowly as a function of √ s, suggesting the ChPT result for this channel can be considered sufficiently reliable for energies in our range of interest. We modify our systematic error determination for √ s ≈ m K to resolve this issue. Notice that the dispersive prediction, whose range of validity is expected to extend above that of ChPT, shows a relation between the phase shift and √ s that is close to linear for a broad range of ππ energy, up to and including the kaon mass [13]. This suggests that the relation between the  According to Eq. 105 from Ref. [55], under the condition where the D-wave phase shift is small, the correction to the s-wave phase shift after including the D-wave phase shift is: where and Here we generalize the numerical recipes for evaluating Lüscher's zeta function given in Ref. [59] to the general cases where l = 0 to evaluate Z d,γ 20 (1, q 2 ) in Eq. (68). The D-wave phase shift, δ 2 (p), can be evaluated by assuming that the energy in the center of mass frame is small in our moving frame calculation so that we can approximate the D-wave phase shift by using the D-wave scattering length: where a I 2 is the D-wave scattering length for isospin-I. The values we used are listed in Tab. XII, which were obtained from Ref. [61], Table VI  In Sec. V we tried different fit ranges to find a balance between minimizing the excited state contamination error and the statistical error. For the preferred fit ranges shown in Tables V and VI, there is no obvious evidence that the neglected excited states give a significant contribution to the fitting result, as can be seen in the energy plots in Fig. 5 and Fig. 6. As discussed in Ref. [30], the presence of an apparent plateau in the fitted energy as However, while such extra assumptions may have been inappropriate for the determination of the preferred central value, their introduction may be a reasonable approach to estimate an excited-state systematic error.
The difficulty associated with an extra state fit can be illustrated by the single-operator analysis used in Ref. [2] to determine the energy of the stationary I = 0, ground state.   One strategy to resolve this problem is to fix one or more of the parameters in the extra-state fit. The first parameter we might fix is the energy of the extra state.
A prediction of this finite-volume energy can be obtained by finding the intersects between Lüscher's formula and a phenomenological model, or a fit to our lattice results, for the scattering phase shift as a function of energy. For simplicity we chose to compute this estimate using the dispersive predictions of Ref. [13]. Unfortunately we found that fixing this extra-state energy alone typically does not solve the problem: sometimes the fitting procedure continues to fail while in those cases where the fit converges, it gives a groundstate energy whose statistical error is several times larger than the statistical error on the result from the preferred fit, which suggests that such an approach may over-estimate the excited state error, conflating it with the statistical error.
Thus, further assumptions are needed to obtain a statistically meaningful estimate of this excited-state contamination error. Because of the very different pattern of operatoreigenstate overlaps we will adopt two different strategies, one for the I = 0 stationary ππ state and the second for the remaining seven cases: the three I = 0 calculations with nonzero total momenta and all four I = 2 calculations. The method applied for the latter we will refer to as "method A", and that applied in the special I = 0 stationary case as "method B".
For the three moving-frame I = 0 and all the I = 2 calculations, the operator-eigenstate overlap matrix is close to diagonal and the two-point function of a given operator with itself is well described by one exponential coming from a single energy eigenstate. We call the operator which couples primarily to the ground state the "ground-state operator", and the operator whose two-point function is dominated by the n th excited state the "n th excited operator". We can then make the reasonable assumption, that because of the small value of the overlap factors between the ground-state operator and those excited states that we include in the preferred fit, the extra-state contribution to these small overlap factors can be neglected. That means that we can focus on the Green's function constructed from the ground-state operator only where the most important effect of an omitted extra state is to change the two parameters (overlap factor and energy) associated with the ground-state.
The argument above suggests that we can perform a fit to the two-point function constructed from the ground-state operator in which the number of states included is one more than the number in the preferred fit, and that we can fix the overlap factors between the ground-state operator and the excited states in that fit to those already determined by the preferred fit. Since these overlap factors are small, we can also fix the energies of those ex-

G. Error budget
We will now combine all of the errors detailed in the earlier sections to provide values for the seven ππ phase shifts that we have computed at specific energies and their corresponding errors. We divide these errors into two categories. The first are the measurement errors I=2 channel P tot = (2, 2, 2) π L P tot = (2, 2, 0) π L P tot = (2, 0, 0) π L P tot = (0, 0, 0) π L Fit range 10-25 12    We refer to the second category of errors as tuning errors. These include the finite volume error, finite lattice spacing error, the unphysical kinematics error, and the higher partial wave correction error. We assign these errors directly to the seven phase shifts. They represent discrepancies between our physical results as presented and those we actually obtain. For example, we describe our results as phase shifts in the continuum limit, but in reality the number we obtain contains finite lattice spacing errors. Similarly we intend our phase shift results to be for the case where m π = 135 MeV. In reality, our calculation was done for a different pion mass and we have made a small correction to the energy at which the phase shift is quoted to compensate for the shifted ππ threshold arising from our incorrect mass and, as was discussed in Section VII D, the remaining errors were estimated using ChPT.
Our first category of errors, the measurement errors, requires special treatment if we are to account for the tight correlation between the way the errors on the measured finitevolume energy shift affect both the implied value of the phase shift and the energy at which that phase shift is quoted. The issue is that we do not know the exact value of the finite volume ππ energy E FV , only that it lies with a certain confidence within the rangē We imagine that the true intersect occurs at some energy E that is close to our measured central value E FV and perform a first-order Taylor expansion in E of both sides of Eq. (70) about this central value Rearranging the above we obtain an expression for the phase shift evaluated at our measured finite-volume energy E FV This approach requires an estimate of the derivative of the phase shift δ(E) with respect to its energy. While this could be obtained from a fit to our data, we find it simplest to use the result from the dispersive analysis [13] which agrees well with our data as can be seen from the comparisons shown in Figures 13 and 14.
Thus, the statistical and excited state contamination errors on the measured finite-volume energies are each converted to errors on the phase shift at the fixed energy E = E FV using the following relations: where ∆E stat and ∆E exc are the statistical and excited state contamination errors that are assigned to the measured finite-volume energy. The three tuning errors ∆δ dis , ∆δ FV , ∆δ unphy coming from non-zero lattice spacing, finite volume and unphysical pion mass are assigned directly to the phase shifts. All of these errors are listed in Table XVII, where the total systematic error is the combination in quadrature of the systematic errors shown in columns 5 through 9.
Another topic that we should discuss is the systematic error on the scattering length of both I = 0 and I = 2 channels. Similar to the discussion above, the systematic error is dominated by the excited state error, so we neglect all other sources of systematic error.
The results are summarized in Tab. XVIII. It can be shown that, despite the fact that the scattering length of I = 0 channel is inconsistent with the experimental value if we only consider the statistical error, as we show in Sec. VI C, the final systematic error for I = 0 channel is three times larger than the central value, and by including this error, the lattice results are consistent with the experimental value. Our large excited state error in moving frame calculation, and the resulting large relative error on the energy difference between E ππ and 2m π , leads to such a large error. of energy. Shown also on the plot are the corresponding dispersive results [13]. Note: the errors are not shown for the dispersive results.

VIII. CONCLUSIONS
In this paper we have presented in detail a lattice calculation of the ππ scattering phase shift for both the I = 0 and I = 2 channels. Our final results are presented in Table XVII and illustrated in Figures 13 and 14. This calculation is performed using a physical pion mass with G-parity boundary conditions in all three spatial directions. This results in a single-pion, finite-volume ground state with non-zero momentum, which is crucial for the closely related calculation presented in Ref. [1] for the I = 0 K → ππ decay amplitude and results [13] that are shown in Fig. 13, but here with an expanded scale.
interpolating operators, one of which is a four-quark operator constructed from two pion interpolating operators which each carry larger-than-minimum momenta, while the other is a scalar two-quark operator (the sigma operator). With these improvements we obtain an error is improved by a factor of 5. ii) We are able to provide a more reliable and detailed systematic error analysis. iii) We have been able to resolve the 3σ discrepancy between our earlier result for this phase shift and that predicted by a dispersive analysis [13] so that our current results agree well with the dispersive prediction (cf. Fig. 13). The discrepancy is now 4 Note this number is slightly different from the number given in Ref. [1]. This is because we have now included an estimate of the error due to the unphysical pion mass, resulting in a correction to the ππ energy and we have refined our excited state error estimation.  Here a 0 (exp) is the experimental value of the scattering length quoted from Ref. [11] understood to have resulted from excited ππ state contamination, which was underestimated in Ref. [2] and is now under much greater control.
More specifically, we have employed a concrete procedure for estimating the error resulting from a nearby excited state that was not included in our fit. As discussed in Sec. VII, we introduce one additional state into our fit but with an energy fixed to that given by the dispersive calculation [13] and with overlap factors to our operators carefully estimated, so as to avoid introducing instability in the fits or inflating the statistical error. The resulting shift in the ground-state energy then provides a meaningful indicator of the size of the corresponding systematic error.
In addition to computing the I = 0 phase shift for two pions with zero total momentum, we also perform a moving-frame calculation with three different total momenta. The observation [49] that three types of lattice symmetry can be used to significantly reduce the number of contractions was exploited to reduce the contraction time by a factor of seven.
The resulting values of the ππ phase shifts at lower energies not only allow us to perform a comparison with dispersive and chiral perturbation theory predictions but also give us an independent evaluation of the Lellouch-Lüscher correction needed to obtain the K → ππ decay amplitude from a finite-volume lattice QCD calculation. Because of the critical role played by the sigma interpolating operator in the stationary frame calculation, we will include a sigma operator with non-zero total momentum in future work. This operator might be expected to strongly couple to more states in the fit, in contrast to the ππ(311, 311) operator, and may significantly reduce the errors as it did for the stationary case. An additional future goal is to extrapolate the result to the continuum limit, since in this work, the calculation is performed on a single 32 3 lattice. We are now preparing such a scaling calculation using two lattice volumes with spatial extents of 48 3 and 64 3 .
In this paper, we have used a combination of bootstrap and jackknife methods [48] together with correlated fits to determine the ground-state energies, the operator-state overlap amplitudes and the analysis of the excited-state error. This method allows us to estimate the goodness-of-fit, which provides guidance as to which model and fitting range we should choose (see Sec. V). We also compared our multi-state fitting with the GEVP method and found that our fitting procedure gave consistent statistical errors for the I = 0 case than did our implementation of the GEVP method. We did not attempt to estimate the systematic errors resulting from the GEVP approach. The GEVP method may well excel when more states and operators are included.

ACKNOWLEDGMENTS
We would like to thank our RBC and UKQCD collaborators for their ideas and assistance and to specifically acknowledge the contributions of Chulwoo Jung and Taku  As shown in Ref. [6] the allowed quark momenta (for G-parity BCs in 3 directions) are ± π 2L (1, 1, 1) + 2π L n , where n is a vector of integers. While any combination of p and q satisfying this condition result in valid pion interpolating operators, we observed in Ref. [6] that the cubic symmetry breaking manifest in the operator amplitudes between pion states of total momentum related by cubic rotations is dramatically suppressed by averaging over pairs of bilinear operators with the same total momentum but with different assignments of quark momenta. The specific criteria for selecting those momenta are discussed in more detail in that paper; here in Table XIX we list only the two choices for each of the 32 total momenta. The momentum distribution is listed below for all 32 pions we use (in units of π/2L): Recall that in addition to the above, we also symmetrize the momentum between the quark and antiquark by averaging the assignments ( p, q) and ( q, p). Thus in practice our pion interpolating operators comprise an average over a total of four quark field bilinears.

σ operator
In this work we use the σ operator only in the case of zero total momentum, and as a result the momenta assigned to the quark and anti-quark fields must be equal and opposite.
We construct an operator that is symmetric under cubic rotations by averaging over 8 orientations of the quark momentum. The list of momenta assigned to the quark operator are given in Table. XX.

Appendix B: Chpt prediction for phase shift
In this appendix we present the partial wave amplitude t I l=0 results from the next-toleading-order (NLO) ChPT. These amplitudes are connected with the scattering phase shift In this appendix we list the contraction formula for each diagram introduced in Sec. III.
The first four diagrams are associated with the product of two ππ interpolating operators, where the four time slices are the time coordinates of the four single-pion interpolating operators, which are t src − 4, t src , t snk and t snk + 4, respectively. The final four expressions correspond to the cases where at least one of the source or sink operators is a σ operator.
The quantity P ta,t b is the G-parity quark propagator from t a to t b while the flavor-spin matrix S 1 is defined as S 1 = σ 3 γ 5 . These eight amplitudes are obtained from the following contractions: T r {P t 2 ,t 4 S 1 P t 4 ,t 2 S 1 } + (t 3 ↔ t 4 ) (C2) R = 1 2 1 2 T r {P t 1 ,t 2 S 1 P t 2 ,t 3 S 1 P t 3 ,t 4 S 1 P t 4 ,t 1 S 1 } + (t 3 ↔ t 4 ) (C3)