\section{Introduction}
\label{sec:intro}

Understanding the interactions of two pions is an interesting endeavor for practitioners of non-perturbative QCD. Not only do we learn how the fundamental interactions of quarks and gluons give rise to the observable properties of hadrons, these two particle systems play an important role in Standard Model processes under intense investigation, such as $K\to\pi\pi$ decays~\cite{Blum:2015ywa,Bai:2015nea,RBC:2020kdj} and the muon's anomalous magnetic moment $g-2$~\cite{RBC:2018dos,Bruno:2019nzm}. Our focus in this study is primarily on isospin $I=0$ and 2 for the former while $I=1$ is important for the latter. Isospin symmetry and Bose-Einstein statistics constrain the states that appear in these processes.

With L\"uscher's technique \cite{Luscher:1990ux} that relates two-pion energy in a finite box with the corresponding scattering phase shift, there have been many studies of two-pion scattering in lattice QCD at unphysical pion masses \cite{Kuramashi:1993ka,CP-PACS:2004dtj,CP-PACS:2005gzm,Beane:2005rj,Beane:2007xs,Feng:2009ij,NPLQCD:2011htk,Dudek:2012gj,Bulava:2016mks,Briceno:2016mjc,Liu:2016cba,Fu:2017apw,Culver:2019qtx}.  For these studies, a chiral extrapolation was needed to obtain physical results.
The analytic evaluation of two-pion scattering in chiral perturbation theory \cite{Bijnens:1995yn,Bijnens:1997vq} (ChPT) was employed in these works to perform the extrapolation of important parameters of two-pion scattering such as the scattering length.
While we can expect the extrapolation is reasonable for the scattering length, which can be extracted near two-pion threshold, the extrapolation of the scattering amplitudes or phase shifts might not be accurate at high energies.

Now it is possible to perform a lattice calculation at the physical pion mass so that we can directly compute the two-pion phase shifts at relatively large energies without a chiral extrapolation.
There was a study where the $I=2$ scattering length was computed including the physical pion mass \cite{Fischer:2020jzp} but no other except the series of our two-pion scattering studies \cite{RBC:2021acc,Meyer:2023prep}.
Here we report results for phase shifts at various energy levels for $I=2$ and 0 at physical pion mass using $2+1$ flavors of M\"obius domain wall fermions (MDWF). This is one of several related projects undertaken by the RBC and UKQCD collaborations.

The RBC and UKQCD collaborations have reported results for CP violation in kaon decays and pion phase shifts at the physical point for the corresponding $I=0$ and 2 final states~\cite{Bai:2015nea,RBC:2020kdj,RBC:2021acc}. Because the physical kinematics for such decays require pions with back-to-back relative momenta, which is not the ground state achieved in ordinary lattice calculations where the pions are at rest \cite{Maiani:1990ca}, G-parity spatial boundary conditions (GPBC) were employed in the simulations~\cite{Christ:2019sah}. GPBC forbid pions with zero momentum, and if the box size is adjusted appropriately, then the two-pion ground state computed on the lattice will have physical momenta satisfying $E_{\pi\pi} \approx M_K$. In the GPBC two-pion scattering work \cite{RBC:2021acc}, we obtained the phase shifts of the $I=0$ channel as well as the $I=2$ channel at various two-pion energies with non-zero pion momenta that are consistent with the prediction from the dispersion theory \cite{Ananthanarayan:2000ht,Colangelo:2000jc,Colangelo:2001df,Garcia-Martin:2011iqs,Garcia-Martin:2011nna} based on the Roy Equation \cite{Roy:1971tc} with inputs obtained by a combination of chiral perturbation theory for the scattering lengths and experimental data for the high energy regime.
GPBC are not implemented without extra cost, however. They are at least twice as expensive for measurements compared with periodic boundary conditions because the G-parity Dirac operator is explicitly two-flavor, with mixing between the flavors occurring at the boundary, and they also require gauge ensembles with the same boundary conditions to be generated. In addition, an important future step is to include isospin breaking effects in the calculation of $\varepsilon^\prime$ which is expected to be significant but GPBC may not be suitable due to the intrinsic role of the isospin symmetry.

The long-term aim of this study is to explore the use of periodic boundary conditions (PBC) to answer the question whether the decay amplitudes with physical kinematics can be extracted reliably from an {\it excited state} computed on the lattice.  As a first step, we investigate pion scattering in this setup.

Two-pion states with a definite total momentum can vary their total energy not only by a standard excitation of a single pion but also by changing the momenta of individual pions, or equivalently the relative momentum.  While a finite box allows the pions to have discrete momenta, $2\pi n/L$ in the case of PBC, labeled by an integer $n$, the interval of momenta $2\pi/L$ and hence that of two-pion state energies have to be small when we employ a sufficiently large volume $L$ to introduce the physical pion mass with controllable finite volume effects.
It may thus be challenging especially for $I=0$, where there are disconnected diagrams and corresponding operators couple with the vacuum state, with the physical pion mass to extract the signals of an excited state with the statistical and systematic errors under control.
In fact, we learned from our earlier works \cite{Bai:2015nea,RBC:2020kdj,RBC:2021acc} that there is significant higher-state contamination in two-pion correlation functions that also needs to be removed.  The higher-state contamination appears quite significant for the $I=0$ channel at physical pion mass.

In this work we attack the problem with several cutting-edge lattice methods: all-mode-averaging (AMA)~\cite{Blum:2012uh,Shintani:2014vja}, all-to-all (A2A) propagators~\cite{Foley:2005ac}, and the variational method using the generalized eigenvalue problem (GEVP)~\cite{Luscher:1990ck,Blossier:2009kd} to compute correlation functions and extract energy eigenvalues.
A companion paper \cite{Meyer:2023prep} using distillation \cite{HadronSpectrum:2009krc} will also be available soon.

The GEVP method provides us with a systematic procedure to decompose correlation functions into contributions from the several lightest states with the same quantum numbers and have been widely used for hadron spectrum studies.
In our particular case, it turned out from our earlier works that introducing a $\sigma$ operator for $I=0$ in the measurements plays a more important role in removing the contamination from excited states than simply increasing statistics.
Originally \cite{Bai:2015nea} we computed $I=0$ two-point functions using a single two-pion operator with 216 configurations and resulted in unexpectedly large contamination from excited states.
Later \cite{RBC:2021acc} we performed multi-state fits to a matrix of correlation functions calculated with three interpolation operators including a $\sigma$ operator as well as another additional two-pion operator and then the obtained energy of the two-pion ground state dramatically changed with 741 configurations, while increasing statistics to 1,400 configurations with the single two-pion operator still gave a ``fake plateau'' and increased the difference in the phase shift from the dispersive approach.
In addition, obtained two-pion energy from the variational method had much less statistical error.
In this work, we introduce a $\sigma$ operator as well as four two-pion operators with various pion momenta for our measurements and GEVP analysis.

In addition, we propose a variant of the GEVP approach which we call the {\it re-based} GEVP (RGEVP).  The eigenvectors of GEVP obtained at a time slice give us a new basis of operators.  In principle, each operator in the new basis couples well with one of the lowest energy states considered in the GEVP.  With limited statistics, since we could lose the signal from one or more of those states at large time separations, it may be reasonable to exclude such noisy states by removing the corresponding operators from the basis so that all the states included in the GEVP analysis have good statistical precision.  The RGEVP is to reduce the size of GEVP by using fewer operators that couple well with states.
We find that this approach gives us an improvement on statistical precision for the ground and first excited states of the $I=0$ channel.

We perform lattice calculation for two-pion scattering with 258 configurations on the $24^3\times64$ lattice with the lattice cutoff $a^{-1} = 1.023$~GeV and 107 configurations on the $32^3\times64$ lattice with $a^{-1} = 1.378$~GeV
\cite{RBC:2014ntl,Tu:2020vpn}.  Both ensembles are generated with $2+1$-flavor M\"obius domain wall fermions and Iwasaki $ + $ DSDR (dislocation suppressed determinant ratio) gauge action.
While better statistical precision is desired and we will update our results in the near future, we take the continuum limit of the scattering phase shifts and scattering lengths at this point.  Our determination of the scattering length does not need a chiral extrapolation, which assumes the leading order as an input from ChPT and hence gives a precise value.  Our results are meaningful as a pure lattice determination, though they have larger uncertainty.

Although the core of this study is the application to $K \to \pi \pi$ and the direct CP violation parameter, $\varepsilon'$, the experience gained here will provide impetus to other applications
of $\pi$-$K$ scattering phases. Examples that we have in mind so far are
direct CP violation in charm decays~\cite{LHCb2019,AS_LAT19}, possible CP violation in $\tau \to \nu K \pi$~\cite{Babar_tauCP}, and
three body proton decays.

This paper is organized as follows. Section~\ref{sec:th} describes the theoretical framework underlying the calculation. In Section~\ref{sec:lat details} we give the lattice details. Section~\ref{sec:results} gives results for the computed two-pion energies and the corresponding phase shifts. Here we also compare our results to recent data-driven studies~\cite{Ananthanarayan:2000ht}. In Section~\ref{sec:comparison} we compare the PBC calculation to the GPBC one~\cite{RBC:2021acc}.
Section~\ref{sec:conclusion} summarizes the present work and future prospects.

\section{Theoretical framework}
\label{sec:th}

\subsection{Operator Construction}
\label{sec:operator def}

In this subsection we describe the operators and states used in this work.

We start with pion operators with definite spatial momentum
\begin{align}
    \pi^a(t,\vec p)
    = \sum_{\vec x,\vec y}\,
    &e^{-i(\vec p_1\cdot\vec x + \vec p_2\cdot\vec y)}
    f_r(||\vec x - \vec y||)
    \notag\\
    &\times\overline\psi(t,\vec x){\rm{i}}\gamma_5F^a\psi(t,\vec y),
    \label{eq:pionOp}
\end{align}
which is defined with the Coulomb gauge fixing and the center-of-mass momentum $\vec p = \vec p_1 + \vec p_2$.
Although the pion operators defined in this way depend on both quark and anti-quark momenta $\vec p_1$ and $\vec p_2$, respectively, we always give zero momentum to the anti-quark ($\vec p_1 = 0$) so that the quark carries the total momentum ($\vec p_2 = \vec p$) in this work and we add only the total momentum as an argument of the operators.
The exponential smearing function is defined as
\begin{equation}
    f_r(||\vec x-\vec y||) = \exp(-||\vec x-\vec y||/r),
\end{equation}
with smearing radius $r$ and the periodic modulus $||\vec x -\vec y||$, the length of the shortest straight path from $\vec y$ to $\vec x$ in the periodic box. This hydrogen-like wave function has been used in our earlier works ~\cite{RBC:2020kdj,RBC:2021acc}.
While the single-pion operator of course depends on the smearing radius, we drop $r$ from the pion operator on the left hand side of Eq.~\eqref{eq:pionOp} for simplicity.
The quark and anti-quark isospin doublets are defined as
\begin{equation}
    \psi
    = \left(
    \begin{array}{c}
    u \\ d
    \end{array}
    \right),\hspace{10mm}
    \overline\psi = (
    \begin{array}{cc}
    \bar u & \bar d
    \end{array}
    ),
\end{equation}
and
\begin{align}
    F^+ &= \frac{1}{2}(\sigma_1 + i\sigma_2),\\
    F^- &= \frac{1}{2}(- \sigma_1 + i\sigma_2),\\
    F^0 &= \frac{1}{2}\sigma_3,
\end{align}
with the Pauli matrices $\sigma_{1,2,3}$.

The two-pion operators are constructed by multiplying two single-pion ones:
\begin{align}
\widetilde O_{\pi\pi}^{I,I_z}(t_1, t_2,\vec P, \vec p/2)
= \sum_{a,b}c^{I,I_z}_{ab}\, &\pi^a(t_1,(\vec P+\vec p)/2)
\notag\\
\times\, & \pi^b(t_2, (\vec P-\vec p)/2),
\end{align}
where $\vec P$ and $\vec p$ are the center-of-mass and relative momenta of the two pseudoscalar operators, respectively, and $a,b\in \{+,-,0\}$.
The coefficients $c^{I,I_z}_{ab}$ project the two-pion operator to an isospin-definite channel labeled by $(I,I_z)$.  Appendix~\ref{sec:contractions} gives the explicit forms of the $(I,I_z) = (2,0)$ and $(0,0)$ two-pion operators.

The discrete, finite lattice breaks the continuum rotational symmetry of angular momentum SO(3) down to a discrete subgroup, which depends on the center-of-mass momentum.  The irreducible representations of such a discrete subgroup do not give rise to angular momentum eigenstates that appear as irreducible representations of SO(3). Instead, they are mixtures which can be classified in terms of the continuum irreducible representations.

It is fairly straightforward to make this classification based on fundamental group theory for one and two particle systems, both moving and at rest~\cite{Moore:2005dw,Dudek:2012gj}.
Our main target is the $s$-wave two-pion states and their phase shifts. The corresponding interpolating operators are defined as
\begin{equation}
O_{\pi\pi}^{I,I_z}(t_1, t_2, \vec P, \vec p/2)
= \sum_{\hat T\in G}
\chi_{A_1}(\hat T)
\widetilde O_{\pi\pi}^{I,I_z}(t_1, t_2,\vec P, \hat T[\vec p/2]),
\label{eq:pipi_swave}
\end{equation}
where we sum over all elements $\hat T$ in the finite-volume symmetry group $G$ and the normalization factor $\chi_{A_1}(\hat T)$ is the character of the group element $\hat T$ in the representation $A_1$~\cite{Moore:2005dw,Dudek:2012gj}.

While one could consider composite operators of two quark bilinear operators, both located at a single time slice, it has been shown that placing them on slightly different time slices is advantageous for reducing statistical noise.  This is the case especially for $I=0$~\cite{Bai:2015nea,RBC:2020kdj,RBC:2021acc}, where the overlap of the two-particle operator with the vacuum state can be suppressed exponentially by the separation $\Delta\equiv |t_2 - t_1|$.

In addition to two pions, we introduce a $\sigma$ operator, or iso-singlet scalar bilinear operator, for $I=0$:
\begin{equation}
    \sigma(t, \vec p)
    = \sum_{\vec x,\vec y}\,
    e^{-i(\vec p_1\cdot\vec x + \vec p_2\cdot\vec y)}
    f_r(||\vec x - \vec y||)
    \overline\psi(t,\vec x)\psi(t,\vec y),
\end{equation}
with $\vec p = \vec p_1 + \vec p_2$.
Again, we set $\vec p_1 = \vec0$ and $\vec p_2 = \vec p$ in this work.
This operator has been found to play an important role in controlling the contamination from excited states~\cite{RBC:2020kdj,RBC:2021acc}.

In this work, we concentrate on the rest frame $\vec P = \vec 0$ and $s$-wave operators and states.
For $I=2, I_z=0$, we consider four values of relative pion momenta and use the following operator basis:
\begin{equation}
    O^{2,0}(t)
    = \left(
    \begin{array}{c}
    O_{\pi\pi}^{2,0}(t,t+\Delta,\vec 0,(0,0,0)\times2\pi/L) \\
    O_{\pi\pi}^{2,0}(t,t+\Delta,\vec 0,(0,0,1)\times2\pi/L) \\
    O_{\pi\pi}^{2,0}(t,t+\Delta,\vec 0,(0,1,1)\times2\pi/L) \\
    O_{\pi\pi}^{2,0}(t,t+\Delta,\vec 0,(1,1,1)\times2\pi/L)
    \end{array}
    \right).
    \label{eq:basisI2}
\end{equation}
Similarly for $I=0,I_z=0$, we define
\begin{equation}
    O^{0,0}(t)
    = \left(
    \begin{array}{c}
    O_{\pi\pi}^{0,0}(t,t+\Delta,\vec 0,(0,0,0)\times2\pi/L) \\
    \sigma(t, \vec 0) \\
    O_{\pi\pi}^{0,0}(t,t+\Delta,\vec 0,(0,0,1)\times2\pi/L) \\
    O_{\pi\pi}^{0,0}(t,t+\Delta,\vec 0,(0,1,1)\times2\pi/L) \\
    O_{\pi\pi}^{0,0}(t,t+\Delta,\vec 0,(1,1,1)\times2\pi/L)
    \end{array}
    \right).
    \label{eq:basisI0}
\end{equation}
The order of the operators in these bases matters when it comes to $N\times N$ GEVP as described in Section~\ref{sec:GEVP}.
When we do not need to specify the isospin, we simply call the four two-pion operators $\pi\pi(000)$, $\pi\pi(001)$, $\pi\pi(011)$ and $\pi\pi(111)$, respectively.

\subsection{Correlation Functions}

For the GEVP a matrix correlation function is defined in the basis of operators given in the previous sub-section,
\begin{equation}
    C_{ij}^{I,I_z}(t)
    = \left\langle
    O^{I,I_z}_i(t) O^{I,I_z}_j(-\Delta_j)^\dag
    \right\rangle,
    \label{eq:def_2ptCorr}
\end{equation}
where
\begin{equation}
\Delta_j = \Bigg\{
\begin{array}{ll}
\Delta & (I=2)
\\[2mm]
(1-\delta_{2,j})\Delta & (I=0)
\end{array},
\end{equation}
translates the source two-pion operators by $\Delta$ but does nothing for the $\sigma$ operator so that the time variable $t$ always indicates the minimum time separation between a bilinear of the source operator and that of the sink operator.  While the measured correlator matrix is not exactly symmetric with finite statistics, we symmetrize it by averaging with the transposed partner.

The Wick contractions of the two-point functions are shown diagramatically in Figs.~\ref{fig:Wick contractions} and \ref{fig:Wick contractions sigma} .  They are denoted ``direct" (D), ``cross" (C), ``rectangle" (R), and ``vacuum" (V).
By taking an appropriate linear combination of these elemental contractions, we construct a correlation function of operators carrying definite isospin.
The $I=2$ channel is contributed by diagrams D and C, while all the diagrams in Figs.~\ref{fig:Wick contractions} and \ref{fig:Wick contractions sigma} contribute to the $I=0$ channel.
We compute the disconnected diagrams for every time translation and take the translation average, while the connected diagrams are computed at every several source time slices, which are specified in Section~\ref{sec:lat details}.
The complete formulae for the $I=2$ and $I=0$ channels are given in Appendix~\ref{sec:contractions}.

For $I=0$ there is an additional complication:
in the rest frame the ground state in this channel is the vacuum. This contribution dominates the correlation function and must be subtracted,
\begin{align}
    C_{ij}^{0,0}(t)
    &= \left\langle
    O^{0,0}_i(t) {O^{0,0}_j(-\Delta_j)}^\dag
    \right\rangle
    \notag\\
    & - \frac{1}{L_t}\sum_{t_{src}=0}^{L_t-1} \left\langle
    O^{0,0}_i(t+t_{src}) \right\rangle
    \left\langle {O^{0,0}_j(t_{src}-\Delta_j)}^\dag
    \right\rangle,
    \label{eq:vac_subt}
\end{align}
where $L_t$ stands for the time extent of the lattice ensemble in lattice units.
While the second term on the right hand side is independent of $t$ in the limit of infinite statistics, we perform this subtraction time-slice-by-time-slice as we found it provides a minor statistical advantage~\cite{RBC:2020kdj,RBC:2021acc}.

\begin{figure}[tbp]
    \centering
    \includegraphics[width=.22\textwidth]{figs/pipi_D}~~~~~~~

    \includegraphics[width=.22\textwidth]{figs/pipi_C}
    \\[6mm]
    \includegraphics[width=.22\textwidth]{figs/pipi_R}~~~~~~~

    \includegraphics[width=.22\textwidth]{figs/pipi_disc}
    \caption{Elemental Wick contractions. Clockwise from top-left: D (direct), C (cross), V (vacuum), and R (rectangle). Linear combinations of the four diagrams are used to construct pion scattering correlation functions for $I=0$, 1, and 2 states with definite lattice hypercubic symmetry.}
    \label{fig:Wick contractions}
\end{figure}

\begin{figure}[tbp]
    \centering
    \vspace{6mm}
    \includegraphics[width=.22\textwidth]{figs/pipi_sigma_conn.pdf}~~~~~~~

    \includegraphics[width=.22\textwidth]{figs/pipi_sigma_disc.pdf}
    \\[6mm]
    \includegraphics[width=.22\textwidth]{figs/sigma_conn.pdf}~~~~~~~

    \includegraphics[width=.22\textwidth]{figs/sigma_disc}
    \caption{Elemental Wick contractions with one (upper) or two (lower) scalar bilinear operators. These are analoguous to the R (left) and V (right) diagrams in Fig.~\ref{fig:Wick contractions} and relevant for only the $I=0$ channel.}
    \label{fig:Wick contractions sigma}
\end{figure}

\subsection{Thermal effects}

Due to the finite time size of the lattice ($L_t$) and the pions satisfying periodic boundary conditions in time, unwanted contributions contaminate the correlation function. These so-called around-the-world (ATW), or thermal, effects arise when one of the source pions propagates forward in time while the other goes backwards through the boundary to reach the sink time slice.  They can be seen by inserting a complete set of states in the two-point (``thermal") correlation function and translating the source and sink operators to equal times.
\begin{align}
&\braket{O_{\pi\pi}(t)O_{\pi\pi}(0)^\dagger}
\notag\\
=&
\sum_m\sum_n \braket{m|O_{\pi\pi}(t) \ket{n}\bra{n}O_{\pi\pi}^\dagger(0)|m}\nonumber
\\ =&\sum_n (e^{-E_n^{\pi\pi} t}+e^{-E_n^{\pi\pi} (L_t-t)}) \braket{0|O_{\pi\pi} \ket{n}\bra{n}O_{\pi\pi}^\dagger|0}\nonumber
\\ &+e^{-E_{\pi(\vec p)} t}e^{-E_{\pi(\vec p)}(L_t-t)}
\notag\\&\hspace{5mm}\times
\braket{\pi(\vec p) \left| O_{\pi\pi}\right|\pi(\vec p)}\braket{\pi(\vec p) \left| O^\dagger_{\pi\pi}\right|\pi(\vec p)}+\dots,
\label{eq:atw}
\end{align}
where we omit the isospin superscripts and momentum arguments for simplicity and the sum over $m$ gives the thermal expectation value.  The first term on the right hand side contains the zero temperature expectation value, while the last term is the thermal contribution which vanishes as $L_t\to \infty$.  Notice when the rest frame is employed, the leading ATW contribution is time independent while the ATW effects of an excited-state pion and those in a moving frame are time dependent. All thermal effects are suppressed exponentially with $L_t$.

For the rest frame, the leading ATW term, which is constant, can be removed simply by a ``matrix" subtraction
\begin{align}
    C_{ij,subt}^{I,I_z}(t)&\equiv C_{ij}^{I,I_z}(t)-C_{ij}^{I,I_z}(t+\delta_t).
    \label{eq:matdt_subt}
\end{align}
where $\delta_t$ is an arbitrary time shift.
This subtraction removes all constant contributions to the correlation functions and therefore the vacuum subtraction \eqref{eq:vac_subt} is in principle unnecessary if this subtraction is applied.
In this work, we still apply the vacuum subtraction so that we can investigate the significance of the ATW effects with the absence of vacuum effects by analyzing both matrix-subtracted and unsubtracted correlators.

\subsection{Generalized eigenvalue problem method}
\label{sec:GEVP}

For the rest of the section we omit the superscripts $I$ and $I_z$ and simplify our notation of the correlator matrix $C_{ij,subt}^{I,I_z}(t) \to C_{ij}(t)$, or $C(t)$ when the operator indices $i$ and $j$ can be dropped without confusion.

As is well known, Euclidean space correlation functions are a sum of exponential terms, each term corresponding to one state in a tower of states with fixed quantum numbers and increasing energies,
\begin{align}
C_{ij}(t)
& = \sum_n A_{n,i}A_{n,j}^*e^{-E_n t},
\label{eq:state_expansion}
\end{align}
where we have neglected the backward propagating contributions proportional to $e^{-E_n(L_t-t)}$ since $t\ll L_t$ in our setup\footnote{This contribution is exponentially suppressed compared to the ATW effect considered earlier.}. $A_{n,i}=\langle 0|O_i|n\rangle (1-e^{-E_n\delta_t})^{1/2}$ is the overlap of the $i$-th operator acting on the $n$-th state and the vacuum multiplied with with a normalization factor due to the matrix subtraction.

To extract the desired excited states in the correlator, we employ the variational method by solving a generalized eigenvalue problem (GEVP)~\cite{Luscher:1990ck,Blossier:2009kd}.
For an $N\times N$ matrix $C(t)$ we solve the following GEVP:
\begin{eqnarray}
C(t)V_n(t,t_0) &=& \lambda_n(t,t_0) C(t_0)V_n(t,t_0),
\label{eq:GEVP}
\end{eqnarray}
with eigenvalues $\lambda_n(t,t_0)$ and eigenvectors $V_n(t,t_0)$, where in principle we can choose $t_0$ in the range $0 < t_0 < t$.
At asymptotically large time separations the eigenvalue behaves as $\lambda_n(t,t_0) = e^{-E_n (t-t_0)}$ where $E_n$ is $n$-th energy state in the GEVP. In Ref.~\cite{Blossier:2009kd} it was shown that the leading correction behaves like $e^{-(E_{N+1}-E_n)t}$ for $t_0\ge t/2$.  In this work, we use the first $N$ operators of the bases in Eqs.~\eqref{eq:basisI2} and \eqref{eq:basisI0} for the $I=2$ and $I=0$ channels, respectively, and solve GEVP with a $t$-independent value of $t - t_0\equiv\Delta_t$.  While the preferable inequality $t_0\ge t/2$ is violated in the region $t<2\Delta_t$, we do not use data at such short times relative to $\Delta_t$ for our final results.

Effective two-pion energies are defined as~\cite{Luscher:1990ck,Blossier:2009kd}
\begin{equation}
     E^{eff}_n(t,t_0) = \ln{\lambda_n(t,t_0)}- \ln{\lambda_n(t+1,t_0)}.
     \label{eq:def_Eeff}
\end{equation}
The corresponding eigenvector $V_n(t,t_0)$ at asymptotic time separation provides a new operator that couples to the $n$-th state but not with the other states in the GEVP\footnote{The coupling is not perfect: in Ref.~\cite{Blossier:2009kd} it is shown that corrections are $O(e^{-(E_{N+1}-E_n)t_0})$.}~\cite{Blossier:2009kd}
\begin{equation}
    \tilde O_n = \sum_i V_{n.i} O_i.
    \label{eq:op eigenvector}
\end{equation}
These eigenvectors play a key role in isolating the weak operator matrix elements between an excited two-pion state and the kaon in $K\to\pi\pi$ decays, for example.

In practice, the large statistical error of correlation functions at large time separations may cause misordering of eigenvalues and eigenvectors for specific jackknife samples resulting in incorrectly large errors in two-pion energies and GEVP eigenvectors.
We attempt to avoid the misordering as follows.

At short time separations where the statistical errors are small enough, sorting eigenvalues into descending order will give us the ascending order of effective energies.  We first solve the $N\times N$ GEVP \eqref{eq:GEVP} at $t=\Delta_t+1$ and $t_0 = 1$ and obtain $N$ eigenvectors
\begin{equation}
V_1(\Delta_t+1),V_2(\Delta_t+1),\ldots,V_N(\Delta_t+1),
\end{equation}
which are ordered by ensuring the descending order of the corresponding eigenvalues. The second argument of the eigenvectors is omitted for simplicity since it is always $\Delta_t$ smaller than the first argument in this work as explained above.
To ensure the correct order also at larger time separations where one or more states may not have a good signal, we employ a recursive procedure as follows.  For $t\ge\Delta_t+2$, we now suppose the ordering of the eigenvectors $V_n(t-1)$ at $t-1$ is successful and give a recipe to obtain the correct order of the eigenvectors $V_n(t)$ at $t$.
We define an $N\times N$ matrix
\begin{eqnarray}
    T(t-1) = (V_1(t-1)\,\,V_2(t-1)\,\ldots V_N(t-1)),
\end{eqnarray}
using the set of the GEVP eigenvectors $V_n(t-1)$ obtained at one time slice earlier.
Then we can calculate approximately diagonal matrices
\begin{align}
    C'(t, t-1) &= T(t-1)^\dag C(t) T(t-1),
    \label{eq:re-basedCt}\\
    C'(t_0, t-1) &= T(t-1)^\dag C(t_0) T(t-1).
    \label{eq:re-basedCt0}
\end{align}
Here, the second argument of $C'$ in the left hand sides corresponds to the argument of $T$ on the right hand sides.  With these near diagonal matrices, we can express the GEVP \eqref{eq:GEVP} as
\begin{align}
C'(t, t-1)& V_n'(t) = \lambda_n(t,t_0)C'(t_0,t-1)V_n'(t),
\label{eq:GEVPprime}
\end{align}
where $V_n'(t)$ satisfies
\begin{equation}
V_n(t) = T(t-1)V_n'(t).
\label{eq:Vprime}
\end{equation}
Again, the second argument $t_0$ of $V_n$ and $V_n'$ is omitted as noted above.  While $V_n'$ should have another argument $t-1$ that specifies the argument of $T$, we also drop it to avoid the confusion with the omitted second argument $t_0$.
If we can obtain the correct order of $V_n'(t)$ when solving the modified GEVP \eqref{eq:GEVPprime}, we can also obtain the corresponding GEVP eigenvectors $V_n(t)$ with the original basis through Eq.~\eqref{eq:Vprime}.
Since $T(t-1)$ decomposes the correlators into contributions from the lowest few states and should be $t$-independent at sufficiently large $t$ where the contamination from the $N+1$-th and higher states is negligible.
Since $T(t-1)$ and $T(t)$ therefore should be mostly identical at such large time separations, $V_n'(t)$ should be close to a unit vector.  The correct ordering can be achieved by choosing $V_n'(t)$ as the GEVP eigenvector of Eq.~\eqref{eq:GEVPprime} whose $n$-th element is most significantly non-zero.

In this work, we order eigenvalues and eigenvectors of GEVP based on the order of eigenvalues at $t_0 = 1,2$ and then switch to the recursive procedure described above.
The reason why we switch the procedure at a certain value of $t_0$ rather than $t$ is because the GEVP eigenvectors $V_n(t,t_0)$ receive the contamination from the higher states by $O(e^{-(E_{N+1}-E_n)t_0})$~\cite{Blossier:2009kd}.

It is guaranteed that the GEVP eigenvalues and eigenvectors in Eq.~\eqref{eq:GEVP} are real when $C(t_0)$ is a positive-definite real symmetric matrix and some software functions to solve GEVP have this assumption.
While the correlator matrix at large time separations with limited statistics have zero-consistent eigenvalues and it is inevitable that some jackknife samples of correlator matrix have negative eigenvalues, that may not necessarily mean we cannot solve GEVP or obtain any information of two-pion signals at those time separations.
The positivity and real symmetry are sufficient but not necessary to give us real eigenvalues and eigenvectors of GEVP.
As shown in Section~\ref{sec:energies}, we find that GEVP with $C(t_0)$ including a negative eigenvalue can still give us good signals of a few lowest energy states with the ordering strategy explained above.

\subsection{Re-based GEVP}

The GEVP method provides a decomposition based on the $N$ lightest states.  It assumes sufficiently large time separations, so the contamination of higher excited states ($n\ge N+1$) can be ignored.  On the other hand, the GEVP becomes increasingly difficult with increasing time separation due to the exponentially deteriorating signal-to-noise of the correlator matrix. Ideally, significant signal loss occurs after reaching a plateau in time of the effective energies.  For small $N$, we expect the plateau to start at larger time separations, where the signal-to-noise ratio is already poor.  For large $N$, plateaus for the various energies move to earlier times, but the larger statistical errors on the higher energy states in the GEVP might spoil the signal of the lower states.  Thus for any choice of $N$ (or operator set) there is a chance that the signal loss occurs before a clear plateau is observed.

To address the problem we propose a modified version of the GEVP which we call the {\it re-based} GEVP (RGEVP) whereby we choose a new reduced basis of operators at increasing times as explained below.  The idea originates from the fact that the number of resolvable states which dictates the appropriate size of GEVP decreases with increasing time separation  and that the GEVP at short time separations, even before reaching a plateau, provides a set of nearly diagonalized operators.

In this work we consider the RGEVP with fixed $t-t_0 \equiv \Delta_t$, and we omit $t_0$ from the arguments of $V_n$ in this subsection.  The recipe for the RGEVP is as follows:
\begin{enumerate}
    \item Choose the GEVP size $N_t$ at each time slice starting at $\Delta_t+1$
    \begin{eqnarray}
        N_{\Delta_t+1} \ge N_{\Delta_t+2} \ge \ldots.
    \end{eqnarray}
    \item Solve the $N_{\Delta_t+1}\times N_{\Delta_t+1}$ GEVP at the beginning time slice $t=\Delta_t+1$ and order the eigenvectors $V_n(\Delta_t+1)$ ensuring the descending order of the corresponding eigenvalues.
    \item For $t\ge\Delta_t+2$, suppose we have obtained the eigenvectors $V_n(t-1)$ at $t-1$ for $n=1,2,\ldots,N_{t-1}$ and calculate the $N_t\times N_t$ re-based correlator matrices by Eqs.~\eqref{eq:re-basedCt} and \eqref{eq:re-basedCt0} with a modified $N_{\Delta_t+1}\times N_t$ re-basing matrix
    \begin{eqnarray}
        T(t-1) = (V_1(t-1)\,\,V_2(t-1)\,\ldots V_{N_t}(t-1)).
    \end{eqnarray}
    \label{step:rebase}
    \item Solve the GEVP \eqref{eq:GEVPprime} with the re-based correlators in Eqs.~\eqref{eq:re-basedCt} and \eqref{eq:re-basedCt0} and obtain the eigenvectors $V_n(t)$ by Eq.~\eqref{eq:Vprime} for $n = 1,2,\ldots,N_t$ using the ordering procedure described in the previous subsection.
    \label{step:solveGEVP}
\end{enumerate}
Note that, in step~\ref{step:rebase}, the number $N_t$ of columns of $T(t-1)$ despite $N_{t-1}(\ge N_t)$ eigenvectors $V_n(t-1)$ obtained at time slice $t-1$ plays a role in reducing the size of the GEVP at time slice $t$ when $N_{t-1}>N_t$.
While, repeating steps \ref{step:rebase} and \ref{step:solveGEVP} and applying the results to Eq.~\eqref{eq:def_Eeff} give us a new series of effective two-pion energies, they are identical to the normal GEVP results up to the time separation where the first actual re-basing, {\it i.e.} reducing the size of GEVP, is performed.
On the other hand, it is interesting to investigate how the GEVP with the reduced new basis behaves at smaller time slices.  Therefore we perform the following additional steps:
\begin{enumerate}
\setcounter{enumi}{4}
\item Repeat steps \ref{step:rebase} and \ref{step:solveGEVP} to obtain $T(\tilde t)$ with a chosen time $\tilde t$ where the final re-basing is performed.
\item Perform the GEVP analysis with the re-based correlator matrix $T(\tilde t)^\dag C(t)T(\tilde t)$ for all available time slices.
\end{enumerate}
These steps can be parametrized by the time slices where the GEVP size is reduced and the new re-based GEVP sizes at these time slices.  We present the RGEVP results with these parameters in Section~\ref{sec:results}.
In this work, we implement the re-basing with the central values of $T(t-1)$ for all jackknife samples to maintain the configuration independence of the new operator basis.

\subsection{Phase shifts}
\label{sec:phase shifts theory}

The L\"uscher method~\cite{Luscher:1990ux} allows us to extract the scattering phase shifts from finite-volume energies on the lattice. The interaction region is supposed to be confined to a volume well contained inside a box of larger volume.  Outside this region the solution of the wave equation corresponds to free (non-interacting) particles, and the boundary conditions of the box impose a quantization condition on the complete solution which necessarily relates the phase shifts to finite-volume energies.

While the extension to moving frames is straightforward \cite{Rummukainen:1995vs}, we limit our discussion to the case of the rest frame, where one obtains the two-pion $s$-wave scattering phase shift $\delta(E_{\pi\pi})$ corresponding to a given two-pion energy as follows:
  \begin{align}
    k&=\sqrt{\frac{E_{\pi\pi}^2}{4}-m_{\pi}^2},

    \label{eq:disp_rel}
    \\q&= \frac{kL}{2\pi},
    \\ \tan~\phi(q)&=\frac{\pi^{3/2}q}{Z_{00}(1;q^2)},
  \\\delta(E_{\pi\pi}) &= -\phi(q)+\pi n,~n\in \mathbb{Z},\label{eq:phase shift}
  \end{align}
where $E_{\pi\pi}$ is the energy of the two-pion state in a finite box of size $L^3$, $m_\pi$ is the pion mass and $Z_{00}(1,q^2)$ is the L\"uscher zeta function, which we compute via an efficient numerical implementation given in Ref.~\cite{Yamazaki:2004qb}.
Eq.~(\ref{eq:phase shift}) is used to compute all $I=0$ and 2 phase shifts.

The method outlined here for obtaining the phase shifts is strictly valid in the limited region $2m_\pi\le E_{\pi\pi}\le4m_\pi$ and up to neglected higher partial waves. In the present work we apply this method also to energies above the $4\pi$ inelastic threshold, and neglect these sources of systematic errors.  Numerical results presented in Section~\ref{sec:phase shifts res} suggest these systematic errors may not be large within the energy range considered here.

\subsection{The dispersion relation method}

Up to small finite lattice spacing effects, L\"uscher's method, explained above, gives an accurate prescription for obtaining the phase shift.
With finite lattice spacing, Eq.~\eqref{eq:disp_rel} needs modification since the appropriate dispersion relation for finite lattice spacing depends on the type of lattice fermion. While simulations at multiple lattice spacings enable the removal of these effects, we remove some of them for each lattice spacing separately~\cite{RBC:2021acc} .

The method is based on the cancellation of artifacts between interacting and non-interacting two-pion energies.
The non-interacting two-pion energies, $E_n^{0}$, are determined from a product of two expectation values of single-pion correlators, $C^{0}(t)$, which is analogous to the D diagram in the interacting case.
Since the correlation functions involving the $\sigma$ operator do not contain the D diagram, they should be treated separately.
We first explain the method for $I=2$, where the $\sigma$ operator is absent, and then explain its generalization to the case including the $\sigma$ operator.

The non-interacting correlator matrix $C^{0}(t)$ is diagonal with analogous two-pion effective energies, $E_1^{0,eff}(t,t_0)<E_2^{0,eff}(t,t_0)<\ldots<E_N^{0,eff}(t,t_0)$. The energy shift of two-pion states due to interactions is
\begin{equation}
    \Delta E^{eff}_n(t,t_0) = E^{eff}_n(t,t_0) - E^{0,eff}_n(t,t_0),
    \label{eq:energy_shift}
\end{equation}
where single-pion discretization errors largely cancel.
Thus we obtain an improved two-pion energy
\begin{equation}
    E_n^{eff\,\prime}(t,t_0)
    = E^{0,disp}_n + \Delta E^{eff}_n(t,t_0),
    \label{eq:newEeff}
\end{equation}
where we define
\begin{equation}
    E^{0,disp}_n = 2\sqrt{m_\pi^2 + \frac{1}{4}|\vec p_n|^2},
    \label{eq:disp_nonintEpipi}
\end{equation}
with a non-interacting pion momentum in the finite box $\vec p_n = (0,0,0), (0,0,2\pi/L), (0,2\pi/L,2\pi/L),\ldots$.

In addition to reducing the scaling violation in the dispersion relation, this method determines two-pion energies and phase shifts with other improvements as well.
The first term on the right hand side of Eq.~\eqref{eq:newEeff}, which is defined in Eq.~\eqref{eq:disp_nonintEpipi}, is as statistically precise as $m_\pi$.
On the other hand, the second term on the right hand side of Eq.~\eqref{eq:newEeff}, which is given in Eq.~\eqref{eq:energy_shift}, is also expected to be more precise than $E_n^{eff}(t,t_0)$ because of the correlation between the first and second terms on the right hand side of Eq.~\eqref{eq:energy_shift}.
Furthermore, the energy difference in Eq.~\eqref{eq:energy_shift} may also remove excited state effects related to single pions and allow effective energies $E_n^{eff\,\prime}$ to plateau sooner.
We find these to be the case, especially for the $I=2$ channel, as discussed in Section~\ref{sec:results}.

For the $I=0$ channel, the statistical errors are dominated by the disconnected diagram which cannot be improved by this method, though we still find some improvement for the combination of diagrams by applying the following procedure with the $\sigma$ operator.
The $I=0$ channel is more complicated not only because of the inclusion of the $\sigma$ operator but also because the interaction between the two pions makes the finite-volume two-pion energies quite unlike the energies of non-interacting two pions, {\it i.e.} $|\Delta E_n^{eff}(t,t_0)|$ for $I=0$ is much larger than that for $I=2$.  Therefore it is less meaningful to identify the one-to-one correspondence between the interacting and non-interacting two-pion energies as in Eq.~\eqref{eq:energy_shift}.
In this work instead of matching the state label $n$ of the interacting and non-interacting two-pion energies, the first and second terms in Eq.~\eqref{eq:energy_shift}, we choose the non-interacting two-pion energy as the one closest to the interacting $n$-th state energy $E_n^{eff}(t,t_0)$ for the procedure explained above.

\section{Ensemble Details and Computational Setup}
\label{sec:lat details}

Our computations are carried out on two ensembles of 2+1 flavors of M\"obius domain wall fermions (MDWF) with physical masses generated by the RBC/UKQCD collaborations~\cite{RBC:2014ntl,Tu:2020vpn}. Both use the Iwasaki-DSDR gauge action~\cite{Renfrew:2008zfx} and correspond to inverse lattice spacings of about 1.0 and 1.4 GeV, respectively, with similar physical volumes ($L\sim5$ fm). The parameters of each ensemble are listed in Tab.~\ref{tab:ensembles}.

\begin{table}[h]
    \centering
    \begin{tabular}{c|c|c|c|c|c|c}
    \hline
         $m_\pi$& lattice &  & $a^{-1}$ &$L$ & trajectories&  \\
        (MeV) & size & $L_s$ &  (GeV) & (fm) & (MD time units) & configs  \\
        \hline
        142.6(3) & $24^3\times64$ & 24 & 1.023(2) & 4.67 & 250-3860 & 258\\
        143.6(9) & $32^3\times64$ & 12 & 1.378(5) & 4.58 & 200-1320 & 107\\
        \hline
    \end{tabular}
    \caption{Ensemble parameters. 2+1 flavors of M\"obius domain wall fermions, generated by the RBC/UKQCD collaborations~\cite{RBC:2014ntl,Tu:2020vpn}. Trajectories used for measurements are separated by 10 or 20 Monte Carlo time units. The last column refers to the number of configurations in each ensemble used for measurements.}
    \label{tab:ensembles}
\end{table}

Correlation functions are computed in an all-to-all propagator (A2A)~\cite{Foley:2005ac} framework using 2,000 low-modes of the preconditioned, squared Dirac operator and spin-color-time diluted random source propagators for the high modes.

We employ all-mode averaging (AMA)~\cite{Blum:2012uh,Shintani:2014vja} to save the computational cost for the conjugate gradient (CG).  While the traditional AMA is to perform fewer exact measurements ($e.g.$ with fewer source locations) for all the configurations with which sloppy measurements are performed, we rather reduce the number of configurations for the exact measurements keeping the A2A procedure as it is.  We first perform both exact and sloppy measurements with $N_{exact}$ configurations and create corresponding jackknife samples of the difference between exact and sloppy correlators, which would correct the bias due to the sloppy CG.  In the case of bin size 1,
\begin{equation}
\Delta C^{(l)}(t) = \frac{1}{N_{exact}-1}\sum_{k\neq l} C_{exact}^{{\cal E}_k}(t)-C_{sloppy}^{{\cal E}_k}(t),
\end{equation}
where $\cal E$ is the list of configurations with which both exact and sloppy measurements are performed and $C_{exact/sloppy}^{{\cal E}_k}(t)$ stands for an exact/sloppy sample of the correlator matrix calculated with a configuration ${\cal E}_k$.
The average of the difference is defined as
\begin{equation}
\Delta C(t) = \frac{1}{N_{exact}}\sum_{l=1}^{N_{exact}}\Delta C^{(l)}(t).
\end{equation}
We also perform sloppy measurements with $N_{sloppy}$ additional configurations in the list $\cal S$ and create corresponding jackknife samples.  If we set the bin size to 1, they are written as
\begin{equation}
C_{sloppy}^{(l)}(t) = \frac{1}{N_{sloppy}-1}\sum_{k\neq l}C_{sloppy}^{{\cal S}_k}(t),
\end{equation}
and the average is given as
\begin{equation}
C_{sloppy}(t) = \frac{1}{N_{sloppy}}\sum_{l=1}^{N_{sloppy}}C_{sloppy}^{(l)}(t).
\end{equation}
With these jackknife samples and averages, we define super-jackknife samples as
\begin{equation}
C_{\rm AMA}^{(l)}(t) = \Bigg\{\begin{array}{ll}
C_{sloppy}(t) + \Delta C^{(l)}(t) & (1\le l\le N_{exact})
\\[2mm]
C_{sloppy}^{(l-N_{exact})}(t) + \Delta C(t) & (l> N_{exact})
\end{array}.
\end{equation}

In this work, the sloppy (high-mode) propagators are computed with 400 and 330 iterations of CG and the exact propagators are computed to the CG residual of $10^{-8}$ with 14 and 17 configurations for the $24^3$ and $32^3$ lattices, respectively.  For the sloppy part of the measurements on the $24^3$, 1.023~GeV, ensemble we employ the zM\"obius approximation~\cite{Abramczyk:2017oxr} of the M\"obius Dirac operator, reducing the size of the fifth dimension by a factor of two.

As described in Section~\ref{sec:operator def}, the two-pion operators are defined as a product of two single-pion operators separated by a parameter $\Delta$.  We choose $\Delta = 3$ and 4 for the $24^3$ and $32^3$ lattices, respectively, so the separation in physical units is about the same.
The smearing radius $r$ of the single-pion and the sigma operators is set to 1.5 and 2.0 in lattice units for the $24^3$ and $32^3$ lattices, respectively.
We average the correlation functions over time location of the source operator $t_{src}$.
The disconnected diagrams are computed at $t_{src} = 0,1,\ldots,63$, while the connected diagrams are computed at $t_{src} = 0, 8, \dots, 56$ for the $24^3$ lattice and $t_{src} = 0, 10, \dots, 50$ for the $32^3$ lattice.

\section{results}
\label{sec:results}

\subsection{Energies}
\label{sec:energies}

Interacting two-pion energies are tabulated and plotted for several values of Euclidean times $t_0$, $t-t_0$, $\delta_t$, and various types ($N\times N$ or RGEVP) of GEVP, with and without the dispersion relation method, in Appendix~\ref{sec:2pi energies}. The parameters for re-basing are found in captions of respective figures and tables throughout this subsection and tables in Appendix~\ref{sec:2pi energies}.  Some general patterns are apparent. For short times the statistical errors are sub-percent, even down to the per mille level in some cases. The effect of increasing $t-t_0$ is small and likewise for  $\delta_t$. The dispersion relation (DR) method greatly enhances statistical precision, especially for the $I=2$ ground state but gains are seen generally. It also lessens (single-pion) excited state contamination. In the next two subsections we discuss the results in detail, including comparisons in a range of $t$ and for both lattice spacings to assess residual excited state contamination and lattice artifacts, respectively.

Our guiding principal throughout this analysis is to stick to as short times as possible, where statistics are better and ATW effects are smaller, taking advantage of the GEVP and the DR method that reduce excited state contamination.
In this section, the state label $n$ begins with 0 so that the ground state is labeled with $n=0$, the first excited state with $n=1$ and so on.

\subsubsection{I=2}

We begin with the ground state for the $24^3$ ensemble. Fig.~\ref{fig:I=2 n=0 energies} shows the ground state energy for several GEVP types and representative values of $\delta_t=2,5,8$, $t-t_0=1,2,3,4$, and $t=4,6,8$. Similar patterns of behavior emerge for the various GEVP types. Without the DR method, increasing either $\delta_t$ or $t-t_0$ tends to decrease the energy, suggesting smaller excited state effects (the statistical errors are relatively large, so we do not draw a strong conclusion). For larger times the effect is smaller. For fixed $t$ we often observe smaller statistical errors as $\delta_t$ increases, though after $\delta_t=5$ the improvement is slight. Increasing $t-t_0$ has little effect. The DR method, on the other hand, shows a dramatic reduction in statistical error but little change after that for the other variables. In either case, there is little dependence on the GEVP type. This is because, in Eq.~\eqref{eq:newEeff}, $E_0^{0,disp} = 2m_\pi$ dominates the statistical error and $\Delta E_0^{eff}(t,t_0)$, which is dependent on GEVP type, is much more precise for the ground state in our measurements.

Fig.~\ref{fig:I=2 4x4 eff energy n=0} shows the effective ground state energy for the $4\times4$ GEVP (our largest basis for $I=2$) with $t-t_0=1$ and matrix subtractions in the range $4\le \delta_t\le11$ with (lower panel) and without (upper panel) the DR method. There are several interesting features. First a clear and stable plateau sets in between $t=4$ and 5 for all non-zero $\delta_t$. Both statistical uncertainties and excited state contamination are significantly reduced by the DR method.
In the upper panel (no matrix subtraction) there is a systematic downward shift for each time slice which grows with increasing $t$. This is a clear indication of the ATW effect, which is eliminated by the matrix subtraction. Perhaps somewhat surprisingly the shift is also eliminated by the DR method. This can be understood as follows. A very similar ATW effect occurs in the single pion case, so the observed cancellation implies that in the interacting case, the two pions do not interact very much. We also observe a small reduction in statistical error as $\delta_t$ increases in the absence of the dispersion method while there is no difference when it is used as the dominant error on the ground state energy after DR is from the error on $m_\pi$. Since contractions were only computed for $t\le16$ to reduce computational cost, the maximum value of $t$ for the effective energy in each case is $t_{max,\delta_t}=15-\delta_t$, as seen in the figure (15 not 16 appears because the effective energy depends on GEVP eigenvalues at $t$ and $t+1$).
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24/efmI2_4x4_n0.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI2_4x4_n0.pdf}
    \caption{The $I=2$ effective ground state energy with the dispersion relation method (lower) and without (upper). $24^3$, $4\times4$ GEVP, $t-t_0=1$. Without matrix subtraction there is an evident downward shift in the energy, an indication of the ATW effect, for the non-DR result.}
    \label{fig:I=2 4x4 eff energy n=0}
\end{figure}

In Tab.~\ref{tab:fitres_Epipi_I2_n0_L24} we tabulate the fit results for the ground state energy using the dispersion relation method for different fit ranges and GEVP types and display them in Fig.~\ref{fig:I=2 fits} (upper panel). In this work all fits are correlated fit to a constant ($t$-independent) parameter, performed separately for each individual effective energy.  The stability of the fits is quite robust for all values of $t$ and GEVP types considered in the figure.

The behavior of the first excited state is similar to the ground state.
In Fig.~\ref{fig:I=2 4x4 eff energy n=1} we again observe a stable plateau beginning with $t=5$; however an ATW effect cannot be seen. But we do see excited state effects, especially without the DR method, and statistical errors are reduced significantly by the DR method. There is little dependence on $\delta_t$.
The small deviation visible at $t_{min}=4$ in the fit results in Tab.~\ref{tab:fitres_Epipi_I2_n1_L24} and Fig.~\ref{fig:I=2 fits} (lower panel) may be associated with the slight blip observed in the effective energy at $t=4$ in Fig.~\ref{fig:I=2 4x4 eff energy n=1}. The $2\times2$ GEVP may be systematically high, especially for $t=4$.

The second and third excited state energies are shown in Fig.~\ref{fig:I=2 4x4 eff energy n=2,3} and listed for various fit ranges and GEVP types in Tabs.~\ref{tab:fitres_Epipi_I2_n2_L24} and \ref{tab:fitres_Epipi_I2_n3_L24} (see also Fig.~\ref{fig:I=2 2-3 fits}). In the upper panel the energy fluctuates down at $t=4$ which leads to an elevated $\chi^2$ in the fit. The same happens at $t=6$ for the third excited state. Again, there is little dependence on $\delta_t$ or GEVP type, except that the $3\times3$ GEVP is a bit high for small $t_{min}$ in the fit (Fig.~\ref{fig:I=2 2-3 fits}).With four interpolating operators for $I=2$, the third excited state is as far as we can go.

\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24/efmI2_4x4_n1.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI2_4x4_n1.pdf}
    \caption{The $I=2$ effective first excited state energy with (lower) and without (upper) the dispersion relation method. $24^3$, $4\times4$ GEVP, $t-t_0=1$. The ATW effect is not observed, $c.f.$ Fig.~\ref{fig:I=2 4x4 eff energy n=0}.}
    \label{fig:I=2 4x4 eff energy n=1}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI2_4x4_n2.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI2_4x4_n3.pdf}
    \caption{The $I=2$ effective second (upper) and third (lower) excited state energies with the dispersion relation method. $24^3$, $4\times4$ GEVP, $t-t_0=1$. The ATW effect is not observed, $c.f.$ Fig.~\ref{fig:I=2 4x4 eff energy n=0}.}
    \label{fig:I=2 4x4 eff energy n=2,3}
\end{figure}

\input{tabs/tab_fitres_pipi2.L24.n0.bs1}
\input{tabs/tab_fitres_pipi2.L24.n1.bs1}
\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L24.n0.bs1.pdf}
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L24.n1.bs1.pdf}
    \caption{$I=2$ $\pi\pi$ ground (upper) and first excited (lower) state energies obtained from fits to a constant for various fit ranges, $t_{min}-9,10$, and GEVP types. $24^3$ ensemble.}
    \label{fig:I=2 fits}
\end{figure}
\input{tabs/tab_fitres_pipi2.L24.n2.bs1}
\input{tabs/tab_fitres_pipi2.L24.n3.bs1}
\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L24.n2.bs1.pdf}
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L24.n3.bs1.pdf}
    \caption{$I=2$ $\pi\pi$ second (upper) and third excited (lower) state energies obtained from fits to a constant for various fit ranges, $t_{min}-9,10$, and GEVP types. $24^3$ ensemble.}
    \label{fig:I=2 2-3 fits}
\end{figure}

The effective ground state energy computed on the $32^3$ lattice is shown in Fig.~\ref{fig:I=2 4x4 eff energy 32c n=0}. The ATW effect is even more pronounced, as expected, since the physical time extent is smaller compared to the $24^3$ lattice, and it is significantly reduced again, but not completely eliminated, by the DR method. The lower panel indicates that the effective energy even with the matrix subtraction might decrease with increasing time from $t=10$ and it could mean there are also higher-order ATW effects that cannot be removed by the matrix subtraction in Eq.~\eqref{eq:matdt_subt}.  Since this tendency is not statistically significant and it is still possible that this is due to statistical fluctuation, we will conclude this point after increasing statistics but do not use data point at $t\ge10$ for our fits.  Again a large reduction in the statistical error and excited state contamination occurs with the dispersion relation method where a plateau emerges beginning at $t=4$. Fig.~\ref{fig:I=2 L32 fits} shows fits for several ranges and GEVP types. Results are tabulated in Tab.~\ref{tab:fitres_Epipi_I2_n0_L32}. There is little dependence on $t_{min}$ or GEVP type.
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs32/efmI2_4x4_n0.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI2_4x4_n0.pdf}
    \caption{The $I=2$ effective ground state energy with (lower) and without (upper) the dispersion relation method. $32^3$, $4\times4$ GEVP, $t-t_0=1$. Their is a pronounced ATW effect without matrix subtraction (upper panel). The dispersion relation method reduces, but does not entirely eliminate the effect (lower panel).}
    \label{fig:I=2 4x4 eff energy 32c n=0}
\end{figure}
\input{tabs/tab_fitres_pipi2.L32.n0.bs1}
\input{tabs/tab_fitres_pipi2.L32.n1.bs1}
\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L32.n0.bs1.pdf}
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L32.n1.bs1.pdf}
    \caption{$I=2$ $\pi\pi$ ground (upper) and first excited (lower) state energies obtained from fits to a constant for various fit ranges, $t_{min}$ to 12,13, and GEVP types. $32^3$ ensemble.}
    \label{fig:I=2 L32 fits}
\end{figure}

The situation is similar for the first excited state energy (see Fig.~\ref{fig:I=2 4x4 eff energy 32c n=1}) except the ATW effect is not detectable, the same as we saw for the $24^3$ lattice. The plateau begins at $t=4$ or 5 with the DR method. The DR method removes a bit of jitter as well, which is probably due to poor statistics.

The energies for higher excited states are shown in Fig.~\ref{fig:I=2 4x4 eff energy 32c n=2,3}. Fits are summarized in Tabs.~\ref{tab:fitres_Epipi_I2_n2_L32} and~\ref{tab:fitres_Epipi_I2_n3_L32} and shown in Fig.~\ref{fig:I=2 2-3 L32 fits}.
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs32/efmI2_4x4_n1.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI2_4x4_n1.pdf}
    \caption{The $I=2$ effective first excited state energy with (lower) and without (upper) the dispersion relation method. $32^3$, $4\times4$ GEVP, $t-t_0=1$.}
    \label{fig:I=2 4x4 eff energy 32c n=1}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI2_4x4_n2.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI2_4x4_n3.pdf}
    \caption{The $I=2$ effective second (upper) and third (lower) excited state energies with the dispersion relation method. $32^3$, $4\times4$ GEVP, $t-t_0=1$.}
    \label{fig:I=2 4x4 eff energy 32c n=2,3}
\end{figure}
\input{tabs/tab_fitres_pipi2.L32.n2.bs1}
\input{tabs/tab_fitres_pipi2.L32.n3.bs1}
\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L32.n2.bs1.pdf}
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi2.L32.n3.bs1.pdf}
    \caption{$I=2$ $\pi\pi$ second (upper) and third excited (lower) state energies obtained from fits to a constant for various fit ranges for the $4\times4$ GEVP. $32^3$ ensemble.}
    \label{fig:I=2 2-3 L32 fits}
\end{figure}

\subsubsection{I=0}

This case is statistically noisier than $I=2$ due to the disconnected diagrams of the correlation function as well as to the coupling with the vacuum state. In Fig.~\ref{fig:I=0 n=0 energies} the effective energies for the ground state are shown for several GEVP types computed on the $24^3$ ensemble. Values of $t-t_0$ range from 1 to 4 with $3\le t \le 7$ and $\delta_t=2,5,8$. Like $I=2$, there is little dependence on GEVP type. The dispersion relation method reduces excited state contamination and statistical errors as does increasing $\delta_t$. The effect of $t-t_0$ is less clear, though it appears it may also reduce excited state effects.

Fig.~\ref{fig:I=0 5x5 eff energy n=0} shows the effective ground state energy for the $5\times5$ GEVP with $t-t_0=1$ and matrix subtractions in the range $7\le \delta_t\le11$, with and without the dispersion relation method and for the RGEVP. A small ATW effect may be visible in the upper panel when no matrix subtraction is performed, and it is largely absent in the middle and lower panels, showing again that the dispersion relation method largely eliminates it. We also observe that the RGEVP makes a moderate improvement on the statistical errors for larger times (lower panel).
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24/efmI0_5x5_n0.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI0_5x5_n0.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR-RGEVP/efmI0_n0_full_T4S012.pdf}
    \caption{The $I=0$ effective ground state energy. $24^3$, $5\times5$ GEVP, $t-t_0=1$ with (middle, lower) and without (upper) the dispersion relation method. The $3\times3$ RGEVP result is shown in the lower panel, starting from the middle panel. The re-basing uses GEVP eigenvectors from time slice $t_0=4$.}
    \label{fig:I=0 5x5 eff energy n=0}
\end{figure}

Fit results are summarized in Tab.~\ref{tab:fitres_Epipi_I0_n0_L24} and Fig.~\ref{fig:I=0 fits} (upper panel). There appears to be a small systematic shift with the minimum time separation in the fit, $t_{\rm min}$, but it is within the statistical uncertainty.
\input{tabs/tab_fitres_pipi0.L24.n0.bs1}
\begin{figure}[hbt]
    \includegraphics[width=0.975\linewidth]{figs/v2/fig_fitres_pipi0.L24.n0.bs1.pdf}
    \includegraphics[width=0.975\linewidth]{figs/v2/fig_fitres_pipi0.L24.n1.bs1.pdf}
    \includegraphics[width=0.975\linewidth]{figs/v2/fig_fitres_pipi0.L24.n2.bs1.pdf}
    \caption{$I=0$ $\pi\pi$ ground (upper), first (middle), and second (lower) excited state energies obtained from fits to a constant for various GEVP types and fit ranges, $t_{min}$ to 6 or 7. $24^3$ ensemble.}
    \label{fig:I=0 fits}
\end{figure}

The first excited state energies are summarized in Tab.~\ref{tab:efm_I0_n1_L24} and Fig.~\ref{fig:I=0 n=1 energies} for a wide range of parameters. In Fig.~\ref{fig:I=0 5x5 eff energy n=1} we show the first excited state energy for the $5\times5$ GEVP and $t-t_0=1$. A plateau begins at $t=3$ or 4, and an interesting systematic begins to emerge as well which is even more pronounced for the second excited state, as we see shortly. A careful inspection of Fig.~\ref{fig:I=0 5x5 eff energy n=1} shows a systematic downward drift of the energy for each $t_{max,\delta_t}$. It is most visible between $t=6$ and 7 for $\delta_t=8$, but a smaller shift appears at smaller times for large values of $\delta_t$. The effect is absent for the RGEVP (lower panel), and leads us to believe it is due to a breakdown of the GEVP at large $t$ due to large statistical errors in the correlation matrix.
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24/efmI0_5x5_n1.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI0_5x5_n1.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR-RGEVP/efmI0_n1_full_T3S012.pdf}
    \caption{The $I=0$ effective first excited state energy with (lower) and without (upper) the dispersion relation method. $24^3$, $5\times5$ GEVP, $t-t_0=1$.}
    \label{fig:I=0 5x5 eff energy n=1}
\end{figure}

Fit results are given in Tab.~\ref{tab:fitres_Epipi_I0_n1_L24} and displayed in Fig.~\ref{fig:I=0 fits} (middle panel). There is no detectable difference as $t_{min}$ or the GEVP type varies.
\input{tabs/tab_fitres_pipi0.L24.n1.bs1}

The second excited state is less well resolved statistically, so we restrict our focus to relatively small times (see Fig.~\ref{fig:I=0 5x5 eff energy n=2}). There is little difference without or with the DR method (upper and middle panels); however the downward trend observed for the first excited state is even more visible here and is also removed by the RGEVP although the statistical errors increase significantly. Fit results are listed in Tab.~\ref{tab:fitres_Epipi_I0_n2_L24} and plotted in Fig.~\ref{fig:I=0 fits}.
\input{tabs/tab_fitres_pipi0.L24.n2.bs1}

\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24/efmI0_5x5_n2.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI0_5x5_n2.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs24DR-RGEVP/efmI0_n2_full_T1S012.pdf}
    \caption{The $I=0$ effective second excited state energy with (middle, lower) and without (upper) the dispersion relation method. $24^3$, $5\times5$ (middle, upper) and $3\times3$ (lower) GEVP, $t-t_0=1$. Note the significant jump from $t=6$ to 7 for $\delta_t=8$ in upper two panels. Smaller shifts happen at earlier times.}
    \label{fig:I=0 5x5 eff energy n=2}
\end{figure}

Another dramatic difference happens when we reduce the GEVP to $3\times3$ ($4\times4$ is essentially the same as $5\times5$) without re-basing (see Fig~\ref{fig:I=0 3x3 eff energy n=2}). The noise increases dramatically because of the significant coupling between the second excited state and the $\pi\pi(011)$ operator which is excluded in the $3\times3$ analysis. This is similar to the effect observed in Refs.~\cite{RBC:2020kdj} and~\cite{RBC:2021acc} where adding a noisy scalar bilinear operator reduced the statistical error of (and excited state contamination in) the ground state, although it is the $\pi\pi(011)$ operator that is not included in the $3\times3$ but in the $4\times4$ and $5\times5$ analyses in this case.

The GEVP eigenvalues are too noisy even for short times to extract meaningful energies for the higher energy states. They will have to wait for better statistics in the future.

\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs24DR/efmI0_3x3_n2.pdf}
    \caption{The $I=0$ effective second excited state energy for the $3\times3$ GEVP, $t-t_0=1$. The noise has dramatically increased for this state compared to Fig.~\ref{fig:I=0 5x5 eff energy n=2}.}
    \label{fig:I=0 3x3 eff energy n=2}
\end{figure}

On the $32^3$ ensemble, the pattern repeats except the effective energies are even noisier as seen in Fig.~\ref{fig:I=0 5x5 eff energy 32c n=0}. The ATW effect is clearly visible in the top panel and absent for the DR method. It is also interesting to note that the statistical error in the latter case is smaller at large time when no matrix subtraction is performed. Fits are summarized in Tab.~\ref{tab:fitres_Epipi_I0_n0_L32} and Fig.~\ref{fig:I=0 32c fits}.
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs32/efmI0_5x5_n0.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI0_5x5_n0.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR-RGEVP/efmI0_n0_full_T1S0123ave_T2S012ave_T4S01ave.pdf}
    \caption{The $I=0$ effective ground state energy without (upper) and with (middle, lower) the dispersion relation method. $32^3$, $5\times5$ GEVP, $t-t_0=1$. The ATW effect is clear in the top panel. The lower panel shows the $2\times2$ RGEVP result, using $t_0=1$, 2, and 4 to go from $5\times5\to2\times2$ in three steps.}
    \label{fig:I=0 5x5 eff energy 32c n=0}
\end{figure}
\input{tabs/tab_fitres_pipi0.L32.n0.bs1}
\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi0.L32.n0.bs1.pdf}
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi0.L32.n1.bs1.pdf}
    \includegraphics[width=\linewidth]{figs/v2/fig_fitres_pipi0.L32.n2.bs1.pdf}
    \caption{$I=0$ $\pi\pi$ ground (upper), first excited (middle), and second excited (bottom) state energies obtained from fits to a constant for various GEVP types and fit ranges, $t_{min}$ to 8 or 13. $32^3$ ensemble.}
    \label{fig:I=0 32c fits}
\end{figure}

After $t=3$ or 4, the statistical errors get very large for the first excited state even with the DR method (see Fig.~\ref{fig:I=0 5x5 eff energy 32c n=1}).
\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI0_5x5_n1.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR-RGEVP/efmI0_n1_full_T1S0123ave_T2S012ave_T4S01ave.pdf}
    \caption{The $I=0$ effective first excited state energy for the $5\times5$ GEVP (upper) and $2\times2$ RGEVP (lower) with the dispersion relation method. $32^3$, $t-t_0=1$.}
    \label{fig:I=0 5x5 eff energy 32c n=1}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=1.\linewidth]{figs/Effs32DR/efmI0_5x5_n2.pdf}
    \includegraphics[width=1.\linewidth]{figs/Effs32DR-RGEVP/efmI0_n2_full_T1S0123ave_T2S012ave.pdf}
    \caption{The $I=0$ effective second excited state energy for the $5\times5$ GEVP (upper) and $3\times3$ RGEVP (lower) with the dispersion relation method. $32^3$, $t-t_0=1$.}
    \label{fig:I=0 5x5 eff energy 32c n=2}
\end{figure}
However, we observe a large reduction in the statistical error using the RGEVP compared to the standard GEVP in this case. Re-basings to go from  $5\times5$ down to $2\times2$ are implemented at $t_0=1$, $t_0=2$ and $t_0=4$. Presumably the improvement occurs because the overlap of the higher states with the lower states is more and more unresolved with increasing $t$, adding only noise to the GEVP. Fit results are tabulated in Tab.~\ref{tab:fitres_Epipi_I0_n1_L32}.
\input{tabs/tab_fitres_pipi0.L32.n1.bs1}

The second excited state is shown in Fig.~\ref{fig:I=0 5x5 eff energy 32c n=2} where again there is a large reduction in statistical error with the RGEVP. A drift downwards with increasing time is observed for larger values of $\delta_t$ for the $3\times3$ RGEVP, though it is not easy to investigate the systematic error at this point because of the accompanied large statistical errors. Fit results are summarized in Tab.~\ref{tab:fitres_Epipi_I0_n2_L32}.
\input{tabs/tab_fitres_pipi0.L32.n2.bs1}

In preparation for the next subsection, we now summarize in Tab.~\ref{tab: energy phase shift summary} the energies and corresponding time-slice fit ranges and GEVP types that we will use to compute the phase shifts. In all cases except the third excited state for $I=2$, which is only accessible with the $4\times4$ GEVP, we choose the RGEVP since the noise is usually reduced or unchanged from the ordinary GEVP. In addition we believe it is more robust at large times since poorly resolved elements of the correlation function are avoided. As mentioned already, we stick to small time slices where the signal is better resolved, and we use the results of the DR method throughout.  With the DR method we see negligible dependence on $t-t_0$, so we take it to be 1 in all cases.

For the $I=2$ ground state we choose 4--10 ($24^3$) and 4--9 ($32^3$) for the fit range. From Figs.~\ref{fig:I=2 fits}~and~\ref{fig:I=2 L32 fits} one can see that the results are very stable with both the range and GEVP type.

Similarly for the first excited state we choose 5--9 and 4--10 (in the case of $24^3$ a small excited state effect may be visible at $t=4$). In both cases the $2\times2$ GEVP appears to be a little high. For the second excited state we take 3--9 and 3--7. Finally for the third excited state our ranges are 3--5 and 4--7.

For $I=0$ it is not as straightforward to choose central values due to larger statistical errors, especially for the $32^3$ ensemble, and as we have seen already, small times comprise the set of usable time slices.

From Figs.~\ref{fig:I=0 fits}($24^3$) and \ref{fig:I=0 32c fits}($32^3$) one sees small variations with $t_{min}$ for the ground states. We choose fit ranges 4--8 and 5--9, respectively. On the $24^3$ lattice we see no significant variation for the first excited state, so we take 3--6. For $32^3$ the errors are large for all but the RGEVP, and we take 4--8. Finally for the second excited state the fit ranges are chosen to be 3--6 and 4--6. In the latter case we use $\delta_t=3$ since we observed a flatter plateau and smaller statistical errors for this choice. Higher excited states cannot be extracted from our data; improved statistics are needed.

\begin{table*}
\centering
\begin{tabular}{|c|c|c|c|c|c|r|r|r|}
\hline
 $I$& $n$ & GEVP type & fit range & $\delta_t$ & $\chi^2/d.o.f.$& energy & $\delta_0$ (deg) & $m_\pi a_0^I$ \\
\hline
\multicolumn{9}{|c|}{$24^3$ lattice}\\
\hline
2 & 0 & RGEVP & 4--10 & 5[4--8] & 0.6 & 0.28130(34)& $-0.374(13)(4)$ & $-0.0496(11)(5)$ \\
2 & 1 & RGEVP & 5--9 & 5[4--8] & 0.3 & 0.60789(31)& $-12.33(22)(20)$&\\
2 & 2 & RGEVP & 3--9 & 5[4--8] & 1.4 & 0.81753(56)& $-20.18(43)(54)$&\\
2 & 3 & $4\times4$ & 3--5 & 5[4--8] & 0.1 & 0.9658(17)& $-26.5(2.4)(4)$&\\
0 & 0 & RGEVP & 4--8 & 7[5--9] & 1.1 & 0.27069(41)& & 0.2038(70)(160)\\
0 & 1 & RGEVP & 3--6 & 9[6--9] & 0.1 & 0.5308(44)& 45.1(2.9)(2.1)& \\
0 & 2 & RGEVP & 3--5 & 9[7--9] & 0.7& 0.695(12) & 83(11)(15)&\\
\hline
\hline
\multicolumn{9}{|c|}{$32^3$ lattice}\\
\hline
2 & 0 & RGEVP & 4--9 & 8[4--8] & 0.8 & 0.21041(37) & $-0.424(51)(23)$&$-0.0537(42)(22)$\\
2 & 1 & RGEVP & 4--10 & 8[4--8] & 0.3 & 0.45648(36)& $-13.37(32)(11)$&\\
2 & 2 & RGEVP & 3--7 & 8[4--8] & 0.3 & 0.61552(58) & $-22.90(59)(34)$&\\
2 & 3 & $4\times4$ & 4--7 & 8[4--8] & 0.1 & 0.7230(17)& $-24.2(3.2)(2.1)$&\\
0 & 0 & RGEVP & 5--9 & 5[5--8] & 0.4 & 0.20248(56)& & 0.1947(150)(126)\\
0 & 1 & RGEVP & 4--8 & 5[5--8] & 0.6 & 0.4052(41)& 38.6(3.6)(10.1) &\\
0 & 2 & RGEVP & 4--6 & 5[4--6] & 0.4 & 0.5304(91)& 71(12)(16) & \\
\hline
\end{tabular}
\caption{Summary of two-pion energies and phase shifts ($\delta)$. The last column gives the scattering length in units of the pion mass. The first error is statistical, and the second ($\delta$, $m_\pi a_0$) a systematic error, which is estimated by varying fit range as well as $\delta_t$ in the range shown in the square brackets.  See text for more detail}
\label{tab: energy phase shift summary}
\end{table*}
\subsection{Phase shifts}
\label{sec:phase shifts res}

The phase shifts are computed using the fitted energies described in the previous subsections. Specifically, we take fit values corresponding to Eq.~(\ref{eq:newEeff}) and insert them into Eq.~(\ref{eq:phase shift}). For the $I=2$ ground state with pions at rest, the interacting two-pion energy is above the $2\, m_\pi$ threshold since the interaction is repulsive, and a phase shift is readily calculated using Eq.~(\ref{eq:phase shift}). For $I=0$ the ground state is below the threshold, so the phase shift is purely imaginary.
The phase shifts for each state and various fit ranges for the energies are summarized in Tabs.~\ref{tab:pshift_I2_n0_L24}-\ref{tab:pshift_I0_n2_L32}.

\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/phaseshift.pdf}
    \caption{Pion scattering phase shifts for ground and higher excited states. $I=0$ (upper) and $I=2$ (lower). Bars denote statistical and systematic  errors added in quadrature (see Tab.~\ref{tab: energy phase shift summary}. The Roy Equation results~\cite{Ananthanarayan:2000ht} are shown by solid lines and corresponding error bands. Vertical gray lines denote $4 m_\pi$ thresholds for $24^3$ and $32^3$ ensembles, above which the method of determining phase shifts~\cite{Luscher:1990ux} in this work is no longer strictly valid.}
    \label{fig:phase shift}
\end{figure}

Results corresponding to the energies (fit ranges) chosen at the end of the last section are shown in Fig.~\ref{fig:phase shift} and given in Tab.~\ref{tab: energy phase shift summary} for both ensembles. Agreement with the Roy Equation result\footnote{In the figures displaying Roy Equation results we always use the formulae in Ref.~\cite{Ananthanarayan:2000ht}, for technical reasons. In later papers we will switch to the update in ~\cite{Colangelo:2001df}.}~\cite{Ananthanarayan:2000ht,Colangelo:2001df} is observed up to the third (second) excited state for $I=2$(0). Even though, strictly speaking, L\"uscher's method~\cite{Luscher:1990ux} is not valid above the inelastic scattering threshold ($4m_\pi$) which is shown for each ensemble by the vertical gray lines in Fig.~\ref{fig:phase shift}, we see no evidence of a dramatic breakdown above this threshold (or higher pion multiplicities, $K\bar K$, and so on). It will be interesting to see in future works with more precision if a signal of such a breakdown emerges.

In addition to statistical errors, we also estimate a systematic error due to choice of fit range for the energy and $\delta_t$ in the matrix subtraction to remove the ATW effect. The error is estimated by comparing phase shifts for various fit ranges and $\delta_t$ values and taking half the difference between minimum and maximum central values for a given state. For the comparison of the systematic error, we use the fit ranges shown in Tabs.~\ref{tab:pshift_I2_n0_L24}--\ref{tab:pshift_I0_n2_L32} and a range of $\delta_t$ values shown in the square bracket in Tab.~\ref{tab: energy phase shift summary}.  We remove the fit results that have more than twice as large statistical error as the fit result shown in Tab.~\ref{tab: energy phase shift summary} from the estimation of the systematic error.  In all cases we do not vary the GEVP type and choose the RGEVP except for the $I=2$ third excited state where we need to use the full-size $4\times4$ GEVP. This is because only the RGEVP appears to give reasonable results for the $I=0$ excited states on the $32^3$ lattice, while the other results are less dependent on GEVP type than fit range.  The systematic errors are small for the ground state (less than the statistical error) and become comparable, and even larger than, the statistical error with increasing energy. Since the $I=0$ energies are noisier than their $I=2$ counterparts, so too are their phase shifts.

The physical energies corresponding to the various states differ slightly on the two ensembles, so to compare the phase shifts at fixed energy, we do a piece-wise linear interpolation at fixed lattice spacing and then extrapolate to the continuum, $a\to0$, at fixed energy. The extrapolation is performed linearly in $a^2$. The results are shown in Fig.~\ref{fig:phase shift 0 and 2}. In each panel the $a\to0$ extrapolation is shown with a green band. The values are compatible with the Roy Equation, albeit within relatively large uncertainties especially for $I=0$. For this study we calculate error bands for the dispersive results using Ref~\cite{Ananthanarayan:2000ht} and we defer a more comprehensive study, including for example the more accurate results of Ref.~\cite{Colangelo:2001df}, to a future publication.

\begin{figure}[hbt]
    \includegraphics[width=\linewidth]{figs/phaseshiftI0.pdf}
    \includegraphics[width=\linewidth]{figs/phaseshiftI2.pdf}
    \caption{Continuum limit of the pion scattering phase shifts for ground and higher excited states (light green bands) for $I=0$ (upper panel) and $I=2$ (lower).  The continuum extrapolation is performed linearly in $a^2$ after linear interpolation in the energy at fixed lattice spacing. The Roy Equation results~\cite{Ananthanarayan:2000ht} are shown by the light red error bands.}
    \label{fig:phase shift 0 and 2}
\end{figure}

\input{tabs/tab_pshift_pipi2.L24.n0.bs1}
\input{tabs/tab_pshift_pipi2.L24.n1.bs1}
\input{tabs/tab_pshift_pipi2.L24.n2.bs1}
\input{tabs/tab_pshift_pipi2.L24.n3.bs1}
\input{tabs/tab_pshift_pipi2.L32.n0.bs1}
\input{tabs/tab_pshift_pipi2.L32.n1.bs1}
\input{tabs/tab_pshift_pipi2.L32.n2.bs1}
\input{tabs/tab_pshift_pipi2.L32.n3.bs1}
\input{tabs/tab_pshift_pipi0.L24.n1.bs1}
\input{tabs/tab_pshift_pipi0.L24.n2.bs1}
\input{tabs/tab_pshift_pipi0.L32.n1.bs1}
\input{tabs/tab_pshift_pipi0.L32.n2.bs1}

Lastly, we discuss the two-pion scattering length.  Since the relation between the two-pion phase shift $\delta_0$ and scattering length $a_0$ is given as\footnote{The subscript refers to the $s$-wave channel.}
\begin{equation}
k\cot{\delta_0(k)} = \frac{1}{a_0} + \frac{1}{2}r_0k^2 +O(k^4),
\label{eq:delta_expansion}
\end{equation}
with the effective range parameter $r_0$,
we can calculate the scattering length by
\begin{equation}
a_0 = \frac{\tan{\delta_0(k)}}{k} + O(k^2),
\label{eq:a0}
\end{equation}
with a sufficiently small value of $k$, which can be obtained with the ground state straightforwardly for the $I=2$ channel.
For the $I=0$ channel, while $k$ is a pure imaginary number for the ground state, we can analytically continue Eq.~\eqref{eq:delta_expansion} and the formulae given in Section~\ref{sec:phase shifts theory} \cite{Fu:2017apw}. With that, we obtain the scattering length as a real number. The values for the scattering lengths for various ranges of fits to the ground state energies are given in Tabs.~\ref{tab:a0_I2_L24}-\ref{tab:a0_I0_L32}. In Tab.~\ref{tab: energy phase shift summary} central values are displayed; systematic errors are computed as before for the phase shifts. A simple linear extrapolation in $a^2$, after combining statistical and systematic errors in quadrature, yields
\begin{eqnarray}
    m_\pi a_0^0&=& 0.184(47),\\
    m_\pi a_0^2 &=& -0.058(11)
\end{eqnarray}
for $I=0$ and 2, respectively.

\input{tabs/tab_a0mpi_pipi2.L24.bs1}
\input{tabs/tab_a0mpi_pipi2.L32.bs1}
\input{tabs/tab_a0mpi_pipi0.L24.bs1}
\input{tabs/tab_a0mpi_pipi0.L32.bs1}

These numbers are consistent with ChPT prediction and earlier lattice calculations \cite{Colangelo:2000jc,FlavourLatticeAveragingGroupFLAG:2021npn} within 1.3$\sigma$ for $I = 2$ and 1$\sigma$ for $I = 0$, while our error is somewhat larger because the earlier works at unphysical pion masses performed a chiral extrapolation with the assumption
\begin{align}
m_\pi a_0^0 &= \frac{7m_\pi^2}{16\pi f_\pi^2}
\left\{
1-\frac{m_\pi^2}{16\pi^2f_\pi^2}
\left[
9\ln\frac{m_\pi^2}{f_\pi^2}-5-l_{\pi\pi}^0
\right]
\right\},
\\
m_\pi a_0^2 &= -\frac{m_\pi^2}{8\pi f_\pi^2}
\left\{
1+\frac{m_\pi^2}{16\pi^2f_\pi^2}
\left[
3\ln\frac{m_\pi^2}{f_\pi^2}-1-l_{\pi\pi}^2
\right]
\right\},
\end{align}
where $f_\pi$ is the pion decay constant and a combination of the low energy constants $l_{\pi\pi}^I$ is the fit parameter for the extrapolation.  Thus the approach taken in the earlier calculations uses scattering lengths computed at unphysical pion masses as the sub-leading contribution to the scattering length at the physical point, with the leading contribution obtained entirely from $m_\pi$ and $f_\pi$ which can be computed very precisely from independent two-point function measurements.  Our results on the other hand are pure lattice results obtained without assuming the leading contribution.

\section{Cost Comparison with G-Parity boundary conditions}
\label{sec:comparison}

One of the main reasons behind the present study is to compare the effectiveness of G-parity boundary conditions against the conventional periodic boundary conditions. In this section we carry out such a comparison.

In Ref.~\cite{RBC:2021acc} measurements were carried out on an ensemble of 741 gauge field configurations with identical parameters to those used in the $32^3$ ensemble except GPBC were used. Nine hundred low-modes and twenty-four (spin-color-flavor diluted) random source fields on each time slice comprised the A2A~\cite{Foley:2005ac} measurement setup. The high-mode part is the same as used in this work, but the low-mode part differs, 900 $vs.$ 2,000 modes in this work and exact eigenvectors $vs.$ approximate coarse-grained eigenvectors based on local-coherence~\cite{Clark:2017wom}.
Another difference is the time-translation interval of connected diagrams.  Most of connected diagrams are calculated with 6 source locations in this work and 8 source locations in the earlier work.
One exception is the connected diagram of $\sigma$--$\sigma$ two-point function, which is computed, again, with 6 source locations in this work but with 64 source locations in the earlier work.  This difference in $\sigma$--$\sigma$ measurement detail would not affect computational cost as it is very cheap, though it might make an impact on the precision.

We will focus on the zero total momentum rest frame. For $I=2$ Tab. X in Ref.~\cite{RBC:2021acc} lists the energy (in lattice units) as 0.41528(46) with $t_0=11$, $t-t_0=1$. In this study the corresponding value is 0.45549(160). Recall that the momentum, in units of $2\pi/L$, of the pions is slightly smaller in the GPBC case, $\sqrt{3\times \left(\frac{1}{2}\right)^2}$ instead of 1 and therefore we do not expect the central values to agree. In fact both the number in our present study and that of Ref.~\cite{RBC:2021acc} agree with the dispersive results at their corresponding momenta. Accounting for the increased error growth due to this slight mismatch in the energy \cite{Parisi:1983ae} and the difference in number of measurements yields a factor $\approx 1.6 \times \sqrt{741/107}=4.2$, which implies the statistical errors are about the same for the two methods in this case.

The more salient comparison is for the $I=0$ case since that is the main reason for using GPBC. In Tab.~X in Ref.~\cite{RBC:2021acc} the energy of the two-pion state is given as 0.3479(11)[10], where the first error is statistical and the second is excited state systematic. It is determined from a three-operator, two-state fit with fit range 6--15. A GEVP analysis with $t_0=5$, $t-t_0=1$, without matrix subtraction, gives 0.3489(11), where the error is statistical. The matrix subtraction was not performed for this analysis as the ATW effects with GPBC are due to a long-time propagation of moving pions and the estimates of the size of the dominant ATW contribution gave values ten times smaller than the statistical errors on the data, and separate multi-operator fits with the ATW term included gave results consistent with zero.  The corresponding result for PBC is 0.406(16), where a matrix subtraction was done to remove ATW contamination. While accounting again for the mismatch in energies and measurements yields a factor of about 3.5, we do not see significant shift by the matrix subtraction for $I=0$ and it would be interesting to study and understand how significant the ATW effect on $I=0$ is.  In any case the factor is not big enough to explain the larger error which is about four times bigger in the PBC case for equivalent statistics and energies. Of course the argument about exponential growth of errors is not perfect since the pre-factors may differ between the two cases.

In the GPBC case, the boundary conditions lead to cubic symmetry breaking at the quark level which is suppressed by averaging pairs of single-pion interpolating operators with different quark and anti-quark momentum assignments but with the same total pion momentum~\cite{Christ:2019sah,RBC:2021acc}. For the ground state eight pairs (averages) of pion interpolating operators are used to construct sixty-four correlation functions, which most strongly overlap with the s-wave G-parity ground state. Here, we use six single-pion operators with momentum $(\pm 1,0,0)$ (plus permutations), or thirty-six correlators, for the corresponding s-wave excited state. While it is difficult to quantify the improvement gained from averaging over more combinations due to correlations, we do expect some benefit and will study this question in future calculations in both setups.

While not as direct, it is interesting and useful to compare the $24^3$ ensemble to the GPBC case. Here the PBC setup is the same as before, but the statistics are based on 258 configurations. The relevant two-pion energy ($t=4$) is 0.5298(64) in lattice units. Converting the errors to GeV and accounting for the different number of measurements yields a factor $\approx 2.53$. Accounting for the different number of pion momentum combinations as before, this factor is further reduced to somewhere between 1 and 2 (These factors are based on an assumption of statistical independence between the momentum orientations that may not be born out in practice). The exponential factor is roughly one in this case which we ignore.  A more thorough comparison between the GPBC and PBC approaches will be performed when we have measurements performed on an equivalent statistical sample.

\section{Conclusion}
\label{sec:conclusion}

In this work we have carried out a study of pion scattering at the physical point using 2+1 flavor M\"obius domain wall fermion ensembles with inverse lattice spacings of 1.023 and 1.378 GeV using periodic boundary conditions (PBC). The main focus was to extract the first excited state energies in the rest frame, in both $I=0$ and  2 channels, and their corresponding phase shifts, using the finite volume L\"uscher formalism~\cite{Luscher:1990ux}. The first excited state energy (roughly) corresponds to the important case of on-shell $K\to\pi\pi$ decay which is our longer term goal~\cite{Tomii:2022vxb}.

The energies were computed using the GEVP method~\cite{Luscher:1990ck,Blossier:2009kd}. In order to extract the desired first excited state and control excited state contamination, several two-pion interpolating operators were used, including a simple scalar bilinear. The other operators were constructed from pions with equal and opposite momenta. The single pion momentum took values (0,0,0), ($\pm1$,0,0), ($\pm1$,$\pm1$,0), ($\pm1$,$\pm1$,$\pm1$) (and permutations) in units of $2\pi/L$. As found in~\cite{RBC:2020kdj,RBC:2021acc}, the inclusion of the scalar bilinear with the quantum numbers of the vacuum is crucial to disentangle the first and second excited state energies.

The GEVP, or matrix correlation function, size ranged from $2\times2$ to $5\times5$. We found in most cases the size did not have a large effect after $2\times2$ for $I=2$ and $3\times3$ for $I=0$. However in some cases the noise of the higher states in the correlation matrix at large times adversely impacted the lower energies. In these cases the overlap of higher states with lower states was small, and the extra states only contributed noise to the lower states. This problem led to a modification of the method: as time increases the operator basis of the GEVP is changed (``re-based") using the eigenvectors of the GEVP from earlier times. While this method gives similar results in most cases studied here, it had a dramatic improvement for the $I=0$ first excited state on the $32^3$ ensemble where the statistics were relatively poor (see Fig.~\ref{fig:I=0 32c fits}). This RGEVP then allowed a relatively precise phase shift determination (see Fig.~\ref{fig:phase shift}).

We started this series of $K\to\pi\pi$ studies with G-parity boundary conditions (GPBC) \cite{Bai:2015nea,RBC:2020kdj,RBC:2021acc} because we anticipated it would be challenging to extract the signal of an excited two-pion state which is necessary for the $K\to\pi\pi$ study using PBC.
After seeing the successful calculation of two-pion scattering and $K\to\pi\pi$ decay amplitudes with GPBC \cite{RBC:2020kdj,RBC:2021acc}, we launched this PBC project to find a practical alternative, as a further check, and to enable calculations with isospin corrections.
We note once again that we have been successful in extracting signals of multiple states: the four lowest energy states for $I=2$ and the three lowest ones (at least on the $24^3$ lattice with better statistics) for $I=0$, despite the anticipated difficulty of extracting the excited state signals.

In this first study we have focused on two important systematic effects. First, by computing energies and phase shifts on ensembles with different lattice spacings we find no statistically significant discretization errors for $I=0$ and small, but statistically significant, effects for $I=2$. Second, we studied the time-dependence of the effective energies and observed noticeable excited state contamination for short times. We also saw that much of the effects arising from single-pion excited states could be removed with the dispersion relation method, which also removes leading discretization errors. Moreover, we expect the systematic uncertainties for PBC to be similar to GPBC. In that case statistical and systematic errors on the energies were estimated to be roughly equal\footnote{In Ref.~\cite{RBC:2021acc} only one lattice spacing was considered.}~\cite{RBC:2021acc}.

Another important motivation of employing PBC is that it appears difficult to use GPBC to compute QED and strong-interaction isospin breaking effects since they explicitly enforce isospin symmetry, though isospin symmetry breaking is expected to significantly impact the value of the direct CP violation parameter $\varepsilon'$ (20--30\%)\cite{Cirigliano:2019cpi}.  Controlling such effects precisely is important for the next generation of $\varepsilon^\prime$ calculations.

We are currently improving the statistics on both ensembles for our companion $K\to\pi\pi$ calculation, and when that is complete, we will have an even better comparison with G-parity and estimates of systematic errors. These results and the ones so far from the kaon decay project leave us optimistic for the PBC method.

\acknowledgments
We thank the U.S. Department of Energy (DOE) for partial support under awards DE-FG02-92ER40716 (TB, DH), DE-SC0010339 (TB, LJ, MT), DE-SC0046548 (ASM), DE-SC0012704 (PAB, TI, CJ, AS) and DE-SC001704 (AS). LJ and MT were also supported in part by the U.S. DOE under Early Career Award No. DE-SC002114.  The work of CK was supported by the U.S. DOE Exascale Computing Project.  PB, TI, CJ acknowledge support under DOE SciDAC-5 LAB 22-2580.
TI is also supported by the Department of Energy,
Laboratory Directed Research and Development (LDRD No. 23-051) of BNL and RIKEN BNL Research Center.
The research of MB was funded through the MUR program for young researchers, ``Rita Levi Montalcini''.  The work of DH was supported by the Swiss National Science Foundation (SNSF) through grant No. 200020\_208222. The software used for this work is based on \href{https://github.com/RBC-UKQCD/CPS}{CPS},
\href{https://github.com/paboyle/Grid}{Grid}, and \href{https://github.com/aportelli/hadrons}{Hadrons}.  We thank the collaborators in the RBC and UKQCD collaborations, especially N.H. Christ, for fruitful discussions.  Computations were carried out with USQCD resources funded by the US DOE at BNL and JLab.