Interactions of two and three mesons including higher partial waves from lattice QCD

We study two- and three-meson systems composed either of pions or kaons at maximal isospin using Monte Carlo simulations of lattice QCD. Utilizing the stochastic LapH method, we are able to determine hundreds of two- and three-particle energy levels, in nine different momentum frames, with high precision. We fit these levels using the relativistic finite-volume formalism based on a generic effective field theory in order to determine the parameters of the two- and three-particle K-matrices. We find that the statistical precision of our spectra is sufficient to probe not only the dominant $s$-wave interactions, but also those in $d$ waves. In particular, we determine for the first time a term in the three-particle K-matrix that contains two-particle $d$ waves. We use three $N_f=2+1$ CLS ensembles with pion masses of $200$, $280$, and $340\;$MeV. This allows us to study the chiral dependence of the scattering observables, and compare to the expectations of chiral perturbation theory.


Computation of finite-volume energies
The computational methods follow closely those employed in Ref. [2], and we review only the high-level features in this section.

Interpolating operators
In order to extract the finite-volume spectrum reliably from LQCD, it is imperative to include interpolating operators with good overlap onto the states of interest. Operators annihilating a two-pion and three-pion state are given by a sum over products of single-pion annihilation operators, each projected to definite momentum p i , ππ (P,Λ) = c (P,Λ) p 1 ,p 2 π p 1 π p 2 , (2.1) πππ (P,Λ) = c (P,Λ) p 1 ,p 2 ,p 3 π p 1 π p 2 π p 3 . (2. 2) The corresponding creation operators are obtained by taking the Hermitian conjugate of the annihilation operators. The momentum combinations encoded in the Clebsch-Gordan coefficients c (P,Λ) are chosen such that the resulting operators transform according to the irreducible representation (irrep) Λ of the little group of the total momentum P = i p i . In this basis, the finite-volume Hamiltonian assumes a block-diagonal form, reflecting the symmetry of the cubic spatial volume or its boosted deformations, thus greatly simplifying the extraction of the spectrum. In addition to the total momenta d 2 = L 2 /(2π) 2 P 2 ≤ 4 commonly used in previous work, total momenta up to d 2 ≤ 9 are included in the present work (excluding momentum of type d = [122]). In particular, the trivial irrep with total momentum of type d = [003] proves to be useful in the three-particle sector, since it provides an additional data point in the threshold region purely on kinematical grounds. The group-theoretical projection of meson-meson operators proceeds along the lines of Refs. [42][43][44]. Our choices for representation matrices of the elements of the little group of the newly included momentum classes d = [012] and d = [112] are given in Appendix A. Three-pion operators are then obtained by iteratively applying the two-meson coupling procedure [2,45], and all interpolators used in this work are tabulated in Appendix B.
Two-kaon and three-kaon annihilation operators are given by replacing the constituent pion annihilation operators with kaon annihilation operators in Eqs. (2.1) and (2.2), with the Clebsch-Gordan coefficients unchanged.

Correlation functions
Calculating correlation functions of operators of the form of Eqs. (2.1) and (2.2) requires a method to handle quark propagation from all spatial sites on the source time slice to all spatial sites on the sink time slice in order to be able to perform the momentum projection for each constituent hadron individually. Such a method is furnished by distillation [12], which treats quark propagation in a low-dimensional subspace that encodes the information relevant for hadronic physics.
The distillation subspace is spanned by the N ev lowest-lying eigenvectors of the threedimensional gauge-covariant Laplacian on each time slice of the lattice. The projection of quark fields into this subspace is equivalent to a smearing of the quark fields with an approximately Gaussian profile having a characteristic smearing radius. In order to keep the smearing width constant as the physical spatial volume of the lattice is increased, the number of retained eigenvectors needs to be scaled linearly with the volume. The concomitant increase in computational cost can be ameliorated by employing stochastic estimators in the distillation subspace [11]. In this so-called stochastic LapH method, a valence quark line is estimated stochastically, Q aα,bβ (x, y) = lim Nr→∞ 1 N r r,d φ (r,d) aα (x)ρ bβ (y) (r,d) * , (2.3) using N r independent sets of diluted [46] stochastic combinations of LapH eigenvectors as sources ρ(y) and corresponding solutions of the Dirac equation Dφ = ρ, where r and d = 1, . . . , N dil denote the noise and dilution indices, respectively, a, b are color indices, and α, β are Dirac spin indices. This method to treat quark propagation affords great flexibility as the computationallyexpensive solutions of the Dirac equation can be re-used across several spectroscopy projects. In a subsequent step, quark sources and solutions are combined into meson functions with support only on the source or sink time slice [11]. The meson functions are rank-two tensors with two open dilution indices and a compound label identifying their spin and spatial structure, meson momentum, and pair of quark noises.
The final step in computing correlation functions consists of performing tensor contractions over dilution indices of sets of meson functions according to the Wick contractions of quark fields [11], and forming the linear combinations of individual momentum assignments governed by the group-theoretical projections discussed in Section 2.1. This step becomes more computationally expensive as systems of an increasing number of valence quarks are considered, and with a naive implementation the associated cost completely dominates that of the whole calculation already for three-meson systems. A significant speedup can be achieved by systematically eliminating all redundant computation through common subexpression elimination [47]. The application of this idea in the present LQCD context was described in Ref. [2]. The publicly available implementation 2 is not restricted to systems of mesons and was used recently to speed up a calculation of two-baryon systems by nearly three orders of magnitude [48]. In the present work, these improvements enable the computation of up to 20,000 distinct correlation functions per gauge configuration and source time, encompassing the evaluation of up to one billion individual diagrams (as defined in Ref. [2]).

Ensemble details
Calculations in this study are performed on three ensembles at a fixed lattice spacing generated through the CLS effort [49]. The simulations use N f = 2+1 nonperturbatively O(a)improved Wilson fermions and the Lüscher-Weisz gauge action with tree-level coefficients. They are performed along a chiral trajectory keeping the trace of the quark mass matrix  N203 48 3  Overview of the lattice geometry, approximate pseudoscalar masses, number of gauge configurations N cfg , fixed source-time positions t src , number of Laplacian eigenvectors N ev , employed dilution scheme, and number of independent quark noises N r used to estimate light and strange quark propagation for ensembles used in this work. The dilution scheme notation is explained in Ref. [11]. All ensembles share the same lattice spacing a ≈ 0.064 fm.  Table 2. Pion and kaon decay constants determined in Ref. [51].
fixed, so a heavier-than-physical light quark mass implies a lighter-than-physical strange quark mass. An overview of the three ensembles used in this work as well as the computational setup is given in Table 1. The lattice spacing on these ensembles was determined to be a = 0.06426 (76) fm using the linear combination 2 3 (F K + 1 2 F π ) of decay constants to set the scale [50]. In addition, the decay constants used in this work were computed in Ref. [51] and are reproduced in Table 2 for convenience. Open temporal boundary conditions were imposed when generating the ensembles to avoid topological charge freezing at fine lattice spacings [52]. Consequently, translation invariance in the temporal direction is broken, and the position of source operators is fixed to the values given in Table 1 rather than being randomized on each gauge configuration. The effect of the boundary conditions on spectral quantities is expected to decay exponentially with the distance from the boundary, with the decay constant governed by the lightest state with vacuum quantum numbers, which is expected to be a two-pion state for the quark masses used in this work. They are thus expected to be most pronounced on the ensemble with the lightest pion mass, the D200 ensemble in the set used in this work. In a previous study on the same ensemble, t src = 32 was found to have negligible temporal boundary effects [53]. The sources placed in the bulk of the lattice for this study are thus expected to be sufficiently far away from the boundary. Additionally, the correlators are always constructed such that the sink times are toward the center of the lattice with respect to the source position. Therefore, the additional source at t src = 92 for the D200 ensemble required backward-time correlators. The operators used for the backward-time correlators are related to the operators in the forward-time correlators by a parity and charge-conjugation transformation.
Except for a twofold increase in statistics on the D200 ensemble and the use of an improved dilution scheme on the N200 ensemble, solutions of the Dirac equation for the  Table 3. Single-hadron energies at rest, determined from single-exponential fits. The range of time separations included in the fit is given by t fit . The high χ 2 red for N203 is discussed in the text.
light quark are re-used from a previous spectroscopic calculation supporting the lattice determination of the hadronic vacuum polarization contribution to the anomalous magnetic moment of the muon [54].

Analysis of correlation functions
Matrices of correlation functions are computed for a wide range of total momenta and irreps. In a first step, the data is averaged over equivalent momenta, irrep rows, and source times on each gauge-field configuration. The subsequent analysis is performed using jackknife resampling and for each total momentum-squared/irrep separately. We provide the jackknife samples of the resulting spectrum in HDF5 format as ancillary files with the arXiv submission.

Single-hadron energies
Pion and kaon correlators are computed in each momentum frame, and then used to extract the single-hadron energies for each total momentum-squared. We employ a singleexponential correlated-χ 2 fit to these correlation functions. The results for the kaon and pion energies at rest are given in Table 3. The high χ 2 red ≡ χ 2 /dof for the single-hadron fits on N203 could be reduced by rebinning the data (as was done on N200 and D200). However, due to the smaller number of configurations for N203, rebinning of the data leads to an unstable covariance matrix when fitting the full set of multi-hadron energies to the quantization condition (an issue discussed further below).
When fitting the multi-hadron spectrum to the quantization condition, the continuum dispersion relation is assumed to be valid for the pion and kaon, and therefore only the energies for the single hadrons at rest are actually needed. We can test this assumption by studying the dispersion relation of our single-hadron states. The results, shown in Figure 1, show no appreciable disagreement with the continuum dispersion, suggesting that discretization effects are small.

Multi-hadron spectrum extraction
In order to extract not only ground states but also a tower of excited states in each irrep, the eigenvectors v n determined from the generalized eigenvalue problem (GEVP) [55,56]  are used to form the correlation functions of rotated operators with 'optimal' overlap onto the nth state in the spectrum,Ĉ where the parentheses denote an inner product over the operator indices. The diagonalization is performed for a single (t 0 , t d ) only, keeping t 0 t d /2 [57], but we checked for a range of sensible values that the resulting spectrum is independent of that choice. Two-and three-pion finite-volume energies are extracted from single-exponential correlated-χ 2 fits to the ratios , (2.6) where the product in the denominator is over two or three single-pion correlation functions, respectively, with momenta chosen to match the closest noninteracting energy to the nth state. The ratio Eq. (2.6) gives access to the energy splitting between the interacting and noninteracting state at sufficiently large time separation, and the lattice energy is reconstructed using the single-pion mass and dispersion relation. Fits are performed to the ratio data in the time separation range [t min − t max ] with t max = 40 fixed on all ensembles (with the exception of the pion spectra on N200, for which t max = 36 led to a more consistent analysis). The lower bound of the fit window is chosen such that the residual excited-state contamination is subdominant compared to the statistical uncertainty. The analysis of the two-kaon and three-kaon data proceeds analogously with the product in the denominator of Eq. (2.6) replaced with single-kaon correlation functions.
Fitting to the ratio of correlators, rather than directly to the rotated correlator, generally allows for an earlier choice of t min and thus a more precise energy extraction. The  . The extracted energy shift in the lab frame as a function of the smallest time separation included in the fit, t min , for the three-kaon A 2 (9) ground state energy on each ensemble. The fits are to the correlator ratio, and a single-exponential is used for the model. The central value and error of the energy shift for the chosen fit is indicated with the horizontal black solid and dashed lines, respectively. dependence of these fits on the choice of t min is shown in Figure 2 for the three-kaon A 2 (9) ground-state energy on each ensemble. The earlier plateau in the correlator ratio can be understood as coming from a partial cancellation of correlations and excited states between the numerator and denominator. However, one disadvantage is the loss of a monotonic decrease in the effective energy of the correlator ratio. Therefore, we verify consistency between the results from the ratio fit and direct fits to the rotated correlators.
The resulting energy splitting a∆E from the fits to the correlator ratio are converted to absolute energies in the center-of-momentum frame aE * by using the extracted singlehadron energies at rest with the continuum dispersion relation. All of the extracted threekaon and three-pion energies on N200 are shown in Figure 3, along with the energies resulting from a global fit using the two-and three-particle quantization condition, to be described below.
The global fits account for the correlations between all levels, i.e., those with both two and three particles of a given type. Since there are many levels in a given fit (up to 77) one might be concerned about the reliability of our determination of the covariance matrix. To assess this, we compared fits to the spectrum with both bootstrap and jackknife resampling and found them to be consistent. On N200 and D200, we also found consistency for the central values of our fits when rebinning the data by 2 and 3, with the errors increasing due to autocorrelations. On N203, we found that rebinning the data leads to inconsistent fits, which can be explained by the much smaller number of configurations as compared to N200 and D200. In fact, using a smaller subset of energies on N203, such that the covariance matrix is stable with the rebinned data, still gives consistent results with the same subset of energies without rebinning, suggesting the rebinning is not necessary in the final analysis on N203. With our final choices of rebinning (3 for both N200 and D200, and none for N203), the number of bins is roughly consistent across all ensembles. This is approximately an order of magnitude larger than the number of energies going into our largest fits. All of our consistency checks suggest this is sufficient for a reliable estimate of the covariance matrix.  In this section we summarize the formalism that we use to fit the spectrum, describe our overall fitting strategy, and present the parametrizations of the K-matrices that are used to describe the two-and three-particle interactions. In addition we collect various results from ChPT that are needed to analyze the mass dependence of our results.

Quantization conditions
To relate the finite-volume spectrum to infinite-volume scattering parameters, we need both two-and three-particle quantization conditions. The former is standard-the original Lüscher quantization condition [58,59] generalized to moving frames in a relativistic formalism [60,61]. For the latter, we use the results of the relativistic field theory (RFT) approach. This formalism, developed in Ref. [14], holds for three identical, spinless particles within the kinematic range for which only a single three-particle channel can go on shell. It applies for relativistic particles, and leads to Lorentz-invariant scattering amplitudes as long as one uses the relativistic form of the kinematic functions (a point discussed extensively in Ref. [22]). Practical implementation of the quantization condition requires truncating the angular momentum of the interactions between pairs of particles. We follow Refs. [19,20] and include both s and d waves (p waves being forbidden for identical particles). Where we break new ground is the inclusion of s and d waves in moving frames; previous work in the RFT approach with moving frames was restricted to s waves [3].
We present here only a summary of the formalism and its implementation, since most of the details are the same as in Refs. [3,19,20]. The formalism is derived for a generic continuum effective field theory restricted to a cubic spatial box of length L. Thus the allowed total momenta are drawn from the finite volume set: P = (2π/L)d, where d ∈ Z 3 . For a given choice of P and L, the two-particle spectrum is given by the energies where E * 2 = E 2 2 − P 2 is the center-of-mass frame (CMF) energy, F is a kinematical function to be discussed shortly, and K 2 is the two-particle K-matrix given in Eq. (3.5) below. The three-particle spectrum is given by the energies E that solve where E * = E 2 − P 2 , F 3 will be discussed shortly, and K df,3 is a three-particle Kmatrix discussed around Eq. (3.13). The quantization conditions hold up to exponentially suppressed corrections, which should be small for our values of M π L and M K L. Although the two quantization conditions have a similar form, they differ substantially in the details. The determinant in Eq. (3.1) runs over indices { , m}, which denote the angular momentum of the two particles in their CMF, while that in Eq. (3.2) runs over {k m}, where k represents the three-momentum of one of the particles, k, which is drawn from the finite-volume set, while { , m} are the relative angular momenta of the other two particles in their CMF. The formalism has a built-in cutoff, such that only a finite number -10 -of values of k contribute, but the sum over { , m} in both quantization conditions must be truncated by hand. As noted above, here we consider max = 2.
Another difference between Eqs. (3.1) and (3.2) concerns the first entry in the determinants. For two particles, the matrix F is a purely kinematic function (a "Lüscher zeta function"), which encodes the effect of working in finite volume, and, for a general moving frame, couples s and d waves. We use the form given in Appendix A of Ref. [19], extended to moving frames. For three particles, F 3 contains not only kinematical functions (F , together with an additional function, G), but also the two-particle K-matrix. Again, the details are given in Ref. [19], except here extended to moving frames. The extension of the implementation of the s-wave-only formalism to moving frames is described in the supplementary material to Ref. [3]. The generalization to include d waves is straightforward.
The final difference between the two quantization conditions concerns the second term in the determinants. For two particles, the K 2 that appears is a version of the two-particle K-matrix, differing from the standard choice by some additional cutoff terms (discussed below). It is an infinite-volume quantity that is algebraically related to the physical twoparticle scattering amplitude M 2 . For the three-particle quantization condition, what appears instead is a three-particle K-matrix, K df,3 , which is an infinite-volume but cutoffdependent amplitude, related to the physical three-particle amplitude M 3 through integral equations [15]. What matters here is that K df,3 is a real, analytic function of Lorentz invariants, and thus can be simply parametrized, as discussed below.
Both quantization conditions can be block-diagonalized by projecting onto irreps of the appropriate symmetry group for the given choice of P . This is the little group LG(P ), composed of elements of the cubic group O h that leave P invariant. We implement this projection using the formalism developed in Ref. [19], extended in Ref. [3], and here generalized to allow for the inclusion of two-particle d waves in moving frames. For given choices of K 2 and K df,3 , the solutions to the quantization condition are obtained by tracking the small eigenvalues of the matrices in the determinants, irrep by irrep, and numerically determining where they cross zero. This procedure has been independently implemented in two Python codes, and in a (much slower) Mathematica code, so that all numerical results presented below have been cross-checked.
As noted in previous implementations of the RFT approach, there are choices of the functions K 2 and K df,3 , as well as of M π L and M K L, for which there are unphysical solutions to the three-particle quantization condition [17,19]. Their source is not fully understood, but is likely due to a combination of numerically enhanced exponentially suppressed errors (which are not controlled in the derivation) and the use of overly restrictive choices for K 2 and K df,3 . We have checked that, for the fits presented in the next section, there are no unphysical solutions in the energy range considered.
Another source of error that we do not systematically control is due to discretization effects. Strictly speaking, the quantization conditions hold only in the continuum limit, since we are assuming relativistic dispersion relations for intermediate particles, and restrict the two-and three-particle interactions using Lorentz invariance. The potential size of discretization errors is unclear. On the one hand, the dispersion relations shown in Figure 1 show no indication of large discretization effects, consistent with the findings of an investigation of resonant I = 1 pion scattering on ensembles with the same action [53]. On the other hand, a recent study of two-baryon scattering amplitudes from LQCD serves as a cautionary tale about potentially large discretization effects in these types of calculations [62]. What is clear is that more work, involving multiple lattice spacings, is needed to demonstrate control over discretization effects.

Fitting strategy
We now describe the overall strategy used for fitting. An important consideration is that there is not a one-to-one relationship between the energy levels and the K-matrices-this holds only for the two-particle spectrum in the approximation of a single partial wave. Thus one must parametrize the K-matrices and perform a global fit to the whole set of levels. We perform such fits to the two-particle spectrum alone, and to the combined two-and three-particle spectra. We note that, on a given ensemble, these two spectra are correlated, as are the individual levels, and we perform a fully correlated fit. More precisely, our fits to the spectrum are carried out by minimizing the following χ 2 function with respect to the parameters in the K-matrices, {p n }, where E j are the measured energy levels, with covariance matrix C, and E QC i ({p n }) are the energy levels predicted by the quantization conditions for a given choice of {p n }. The model functions yielding the E QC i ({p n }) depend on the data through the particle masses and box length L, so the covariance of the residuals [63]. However, we use C ij = cov(E i , E j ) instead of cov(r i , r j ) since it makes little difference to the final fit parameters but significantly simplifies the analysis. We estimate C by means of jackknife samples, ignoring the sample to sample variations in M π and M K , which are very small effects compared to the variations in the E i . We do not apply the correction discussed in Ref. [64] since any resulting change to the final fit parameters are expected to be insignificant.
In order to estimate the uncertainties of the best fit parameters, we use the derivative method. This uses the matrix of covariances between the parameters p n and p m , given by where the partial derivatives are evaluated numerically at the minimum of χ 2 . The main advantage of this approach is that one avoids refitting the spectrum in each sample, which is computationally very expensive for the three-particle case. To ensure the validity of this procedure, we have checked for two-particle fits that it yields almost identical uncertainty estimates as obtained with the standard jackknife method. We have also checked that using bootstrap instead of jackknife samples does not lead to significant changes in the results. One might have thought that, instead of a global fit, a better procedure would be to fix K 2 from fits to the two-particle spectrum, and then use the result as input into the -12 -three-particle spectrum, so as to determine K df,3 . Such a procedure would, however, fail to make use of the constraints on two-particle interactions arising from the interactions of pairs in the three-particle system. Indeed, the three-particle spectrum provides significant additional constraints on two-particle interactions, since there are three interacting pairs.
As a final comment on fitting, we consider the appropriate range of E * to use. Since we are considering isosymmetric QCD, G parity is an exact symmetry, and there are no transitions between sectors with even and odd numbers of pions. Thus, for two pions, the range of validity of the quantization condition is 0 < E * < 4M π , while for three pions it is M π < E * < 5M π . For kaons, by contrast, there is no constraint from G parity, and the process K + K + → K + K + π 0 is allowed, for example if the kaons are in a relative d wave. Also allowed is K + K + → K + K 0 π + . Because of this, the upper bounds on the validity of the quantization conditions for two and three kaons are, strictly speaking, 2M K + M π and 3M K + M π , respectively. We expect, however, the coupling to the channels with an additional pion to be very small. This is because these transitions are induced by the chiral anomaly, and thus, in ChPT, by the Wess-Zumino-Witten (WZW) term. However, the WZW term itself does not lead to a K + K + → K + K 0 π + transition, due to its antisymmetry; the closest transition that is induced is K + K 0 → K + K 0 π 0 . To obtain the desired transition requires an additional loop (which attaches a K 0 K + → K 0 K + vertex). This implies that the process is at least of next-to-next-to-leading order (NNLO) in ChPT, since the WZW term is itself of next-to-leading order (NLO). Thus we expect the coupling to an additional pion to turn on only far above threshold, where the p 4 suppression is less significant. On the other hand, we do not include mixed-flavor operators in our calculations to capture these extra levels above threshold. Although this could result in unreliable energy extractions in this region, we expect that the overlap of our operators with these mixed-flavor states is small enough to still obtain an accurate determination of the desired energies. Because of this we have included in the fits to 2K + and 3K + some levels that lie slightly above the strict threshold. Details are given in the next section.

Parametrizations of K-matrices
We now summarize the parametrizations that we use for K 2 and K df,3 . The former is given by where the subscript on E * 2 emphasizes that this is the CMF energy of a two-particle system, and Here q is the magnitude of the three-momentum of each of the two particles in their CMF, and H(q 2 ) is a smooth cutoff function that equals unity for physical scattering (q 2 ≥ 0). We follow all previous work implementing the RFT formalism and use the form of H(q 2 ) from Ref. [14]. We stress that the term proportional to [1 − H(q 2 )] vanishes identically -13 -for physical scattering. It is not necessary for the two-particle quantization condition, although it can be included (H dependence for subthreshold two-particle finite-volume states is cancelled by the corresponding dependence of F ), but it is essential for the threeparticle case, where subthreshold momenta are included down to q 2 = −M 2 . We also note that the additional scheme dependence in K 2 introduced in Ref. [20] is not needed here as there are no resonances.
To complete the parametrization of the two-particle phase shifts we need to specify the forms we use for cot δ . For s-wave scattering we consider two choices. The first is motivated by the ChPT expressions for the scattering of identical pseudo-Goldstone bosons (e.g., π + π + and K + K + ), which includes the Adler zero below threshold: This contains the dimensionless parameters z 2 , B 0 , B 1 , and B 2 . We stress that, here and below, the full expressions contain an infinite set of higher order terms in q 2 that we are setting to zero. We use this expression for both pions and kaons, with M = M π and M K , respectively. At lowest order in ChPT, z 2 = 1, and we perform both fits with z 2 fixed to this value, and others allowing it to vary. We also do fits with and without the parameter The second form is the effective-range expansion (ERE) given in terms of the scattering length a 0 , effective range r 0 , and quadratic parameter P 0 . We use the sign convention in which the scattering length is positive for repulsive interactions. In principle, the radius of convergence of the ERE is given by the position of the Adler zero, although this form has been used beyond this radius in many previous LQCDbased analyses of above-threshold scattering. The parameters in Eq. (3.8) are related to those of the ERE through We use the combination M 2 a 0 r 0 rather than M r 0 alone, since the former does not diverge in the chiral limit, whereas the latter can [see Eq. (3.16) below]. For d-wave scattering we use a simpler, one-parameter form, as we find that this provides an adequate description of our data. The overall factor of E * 2 = √ s is adopted from standard continuum analyses (see, e.g., Refs. [65,66]) and implies that higher-order terms in the ERE for δ 2 are present, although with fixed coefficients. The factor of −1 is chosen to avoid unphysical poles in the subthreshold scattering amplitude, which is given by Such poles arise for |q| ≈ M from the D 0 term alone if the interactions are repulsive (as they are here). The use of this ad hoc factor is adapted from that used in continuum analyses of s-wave scattering [65]. The d-wave scattering length is then given by M 5 a 2 = −1/(D 0 − 1), using the conventional definition in which a 1/5 2 has dimensions of length, and the same sign convention as for a 0 . For the three-particle K-matrix we use the threshold expansion worked out in Ref. [19]. This expands K df,3 (p 1 , p 2 , p 3 ; p 1 , p 2 , p 3 ) in powers of ∆ = (E * 2 − 9M 2 )/(9M 2 ) and related quantities, where p i (p i ) are the initial (final) momenta for on-shell three-to-three scattering. We use the expansion through quadratic order, where and K iso,0 df,3 , K iso,1 df,3 , K iso,2 df,3 , K A , and K B are real constants. The first three terms depend only on the overall CMF energy, but not otherwise on the momenta of the particles, and are referred to as "isotropic" contributions. The final two terms do depend on momenta through the dimensionless quantities ∆ A and ∆ B , which are defined in Ref. [19], and are of quadratic order in ∆. Eq. (3.13) is the most general form consistent with Lorentz symmetry, time-reversal, parity, and particle exchange symmetry. To use this result in the quantization condition, Eq. (3.2), one must decompose K df,3 into the {k m} basis discussed above. How to do so for the rest frame (P = 0) is explained in Ref. [19], and the generalization to moving frames is straightforward. We also note that, when decomposed into irreps of the little group of the various frames, the K iso df,3 and K A terms contribute only to trivial irreps, while the K B term can also contribute to nontrivial irreps, and does so to all the nontrivial irreps included in our fits. It is for this reason that K B can be determined more easily than K A . In fact, in our minimal fits we use only the three parameters K iso,0 df,3 , K iso,1 df,3 , and K B . In the two-particle context, it is standard to consider the total angular momentum J of the system, and we have done that above when distinguishing s-and d-wave choices for K 2 , which have J = 0 and 2, respectively. The finite-volume irreps do not uniquely pick out values of J, but do restrict them. For example, the two-particle E + g irrep in the rest frame contains J = 2, 4, . . . but not J = 0, while trivial irreps in all frames contain J = 0 as well as higher values. The relation between irreps and J or helicity components is given, for example, in Table II of Ref. [67].
In infinite volume, one can also classify the three-particle interactions by their total angular momentum. Although this is not discussed in Ref. [19], it is straightforward to see from the explicit forms of ∆ A and ∆ B that the K A term contains only J = 0, as do the isotropic terms, while ∆ B contains J = 0, 1 and 2 components. The presence of these higher angular momenta provides an alternative way of understanding why only the K B term contributes to nontrivial reps, since these require J > 0. We note that a total J = 2, say, can arise both from an = 2 contribution to the pair sub-interaction combined with a relative s-wave between the pair and the remaining particle, and from an = 0 contribution combined with a pair-spectator relative d-wave, as well as from other options.

Results from chiral perturbation theory
We collect in this subsection the results from ChPT that we need when fitting the dependence of the quantities derived from fits versus M 2 π . We use both SU(2) and SU (3) ChPT. For SU (2) ChPT to be valid, we need both M 2 π /(4πF π ) 2 1 and M 2 π /M 2 K 1. The former is indeed small, ranging up to 0.075 on our ensembles. The latter is significantly larger, with the values {0.18, 0.38, 0.61} on the D200, N200, and N203 ensembles, respectively. Thus SU(2) ChPT is of borderline applicability for the N203 ensemble. The rapid rise in this ratio is due to the path along which we take the chiral limit, namely with m av = (2m q +m s )/3 fixed (with m q the common up-and down-quark mass and m s the strange-quark mass). At leading order (LO) in ChPT, this implies that In the following, we use the abbreviations where F π and F K are the physical pion and kaon decays constants, in the F π 92 MeV convention. We work mostly at NLO, and do not indicate the higher-order terms that are dropped.
We begin with the expression for the pion scattering length in SU (2) ChPT [68] (choosing the form commonly used in the lattice community and given in Ref. [69]) Here ππ is a low-energy coefficient (LEC) evaluated at the scale √ 2F π . A subtlety here concerns the chiral trajectory that we follow. In SU(2) ChPT, ππ depends implicitly on m s , and so has, in principle, different values for our three ensembles, none of which equal the standard value, which has m s = m phys s . Since ππ is a smooth function of m s , the difference between the values on our ensembles and the standard value is proportional to m s − m phys s = 2(m phys q − m q ), and thus proportional to M 2 π . As the ππ term is itself of NLO, the part proportional to M 2 π is effectively of NNLO (or higher order), and can be consistently dropped.
The expression for the effective range in SU (2) ChPT is [70] where C 3 is an LEC, and is evaluated at the scale √ 2F π . Similarly to ππ , C 3 can be treated as a constant in an NLO expression.
One can also obtain the NLO expression for the position of the Adler zero in pion scattering. Rewriting the result of Ref. [71], we have where z is an LEC evaluated at √ 2F π , which is related to the standard SU(2) ChPT LECs (see, e.g., Ref. [69]) by The d-wave I = 2 pion amplitude vanishes at LO in ChPT. Although a NNLO expression exists in the literature, our results are insufficient to attempt a fit to this form. Furthermore, there is a subtlety related to a change in sign of the phase shift at physical pion masses just above threshold, an issue we discuss below when we analyze our results.
Expressions in SU (3) ChPT for the pion and kaon scattering lengths can be obtained from Refs. [72,73]. The NLO terms depend not only on the pion and kaon masses, but also on that of the η meson. Within these terms, it is consistent to rewrite the η mass using the LO expression 3M 2 η = 4M 2 K − M 2 π , and also to treat F π and F K as interchangeable. In this way, we can write the results as functions of x π and x K alone, finding (3.20) These two expressions contain only a single LEC, L ππ , which is here evaluated at the scale 4πF π . We stress that, in SU (3) ChPT, there is no subtlety due to our choice of chiral trajectory, since all dependence on m s is explicit. In other words, the trajectory is encoded into the values of x π and x K . Expressions in SU (3) ChPT for the kaon effective range and Adler zero position could, in principle, be extracted from the results given in Ref. [8] for the K − K − scattering amplitude. We have not done so, however, as our simulation results for M 2 K a KK 0 r KK 0 presented below lie close to unity, and thus very far from the LO chiral prediction of 3. This indicates very large NLO corrections, and a breakdown in convergence.
An alternative would be to use SU (2) ChPT, treating the kaon as a heavy source for pions, following Ref. [74]. However, the analysis for kaon scattering has not been carried out (and would require methods similar to that used to study N N scattering in EFT). Thus we use simple analytic parametrizations of Finally, we turn to the chiral predictions for the three-particle K-matrix. At tree level (LO) K df,3 is purely isotropic, with only K iso,0 df,3 and K iso,1 df,3 being nonzero [3] These expressions hold both for three-π + and three-K + systems, using the corresponding masses and scattering lengths. The other constants in Eq. (3.13) (K iso,2 df,3 , K A , and K B ) can appear first at NLO, and thus are suppressed by at least one additional power of 1/F 2 ∝ M a 0 .

Fitting the spectrum
In this section, we present the results of fits to the two-and three-particle spectra, describing in turn the results for pions and kaons. The parametrizations explained in Section 3.3 are used. We comment on the main features of the fit, such as goodness of fit and which parameters are needed, but leave the interpretation of the results for the fit parameters themselves to the following section. An example of these fits on the N200 ensemble was already shown in Figure 3. We also note that the precise set of levels used in the fits of this section is given in Appendix C.

Fits to the spectrum of two and three pions
We start with fits to the energy levels in the 2π + and 3π + sectors. For the two-pion s-wave interaction, we use the Adler-zero form, given in Eq. (3.8), rather than the ERE form. The former has been found previously to provide a better description for light pions [6,19], and we confirm this result below. When d waves are included, we use the expression in Eq. (3.11). For K df,3 , we consider only three of the terms in Eq. (3.13): the leading two isotropic ones, K iso,0 df,3 , K iso,1 df,3 (referred to below as s-wave terms), and K B (referred to as a d-wave term). We find that this choice provides a good description of the levels, whereas the other parameters in K df,3 are poorly determined if included.
We present results for the following set of representative fits: 1. A fit solely to two-particle energies, including both s-and d-wave interactions. This fit shows the information that can be extracted from the two-particle spectrum alone, and thus is a useful point of comparison for three-particle fits.

2.
A combined fit to two-and three-pion energies including only levels in the trivial irreps in each frame, and with only s-wave contributions for both two-and threeparticle interactions. Furthermore, the position of the Adler zero is fixed to its LO ChPT value (z 2 = 1). This is the type of fit used in previous work [6,19], albeit to fewer frames than available here.

3.
A combined fit to all two-and three-pion levels, including all irreps, and including both s-and d-wave interactions. We again fix z 2 = 1, and do not include the quadratic parameter B 2 from Eq. (3.8). We call this the "minimal" complete fit.
4. An extension of the minimal fit in which we do not fix the position of the Adler zero. The goal is to check the validity of the z 2 = 1 hypothesis.
5. An extension of the minimal fit in which we keep z 2 = 1, but allow B 2 to vary in order to test whether this can improve the fit. This fit also gives information on the convergence of the q 2 expansion in the Adler-zero form.
The fits can, in principle, extend in the CMF energy up to the inelastic thresholds, which are 4M π and 5M π , respectively, for two-and three-pion spectra. In practice, we quote results from fits to a smaller energy range on two of the three ensembles to avoid possible -18 -issues with the estimate of the covariance matrix. Fits to levels within varying energy ranges yield compatible results, even though χ 2 red increases as more levels are included. This χ 2 red behavior is to be expected; we parameterize the K-matrices by expanding about the appropriate thresholds, so it is no surprise that our models become less accurate if too many high-energy levels are included in fits. This may indicate that the threshold expansions used in our parametrizations are breaking down for higher values of the CMF energy.
The results of the fits are shown in Tables 4 to 6, in each case ordered from left to right as in the list above. The tables give the number of degrees of freedom, from which the number of levels included in the fits can be seen. For example, on the N203 ensemble, we fit to 38 two-particle levels and 35 three-particle levels. In addition to quoting results for the fit parameters themselves, we also quote, in the last two rows, the results for the s-wave scattering length and effective range, obtained using Eq. (3.10), to facilitate a more direct comparison of the results of the fits.
We draw several global conclusions from the results in the tables.
a. The values of χ 2 red for the best fits are generally reasonable, although always larger than unity. Such values are typical in the analysis of lattice spectra [6,7,10], and may indicate the relevance of neglected systematic uncertainties, e.g., discretization effects, or exponentially suppressed effects in the quantization condition. Goodness of fit becomes poorer for lighter pion masses, possibly indicating a reduction in the range of validity of our truncated threshold expansions.
b. The inclusion of two-and three-particle d-wave interactions yields a better description of the data, as shown by the smaller χ 2 red compared to fits of type 2. This result is particularly striking as the levels included in this fit type are those in the trivial irreps, which are the least sensitive to d-wave interactions. Were we to attempt a fit to all irreps without including the d-wave terms, a significantly higher χ 2 red would result, as shown by the significance of the parameter D 0 in the fits in which it is included. We will further elaborate on these points when discussing fits to the kaon spectra.
c. The results for M π a ππ 0 and M 2 π a ππ 0 r ππ 0 are consistent across all five fits on all ensembles (with the fit 1 result for M 2 π a ππ 0 r ππ 0 on the N203 ensemble being the only outlier), indicating that we are obtaining a consistent description of the two-particle interactions from two-and three-pion levels. d. A comparison of the results of fits 1 and 4 suggests that the addition of the threeparticle levels improves the precision with which we can extract two-particle interactions, most notably for the d-wave parameter D ππ 0 .
e. The minimal fits (fit 3) have essentially the same χ 2 red as those of fits 4 and 5, and the results for z 2 (in fit 4) and B 2 (in fit 5) are consistent with unity and zero, respectively. Thus we conclude that the minimal description is adequate for our dataset, and use -19 -the parameters from this fit in our investigation of quark-mass dependence in the following subsection.
We close this section by showing the results of two additional fits to the N203 ensemble that motivate the choice of fits presented above. First, we test whether the s-wave parameters in K df,3 that we previously set to zero, K iso,2 df,3 and K A , are relevant for a description of our data. We use the ensemble with the heaviest pion mass for this test, since these parameters are of higher order in the chiral expansion than those we keep in our standard fits. We show in Table 7 the result of a fit in which all ten parameters in both the twoand three-particle interactions are turned on. We observe that the additional parameters do not lead to a reduction in χ 2 red , but do lead to much larger uncertainties in the fit parameters associated with s-wave interactions. Since we have not added additional d-wave parameters, we expect the results for these parameters to be unchanged, which is confirmed by the fit. Not visible from the table is the fact that the correlation between the additional fit parameters is substantial. We conclude that there is no need to include the additional parameters in order to represent our data.
Our second goal is to check whether the use of the Adler-zero form remains appropriate at heavier pion masses. To study this, we perform a minimal fit to the N203 data using the ERE form for the two-particle s-wave interaction, Eq. (3.9). We find, as shown in Table 8, that this fit is significantly disfavored compared to the Adler-zero form.    Table 5. Fit results for pions on the N200 ensemble using M π L = 4.419849, with notation as in Table 4. Each fit contains all 2π + and 3π + energy levels (in the appropriate irreps) below E * 2,max = 4M π and E * max = 5M π , respectively.

Fits to two and three kaons
We now turn to the fits of the 2K + and 3K + levels. We first use the same five fits as for pions, listed in the previous section, with results given in Tables 9 to 11. Note that we choose E * max to be slightly above the relevant inelastic thresholds, which are 2M K + M π and 3M K + M π , respectively, for two and three kaons. The values of M π /M K for our ensembles can be determined from Table 3.
Our observations concerning the fits are similar to those described above for the pion fits. In particular, the inclusion of d-wave interactions remains crucial to obtain reasonable fits, even if considering only trivial irreps. To illustrate this, we do a simple exercise using the N203 ensemble of Table 9. Taking the best fit parameters of fit 3 in the set of energy levels of fit 2, we get χ red = 1.46. This is much lower than that of fit 2, and thus a significantly better description of the levels in trivial irreps is achieved when D 0 and K B are included.
The trend of χ 2 red with the pion mass is opposite to that for the pion fits, but this can be perhaps understood because the kaon mass increases as the pion mass decreases. Another small difference is that the significance of the difference z 2 − 1, while still less than 2σ, is greater on the N203 and D200 ensembles than for pions. For this reason we choose our canonical fit in the next section to be fit 4, i.e. with z 2 left free.
As noted in the introduction, the dominant contribution to the shift in energies from -22 -their noninteracting values is due to two-particle interactions. Although our fits indicate that including the terms in K df,3 , in particular K B , leads to improved fits, it is interesting to have a visual representation of the contribution of K df,3 to the shifts. To show this, we have used the quantization conditions to determine the spectrum on the N200 ensemble taking all the parameters from fit 4 except that we set K df,3 to zero. The resulting energy shifts for the 23 three-kaon levels included in the fits are shown in Fig. 4, along with those for fits 2 and 4, as well as the results from the simulations. Note that the four levels in nontrivial irreps are not included in fit 2. We draw several conclusions. First, the shifts due to K df,3 are comparable to our present statistical errors. Second, one can see, particularly  ). For each level, the shifts determined from the lattice simulation are given by the open circles with error bars, followed to the right by the predictions of fit 2 (teal), fit 4 (orange), and fit 4 with K df,3 = 0 (red), respectively. The predictions for fit 2 are absent for the four nontrivial irreps on the right, as these levels were not included in this fit. See text for further discussion.
-23 -from the levels on the right half of the plot, that including K df,3 improves the fit more often than not. Third, for the trivial irreps (the first 19 in the plot), the difference between fits 2 and 4 is again comparable to our errors. Finally, although one cannot determine the goodness of fit from this figure, since only diagonal errors are shown, we can obtain quantitative comparisons by calculating χ 2 for the complete two-and three-kaon fits for various choices of level. We find that, if we keep only the 41 levels of fit 2, then fit 4 yields χ 2 = 52.5, while fit 4 with K df,3 = 0 gives χ 2 = 180, to be compared to 101 for fit 2. We also note that fit 4 with K df,3 = 0 yields χ 2 = 192 on the complete set of levels, to be compared to χ 2 = 64 from the full fit 4.
In ChPT, kaon interactions are generally stronger than those of pions due to their heavier mass. Thus the higher-order parameters in K df,3 might have more impact here. We have tested this with a full 10-parameter fit on the N200 ensemble, chosen as it has the greatest statistical significance for the lower-order terms in K df,3 . The resulting fit parameters are given in Table 12. Comparing to the standard fits in Table 10, we see that adding more parameters does lead to a better fit (unlike for pions), with the result for common parameters being consistent with those from the standard fits. However, the large diagonal errors, and the large correlations (not shown), indicate that there are redundant directions in parameter space, so it is difficult to extract conclusions. The only exception is that there is an indication of a nonzero K A .   (7) --50 (7) -50 (7) -51 (7 Table 10. Fit results for kaons on the N200 ensemble using M K L = 7.225143, with notation as in Table 9. Each fit contains all 2K + energy levels (in the appropriate irreps) below E * 2,max = 2.75M K and all 3K + levels (in the appropriate irreps) below E * max = 3.75M K .   Table 11. Fit results for kaons on the D200 ensemble using M K L = 9.994083, with notation as in Table 9. Each fit contains all 2K + and 3K + energy levels (in the appropriate irreps) below E * 2,max = 2.53M K and E * max = 3.53M K , respectively.

N200(K)
(2K/3K, ≤ 2) Full fit   Finally, a natural question to ask is whether is it justified to include the Adler zero in the parametrization of two-kaon interactions. This is the case at the SU(3)-flavor symmetric point, since the kaon interactions there are identical to those of pions. However, the situation when M K M π is unclear. To study this, we have performed fits on all ensembles using the standard ERE parametrization of the s-wave phase shift, both for the two-kaon spectrum alone and for the combined two-and three-kaon spectra. The results are given, respectively, in Tables 13 and 14.
The results indicate that the ensemble N203, closest to the SU(3) point, seems to be better fit by an Adler-zero parametrization. By contrast, results from N200 and D200 are equally well described by the ERE form, with the resulting fit parameters being similar. This differs from what we found with pions, where the ERE form was disfavored. One reason for this difference might be that the kaon fits are over a smaller energy range, and thus lie within the range of convergence of the ERE.

Discussion of results
We now turn to an analysis of the K-matrix parameters obtained from the fits presented in the preceding section. We compare their dependence on M π and M K to the predictions from ChPT described in Section 3.4. We discuss the parameters for pion and kaon scattering in two separate subsections.
We begin with some general comments. First, in all fits, we consider only the uncertainty in the scattering parameters (the y coordinates), treating the x coordinates (which are related to M π /F π and M K /F K ) as error-free. This is justified because the relative uncertainty in x coordinates is about an order of magnitude smaller than those in the y coordinates (see Table 2). Second, in our fit results and chiral extrapolations, we quote only the statistical uncertainty and the systematic uncertainty due to the choice of fit. We do not account for discretization errors, which are of O(a 2 ) and thus could be ∼ 5%, as is the case for the decay constants [50]. We also neglect volume dependence of the form exp(−M π L) (likely to be at or below the percent level). Thus our results have systematic errors that are not fully controlled, but break new ground by providing the first lattice results for some quantities. Finally, when we present extrapolations to the physical point, we use the physical charged pion and kaon masses (139.57 and 493.677 MeV, respectively [75]) and the corresponding decay constants (130.2 and 155.7 MeV, respectively [75]).

Results for multi-pion systems
We start with the quantities that describe two-and three-pion scattering amplitudes. In the two-pion sector, we will discuss the threshold parameters, the mass dependence of the position of the Adler zero, and d-wave interactions. We compare our results for three pions with previous work that included only s-wave interactions, and discuss the chiral behavior of the components of K df,3 .

Two-pion threshold parameters
Here we analyze the threshold parameters of two pions at maximal isospin, the scattering length, a ππ 0 , and the effective range, r ππ 0 . This system has been studied using LQCD in many earlier works [3,6,70,[76][77][78][79][80][81][82][83][84][85][86][87]. Threshold parameters have been obtained using two approaches: (i) applying the threshold expansion to order L −5 to extract the scattering length; and (ii) using the full spectrum to determine the phase shift as a function of q 2 . We use the latter approach, which is required to obtain the effective range, and uses the input from many more spectral levels.
As explained in the previous section, we take the results from the minimal fit (fit 3). We add a systematic uncertainty associated with the dependence on choice of fit, taken to be the standard deviation of the results of fits 3-5. The inputs to the chiral fits are summarized in Table 15, where we also include the leading-order ChPT predictions for comparison.
Ensemble  Table 15. Two-pion threshold parameters used for chiral fits, obtained as described in the text, along with the leading-order predictions of ChPT.
We fit these results to the SU(2) ChPT expressions, given in Eqs.    while for the effective range we obtain We plot these fits in Figure 5, which show the complete consistency with the ChPT expressions. The extrapolation to the physical point yields (M π a ππ 0 ) phys = 0.0429(1) , (M 2 π r ππ 0 a ππ 0 ) phys = 2.68(2) , (5.3) also shown in the figures. We stress again that the errors do not include all sources of systematic uncertainty, so that these numbers cannot be quantitatively compared to those from other works. However, we do note that the result for the scattering length agrees within ∼ 3% with the only fully controlled result (according the FLAG review [69]), (M π a ππ 0 ) phys = 0.0442(2)( 4 0 ) , [Ref. [77]] . (5.4)

Position of the Adler zero
In the context of LQCD determinations of the isospin-2 ππ s-wave phase shift, it was proposed in Ref. [3] to use the Adler-zero parametrization, Eq. (3.8), in place of the ERE form used previously, Eq. (3.9). Explicit inclusion of the Adler zero extends the range of validity in q 2 , and is particularly important when including subthreshold amplitudes (q 2 < 0), which is unavoidable when fitting the three-particle spectrum [13]. Indeed, it was found in Ref. [3], and subsequently supported by the results of Ref. [6], that the Adler-zero form led to better fits than those using the ERE. However, in most fits done to date, the position of the Adler zero has been fixed to its leading-order value in ChPT, z 2 = 1. Since our fits are more precise than those done previously, we can attempt to study the chiral behavior of the Adler zero, and compare to the NLO prediction from ChPT, Eq. (3.17).
We use the results from our fits in which z 2 is a free parameter (fit 4 in Tables 4 to 6 above). We find z = 8.5(2.2), χ 2 /dof = 0.6/(3 − 1) , (5.5) and show the fit in Figure 6. The extrapolation to the physical point yields (z 2 ) phys = 0.94 (2) . (5.6) Nevertheless, at the level of precision achieved here, it is reasonable to set z 2 = 1, as the magnitude of the NLO correction is smaller than the statistical uncertainty in our determination of z 2 .

Two-pion d-wave scattering length
We next analyze our results for d-wave interactions. While the s-wave interaction in the I = 2 channel has been studied extensively on the lattice, much less is known for the d-wave interaction, which to our knowledge has only been extracted in Refs. [6,67,88]. At physical quark masses, the isospin-2 d-wave phase shift exhibits an interesting feature: dispersion relations, Roy equations, and chiral perturbation theory predict a change of sign near threshold. It starts out positive (attractive) for very small q 2 and then passes through zero and becomes increasingly negative (repulsive) [66,[89][90][91][92]. In terms of k 5 cot δ 2 , this behavior implies a pole slightly above threshold. It is unclear, however, whether this phenomenon persists for heavier pion masses. Indeed, Ref. [93] studied δ 2 using ChPT, and found poor convergence for the region of pion masses of relevance here.
In practice, probing the near-threshold dependence of δ 2 using the Lüscher method is difficult for two reasons. First, δ 2 is small close to threshold due to its scaling as q 5 , resulting in tiny shifts of the energy spectrum. Second, the energy levels most sensitive to the d-wave interaction lie well above threshold for typical values of M π L currently used, and are thus sensitive to δ 2 only in the region where the phase shift is expected to be negative. For this reason we opt to use the simple one-parameter form Eq. (3.11), which implies a uniformly negative phase shift in all our fits, with D 0 < 0. As noted in Section 3.4, δ 2 -30 -vanishes at LO in ChPT. We therefore expect the d-wave scattering length to behave as up to logarithmic corrections, and choose to fit to this simple power-law dependence. As shown in Figure 7, our results are in excellent agreement with this behavior. Given the possibility discussed above of rapid changes in δ 2 near threshold, we refrain from quoting a value extrapolated to the physical point. We can make a qualitative comparison to the results of Ref. [6]. The values of the scattering lengths in the two works are of the same order of magnitude, but those from Ref. [6] do not show the same monotonic chiral behavior. The largest tension is at our heaviest pion mass, which is comparable to the heaviest pion mass used in Ref. [6]. We note that our results are based on a global fit including many two-and three-pion levels in irreps that are sensitive to the d-wave amplitude, whereas those from Ref. [6] are from just a few two-pion levels in the nontrivial irreps. Another difference is that Ref. [6] employs an N f = 2 simulation, compared to ours that uses N f = 2 + 1 sea quarks. Relative cutoff effects may also be significant, given that our lattice spacing is significantly finer. Thus we view our results as more reliable, although further work will be needed to understand the differences.

Comparison to previous s-wave three-pion fits
The three-pion coupling at maximal isospin has been extracted in the RFT approach [3,6,7], the FVU approach [10,35], and via the threshold expansion [1,9]. All these studies have included only s-wave contributions to the two-and three-pion interactions. In this section, we compare results that we have obtained from fits that closely match those used in these previous studies, i.e., restricted to s-wave interactions and using only a subset of the moving frames and levels. Specifically, we compare to the RFT results of Refs. [3,6], which use the same fit model as this work, allowing for a direct comparison. We stress, however, that, as seen in Section 4, fitting without including d-wave terms leads to poorer fits (even if only trivial irreps are included). Thus, the present comparison must be understood as a consistency check.  Table 16. Global fits to the two-and three-pion spectrum of D200 including only s-wave interactions. We use the 22 energy levels as in Ref. [3]: eleven 2π + levels (all in trivial irreps), and eleven 3π + levels, including three in nontrivial irreps.
We start with a direct comparison to Ref. [3]. That work used the two-and threepion spectra of Ref. [2], which were calculated on the D200 ensemble also used here. The present determination of the spectrum differs in several ways: (i) we increased statistics by including measurements on more gauge configurations as well as adding a second source time per configuration; (ii) we rebin data in blocks of three configurations to ameliorate autocorrelation; (iii) following these changes, we re-assessed the excited-state systematics of the spectrum extraction. The fit results are shown in Table 16, along with those from Ref. [3]. Note that this is a new fit, different from those presented in Section 4.1. We stress that we are fitting to the same set of levels, so differences in fit parameters result solely from updates in the energies of the levels. We observe a substantial difference between the two fits, which, including correlations, is about 3σ. In particular, the determination of M 2 π K iso,0 df,3 , which was found to be different from zero at the 2σ level in Ref. [3], is called into question in view of the change of sign in the new fit.
Our interpretation of the discrepancy is the following. First, a larger rebinning reduces the impact of the autocorrelation between samples, and leads to more reliable determinations of the energy levels. Second, the increase in statistics enables the use of larger values of t min , which reduces the contamination from excited states. In summary, we now think that there may be systematic errors in the spectrum of Ref. [2] that have not been accounted for. Since the impact of three-particle interactions on the spectrum is small, such effects can lead to significant systematic errors in the parameters in K df,3 . Note that in our preferred fit which includes d-wave interactions, the value of M 2 π K iso,0 df,3 is shifted by more than one standard deviation-compare Table 16 and fit 3 in Table 6. This is another indication of the challenge in determining these parameters.
We end this discussion with a plot comparing the determinations of the s-wave components of K df,3 from this work with those from Ref. [6], as well as Ref. [3]. Here we use the results of the full s-wave-only fits (fit 2), since these use the same fit functions as Refs. [3,6]. There are however some differences between our fit 2 and Refs. [3,6]: (i) fit 2 includes only trivial irreps, while Refs. [3,6] include both trivial and nontrivial irreps in the three-pion sector, and (ii) here we use more moving frames than in previous work. The comparison is shown in Figure 8, which also shows the LO predictions from ChPT from Eq. (3.21). For this figure alone we use M π a ππ 0 for the x axis, as this has been used in the prior works. As a consequence the error bars now become ellipses.
-32 -We first note that, for the two higher-mass ensembles, most of our results (red ellipses) are statistically different from zero, unlike for the D200 ensemble. The figure also shows that the just-discussed tension with Ref. [3] (blue ellipses) remains when we use a fit to levels in the trivial irrep in all frames. What we see in addition is a significant tension at the heaviest pion mass with the results from the ETMC collaboration [6] (orange ellipses). Just as for the difference in the d-wave scattering length, it is difficult to understand this tension. Possible sources of difference are that the two calculations correspond to different lattice regularizations, and thus have different O(a 2 ) discretization errors, and the fact that the ETMC result is for N f = 2, while our N203 ensemble is close to the SU(3)-symmetric point. More investigation will be needed to understand how these effects enter in K df,3 .  Tables 4 to 6. Results for Refs. [3,6] are also shown.

Results for K df,3 from full fits including d waves
We conclude the discussion for pions by presenting our final results for the terms in K df,3 , which are obtained from fits to levels in all irreps including d-wave interactions. We recall that such fits provide a much better description of the spectrum than fit 2, the results from which are discussed in the previous section for comparative purposes. We have results for K iso,0 df,3 , K iso,1 df,3 , and K B , and, as before, use fit 3 for our central values. The chiral dependence of the isotropic parameters is shown in Figure 9. We observe statistically significant deviations from zero at heavier pion masses. Our results are consistent with a linear chiral dependence on (M π /F π ) 4 as expected from ChPT, see Eq. (3.21). However, the respective slopes are in significant tension with the ChPT predictions, hinting at a possible breakdown in the convergence of ChPT. Indeed, as noted already in Ref. [3], there are reasons to expect large NLO corrections to K df,3 . A final comment concerns the results extrapolated to the physical point: since K df,3 vanishes rapidly towards the chiral limit, it will be very difficult in practice to extract the three-pion interaction directly at this point.  linear fit This work Physical point Figure 10. Chiral scaling of K B for pions as a function of (M π /F π ) 6 , including a linear fit. Notation as in Figure 5.
We now turn to K B . As explained in Section 3.4, while no ChPT prediction is available, we expect the chiral scaling up to logarithms. This expectation is borne out by our results shown in Figure 10. The extraction of K B , despite being of quadratic order in the threshold expansion, is aided by it being the only contribution from K df,3 to a set of energy levels in nontrivial irreps. Similarly to the isotropic parameters, the extrapolation indicates a very small value for K B at the physical point, which will presumably be difficult to extract directly from simulations with close-to-physical quark masses.
Ensemble  Table 17. Two-kaon threshold parameters from fit 4 in Tables 9 to 11, with systematic errors due to choice of fit obtained as discussed in the text.  Figure 11. Results for the s-wave two-kaon scattering length and effective range. Notation as in Figure 5, except that in the left panel we include, with an empty triangle, the result for the chiral extrapolation from an SU(3) ChPT fit to both kaon and pion scattering lengths.

Results for multi-kaon systems
Unlike for the isospin-2 ππ system, there are relatively few studies of two kaons at maximal isospin [8,86,94,95]. Furthermore, the three-kaon K-matrix has yet to be explored, as the only other study of three kaons at maximal isospin set it to zero [8]. Thus we provide here the first exploration of this quantity. The remainder of this section has a similar structure to that for pions. As discussed in Section 4.2, we use fit 4, which includes the position of the Adler zero as a free parameter, as our reference fit.

Two-kaon threshold parameters
As in the two-pion case, we start by looking at the two-kaon threshold parameters. The results are summarized in Table 17, along with the LO chiral predictions. The central values are from fit 4 in Tables 9 to 11, while the systematic uncertainty is given by the standard deviation of the results from all the fit models, including the ERE fits on the D200 and N200 ensembles. The deviation from LO ChPT is more pronounced for kaons than for pions (Table 15), indicating poorer convergence of the chiral expansion.
The results for the scattering length and effective range are shown in Figure 11. For both quantities we perform a linear fit in (M π /F π ) 2 , motivated by SU (2)  where the quoted uncertainties are only statistical. We can also perform a simultaneous SU(3) ChPT fit to the kaon and pion scattering lengths using the NLO results Eqs. (3.19) and (3.20). We stress that a single, common LEC enters in both expressions, so that we have one free parameter to describe six data points. In the fit, we ignore the (small) uncertainties in (M π /F π ) 2 -as we do for all fits presented in this section-as well as the correlations between a KK 0 and a ππ 0 evaluated on the same ensemble. The resulting percent-level determination of L ππ (µ = 4πF π ) = −1.13(3) · 10 −3 , χ 2 /dof = 0.98/(6 − 1), (5.10) does not reflect that there is a family of NLO ChPT expressions, which differ only by higher-order terms, but lead to slightly different results for L ππ . Here we use the form in which both equations are written as functions of just x π and x K , but it would be equally valid, for example, to replace F K → F π in the NLO terms. Comparing the results from fitting with different allowed NLO forms, we find that the resultant L ππ varies at the 10% level. We thus assign a systematic error of this size, which we interpret as due to missing NNLO terms. Including this uncertainty, the results of extrapolating to the physical point are The kaon result is also shown in Figure 11a. We see that the statistical error is dominated by the systematic error from NNLO effects. The result for the physical pion scattering length in Eq. (5.11) is in complete agreement with that from the SU(2) ChPT extrapolation, given in Eq. (5.3). For the kaon scattering length, however, the result from the linear fit, Eq. (5.9), disagrees significantly with that based on SU(3) ChPT, Eq. (5.11). This difference can be seen also in Figure 11a.
We think that this roughly 10% difference is mainly due to discretization effects. In particular, we note that, along our chiral trajectory, the ratio M K /F K is expected to increase monotonically towards physical quark masses, yet its value on our most chiral ensemble 4.513(9) (see Table 2) is larger than the physical value 4.472. Indeed, sizeable discretization effects of 3 − 4% were observed in similar ratios in Ref. [50]. Since the expression for M K a KK 0 , Eq. (3.20), is proportional to (M K /F K ) 2 , the discretization errors in the SU(3) physical-point prediction could be as big as 6 − 8%, and thus largely explain the discrepancy between the extrapolations based on SU(2) and SU(3) ChPT.

Two-kaon d-wave interactions
We now turn to the d-wave two-kaon interaction, which to our knowledge has not been previously studied in lattice calculations. In our fits, this is described by a single parameter, D 0 , which is constrained quite well by both two-and three-kaon levels in nontrivial irreps.
-36 -Our results for the d-wave scattering length, M 5 K a KK 2 , are shown in Figure 12. Since we do not know of an SU(2) ChPT expression for this quantity, we perform a simple linear fit, which represents the data rather poorly. The conclusion is that there is some evidence that the d-wave scattering length increases as one approaches the physical point. . Notation as in Figure 5.

Results for K df,3 for three kaons
Our final results are for the parameters of the three-kaon K-matrix: K iso,0 df,3 , K iso,1 df,3 , and K B , which have not previously been determined. In order to simplify the notation, we use the same symbols for these parameters as for pions, although they are different physical quantities.
The results for K iso,0 df,3 and K iso,1 df,3 are shown in Figure 13. From the left panel, we see that K iso,0 df,3 is consistent with zero for all masses, and also after linear extrapolation to the physical point. For K iso,1 df,3 , by contrast, nonzero values are found for the two heaviest pion masses, and a linear extrapolation to a nonzero value at the physical point is reasonable. The only theoretical guidance we have is from the LO ChPT result of Eq. (3.21), which predicts proportionality to M 4 K /F 4 K . Given our chiral trajectory, this would lead to an increase as we move to smaller values of M 2 π /F 2 π . However, the values predicted by LO ChPT are very far from those we find. For instance, on the N203 ensemble the results are which do not even appear in the plot ranges.
Our results for K B are shown in Figure 14. As can be seen, the result on the D200 ensemble has by far the largest statistical errors. This is because the large value of M K L ∼ 10 suppresses the contribution to shifts in the spectrum from interactions that are of higher order in the threshold expansion. Our results are consistent with an increase in K B as we approach the physical point.

Conclusion
In this work we have extended the application of LQCD to multihadron systems by utilizing state-of-the-art numerical methods to determine an order of magnitude more two-and three-particle spectral levels than in previous work. On each ensemble we have determined, for both pions and kaons, 50 − 80 levels below the relevant inelastic thresholds, roughly equally split between those for two and three particles. The energies of the levels have been determined to a precision of 1 − 5%, with the shifts from the corresponding free energies determined with errors of 5 − 15%. This unprecedented number of levels with such high precision has been made possible with stochastic LapH and advanced contraction algorithms. The jackknife samples of all energies extracted in this work are provided in HDF5 format as ancillary files with the arXiv submission We have found that these levels can be described well by the two-and three-particle quantization conditions using a relatively small number of underlying parameters. Since the major contribution to the energy shifts to three-particle levels arises from two-particle interactions, our strategy has been to do a simultaneous fit to two-and three-particle levels, including all correlations. We find that this leads to a better determination of the twoparticle scattering parameters than a fit to the two-particle levels alone. The goodness of the fits is reasonable, with χ 2 red = 1.2 − 1.7. The 3π + and 3K + systems are nonresonant, with repulsive two-particle interactions that are expected to have mild energy dependence. This allows the two-particle interactions to be described by relatively few parameters in the relevant energy range. Our fits have been able to determine these parameters accurately. For example, the pion and kaon scattering lengths are determined with statistical errors of 1 − 2%. These systems thus provide a good testbed for studying the significance of three-particle interactions, which themselves make only a small impact on the energy levels. Previous work has either found no statistical evidence for such interactions, or a barely significant signal.
Our main conclusions are as follows. First, we find that a parametrization of the two-pion K-matrix that includes the Adler zero expected from chiral symmetry is favored, compared to an effective range expansion. For our lightest kaons, the presence of the Adler zero is also preferred. In this regard, we stress that the three-particle quantization condition involves contributions from subthreshold two-particle scattering, and thus is more sensitive to the Adler zero. We also find that the position of the Adler zero is consistent with the expectations of ChPT, albeit with relatively large errors.
Our second conclusion is that reasonable fits require the presence of d-wave parameters, both for two-and three-particle interactions. In the two-particle sector, we determine the d-wave scattering length with 10 − 20% statistical errors, and observe the expected chiral behavior. For three particles, our results provide the first determination of the d-wave three-particle parameter, K B . We note that, once d-wave parameters are included, there is no longer a one-to-one connection between levels and scattering parameters, even in the two-particle sector. Thus a global fit is required.
Our final conclusion is that the determination of three-particle interactions requires the use of many energy levels, in a variety of irreps, determined with sufficient accuracy. In particular, to extract d-wave contributions it is important to use levels in nontrivial irreps. We find nonzero results not only for K B , but also the s-wave parameters K iso,0 df,3 and K iso,1 df,3 in most of the cases. The statistical significance of these results vary, but is greater than in previous studies, exceeding 3σ in several cases. For pions, the s-wave part of K df,3 has the expected linear dependence on M 4 π , but the coefficient is in disagreement with leading order ChPT. The d-wave part K B also has the expected M 6 π dependence, but in this case a ChPT prediction for the coefficient is not available. For both quantities our chiral extrapolations indicate that it will be very difficult to determine the three-particle interaction for physical quark masses. We also observe some tensions between our results for K iso,0 df,3 and K iso,1 df,3 and those of previous work [6]. A noteworthy feature of our fits using the quantization conditions is that they continue to work quite well above the energies that we have included in the fits (see Figure 3). This shows that we have not forced the fits to work in a limited energy range only to have them quickly fail outside that range, and thus supports the applicability of the threshold expansion that we have used for the K-matrices. In particular, for three kaons, these extra levels lie outside the range that can be rigorously described by the quantization condition, indicating that, as expected, anomaly-induced transitions and overlaps onto mixed-flavor states are suppressed.
There are many exciting directions in which this work can be extended. Most straightforward is to the maximal charge systems with hybrid flavor content, e.g. π + π + K + and π + K + K + . The formalism for analyzing such "2 + 1" systems has very recently been developed [25]. More challenging are systems of three pions or kaons with isospin less than the maximal value, which involve quark annihilation diagrams and, in some cases, resonant behavior. The relevant formalism has been derived in Ref. [21]. We hope to study such systems in the near future.
Another future direction is to extend the study of 3π + and 3K + interactions to ensembles with different lattice spacings and volumes. This will provide a direct check on the importance of discretization errors, and on whether our fits can accurately describe the volume dependence. The former have been shown to be relevant in two-baryon systems [62]. In the latter regard, we note that the formalism that we use drops corrections proportional to e −mπL , and there is no expectation in the short term that such terms could be included in a rigorous way. Thus we must assume that they are small, and the only consistency check that we have on this is to compare results at different volumes.
Finally, we note that, to obtain the three-particle S-matrix elements from the Kmatrices determined here, one needs to solve the integral equations presented in Ref. [15]. To do so in the presence of d-wave interactions requires a generalization of the methods developed in Refs. [7,96]. also acknowledges financial support from Generalitat Valenciana through the plan GenT program (CIDEGENT/2019/040).
Calculations for this project were performed on the HPC clusters "HIMster II" at the Helmholtz-Institut Mainz and "Mogon II" at JGU Mainz. We are grateful to our colleagues in the CLS initiative for sharing ensembles.

A Little-group representation matrices
Our approach to designing the single-hadron and multi-hadron operators used in this work has been described in detail in Ref. [44]. However, single-hadron operators with momenta in directions such as (0, 1, 2) and (1, 1, 2), where Cartesian components are used, were not treated in Ref. [44]. Hence, we provide specific details concerning only these additional operators in this appendix.
For each class of momenta, we choose one representative reference momentum direction p ref . We then construct operators that transform irreducibly under the little group of p ref .
Recall that the little group of p ref is the subset of the symmetry operations that leave the reference momentum p ref invariant. For each momentum direction p in a class of momenta, we select one reference rotation R p ref that transforms p ref into p. As long as the selected rotation transforms p ref into p, it does not matter which group element is chosen, but a choice must be made in order to specific all phases between the single-hadron operators of different momentum directions. These phases must be known in order to properly construct the multi-hadron operators. All single-hadron operators having a momentum in the direction of p are then obtained by applying the reference rotation to the corresponding operators constructed using the momentum in the direction of p ref . Our choices of reference momenta directions and reference rotations for the additional operators used in this work are listed in Table 18.
The little groups associated with the additional momentum directions used in this work are listed in Table 19. Although this work involves only mesons, we wish to provide details concerning both the single-valued and double-valued representations for possible future calculations with baryons. The double-valued representations of a group G are constructed by extending the group elements to form the so-called "double group" G D . The elements of the double groups associated with our choices of additional reference momentum directions are explicitly presented in Table 19, grouped into their conjugacy classes.
The characters of the irreducible representations (irreps) of the little group C s for the additional momenta directions are presented in Table 20. The one-dimensional singlevalued irreps are labeled by A, and the one-dimensional double-valued irreps are denoted by F . Since all of the irreps are one dimensional, the characters are the representation matrices.

B Tables of interpolating operators
The three-and two-hadron operators used in this work are listed in Tables 21 to 31 below. The notation for the irreps follows the conventions in Ref. [44]. The subscripts g/u denote  [44] for definitions of the rotation operators C nj below. E is the identity element. even/odd parity, and the superscripts +/− denote even/odd G-parity. The Clebsch-Gordan coefficients that fully define each operator are not given, but are available upon request. Table 19. The little groups corresponding to reference momentum directions (0, 1, 2) and (1,1,2) are isomorphic to C s . The elements of the double groups C D s for these momenta directions are listed above, grouped into conjugacy classes. E is the identity element, E represents a rotation by 2π about any axis, and G = EG for each element G in C s . Spatial inversion is denoted by I s .  Table 19 for the definitions of the classes C n . Since all of the representations are one dimensional, the characters are the representation matrices.            -51 -  Table 32. Energy levels used in the fits of this work. Notation is as follows: "2A − 1u + E − u " means the lowest two levels in the A − 1u irrep, and the lowest in the E − u irrep.

C Energy levels used in fits
-52 -