Pion scattering in the isospin I=2 channel from elongated lattices

TL;DR

This study uses lattice QCD to analyze pion-pion scattering in the I=2 channel, extracting phase shifts and scattering lengths that agree with experimental data, across different lattice configurations and quark masses.

Contribution

First lattice QCD calculation of I=2 pion scattering using elongated lattices with a combined analysis at multiple quark masses.

Findings

Finite-volume spectra match experimental phase shifts.

Scattering lengths agree with chiral perturbation theory.

Results are consistent across different lattice geometries.

Abstract

Pion-pion elastic scattering in the isospin I=2 channel is investigated in two-flavor dynamical lattice QCD. Six ensembles are used with lattices elongated in one of the spatial dimensions at two quark masses corresponding to a pion mass of 315 MeV and 226 MeV. The energy of the low-lying states below the inelastic threshold are extracted in each case using the standard variational method.The extracted finite-volume spectrum is fitted by the inverse amplitude method simultaneously for both quark masses and extrapolated thereafter to the physical point. The resulting phase-shifts and scattering length are compared with those from experiment, leading-order chiral perturbation theory and other lattice studies. Our calculations match the experimental results.

Tables

Table 1. Table 1: Resolution of angular momentum in terms of irreps of the O h subscript 𝑂 ℎ O_{h} and the D 4 ​ h subscript 𝐷 4 ℎ D_{4h} group.

$\ell$	$O_h$	$D_{4h}$
0	$A_1^{+}$	$A_1^{+}$
1	$F_1^{-}$	$A_2^{-}\oplus E^{-}$
2	$E^{+}\oplus F_2^{+}$	$A_1^{+}\oplus B_1^{+}\oplus B_2^{+}\oplus E^{+}$
3	$A_2^{-}\oplus F_1^{-}\oplus F_2^{-}$	$A_2^{-}\oplus B_1^{-}\oplus B_2^{-}\oplus 2E^{-}$
4	$A_1^{+}\oplus E^{+}\oplus F_1^{+}\oplus F_2^{+}$	$2A_1^{+}\oplus A_2^{+}\oplus B_1^{+}\oplus B_2^{+}\oplus 2E^{+}$

Table 2. Table 2: The “seed” momentum 𝐩 𝐩 \mathbf{p} used to create interpolating fields for the zero momentum states and the moving states. The momentum components are indicated in units of the smallest non-zero momentum allowed in the corresponding direction; for the elongated boxes the smallest momentum in the z 𝑧 z -direction is reduced proportional to the elongation. Each of these operators is projected onto the A 1 + superscript subscript 𝐴 1 A_{1}^{+} irrep for the D 4 ​ h subscript 𝐷 4 ℎ D_{4h} group in the elongated case, or O h subscript 𝑂 ℎ O_{h} group for the cubic case (indicated with an star in the enesemble label).

		$𝐏=\left[000\right]$				$𝐏=\left[001\right]$
${ℰ}_1^{\ast }$		$\left[000\right]$	$\left[001\right]$	$\left[002\right]$		$\left[001\right]$	$\left[101\right]$	$\left[002\right]$
${ℰ}_2$		$\left[000\right]$	$\left[001\right]$	$\left[100\right]$		$\left[001\right]$	$\left[101\right]$	$\left[002\right]$
		$\left[101\right]$				$\left[111\right]$
${ℰ}_3$		$\left[000\right]$	$\left[001\right]$	$\left[100\right]$		$\left[001\right]$	$\left[101\right]$	$\left[002\right]$
		$\left[101\right]$	$\left[002\right]$	$\left[102\right]$		$\left[111\right]$	$\left[003\right]$	$\left[102\right]$
${ℰ}_4^{\ast }$		$\left[000\right]$	$\left[001\right]$			$\left[001\right]$	$\left[101\right]$	$\left[002\right]$
		$\left[000\right]$	$\left[001\right]$	$\left[100\right]$		$\left[001\right]$	$\left[101\right]$	$\left[002\right]$
		$\left[000\right]$	$\left[001\right]$	$\left[100\right]$		$\left[001\right]$	$\left[101\right]$	$\left[002\right]$

Table 3. Table 3: A summary of lattice details used including the lattice spacing a 𝑎 a , number of gauge configurations. The nucleon mass, pion decay constant and kaon decay constant are represented.

ensemble	$N_t\times N_{x,y}^2\times N_z$	$\eta$	$a\left[\text{fm}\right]$	$N_{\text{cfg}}$	$a{M}_{\pi }$	$a{m}_{u/d}^{pcac}$	$a{f}_{\pi }$
${ℰ}_1$	$48\times {24}^2\times 24$	$1.00$	$0.1210\left(2\right)\left(24\right)$	$300$	$0.1931\left(4\right)$	$0.01226\left(5\right)$	$0.0648\left(8\right)$
${ℰ}_2$	$48\times {24}^2\times 30$	$1.25$	$-$	$-$	$0.1944\left(3\right)$	$0.01239\left(4\right)$	$0.0651\left(6\right)$
${ℰ}_3$	$48\times {24}^2\times 48$	$2.00$	$-$	$-$	$0.1932\left(3\right)$	$0.01227\left(5\right)$	$0.0663\left(6\right)$
${ℰ}_4$	$64\times {24}^2\times 24$	$1.00$	$0.1215\left(3\right)\left(24\right)$	$400$	$0.1378\left(6\right)$	$0.00612\left(5\right)$	$0.0600\left(10\right)$
${ℰ}_5$	$64\times {24}^2\times 28$	$1.17$	$-$	$378$	$0.1374\left(5\right)$	$0.00620\left(4\right)$	$0.0600\left(8\right)$
${ℰ}_6$	$64\times {24}^2\times 32$	$1.33$	$-$	$400$	$0.1380\left(5\right)$	$0.00619\left(4\right)$	$0.0599\left(10\right)$

Table 4. Table 4: The fitted LECs ( l r i ⋅ 10 3 ⋅ superscript subscript 𝑙 𝑟 𝑖 superscript 10 3 l_{r}^{i}\cdot 10^{3} ) and m π subscript 𝑚 𝜋 m_{\pi} and f π subscript 𝑓 𝜋 f_{\pi} from the IAM analysis.

$l_r^1=+{11.6}_{-14.4}^{+16.2}$     $l_r^2=-{0.7}_{-7.4}^{+6.9}$     $l_r^4=+{52.4}_{-25.5}^{+25.1}$
$m_{\pi }^{\mathrm{light}}={223.83}_{-0.18}^{+0.19}$ MeV	$m_{\pi }^{\mathrm{heavy}}={315.25}_{-0.09}^{+0.06}$ MeV
$f_{\pi }^{\mathrm{light}}={97.45}_{-0.30}^{+0.56}$ MeV	$f_{\pi }^{\mathrm{heavy}}={107.44}_{-0.25}^{+0.16}$ MeV

Table 5. Table 5: Scattering length determined from the effective-range expansion.

	heavy	light	phys
$m_{\pi }$	315 MeV	226 MeV	139 MeV
$m_{\pi }{a}_0$	$-{0.1673}_{-0.0191}^{+0.0177}$	$-{0.1001}_{-0.0064}^{+0.0067}$	$-{0.0436}_{-0.0012}^{+0.0013}$

Table 6. Table 6: Energy levels in isospin 2 with fitting details. η 𝜂 \eta is the elongation and t 0 subscript 𝑡 0 t_{0} is the variational time.

$M_{\pi }=315\text{MeV}$
$𝐏$	$\eta$	$a{m}_{\pi }$	$n$	$t_0$	fit window	aE	${\chi }^2$/dof	Q
$(0,0,0)$	$1.0$	$0.1931\left(4\right)$	$1$	$3$	$5-15$	0.391(1)	1.36	$0.23$
			$2$	$3$	$7-17$	0.669(2)	1.24	$0.28$
	$1.25$	$0.1944\left(3\right)$	$1$	$3$	$4-13$	0.393(1)	1.23	$0.29$
			$2$	$3$	$6-14$	0.577(1)	1.10	$0.36$
			$3$	$3$	$5-13$	0.663(1)	0.98	$0.42$
	$2.0$	$0.1932\left(3\right)$	$1$	$3$	$4-13$	0.3883(6)	0.90	$0.48$
			$2$	$3$	$6-13$	0.4712(7)	0.76	$0.52$
			$3$	$3$	$2-12$	0.6509(6)	0.97	$0.44$
			$4$	$3$	$5-11$	0.6579(8)	0.76	$0.47$
			$5$	$3$	$5-12$	0.715(1)	0.65	$0.58$
$(0,0,1)$	$1.0$	$0.1931\left(4\right)$	$1$	$3$	$6-13$	0.526(1)	1.02	$0.38$
			$2$	$3$	$6-15$	0.762(4)	1.00	$0.42$
	$1.25$	$0.1944\left(3\right)$	$1$	$3$	$4-10$	0.487(1)	0.91	$0.40$
			$2$	$3$	$6-14$	0.729(1)	0.61	$0.66$
			$3$	$3$	$4-16$	0.754(2)	1.40	$0.19$
	$2.0$	$0.1932\left(3\right)$	$1$	$3$	$3-12$	0.4311(6)	1.06	$0.38$
			$2$	$3$	$6-12$	0.5630(8)	0.51	$0.60$
			$3$	$3$	$3-15$	0.6880(7)	1.04	$0.40$
			$4$	$3$	$2-16$	0.764(1)	0.68	$0.74$
			$5$	$3$	$4-10$	0.782(1)	1.01	$0.36$

Table 7. Table 7: Covariance matrices for each ensemble. The ordering of the levels is consistent with the order in Table  6 .

Equations

$$C_{ij}(t)=⟨{\mathcal{O}}_i(t){\mathcal{O}}_j^{†}(0)⟩.$$

$$C(t_0{)}^{-\frac{1}{2}}C(t)C(t_0{)}^{-\frac{1}{2}}{\psi }^{(n)}(t,t_0)={\lambda }^{(n)}(t,t_0){\psi }^{(n)}(t,t_0),$$

$${\lambda }^{(n)}(t,t_0)\propto e^{-E_{n}t}\left[1+\mathcal{O}(e^{-\Delta E_{n}t})\right],n=1,\dots ,N.$$

$$\pi \pi (\mathbf{P},\mathbf{p},t)\equiv {\pi }^{+}(\Gamma (\mathbf{p}),t){\pi }^{+}(\Gamma (\mathbf{P}-\mathbf{p}),t),$$

$${\pi }^{+}(\Gamma (\mathbf{p}),t)=\bar{d}(t)\Gamma (\mathbf{p})u(t).$$

$$\pi \pi (\mathbf{P},\mathbf{p})=\frac{1}{|G|}\sum _{g\in G}{\chi }_{A_1^{+}}(g)\pi \pi (R(g)\mathbf{P},R(g)\mathbf{p}).$$

$$\begin{array}{rl}C(\mathbf{P},& \mathbf{p},{\mathbf{p}}^{\prime },t)\equiv ⟨\pi \pi (\mathbf{P},{\mathbf{p}}^{\prime },t)\pi \pi (\mathbf{P},\mathbf{p},0{)}^{†}⟩\\ =& -[5(\mathbf{P}-{\mathbf{p}}^{\prime })t\mid 5(-\mathbf{P}+\mathbf{p})0\mid 5{\mathbf{p}}^{\prime }t\mid 5(-\mathbf{p})0]\\ & +[5(\mathbf{P}-{\mathbf{p}}^{\prime })t\mid 5(-\mathbf{P}+\mathbf{p})0][5{\mathbf{p}}^{\prime }t\mid 5(-\mathbf{p})0]\\ & -[5(\mathbf{P}-{\mathbf{p}}^{\prime })t\mid 5(-\mathbf{p})0\mid 5({\mathbf{p}}^{\prime })t\mid 5(-\mathbf{P}+\mathbf{p})0]\\ & +[5(\mathbf{P}-{\mathbf{p}}^{\prime })t\mid 5(-\mathbf{p})0][5({\mathbf{p}}^{\prime })t\mid 5(-\mathbf{P}+\mathbf{p})0].\end{array}$$

$$\left[i_1{\mathbf{p}}_1{j}_1\mid \dots \mid i_k{\mathbf{p}}_k{j}_k\right]\equiv \text{Tr}\prod _{\alpha =1}^k\Gamma ({\mathbf{p}}_{\alpha })M^{-1}(t_{j_{\alpha }},t_{j_{\alpha +1}}),$$

$$\lambda (t)=A_1{e}^{-E_1(t-t_0)}+A_1^{\prime }{e}^{-E_1^{\prime }(t-t_0)}+B{e}^{-\Delta E(t-t_0)}.$$

$$C(t)=\frac{1}{Z}\sum _n⟨n\left|e^{-H(T-t)}\pi \pi (\mathbf{P},{\mathbf{p}}^{\prime })e^{-Ht}\pi \pi (\mathbf{P},\mathbf{p}{)}^{†}\right|n⟩,$$

$$\delta C(t)=\frac{1}{Z}\sum _{\mathbf{p}}{e}^{-E_{\pi }(\mathbf{p})T}{e}^{-\Delta E(\mathbf{p})t},$$

$$\cot {\delta }^{20}={\mathcal{W}}_{00}=\frac{{\mathcal{Z}}_{00}(1,q^2;\eta )}{{\pi }^{3/2}\eta q}.$$

$$p\cot \delta (p)=\frac{1}{a_0}+\frac{1}{2}{r}_0{p}^2.$$

$$m_{\pi }{a}_0=-\frac{m_{\pi }^2}{16\pi f_{\pi }^2}\left[1+\frac{m_{\pi }^2}{32{\pi }^2{f}_{\pi }^2}\left(3\ln \frac{m_{\pi }^2}{2f_{\pi }^2}-1-l_{\pi \pi }\right)\right],$$

$$\displaystyle\cot\,\delta_{\rm IAM}(s)=\frac{\sqrt{s}}{2p}\Bigg{(}\frac{T_{2}(s)-\mathop{\hbox{Re}}{T_{4}(s)}}{(T_{2}(s))^{2}}\Bigg{)}\,,$$

$$m_{\pi }^2=m^2\left(1-\frac{m^2}{32{\pi }^2{f}_0^2}{\bar{l}}_3\right)\text{ \& }{f}_{\pi }=f_0\left(1+\frac{m^2}{16{\pi }^2{f}_0^2}{\bar{l}}_4\right).$$

$$l_i^r=\frac{{\gamma }_i}{32{\pi }^2}\left({\bar{l}}_i+\log \frac{m^2}{{\mu }^2}\right),$$

Full text

I Introduction

The prediction of scattering phase-shifts of strongly interacting systems directly from quark-gluon dynamics has become possible through rapid advances in lattice QCD calculations. In such ab-initio calculations one can also vary parameters that are inaccessible in experiment like the number of flavors $N_{f}$. In addition, they serve to test models that extrapolate in hadron masses, allowing for deeper insights into the QCD dynamics at the hadronic scale, including resonances and their dependence on quark masses and flavors.

Lattice QCD calculations are performed in a small cubic volume and in imaginary time, so that a direct evaluation of phase-shifts is not possible. However, the discrete energy eigenvalues of the QCD Hamiltonian in a cube with periodic boundary conditions can still be put into relation with phase-shifts as shown by Lüscher Lüscher (1991, 1986). Varying the lattice size, through, e.g., elongated boxes Feng et al. (2004); Lee and Alexandru (2017), allows one to produce a number of phase-shifts at the eigenergies; additional phase-shifts can be determined by projecting the particles to moving frames Rummukainen and Gottlieb (1995). This enabled lattice QCD calculations of the $\rho$ meson phase-shifts in $N_{f}=2$ Aoki et al. (2007); Feng et al. (2011); Göckeler et al. (2008); Lang et al. (2011); Pelissier and Alexandru (2013); Bali et al. (2016); Guo et al. (2016) and $N_{f}=2+1$ Aoki et al. (2011); Dudek et al. (2013); Feng et al. (2015); Metivet (2015); Wilson et al. (2015); Andersen et al. (2019); Fu and Wang (2016); Alexandrou et al. (2017). The isoscalar sector is particularly interesting because of the presence of the broad $f_{0}(500)$ “$\sigma$” resonance. Only recently, phase-shifts in this channel have been calculated Briceño et al. (2017a, 2018); Guo et al. (2018a) (for calculations of the scattering lengths, see, e.g., Refs. Fu and Chen (2018); Liu et al. (2017)).

The infinite-volume extrapolation of the $\pi^{+}\pi^{+}$ system was the first physical application of the original Lüscher formalism, for the scattering length Sharpe et al. (1992); Kuramashi et al. (1993); Gupta et al. (1993); Yamazaki et al. (2004); Aoki et al. (2005); Beane et al. (2006, 2008a); Feng et al. (2010); Yagi et al. (2011); Fu (2013); Sasaki et al. (2014); Kawai et al. (2018) and extended to higher energies Dudek et al. (2012); Bulava et al. (2016); Akahoshi et al. (2019). See also Refs. Helmes et al. (2015, 2019) for a comprehensive calculation of meson-meson scattering lengths. Systems of more than two pions with maximal isospin have been calculated in lattice QCD Hörz and Hanlon (2019); Detmold et al. (2008) and serve as a test ground for infinite-volume mapping techniques that are currently being developed Blanton et al. (2019); Pang et al. (2019); Briceño et al. (2019); Mai and Döring (2019); Romero-López et al. (2018); Guo et al. (2018b); Hammer et al. (2017); Briceño et al. (2017b); Hansen and Sharpe (2015); Detmold and Savage (2008). In Ref. Mai and Döring (2019) such a formalism was, for the first time, applied to the above-threshold $3\pi^{+}$ system.

In this study we determine the $\pi^{+}\pi^{+}$ isospin $I=2$ phase-shifts using elongated boxes. Our study is carried out using $N_{f}=2$ dynamical configurations with nHYP fermions Hasenfratz et al. (2007). We analyze two sets of ensembles with different sea quark masses: one corresponding to $m_{\pi}=315$ MeV and the other one to $m_{\pi}=226$ MeV. For each pion mass we use three ensembles with different lattice geometry. For each ensemble we analyze states at rest, ${\bf P}=(0,0,0)$, and states moving along the elongated direction with momentum ${\bf P}=(0,0,1)$. For each case, we use two-hadron $\pi\pi$ interpolators with different back-to-back momenta in the variational basis. We extract the lowest states energy states using the variational method Luscher and Wolff (1990).

The scattering length is computed at the two pion masses and the results are extrapolated to the physical point using chiral perturbation theory. In the wider energy range the inverse amplitude method (IAM) Truong (1988); Gómez Nicola et al. (2008) is utilized to fit the data and obtain predictions for the phase-shifts at the physical quark mass. This model is unitary and matches the chiral pion-pion amplitude Gasser and Leutwyler (1985, 1984) up to the next-to-leading order. The obtained predictions overlap nicely with the experimentally obtained values.

In comparison to previous studies in which this method was applied Hu et al. (2017); Döring et al. (2018); Hu et al. (2016); Bolton et al. (2016) we allow here the pion mass and decay constant to vary and include their full correlations with the energy eigenvalues in the fit.

In an upcoming paper Mai et al. (2019) we will use the $I=2$ energy eigenvalues determined here, together with the corresponding results of the isovector and isoscalar channels Guo et al. (2016, 2018a), to perform a global analysis of pion-pion-scattering with IAM. The full energy dependence of the $I=2$ partial-wave amplitude is also needed as input for upcoming calculations of the of $3\pi^{+}$ system above threshold along the lines of Ref. Mai and Döring (2019), as required by three-body unitarity Mai et al. (2017).

This paper is organized as follows. In Section II details of the lattice calculation are provided, including the variational basis and interpolating operators, followed by a description of the extraction of energy eigenvalue in Section III. The determination of phase-shifts, effective range expansion fit, IAM fit, and chiral extrapolation is discussed in Section IV.

II Lattice setup

For the $I=2$ channel we cannot use $\bar{q}q$ interpolating fields since the maximum isospin for such operators is $1$. We will construct our variational basis out of two-pion interpolating fields, two $\pi^{+}$ pions to be precise, projected to the appropriate momentum combinations. Different interpolating fields will only differ by the choice of pion momenta. Using the interpolating fields we construct the correlator matrix:

[TABLE]

The eigenvalues of this matrix are obtained by solving the generalized eigenvalue problem,

[TABLE]

where $t_{0}$ is a parameter. The energies of the $\pi^{+}\pi^{+}$ system can be extracted from the long-time behavior of the eigenvalues, Luscher and Wolff (1990); Blossier et al. (2009)

[TABLE]

When $t<2t_{0}$ the correction term vanishes at a rate given by $\Delta E_{n}=E_{N+1}-E_{n}$, the difference between the energy level of interest and the lowest level excluded from the variational basis Blossier et al. (2009). The variational method filters out contaminations in the spectrum from states included in the basis.

For each ensemble we include in the variational basis all interpolating fields with momenta corresponding to two-hadron states that in the absence of interactions will be below the inelastic thresholdat $E=4m_{\pi}$, and the next momentum just above it. In our analysis we include only the energies extracted from the variational method that lay below the inelastic threshold.

In this work we will use two cubic lattices, and four lattices elongated in one spatial direction. The relevant symmetry groups are $O_{h}$ and $D_{4h}$ for the cubic and elongated boxes, respectively. The elongation direction is chosen to be $z$. The angular momentum labels used for the irreducible representations (irreps) of $SO(3)$ are split into multiplets of the irreps of the lattice symmetry groups. The splitting of angular momentum for $O_{h}$ and $D_{4h}$ can be found in Table 1.

To get additional states in this scattering region, we will also use interpolators for non-zero momentum states. By constructing operators with a total momentum $\mathbf{P}$ we will access additional energy levels in the elastic region. The relativistic effects cause the box to shrink in the direction of the total momentum, changing the symmetry group. We align the boost direction with the direction of elongation so that the relevant symmemtry group remains $D_{4h}$.

The $\pi^{+}\pi^{+}$ interpolating fields are constructed using two $\pi^{+}$ interpolators with the appropriate momenta:

[TABLE]

where

[TABLE]

Here the quark fields $u(t)$ and $d(t)$ represent a three-dimensional slice of the field, and they can be viewed as $N=4\times 3\times V_{3}$ vectors, where $V_{3}$ is the number of points in a time slice. The matrix $\Gamma$ is an $N\times N$ matrix that acts in the spin, color and position space. The only structure used in this study is $\Gamma(\mathbf{p})=\gamma_{5}e^{i\mathbf{p}}$. This matrix acts trivially in the color space and in position space we have $[e^{i\mathbf{p}}]_{\mathbf{x},\mathbf{y}}=e^{i\mathbf{p}\mathbf{x}}\delta(\mathbf{x}-\mathbf{y})$. For more details see Refs. Pelissier and Alexandru (2013); Guo et al. (2016).

To access the zero angular momentum, $l=0$, phase-shifts in $\pi^{+}\pi^{+}$ we need to project our operators to the $A_{1}^{+}$ irrep of $O_{h}$ and $D_{4h}$. According to Table 1 the lowest contribution to this irrep comes from $l=0$ states and the corrections come from $l=4$ states for $O_{h}$ and $l=2$ states for $D_{4h}$. These higher phase-shifts are expected to be small in the kinematic region we explore and can be safely neglected.

To construct interpolating fields with the right symmetry properties, we start with a “seed” interpolating field and project it on the $A_{1}^{+}$ irrep using:

[TABLE]

Here $G$ is the group ($O_{h}$ or $D_{4h}$), $g$ is an element of $G$, $\chi$ is the character of $g$ in the irrep $A_{1}^{+}$ and $R(g)$ is the rotation corresponding to the element $g$.

The collection of “seed” momenta used in this study are presented in Table 2. As mentioned earlier, we used all momenta that in the non-interacting case had energy less than $4m_{\pi}$ and the next just above. Then the operator is plugged into Eq. (6) to generate the linear combination that overlaps with $A_{1}^{+}$. In Fig. 1 we highlight which operators were used in the $315\mathop{\hbox{MeV}}$ ensembles. For example, at elongation $\eta=1.25$ there are three energies in the elastic scattering region, thus four operators are selected.

These operators are evaluated on a set of six ensembles. The parameters for these ensembles are listed in Table 3. To compute the correlation functions we need to perform the Wick contractions and the correlation functions become functions of quark propagators. All correlation functions used in this study are linear combinations of

[TABLE]

Above we introduced the following notation to denote the quark propagator traces:

[TABLE]

where $M^{-1}(t,t^{\prime})$ represents the quark propagator $N\times N$ matrix between two time-slices.

Evaluating these diagrams requires the all-to-all quark propagator. This is a very expensive calculation and to avoid it we use the Laplacian-Heaviside method (LapH) Peardon et al. (2009). The basic idea is to replace the quark interpolating fields with smeared quark interpolators that have the same symmetries. The smearing is introduced by truncating the three-dimensional Laplacian operator on each time-slice by keeping the lowest-lying $N_{v}$ operator modes. The net effect is that the point-quark all-to-all propagator $M^{-1}$ is replaced with the propagator for smeared quarks ${\widetilde{M}}^{-1}$. The smeared propagator can be computed efficiently since we only have to invert the Dirac matrix for the LapH modes, rather than for each point on the lattice. We stress that this replacement is not an approximation, but rather creates operators with the right quantum numbers which have a different overlap with the relevant states. The number of LapH modes $N_{v}$ controls the smearing radius and the overlap; if the truncation is aggressive, keeping only a few LapH modes, the overlap decreases and correlation functions need to be computed accurately at large time separations to resolve the relevant states. For our study we used $N_{v}=100$ which corresponds to a smearing radius of roughly $0.5\mathop{\hbox{fm}}$ Guo et al. (2016). To compute the smeared propagator efficiently we use GPU inverters Alexandru et al. (2012). In our calculations all steps remain unchanged except that the quark traces in Eq. (8) are computed using the smeared propagator ${\widetilde{M}}^{-1}$.

III Extracting finite-volume spectrum

The variational method gives us the generalized eigenvalues as a function of time separation between sink and source. To extract the energy levels we have to fit these functions. When the correlation functions are saturated by a single state, the mass can be extracted using a single-exponential fit. To determine the appropriate range for this fit we plot the effective mass, $m_{\text{eff}}\equiv-\log\lambda(t+1)/\lambda(t)$, and look for a plateau. As we can see from Fig. 2, the data is very precise and the effect of the corrections is statistically significant for all time ranges (see the inset). As such, we have to fit the correlators to include the effects of the excited states and thermal corrections:

[TABLE]

Above the fitting parameters are the energy levels $E_{1}$ and $E_{1}^{\prime}$, the spectral weights $A_{1}$ and $A_{1}^{\prime}$, and $B$ the coefficient of the leading thermal correction, to be discussed below. The inclusion of the excited level allows us to fit the correlator at earlier times and $E^{\prime}_{1}$ should have a value corresponding to the lowest states excluded from the variational basis. This is indeed the case in our analysis.

One unusual feature for two-pion correlation functions is that the wrap-around effects, due to thermal fluctuations, lead to constants or slowly decaying exponential contributions to the correlation function Dudek et al. (2012). To see how the thermal effects contribute, consider the correlation function

[TABLE]

where $T$ is the time-extent of the lattice. The leading order contribution comes from $\left|n\right\rangle=\left|0\right\rangle$, which leads to the usual superposition of exponentials for the correlation function. The next order term for $\pi\pi$ comes from $\left|n\right\rangle=\left|\pi(\mathbf{p})\right\rangle$, single pion states with momentum $\mathbf{p}$. The effect of these states is to generate a set of slowly decaying exponentials

[TABLE]

with $E_{\pi}(\mathbf{p})$ the energy of a single pion with momentum $\mathbf{p}$ and $\Delta E=E_{\pi}(\mathbf{P}+\mathbf{p})-E_{\pi}(\mathbf{p})$. Note that these contributions are suppressed by coefficients that vanish exponentially fast with $T$. However, when fitting our correlators for $t\sim T/2$, their effects are important, as can be seen from Fig. 2, and this correction needs to be included. In our fits we fix the exponential term $\Delta E$ to the lowest value generated by varying $\mathbf{p}$ over the allowed momenta: for $\mathbf{P}=0$ we fix $\Delta E=0$ and for $\mathbf{P}=[001]$ we fix $\Delta E=\sqrt{m_{\pi}^{2}+(2\pi/L/\eta)^{2}}-m_{\pi}$.

For elongated boxes an interesting feature arises: the non-interacting levels cross as we vary the elongation. This can be easily seen in Fig. 1. In particular, we draw the reader’s attention to the intersection between the $\pi\pi([000],[100])$ and $\pi\pi([000],[002])$ that occurs at $\eta=2$. This has consequences also for the energy levels in the interacting case. Even in the presence of interactions, states appear with energies at (or very close to) non-interacting levels. To understand this note that the interacting energy levels fall in between two non-interacting ones. This can be understood in the case that the only partial wave contributing is $\ell=0$. In this case the energy levels satisfy Eq. (12) below and the ${\cal Z}_{00}$ function has poles at the energies corresponding to non-interacting levels, bracketing the interacting solutions. In Fig. 3 we show the solutions of Eq. (12) as we vary the elongation of the box past $\eta=2$. Notice the vertical lines that are the poles of the ${\cal Z}_{00}$ functions: one corresponding to $\pi\pi([000],[100])$ state that remains fixed since the momentum in the transverse directions does not change as we change elongation, and one corresponding to $\pi([000],[002])$ that moves as we vary $\eta$. The two lines bracket one of the solutions of Eq. (12) and when $\eta=2$ the poles merge and the corresponding solution has exactly the same energy as the non-interacting two-pion state.

Note that this argument relies on Eq. (12) being exact, which only works in the limit where the higher partial waves are zero. We expect that when the other partial waves are taken into account, the energies will be shifted slightly away from this position. This is for example the case when the poles for $\pi\pi([000],[100])$ and $\pi\pi([000],[001])$ merge, which happens for $\eta=1$, corresponding to a cubic box. In that case, the energy level that is pinched by the poles ends up belonging to the $E^{+}$ representation of the symmetry group for cubic boxes $O_{h}$. The lowest angular momentum that this irrep overlaps with is $\ell=2$. In this case it is then natural that the energy for this level coincides with the non-interacting case if we assume that all partial waves above $\ell=0$ are small: the shift will be proportional to the $\delta^{I=2,l=2}$.

It is important to emphasize that the presence of these levels do not imply that the phase-shift $\delta^{20}$ is zero for this energy, as one would expect. For example in Fig. 3 the phase-shift is clearly non-zero at the pole. To understand this note that the connection between the phase-shift and energy is controlled by the ${\cal Z}_{00}$ function which is infinitely steep when the poles merge. Thus, finite changes in the phase-shift away from zero lead to infinitesimal (zero) changes in the energy away from the pole. This is also important when analyzing the levels extracted from numerical simulations. The energy levels will be determined with some finite stochastic error, which is mapped through this infinitely steep function into infinite errors in the phase-shift space. As such, these energy levels offer no real constrain on the phase-shifts and we do not include them in our analysis. Note that large error bars in phase-shifts always arise when the error-bars cross a pole, which happens even without pinching (see for example the left panel of Fig. 5).

IV Extracting physical quantities

Having computed the finite-volume energy levels of $\pi\pi$ scattering in maximal isospin, we now connect them to physical quantities. For energies below the inelastic threshold we use Lüscher’s formula Lüscher (1991) and the extensions to elongated boxes Feng et al. (2004) and boosted frames Lee and Alexandru (2017). The s-wave is the lowest partial wave that contributes to $\pi\pi$ scattering in isospin 2. For cubic boxes, elongated boxes, and boosted systems, the irrep that overlaps with $\delta^{20}$ is $A_{1}^{+}$. Higher partial waves ($l=2,4,6\ldots$) also overlap with this irrep in a finite volume, however, we neglect these from further analysis since they are known to be negligible in this channel, see, e.g., Refs. Pelaez (2016); Dudek et al. (2012). Thus, for the extraction of the phase-shifts in this channel the relevant Lüscher formula is

[TABLE]

Boosting the system has the effect of changing the value of $\eta$, see Ref. Lee and Alexandru (2017) for further details and explicit form of $\mathcal{Z}_{00}(1,q^{2};\eta)$. In the following we will use two distinct parameterizations of the phase-shifts, which also will allow us to interpolate them in energy as well as to extrapolate them to the physical point.

IV.1 Effective range expansion

The effective range expansion (ERE) is expected to hold in channels without resonances in vicinity of the production threshold. The drawback is that with a finite number of terms we can only describe low momentum data. We will use the first two terms of the ERE, i.e.,

[TABLE]

The parameters $a_{0}$ and $r_{0}$ are the scattering length and effective range, respectively. The lattice energy spectrum can be used to compute phase-shifts, and we can find an $a_{0}$ and $r_{0}$ from a fit to this results. To perform this fit, the above-given phase-shift is related to the finite-volume spectrum via the Lüscher’s formula. The corresponding correlated $\chi^{2}$ with the energy eigenvalues is minimized then with respect to $a_{0}$ and $r_{0}$.

For both pion masses we restrict the fit to the lowest two energy levels in each ensemble. For the $315\mathop{\hbox{MeV}}$ data the scattering length extracted is $m_{\pi}a_{0}=-0.20(3)$ with a $\chi^{2}/\text{dof}=5.6/(6-2)$. The $226\mathop{\hbox{MeV}}$ ensemble yields a scattering length $m_{\pi}a_{0}=-0.10(2)$ with a $\chi^{2}/\text{dof}=3.3/(6-2)$. The $315\mathop{\hbox{MeV}}$ scattering length is within two sigma of the leading order (LO) ChPT value for that pion mass, while the $226\mathop{\hbox{MeV}}$ scattering length is consistent with the LO ChPT value within one sigma.

We can extrapolate these scattering lengths to the physical pion mass using ChPT at NLO Gasser and Leutwyler (1984). The expansion of $m_{\pi}a_{0}$ reads

[TABLE]

where $l_{\pi\pi}$ is a combination of the usual low-energy constants (LECs). This function is fit to our two extracted scattering lengths at unphysical pion masses, with respect to $l_{\pi\pi}$. With a $\chi^{2}/\text{dof}=0.74/(2-1)$ we obtain a value of $l_{\pi\pi}=-1.09(2.52)$. Evaluating the ChPT expansion at the physical pion mass we obtain $m_{\pi}a_{0}=-0.0455(16)$. We compare this result to the physical value and other lattice studies in Fig. 4 indicated as “ERE.” The gray band represents the result of the estimation from the Roy equation Ananthanarayan et al. (2001). The references for other determinations of the scattering length are, in order: NPLQCD2006 Beane et al. (2006), NPLQCD2008 Beane et al. (2008a), ETMC2010 Feng et al. (2010), ETMC2015 Helmes et al. (2019), Yagi2011 Yagi et al. (2011), Fu2013 Fu (2013), and PACS-CS2013 Sasaki et al. (2014). This figure also contains a second scattering length from the lattice data using the inverse amplitude method as explained below, indicated as “IAM.”

IV.2 Inverse amplitude method

The available data on the finite-volume spectrum of $\pi\pi$ scattering at maximal isospin covers a large region in energy and pion mass. Therefore, to connect those data to each other but also to the physical point one requires a framework that can extrapolate in energy and $m_{\pi}$. The so-called Inverse Amplitude Method (IAM) Truong (1988); Gómez Nicola et al. (2008, 2008) reconciles both these requirements. In particular, it provides a scattering amplitude, which fulfills two-body unitarity exactly and has correct chiral behaviour up to next-to-leading chiral order Gasser and Leutwyler (1985, 1984). Furthermore, it fulfills the general requirements on the chiral trajectory for resonances, derived in Ref. Bruns and Mai (2018) to all chiral orders. The corresponding cotangent of the phase-shift reads

[TABLE]

where $T_{2}(s)$ and $T_{4}(s)$ are the leading and next-to-leading order chiral amplitudes, respectively, projected to isospin $I=2$ and angular momentum $l=0$. The total energy squared of the system is denoted by the Mandelstam $s=E^{2}$.

The leading order chiral amplitude $T_{2}(s)$ is a function of energy, Goldstone-boson mass, $m^{2}=B(m_{u}+m_{d})$ and pion decay constant in the chiral limit, $f_{0}$. The amplitude $T_{4}$ involves in the two-flavor case two low-energy constants (LECs) $\bar{l}_{1}$ and $\bar{l}_{2}$. Two additional low-energy constants $\bar{l}_{3}$, $\bar{l}_{4}$ enter the NLO chiral amplitude when replacing the above mass and decay constants by their physical (lattice) values using one-loop results Gasser and Leutwyler (1984),

[TABLE]

The constants $\bar{l}_{i}$ do not depend on the regularization scale, but only on the parameters of the underlying theory - the quark masses. However, they are related to the scale-dependent, but quark-mass independent renormalized LECs via

[TABLE]

where $\gamma_{1}=\nicefrac{{1}}{{3}}$, $\gamma_{2}=\nicefrac{{2}}{{3}}$, $\gamma_{3}=-\nicefrac{{1}}{{2}}$, $\gamma_{4}=2$. For a fixed scale $\mu$ one can, thus, determine the renormalized LECs and then make predictions for the two-particle scattering for a setup with a different pion mass. In the course of this work we use dimensional regularization with $\mu=770$ MeV, but emphasize that the expression (15) is manifestly scale independent.

As discussed before, the finite volume spectrum consists of a large set of energy eigenvalues as well as $m_{\pi}$ and $f_{\pi}$. In the past it has been noted Bulava et al. (2016) that variations of the pion mass and decay constant can lead to non-negligible effects in the energy spectrum. Therefore, it is important to asses this source of uncertainty in a systematic way. In relation to this we noted some systematic effects in $\mathcal{E}_{2}$, c.f., different central value of the pion mass recorded in Table 3. Thus, we exclude this set from further fits, which leaves us with 21 energy eigenvalues, five pion masses and five decay constants to fit. The latter is performed plugging in Eq. (15) into Eq. (12) while minimizing the correlated $\chi^{2}$ with respect to 3 dynamical parameters 111We found that the value of $l_{r}^{3}$ does not lead to any notable improvement of the fit and fix this value to the one reported by FLAG Aoki et al. (2019): $l_{r}^{3}=8.94\cdot 10^{-6}$. $l_{r}^{1}$, $l_{r}^{2}$ and $l_{r}^{4}$. We let the fit also determine the values of the pion mass and decay constants that are common for the three light and three heavy ensembles, respectively. Note that all correlations between pion mass, decay constant and energy eigenvalues are taken into account. The best parameters of the overall fit to light and heavy data produce $\chi^{2}/{\rm dof}=75.4/(31-7)$ with the LECs given in Table 4.

The errors on the parameter have been estimated in a re-sampling procedure. We have explored the origin of the large $\chi^{2}$ of this combined fit extensively. As a matter of fact, the fit to the results of only heavy or light ensembles gives $\chi^{2}/{\rm dof,heavy}=21.2/(17-5)$ and $\chi^{2}/{\rm dof,light}=7.0/(14-5)$. Thus, it appears that IAM is flexible enough to address all lattice QCD results up to the inelastic threshold as a function of energy. However, the simultaneous parametrization of the pion mass dependence is less reliable for such a large pion mass range. To some extent this is expected, since IAM coincides with the chiral expansion only up to next-to-leading chiral order. Furthermore, we refer at this point to the recent study of systematic uncertainties tied to the use of lattice data at unphysical pion mass in the context of ChPT Dürr (2015).

We note that, for the individual fits where the model fits well the data, a significant part of the $\chi^{2}$ comes from the cross-correlation between the data points. This demonstrates that the cross-correlations should be taken into account.

Table 5 lists the scattering lengths for the obtained set of parameters.

The result at the physical pion mass is depicted together with other lattice QCD based estimations in Fig. 4. It shows that the extrapolated value is well in agreement with the value extracted from the Roy equations Ananthanarayan et al. (2001) as well as overlaps with most recent lattice QCD determinations.

The compilation of phase-shifts for the combined fit to light and heavy energy eigenvalues is depicted in Fig. 5. It also contains the lattice results after extrapolating them to the physical point (right panel). Recall that fitting is performed at the level of energy eigenvalues including full information about cross correlations. Overall we note, that while the LECs of the model have sizable statistical uncertainties, the corresponding error bands on phase-shifts are quite narrow. Clearly, the prediction of the phase-shifts at the physical point depicted in Fig. 5 shows a nice agreement with the available experimental data even far beyond the elastic region. This observation has been also noted in Ref. Mai and Döring (2019) where the two and three pion lattice QCD results of the NPLQCD collaboration Beane et al. (2008b); Detmold et al. (2008) were discussed.

V Conclusion

We have performed a calculation of the pion-pion elastic scattering in the isopsin $I=2$ channel in two-flavor dynamical lattice QCD. One novel feature compared with other $I=2$ studies is the use of elongated lattices which has proven a cost-effective approach in mapping out the momentum dependence in scattering processes. We considered six ensembles (see Table 3) with elongation up to a factor of 2 in one of the spatial dimensions, and two pion masses at 315 MeV and 226 MeV. Boosting of the $\pi\pi$ system in the elongated direction is also considered for enhanced energy coverage. In each case, we extract multiple low-lying states using the standard variational method.

The data is analyzed simultaneously at both pion masses making use of the inverse amplitude method. This unitary model matches the chiral $\pi\pi$ scattering amplitude to next-to-leading order and allows, thus, for chiral extrapolation to the physical point in a wide energy range. In the current application of the method, we allow the pion mass and decay constant to vary and include their full correlations with the energy eigenvalues in the corresponding fit. The scattering length extrapolated to the physical point reads $m_{\pi}{}a_{0}=-0.0436(13)$. As an additional check we also perform an effective range expansion, which in combination with the perturbative chiral form of the scattering length gives a consistent value of $m_{\pi}{}a_{0}=-0.0455(16)$ at the physical point.

A closer analysis of the $\chi^{2}$ reveals a slight tension in the model description between light and heavy ensembles if the full energy range is fitted. Overall, this suggests that more reliable extraction of low-energy parameters requires results at lower pion mass as well. Another possibility, is to include cross channels, $I=0,1$, in the analysis, which might restrict these constants more strongly. In any case, we note that the phase-shifts extrapolated to the physical point match the experimental phase-shifts.

Overall our results demonstrate the efficacy of using elongated lattices and IAM analysis as an effective tool in studying hadron-hadron scattering processes from first principles. Work is under way to carry out a global IAM analysis across all three isospin channels of $\pi\pi$ scattering on elongated lattices. Such a study can pave the way for extending to systems that include three pions above threshold.

Acknowledgements.

CC, AA, and FXL are supported in part by the National Science Foundation CAREER grant PHY-1151648 and by U.S. DOE Grant No. DE-FG02-95ER40907. MD and MM acknowledge support by the National Science Foundation CAREER grant PHY-1452055. MD is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177 and grant de-sc001658MM. AA gratefully acknowledges the hospitality of the Physics Department at the University of Maryland where part of this work was carried out. The computations were performed on the GWU Colonial One computer cluster and the GWU IMPACT collaboration clusters.

Appendix A Extracted energies and correlation matrices

In this section we summarize results from fitting the correlation function to extract the finite volume energies. The energies from non-boosted systems were fit with the functional form $A_{1}e^{-E_{1}t}+A_{1}^{\prime}e^{-E_{1}^{\prime}t}+B$, while boosted systems were fit with $A_{1}e^{-E_{1}t}+A_{1}^{\prime}e^{-E_{1}^{\prime}t}+Be^{-\Delta Et}$ with $\Delta E=\sqrt{m_{\pi}^{2}+[2\pi/(\eta L)]^{2}}-m_{\pi}$. The results of all extractions are reported in Table 6. The fit window was chosen to minimize the $\chi^{2}$ per degree of freedom. The parameter Q is the confidence level of the fit corresponding to the probability that $\chi^{2}$ is larger then the fit result. The value of $am_{\pi}$ is extracted from the two point correlation function. We also include tables of the cross-correlations between the extracted energies for all $\pi\pi$ scattering channels in Table 7.