A cross-channel study of pion scattering from lattice QCD

TL;DR

This study performs a comprehensive analysis of pion-pion scattering using lattice QCD data across multiple isospin channels, employing a chiral model to improve parameter determination and physical extrapolation.

Contribution

It introduces a simultaneous multi-channel analysis of lattice QCD pion scattering data with a chiral model, enhancing parameter accuracy and physical extrapolation over single-channel methods.

Findings

Precise determination of scattering lengths: $m_ a_0^{I=0}=0.2132(9)$ and $m_ a_0^{I=2}=0.0433(2)$.

Extraction of resonance parameters: $m_ \sigma=443(3)-i221(6)$ MeV, $m_ ho=724(4)-i67(1)$ MeV.

Improved parameter constraints through cross-channel fitting.

Abstract

We use a chiral model for pion interactions, in the inverse amplitude formalism, to perform a simultaneous analysis of lattice QCD results for pion-pion scattering in all three isospin channels. The input is the finite-volume two-pion spectrum computed using lattice QCD from six ensembles on lattices elongated in one of the spatial dimensions. A two-flavor dynamical lattice QCD action is used with two quark masses corresponding to a pion mass of 315 MeV and 224 MeV. The spectrum in the elastic region is subjected to a global fit which takes into account full correlations across isospin, pion mass and decay constant. The parameters from the fit are used to perform a chiral extrapolation to the physical point. The cross-channel fit results in a more precise determination of the parameters of the model when compared with single channel fits. We obtain $m_\pi a_0^{I=0}=0.2132(9)$, and $m_\pi a_0^{I=2}=0.0433(2)$ as well as $m_\sigma=443(3)-i221(6)$ MeV and $m_\rho=724(4)-i67(1)$ MeV. Several aspects of scale setting and consistency with previous analyses of lattice QCD results are discussed as well.

Tables

Table 1. Table 1: Summary of the ensembles used in all isospin channels. Included for reference are the lattice geometry and elongation  η 𝜂 \eta , the lattice spacing  a 𝑎 a , the number of configurations in each ensemble, the pion mass, the current quark mass (see Eq.  5 ), and the pion decay constant.

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$		$300$	$0.1944\left(3\right)$	$0.01239\left(4\right)$	$0.0651\left(6\right)$
${ℰ}_3$	$48\times {24}^2\times 48$	$2.00$		$300$	$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 2. Table 2: 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	$T_1^{-}$	$A_2^{-}\oplus E^{-}$
2	$E^{+}\oplus T_2^{+}$	$A_1^{+}\oplus B_1^{+}\oplus B_2^{+}\oplus E^{+}$
3	$A_2^{-}\oplus T_1^{-}\oplus T_2^{-}$	$A_2^{-}\oplus B_1^{-}\oplus B_2^{-}\oplus 2E^{-}$
4	$A_1^{+}\oplus E^{+}\oplus T_1^{+}\oplus T_2^{+}$	$2A_1^{+}\oplus A_2^{+}\oplus B_1^{+}\oplus B_2^{+}\oplus 2E^{+}$

Table 3. Table 3: The fitted LECs ( l i r ⋅ 10 3 ⋅ subscript superscript 𝑙 𝑟 𝑖 superscript 10 3 l^{r}_{i}\cdot 10^{3} ) and m π subscript 𝑚 𝜋 m_{\pi} and f π subscript 𝑓 𝜋 f_{\pi} from the mIAM fit ( χ dof 2 = 218 / ( 88 − 7 ) subscript superscript 𝜒 2 dof 218 88 7 \chi^{2}_{\rm dof}=218/(88-7) ) to all lattice QCD results  Culver  et al. ( 2019 ); Guo  et al. ( 2018 , 2018 ) .

$l_1^r=-{4.07}_{-0.13}^{+0.12}$     $l_2^r=+{5.14}_{-0.19}^{+0.23}$     $l_4^r=+{9.05}_{-0.70}^{+0.54}$
$m_{\pi }^{\mathrm{light}}={223.98}_{-0.13}^{+0.07}$ MeV	$m_{\pi }^{\mathrm{heavy}}={315.40}_{-0.05}^{+0.09}$ MeV
$f_{\pi }^{\mathrm{light}}={98.29}_{-0.12}^{+0.22}$ MeV	$f_{\pi }^{\mathrm{heavy}}={107.48}_{-0.12}^{+0.08}$ MeV

Table 4. Table 4: Scattering lengths and pole positions from the global mIAM fit to lattice QCD data  Culver  et al. ( 2019 ); Guo  et al. ( 2018 , 2018 ) . The last column shows the extrapolation to the physical point. Error bars are determined in a re-sampling routine corresponding to the parameters quoted in Table  3 . The couplings g 𝑔 g quoted in parentheses are defined as residua at the complex pole positions at s ∗ = m σ 2 superscript 𝑠 superscript subscript 𝑚 𝜎 2 s^{*}=m_{\sigma}^{2} or s ∗ = m ρ 2 superscript 𝑠 superscript subscript 𝑚 𝜌 2 s^{*}=m_{\rho}^{2} via g 2 = lim s → s ∗ | ( s − s ∗ ) ​ T ​ ( s ) | superscript 𝑔 2 subscript → 𝑠 superscript 𝑠 𝑠 superscript 𝑠 𝑇 𝑠 g^{2}=\lim_{s\to s^{*}}|(s-s^{*})T(s)| .

$m_{\pi }$ [MeV]	$\sim 315$	$\sim 224$	$139$
$m_{\pi }{a}_0^{I=0}$	$+{1.9008}_{-0.0593}^{+0.0521}$	$+{0.6985}_{-0.0015}^{+0.0010}$	$+{0.2132}_{-0.0009}^{+0.0008}$
$m_{\pi }{a}_0^{I=2}$	$-{0.1538}_{-0.0018}^{+0.0021}$	$-{0.0952}_{-0.0009}^{+0.0010}$	$-{0.0433}_{-0.0002}^{+0.0002}$
$m_{\sigma }$ [MeV]	$+{591}_{-5}^{+6}-i{109}_{-4}^{+4}$	$+{502}_{-4}^{+4}-i{175}_{-5}^{+6}$	$+{443}_{-3}^{+3}-i{221}_{-6}^{+6}$
$g_{\sigma \pi \pi }$ [MeV]	${533}_{-2}^{+2}$	${426}_{-2}^{+2}$	${397.8}_{-0.6}^{+0.6}$
$m_{\rho }$ [MeV]	$+{789}_{-1}^{+1}-i{20}_{-0}^{+0}$	$+{738}_{-1}^{+2}-i{43}_{-1}^{+1}$	$+{724}_{-4}^{+2}-i{67}_{-1}^{+1}$
$g_{\rho \pi \pi }$ [MeV]	${226}_{-2}^{+2}$	${282}_{-2}^{+3}$	${323}_{-3}^{+5}$

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].$$

$$⟨P(t)P(0{)}^{†}⟩\stackrel{t\to \infty }{\to }\frac{Z^2}{2m_{\pi }}{e}^{-m_{\pi }t}.$$

$$m_{\text{PCAC}}\equiv \frac{1}{2}\frac{⟨{\partial }_t{A}_4(t)P(0{)}^{†}⟩}{⟨P(t)P(0{)}^{†}⟩},$$

$$f_{\pi }\equiv \sqrt{2Z^2}\frac{m_{\text{PCAC}}}{m_{\pi }^2}.$$

$$T_{mIAM}^{Il}(s)=\frac{(T_2^{Il}(s){)}^2}{T_2^{Il}(s)-T_4^{Il}(s)+A_m^{Il}(s)}.$$

$$A_m^{Il}(s)$$

$$-\frac{(s_2-s_A)(s-s_2)}{s-s_A}{\left(T_2^{Il}(s_2)-T_4^{Il}(s_2)\right)}^{\prime },$$

$$m_{\pi }^2=m^2\left(1-\frac{m^2}{32{\pi }^2{f}_0^2}{\bar{l}}_3\right),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)$$

$$T_{\mathrm{mIAM}}^{Il}(s)=\frac{\sqrt{s}}{2p(\cot {\delta }_{\mathrm{mIAM}}^{Il}(s)-i)}.$$

$$\cot {\delta }_{\mathrm{mIAM}}^{Il}(s)$$

$$\frac{\sqrt{s}}{2p}$$

$$p\cot {\delta }^{00}(s)$$

$$\cot {\delta }^{11}(s)$$

$$\cot {\delta }^{20}(s)$$

$${\chi }^2\equiv \sum _{i=1}^6(X_i-Y_i{)}^T\cdot {\mathrm{Cov}}_{{\mathcal{E}}_i}^{-1}\cdot (X_i-Y_i)$$

Full text

I Introduction

Lattice QCD calculations provide an ab-initio access to particle scattering subject to the strong interaction. Calculations of two-hadron states are performed in a small cubic volume with periodic boundary conditions and they are connected analytically to infinite-volume scattering amplitudes via mapping established in Refs. Luscher (1986); Lüscher (1986, 1991). In this context, pion-pion scattering has been a prime subject for lattice QCD calculations, in isospin $I=2$ Sharpe et al. (1992); Kuramashi et al. (1993); Gupta et al. (1993); Yamazaki et al. (2004); Aoki et al. (2005); Beane et al. (2006, 2008); Feng et al. (2010); Yagi et al. (2011); Fu (2013); Sasaki et al. (2014); Kawai et al. (2018); Dudek et al. (2012); Bulava et al. (2016); Akahoshi et al. (2019); Helmes et al. (2015, 2019); Hörz and Hanlon (2019), $I=1$ Aoki et al. (2007); Feng et al. (2011); Göckeler et al. (2008); Lang et al. (2011); Bali et al. (2016); Guo et al. (2016); Pelissier and Alexandru (2013); Aoki et al. (2011); Dudek et al. (2013); Feng et al. (2015); Metivet (2015); Wilson et al. (2015); Alexandrou et al. (2017); Andersen et al. (2019); Fu and Wang (2016); Alexandrou et al. (2017); Werner et al. (2019), and $I=0$ Briceño et al. (2017, 2018); Guo et al. (2018); Fu and Chen (2018); Liu et al. (2017).

Since most lattice QCD calculations are carried out, for technical reasons, using quark masses heavier than the physical ones, extrapolation of lattice QCD results to the physical point requires a model for energy and pion mass dependence. Chiral Perturbation Theory (ChPT) allows for the controlled expansion of the $\pi\pi$ scattering amplitude in quark masses and meson momenta. The ChPT expansion is valid only in the non-resonant region around the threshold. Resonances such as $\rho(770)$ and $f_{0}(500)$ (or ‘$\sigma$’), but also the high energy behavior of the scattering amplitudes, require a non-perturbative treatment, usually guided by imposing constraints from analyticity and two-body unitarity. The inverse amplitude method (IAM) Truong (1988); Pelaez and Rios (2006); Gómez Nicola et al. (2008); Pelaez and Rios (2010) provides a well-established framework in this sense. It is unitary and matches the chiral pion-pion amplitude Gasser and Leutwyler (1985, 1984) up to the next-to-leading order (NLO).

The IAM describes different isospin channels within the same approach; the fit parameters in form of low-energy constants, pion masses, and pion decay constants enter different channels and connect them. Here, for the first time we apply this method to lattice QCD data in different isospin channels determined from the same ensembles Guo et al. (2018, 2016); Culver et al. (2019), thus minimizing the effect of data inconsistencies.

In comparison to previous studies in which the IAM 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 two-hadron energy levels (ELs) in the fit; for the isospin $I=2$ channel alone, this was recently achieved in Ref. Culver et al. (2019). In addition, the correlations of ELs between different isospin channels are taken into account, for the first time. We work in lattice units throughout the paper up to the last step, when evaluating the amplitude at the physical pion mass, which involves a scale setting.

The global fit allows one to address issues that remained unresolved in previous studies. For example, using SU(3) unitarized ChPT with contact terms Oller et al. (1999), which can only be approximately matched to NLO SU(3) CHPT, it was found that data from $N_{f}=2$ lattice QCD calculations by different groups Bali et al. (2016); Guo et al. (2016); Göckeler et al. (2008); Lang et al. (2011) extrapolate to consistently low physical $\rho$ masses Hu et al. (2016) (see Refs. Gomez Nicola and Pelaez (2002); Nebreda and Peláez. (2010); Guo et al. (2012, 2015) for more recent one-loop SU(3) versions). Within the model, the $K\bar{K}$ channel could explain the discrepancy, but also missing higher orders in the model could be responsible, or systematic effects on the side of lattice QCD calculations: Small volumes could induce large exponential corrections that cannot be systematically represented in models for the scattering region. Scale setting can be another source of discrepancy in lattice QCD data. See also Ref. Hu et al. (2017) for a corresponding analysis of $N_{f}=2+1$ lattice QCD calculations. Here, we have the opportunity to analyze different isospin channels at the same time with a model that matches to NLO two-flavor ChPT. This will allow us to avoid some of the potential pitfalls at least for the results computed by the GW lattice QCD group Guo et al. (2018, 2016); Culver et al. (2019).

II Lattice QCD details

The input used in this study are two-hadron state energies in a finite box with periodic boundary conditions computed using lattice QCD. Our results use two different quark masses corresponding to a pion mass of $m_{\pi}=224\mathop{\hbox{MeV}}$ and $m_{\pi}=315\mathop{\hbox{MeV}}$. For each quark mass we use three different box geometries to scan the energy region below the inelastic threshold in each channel. The parameters for each of the six lattice ensembles are listed in Table 1.

All lattice results are computed using QCD with two mass-degenerate quark flavors ($N_{f}=2$ QCD). This is a good approximation of the real world when focusing on the lightest hadrons which are composed mainly of up and down quarks, as is the case in this study. Moreover, $N_{f}=2$ QCD is also a very interesting theoretical model with a minimal number of parameters: one quark mass and $\Lambda_{\text{QCD}}$. This can be used as a precise testbed for non-perturbative aspects of QCD. To generate the ensembles, both gauge and quark actions use improved discretizations. For the gauge action we use Lüscher-Wise action Luscher and Weisz (1985) and for the quarks we use nHYP discretization Hasenfratz et al. (2007).

The data analysed here correspond to two-pion states for all three possible isospin combinations. Specific details can be found in Refs. Pelissier and Alexandru (2013); Guo et al. (2018) for $I=1$, Ref. Guo et al. (2016) for $I=0$, and Ref. Culver et al. (2019) for $I=2$. Here we review briefly the lattice methods used to compute the two-pion state energies and the other relevant observables used in this study.

II.1 Two-pion finite volume spectrum

The spectrum of two-hadron states in a box with periodic boundary conditions is quantized. The energy levels, in the elastic region, can be related with the scattering amplitude in the infinite volume. To extract the energy levels we use the standard variational method Luscher and Wolff (1990). For each isospin channel we construct a set of interpolating fields that are expected to have large overlap with the lowest lying states in the spectrum. We include sufficient interpolators to resolve all the states with energies below the inelastic threshold. The interpolators $\mathcal{O}_{i}$ are used to construct a correlator matrix,

[TABLE]

The eigenvalues are extracted by solving the generalized eigenvalue problem,

[TABLE]

Here $t_{0}$ is a parameter that is adjusted for each isospin to help dampen the effects of excited state contributions. It was shown in Refs. Luscher and Wolff (1990); Blossier et al. (2009), that the energies of the system can be extracted from the long-time behavior of the generalized eigenvalues

[TABLE]

Finite volume states will inherit the symmetries of the box and they can be classified according to the irreducible representations (irreps) of the box symmetry group. To properly identify the states corresponding to the energy levels extracted from the variational analysis, we need to design interpolating fields that have the appropriate transformation properties. For this study we use both cubic boxes and boxes that are elongated in one of the dimensions. The relevant lattice symmetry group for the cubic box is $O_{h}$ whose 10 irreps are conventionally named as $A^{\pm}_{1},A^{\pm}_{2},E^{\pm},T^{\pm}_{1},T^{\pm}_{2}$, and for the elongated box $D_{4h}$ whose 10 irreps are $A^{\pm}_{1},A^{\pm}_{2},E^{\pm},B^{\pm}_{1},B^{\pm}_{2}$. In the energy range we study, scattering in the $I=0$ and $I=2$ is dominated by the $\ell=0$ partial wave, and for $I=1$ the $\ell=1$ partial wave. In a finite box, different partial waves are mixed by the finite volume effects. The rotationally-symmetric SO(3) angular momentum multiplets in the continuum are split into multiplets transforming under the symmetry group of the box. The splitting is shown in Table 2 (for details see Ref. Lee and Alexandru (2017)). We see that the relevant lattice irreps for $I=0,2$ channels are $A_{1}^{+}$ in both box symmetries. For the $I=1$ channel we use the $A_{2}^{-}$ irrep which is sensitive to the elongation.

Once the appropriate operators are determined, correlation functions are computed using Wick contractions on the quark fields. The resulting quark-diagrams depend on the isospin channel – the details are included in the references listed above. The correlation functions require the evaluation of the all-to-all quark propagator, that is the quark propagator from every point on the lattice to any other point. To avoid the full expense of this calculation we use the Laplacian-Heaviside method (LapH) Peardon et al. (2009). The idea is to truncate the quark interpolating fields by dropping out the high-frequency modes of the three-dimensional Laplacian on each time slice while preserving the symmetry of the resulting “smeared” fields. The interpolating fields constructed out of these quarks excite the same QCD states, but overlap better with the low-energy states. An additional advantage is that we only need to invert the Dirac matrix for the LapH modes, reducing the numerical cost. For all isospin channels we used the $N_{v}=100$ lowest eigenvectors corresponding to a smearing radius of approximately $0.5\mathop{\hbox{fm}}$ Guo et al. (2016). The smeared quark propagators were computed efficiently using a set of GPU inverters Alexandru et al. (2012).

II.2 Other observables

Besides the energy of the two-pion states, the other lattice QCD inputs for the analysis are the pion mass and the pion decay constants. These parameters for each lattice ensemble are listed in Table 1. We note that all ensembles are generated with the same coupling, which should generate the same lattice spacing (or cutoff). For ensemble $\mathcal{E}_{1,2,3}$ the quark mass is the same and similarly for ensembles $\mathcal{E}_{4,5,6}$. The differences in pion mass within these sets is thus expected to be the result of statistical fluctuation and/or finite volume effects. The pion mass was computed by evaluating the two-point function of the pion using LapH with a $\bar{u}\gamma_{5}d$ interpolating field.

The value of the pion decay constant $f_{\pi}$ was computed using standard methods (see for example Ref. Brandt et al. (2013)). We use two two-point correlation functions: $\left\langle A_{4}(t)P(0)^{\dagger}\right\rangle$ and $\left\langle P(t)P(0)^{\dagger}\right\rangle$, where $P$ is the pseudo-scalar density $\bar{q}\gamma_{5}q$ and $A_{4}$ is the axial current density $\bar{q}\gamma_{4}\gamma_{5}q$. Both $P(t)$ and $A_{4}(t)$ are projected to zero spatial momentum. From the $\left\langle P(t)P(0)^{\dagger}\right\rangle$ correlation function the pion mass $m_{\pi}$ and overlap factor $Z$ are extracted:

[TABLE]

The ratio of the two correlation functions is used to calculate the current quark mass

[TABLE]

where $\partial_{t}A_{4}(t)\equiv[A_{4}(t+1)-A_{4}(t-1)]/2$. Using these values the pion decay constant is defined as

[TABLE]

The decay constant needs to be renormalized, but the renormalization is expected to introduce only a couple of percent shift Hoffmann et al. (2007).

For each pion mass we use one cubic box, and two elongated boxes. The elongated boxes are important for getting a good scan of the relevant scattering region in each isospin channel. The lowest momentum of a particle moving in the elongated direction is $\frac{2\pi}{L\eta}$. The energy of the particle thus changes as we vary $\eta$. For two-pion states with total momentum $\mathbf{P}=[000]$ the energy changes little with $\eta$ unless the state corresponds to two pions with non-zero relative momentum, which tends to be the case for the excited levels in each channel. For states with total momentum $\mathbf{P}=[100]$ in the elongated direction, even for the lowest state the energy usually varies with $\eta$ and the energy levels cover well the kinematic region of interest. This is indeed observed in the extracted spectrum shown in Figure 1.

III Global study of $\pi\pi$ scattering

The determination of the $I=2$ finite volume spectrum from lattice QCD calculations concludes a multi-year program Guo et al. (2018, 2016); Culver et al. (2019) of the GW lattice group in obtaining comprehensive information on $\pi\pi$ scattering from first principles. The obtained set of energy eigenvalues covers a large energy region from below the production threshold to and beyond the resonance region. At the same time these results describe $\pi\pi$ scattering at two unphysical pion masses ($\sim 1.5$ and $\sim 2.5~{}m_{\pi}^{\rm phys}$). Thus, they provide a unique opportunity for mapping out the $m_{\pi}$ vs. $E$ plane with respect to $\pi\pi$ interactions in all three isospins, which is explored in the following.

To make full use of the available information, a scattering amplitude is required which not only takes into account the analytic properties in $E$ but also the chiral behavior, consistent with constraints from perturbative ChPT at NLO as well as those from Ref. Bruns and Mai (2018). A method reconciling both demands is the so-called modified inverse amplitude method (mIAM) Truong (1988); Dobado and Pelaez (1997); Gómez Nicola et al. (2008, 2008). In the past, it has been shown to be very successful in describing experimental data on $\pi\pi$ scattering Hanhart et al. (2008); Pelaez (2016); Nebreda and Peláez. (2010) in all three isospin channels, while also having correct chiral behavior up to the next-to-leading chiral (NLO) order. Likewise it fulfills the general requirements on the chiral trajectory for resonances, derived in Ref. Bruns and Mai (2018) to all chiral orders.

III.1 Infinite volume spectrum

The modified inverse amplitude method is based on the leading ($T_{2}^{Il}(s)$) and NLO ($T_{4}^{Il}(s)$) chiral amplitudes projected to a specific isospin ($I$) and angular momentum ($l$). A unitary scattering amplitude $T_{mIAM}^{Il}(s)$ can then be derived Truong (1988) using dispersion relations, namely

[TABLE]

The term $A^{Il}_{m}(s)$ in the denominator does not arise for dynamical reasons, but has been introduced Fernandez-Fraile et al. (2007); Gómez Nicola et al. (2008) to avoid the appearance of an unphysical pole. Such a pole is associated with appearance of the so-called Adler zero – a sub-threshold energy at which the amplitude vanishes as demanded by chiral symmetry. Explicitly it reads

[TABLE]

where $s_{A}$ and $s_{2}$ are the zeros of $T_{2}(s)-T_{4}(s)$ and $T_{2}(s)$, respectively. With this the $T_{mIAM}$ is an analytic, unitary scattering amplitude, which indeed reproduces the usual chiral expansion and is crossing symmetric up to the next-to-leading chiral order.

The leading order chiral amplitude 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}^{Il}$ 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 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$. Hence, for a fixed scale $\mu$ one can determine the renormalized LECs and then make predictions for two-particle scattering at a different pion mass. In the course of this work we use dimensional regularization with $\mu=770$ MeV, but emphasize that the amplitude (7) is manifestly scale independent.

III.2 Finite volume spectrum

The scattering amplitude introduced in the previous section describes the scattering of two pions in infinite volume in terms of a complex-valued function of continuous energy/momentum variables. In finite volume, momenta are discretized leading to a discretized interaction spectrum. The way to convert the latter into phase-shifts is given by Lüscher’s method Lüscher (1986); Luscher (1986), see also Refs. Lee and Alexandru (2017); Briceno et al. (2018). In the context of mIAM its implementation has been used and is described in Refs. Mai and Döring (2019); Mai (2018).

In infinite volume, the scattering amplitude (7) is related to the phase-shifts via

[TABLE]

We use Eq. (7) to compute the phase-shifts predicted by mIAM:

[TABLE]

where $\bar{T}_{4}^{Il}$ denotes the NLO chiral amplitude without $s$-channel loop diagrams and $J(s)$ denotes the meson-meson loop in dimensional regularization for $\mu=770$ MeV. The determination of the corresponding finite-volume spectrum amounts to finding the roots of the following set of equations

[TABLE]

for the irreps $A_{1}^{+}$, $A_{2}^{-}$, and $A_{1}^{+}$, respectively. In every interaction channel ($I=0,1,2$ from top to bottom, respectively) the right-hand side carries the information on the geometry of the lattice (size $L$, elongation $\eta$), and kinematics via $q=L/(2\pi)p(s)$, where $p$ is the magnitude of the pion relative three-momentum in the center of the mass system. The required Lüscher functions for elongated boxes as well as corresponding formulas for moving frames are quoted in Ref. Lee and Alexandru (2017). The left hand side of the above equations contains only the information on the pion interaction in the corresponding channel.

We note that Eqs. (12) are valid when neglecting higher partial waves and only below the inelastic threshold (4$m_{\pi}$). The cutoff in the angular momentum space is justified by the smallness of higher partial waves in all channels. Additionally, note that the factor $p$ on both sides of the first equation makes it well-defined below $\pi\pi$ threshold, where the lattice QCD result also exist, see also previous studies Döring et al. (2018); Guo et al. (2018); Briceño et al. (2017) of this channel.

III.3 Fit to finite volume spectrum

As discussed above there are four parameters in the model – the LECs $l_{i}^{r}$. Multiple work-flows are possible when confronting lattice QCD results with the model, i.e. with respect of how many interaction channels are included, or which of the experimental or lattice data are included. After evaluating several such options, we decided to fit to lattice QCD results only and predict $\pi\pi$ scattering at the physical point. This is possible due to the correct chiral properties of the mIAM also in the vicinity of resonances Bruns and Mai (2018). Additionally, the fitted lattice QCD eigenvalues include correlations111The covariance matrices as well as the energy eigenvalues are collected for all channels for convenience at https://github.com/chrisculver/Pipi_Energies_Covariances., which as we found, make considerable contribution to the total $\chi^{2}$.

Overall, the data consist of 95 energy eigenvalues in the $I=0,1,2$ channels, extracted from 6 ensembles as in Table 1 and depicted in Figure 1. For each of the ensembles the $m_{\pi}$ and $f_{\pi}$ have been recorded including the corresponding correlations. Varying the pion mass can lead to non-negligible effects as noted in Ref. Bulava et al. (2016). In relation to the recorded pion mass we noted a systematic effects in $\mathcal{E}_{2}$, i.e., its different central value in Table 1. Thus, we exclude this ensemble from further fits. This leaves us with 88 data points and correlations between different channels, energy levels, pion masses and decay constants within each ensemble. In exploratory fits we have found that the value of $l^{r}_{3}$ does not lead to any notable improvement of the $\chi^{2}$ so that we fix it to the value reported by FLAG Aoki et al. (2019), i.e., $l^{r}_{3}=8.94\cdot 10^{-6}$. This leaves us with three fit parameters ($l^{r}_{1}$, $l^{r}_{2}$, $l^{r}_{4}$) of the model as well two pairs $(m_{\pi},f_{\pi})_{\rm heavy/light}$ for heavy and light ensembles, respectively. We allow the pertinent four parameters to vary instead of using their central values due to the high precision of the data especially in the threshold proximity. In each fit scenario the correlated $\chi^{2}$,

[TABLE]

is minimized with respect to the seven fit parameters. Here $X_{i}$ denotes the vector containing central values of $m_{\pi}$, $f_{\pi}$ and energy eigenvalues of a given lattice ensemble. The vector $Y_{i}$ contains the corresponding result obtained from the model as in Eq. (12). The covariance matrix ${\rm Cov}_{\mathcal{E}_{i}}$ contains uncertainties and correlations of all data, including correlations between different isospin channels.

III.4 Fit results

The best fit gives $\chi^{2}_{\rm dof}=218/(88-7)~{}(\approx 2.7)$ with the parameters collected in Table 3. We explored various other fit scenarios to determine the root of the relatively large $\chi^{2}$ value. For instance, we note that a sizable contribution to the total $\chi^{2}$ is generated by the correlations within each ensemble across different isospin channels (up to 50%). Furthermore, we have performed separate fits for light and heavy ensembles only. The best fit results of those are $\chi^{2}_{\rm light}=63.8/(47-5)\approx 1.5$ and $\chi^{2}_{\rm heavy}=107.2/(43-5)\approx 2.8$. The tension between the heavy and light mass fits, noted in Ref. Culver et al. (2019) (i.e. $\chi^{2}_{{\rm light\&heavy}}/(\chi^{2}_{{\rm light}}+\chi^{2}_{{\rm heavy}})\approx 2$), is reduced when all channels are considered. Moreover, the fact of $\chi^{2}_{\rm heavy}>\chi^{2}_{\rm light}$ is easy to understand since the model used relies on the chiral amplitudes up to the fourth chiral order. Thus, it is not too surprising that the corresponding description begins to fail at the heavier pion mass. In this respect, see the discussion in Ref. Dürr (2015). Also, as depicted in Figure 1, large contributions to the $\chi^{2}$ come from energy eigenvalues at high energies in $I=1$ and $I=2$ channels, where the error bars are especially narrow. A similar observation was also made in Ref. Guo et al. (2016) dealing with the $I=1$ channel only, supported by two simpler (Breit-Wigner and Lippmann-Schwinger type) models.

The probability distribution function of the low-energy constants is depicted in Figure 2 together with that of the re-sampled $\chi^{2}_{\rm res}$ in the inset. The latter should follow a non-central $\chi^{2}$ distribution function with non-centrality parameter $\lambda=\chi^{2}=218$ of the fit to the original data and $k=88-7$ degrees of freedom, peaking at a re-sampled $\chi_{\rm res}^{2}$ of around $\chi_{\rm res}^{2}=300$ shown as the blue curve overlaying the histogram in the inset of Figure 2. The difference between histogram and theoretical expectation for the $\chi_{\rm res}^{2}$ distribution is explained by the mentioned shortcoming of the fit at higher energies / pion mass and across different isospin channels. The errors on LECs are on the order of 5-10%, which is an order of magnitude smaller than those from the fit of the same model to one channel ($I=2$) only Culver et al. (2019). This shows that inclusion of all isospin channels indeed restricts the model and thus the extraction of physically relevant information strongly.

Finally, the results of the fit in all three channels at two pion masses are also presented in terms of infinite-volume quantities in the lower part of Figure 1. There the energy eigenvalues obtained from a lattice calculation are also mapped to the phase-shifts using Lüscher’s method. In terms of the phase-shifts some of the results produce error bars wrapping through the whole codomain of this mapping. This occurs when the error on energy eigenvalue overlaps with a non-interacting value.

III.5 Predictions at the physical point

The parameters of the model $(l^{r}_{1},l^{r}_{2},l^{r}_{4})$ have been determined in a fit to the lattice results in a model with correct (up to second chiral order) pion mass dependence. This allows us to extrapolate the amplitude to the physical point to confront the pertinent phase-shifts with the phenomenological results. The extrapolations are shown in Figure 3 together with the phase-shifts extracted from experiment. The bands show the $1\sigma$ region, originating from the re-sampling of the fits. To emphasize that unitarity is strictly fulfilled only up to the first inelastic threshold ($4m_{\pi}$), we use a different color for the predicted curves above this region.

We observe that in the even isospin channels the prediction agrees with the experimental data very well in the elastic region and even beyond it. In the $I=1$ channel the functional behaviour is very similar to the one suggested by experiment, but is shifted to the left. This suggests a lighter mass of the $\rho$ resonance and we will return to this discussion point below.

The physical parameters, such as scattering lengths, resonance pole positions and couplings have been discussed in length in the previous papers, dealing with single channels Hu et al. (2016); Guo et al. (2018); Hu et al. (2017, 2016); Guo et al. (2016). There the dependence on the utilized model has been discussed using a broad class of chiral unitary models, Breit-Wigner type, and models based on conformal mapping. Given that the mIAM is a better compromise between constraints from chiral symmetry Bruns and Mai (2018) and analytic properties of scattering amplitudes we simplify the discussion here by discussing the results of this approach only. The collection of observables at physical and both unphysical pion masses can be found in Table 4.

The sizes of the even-isospin scattering lengths are slightly smaller than the phenomenological values $m_{\pi}\,a_{0}^{I=0}=0.2198(46)_{\rm stat}(16)_{\rm syst}(64)_{\rm th}$ and $m_{\pi}\,a_{0}^{I=2}=-0.0445(11)_{\rm stat}(4)_{\rm syst}(8)_{\rm th}$ of Ref. Caprini et al. (2012). Interestingly, most lattice QCD based determinations of the latter tend to be smaller in magnitude than the results based on experimental results, see the discussion in the FLAG report Aoki et al. (2019).

The extrapolated pole position of the isoscalar resonance agrees well in the real part but is too small in the imaginary part when compared to the phenomenologically driven analysis of Ref. Pelaez (2016): $m_{\sigma}=449_{-16}^{+22}-i275_{-12}^{+12}$ MeV which is an average of analyses based on Roy equations and related methods Colangelo et al. (2001); Caprini et al. (2006); Moussallam (2011); Garcia-Martin et al. (2011). For the chiral extrapolations of the HadronSpectrum data Briceño et al. (2017) in Ref. Döring et al. (2018) and the GW lattice QCD data Guo et al. (2018) similarly narrow $\sigma$ pole positions were found; also, in fits to experimental data only, e.g., in Refs. Oller (2003); Dobado and Pelaez (1997); Döring and Meißner (2012) based on the mIAM or unitarized ChPT with contact terms only, the $\sigma$ was rather narrow compared to the averaged value of Ref. Pelaez (2016) from Roy equations. These observations suggest that methods based only on s-channel unitarization (with up to one loop in the $t,u$-channels) tend to produce slightly narrower $\sigma$ resonances compared to the Roy equations which provide better sub-threshold analytic properties.

The extrapolated result on the isovector resonance agrees well in its width with the phenomenological value Tanabashi et al. (2018) of $\Gamma\approx 150$ MeV. However, its mass is too small by $\sim 40$ MeV. This corroborates a similar finding for the chiral extrapolation of several $N_{f}=2$ lattice QCD calculations in Ref. Hu et al. (2016) based on a simpler model and using only data from the $I=1$ channel. In this model only NLO contact interactions were used in the unitarization, based on Ref. Oller et al. (1999). As shown here, the light $\rho$ is not a consequence of this simplification.

Another potential reason for the discrepancy between the physical $\rho$ mass and our extrapolation can be attributed to the ambiguity in determining the lattice spacing. We argue here that this issue is a bit more subtle: the lattice spacing is not directly relevant while the definition of the physical point is. To see this note that all fits are carried using inputs in lattice units: energy levels, pion masses and decay constants, and their correlations are all dimensionless. The lattice data together with the mIAM predicts the phase-shifts $\delta^{IJ}(l_{i}^{r},E/m_{\pi},f_{\pi}/m_{\pi},\mu/m_{\pi})$, all parameters being dimensionless ratios. The LEC’s $l_{i}^{r}$ are fixed by the fit to the dimensionless lattice data. The dependence on $\mu/m_{\pi}$ is very small and can be disregarded. The prediction of $\delta^{IJ}$ versus $E/m_{\pi}$ is completely independent from the lattice spacing $a$. The only relevant parameter is $f_{\pi}/m_{\pi}$, which needs to be set to the value corresponding to the physical point. For the results in Fig. 3 we define the physical point by setting $f_{\pi}/m_{\pi}=92.4/139$, the physical values for $\pi^{+}$. Note that there is an ambiguity in defining this point for $N_{f}=2$ simulations due to the absence of the strange quark and isospin breaking effects.

Finally, we have also checked the effects on the pion mass due to the elongation of the box. We could not detect any systematic shift of the pion mass when allowing for different fit parameters for the pion masses in the different spatial elongations. In particular, the pertinent independent-mass fit leads to a very similar chiral extrapolation. Other possible explanations for the discrepancy could be related to the missing $K\bar{K}$ channel Hu et al. (2016) or higher orders effects in the chiral expansion Niehus et al. (2019).

IV Summary and conclusions

We study pion-pion elastic scattering across all three isospin channels in two-flavor dynamical lattice QCD. It is the first time such a cross-channel analysis of lattice QCD scattering data is attempted. We use elongated lattices which offer a cost-effective alternative to cubic lattices in mapping out the momentum dependence in scattering processes. We consider six ensembles with elongations up to a factor of two in one of the spatial dimensions, and two quark masses corresponding to pion masses at 315 MeV and 224 MeV. The lattice input include two-pion states at rest and also states boosted along the elongated direction for enhanced momentum coverage.

To make contact with phenomenology, we put the finite-volume spectrum across all three channels through a correlated global analysis using the inverse amplitude method. This method unitarizes chiral perturbation theory in the s-channel, matching the chiral $\pi\pi$ scattering amplitude to next-to-leading order, and it allows for chiral extrapolation of lattice data to the physical point over a wide energy range. We treat the pion mass and decay constant as fit parameters and include their full correlations with the energy eigenvalues in the corresponding analysis; likewise, full correlations across different isospin channels are included.

The model qualitatively captures the scattering phase-shifts in the entire elastic region. However, the $\chi^{2}$ indicates a tension between the model and the lattice QCD data. Remarkably, a large contribution to the $\chi^{2}$ originates from correlations between different isospin channels, which have been taken into account for the first time in the present study. Quantitatively the agreement is better for the lower quark mass and lower energies as expected from a ChPT-based model. The cross-channel fit provides better constraints on the model parameters and this in turn leads to a tighter extrapolation to the physical point. The extrapolation agrees well with the experimental phase-shifts in the $I=0$ and $I=2$ channel. In the $I=1$ channel, the extrapolations favor a lower mass for the $\rho$ resonance, as was found in previous studies. We argued here that this disagreement is not related to the determination of lattice spacing, since this has almost no effect on the determination of the model parameters. Among other possible explanations for the discrepancy, there is, however, an ambiguity in defining the physical point since the experimental data is not directly comparable with $N_{f}=2$ simulations.

Overall our results demonstrate that elongated lattices combined with a global analysis in the inverse amplitude method can be an effective tool in probing hadron-hadron scattering processes from first principles. We plan to extend the approach to more systems, such as three pions above threshold, and pion-baryon scattering.

Acknowledgements.

We thank Dehua Guo for helping us compute the correlation between the energies in the $I=2$ channel and the $I=0,1$ channels. 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. MM and MD are supported in part by the National Science Foundation CAREER grant PHY-1452055 and by U.S. DOE Grant No. DE-AC05-06OR23177. 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.