Scattering of two and three physical pions at maximal isospin from lattice QCD

We present the first direct Nf=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_f=2$$\end{document} lattice QCD computation of two- and three-π+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi ^+$$\end{document} scattering quantities that includes an ensemble at the physical point. We study the quark mass dependence of the two-pion phase shift, and the three-particle interaction parameters. We also compare to phenomenology and chiral perturbation theory (ChPT). In the two-particle sector, we observe good agreement to the phenomenological fits in s- and d-wave, and obtain Mπa0=-0.0481(86)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_\pi a_0 = -0.0481(86)$$\end{document} at the physical point from a direct computation. In the three-particle sector, we observe reasonable agreement at threshold to the leading order chiral expansion, i.e. a mildly attractive three-particle contact term. In contrast, we observe that the energy-dependent part of the three-particle quasilocal scattering quantity is not well described by leading order ChPT.


Introduction
Quantum chromodynamics (QCD) describes the interaction of quarks and gluons, while only hadrons (mesons and baryons) are experimentally observable. They are low energy bound states, or resonances of the former fundamental particles. Understanding the interactions of two or more hadrons is highly relevant for several reasons. For instance, resonances become visible only when studying the interaction of other hadrons. And for understanding experimental signatures of particle decays, the interactions of the final states need to be understood. a e-mail: fernando.romero@uv.es (corresponding author) Lattice QCD, the formulation of QCD on a spacetime lattice, offers the opportunity of first principles, numerical explorations of few-particle scattering amplitudes. Maybe the most obvious example for the importance of three-particle interactions is the ω-meson, which decays predominantly into three pions with J P = 1 − [1]. Another one would be the Roper resonance [2], with both N π and N ππ decay channels. However, since the investigation of three-particle interactions from lattice QCD is in its infancy, three weakly interacting pions with isospin I = 3 is an interesting and important benchmark system.
The extraction of two-particle scattering amplitudes in Lattice QCD is by now well established for 2 → 2 systems, both theoretically [4][5][6][7][8][9][10][11][12][13][14][15][16], and in practice [3, (see Ref. [41] for a review). One of the most studied systems is isospin-2 ππ scattering. To illustrate the state-of-the-art, we show in Fig. 1 the ππ I = 2 scattering length M π a 0 as a function of M π / f π comparing this work's result to the N f = 2 + 1 + 1 results of Ref. [3]. The new N f = 2 point at a slightly less than physical value of M π / f π as well as the other two new points are compatible within errors with leading order (LO) ChPT (dashed line).
In this article, we present results for scattering quantities of two and three-pion systems with maximal isospin, including for the first time an ensemble at the physical point. This work breaks new ground on several fronts: the first direct computation at the physical point of the I = 2 s-and dwave phase shift, and the chiral dependence of the three-π + quasilocal interaction.

Scattering amplitudes from lattice QCD
The calculation of scattering amplitudes from lattice simulations proceeds in an indirect way. The required physical quantities from the lattice are the finite-volume interacting energies of two and three particles -the finite volume spectrum. The mapping between the finite volume spectrum and infinite-volume scattering quantities -the so-called quantization condition -is known but highly nontrivial. It is valid up to effects that vanish exponentially with the pion mass, ∼ exp(−M π L).
The two-particle quantization condition (QC2) takes the form of a determinant equation [4][5][6] (we assume two identical scalars): Here, F 2 and K 2 are both matrices in angular momentum space , m. The matrix elements of F 2 are kinematical functions (Lüscher zeta function) that depend on the three-momentum of the system, P and the center-of-mass (CM) energy, E * . (K 2 ) m, m = δ m, m (K 2 ) is simply the infinite-volume scattering K-matrix projected to the corresponding partial wave. In order to render the matrices finitedimensional, a truncation must be applied in , by assum-ing that K 2 vanishes for higher partial waves. Furthermore, the relations between K 2 , the phase shift (δ ), and the scattering amplitude (M 2 ) are trivial. More details can be found in Ref. [41]. The three-particle quantization condition (QC3) for identical (pseudo)scalars in the RFT approach reads (G-parity is assumed) [42]: Even though this looks formally identical to Eq. 1, there are some distinct features. First, the matrices in Eq. 2 live in a larger k m space, where , m are the angular momentum indices of the interacting pair, and k labels the threemomentum of the third particle -the spectator. Next, F 3 depends on geometric functions (like F 2 itself), but also on K 2 . Thus, two-particle interactions are a necessary ingredient for three-particle scattering. Note that an analytical continuation of K 2 below threshold is needed for the QC3. Finally, K 3,df is a real, singularity-free, quasilocal, intermediate three-particle scattering quantity -which we aim to determine. As in the case of the QC2, Eq. 2 is infinite-dimensional, and must be truncated. The truncation in k is due to a cut-off function, whereas for , m one assumes that K 3,df vanishes above some value of , see Refs. [42,65] for details. Establishing the connection between K 3,df and the physical scattering amplitude, M 3 requires a set of integral equations, derived in Ref. [43] and solved in Ref. [47]. In this work, we focus only on the extraction of K 3,df . In a finite volume, partial waves mix and, thus, F 2 and F 3 are nondiagonal in , m. The correct labels are then irreducible representations (irreps) of the discrete symmetry group, which we label as . The subduction of angular momenta into irreps is known [76, Table 2]. Therefore, one block-diagonalizes the quantization conditions into irreps, see Refs. [16,49,58,69].

Lattice computation
This work uses N f = 2 flavour lattice QCD ensembles generated by the Extended Twisted Mass collaboration (ETMC) [80], including one ensemble at the physical pion masssee Table 1. For the ensemble generation the Iwasaki gauge action [81] was used together with Wilson clover twisted mass fermions at maximal twist [82]. The latter guarantees scaling towards the continuum with only O(a 2 ) artefacts in the lattice spacing a [83]. The presence of the clover term (with coefficient c sw ) has been shown to further reduce the O(a 2 ) artefacts, in particular isospin-breaking effects of the twisted-mass formulation, which have been empirically found to be very small for masses and decay constants [80]. For the two-pion scattering length with I = 2, discretisation artefacts are only of order O(am q ) 2 , with m q the up/down quark mass [84]. Another possible source of O(a 2 ) effects that should be mentioned is the π 0 contamination in the correlation functions due to the breaking of parity in twisted mass. However, it is also important to realise that at maximal isospin there is no mixing with other flavour states due to broken isospin symmetry. Parametrically, O(a 2 ) artefacts are ∼ 2.5% and O(am q ) 2 ≤ 0.4% for this lattice spacing, and thus well below our statistical uncertainty. The two-and three-π + energy spectrum is measured from Euclidean correlation functions of operators with the corresponding quantum numbers. By means of the single pion operators (π + = −ūγ 5 d), we construct two-particle operators as where p i labels the momentum of each single pion, and similarly for three pions x,y,z e i p 1 x+i p 2 y+i p 3 z × π + (x) π + (y) π + (z).
Correlation functions are computed using the stochastic Laplacian-Heaviside smearing [85,86] with algorithmic parameters as in Ref. [87]. In addition, operators that transform under a specific irrep of a discrete symmetry group are constructed following Ref. [32]. In the two-pion case we use the irreps A (+) 1 , E (+) , B 1 and B 2 , in the three pion channel 1 , E (−) , A 2 , B 1 and B 2 , for all P 2 ≤ 4 with P the centreof-mass momentum. We refer to Table 9 in the appendix for an overview. We extract the spectrum in each irrep independently using the generalized eigenvalue method (GEVM) [6,88,89] and also the GEVM/PGEVM method [90], see the appendix for more details.
A technical issue of lattice calculations with (anti)periodic boundary conditions in the time direction is the presence of so-called thermal states, i.e. effects from states that propagate backwards in time across the boundary. They vanish with M π T → ∞, but at finite values of T , these effects are significant and need to be treated accordingly. In fact, thermal pollutions are one of the major systematic uncertainties in our calculation. We deal with them as follows: using the operators discussed above we build correlator matrices which are input to the GEVM/PGEVM which in turn have so-called principal correlators as output. From the latter energy levels and corresponding error estimates are extracted from bootstrapped, fully correlated fits to the data with fit ranges chosen by eye. We use five different treatments to arrive from a correlator matrix at an energy level. Details of those five treatments are explained in Appendix A1.
As also explained in Appendix A1, the different energy levels per principal correlator (up to five) are then combined using a correlated weighted average. However, to account for the spread between the different methods we use a procedure discussed in Ref. [32] to widen the resampling distribution: for energy level E we compute the scaling factor where δ E is the statistical uncertainty of the weighted average and E Y is the difference between method Y and the weighted average. By scaling the resampling distribution of the weighted average with w, we obtain a distribution that reflects both the statistical and the systematic uncertainties, while still being usable in the bootstrap analysis chain. The energy levels are publicly available [91]. The finite-volume scattering formalism is applicable under the assumption that exponential finite volume effects are negligible. On the physical point ensemble, we have M π L ≈ 3, which implies e −M π L ∼ 5% and might be considered to be at the edge of feasibility. However, based on a ChPT analysis, finite-volume effects are also proportional to [M π /(4π F π )] 2 , which at the physical point reduces finitevolume effects sizably. Moreover, as argued in Ref. [74], if the volume-dependent mass is used to analyze the multiparticle energy levels, the leading finite-size effects cancel. For the other two ensembles we have M π L > 5, which is safe concerning finite volume effects.

Results
In the case of two pions, by keeping only s-wave interactions in A 1 irreps, the projected QC2 becomes a one-to-one corre- spondence of an energy level to a phase shift point [6,7]. For the analysis, we need an appropriate phase shift parametrization. We use a model that incorporates the expected Adler zero [69,94]: with s the center-of-mass energy squared and k 2 = s/4−M 2 π . We will fix the position of the Adler zero to its leading order chiral perturbation theory (LO ChPT) value: z 2 = M 2 π . Even though higher order corrections are to be expected, its value has been seen to be compatible with LO ChPT when left free [92,95,96]. Note that in Eq. 6 with fixed Adler zero, we have M π a 0 = 1/B 0 .
We perform a correlated two-parameter fit to the energy levels. The results for the three ensembles are shown in Table  2. In all cases, the magnitude of the B i coefficients decreases with increasing order, indicating that the expansion converges quickly enough even at the heaviest pion mass. Still, for the heaviest ensemble (cA2.60.32), we also attempt a fit with a quadratic term in k 2 , B 2 and observe a small, barely significant value for B 2 and no substantial change in B 0 and B 1 . Based on ChPT, better convergence is expected for lighter pions.
The s-wave phase shift is visualised for the physical point ensemble in the left panel of Fig. 2. In this plot we also compare to other results in the literature. For the other two ensembles the corresponding plots can be found in the left panels of Figs. 9 and 10, respectively, in the appendix.
One interesting point to discuss is the suitability of the δ 0 parametrization. It has been customary to use a standard effective range expansion parametrization (ERE) for isospin-2 ππ scattering: However, the presence of the Adler zero limits the radius of convergence to k 2 ∼ 0.5M 2 π . For this reason, explicitly incorporating the Adler zero must improve the radius of convergence, and has been shown to provide a better description of the data [69]. Here, we compare again the two fit models. The ERE results are shown in Table 3. As can be seen, the values of χ 2 in the case of the ERE fits are always larger than their Adler-zero counterparts given in Table 2. This further supports the usage of the Adler-zero parametrization for I = 2 ππ scattering.
Similarly, the d-wave phase shift can be obtained from most of the nontrivial irreps when neglecting > 2 waves [15,16]. Since we have few data points, we attempt the following fit (see Table 4): The best fit curve for the physical point ensemble is show in the right panel of Fig. 2 and compared to Ref. [92]. Again, for the other two ensembles the corresponding plots can be found in the appendix in the right panels of Figs. 9 and 10, respectively.
In the three pion case we need to parametrize K 3,df . For this, we expand K 3,df about threshold up to linear terms of relativistic invariants [49]: where K iso,0 df,3 and K iso,1 df,3 are the numerical constants to be determined. This parametrization has no momentum dependence, and thus receives the name "isotropic". It is the threeparticle equivalent of keeping only s-wave interactions. At the next order in the expansion, O( 2 ), three new parameters arise, for which also the d-wave must be included [49]. This is beyond the scope of the present analysis (Table 5).
Following the strategy outlined in Ref. [69], we perform a simultaneous s-wave only fit to two-π + A 1 levels, and all three-π + levels. For this, we use the δ 0 model in Eq. 6 and the K 3,df parametrization in Eq. 9 -four parameters in total, see Table 6. As can be seen the best fit values for B 0 and B 1 agree well between the two-particle and the global fit, with even smaller errors in the case of the latter. For convenience, we provide the full covariance matrices of the fits in Table 6 in the appendix, see Eqs. (B2) to (B4).
We have also performed fits including only the constant term K iso,0 df,3 , the results of which can be found in the appendix. We observe that for the ensembles with larger than physical pion mass value the inclusion of the linear term seems necessary.
In Fig. 11 in the appendix we provide as an example for the physical point ensemble the measured energy spectrum in the two-and three particle sectors separately. In that figure we also compare to the noninteracting energy levels. Moreover, we give the energy levels predicted by our fits, see Tables 2, 4, 6
In Fig. 1 we also compare to results from N f = 2 + 1 + 1 calculations from Ref. [3] and with LO ChPT. Within the uncertainties we do not observe a significant difference between N f = 2 and N f = 2 + 1 + 1 results. Moreover, as was found in all previous investigations of two pions at maximal isospin, LO ChPT describes the mass dependence extraordinarily well. At the physical point, LO ChPT predicts M π a 0 −0.04438, which agrees within error bars with the value we report here, see above. Unfortunately, our determination here suffers from relatively large statistical uncertainties and, thus, cannot compete with determinations based on chiral extrapolations. A summary of various determinations from the literature is compiled in Table 5.
Regarding the d-wave phase shift, we have mild statistical evidence that it is repulsive at the physical point in the considered energy region. We observe agreement within 1σ with Ref. [92], as shown in Fig. 2b. An interesting feature of the phenomenological fits to δ 2 is that there is a sign change near threshold, which yields an attractive phase shift at threshold [92,95,96,106]. We cannot confirm or deny such behaviour, as the explored energy region is too far above threshold. For larger pion mass values, we obtain a similar behaviour. The d-wave phase shift is more repulsive for the two larger pion mass values -see Table 4 and the appendix.
We show our results in the three-particle sector in Fig. 3. As can be seen in Fig. 3a, there is significant evidence that K 3,df at threshold (K iso,0 df,3 ) is positive (attractive). Even though we find reasonable agreement with the LO ChPT [69] prediction, the data suggests that NLO effects can be significant, and it may be worth to extend the ChPT result to one loop in future work. For K iso,1 df,3 , the situation is somewhat different. All evidence points to a negative value, very far from the ChPT results. While one could conclude that a NLO ChPT description is required, there is a subtlety in the LO ChPT prediction: it assumes that the connection between K 3,df and M 3 -which involves integral equations -is trivial in LO ChPT [69]     where M 3,df is the divergence-free three-to-three amplitude [43]. As argued in Ref. [69], this induces large errors in K iso,1 df,3 (up to 50% for 200 MeV pions). The situation is expected to be more dramatic for heavier pions, like our two results at 242 and 340 MeV, for which the largest difference is seen.
In order to address this rigorously, the integral equation must be systematically solved, which is beyond the scope of this work.

Conclusion
We have presented the first N f = 2 lattice calculation of twoand three-π + scattering at the physical point. In the two pion channel we observe very good agreement with other lattice calculations and ChPT or ChPT combined with Roy-Steiner equations for the s-wave phase shift. In particular, for the whole range of pion mass values we have available here we do not observe a significant deviation from LO ChPT or a significant difference to N f = 2 + 1 + 1 lattice results. For the d-wave our uncertainties are relatively large. However, thanks to the physical point ensemble we can directly compare to phenomenology and observe reasonable agreement. For the d-wave phase shift smaller scattering momenta would Table 6 Two-and three-pion fits using the Adler-zero form (z 2 = M 2 π , fixed). Since we only include s-wave interactions, we use two-pion levels in the A 1 irrep, and all irreps for three-pions. Recall that 1  Fig. 3 Constant(left) and linear(right) terms of K 3,df as a function of the s-wave scattering length. We also include the results of Ref. [69] be desirable in order to be able to shed light on a possible sign change at small k 2 -values.
For the three pion case, we observe reasonable agreement with other lattice calculations, phenomenology, and ChPT. By including two ensembles at heavier pion masses, we have gained insight on the chiral dependence of three-π + scattering quantities for the first time. We use an isotropic parametrisation of K 3,df , the real, singularity free, quasilocal, intermediate three particle scattering quantity. Here we find good agreement to LO ChPT for the constant term in K 3,df in an expansion about threshold, but an opposite sign compared to LO ChPT for the next-to-leading term. We have discussed possible explanations for this. On the other hand, qualitative agreement is found for both terms with the other available lattice calculation of these quantities.
This letter represents a step towards exploring and understanding the hadronic spectrum of QCD, and shows that three-particle quantities can be extracted with current techniques. In the very near future we expect more lattice calculations of three-body observables with increasing accuracy and describing systems with growing complexity -e.g. threeparticle resonances such as the ω.  [110][111][112], Lemon [113], QUDA [114][115][116], R [117], hadron [118] and paramvalf [119] have been used.

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: The energy levels of this work can be found in a GitHub repository [91].] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Extraction of the energy levels
In this section, we provide more details regarding the extraction of energy levels from the correlation functions of one, two and three charged pions. All the required quark contraction diagrams are shown in Fig. 4. For the observables in question we have determined the integrated autocorrelation times using the method put forward in [120] and found that we can treat our measurements as decorrelated. The statistical analysis is performed via bootstrap.

A.1 Thermal pollutions
Given the individual pion momenta p i , i = 1, 2, 3, we adopt the following convention to express the total momentum P = i p i and the relative momenta q j , j = 1, 2 The spectral composition of a three-pion correlation function (with periodic boundary conditions) reads The double sum is over all states m, n with the correct quantum numbers. The desired signal arises when m is the vacuum, and n the three-pion state, or vice versa. Usually one would expect that all other contributing to the spectral decomposition are exponentially suppressed compared to this ground state. Here this is not the case, because there are nonzero contributions to the spectral decomposition for finite T for instance when m is an intermediate two-pion state and n is a one-pion state. Such so-called thermal pollution states have a time dependence proportional to exp(− E t), with E = E 2π − E π , which can dominate the correlation function for large enough t when E < E 3π . There is an additional backward propagating part as well which goes as exp(− E · (T − t)). Together they either form a cosh (sum) or a sinh (difference). For three pions we only have time-even operators and therefore everything will have a cosh-shape. The amplitude of the cosh will be proportional to exp(−(E 2π + E π )T ), which vanishes for T → ∞. The size of the pollution will depend on the individual momenta of the involved pions through the energy E 2π and E π . The most significant pollution will be the one leading to the smallest E, which usually corresponds to the smallest involved momenta.
The thermal pollutions depend also on the frame and irrep. Let us illustrate this for a specific example: assume that n is a one-pion state |p 1 and m a two-pion state |p 2 , p 3 with free energies given by the dispersion relation. In this specific case only summands where p 1 | O |p 2 , p 3 = 0 contribute, i.e. the three-pion operator O must couple to the momenta p 1 , p 2 and p 3 .
The individual particle momenta that couple to a multiparticle operator can be inferred from group theory. Consider the frame P 2 = 0, then the three-pion operator will be in some irrep − , and the single pion always in the A − 1 . Therefore, the two-pion system needs to be in the opposite parity irrep + such that A − 1 ⊗ + = − . Note that only the irreps for P 2 = 0 have a parity index, that is, in moving frames parity is not a good quantum number. In this situation, the momenta of the two-pion system can only take the values that actually couple to the irrep of the three-particle operator.
The allowed contributions are generated from all permutations of the three-pion individual momenta. Using the measured pion rest mass M π and the free particle dispersion relation (assuming weak interactions between the two pions) we can thus estimate the relevant energies E π (p 1 ) and E 2π (p 2 , p 3 ). Using these together with the T -values we can now estimate for every ensemble, irrep and total momentum which thermal contribution is -up to unknown matrix elements -largest. Since we are able to remove only a single thermal state, this is the only way to single out the relevant parameters for the possible subtraction of these polluting states. Figure 5 shows the contributing thermal states for two example cases, left the A − 1 irrep in the P 2 = 0 frame, right the B 1 irrep in the P 2 = 2 frame. The different correlators shown correspond to different combinations of single and two pion momenta. For these cases the largest contribution is coming from (p 2 1 = 0, p 2 2 = 0, p 2 3 = 0) and (p 2 1 = 1, p 2 2 = 2, p 2 3 = 1), respectively. The other possible contributions are suppressed by two orders of magnitude or even exponentially.
To be precise, in order to find the dominating contribution for each irrep, ensemble and frame, we take the largest thermal contribution at t = 10, from which we can estimate E. To illustrate this procedure further, we will look at irrep = B 1 with P 2 = 2. The three-particle momenta that couple to the operator below our threshold are listed in Table 7. As the three particles are indistinguishable, we can partition them at will into a one-pion and two-pion state. The twoparticle momenta must again be a valid two-particle system, otherwise they cannot be an intermediate thermal state. Table  8  Each line corresponds to a particular combination of the individual particle momentum magnitudes of the one-pion (p 1 ) and the two-pion system (p 2 and p 3 ) Table 7 Possible three-pion individual momenta in the = B 1 irrep with total momentum P 2 = 2 Table 8 Possible two-pion individual momenta in the B 1 irrep for different values of P 2 of the two-pion subsystem Thus, again for the example of the B 1 irrep, we have to go through the following possibilities: • We take (0, 1, −1) for the one pion and (0, 0, 1) and (1, 0, 0) for the other two. The two-pion system has P 2 = 2, but the lowest contribution in that irrep has larger momenta. So this does not contribute. • The single pion has p 1 = (1, 1, 0) and the two-pion system gets p 2 = (1, 0, 0) and p 3 = (−1, 0, 0). The two-pion system therefore has total momentum P 2 = 0, but there is no contribution to B 1 in that moving frame. Therefore this example does not contribute to the thermal states. • A contribution is obtained using (1, 0, 0) for the one pion momentum, and (−1, 1, 0) and (1, 0, 0) for the two-pion system. In the latter, we have P 2 = 1, which corresponds to the first entry in Table 8 (albeit after an inconsequential global rotation). This contributes as a thermal state, incidentally it is the largest one as shown in Fig. 5b.
Of course, there are many more possibilities to check for. Using this method we determine the leading thermal state for every correlator matrix and can use this as input for thermal state treatments, detailed below.

A.2 General technicalities
Multi-particle correlators in general are contaminated with excited states at early times, and with thermal pollution at late time slices. Fitting too early will overestimate the energy, while fitting too late may underestimate it. In order to obtain a robust energy estimate, we use combinations of different methods to attenuate these issues. The order of application of these methods is illustrated with a flow chart in Fig. 6. The detour arrows indicate optional parts of the chain. We will explain the different methods in order. First the correlator matrices can optionally be treated with weight-shift-reweight [76] to suppress thermal states at the cost of larger statistical uncertainty. Then we independently use the original and treated correlator matrix and apply the GEVM, which yields the principal correlators. These principal correlators can be used to build ratios [17,31] or left as-is. All variants can optionally be fed into the Prony generalized Eigenvalue method (PGEVM) [90] with t 0 = 2 fixed to suppress excited states (The PGEVM with δ 0 fixed, see Ref. [90] for details, turned out to not be reliable).
The resulting treated correlators are evaluated by looking at the so called effective mass. The simplest definition of it which assumes a signal proportional to exp(−Et) only. There are generalizations that take back-propagation, shifting or weighting into account. Depending on the treatment of the correlator we choose the appropriate effective mass. We do not use all of the possible treatments in our analysis, but only the following five: no treatment (i.e. all optional parts are left out), only PGVM, only ratio, only weight and shift and finally the combination of weight and shift with PGEVM. In more detail this means: No treatment When no thermal states contribute (like in E irreps in the two pion channel), a simple cosh-like model is fitted: If thermal states are present in the given irrep, a two-state model with constrained second energy E 1 will be fitted to the data (for how E 1 and its error is determined, see Appendix A1) The constraint is implemented by augmenting the χ 2 function to be minimized by a term where E 1 is the fit parameter,Ē 1 is the determined central value for the thermal energy and δ E 1 the statistical uncertainty on E 1 . PGEVM This method works well when there are no significant thermal state contributions. We fit a simple exponential model at early times. Ratio We take the ratio of the principal correlator obtained from the GEVP (no weight-and-shift applied) and form ratios with the one-pion correlation function: The ratio R 3 is chosen as a double ratio such that in the numerator, thermal state contributions ∝ exp(− E t) are removed, since E ≈ E π . The resulting sinhlike correlator needs to be divided by another sinh-like expression, that's why we take the difference also in the denominator. Among different ratio expressions we have tested, this one works best in the sense that the plateau is longest. An exponential model is fitted to the ratios where the signal behaves like Note that for the ratios we do not include backwards propagating parts and thus do not extend fit ranges too far towards T /2.

Weight-shift
The correlator matrix has the leading thermal state removed [76] and, therefore, the principal correlators can be fitted with a cosh-like model which incorporates the weight-shift-reweight procedure.
Weight-shift and PGEVM In general the additional suppression of excited states by the application of the PGEVM works well after weight-shift has been applied beforehand. The resulting correlator is fitted with an exponential model. Fit ranges can be chosen early enough such that the neglect of backwards propagating parts is not significant. Figure 7 shows a comparison between no treatment, weight-shift and the ratio R 3 for a case with heavy thermal pollution. One can see how the effective mass of the plain correlator does not show any plateau due to the high degree of thermal pollution. The effective mass of the weighted correlator, however, exhibits a plateau between t 1 = 12 and t 2 = 14, but still shows a drop beyond. However, this three time slice plateau can only be identified when compared to the effective mass given by the ratio R 3 . This likely stems from the second leading thermal state as visible in Fig. 5a. The ratio however has a long plateau that is compatible with the weight-shift method a posteriori. In general we see that with the ratio method it is possible to fit energy levels with strong thermal pollution when other methods fail to produce a plateau. The statistical uncertainty from the energy determination with the ratio is also lower than with other methods in most cases.
In some cases the thermal states are so pronounced that no plateau can be identified, even after applying the PGEVM. In these cases the method is not used for that particular level.  These cases work much better with either the multi-state model, weight-shift-reweight or the combination of weightshift-reweight and the PGEVM. The ratio method seems to be the most robust one, it shows plateaus even when other methods fail to produce one. Also, the statistical uncertainty seems to be lower compared to the other methods in general.
For every principal correlator we attempt to extract the energy with all the five methods detailed above. If a plateau can be identified, we use the extracted energy level. All such determinations per principal correlator are combined with a correlated weighted average. In order to incorporate the systematic spread between the central values, we also compute a systematic error scaling factor as introduced in Ref. [32]: for energy level E we compute the scaling factor w Eq. (5), as mentioned in the main text.
To illustrate this method to incorporate the systematic error into the resampling distribution, we use two artificially generated data points with central values X 1 and X 2 and corresponding standard errors generated in four ways, where either the central values and/or errors are chosen to be the same or different. All combinations thus give four cases, which are shown in the quadrants of Fig. 8 (upper left: different mean, different errors; upper right: same mean, different errors; lower left: different mean, same errors; lower right: all the same). The central values with standard errors for X 1 and X 2 are shown as the first two pairs of points in each quadrant. The third pair shows the weighted average of the two estimates and the fourth pair the result after the rescaling. One can nicely see how the weighted average gravitates toward the data point with the smaller uncertainty (hence higher weight) and how the rescaling incorporates the spread between the central values. The method works well for both bootstrap and jackknife resampling.
In order to choose appropriate fit ranges for the different methods, we proceed iteratively, selecting fit ranges by eye guided by the p-value of the fit. Energy levels are included in the further analysis only if a plateau of at least five time slices length could be identified for the T = 96 lattices and of at least four time slices for the T = 64 lattice. Some energy levels show significant tension between the different fitting methods after this first iteration. In these cases, we re-evaluate the plateaus to arrive at our final choices.

Appendix B: Fitting the spectrum
Here, we aim to extend the discussion of the fitting procedure of the spectrum to the quantization condition. The summary of the frames, irreps and energies used in this work is shown in Table 9.  In both, the two and three-particle sector, we define the χ 2 as: where C is the covariance matrix of the energy levels, estimated from the bootstrap samples. Best fit parameters are obtained using the Levenberg-Marquardt algorithm. The range of validity of the quantization conditions is limited by the first inelastic threshold. This is E * = 4M π (5M π ) for the two-particle (three-particle) quantization condition. We generally include levels up to that threshold, however, for the physical point ensemble (cA2.09.48), we have included levels higher up in energy. Since the 2π → 4π , and 3π → 5π couplings are very small, we expect this to be a valid approximation. In fact, phenomenological studies set the first relevant inelasticity to be the ρππ channel (E * ∼ 8M π for physical kinematics) [92,95,96].
As mentioned in the main text, we show here additional two-pion phase shift plots: Fig. 9 for cA2.30.48, and Fig. 10 for cA2.60.32. In the case of the s-wave phase shift, we also compare to LO ChPT. As can be seen, the ChPT prediction describes less accurately the data at heavier pion massescompare to Fig. 2.
B.2 Additional discussion on three-pion fits First, we perform a global fit to two-and three-particle levels that includes only a constant term in K 3,df . This is shown in Table 10. As can be seen, the quality of the fit is significantly worse for the heavier ensembles than in the linear fits of Table 6 in the main text. For the ensemble at the physical point (cA2.09.48), the value of χ 2 is basically the same, but in both cases K 3,df is compatible with zero. We thus conclude that the linear model of K 3,df in Eq. 6 in the main text is more appropriate for this system.
Next, the full covariance matrices of the fits in Table  6 in the main text are provided. We use the form C = D R D, with D being a diagonal matrix with the standard errors of the parameters. We ordered the entries as: 1/B 0 , B 1 , M 2 π K iso,0 df,3 , M 2 π K iso,1 df,3 . For s-wave we use a model that incorporates the Adler-zero, whereas for d-wave we fit to a constant in the region for which we have data. Two points have been omitted in the plot due to the very large errorbars Table 10 Two-and three-pion fits using the Adler-zero form (z 2 = M 2 π , fixed). Here we assume that K 3,df is given by a constant: K 3,df = K iso,0 We observe a large correlation within the two and threeparticle sectors separately -the pairs 1/B 0 , B 1 , and M 2 π K iso,0 df,3 , M 2 π K iso,1 df,3 are highly correlated. In contrast, the correlation between the two-and three-particle sectors is milder.
(a) (b) Fig. 11 The center-of-mass spectrum for two and three pions on the physical point ensemble (cA2.09.48). The red data points are the energy levels determined from the correlator. The black lines denote the prediction from the quantization condition. For the two-pion A 1 levels, and all three-pion levels, we use the fit in Table 6 in the main text. For the non-A 1 two-pion levels, which are dominated by d-wave interactions, we use the fit in Table 4 in the main text. The short dashed gray lines denote the noninteracting energy levels. We also include the relevant inelastic thresholds as long dotted gray lines B.3 Two-and three-pion spectrum We conclude the discussion by comparing the spectrum from the lattice to the one predicted by the quantization conditions using the best fits. This is shown in Fig. 11 for the ensemble at the physical point.