for the Hadron Spectrum Collaboration

#  The energy-dependent $\pi^{+}\pi^{+}\pi^{+}$ scattering amplitude from QCD

Maxwell T. Hansen<EMAIL_ADDRESS>Theoretical Physics Department,
CERN, 1211 Geneva 23, Switzerland  Raul A. Briceño<EMAIL_ADDRESS>Thomas
Jefferson National Accelerator Facility, 12000 Jefferson Avenue, Newport News,
VA 23606, USA  Department of Physics, Old Dominion University, Norfolk, VA
23529, USA  Robert G. Edwards<EMAIL_ADDRESS>Thomas Jefferson National
Accelerator Facility, 12000 Jefferson Avenue, Newport News, VA 23606, USA
Christopher E. Thomas<EMAIL_ADDRESS>DAMTP, University of
Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK  David J. Wilson
<EMAIL_ADDRESS>DAMTP, University of Cambridge, Wilberforce Road,
Cambridge, CB3 0WA, UK

###### Abstract

Focusing on three-pion states with maximal isospin ($\pi^{+}\pi^{+}\pi^{+}$),
we present the first non-perturbative determination of an energy-dependent
three-hadron scattering amplitude from first-principles QCD. The calculation
combines finite-volume three-hadron energies, extracted using numerical
lattice QCD, with a relativistic finite-volume formalism, required to
interpret the results. To fully implement the latter, we also solve integral
equations that relate an intermediate three-body K matrix to the physical
three-hadron scattering amplitude. The resulting amplitude shows rich analytic
structure and a complicated dependence on the two-pion invariant masses,
represented here via Dalitz-like plots of the scattering rate.

††preprint: CERN-TH-2020-147, JLAB-THY-20-3242

_Introduction_ — The three-body problem lies at the core of a broad range of
outstanding questions in quantum chromodynamics (QCD). The largest uncertainty
in QCD-based structure calculations of light nuclei, for example, is the
estimate of the three-nucleon force (see Ref. Piarulli  _et al._ (2018)). In
addition, many QCD resonances have significant branching fraction to channels
with three or more hadrons. The Roper resonance, for example, has defied
simple quark-model descriptions, due in part to its nature as a broad
resonance with a $\sim 30\%$ branching fraction to $N\pi\pi$. A rigorous QCD
calculation would elucidate the role of non-perturbative dynamics in the
Roper’s peculiar properties, e.g. the fact that it has a lower mass than the
negative-parity ground state, which seems unnatural from the perspective of
the quark model Isgur and Karl (1977, 1979).

As a necessary step towards studying a broad class of three-hadron systems, in
this work we present the first study of an energy-dependent three-body
scattering amplitude from QCD. This non-perturbative result is achieved by the
coalescence of three novel techniques: a calculation of finite-volume three-
hadron energies based in numerical lattice QCD, a relativistic finite-volume
formalism to relate the energies to K matrices, and a numerical evaluation of
corresponding integral equations to convert the latter into the three-hadron
scattering amplitude. The theoretical basis required to achieve these final
two steps was derived in Refs. Hansen and Sharpe (2014, 2015).111A large body
of work has investigated general methods for relating finite-volume energies
to scattering amplitudes for both two- and three-body states. See Refs.
Lüscher (1991); Rummukainen and Gottlieb (1995); Bedaque (2004); Kim  _et al._
(2005); Fu (2012); Leskovec and Prelovsek (2012); Göckeler  _et al._ (2012);
He  _et al._ (2005); Lage  _et al._ (2009); Bernard  _et al._ (2011); Döring
_et al._ (2011); Döring and Mei$\ss$ner (2012); Agadjanov  _et al._ (2014);
Döring  _et al._ (2012); Hansen and Sharpe (2012); Briceño and Davoudi
(2013a); Guo  _et al._ (2013); Briceño  _et al._ (2018a); Briceño (2014) and
Refs. Romero-López  _et al._ (2019); Briceño  _et al._ (2019a, 2017a); Hansen
and Sharpe (2019); Briceño  _et al._ (2018b); Blanton  _et al._ (2019); Hammer
_et al._ (2017a, b); Briceño and Davoudi (2013b); Polejaeva and Rusetsky
(2012); Mai and Döring (2017); Guo and Gasparian (2017, 2018); Guo and Morris
(2019); Romero-López  _et al._ (2018); Jackura  _et al._ (2019a); Briceño  _et
al._ (2019b); Hansen  _et al._ (2020); Blanton and Sharpe (2020a, b); Beane
_et al._ (2020), respectively.

This work considers the scattering of three-pion states with maximal isospin
$(I=3)$ in QCD with three dynamical quarks ($N_{f}=2+1$): two degenerate light
quarks, with heavier-than-physical mass corresponding to a pion mass
$m_{\pi}\approx 391\,\text{MeV}$, and a strange quark. This channel offers an
optimal benchmark case, since both the maximal-isospin three-pion system and
its two-pion subsystem are expected to be weakly interacting and non-resonant.

Many numerical studies of three-hadron states have been published over the
last decade, ranging from early work deriving and fitting large-volume
expansions of the three-pion ground state Beane  _et al._ (2008); Detmold  _et
al._ (2008, 2012) to more recent results using quantization conditions to
study ground Mai and Döring (2019) and excited states Blanton  _et al._
(2020); Mai  _et al._ (2020); Guo and Long (2020), with the latter set each
analyzing the lattice QCD spectrum published in Ref. Hörz and Hanlon (2019).
Independent sets of finite-volume energies have also been calculated and
analyzed in Refs. Culver  _et al._ (2020) and Fischer  _et al._ (2020). The
present investigation goes beyond this previous work, by providing the first
complete numerical determination of physical scattering amplitudes for three-
body systems.

In the following, we first discuss our determination of two- and three-pion
finite-volume energies, before describing the fits used to relate these to
infinite-volume K matrices. The latter then serve as inputs to known integral
equations, which we solve numerically to extract the $3\pi^{+}\to 3\pi^{+}$
scattering amplitude. Additional details of the analysis are discussed in the
supplemental material.

Figure 1: The $\pi^{+}\pi^{+}$ and $\pi^{+}\pi^{+}\pi^{+}$ finite-volume
spectra in the center-of-momentum frame for the relevant finite-volume irreps
with various overall momenta, as explained in the text. Points are computed
energy levels on the two volumes with error bars showing statistical
uncertainties. Each rectangular insert shows a vertical zoom of the region
indicated by the small neighboring rectangle. Grey curves are the “non-
interacting” finite-volume energies, i.e. the energies in the absence of any
interactions between pions. Orange curves are predictions from the finite-
volume formalism based only on the two-particle scattering length, given in
Eq. (4) (here with the local three-body interaction set to zero).

_Spectral Determination_ — Figure 1 summarizes the two- and three-pion finite-
volume spectra calculated in this work.222Two-pion energies on the larger
volume have already appeared in Ref. Dudek  _et al._ (2012).

Computations were performed on anisotropic lattices which have a temporal
lattice spacing, $a_{t}$, finer than the spatial lattice spacing, $a_{s}$
($a_{t}=a_{s}/\xi$ with $\xi=3.444(6)$ Dudek  _et al._ (2012)). Two lattice
ensembles were used, differing only in the volume:
$(L/a_{s})^{3}\times(T/a_{t})=20^{3}\times 256$ (with 256 gauge-field
configurations) and $24^{3}\times 128$ (with 512 configurations). We use $2+1$
flavors of dynamical clover fermions, with three-dimensional stout-link
smearing in the fermion action, and a tree-level Symanzik-improved gauge
action. The bare parameters and basic lattice properties are detailed in Refs.
Edwards  _et al._ (2008); Lin  _et al._ (2009). Setting the scale via
${a_{t}^{-1}=m_{\Omega}^{\mathrm{exp}}\,(a_{t}m_{\Omega}^{\mathrm{latt}}})^{-1}$,333where
$a_{t}m_{\Omega}^{\mathrm{latt}}=0.2951(22)$ was measured in Ref. Edwards  _et
al._ (2011) and $m_{\Omega}^{\mathrm{exp}}$ is the experimentally determined
$\Omega$ baryon mass from Ref. Zyla  _et al._ (2020) and combining with
$a_{t}m_{\pi}=0.06906(13)$ Dudek  _et al._ (2012) and $a_{t}m_{K}=0.09698(9)$
Wilson  _et al._ (2015a), yields $m_{\pi}\approx 391\,\text{MeV}$ and
$m_{K}\approx 550\,\text{MeV}$. The values of $a_{t}m_{\pi}$ and $\xi$
translate into spatial extents of $m_{\pi}L=4.76$ and $m_{\pi}L=5.71$ for the
two ensembles.

The spectrum of energies in a finite volume is discrete and each energy level
provides a constraint on the scattering amplitudes at the corresponding
center-of-momentum energy. To obtain more constraints, we compute spectra for
systems with overall zero and non-zero momentum, $\bm{P}$. Momenta are
quantized by the cubic spatial boundary conditions,
$\bm{P}=\tfrac{2\pi}{L}(n_{1},n_{2},n_{3})$, where $\\{n_{i}\\}$ are integers,
and we write this using a shorthand notation as $[n_{1}n_{2}n_{3}]$.

In this work we restrict attention to $S$-wave scattering. The reduced
symmetry of a cubic lattice means that total angular momentum, $J$, is not a
good quantum number and instead channels are labelled by the irreducible
representation (irrep, $\Lambda$) of the octahedral group with parity for
$\bm{P}=\bm{0}$ or the relevant subgroup that leaves $\bm{P}$ invariant for
$\bm{P}\neq\bm{0}$ Johnson (1982); Moore and Fleming (2006). We consider the
relevant irreps which contain $J=0$: $A_{1}^{-}(A_{1}^{+})$ for $\pi\pi\pi$
($\pi\pi$) at rest and $A_{2}(A_{1})$ for $\pi\pi\pi$ ($\pi\pi$) with non-zero
$\bm{P}$. Isospin, $I$, and $G$-parity, $G$, are good quantum numbers in our
lattice formulation; these distinguish the two-pion ($I^{G}=2^{+}$) and three-
pion ($I^{G}=3^{-}$) channels. We neglect higher partial waves here, in
particular the two-particle $D$-wave which mixes with the $S$-wave in the
finite-volume energies. As described in Ref. Dudek  _et al._ (2012), a nonzero
$D$-wave interaction can be extracted, in particular if aided by the
consideration of other, non-trivial finite-volume irreps, but has a small
influence on the two-pion energies considered here.444There is, in principle,
a systematic uncertainty associated with neglecting the $D$-wave contribution.
Given the consistency of our results with Ref. Dudek  _et al._ (2012), this
appears to be below the statistical uncertainty in the present fits. See also
Secs. VIII A and B of that work for more discussion.

To reliably extract the finite-volume energies we have computed two-point
correlation functions using a large basis of appropriate interpolating
operators. From these, the spectra are determined using the variational method
Michael (1985); Lüscher and Wolff (1990); Blossier  _et al._ (2009), with our
implementation described in Refs. Dudek  _et al._ (2008, 2010). This amounts
to calculating a matrix of correlation functions,

$G_{ij}(t)=\langle\mathcal{O}^{\vphantom{\dagger}}_{i}(t)\mathcal{O}^{\dagger}_{j}(0)\rangle\,,$
(1)

and diagonalizing $M(t,t_{0})=G(t_{0})^{-1/2}\cdot G(t)\cdot G(t_{0})^{-1/2}$
for a fixed $t_{0}$. One can show that the corresponding eigenvalues satisfy
$\lambda_{n}(t,t_{0})\to e^{-E_{n}(L)(t-t_{0})}$, where $E_{n}(L)$ is the
$n$’th energy level with overlap to some of the operators in the basis. This
basic methodology has been applied to a wide range of two-hadron scattering
observables for several phenomenologically interesting channels Dudek  _et
al._ (2012, 2013); Briceño  _et al._ (2015, 2016, 2017b); Dudek  _et al._
(2016); Moir  _et al._ (2016); Woss and Thomas (2016); Briceño  _et al._
(2018c); Woss  _et al._ (2018); Wilson  _et al._ (2015b); Cheung  _et al._
(2016); Wilson  _et al._ (2019). See Sec. .1 of the supplemental material for
some example plots of $\lambda_{n}(t,t_{0})$.

In order to robustly interpolate the two- and three-pion energy eigenstates we
use operators with two- and three-meson-like structures in the appropriate
irrep, constructed from products of single-meson-like operators projected to
definite spatial momentum. The latter are built from linear combinations,
chosen to optimize overlap to the single-pion states, of fermion bilinears of
the form, $\bar{\psi}\Gamma{D}\dots D\psi$, where $\psi$ is a quark field and
$D$ is a discretized covariant derivative. Details of these operator
constructions are given in Sec. .5 of the supplemental material with further
details relevant to the three-meson-like operators presented in Ref. Woss  _et
al._ (2019). Using such a wide variety of optimized operators, and especially
multi-hadron operators with momentum-projected single-hadron components,
allows one to minimize excited state contamination and extract the energies
reliably and precisely from small values of $t$. This approach is made
feasible due to the distillation method Peardon  _et al._ (2009) which we
employ to efficiently compute the numerous quark-field Wick contractions that
are required. We use 128 distillation vectors for the $20^{3}$ ensemble and
162 for the $24^{3}$.

Returning to the two- and three-pion spectra summarized in Fig. 1, we observe
a one-to-one correspondence between the computed energy levels and the non-
interacting energies in all panels, with the computed values slightly higher
in energy than the non-interacting levels. This suggests that the system is
weakly interacting and repulsive in both the two- and three-hadron sectors.

Figure 2: Example of data and fits for $\mathcal{K}_{2}$ and
$\mathcal{K}_{3,{\sf iso}}$, as described in the text. The red points are
given by substituting finite-volume energies into $-1/F(E_{2},\bm{P},L)$ and
$-1/F_{3,{\sf iso}}(E_{3},\bm{P},L)$ for the two- and three-particle energies,
respectively, with the volume and $\bm{P}$ indicated in the legend. A symbol
appearing at the very top or bottom represents a case where the central value
falls outside the plotted region. The dark cyan bands represent the fit shown
in Eq. (4) while the lighter bands show the spread covered by the various fits
described in the supplemental material. For the bottom panel we normalize to
$m_{\pi}^{2}\mathcal{K}_{3,{\sf
iso}}^{\text{LO}}=4608\,\pi^{2}(m_{\pi}a_{0})^{2}$, with $m_{\pi}a_{0}$ taken
from Eq. (4). This simple relation between $\mathcal{K}_{3,{\sf iso}}$ and the
two-particle scattering length holds at leading order in chiral perturbation
theory at threshold, as was first derived in Ref. Blanton  _et al._ (2020).
The grey curve gives the full leading-order prediction, which is linear in
$E_{3}^{\star 2}$.

_Analyzing the finite-volume spectra_ — We now describe our method for
determining two- and three-body K matrices from the extracted finite-volume
energies, beginning with an overview of scattering observables:

The two-pion scattering amplitude is defined as the connected part of the
overlap between an incoming $\pi^{+}\pi^{+}$ asymptotic state (with momenta
$\bm{p},-\bm{p}$) to an outgoing $\pi^{+}\pi^{+}$ state (with
$\bm{p}^{\prime},-\bm{p}^{\prime}$). Without loss of generality, here we have
assumed the center-of-momentum frame. We also define
$p=|\bm{p}|=|\bm{p}^{\prime}|$, where we have used that the magnitudes must be
equal to satisfy energy conservation. In addition, $s_{2}\equiv E_{2}^{\star
2}\equiv 4(p^{2}+m_{\pi}^{2})$ defines the squared center-of-momentum frame
energy. The only remaining degree of freedom is the scattering angle between
$\bm{p}$ and $\bm{p}^{\prime}$. In this work we focus on the $S$-wave
scattering amplitude, denoted $\mathcal{M}_{2}$, in which this angle is
integrated to project onto zero-angular-momentum states. Finally we recall the
simple relation between $\mathcal{M}_{2}$ and the K matrix in the elastic
region, $\mathcal{K}_{2}^{-1}=\text{Re}\,\mathcal{M}_{2}^{-1}$.555The
imaginary part of $\mathcal{M}_{2}^{-1}$ is completely fixed by unitarity so
that $\mathcal{K}_{2}$ is the only part free to depend on the dynamics of the
system. We work with the simple phase space factor, proportional to the
momentum magnitude. See, e.g., Ref. Briceño  _et al._ (2018a) for more
details. In contrast to $\mathcal{M}_{2}$, $\mathcal{K}_{2}$ is real for real
$s_{2}$ and is meremorphic in a region of the complex $s_{2}$ plane around
$s_{2}=4m_{\pi}^{2}$. In this work we also consider an analogous, three-body K
matrix, introduced in Ref. Hansen and Sharpe (2014) and denoted by
$\mathcal{K}_{\text{df},3}$.

In the two-pion sector, in the case that the $S$-wave interactions are
dominant, the scalar-irrep finite-volume energies satisfy the quantization
condition Lüscher (1991); Rummukainen and Gottlieb (1995); Kim  _et al._
(2005),

$\displaystyle\mathcal{K}_{2}(E^{\star}_{2})+F^{-1}(E_{2},\bm{P},L)=0\,,$ (2)

where $E_{2}^{\star}\equiv\sqrt{E_{2}^{2}-\bm{P}^{2}}$ is the center-of-
momentum energy and $F(E_{2},\bm{P},L)$ is a known geometric function. For the
three-body sector, we use the isotropic approximation of the general formalism
derived in Ref Hansen and Sharpe (2014), which takes an analogous form, now
for pseudoscalar-irrep energies

$\displaystyle\mathcal{K}_{3,{\sf iso}}(E^{\star}_{3})+F^{-1}_{\rm
3,{\sf{iso}}}[\mathcal{K}_{2}](E_{3},\bm{P},L)=0\,,$ (3)

where the notation is meant to stress that $F_{\rm
3,{\sf{iso}}}[\mathcal{K}_{2}](E_{3},\bm{P},L)$ is a functional of
$\mathcal{K}_{2}(E_{2}^{\star})$. $F_{\rm 3,{\sf{iso}}}$ is defined in Eq.
(39) of Ref. Hansen and Sharpe (2014). Here $\mathcal{K}_{3,{\sf iso}}$ is the
component of $\mathcal{K}_{\mathrm{df},3}$ that only depends on the total
three-hadron energy, i.e. is “isotropic”. Equation (3) holds only when
$\mathcal{K}_{\mathrm{df},3}$ is well approximated to be isotropic and our
fits give evidence that this is a good approximation for this system.

Combining these two conditions with the energies plotted in Fig. 1 allows one
to constrain both the two- and three-hadron K matrices. One strategy is to fit
a parametrization of $\mathcal{K}_{2}$ and use this to determine the energy
dependence of $\mathcal{K}_{3,{\sf iso}}$ as summarized in Fig. 2. An
alternative approach is to parametrize both K matrices and fit these
simultaneously to the entire set of finite-volume energies. A detailed
discussion with a wide range of fits is given in Sec. .2 of the supplemental
material. Both strategies give consistent results and the key message is that
the full data set is well described by a constant $\mathcal{K}_{3,{\sf iso}}$
that is consistent with zero, together with the leading-order effective range
expansion: $\tan\delta(p)=-a_{0}p$ with $\mathcal{K}_{2}(E^{\star}_{2})=-16\pi
E_{2}^{\star}\tan\delta(p)/p$. Here the second equation defines the $S$-wave
scattering phase shift, $\delta(p)$, and the first defines the scattering
length, $a_{0}$. Our best fit, performed simultaneously to all spectra shown
in Fig. 1 but with a cutoff in the center-of-momentum frame energies
included,666This fit is denoted by B2+3 in Sec. .2 of the supplemental
material. As explained there, the fitted data includes all two-pion energies
below $E^{\star}_{2,{\sf cut}}=3.4m_{\pi}$ and all three-pion energies below
$E^{\star}_{3,{\sf cut}}=4.4m_{\pi}$, with both cutoffs applied to energies in
the center-of-momentum frame. yields

$\begin{aligned} m_{\pi}a_{0}&=0.296\pm 0.008\\\\[-3.0pt]
m_{\pi}^{2}\mathcal{K}_{3,{\sf iso}}&=-339\pm 770\end{aligned}\
\bigg{[}\begin{aligned} \\!1.0\ &0.6\\\\[-3.0pt] &1.0\end{aligned}\bigg{]}\,,$
(4)

with a $\chi^{2}$ per degree-of-freedom of $64.5/(37-2)=1.84$. The square-
bracketed matrix gives the correlation between the two fit parameters. This is
consistent with the previous determination of the scattering length at this
pion mass, presented in Ref. Dudek  _et al._ (2012), and is also the value
used to generate the orange curves in Fig. 1 (together with
$\mathcal{K}_{3,{\sf iso}}=0$). In Fig. 2 we illustrate the same fit using the
darker cyan curves. In addition, we include the lighter bands as a systematic
uncertainty, estimated from the spread of various constant and linear fits, as
detailed in Sec. .2 of the supplemental material.

_$3\pi^{+}$ scattering amplitude_ — Following the relativistic integral
equations presented in Ref. Hansen and Sharpe (2015), we can write the $J=0$
and $\mathcal{K}_{3,{\sf iso}}=0$ amplitude as follows:

$\mathcal{M}^{(u,u)}_{3}(p,k)=-\mathcal{M}_{2}(E_{2,p}^{\star})G_{\sf
s}(p,k)\mathcal{M}_{2}(E_{2,k}^{\star})\\\\[3.0pt]
-\mathcal{M}_{2}(E_{2,p}^{\star})\int_{k^{\prime}}\\!G_{\sf
s}(p,k^{\prime})\mathcal{M}^{(u,u)}_{3}(k^{\prime},k)\,,$ (5)

where $\int_{k}\equiv\int dk\,k^{2}/[(2\pi)^{2}\omega_{k}]$ and we have
introduced

$\displaystyle G_{\sf s}(p,k)$
$\displaystyle\equiv-\frac{H(p,k)}{4pk}\log\left[\frac{\alpha(p,k)-2pk+i\epsilon}{\alpha(p,k)+2pk+i\epsilon}\right]\,,$
(6) $\displaystyle\alpha(p,k)$
$\displaystyle\equiv(E_{3}-\omega_{k}-\omega_{p})^{2}-p^{2}-k^{2}-m^{2}\,.$
(7)

$\mathcal{M}_{2}$ is the $S$-wave two-particle scattering amplitude,
introduced above, which depends on the invariant $E_{2,k}^{\star
2}\equiv(E_{3}-\omega_{k})^{2}-k^{2}$, with $\omega_{k}=\sqrt{k^{2}+m^{2}}$.
The $(u,u)$ superscript emphasizes that specific spectator momenta, $k$ and
$p$, are singled out in the initial and final states respectively. The
function $G_{\sf s}$ encodes the spectator exchange, projected to the
$S$-wave. It inherits a scheme dependence through the smooth cutoff function
$H$, defined in Eqs. (28) and (29) of Ref. Hansen and Sharpe (2014). This
scheme dependence is matched by that inside of $\mathcal{K}_{3,{\sf iso}}$
such that the resulting scattering amplitude is universal.

Figure 3: _Top:_ Dalitz-like plot of $m_{\pi}^{4}|\mathcal{M}_{3}|^{2}$ for
$\sqrt{s_{3}}=3.7m$ with final kinematics fixed to $\\{\bm{p}^{\prime
2}_{1},\,\bm{p}^{\prime
2}_{2}\\}=\\{0.01m_{\pi}^{2},\,0.7m_{\pi}^{2}\\}\,\Longrightarrow\,\\{m^{\prime}_{12},\,m^{\prime}_{13}\\}=\\{2.1m_{\pi},\,2.25m_{\pi}\\}$.
_Bottom:_ Same total energy, now with incoming and outgoing kinematics set
equal, as discussed in the text.

To use Eq. (5) in practice, one requires a parameterization for
$\mathcal{M}_{2}$. As described in the previous section, the $\pi^{+}\pi^{+}$
system is well described using the leading order effective range expansion for
$\mathcal{M}_{2}$,

$\displaystyle\mathcal{M}_{2}(E_{2}^{\star})$ $\displaystyle=\frac{16\pi
E_{2}^{\star}}{-1/a_{0}-i\sqrt{E_{2}^{\star 2}/4-m_{\pi}^{2}}}\,.$ (8)

Following the derivation of Ref. Hansen and Sharpe (2015), the final step is
to symmetrize with respect to the spectators, to reach

$\mathcal{M}_{3}(s_{3},m^{\prime 2}_{12},m^{\prime
2}_{13},m^{2}_{12},m^{2}_{13})=\sum_{\bm{p}_{i}\in\mathcal{P}_{p}}\sum_{\bm{k}\in\mathcal{P}_{k}}\mathcal{M}^{(u,u)}_{3}(p,k)\,,$
(9)

where $\mathcal{P}_{p}=\\{\bm{p},\bm{a}^{\prime},-\bm{p}-\bm{a}^{\prime}\\}$
and $\mathcal{P}_{k}=\\{\bm{k},\bm{a},-\bm{k}-\bm{a}\\}$. We have presented
the left-hand side as a function of the five Lorentz invariants that survive
after truncating to $J=0$ in both the two and three particle sector: the
squared three-hadron center-of-momentum frame energy, $s_{3}$, as well as two
pion-pair invariant masses for each of the initial and final states. These are
defined by introducing the notation
$\\{\bm{k},\bm{a},-\bm{k}-\bm{a}\\}=\\{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}\\}$,
then for example

$m_{12}^{2}=(p_{1}+p_{2})^{2}=(E_{3}^{\star}-[{m_{\pi}^{2}+\bm{p}_{3}^{2}}]^{1/2})-\bm{p}_{3}^{2}\,,$
(10)

where the middle expression depends on on-shell four-vectors with
$p_{1}^{2}=m_{\pi}^{2}$.

In the top panel of Fig. 3 we show a Dalitz-like plot of
$|\mathcal{M}_{3}|^{2}$ as a function of $(m_{12},m_{13})$, with all other
kinematics fixed as indicated in the caption. In a usual Dalitz description,
the incoming energy is fixed by the decaying particle so that only the
outgoing kinematics can vary, whereas here we simply fix the other kinematics.
The inputs to this plot are the best-fit scattering length, given in Eq. (4),
together with $\mathcal{K}_{3,{\sf iso}}=0$. The bottom panel of Fig. 3 shows
the same $\sqrt{s_{3}}$ but varies incoming and outgoing kinematics according
to $m_{12}=m_{12}^{\prime}$ and $m_{13}=m_{13}^{\prime}$.

Additional details concerning the $S$-wave integral equations are presented in
Secs. .3 and .4 of the supplemental material, where we also describe the
propagation of the uncertainties of $m_{\pi}a_{0}$ and $\mathcal{K}_{3,{\sf
iso}}$ into the predicted amplitude.777See also Ref. Jackura  _et al._ (2019b)
for more details on expressing the three-particle amplitude via a truncated
partial wave series and Ref. Sadasivan  _et al._ (2020) for a discussion of
integral equations and their solutions in a resonant three-hadron channel.

_Summary_ — In this work we have presented the first lattice QCD determination
of the energy-dependent three-to-three scattering amplitude for three pions
with maximal isospin. The calculation proceeded in three steps: _(i)_
determining finite-volume energies with $\pi^{+}\pi^{+}\pi^{+}$ quantum
numbers, _(ii)_ using the framework of Ref. Hansen and Sharpe (2014) to
extract two- and three-body K matrices from these, and _(iii)_ applying the
results of Ref. Hansen and Sharpe (2015) to convert these to the three-hadron
scattering amplitude, by solving known integral equations. The three steps are
summarized, respectively, by Figs. 1, 2 and 3 of the text.

Having established this general workflow, it is now well within reach to
rigorously extract three-hadron resonance properties from lattice QCD
calculations. In particular the formalism has recently been extended to three-
pion states with any value of isospin in Ref. Hansen  _et al._ (2020). This
should enable studies, for example, of the $\omega$, $h_{1}$ and $a_{1}$
resonances. The main outstanding challenges here include rigorous resonant
parametrizations of the intermediate three-body K matrix, as well as a better
understanding of the analytic continuation required to identify the resonance
pole position.

###### Acknowledgements.

The authors would like to thank M. Bruno, J. Dudek, A. Jackura, L. Leskovec
and A. Rodas for useful conversations, as well as our other colleagues within
the Hadron Spectrum Collaboration. The authors would also like to thank W.
Detmold, F. Romero-López and S. R. Sharpe for useful discussions and for
valuable feedback on a previous version of this manuscript. RAB, RGE and CET
acknowledge support from the U.S. Department of Energy contract DE-
AC05-06OR23177, under which Jefferson Science Associates, LLC, manages and
operates Jefferson Lab. RAB also acknowledges support from the USDOE Early
Career award, contract de-sc0019229. CET and DJW acknowledge support from the
U.K. Science and Technology Facilities Council (STFC) [grant number
ST/P000681/1]. DJW acknowledges support from a Royal Society University
Research Fellowship. MTH, CET and DJW acknowledge the MITP topical workshop
“Scattering Amplitudes and Resonance Properties for Lattice QCD” for
stimulating this project and for hospitality during the initial discussions.
MTH and CET also acknowledge the CERN-TH Institute “Advances in Lattice Gauge
Theory”, which provided the opportunity to make significant progress on this
work. CET further acknowledges CERN TH for hospitality and support during a
visit in January and February of this year. The software codes Chroma Edwards
and Joo (2005) and QUDA Clark  _et al._ (2010); Babich  _et al._ (2010); Clark
_et al._ (2016) were used. The authors acknowledge support from the U.S.
Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research and Office of Nuclear Physics, Scientific Discovery through
Advanced Computing (SciDAC) program. Also acknowledged is support from the
U.S. Department of Energy Exascale Computing Project. This work was also
performed on clusters at Jefferson Lab under the USQCD Collaboration and the
LQCD ARRA Project. This research was supported in part under an ALCC award,
and used resources of the Oak Ridge Leadership Computing Facility at the Oak
Ridge National Laboratory, which is supported by the Office of Science of the
U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research
used resources of the National Energy Research Scientific Computing Center
(NERSC), a DOE Office of Science User Facility supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
The authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin for providing HPC resources. Gauge
configurations were generated using resources awarded from the U.S. Department
of Energy INCITE program at the Oak Ridge Leadership Computing Facility, the
NERSC, the NSF Teragrid at the TACC and the Pittsburgh Supercomputer Center,
as well as at Jefferson Lab.

## References

  * Piarulli  _et al._ (2018) M. Piarulli _et al._ , Phys. Rev. Lett. 120, 052503 (2018), arXiv:1707.02883 [nucl-th] .
  * Isgur and Karl (1977) N. Isgur and G. Karl, Phys. Lett. B 72, 109 (1977).
  * Isgur and Karl (1979) N. Isgur and G. Karl, Phys. Rev. D 19, 2653 (1979), [Erratum: Phys.Rev.D 23, 817 (1981)].
  * Hansen and Sharpe (2014) M. T. Hansen and S. R. Sharpe, Phys. Rev. D90, 116003 (2014), arXiv:1408.5933 [hep-lat] .
  * Hansen and Sharpe (2015) M. T. Hansen and S. R. Sharpe, Phys. Rev. D92, 114509 (2015), arXiv:1504.04248 [hep-lat] .
  * Lüscher (1991) M. Lüscher, Nucl. Phys. B354, 531 (1991).
  * Rummukainen and Gottlieb (1995) K. Rummukainen and S. A. Gottlieb, Nucl. Phys. B450, 397 (1995), arXiv:hep-lat/9503028 [hep-lat] .
  * Bedaque (2004) P. F. Bedaque, Phys. Lett. B593, 82 (2004), arXiv:nucl-th/0402051 [nucl-th] .
  * Kim  _et al._ (2005) C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl. Phys. B727, 218 (2005), arXiv:hep-lat/0507006 [hep-lat] .
  * Fu (2012) Z. Fu, Phys. Rev. D85, 014506 (2012), arXiv:1110.0319 [hep-lat] .
  * Leskovec and Prelovsek (2012) L. Leskovec and S. Prelovsek, Phys. Rev. D85, 114507 (2012), arXiv:1202.2145 [hep-lat] .
  * Göckeler  _et al._ (2012) M. Göckeler, R. Horsley, M. Lage, U. G. Meißner, P. E. L. Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti, Phys. Rev. D86, 094513 (2012), arXiv:1206.4141 [hep-lat] .
  * He  _et al._ (2005) S. He, X. Feng, and C. Liu, JHEP 07, 011 (2005), arXiv:hep-lat/0504019 [hep-lat] .
  * Lage  _et al._ (2009) M. Lage, U.-G. Mei$\ss$ner, and A. Rusetsky, Phys. Lett. B681, 439 (2009), arXiv:0905.0069 [hep-lat] .
  * Bernard  _et al._ (2011) V. Bernard, M. Lage, U. G. Mei$\ss$ner, and A. Rusetsky, JHEP 01, 019 (2011), arXiv:1010.6018 [hep-lat] .
  * Döring  _et al._ (2011) M. Döring, U.-G. Mei$\ss$ner, E. Oset, and A. Rusetsky, Eur. Phys. J. A47, 139 (2011), arXiv:1107.3988 [hep-lat] .
  * Döring and Mei$\ss$ner (2012) M. Döring and U. G. Mei$\ss$ner, JHEP 01, 009 (2012), arXiv:1111.0616 [hep-lat] .
  * Agadjanov  _et al._ (2014) D. Agadjanov, U.-G. Mei$\ss$ner, and A. Rusetsky, JHEP 01, 103 (2014), arXiv:1310.7183 [hep-lat] .
  * Döring  _et al._ (2012) M. Döring, U. G. Mei$\ss$ner, E. Oset, and A. Rusetsky, Eur. Phys. J. A48, 114 (2012), arXiv:1205.4838 [hep-lat] .
  * Hansen and Sharpe (2012) M. T. Hansen and S. R. Sharpe, Phys. Rev. D86, 016007 (2012), arXiv:1204.0826 [hep-lat] .
  * Briceño and Davoudi (2013a) R. A. Briceño and Z. Davoudi, Phys. Rev. D88, 094507 (2013a), arXiv:1204.1110 [hep-lat] .
  * Guo  _et al._ (2013) P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak, Phys. Rev. D88, 014501 (2013), arXiv:1211.0929 [hep-lat] .
  * Briceño  _et al._ (2018a) R. A. Briceño, J. J. Dudek, and R. D. Young, Rev. Mod. Phys. 90, 025001 (2018a), arXiv:1706.06223 [hep-lat] .
  * Briceño (2014) R. A. Briceño, Phys. Rev. D89, 074507 (2014), arXiv:1401.3312 [hep-lat] .
  * Romero-López  _et al._ (2019) F. Romero-López, S. R. Sharpe, T. D. Blanton, R. A. Briceño, and M. T. Hansen, JHEP 10, 007 (2019), arXiv:1908.02411 [hep-lat] .
  * Briceño  _et al._ (2019a) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D99, 014516 (2019a), arXiv:1810.01429 [hep-lat] .
  * Briceño  _et al._ (2017a) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D95, 074510 (2017a), arXiv:1701.07465 [hep-lat] .
  * Hansen and Sharpe (2019) M. T. Hansen and S. R. Sharpe, Ann. Rev. Nucl. Part. Sci. 69, 65 (2019), arXiv:1901.00483 [hep-lat] .
  * Briceño  _et al._ (2018b) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D98, 014506 (2018b), arXiv:1803.04169 [hep-lat] .
  * Blanton  _et al._ (2019) T. D. Blanton, F. Romero-López, and S. R. Sharpe, JHEP 03, 106 (2019), arXiv:1901.07095 [hep-lat] .
  * Hammer  _et al._ (2017a) H.-W. Hammer, J.-Y. Pang, and A. Rusetsky, JHEP 09, 109 (2017a), arXiv:1706.07700 [hep-lat] .
  * Hammer  _et al._ (2017b) H. W. Hammer, J. Y. Pang, and A. Rusetsky, JHEP 10, 115 (2017b), arXiv:1707.02176 [hep-lat] .
  * Briceño and Davoudi (2013b) R. A. Briceño and Z. Davoudi, Phys. Rev. D87, 094507 (2013b), arXiv:1212.3398 [hep-lat] .
  * Polejaeva and Rusetsky (2012) K. Polejaeva and A. Rusetsky, Eur. Phys. J. A48, 67 (2012), arXiv:1203.1241 [hep-lat] .
  * Mai and Döring (2017) M. Mai and M. Döring, Eur. Phys. J. A53, 240 (2017), arXiv:1709.08222 [hep-lat] .
  * Guo and Gasparian (2017) P. Guo and V. Gasparian, Phys. Lett. B774, 441 (2017), arXiv:1701.00438 [hep-lat] .
  * Guo and Gasparian (2018) P. Guo and V. Gasparian, Phys. Rev. D97, 014504 (2018), arXiv:1709.08255 [hep-lat] .
  * Guo and Morris (2019) P. Guo and T. Morris, Phys. Rev. D99, 014501 (2019), arXiv:1808.07397 [hep-lat] .
  * Romero-López  _et al._ (2018) F. Romero-López, A. Rusetsky, and C. Urbach, Eur. Phys. J. C78, 846 (2018), arXiv:1806.02367 [hep-lat] .
  * Jackura  _et al._ (2019a) A. W. Jackura, S. M. Dawid, C. Fernández-Ramírez, V. Mathieu, M. Mikhasenko, A. Pilloni, S. R. Sharpe, and A. P. Szczepaniak, Phys. Rev. D100, 034508 (2019a), arXiv:1905.12007 [hep-ph] .
  * Briceño  _et al._ (2019b) R. A. Briceño, M. T. Hansen, S. R. Sharpe, and A. P. Szczepaniak, Phys. Rev. D100, 054508 (2019b), arXiv:1905.11188 [hep-lat] .
  * Hansen  _et al._ (2020) M. T. Hansen, F. Romero-López, and S. R. Sharpe, JHEP 07, 047 (2020), arXiv:2003.10974 [hep-lat] .
  * Blanton and Sharpe (2020a) T. D. Blanton and S. R. Sharpe, (2020a), arXiv:2007.16190 [hep-lat] .
  * Blanton and Sharpe (2020b) T. D. Blanton and S. R. Sharpe, (2020b), arXiv:2007.16188 [hep-lat] .
  * Beane  _et al._ (2020) S. Beane _et al._ , (2020), arXiv:2003.12130 [hep-lat] .
  * Beane  _et al._ (2008) S. R. Beane, W. Detmold, T. C. Luu, K. Orginos, M. J. Savage, and A. Torok, Phys. Rev. Lett. 100, 082004 (2008), arXiv:0710.1827 [hep-lat] .
  * Detmold  _et al._ (2008) W. Detmold, M. J. Savage, A. Torok, S. R. Beane, T. C. Luu, K. Orginos, and A. Parreno, Phys. Rev. D78, 014507 (2008), arXiv:0803.2728 [hep-lat] .
  * Detmold  _et al._ (2012) W. Detmold, K. Orginos, and Z. Shi, Phys. Rev. D 86, 054507 (2012), arXiv:1205.4224 [hep-lat] .
  * Mai and Döring (2019) M. Mai and M. Döring, Phys. Rev. Lett. 122, 062503 (2019), arXiv:1807.04746 [hep-lat] .
  * Blanton  _et al._ (2020) T. D. Blanton, F. Romero-López, and S. R. Sharpe, Phys. Rev. Lett. 124, 032001 (2020), arXiv:1909.02973 [hep-lat] .
  * Mai  _et al._ (2020) M. Mai, M. Döring, C. Culver, and A. Alexandru, Phys. Rev. D 101, 054510 (2020), arXiv:1909.05749 [hep-lat] .
  * Guo and Long (2020) P. Guo and B. Long, Phys. Rev. D 101, 094510 (2020), arXiv:2002.09266 [hep-lat] .
  * Hörz and Hanlon (2019) B. Hörz and A. Hanlon, Phys. Rev. Lett. 123, 142002 (2019), arXiv:1905.04277 [hep-lat] .
  * Culver  _et al._ (2020) C. Culver, M. Mai, R. Brett, A. Alexandru, and M. Döring, Phys. Rev. D 101, 114507 (2020), arXiv:1911.09047 [hep-lat] .
  * Fischer  _et al._ (2020) M. Fischer, B. Kostrzewa, L. Liu, F. Romero-López, M. Ueding, and C. Urbach, (2020), arXiv:2008.03035 [hep-lat] .
  * Dudek  _et al._ (2012) J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D86, 034031 (2012), arXiv:1203.6041 [hep-ph] .
  * Edwards  _et al._ (2008) R. G. Edwards, B. Joo, and H.-W. Lin, Phys. Rev. D78, 054501 (2008), arXiv:0803.3960 [hep-lat] .
  * Lin  _et al._ (2009) H.-W. Lin _et al._ (Hadron Spectrum), Phys. Rev. D79, 034502 (2009), arXiv:0810.3588 [hep-lat] .
  * Edwards  _et al._ (2011) R. G. Edwards, J. J. Dudek, D. G. Richards, and S. J. Wallace, Phys. Rev. D84, 074508 (2011), arXiv:1104.5152 [hep-ph] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Wilson  _et al._ (2015a) D. J. Wilson, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D91, 054008 (2015a), arXiv:1411.2004 [hep-ph] .
  * Johnson (1982) R. C. Johnson, Phys. Lett. 114B, 147 (1982).
  * Moore and Fleming (2006) D. C. Moore and G. T. Fleming, Phys. Rev. D73, 014504 (2006), [Erratum: Phys. Rev.D74,079905(2006)], arXiv:hep-lat/0507018 [hep-lat] .
  * Michael (1985) C. Michael, Nucl. Phys. B259, 58 (1985).
  * Lüscher and Wolff (1990) M. Lüscher and U. Wolff, Nucl. Phys. B 339, 222 (1990).
  * Blossier  _et al._ (2009) B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer, JHEP 04, 094 (2009), arXiv:0902.1265 [hep-lat] .
  * Dudek  _et al._ (2008) J. J. Dudek, R. G. Edwards, N. Mathur, and D. G. Richards, Phys. Rev. D77, 034501 (2008), arXiv:0707.4162 [hep-lat] .
  * Dudek  _et al._ (2010) J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards, and C. E. Thomas, Phys. Rev. D82, 034508 (2010), arXiv:1004.4930 [hep-ph] .
  * Dudek  _et al._ (2013) J. J. Dudek, R. G. Edwards, and C. E. Thomas (Hadron Spectrum), Phys. Rev. D87, 034505 (2013), [Erratum: Phys. Rev.D90,no.9,099902(2014)], arXiv:1212.0830 [hep-ph] .
  * Briceño  _et al._ (2015) R. A. Briceño, J. J. Dudek, R. G. Edwards, C. J. Shultz, C. E. Thomas, and D. J. Wilson, Phys. Rev. Lett. 115, 242001 (2015), arXiv:1507.06622 [hep-ph] .
  * Briceño  _et al._ (2016) R. A. Briceño, J. J. Dudek, R. G. Edwards, C. J. Shultz, C. E. Thomas, and D. J. Wilson, Phys. Rev. D93, 114508 (2016), arXiv:1604.03530 [hep-ph] .
  * Briceño  _et al._ (2017b) R. A. Briceño, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. Lett. 118, 022002 (2017b), arXiv:1607.05900 [hep-ph] .
  * Dudek  _et al._ (2016) J. J. Dudek, R. G. Edwards, and D. J. Wilson (Hadron Spectrum), Phys. Rev. D93, 094506 (2016), arXiv:1602.05122 [hep-ph] .
  * Moir  _et al._ (2016) G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, and D. J. Wilson, JHEP 10, 011 (2016), arXiv:1607.07093 [hep-lat] .
  * Woss and Thomas (2016) A. J. Woss and C. E. Thomas,  _Proceedings, 34th International Symposium on Lattice Field Theory (Lattice 2016): Southampton, UK, July 24-30, 2016_ , PoS LATTICE2016, 134 (2016), arXiv:1612.05437 [hep-lat] .
  * Briceño  _et al._ (2018c) R. A. Briceño, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D97, 054513 (2018c), arXiv:1708.06667 [hep-lat] .
  * Woss  _et al._ (2018) A. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards, and D. J. Wilson, JHEP 07, 043 (2018), arXiv:1802.05580 [hep-lat] .
  * Wilson  _et al._ (2015b) D. J. Wilson, R. A. Briceño, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D92, 094502 (2015b), arXiv:1507.02599 [hep-ph] .
  * Cheung  _et al._ (2016) G. K. C. Cheung, C. O’Hara, G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, and D. Tims (Hadron Spectrum), JHEP 12, 089 (2016), arXiv:1610.01073 [hep-lat] .
  * Wilson  _et al._ (2019) D. J. Wilson, R. A. Briceño, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. Lett. 123, 042002 (2019), arXiv:1904.03188 [hep-lat] .
  * Woss  _et al._ (2019) A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D100, 054506 (2019), arXiv:1904.04136 [hep-lat] .
  * Peardon  _et al._ (2009) M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), Phys. Rev. D80, 054506 (2009), arXiv:0905.2160 [hep-lat] .
  * Jackura  _et al._ (2019b) A. Jackura, C. Fernández-Ramírez, V. Mathieu, M. Mikhasenko, J. Nys, A. Pilloni, K. Saldaña, N. Sherrill, and A. P. Szczepaniak (JPAC), Eur. Phys. J. C79, 56 (2019b), arXiv:1809.10523 [hep-ph] .
  * Sadasivan  _et al._ (2020) D. Sadasivan, M. Mai, H. Akdag, and M. Döring, Phys. Rev. D 101, 094018 (2020), arXiv:2002.12431 [nucl-th] .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),  _Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004, Batavia, USA, June 21-26, 2004_ , Nucl. Phys. Proc. Suppl. 140, 832 (2005), arXiv:hep-lat/0409003 [hep-lat] .
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2010) R. Babich, M. A. Clark, and B. Joo, in  _SC 10 (Supercomputing 2010) New Orleans, Louisiana, November 13-19, 2010_ (2010) arXiv:1011.0024 [hep-lat] .
  * Clark  _et al._ (2016) M. A. Clark, B. Joó, A. Strelchenko, M. Cheng, A. Gambhir, and R. C. Brower, in  _SC ’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ (2016) pp. 795–806.

## Supplemental Material

###  .1 Spectra

In this section we give details concerning the finite-volume spectra described
in the main text. We focus here on two representative examples for the
correlators used to extract the three-pion energies. The quality of two-pion
correlators can be inferred from the earlier work presented in Ref. Dudek  _et
al._ (2012), which includes a partially overlapping data set.

Figure 4: GEVP eigenvalues used to determine the $\bm{P}=[000]$,
$\Lambda=A_{1}^{-}$, three-pion spectrum on the $20^{3}\times 256$ ensemble.
As explained in the text, the diagonalization is performed with
$t_{0}/a_{t}=10$ and the eigenvalues are rescaled by their expected large-$t$
fall-off.

As a first example, consider the three-pion spectrum for the $\bm{P}=[000]$,
$\Lambda=A_{1}^{-}$ irrep on the $20^{3}\times 256$ ensemble. In this case
$G_{ij}(t)=\langle\mathcal{O}_{i}(t)\mathcal{O}^{\dagger}_{j}(0)\rangle$ is an
$8\times 8$ matrix of correlators, built from the first 8 operators listed in
Table 6 of Sec. .5. In Fig. 4 we plot the corresponding eigenvalues,
$\lambda_{n}(t,t_{0})$, of the matrix

$M(t,t_{0})\equiv G(t_{0})^{-1/2}\cdot G(t)\cdot G(t_{0})^{-1/2}\,,$ (11)

entering the generalized eigenvalue problem (GEVP). These are determined for
$t_{0}/a_{t}=10$ and are plotted vs. $t/a_{t}$ for a range of values both
before and after the reference time. To display the eigenvalues in a useful
manner, we plot the combination $e^{E_{n}(t-t_{0})}\lambda_{n}(t,t_{0})$,
where $E_{n}$ has been determined from a two-state fit to
$\lambda_{n}(t,t_{0})$

$\lambda_{n}(t,t_{0})=(1-A_{n})e^{-E_{n}(t-t_{0})}+A_{n}e^{-E^{\prime}_{n}(t-t_{0})}\,.$
(12)

The quality of the fit is indicated by the $\chi^{2}/\text{dof}$ in each
panel. The plotted combination behaves as expected for a successful GEVP,
showing a reasonable plateaux over a range of $t/a_{t}$.

This result also exhibits no evidence for thermal states on this lattice, as
expected given the length of the temporal extent, $m_{\pi}T\approx 17.7$. We
detour slightly, to explain this point in more detail:

In general, for multi-pion systems, the leading finite-$T$ effects are given
by

$G_{ij}(t)=\langle
0|\mathcal{O}_{i}(0)\,e^{-\hat{H}t}\,\mathcal{O}^{\dagger}_{j}(0)|0\rangle+e^{-m_{\pi}(T-t)}\langle\pi^{-}|\mathcal{O}_{i}(0)\,e^{-\hat{H}t}\,\mathcal{O}^{\dagger}_{j}(0)|\pi^{-}\rangle+\cdots\,,$
(13)

where $\hat{H}$ is the Hamiltonian and the ellipsis represents thermal effects
falling faster than $e^{-m_{\pi}T}$. For concreteness, we have assumed that
$O^{\dagger}(0)$ creates three-$\pi^{+}$ quantum numbers, so that
$\mathcal{O}^{\dagger}(0)|\pi^{-}\rangle$ has the quantum numbers of two pions
with isospin two. A spectral decomposition of Eq. (13) then yields

$G_{ij}(t)=\sum_{n}c^{(n)}_{i}c^{(n)*}_{j}e^{-E_{n}^{\pi\pi\pi}t}+e^{-m_{\pi}(T-t)}\sum_{n}b^{(n)}_{i}b^{(n)*}_{j}e^{-E_{n}^{\pi\pi}t}+\cdots\,,$
(14)

where the sum in the first term (second term) runs over all maximal-isospin
three-pion (two-pion) finite-volume states with specified $\bm{P}$. In the
case of $\bm{P}=[000]$, the two- and three-pion ground states take the form
$Nm_{\pi}+\Delta E_{N}$ where $N=2,3$ and $\Delta E_{N}\sim 1/L^{3}$ for
weakly-interacting systems. Taking the leading ($n=0$) terms of Eq. (14) and
substituting this scaling then yields

$G_{ij}(t)=c^{(0)}_{i}c^{(0)*}_{j}e^{-(3m_{\pi}+\Delta
E_{3})t}+b^{(0)}_{i}b^{(0)*}_{j}e^{-m_{\pi}(T+t)}e^{-\Delta E_{2}t}+\cdots\,.$
(15)

When this same exercise is performed for a _two-pion_ correlator, again
focusing on the case of $\bm{P}=[000]$, one finds that the leading thermal
contamination is a constant in $t$ in the non-interacting limit. As discussed
in Ref. Dudek  _et al._ (2012) this can thus be removed by applying a shift to
the correlator $G_{ij}(t)\to G_{ij}(t)-G_{ij}(t+\delta t)$. In the present
case, however, the leading contaminations are $t$-dependent. One option is to
reweight and shift, i.e.

$G_{ij}(t)\ \to\ e^{-m_{\pi}t}\big{(}G_{ij}(t)e^{m_{\pi}t}-G_{ij}(t+\delta
t)e^{m_{\pi}(t+\delta t)}\big{)}\,.$ (16)

This approach, already used in Ref. Dudek  _et al._ (2012) for $\pi\pi$
systems with non-zero total momentum, reduces thermal effects at the cost of
generally degrading the signal quality. Fortunately, for the $20^{3}\times
256$ lattice, this is not required. Comparing the leading and subleading terms
of Eq. (15), and neglecting the interactions, one finds that the relative size
of the three-pion thermal contamination is $e^{-m_{\pi}(T-2t)}$. Thus,
assuming that the relevant matrix elements have the same order of magnitude,
for the range of $t$ considered these effects are $\sim 10^{-7}$ and are
safely negligible, despite the high statistical precision of the extracted
energies. This concludes our general comments on thermal effects.

Figure 5: GEVP eigenvalues used to determine the $\bm{P}=[000]$,
$\Lambda=A_{1}^{-}$, three-pion spectrum on the $24^{3}\times 128$ ensemble.
The diagonalization is performed with $t_{0}/a_{t}=9$ and the eigenvalues are
rescaled as in Fig. 4.

As a second example, in Fig. 5 we consider the same three-pion quantum numbers
($\bm{P}=[000]$, $\Lambda=A_{1}^{-}$) on the $24^{3}\times 128$ ensemble.
Because the larger spatial volume lowers the value of the $n$th level, here we
include 2 additional operators to better absorb the excited states. As with
the previous example, the $\chi^{2}/\text{dof}$ and the plotted curves provide
strong evidence of a successful GEVP extraction. For this case,
$e^{-m_{\pi}(T-2t)}\sim 10^{-3}$ so that finite-$T$ effects potentially
present a more significant issue. As a result we have also considered shifting
and reweighting, as summarized by Eq. (16), in our various fits. However,
across all values of $\bm{P}$, we find that more stable fits are achieved via
the unmodified correlators, relying on the basis of operators to push the
extraction to earlier times and examining the resulting GEVP eigenvalues. This
is in contrast to the $24^{3}\times 128$ two-pion fits, where the extractions
are improved by shifting in certain cases, as described in Ref. Dudek  _et
al._ (2012).

###  .2 K-matrix fits

In this subsection we give additional details concerning the K-matrix fits,
summarized in the main text. We present four basic types of fits:

  1. 1.

Fitting only the $\pi^{+}\pi^{+}$ spectra to various choices of
$p\cot\delta(p)$ [Table 1].

  2. 2.

Fitting only the $\pi^{+}\pi^{+}\pi^{+}$ spectra to various choices of
$p\cot\delta(p)$, with $\mathcal{K}_{3,{\sf iso}}=0$ [Table 2].

  3. 3.

Fitting only the $\pi^{+}\pi^{+}\pi^{+}$ spectra to various choices of
$\mathcal{K}_{3,{\sf iso}}$, with $p\cot\delta(p)$ fixed by independent
$\pi^{+}\pi^{+}$ fits [Table 3].

  4. 4.

Fitting all spectra simultaneously to various choices of $p\cot\delta(p)$ and
$\mathcal{K}_{3,{\sf iso}}$ [Table 4].

Here $\delta(p)$ is the $S$-wave, $\pi^{+}\pi^{+}$ scattering phase shift,
related to the scattering amplitude via

$\mathcal{M}_{2}(E_{2}^{\star})=\frac{16\pi
E_{2}^{\star}}{p\cot\delta(p)-ip}\,,$ (17)

where $p^{2}=E_{2}^{\star 2}/4-m_{\pi}^{2}$. One standard parametrization of
the scattering amplitude follows from the effective range expansion

$p\cot\delta(p)=-\frac{1}{a_{0}}+\frac{1}{2}r_{0}p^{2}+\mathcal{O}(p^{4})\,,$
(18)

and below we present fits to the leading term as well as to the leading two
terms.

Fit |  $E^{\star}_{2,\text{cut}}$ |  $p\,\text{cot}\delta(p)$ |  fit result |  $\chi^{2}/\text{dof}$
---|---|---|---|---
A2 | $4.0m_{\pi}$ |  $-1/a_{0}$ |  $m_{\pi}a_{0}=0.278\pm 0.007$ |  $89.8/(32-1)=2.90$
B2 | $3.4m_{\pi}$ |  $-1/a_{0}$ |  $m_{\pi}a_{0}=0.292\pm 0.010$ |  $26.9/(21-1)=1.35$
C2 | ${4.0m_{\pi}}$ |  $-1/a_{0}+r_{0}p^{2}/2$ |  $\begin{aligned} m_{\pi}a_{0}&={0.317\pm 0.015}\\\\[-3.0pt] m_{\pi}r_{0}&={-0.39\pm 0.12}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{-0.9}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${79.6/(32-2)=2.65}$
D2 | ${3.4m_{\pi}}$ |  $-1/a_{0}+r_{0}p^{2}/2$ |  $\begin{aligned} m_{\pi}a_{0}&={0.258\pm 0.018}\\\\[-3.0pt] m_{\pi}r_{0}&={0.68\pm 0.33}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{-0.9}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${22.7/(21-2)=1.20}$
E2 | ${4.0m_{\pi}}$ |  $-A(c_{0},p)/a_{0}$ |  $\begin{aligned} m_{\pi}a_{0}&={0.355\pm 0.021}\\\\[-3.0pt] m_{\pi}c_{0}&={11.2\pm 2.1}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{-0.8}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${96.7/(32-2)=3.22}$
F2 | ${3.4m_{\pi}}$ |  $-A(c_{0},p)/a_{0}$ |  $\begin{aligned} m_{\pi}a_{0}&={0.260\pm 0.035}\\\\[-3.0pt] m_{\pi}c_{0}&={3.7\pm 1.1}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{-0.9}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${22.7/(21-2)=1.20}$
G2 | ${4.0m_{\pi}}$ |  $A(1,p)(-1/a_{0}+c_{0}p^{2})$ |  $\begin{aligned} m_{\pi}a_{0}&={0.223\pm 0.019}\\\\[-3.0pt] m_{\pi}c_{0}&={-2.88\pm 0.19}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{\phantom{+}0.9}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${70.1/(32-2)=2.34}$
H2 | ${3.4m_{\pi}}$ |  $A(1,p)(-1/a_{0}+c_{0}p^{2})$ |  $\begin{aligned} m_{\pi}a_{0}&={0.184\pm 0.022}\\\\[-3.0pt] m_{\pi}c_{0}&={-2.2\pm 0.4}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{\phantom{+}0.96}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${26.5/(21-2)=1.40}$
I2 $[r=0.01\ (9/32)]$ | $4.0m_{\pi}$ |  $-1/a_{0}$ |  $m_{\pi}a_{0}=0.292\pm 0.008$ |  $40.5/(32-1)=1.31$
J2 $[r=0.01\ (9/32)]$ | ${4.0m_{\pi}}$ |  $-1/a_{0}+r_{0}p^{2}/2$ |  $\begin{aligned} m_{\pi}a_{0}&={0.300\pm 0.016}\\\\[-3.0pt] m_{\pi}r_{0}&={-0.08\pm 0.14}\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &{-0.8}\\\\[-3.0pt] &\phantom{+}1.0\end{aligned}\bigg{]}$ | ${40.3/(32-2)=1.34}$
Table 1: Summary of fits to $\pi^{+}\pi^{+}$ finite-volume energies, for
various choices of $p\cot\delta(p)$. All values of total momentum $\bm{P}$
(from $[000]$ to $[002]$) and both volumes ($20^{3}$ and $24^{3}$) are used in
each fit. The entries below the lower double horizontal line are determined
using a regularized covariance matrix as explained in the text. The function
$A(c_{0},p)=c_{0}m_{\pi}\sqrt{p^{2}+m_{\pi}^{2}}/(2p^{2}+c_{0}m_{\pi}^{2})$
encodes the effect of the Adler zero, with $c_{0}=1$ corresponding to the pole
position from leading-order chiral perturbation theory. The columns are
understood as follows: “Fit” gives a label to the fit (and defines the
regularized covariance matrix for the final two fits);
“$E^{\star}_{2,\text{cut}}$” gives the two-particle center-of-momentum frame
energy cutoff (i.e. only points with central values below this threshold enter
the fit); “$p\,\text{cot}\delta(p)$” indicates the fit function; “fit result”
displays the extracted parameters and their correlation;
“$\chi^{2}/\text{dof}$” gives the value of $\chi^{2}(\\{\eta_{i}\\})$
(evaluated at the best fit parameters) divided by the number of degrees of
freedom.

In the case of $\pi^{+}\pi^{+}$ scattering, chiral perturbation theory
predicts the Adler zero, which leads to a pole in $p\cot\delta(p)$, limiting
the range of convergence for the effective range expansion. This motivates the
alternative form

$p\cot\delta(p)=A(c,p)\Big{[}\\!-\frac{1}{a_{0}}+c^{\prime}p^{2}+\mathcal{O}(p^{4})\Big{]}\,,$
(19)

where we have introduced

$A(c,p)\equiv\frac{c\,m_{\pi}\sqrt{p^{2}+m_{\pi}^{2}}}{2p^{2}+c\,m_{\pi}^{2}}\,.$
(20)

The function $A(c,p)$ is chosen to match the analytic structure predicted by
leading-order chiral perturbation theory (both the energy numerator and the
pole in the denominator) and is normalized so that $A(c,0)=1$. The leading-
order prediction for the pole position corresponds to $c=1$ and in the
following we present fits both with $c$ fixed and allowed to vary.

As we explain in detail in Sec. .3 below, the three-particle scattering
amplitude, $\mathcal{M}_{3}$, is determined from the two-particle scattering
amplitude together with a local three-particle K matrix, first introduced in
Ref. Hansen and Sharpe (2014) and denoted by $\mathcal{K}_{\text{df},3}$. As
already described in the main text, we work here in the isotropic
approximation, for which this three-particle K matrix reduces to a simple
function of the total center-of-momentum frame energy, denoted
$\mathcal{K}_{3,{\sf iso}}(E_{3}^{\star})$. This quantity admits an expansion
similar to the effective range expansion

$\mathcal{K}_{3,{\sf
iso}}(E_{3}^{\star})=c_{1}/m_{\pi}^{2}+c_{2}\Delta/m_{\pi}^{4}+\mathcal{O}(\Delta^{2})\,,$
(21)

where $\Delta\equiv E_{3}^{\star 2}-9m_{\pi}^{2}$. The fits presented below
take either the first or else the first two terms in this expansion.

In each case, the fit is performed by minimizing the
$\chi^{2}(\\{\eta_{i}\\})$, where

$\chi^{2}(\\{\eta_{i}\\})\equiv[\bm{\mathcal{E}}_{\sf
d}-\bm{\mathcal{E}}(\\{\eta_{i}\\})]\cdot C^{-1}\cdot[\bm{\mathcal{E}}_{\sf
d}-\bm{\mathcal{E}}(\\{\eta_{i}\\})]^{\text{T}}\,.$ (22)

Here $\bm{\mathcal{E}}_{\sf d}$ is a vector built from two and three-particle
energies extracted from the lattice calculation, $C$ is the covariance matrix,
and $\bm{\mathcal{E}}(\\{\eta_{i}\\})$ is a vector of solutions to the two-
and three-particle quantization conditions. The solved energies depend on
$\\{\eta_{i}\\}$, which stands for all two- and three-particle K-matrix
parameters, over which the minimization is performed.

In certain cases the low-lying eigenvalues of $C$ cannot be estimated reliably
and, if underestimated, can lead to artificially enhanced eigenvalues in
$C^{-1}$, and therefore inflated values for $\chi^{2}(\\{\eta_{i}\\})$. To
study this problem we have also considered an alternative method in which the
low lying eigenvalues of $C$ are adjusted. To do so, one first diagonalizes
$C$

$C=R^{\text{T}}\cdot D\cdot R\,,$ (23)

where $R$ is an orthogonal matrix of eigenvectors and $D$ a diagonal matrix of
eigenvalues. We label the eigenvalues (ordered from smallest to largest) by
$\lambda_{1},\cdots,\lambda_{N}$ and note that these are positive. We assume
also that the rows and columns of $R$ are organized such that

$D=\text{diag}[\lambda_{1},\lambda_{2},\cdots,\lambda_{N}]\,.$ (24)

$\chi^{2}(\\{\eta_{i}\\})$ may be poorly estimated if there is a large
hierarchy between the smallest and largest eigenvalues, $\lambda_{1}$ and
$\lambda_{N}$, respectively. This motivates the definition

$C_{r}=R^{\text{T}}\cdot D_{r}\cdot R\,,$ (25)

where $D_{r}$ is a diagonal matrix defined as

$D_{r}=\text{diag}\Big{[}\text{max}[\lambda_{1},r\lambda_{N}],\
\text{max}[\lambda_{2},r\lambda_{N}],\ \cdots,\
\text{max}[\lambda_{N},r\lambda_{N}]\Big{]}\,.$ (26)

Note, if $r=0$, then $D=D_{r}\ \Longrightarrow\ C=C_{r}$. As this parameter is
increased, the lowest eigenvalues are adjusted to some fixed fraction of the
largest value. This approach defines a new test statistic and, in principle,
one can sample its corresponding distribution to define $p$-values and assess
the quality of fits. This goes beyond the scope of this work and we only
perform the modified fits as a cross check to show that the extracted fit
parameters are robust under these regularizations of the covariance matrix.
Such fits are reported in the tables of this section with labels of the form
$[r=0.01\ (n/m)]$ where $r$ indicates the adjustment parameter and $(n/m)$
indicates the number of eigenvalues that have been changed versus the total
number. We have additionally re-done many of the fits summarized Tables 1-4
with the pion mass, $m_{\pi}$, and the anisotropy, $\xi$, varied by one
standard deviation. We find in all cases that the effect of this shift is well
below the statistical uncertainties on the extracted fit parameters.

Fit |  $E^{\star}_{3,\text{cut}}$ |  $p\cot\delta(p)$ |  fit result |  $\chi^{2}/\text{dof}$
---|---|---|---|---
A3(K2) | $4.4m_{\pi}$ |  $-1/a_{0}$ |  $m_{\pi}a_{0}=0.293\pm 0.011$ |  $31.5/(16-1)=2.10$
B3(K2) $[r=0.01\ (5/16)]$ | $4.4m_{\pi}$ |  $-1/a_{0}$ |  $m_{\pi}a_{0}=0.298\pm 0.014$ |  $24.5/(16-1)=1.63$
Table 2: Summary of fits to $\pi^{+}\pi^{+}\pi^{+}$ finite-volume energies, for $p\cot\delta(p)=-1/a_{0}$ with $\mathcal{K}_{3,{\sf iso}}=0$ fixed. All values of total momentum $\bm{P}$ (from $[000]$ to $[111]$) and both volumes ($20^{3}$ and $24^{3}$) are used in each fit. Columns as in Table 1, with “$E^{\star}_{3,\text{cut}}$” indicating the three-particle center-of-momentum frame energy cutoff.  Fit |  $p\cot\delta(p)$ |  $E^{\star}_{3,\text{cut}}$ |  $\mathcal{K}_{3,{\sf iso}}$ |  fit result | $\chi^{2}/\text{dof}$
---|---|---|---|---|---
A3(K3) |  $-1/a_{0}$ $($B${}_{2})$ | $4.4m_{\pi}$ |  $c_{1}/m_{\pi}^{2}$ |  $c_{1}=-253\pm 874$ |  $31.4/(16-1)=2.10$
B3(K3) |  $-1/a_{0}$ $($B${}_{2})$ | $4.4m_{\pi}$ |  $c_{1}/m_{\pi}^{2}+c_{2}\Delta/m_{\pi}^{4}$ |  $\begin{aligned} c_{1}&=5039\pm 1731\\\\[-3.0pt] c_{2}&=-637\pm 92\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0&\,-0.8\\\\[-3.0pt] \\!\phantom{1.0}&\phantom{++}\\!1.0\end{aligned}\bigg{]}$ |  $25.9/(16-2)=1.85$
C3(K3) $[r=0.01\ (5/16)]$ |  $-1/a_{0}$ $($B${}_{2})$ | $4.4m_{\pi}$ |  $c_{1}/m_{\pi}^{2}$ |  $c_{1}=-524\pm 892$ |  $24.4/(16-1)=1.63$
Table 3: Summary of fits to $\pi^{+}\pi^{+}\pi^{+}$ finite-volume energies, for various choices of $\mathcal{K}_{3,{\sf iso}}$ (with $p\cot\delta(p)$ given by fit B2 of Table 1). Columns as in Tables 1 and 2, with “$\mathcal{K}_{3,{\sf iso}}$” indicating the fit function used and $\Delta\equiv E_{3}^{\star 2}-9m_{\pi}^{2}$ encoding a linear-dependence in the squared center-of-momentum frame energy.  Fit |  $E^{\star}_{2,\text{cut}}$ |  $E^{\star}_{3,\text{cut}}$ |  $p\cot\delta(p)$ |  $\mathcal{K}_{3,{\sf iso}}$ |  fit result |  $\chi^{2}/\text{dof}$
---|---|---|---|---|---|---
A2+3 | $3.4m_{\pi}$ | $4.4m_{\pi}$ |  $-1/a_{0}$ |  $0$ |  $m_{\pi}a_{0}=0.300\pm 0.007$ |  $64.7/(37-1)=1.80$
B2+3 | $3.4m_{\pi}$ | $4.4m_{\pi}$ |  $-1/a_{0}$ |  $c_{1}/m_{\pi}^{2}$ |  $\begin{aligned} m_{\pi}a_{0}&=0.296\pm 0.008\\\\[-3.0pt] c_{1}&=-339\pm 770\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &0.6\\\\[-3.0pt] &1.0\end{aligned}\bigg{]}$ | $64.5/(37-2)=1.84$
C2+3 $[r=0.005\ (11/37)]$ | $3.4m_{\pi}$ | $4.4m_{\pi}$ |  $-1/a_{0}$ |  $0$ |  $m_{\pi}a_{0}=0.297\pm 0.008$ |  $50.9/(37-1)=1.42$
D2+3 $[r=0.005\ (11/37)]$ | $3.4m_{\pi}$ | $4.4m_{\pi}$ |  $-1/a_{0}$ |  $c_{1}/m_{\pi}^{2}$ |  $\begin{aligned} m_{\pi}a_{0}&=0.293\pm 0.010\\\\[-3.0pt] c_{1}&=-426\pm 814\end{aligned}\ \bigg{[}\begin{aligned} \\!1.0\ &0.7\\\\[-3.0pt] &1.0\end{aligned}\bigg{]}$ | $50.7/(37-2)=1.45$
Table 4: Summary of combined fits to both two- and three-pion energies.
Columns as in Tables 1-3.

###  .3 Details of the three-particle integral equations

In this subsection we give additional details concerning Eq. (5) of the main
text and give details on its numerical implementation. We begin by reviewing
the integral equations presented in Ref. Hansen and Sharpe (2015). As the
results of the fits summarized in the main text (and detailed in the previous
subsection) are consistent with $\mathcal{K}_{3,{\sf iso}}=0$, we focus here
on the case of a weak three-body interaction, keeping only the linear
contribution in this term. We begin with the unsymmetrized three-body
scattering amplitude

$\displaystyle\mathcal{M}^{(u,u)}_{3}(\bm{p},\bm{k})$
$\displaystyle=\mathcal{D}^{(u,u)}(\bm{p},\bm{k})+\mathcal{E}^{(u)}(\bm{p})\mathcal{K}_{3,{\sf
iso}}\mathcal{E}^{(u)}(\bm{k})+\mathcal{O}(\mathcal{K}_{3,{\sf iso}}^{2})\,,$
(27)

where the superscripts indicate the lack of exchange symmetry, and $\bm{k}$
and $\bm{p}$ specify the momenta of the so-called spectator particles in the
initial and final state, respectively. In general, the factors appearing in
Eq. (27) carry angular momentum indices. However, as discussed in the main
text, the $\pi^{+}\pi^{+}$ system at low energies is dominated by the $S$-wave
component. Thus, we restrict attention here to $\pi^{+}\pi^{+}$ with zero
angular momentum, such that $\mathcal{M}^{(u,u)}_{3}(\bm{p},\bm{k})$ and
$\mathcal{D}^{(u,u)}(\bm{p},\bm{k})$ are simple functions, with no implicit
indices.

In this limiting case, $\mathcal{D}^{(u,u)}(\bm{p},\bm{k})$ satisfies the
implicit equation

$\displaystyle\mathcal{D}^{(u,u)}(\bm{p},\bm{k})$
$\displaystyle=-\mathcal{M}_{2}(E_{2,p}^{\star})G^{\infty}(\bm{p},\bm{k})\mathcal{M}_{2}(E_{2,k}^{\star})-\mathcal{M}_{2}(E_{2,p}^{\star})\int\\!\frac{d^{3}\bm{k}^{\prime}}{(2\pi)^{3}2\omega_{k^{\prime}}}G^{\infty}(\bm{p},\bm{k}^{\prime})\mathcal{D}^{(u,u)}(\bm{k}^{\prime},\bm{k})\,,$
(28)

where $E_{2,k}^{\star 2}\equiv(E-\omega_{k})^{2}-k^{2}$ is the center-of-
momentum energy for the non-spectator pair (with $k=|\bm{k}|$ defined in the
three-particle zero-momentum frame). In words, the unsymmetrized amplitude
$\mathcal{D}^{(u,u)}$ can be evaluated by solving an integral equation
depending only on the two-particle scattering amplitude
$\mathcal{M}_{2}(E_{2,p}^{\star})$ and the exchange propagator

$G^{\infty}(\bm{p},\bm{k})\equiv\frac{H(p,k)}{b_{pk}^{2}-m^{2}+i\epsilon}\,.$
(29)

Here we have suppressed the $\pi$ subscript on the mass and have introduced
$b^{2}_{pk}\equiv(E^{\star}-\omega_{p}-\omega_{k})^{2}-(\bm{p}+\bm{k})^{2}$ as
well as $H(p,k)$, a cut-off function built into the relation between finite-
volume energies and $\mathcal{K}_{3,{\sf iso}}$, as well as that between
$\mathcal{K}_{3,{\sf iso}}$ and the physical scattering amplitude. The
definition used here is

$H(p,k)\equiv J\big{(}E_{2,k}^{\star
2}/[2m]^{2}\big{)}\,J\big{(}E_{2,p}^{\star 2}/[2m]^{2}\big{)}\,,\qquad\qquad
J(x)\equiv\begin{cases}0\,,&x\leq 0\,;\\\
\exp\left(-\frac{1}{x}\exp\left[-\frac{1}{1-x}\right]\right)\,,&0<x\leq
1\,;\\\ 1\,,&1<x\,.\end{cases}$ (30)

(See also Refs. Hansen and Sharpe (2014, 2015) for more discussion on this
technical aspect.) Finally $\mathcal{E}^{(u)}(\bm{p})$ is a function closely
related to $\mathcal{D}^{(u,u)}(\bm{p},\bm{k})$ and defined, in a specific
limiting case, in Eq. (38) below. In this work it is only relevant to
demonstrate that $\mathcal{K}_{3,{\sf iso}}$ contributes negligibly to both
the central value and uncertainty of the three-hadron scattering amplitude, as
we describe in the following section.

The next step is to project the remaining directional freedom within
$\mathcal{D}^{(u,u)}$ onto vanishing three-particle angular momentum, i.e. to
$J=0$. Defining

$\displaystyle\mathcal{D}^{(u,u)}_{{\sf s}}(p,k)$
$\displaystyle\equiv\int\frac{d\Omega_{\hat{k}}}{4\pi}\frac{\,d\Omega_{\hat{p}}}{4\pi}\mathcal{D}^{(u,u)}(\bm{p},\bm{k})\,,$
(31)

one can show that $\mathcal{D}_{\sf s}$ satisfies a one-dimensional integral
equation of the form

$\displaystyle\mathcal{D}^{(u,u)}_{{\sf s}}(p,k)$
$\displaystyle=-\mathcal{M}_{2}(E_{2,p}^{\star})G_{\sf
s}(p,k,\epsilon)\mathcal{M}_{2}(E_{2,k}^{\star})-\mathcal{M}_{2}(E_{2,p}^{\star})\int_{0}^{k_{\sf
max}}\frac{k^{\prime 2}\,dk^{\prime}}{(2\pi)^{2}\omega_{k^{\prime}}}\,G_{\sf
s}(p,k^{\prime},\epsilon)\mathcal{D}^{(u,u)}_{\sf s}(k^{\prime},k),$ (32)

where

$\displaystyle G_{\sf s}(p,k,\epsilon)$
$\displaystyle\equiv\int\frac{d\Omega_{\hat{p}}}{4\pi}\frac{\,d\Omega_{\hat{k}}}{4\pi}G^{\infty}(\bm{p},\bm{k})=-\frac{H(p,k)}{4pk}\log\left[\frac{2pk-(E-\omega_{k}-\omega_{p})^{2}+p^{2}+k^{2}+m^{2}-i\epsilon}{-2pk-(E-\omega_{k}-\omega_{p})^{2}+p^{2}+k^{2}+m^{2}-i\epsilon}\right]\,.$
(33)

In Fig. 6 we plot $G_{\sf s}$ for a range of kinematic values. Setting
$\mathcal{K}_{3,{\sf iso}}=0$ and combining Eqs. (27), (32) and (33), we
arrive at Eq. (5) of the main text. Note that, in Eq. (32), we have included
an explicit cutoff at $k_{\sf max}=(s+m^{2})/(2\sqrt{s})$. This is done
without any additional approximation as $H(p,k)$ has vanishing support for
$k>k_{\sf max}$.

We comment here that, in general, projecting $\mathcal{M}_{3}^{(u,u)}$ to
definite partial waves is not equivalent to doing the same for the symmetrized
three-hadron scattering amplitude. However, as discussed around Eqs. (15) and
(16) of Ref. Jackura  _et al._ (2019b) any definite-angular-momentum component
of the the fully symmetric $\mathcal{M}_{3}$ can be assembled from known
combinations of the various angular-momentum components of
$\mathcal{M}_{3}^{(u,u)}$. Thus the latter form a basis for constructing the
physical three-hadron amplitudes.

Figure 6: Plots of the kernel function $G_{\sf s}(p,k,\epsilon)$ [defined Eq.
(33)] vs $E_{2,k}^{\star}$, here for four different energies as shown in the
legend and with $p=0.1m_{\pi}$ and $\epsilon\to 0^{+}$. For
$\text{Im}[p\,k\,G_{\sf s}(p,k,\epsilon)]$, in the right panel, we have
vertically off-set the curves for improved readability.

To solve Eq. (32) numerically we replace the integral $\int_{0}^{k_{\sf
max}}dk^{\prime}$ with a discrete sum $\sum_{k^{\prime}}\Delta k$ containing
$N$ terms. Then a discretized version of the equation can be written in a
matrix form

$\bm{D}(N,\epsilon)=-\bm{\mathcal{M}}\cdot\bm{G}(\epsilon)\cdot\bm{\mathcal{M}}-\bm{\mathcal{M}}\cdot\bm{G}(\epsilon)\cdot\bm{P}\cdot\bm{D}(N,\epsilon)\,,$
(34)

where we have introduced the following $N\times N$ matrix representations

$\displaystyle\bm{G}_{pk}(\epsilon)=G_{\sf
s}(p,k,\epsilon)\,,\qquad\bm{\mathcal{M}}_{k^{\prime}k}=\delta_{k^{\prime}k}\,\mathcal{M}_{2}(E_{2,k}^{\star})\,,\qquad\bm{P}_{k^{\prime}k}=\delta_{k^{\prime}k}\,\frac{k^{2}\Delta
k}{(2\pi)^{2}\omega_{k}}\,,$ (35)

as well as $\bm{D}(N,\epsilon)$, which becomes our target quantity in the
ordered double limit

$\mathcal{D}^{(u,u)}_{\sf s}(p,k)=\lim_{\epsilon\to
0}\lim_{N\to\infty}\bm{D}_{pk}(N,\epsilon)\,.$ (36)

Here, in a slight abuse of notation, the indices $p\,k$ represent the choices
that are closest to physical momenta $p$ and $k$ for a given $N$ value. Eq.
(34) can then be solved through a matrix inverse to yield

$\bm{D}(N,\epsilon)=-\big{[}\mathbb{I}+\bm{\mathcal{M}}\cdot\bm{G}(\epsilon)\cdot\bm{P}\big{]}^{-1}\cdot\bm{\mathcal{M}}\cdot\bm{G}(\epsilon)\cdot\bm{\mathcal{M}}\,.$
(37)

Finally we return to the endcap factors appearing on either side of
$\mathcal{K}_{3,{\sf iso}}$ in Eq. (27). These are defined in Eqs. (105) and
(106) of Ref. Hansen and Sharpe (2014) and in the overall $S$-wave
approximation they take the form

$\mathcal{E}^{(u)}_{\sf
s}(p)=\frac{1}{3}-\mathcal{M}_{2}(E_{2,p}^{\star})\rho(E_{2,p}^{\star})-\int\\!\frac{dp^{\prime}p^{\prime
2}}{(2\pi)^{2}\omega_{p^{\prime}}}\mathcal{D}_{\sf
s}^{(u,u)}(p,p^{\prime})\rho(E_{2,p^{\prime}}^{\star})\,,$ (38)

where

$\displaystyle\rho(E_{2,p}^{\star})$
$\displaystyle\equiv-i\frac{J\big{(}E_{2,k}^{\star 2}/[2m]^{2}\big{)}}{16\pi
E_{2,k}^{\star}}\sqrt{E_{2,p}^{\star 2}/4-m^{2}}\,,$ (39)

and it is understood that the the $i\sqrt{-x}$ branch is taken for $x<0$.

As $\mathcal{K}_{3,{\sf iso}}$ is consistent with zero in all fits we have
considered, the main purpose in keeping track of these quantities is to
estimate the propagation of uncertainties into $\mathcal{M}_{3}$. In
particular we note

$\Delta\mathcal{M}_{3}^{2}=\sum_{\eta,\eta^{\prime}}[\partial_{\eta}\mathcal{M}_{3}]\,C_{\eta\eta^{\prime}}\,[\partial_{\eta^{\prime}}\mathcal{M}_{3}]\,,$
(40)

where $C_{\eta\eta^{\prime}}$ represents the fit-parameter covariance matrix
and the sums run over all inputs to the scattering amplitude, in particular
the scattering length $a_{0}$ as well as $\mathcal{K}_{3,{\sf iso}}$. In the
next section we present numerical solutions for $\mathcal{M}_{3}$ based on the
parameters extracted from the finite-volume energies.

###  .4 Solutions for $\mathcal{M}_{3}$

Figure 7: Real (left) and imaginary (right) parts of
$p\,k\,\mathcal{M}^{(u,u)}_{\sf s}$, determined by solving Eq. (37) using
$m_{\pi}a_{0}=0.296$ and $\mathcal{K}_{3,{\sf iso}}=0$ for the central values.
As we explain in the text, the uncertainties here follow from propagating the
uncertainties on $\mathcal{K}_{3,{\sf iso}}$ and $m_{\pi}a_{0}$, taken from
Fit B2+3 in Table 4.

Figure 7 shows the result of solving Eq. (37) and varying $N$ and $\epsilon$
in order to ensure that the ordered double limit has been saturated.
Symmetrizing this function over the three incoming and three outgoing momenta
yields the full three-hadron scattering amplitude, plotted in Fig. 3 of the
main text.

To assign an uncertainty estimate to $\mathcal{M}^{(u,u)}_{\sf s}(p,k)$ we
have applied a straightforward adaptation of Eq. (40). Specializing to the
case of only $m_{\pi}a_{0}$ and $\mathcal{K}_{3,{\sf iso}}$ as input
parameters, and neglecting the correlations between them, this becomes

$\big{(}\Delta\mathcal{M}^{(u,u)}_{\sf
s}(p,k)\big{)}^{2}=\big{(}\Delta[m_{\pi}a_{0}]\big{)}^{2}\bigg{(}\frac{\partial\mathcal{M}^{(u,u)}_{\sf
s}(p,k)}{\partial(m_{\pi}a_{0})}\bigg{)}^{2}+\big{(}\Delta\mathcal{K}_{3,{\sf
iso}}/9\big{)}^{2}+\mathcal{O}(\Delta\mathcal{K}_{3,{\sf
iso}}\rho\mathcal{M}_{2})^{2}\,,$ (41)

where $\Delta(m_{\pi}a_{0})$ and $\Delta\mathcal{K}_{3,{\sf iso}}$ denote the
uncertainties on the input parameters and $\Delta\mathcal{M}^{(u,u)}_{\sf
s}(p,k)$ is the resulting amplitude uncertainty plotted in the figure. In
practice one finds that $\Delta\mathcal{K}_{3,{\sf iso}}$ contributes
negligibly to the overall uncertainty, simply because the series of
$\mathcal{M}_{2}$ and $G_{\sf s}$ insertions dominates the value of the
amplitude for weakly-interacting systems. For this reason we have taken the
leading part of $\mathcal{E}_{\sf s}^{(u)}$, resulting in the
$\Delta\mathcal{K}_{3,{\sf iso}}/9$ term. The additional corrections,
indicated by the final term, are negligible for this system.

The dominant source of uncertainty enters through the scattering length, i.e.
the first term in Eq. (41). To estimate this, we have numerically performed
the derivative with respect to $m_{\pi}a_{0}$. Since $\mathcal{M}^{(u,u)}_{\sf
s}(p,k)$ is dominated by the contribution proportional to
$\mathcal{M}_{2}^{2}$, in this case the overall uncertainty is well
approximated by

$\frac{\Delta\mathcal{M}^{(u,u)}_{\sf s}(p,k)}{\mathcal{M}^{(u,u)}_{\sf
s}(p,k)}=2\frac{\Delta[m_{\pi}a_{0}]}{m_{\pi}a_{0}}\,.$ (42)

To produce Fig. 7 we have used $\Delta[m_{\pi}a_{0}]=0.016$ and
$\Delta\mathcal{K}_{3,{\sf iso}}=770$. These are taken from B2+3 in Table 4
with the uncertainty on $m_{\pi}a_{0}$ doubled to account for systematic
variations between the various fits that have been performed.

###  .5 Operator construction and lists

Following Ref. Dudek  _et al._ (2012), to determine the $\pi\pi$ $I=2$ finite-
volume energies we compute correlation functions featuring operators
constructed to resemble a $\pi\pi$ structure. Schematically, for an operator
in lattice irrep $\Lambda_{12}$ and row $\mu_{12}$ with overall momentum
$\bm{k_{12}}$,

$(\pi\pi)^{[\bm{k_{1}},\bm{k_{2}}]\dagger}_{\Lambda_{12}\,\mu_{12}}(\bm{k_{12}})=\sum\limits_{\begin{subarray}{c}\bm{k_{1}},\,\bm{k_{2}}\\\
\bm{k_{1}}+\bm{k_{2}}=\bm{k_{12}}\end{subarray}}\mathcal{C}(\bm{k_{12}},\Lambda_{12},\mu_{12};\bm{k_{1}},\Lambda_{1};\bm{k_{2}},\Lambda_{2})\;\pi^{\dagger}_{\Lambda_{1}}(\bm{k_{1}})\;\pi^{\dagger}_{\Lambda_{2}}(\bm{k_{2}})\,,$
(43)

where the sum is over all momenta related to $\bm{k_{1}}$ and $\bm{k_{2}}$ by
allowed lattice rotations, $\mathcal{C}$ is an appropriate generalized
Clebsch-Gordan coefficient for $\Lambda_{1}\otimes\Lambda_{2}\to\Lambda_{12}$,
and flavor indices and the projection onto $I=2$ are not written explicitly.
Here $\pi_{\Lambda_{i}}^{\dagger}(\bm{k_{i}})$ is the optimal linear
combination of operators to interpolate a $\pi$ with momentum $\bm{k_{i}}$ in
irrep888These are all one dimensional and so we omit the irrep row index.
$\Lambda_{i}$ using a basis of fermion-bilinear operators featuring various
Dirac $\gamma$ matrices and gauge-covariant derivatives – see Ref. Dudek  _et
al._ (2012) for details. In this work the basis of fermion-bilinear operators
used for a $\pi$ operator has up to three derivatives for $\pi$ at rest and up
to one derivative for $\pi$ at non-zero momentum, except we use up to two
derivatives for $1\leq|\bm{k_{i}}|^{2}\leq 4$ on the $24^{3}$ volume. The
operators used to compute the $\pi\pi$ spectra shown in Fig. 1 of the main
text are listed in Table 5.

$\bm{P}$ $\Lambda$ | $L/a_{s}=20$ | $L/a_{s}=24$
---|---|---
[000] $A_{1}^{+}$ |  $\pi_{[000]}\pi_{[000]}$, $\pi_{[100]}\pi_{[100]}$, $\pi_{[110]}\pi_{[110]}$, $\pi_{[111]}\pi_{[111]}$, |  $\pi_{[000]}\pi_{[000]}$, $\pi_{[100]}\pi_{[100]}$, $\pi_{[110]}\pi_{[110]}$, $\pi_{[111]}\pi_{[111]}$,
| $\pi_{[200]}\pi_{[200]}$ | $\pi_{[200]}\pi_{[200]}$
[100] $A_{1}$ |  $\pi_{[000]}\pi_{[100]}$, $\pi_{[100]}\pi_{[110]}$, $\pi_{[110]}\pi_{[111]}$, $\pi_{[100]}\pi_{[200]}$ |  $\pi_{[000]}\pi_{[100]}$, $\pi_{[100]}\pi_{[110]}$, $\pi_{[110]}\pi_{[111]}$, $\pi_{[100]}\pi_{[200]}$,
|  |  $\pi_{[110]}\pi_{[210]}$, $\pi_{[200]}\pi_{[210]}$, $\pi_{[111]}\pi_{[211]}$
[110] $A_{1}$ |  $\pi_{[000]}\pi_{[110]}$, $\pi_{[100]}\pi_{[100]}$, $\pi_{[100]}\pi_{[111]}$, $\pi_{[110]}\pi_{[110]}$ |  $\pi_{[000]}\pi_{[110]}$, $\pi_{[100]}\pi_{[100]}$, $\pi_{[100]}\pi_{[111]}$, $\pi_{[110]}\pi_{[110]}$,
|  |  $\pi_{[110]}\pi_{[200]}$, $\pi_{[100]}\pi_{[210]}$, $\pi_{[111]}\pi_{[210]}$, $\pi_{[110]}\pi_{[211]}$
[111] $A_{1}$ |  $\pi_{[000]}\pi_{[111]}$, $\pi_{[100]}\pi_{[110]}$ |  $\pi_{[000]}\pi_{[111]}$, $\pi_{[100]}\pi_{[110]}$, $\pi_{[111]}\pi_{[200]}$, $\pi_{[110]}\pi_{[210]}$,
|  | $\pi_{[100]}\pi_{[211]}$
[200] $A_{1}$ |  $\pi_{[100]}\pi_{[100]}$, $\pi_{[000]}\pi_{[200]}$, $\pi_{[110]}\pi_{[110]}$, $\pi_{[111]}\pi_{[111]}$ |  $\pi_{[100]}\pi_{[100]}$, $\pi_{[000]}\pi_{[200]}$, $\pi_{[110]}\pi_{[110]}$, $\pi_{[100]}\pi_{[210]}$,
|  |  $\pi_{[111]}\pi_{[111]}$, $\pi_{[110]}\pi_{[211]}$, $\pi_{[210]}\pi_{[210]}$
Table 5: The $\pi\pi$ $I=2$ operators, $\pi_{\bm{k_{1}}}\pi_{\bm{k_{2}}}$,
used to compute the finite-volume energy levels shown in Fig. 1 of the main
text (upper plots) in irrep $\Lambda$ with overall momentum $\bm{P}$. These
are constructed from optimized $\pi$ operators with momentum types
$\bm{k_{1}}$ and $\bm{k_{2}}$; different momentum directions are summed over
as in Eq. (43). Momenta are displayed using the shorthand notation
$[ijk]=\frac{2\pi}{L}(i,j,k)$.

In a similar way, operators used to compute $\pi\pi\pi$ $I=3$ energies
resemble a $\pi\pi\pi$ structure and are formed by combining a $\pi\pi$ $I=2$
operator with a $\pi$ operator, as detailed in Ref. Woss  _et al._ (2019).
Schematically, for an operator in lattice irrep $\Lambda$ and row $\mu$ with
overall momentum $\bm{P}$,

$(\pi\pi\pi)^{[\bm{k_{12}}[\bm{k_{1}},\bm{k_{2}}],\bm{k_{3}}]\dagger}_{\Lambda\,\mu}(\bm{P})=\sum\limits_{\begin{subarray}{c}\bm{k_{12}},\,\bm{k_{3}}\\\
\bm{k_{12}}+\bm{k_{3}}=\bm{P}\end{subarray}}\mathcal{C}(\bm{P},\Lambda,\mu;\bm{k_{12}},\Lambda_{12},\mu_{12};\bm{k_{3}},\Lambda_{3})\;(\pi\pi)^{[\bm{k_{1}},\bm{k_{2}}]\dagger}_{\Lambda_{12}\,\mu_{12}}(\bm{k_{12}})\;\pi^{\dagger}_{\Lambda_{3}}(\bm{k_{3}})\,,$
(44)

where the sum is over all momenta related to $\bm{k_{12}}$ and $\bm{k_{3}}$ by
allowed lattice rotations and, again, flavor indices and the projection onto
$I=3$ are not written explicitly.

From Bose symmetry, a $\pi\pi\pi$ system must be symmetric under the
interchange of any pair of pions. The pions have no intrinsic spin and we are
considering $I=3$ which means that the flavor structure is symmetric under
interchange of any pair; therefore, the spatial structure must also be
symmetric under the interchange of any pair of pions. The operator
construction in Eq. (44) gives operators with the correct symmetry properties
because the three pions are identical, but it does not make this symmetry
manifest and two different sets of ($|\bm{k_{1}}|$, $|\bm{k_{2}}|$,
$|\bm{k_{3}}|$, $|\bm{k_{12}}|$, $\Lambda_{12}$) may lead to equivalent
operators, or a number of different sets may give linearly-dependent
operators.999Strictly we mean the types of momenta (i.e. the equivalence class
of momenta related by rotations in the octahedral group or little group)
rather than the magnitudes, but there is no distinction for the momenta we are
considering here. To ensure we have an appropriate set of independent
operators, for each $\bm{P}$ and $\Lambda$, we construct an operator for every
possible set ($|\bm{k_{1}}|$, $|\bm{k_{2}}|$, $|\bm{k_{3}}|$, $|\bm{k_{12}}|$,
$\Lambda_{12}$) with $|\bm{k_{1}}|^{2}+|\bm{k_{2}}|^{2}+|\bm{k_{3}}|^{2}$ less
than some cutoff. We write Eq. (44) schematically as,

$(\pi\pi\pi)_{\Lambda\mu}(\bm{P})=\sum\widetilde{\mathcal{C}}(\bm{k_{1}},\bm{k_{2}},\bm{k_{3}})\,\pi^{\dagger}(\bm{k_{1}})\pi^{\dagger}(\bm{k_{2}})\pi^{\dagger}(\bm{k_{3}})=\widetilde{\mathcal{C}}\cdot\mathcal{V}_{[\pi\pi\pi]}\,,$
(45)

where $\widetilde{\mathcal{C}}$ with no labels represents a row vector of
coefficients and $\mathcal{V}_{[\pi\pi\pi]}$ a column of operators, such that
the dot-product reproduces the sum. We then introduce the matrix $R(ijk)$
which acts on $\mathcal{V}_{[\pi\pi\pi]}$ by mapping a given entry
$\pi^{\dagger}(\bm{p}_{1})\pi^{\dagger}(\bm{p}_{2})\pi^{\dagger}(\bm{p}_{3})$
into
$\pi^{\dagger}(\bm{p}_{i})\pi^{\dagger}(\bm{p}_{j})\pi^{\dagger}(\bm{p}_{k})$.
This allows us to define the symmetrized vector of Clebsch-Gordan coefficients

$\widetilde{\mathcal{C}}\cdot R(123)+\widetilde{\mathcal{C}}\cdot
R(231)+\widetilde{\mathcal{C}}\cdot R(312)+\widetilde{\mathcal{C}}\cdot
R(132)+\widetilde{\mathcal{C}}\cdot R(213)+\widetilde{\mathcal{C}}\cdot
R(321)\,,$ (46)

where $R(123)$ is just the identity matrix. The final step is to check that
the resulting vector is non-zero and linearly independent from the analogous
expressions for the already considered operators. The resulting sets of
independent operators used in this work are listed in Tables 6, 7, 8, 9 and
10.101010This is not always a unique choice and any independent set of
operators could be used to achieve the same results. In extracting the
energies, subsets of the full operator set were also considered, in order to
investigate sensitivity to the detailed choice of operator basis.

Finally, in order to give further intuition into the operators used, we also
include a diagrammatic representation of the individual pion momentum
assignments in the tables. The diagrams portray the integer vectors
$\bm{d}_{1}$, $\bm{d}_{2}$, $\bm{d}_{3}$, each given by
$\bm{d}_{i}=L\bm{k}_{i}/(2\pi)$. The vectors are assigned a color (orange and
green for the first two pions, and blue for the third) and the absence of any
given color corresponds to a vector of magnitude zero. As summarized by Eq.
(44), our operator construction is based on combining two-pion operators in a
definite irrep with the third pion. An alternative basis is given by summing a
given momentum assignment, represented by a given set
$\bm{d}_{1}\bm{d}_{2}\bm{d}_{3}$, over all rotations in the octahedral group
(in the case of $\bm{P}=[000]$) or else a little group thereof (for non-zero
total momentum) weighted by the appropriate Clebsch-Gordan coefficients. The
operators reached via this alternative construction are equal to a linear
combination of those given by Eq. (44).

For example, on the third line of Table 7, two distinct momentum assignments
arise from combining the $\pi\pi_{[000]A_{1}^{+}}$ with the third
$\pi_{[100]}$ operator. In this case the $\pi\pi_{[000]A_{1}^{+}}$ is built
from individual pions with a unit of back-to-back momentum. When one sums over
the coefficients projecting onto $A_{1}^{+}$, contributions arise with the
back-to-back axis both aligned and perpendicular to the total momentum
direction, $\bm{P}=[001]$. Thus, momentum assignments corresponding to both
diagrams shown in line 3 of Table 7 contribute to the operator on that line.
By contrast, on line 4 of Table 7 only a single momentum configuration
contributes, as indicated. This implies that the two operators are
independent, since they are built from independent linear combinations of the
two momentum configurations.

Operators 8 through 11 of Table 7 give a more complicated example. The first
two (8 and 9) correspond to two linear combinations of two configuration
types, labeled with subscripts 1 and 3, and the next two (10 and 11) are equal
to linear combinations of the operators labeled 1, 2, and 4.

The diagrams in Tables 6, 7, 8, 9 and 10 provide a cross check on the linear-
independence of the operators. Each row corresponds to a linear combination of
the displayed momentum configurations and the number of linearly independent
operators is equal to the number of distinct diagrams.

|  $\bm{d_{1}^{2}\,d_{2}^{2}\,d_{3}^{2}}$ | $\pi_{\bm{k_{1}}}$ |  $\pi_{\bm{k_{2}}}$ |  $\pi\pi_{\bm{k_{12}}\Lambda_{12}}$ ($I=2$) | $\pi_{\bm{k_{3}}}$ |  momentum configurations
---|---|---|---|---|---|---
1 |  000 | $\ \pi_{[000]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[000]}\ $ |  |
---
2 |  011 | $\ \pi_{[100]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |
---
3 |  022 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
4 |  112 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
5 |  033 | $\ \pi_{[111]}\ $ | $\ \pi_{[111]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[000]}\ $ |  |
---
6 |  114 | $\ \pi_{[200]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |
---
7 |  123 | $\ \pi_{[111]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
8 |  222 | $\ \pi_{[110]}\ $ | $\ \pi_{[110]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
9 |  044 | $\ \pi_{[200]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[200]A_{1}}\ $ | $\ \pi_{[200]}\ $ |  |
---
10 |  125 | $\ \pi_{[210]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
Table 6: The $\pi\pi\pi$ $I=3$ operators used to compute the finite-volume energy levels shown in Fig. 1 of the main text (lower left plot) in irrep $A_{1}^{-}$ with overall momentum $\bm{P}=[000]$, labeled by $\bm{d}_{1}^{2}\,\bm{d}_{2}^{2}\,\bm{d}_{3}^{2}$ where $\bm{d_{i}}=\bm{k_{i}}(L/2\pi)$. Different momentum directions for the momentum types $\bm{k_{1}}$, $\bm{k_{2}}$, $\bm{k_{3}}$ and $\bm{k_{12}}$ are summed over as in Eq. (44). Operators 1 to 8 are used on the $20^{3}$ volume and operators 1 to 10 are used on the $24^{3}$ volume. The momentum configuration diagrams in the rightmost column are explained in the text. Operators separated by a single horizontal line correspond to states that are degenerate in the non-relativistic, non-interacting theory. |  $\bm{d_{1}^{2}\,d_{2}^{2}\,d_{3}^{2}}$ | $\pi_{\bm{k_{1}}}$ |  $\pi_{\bm{k_{2}}}$ |  $\pi\pi_{\bm{k_{12}}\Lambda_{12}}$ ($I=2$) | $\pi_{\bm{k_{3}}}$ |  momentum configurations
---|---|---|---|---|---|---
1 |  001 | $\ \pi_{[000]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[100]}\ $ |  |
---
2 |  012 | $\ \pi_{[100]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
3 |  111 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[100]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!\\!2}$
---
4 |  111 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
5 |  014 | $\ \pi_{[100]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[200]}\ $ |  |
---
6 |  023 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[111]}\ $ |  |
---
7 |  113 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[111]}\ $ |  |
---
8 |  122 | $\ \pi_{[110]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[111]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  ${}^{\\!\\!\\!\\!\\!\\!1}$ ${}^{\\!\\!\\!\\!\\!\\!3}$
9 |  122 | $\ \pi_{[110]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[111]E_{2}}\ $ | $\ \pi_{[110]}\ $
10 |  122 | $\ \pi_{[110]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  ${}^{\\!\\!\\!\\!\\!\\!1}$ ${}^{\\!\\!\\!\\!\\!\\!2}$ ${}^{\\!\\!\\!\\!\\!\\!4}$
11 |  122 | $\ \pi_{[110]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[100]B_{1}}\ $ | $\ \pi_{[110]}\ $
12 |  025 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[210]}\ $ |  |
---
13 |  115 | $\ \pi_{[210]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[100]}\ $ |
14 |  115 | $\ \pi_{[210]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]B_{1}}\ $ | $\ \pi_{[100]}\ $
15 |  124 | $\ \pi_{[200]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[210]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
16 |  124 | $\ \pi_{[200]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
Table 7: As Table 6 but for the $A_{2}$ irrep with overall momentum $\bm{P}=[001]$. Operators 1 to 8 are used on the $20^{3}$ volume and operators 1 to 16 are used on the $24^{3}$ volume. Operators separated by a single horizontal line correspond to states that are degenerate in the non-relativistic, non-interacting theory. Operators with no horizontal line correspond to states that are degenerate in the relativistic, non-interacting theory but are split by the interactions. |  $\bm{d_{1}^{2}\,d_{2}^{2}\,d_{3}^{2}}$ | $\pi_{\bm{k_{1}}}$ |  $\pi_{\bm{k_{2}}}$ |  $\pi\pi_{\bm{k_{12}}\Lambda_{12}}$ ($I=2$) | $\pi_{\bm{k_{3}}}$ |  momentum configurations
---|---|---|---|---|---|---
1 |  002 | $\ \pi_{[000]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[110]}\ $ |  |
---
2 |  011 | $\ \pi_{[100]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |
---
3 |  013 | $\ \pi_{[100]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[111]}\ $ |  |
---
4 |  022 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
5 |  112 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!2}$
---
6 |  112 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
7 |  112 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[200]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
8 |  112 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[000]E^{+}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!2}$
---
Table 8: As Table 6 but for the $A_{2}$ irrep with overall momentum $\bm{P}=[011]$ (continued in Table 9). Operators 1 to 11 are used on the $20^{3}$ volume and operators 1 to 20 are used on the $24^{3}$ volume. |  $\bm{d_{1}^{2}\,d_{2}^{2}\,d_{3}^{2}}$ | $\pi_{\bm{k_{1}}}$ |  $\pi_{\bm{k_{2}}}$ |  $\pi\pi_{\bm{k_{12}}\Lambda_{12}}$ ($I=2$) | $\pi_{\bm{k_{3}}}$ |  momentum configurations
---|---|---|---|---|---|---
9 |  015 | $\ \pi_{[100]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[210]}\ $ |  |
---
10 |  024 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[200]}\ $ |  |
---
11 |  114 | $\ \pi_{[200]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[210]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |
---
12 |  123 | $\ \pi_{[111]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[211]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |
---
13 |  123 | $\ \pi_{[111]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!\\!\\!2}$
---
14 |  123 | $\ \pi_{[111]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]B_{2}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!\\!\\!2}$
---
15 |  222 | $\ \pi_{[110]}\ $ | $\ \pi_{[110]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
  ${}^{\\!\\!\\!\\!\\!\\!2}$
---
  ${}^{\\!\\!\\!\\!\\!\\!3}$
---
16 |  222 | $\ \pi_{[110]}\ $ | $\ \pi_{[110]}\ $ | $\ \pi\pi_{[211]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
17 |  222 | $\ \pi_{[110]}\ $ | $\ \pi_{[110]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!\\!\\!4}$
---
18 |  222 | $\ \pi_{[110]}\ $ | $\ \pi_{[110]}\ $ | $\ \pi\pi_{[200]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!4}$
---
${}^{\\!\\!\\!\\!\\!\\!3}$
---
19 |  026 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[211]}\ $ |  |
---
20 |  035 | $\ \pi_{[111]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[111]A_{1}}\ $ | $\ \pi_{[210]}\ $ |  |
---
Table 9: Continuation of Table 8. |  $\bm{d_{1}^{2}\,d_{2}^{2}\,d_{3}^{2}}$ | $\pi_{\bm{k_{1}}}$ |  $\pi_{\bm{k_{2}}}$ |  $\pi\pi_{\bm{k_{12}}\Lambda_{12}}$ ($I=2$) | $\pi_{\bm{k_{3}}}$ |  momentum configurations
---|---|---|---|---|---|---
1 |  003 | $\ \pi_{[000]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[111]}\ $ |  |
---
2 |  012 | $\ \pi_{[110]}\ $ | $\ \pi_{[000]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |
---
3 |  111 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |
---
4 |  113 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[000]A_{1}^{+}}\ $ | $\ \pi_{[111]}\ $ |  |
---
5 |  113 | $\ \pi_{[100]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[200]A_{1}}\ $ | $\ \pi_{[111]}\ $ |  |
---
6 |  122 | $\ \pi_{[110]}\ $ | $\ \pi_{[100]}\ $ | $\ \pi\pi_{[100]A_{1}}\ $ | $\ \pi_{[110]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!1}$
---
${}^{\\!\\!\\!\\!\\!\\!2}$
---
7 |  122 | $\ \pi_{[110]}\ $ | $\ \pi_{[110]}\ $ | $\ \pi\pi_{[110]A_{1}}\ $ | $\ \pi_{[100]}\ $ |  |  ${}^{\\!\\!\\!\\!\\!\\!2}$
---
Table 10: As Table 6 but for the $A_{2}$ irrep with overall momentum
$\bm{P}=[111]$. Operators 1 to 7 are used on both the $20^{3}$ and $24^{3}$
volumes.