Interactions of $\pi K$, $\pi \pi K$ and $KK\pi$ systems at maximal isospin from lattice QCD

We study the interactions of systems of two and three nondegenerate mesons composed of pions and kaons at maximal isospin using lattice QCD, specifically $\pi^+K^+$, $\pi^+\pi^+K^+$ and $K^+K^+\pi^+$. Utilizing the stochastic LapH method, we determine the spectrum of these systems on two CLS $N_f=2+1$ ensembles with pion masses of $200$ MeV and $340$ MeV, and include many levels in different momentum frames. We constrain the K matrices describing two- and three-particle interactions by fitting the spectrum to the results predicted by the finite-volume formalism, including up to $p$ waves. This requires also results for the $\pi^+\pi^+$ and $K^+ K^+$ spectrum, which have been obtained previously on the same configurations. We explore different fitting strategies, comparing fits to energy shifts with fits to energies boosted to the rest frame, and also comparing simultaneous global fits to all relevant two- and three-particle channels to those where we first fit two-particle channels and then add in the three-particle information. We provide the first determination of the three-particle K matrix in $\pi^+\pi^+K^+$ and $K^+ K^+ \pi^+$ systems, finding statistically significant nonzero results in most cases. We include $s$ and $p$ waves in the K matrix for $\pi^+ K^+$ scattering, finding evidence for an attractive $p$-wave scattering length. We compare our results to Chiral Perturbation Theory, including an investigation of the impact of discretization errors, for which we provide the leading order predictions obtained using Wilson Chiral Perturbation Theory.

multihadron systems. As the complexity of the system grows, the number of quantities to constrain from LQCD becomes larger. Indeed, a qualitative step can already be seen in this work compared to that for identical particles: we require more terms in K df, 3 , and two different two-meson amplitudes. In particular, the question arises whether it is preferable to determine the two-meson amplitudes from results for the two-particle spectrum (e.g. for π + π + and π + K + ), and then use the results in a fit to the three-particle spectrum (e.g. for π + π + K + ), or, instead, do a global fit to all channels at once. This is but one example of a general issue. An extreme case is provided by the determination of the γ * γ * → ππ amplitude using finite-volume methods: the formalism requires also input from the finitevolume processes ππ → ππ, πγ → π, γ → ππ, and ππγ → ππ [55]. To determine which fitting procedures are preferable for more complex fits, we investigate and compare several fitting strategies.
Byproducts of this work are well-determined two-particle amplitudes. In particular, we compute the isospin I = 3/2 πK scattering amplitude in both s and p waves. By comparing our results for the s-wave scattering length to ChPT, we are able to extract low energy constants (LECs) in a threshold expansion. Our results for the I = 3/2 p-wave scattering length are at lower pion masses than that obtained previously in Ref. [56]. We compare our result with this LQCD calculation, with the ChPT prediction, and with the results from a dispersive analysis.
This work also contains the first estimates of discretization errors in the three-particle K matrix. To achieve this, we present a new calculation of K df,3 using an extension of continuum chiral perturbation theory (ChPT) in which the effects of discretization errors are included, so-called Wilson ChPT (WChPT) [57,58]. This allows us to determine the dependence of K df,3 on the lattice spacing in terms of low energy constants. We also calculate the corresponding dependence for two-particle amplitudes, which allows a determination of the relevant LECs by fitting to our results for these amplitudes. Thus we can estimate the magnitude of discretization effects in K df, 3 . This paper is organized as follows. In section 2 we provide details of the LQCD simulation, describe our choice of interpolating operators, and discuss the fitting and results for the single and multiple-hadron spectra. Then, in section 3, we review the finite-volume formalism required for this work, present the parametrizations of the K matrices that we use, and collect the results from ChPT that we need. Next, in section 4, we describe the strategies that we use to fit to the spectra, and the results obtained, comparing different approaches. We collect our final results for infinite-volume scattering parameters in section 5, and compare them to ChPT, extracting several low-energy coefficients. We conclude in section 6. We include four appendices. Appendix A displays a table with the energy levels used in fits, appendix B provides details on the calculation of discretization effects using WChPT, appendix C sketches the derivation of the analytical result for the p-wave πK scattering length at NLO in ChPT, and appendix D summarizes the sets of operators used in this work.

Finite-volume Spectrum Extraction
Here we present the details regarding the extraction of the finite-volume spectrum from the two-point correlation functions computed using LQCD. The computational methods and strategies used follow those laid out in Ref. [49]. However, for convenience, we review the pertinent details.

Computation of Correlators
The extraction of finite-volume energies proceeds by first calculating two-point temporal correlation functions, which can be seen to contain all the information on the spectrum through their spectral decomposition where O i (t) and O j (t) are annihilation and creation operators, 2 respectively, |Ω corresponds to the vacuum state, and E n is the energy of the nth eigenstate |n of the Hamiltonian. The indices on the operators indicate the possibility of a set of linearly independent interpolators that all share the same quantum numbers. Details on the procedure for extracting the spectrum from these correlators will be given in later sections, and we now turn to how the correlators themselves are computed. The multi-hadron interpolating operators we consider in this work (see section 2.2 and appendix D) involve definite-momentum projections for the individual hadrons, which in turn necessitates quark propagators from all spatial sites at a given source time to all other lattice sites. One such method, referred to as distillation, acheives this with a smaller computational cost by utilizing a particular smearing based on the covariant Laplacian that cuts off higher-lying modes within the quark fields [59]. The smeared quark propagators can then be obtained by performing inversions within the much smaller subspace spanned by the N ev retained eigenvectors of the covariant Laplacian. However, the number of required eigenvectors needed to keep a constant smearing radius grows proportionally with the volume L 3 and can quickly become prohibitively expensive for large volumes. To mitigate this issue, rather than use distillation directly, we use a stochastic variant, referred to as stochastic LapH which has a better cost scaling [60]. This strategy combines stochastic sources within the distillation subspace and the method of dilution [61] to estimate the smeared quark propagators in terms of source (r,d) and sink ϕ (r,d) functions where r labels the N r stochastic sources, d labels the dilution partition, the columns of V s contain the N ev retained eigenvectors, P [d] is a dilution projector, S = V s V † s is the smearing kernel, D is the Dirac matrix, and (r) is a noise vector within the distillation subspace. -4 -From the quark sinks and sources we form the meson sinks and sources, respectively, which are rank-2 tensors in the dilution indices. The final correlators are constructed from contractions of these tensors over the dilution indices. We make use of common subexpression elimination [62] to reduce the number of needed contractions and diagram consolidation to speed up the optimization. These algorithms have also been utilized for two-baryon [63] and meson-baryon [64] systems. Further details on our implementation of these contraction speedups can be found in Ref. [41], the code for which has been made publicly available [65].

Interpolating Operators
The single-hadron interpolators we use, which annihilate the single-hadron states, are (2.5) where u, d, s are up, down, and strange quark fields, respectively. These are then substituted into our two-and three-hadron interpolating operators, which are of the form

4)
respectively, where P = i p i is the total momentum, Λ is the irrep of the little group of P, f i labels the flavor of the individual hadrons, and c P,Λ are Clebsch-Gordan coefficients. The Clebsch-Gordan coefficients are determined by requiring the overall operator transform according to the irrep Λ. The sum on the right-hand-side includes all momenta related via rotations within the little group of P, with the constraint that i p i = P. A table listing the multi-hadron interpolating operators used in this work is given in appendix D.

Lattice details
Our calculations are performed on two ensembles generated by the CLS (Coordinated Lattice Simulations) consortium [66]. The ensembles use N f = 2 + 1 flavors of nonperturbatively O(a)-improved Wilson fermions and the tree-level O(a 2 )-improved Lüscher-Weisz gauge action. The bare quark masses are tuned such that they follow a chiral trajectory in which the sum of the quark masses is held fixed. The practical effect of this choice is that the kaon mass approaches its physical value from below as the pion mass approaches its physical value from above. The two ensembles used in this study, named N203 and D200, both share a lattice spacing of a ≈ 0.06426 (76) fm, which was determined from the linear combination 2 3 (F K + 1 2 F π ) of decay constants [67]. 3 Other details, including the stochastic LapH setup, on these two ensembles can be found in table 1.
(L/a) 3 × (T /a) M π [MeV] M K [MeV] N cfg t src /a N ev dilution N r ( /s) N203 48 3 Table 1. Specific details on the ensembles used in this work, including the name, geometry, approximate pseudoscalar masses, number of configurations N cfg , source positions t src used, number of eigenvectors N ev of the covariant Laplacian retained, dilution scheme (see Ref. [60] for details), and number of noises N r used for the light (l) and strange (s) quark sources. Both ensembles have the same lattice spacing a ≈ 0.063 fm.
In order to prevent topological charge freezing at fine lattice spacings, open temporal boundary conditions are employed [70]. This requires care when choosing the temporal locations of the source and sink times used for the correlators, as we must make sure that they are sufficiently far from the temporal boundaries in order to suppress any effects from the boundary. As there was no need to produce additional quark sinks beyond what was used in our previous study [49], the arguments used there to justify the source and sink positions carry over here. Essentially, evidence for sufficient suppression of boundary effects on the D200 ensemble was given in Ref. [71], where it was found that a temporal distance of ∼ 32a from the boundary was enough for the exponentially decaying boundary effects to be negligible. Further, it is expected that the boundary effects are more severe on D200 than N203, as the leading contribution comes from the lowest state with quantum numbers of the vacuum, which should be a two-pion state for the quark masses considered here, and therefore has a smaller energy on D200. Thus, as the source positions considered for N203 are even further from the boundary than D200, our choices should be safe from the effects of the open boundary conditions. Note that the source position of t src = 92a for D200 only has sink times smaller than 92a associated with it (i.e. the correlators go backward in time, see Ref. [49] for more details).
Finally, autocorrelations, which lead to underestimated errors, can be checked for by observing dependence on the error estimates from averaging N rebin successive configurations across all the original measurements into N cfg /N rebin new bins. While there is evidence that values as high as N rebin = 20 are needed for D200 to completely remove autocorrelations [64], this is not plausible for our use-case, as the number of energies used in our fits in section 4.2 is too high to reliably estimate the covariance matrix with so few bins. However, we have found little to no dependence on the final results for N203 when using N rebin = 1 or N rebin = 3, suggesting the observables of interest here are not affected significantly by autocorrelations. We therefore use N rebin = 1 for N203, while using N rebin = 3 for D200 in order to still obtain reliable estimates for the covariance matrix while removing some autocorrelation. Additionally, we note that the configurations used on N203 are separated in Markov time by twice the distance used for D200, which is why we use a conservative choice for the rebinning on D200.  Table 2. Pion and kaon masses and decay constants for the two ensembles considered in this work. The masses were determined in our previous work [49]. The decay constants were determined in Ref. [72].

Finite-volume energies from correlators
As can be seen from the spectral decomposition in eq. (2.1), in principle one can extract any energy so long as the operators used have non-zero overlap onto the corresponding eigenstate. However, with finite statistics, reliably determining the states beyond the first few terms from fits to a single correlator is difficult. As we are only interested in the ground states for the single-hadron particles, we can obtain these from single-exponential fits to the correlators starting after all higher-lying states are exponentially suppressed. Fortunately, for the single hadron correlators, the signal-to-noise ratio is either constant or slowly decaying which allows for a good signal after all higher-lying states have decayed away. The needed single-hadron masses were determined in our previous work [49] and are reproduced in table 2 for convenience, along with the decay constants needed for the chiral extrapolations performed later. Our eventual goal is to obtain the multi-hadron interactions, which are constrained by the multi-hadron finite-volume energies. Each of these energies provides a constraint on these interactions, and, therefore, including as many energies as possible will typically improve the reliability of the extracted interactions. In general, the gaps between the multi-hadron states are smaller than those of the single-hadron states, making a reliable extraction more challenging, especially coupled with the increased statistical errors of the multi-hadron fits. Thus, we need a method to reliably determine several energies from the multi-hadron correlators. We utilize a variational method [73,74], which relies on solving a generalized-eigenvalue problem (GEVP) using a correlator matrix built from sets of interpolating operators that all transform in the same way. This is the same method used in Ref. [49], but we repeat most of the details here to make the discussion self-contained.
For each overall flavor, irrep, and total momentum squared, we construct a set of N operators and use them to calculate a correlation matrix C ij (t sep ) as in eq. (2.1). We then form a GEVP as where t 0 is referred to as the metric time. As long as t 0 ≥ t/2, one can show that the generalized eigenvalues behave as where n = 0, . . . , N − 1, E n is the nth eigenenergy, and ∆ n ≡ E N − E n . Thus, the method provides a straightforward way of extracting the excited-state energies without having to -7 -perform a fit which includes many exponentials. In fact, the eigenvalues of C(t) also share this advantage, but with a gap ∆ n ≡ min m =n |E n − E m | that is smaller than or equal to the gaps found in the case of the GEVP. Therefore, the gap from the GEVP helps to further suppress unwanted contributions in the generalized eigenvalues. One technical complication involves matching the eigenvectors at one time separation to another or one resampling to the next (i.e. eigenvector pinning), which can easily become problematic from ambiguous choices, especially at large time separations where the noise grows. Instead, we solve the GEVP on the mean and at a single time separation t d , where t 0 < t d 2t 0 , and then use the eigenvectors to rotate the original correlator matrix on all resamplings and at all other times not equal to t d where the parentheses indicate an inner product. To avoid any systematics associated with only solving the GEVP for one time separation, we look for stability in the extracted spectrum as the diagonalization time t d and metric time t 0 are varied.
To obtain the two-and three-hadron finite-volume energies, we use single-exponential correlated-χ 2 fits to the ratios is an average of single-hadron correlators with flavor f i over all rotationally equivalent momentum with magnitude p 2 i . The product of single-hadron correlators in the denominator is chosen based on the expected non-interacting energy level associated with the nth eigenstate, in which case the asymptotic behavior is where ∆E n lab is the energy shift of the nth eigenstate from its non-interacting value in the lab frame. The use of the ratio has a few advantages over fitting directly toĈ n (t): there is a strong cancellation of correlated fluctuations between the numerator and denominator, and the plateau in the effective energy of the ratio begins at earlier time separations and is more stable across several time separations. However, despite the earlier plateau seen from the ratio, in order to not introduce any further systematics, we typically make a conservative choice for the beginning of our fit range t min such that the single-hadron correlators have already attained their asymptotic behavior.
In figure 1, we show the dependence of several extracted energy shifts on t min and the GEVP parameters t 0 and t d . The energies shown correspond to the first excited state in an irrep with P 2 = 2(2π/L) 2 and the ground state in an irrep with the largest momentum-squared used for a given flavor. We include all three mixed-flavor systems and both ensembles. As can be seen in these t min plots, typically there is a wide region of t min with consistent energy shifts before correlated fluctuations take over when the signal starts to be lost. The ability of the variational method to suppress contributions from unwanted nearby states is illustrated by the results with P 2 = 2(2π/L) 2 , for in these cases -8 -  Figure 1. The dependence of the energy shift in the lab frame extracted from fits to the ratio R n (t) on the smallest time separation t min included in the fit. The colors correspond to different values of the GEVP metric time t 0 and diagonalization time t d : t 0 /a = 4, t d /a = 8 (blue), t 0 /a = 6, t d /a = 12 (orange), t 0 /a = 12, t d /a = 24(green). A small horizontal offset is applied to separate these three choices. Each plot indicates the ensemble, irrep, total momentum-squared (in units of (2π/L) 2 , energy level, and flavor. For example, in the bottom-right plot, the notation A 2 (2) indicates the A 2 irrep with momentum squared being 2(2π/L) 2 , while E ππK 1 indicates that this is the first excited π + π + K + level in this irrep. In each panel, the black horizontal solid and dashed lines indicate the mean and error, respectively, of the energy shift for the chosen fit, for which the value of t min is indicated by the vertical dashed line. there are usually several nearby energy levels. Additionally, the dependence on the choices for (t 0 , t d ) is either very mild or not visible. 4 The value of t min is chosen such that it is much larger than the onset of the stability in the extracted energy, but not so large that correlated fluctuations begin to arise. Among the fits satisfying this critera, our final value is based on making a conservative choice to ensure that any systematics are smaller than the statistical errors, while also making sure the fit quality is reasonable. Examples of the choices of t min are shown in the figure.
To illustrate the number of levels that we are able to determine, and the errors that we obtain, we show in figure 2 the energy levels for the ππK system on the N203 ensemble. To better compare the levels from different momentum frames, we show the center-of-mass frame (CMF) energies E * = E 2 − P 2 . Also shown (as teal dots) are the result of our standard fit to these levels, to be discussed below, in which we fit only to levels that lie below the first inelastic threshold. Although the formalism is not, strictly speaking, valid above this threshold, we also display its predictions for some of the higher levels (as orange dots). Analogous plots for the other three-meson systems that we consider are shown in appendix A. Detailed discussion of these and other fits are provided below.

Theoretical background
In this section we collect theoretical results needed to constrain the two-and three-particle K matrices from the two-and three-particle spectra obtained using lattice QCD. As discussed in the previous section, we have results for both degenerate and nondegenerate twoand three-particle channels. The formalism for the degenerate cases has been reviewed in Ref. [49], so we focus on the formalism needed for the new channels considered here, i.e. πK for the two-particle case, and ππK/KKπ for three particles. We first provide a brief recapitulation of the quantization conditions and our fitting strategy, then describe the parametrizations of K matrices that we use, and finally collect the predictions of ChPT.

Quantization conditions and fitting strategy
To extract infinite-volume scattering parameters from a finite-volume energy spectra, we make use of both two-and three-particle quantization conditions. The two-particle finitevolume formalism was first developed by Lüscher [75,76]. We will need the generalizations to moving frames and to nondegenerate particles given in Refs. [77][78][79].
The three-particle formalism has been developed using three approaches: the relativistic field theory (RFT) approach [14,15], which will be the basis of this analysis, the non-relativistic effective field theory (NREFT) approach [17,18] (subsequently relativized in Ref. [34]), and the finite-volume unitarity (FVU) approach [40]. The formal equivalence (up to technical differences) of the RFT and FVU formalisms was established in Ref. [28], and the equivalence of FVU and relativized NREFT formalisms is noted in Ref. [34]. Reviews of these approaches and comparisons between them can be found in Refs. [3,4,6,8], and a direct comparison of their application is given in Ref. [50]. We will use the RFT result for three identical particles given in Refs. [14,15], and the extension given in Refs. [27,31,33] to so-called "2+1"-systems, i.e. those involving one distinct and two identical spinless particles.
We will not recap the derivations here. We only note that the RFT method is based on an all-orders diagrammatic analysis in a generic relativistic effective field theory, and applies in a kinematic regime in which only three-particle channels can go on shell. As one might expect from the name, this approach applies for relativistic particles; and if the relativistic forms of kinematic functions are used (see Ref. [24]), the formalism leads to Lorentz-invariant scattering amplitudes.
An extensive discussion of the implementation of the quantization condition for 2+1 systems is given in Ref. [54], and thus we present only an overview here. In particular, quantities not defined here can be found in that work. We also make use of the formalism for identical particles derived in Ref. [14], and in particular the implementation presented in Ref. [24]. A general feature of the three-particle formalism is a division into a spectator particle and the remaining pair or dimer. For practical applications, one must impose a cutoff on the angular momentum, , of the pairs. In our previous work analyzing 3π + and 3K + spectra, we used max = 2, thus including s and d waves (p wave being forbidden for identical particles) [49]. However, working up to max = 2 is not possible here, due to the proliferation of fit parameters, as discussed in Ref. [54] and recapitulated below. Instead we use max = 1, with s and p waves in the π + K + subsystems, and only s waves for subsystems with identical particles (2π + and 2K + ).
We begin by considering the two-particle quantization condition. The inputs are the total momentum, P , the box size, L, a kinematic function denoted F , which contains the effects of finite-volume physics, and the two-particle K matrix, K 2 . By solving the following two-particle quantization condition, we can determine the two-particle energies, E 2 , and, from these, the corresponding centerof-mass frame (CMF) energies, E * 2 = E 2 2 − P 2 , up to corrections suppressed by factors of exp(−M π L) and exp(−M K L). As in the lattice simulations, we work in a cubic spatial box with side length L, which restricts the allowed total momenta to the set P = (2π/L)d, where d ∈ Z 3 .
The three-particle quantization condition determines the energies of three-particle states in a finite volume, and is given by Here E is the lab-frame energy of the three-particle state, E * = E 2 − P 2 is the corresponding CMF energy, and K df,3 is a three-particle K matrix discussed further below. Although, superficially, the three-particle quantization condition may look similar to that for two particles, the former hides significant complexity. In particular, we note that whereas the quantity F that appears in the two-particle quantization condition is a purely kinematic function, F 3 depends on F , on an additional kinematic function G, as well as on the two-particle K matrix, K 2 . The detailed form is given in Ref. [54]. We also note that while the three-particle quantization in several different contexts takes exactly the same form as in eq. (3.2) (see, e.g. Ref. [80] for the case of degenerate but nonidentical scalars), the indices over which the determinant is taken differ. We now describe the meaning of the determinants in eq. (3.1) and eq. (3.2). For the two-particle quantization condition, F and K 2 are matrices in which the indices { , m} label the angular momentum of the two-particle state, and the determinant runs over these two labels. For the three-particle quantization condition, additional indices are required. These are given by the spectator momentum k, which is constrained to lie in the finite-volume set, and the angular-momentum indices for the pair, { , m}. In addition, we need an index, i, to specify whether the spectator is identical to one of the other two particles (i = 1) or whether it is distinct (i = 2). For such a system, we therefore have a determinant which runs over the indices {k i , , m, i}, where k i labels the momentum of a spectator of particle species i. The quantities F 3 and K df,3 are matrices with these indices.
The derivation of Refs. [14,33] requires a cutoff on the spectator momenta, which must be implemented by a smooth function (rather than a sharp cutoff). One way to understand the need for the cutoff is that, for a given total energy E, the pair is driven below threshold as the spectator momentum increases. Far enough below threshold the two-particle interaction has a left-hand cut, and this must be avoided as it introduces power-law volume dependence that is not controlled in the derivation. The details of the required cutoff are discussed in Ref. [54]. It implies that only a finite number of values of k i contribute. Combined with the cutoff on discussed above, this implies that all matrices in eq. (3.2) are of finite dimension. 5 We now return to the K matrices appearing in the quantization conditions. Both K 2 and K df,3 are infinite-volume Lorentz-invariant quantities, but have important differences. In particular, K 2 is algebraically related to the two-particle scattering amplitude M 2 , while K df,3 is related to the three-particle scattering amplitude M 3 through integral equations [15]. Furthermore, K df,3 depends on the cutoff function, whereas K 2 is cutoff independent above threshold. However, both K 2 and K df,3 share the property of being real and smooth functions of Lorentz invariants below the relevant inelastic thresholds. 6 We make use of these properties to parametrize the K matrices in section 3.2.
For given choices of K 2 and K df,3 , the quantization conditions predict the two-and three-particle spectra. Details of our numerical methods are given in Ref. [54], and a basic python implementation is publicly available [81]. We first block-diagonalize the quantization conditions by projecting them onto irreducible representations of the subgroup of the cubic group that leaves the overall momentum, P , invariant. This is the little group LG(P ). Within each block, we track the smallest eigenvalues of the matrix lying within the determinant, and find zero crossings using a root-finding algorithm. A finite-volume energy level in the given irrep is predicted to occur for each such crossing. We find that it is useful to have a good initial guess for the energy levels, and, for the most part, this is provided by either the measured energy levels or the free levels. We have implemented this methodology in three independent python codes, and all results presented below have been obtained with at least two of these codes. We find that, when running on about ∼ 50 cores, the convergence for the most challenging cases (simultaneous fits to 2 two-particle and 1 three-particle channels) takes 2 − 3 days. The use of the Python compiler numba [82] to speed up core routines is essential to achieve this speed.
As discussed in Refs. [21,24,26], one must ensure that the zero crossings are physical. In particular, the eigenvalue must cross from negative to positive values as the energy is increased, and no higher-order zeros are allowed. The exception to the latter restriction is that there can be higher-order zeros at noninteracting energies, but these are present only because of the truncation of the K matrices, as discussed extensively in Ref. [24]. Such solutions can be dropped. We find that all the crossings are physical. 5 We note that in the NREFT approach a hard cutoff can be used, and its value is not constrained by the position of the left-hand cut [17,18,34]. 6 The only exceptions are that K2 and K df,3 can have poles in the presence of two-and three-particle resonances, respectively. These exceptions do not occur in this work.

Parametrizing K matrices
We now turn to the parametrizations that we use for the matrices K 2 and K df,3 . The former is diagonal in angular momenta, We have written the expression so that it holds for both degenerate (ππ and KK) and nondegenerate (πK) channels: in the former case, η = 1/2, while η = 1 in the latter. In both cases, q is the the magnitude of the three-momentum for each of the two particles in the CMF of the pair. The function H(q 2 ) plays the role of a cutoff. It equals unity for q 2 ≥ 0 [so that the 1 − H term in eq. (3.4) vanishes above threshold], and smoothly drops to zero well below threshold (so that K ( ) . The form of this function depends on the particle masses, as explained in Ref. [54]. The H function does not play a role in the two-particle quantization condition, eq. (3.1), because it also appears in the quantity F , in such a way that the H dependence cancels. It is, however, essential in the three-particle quantization condition (in which K 2 enters through F 3 ), as it cuts off the sum over the spectator momenta.
As in Ref. [49], we explore two choices for the parametrization of the phase shift: an "Adler zero" form, and the effective-range expansion (ERE). The former is motivated by chiral perturbation theory, which predicts that the scattering amplitude vanishes below threshold at the position denoted the Adler zero. For = 0, the Adler zero form is where z 2 and the B n are dimensionless parameters. Two masses appear in this parameterization: M 1 , the mass which we use to set the units, e.g. of q, and M 2 , the mass of the other particle. M 1 and M 2 are chosen from the set {M π , M K } based on the particular scattering process being considered, and may or may not be distinct. At leading order in ChPT z 2 = 1, while B 0 and B 1 take nonzero values to be discussed in the next section, with all other parameters vanishing. In practice, we use two choices of parametrization: 1. ADLER2, in which z 2 = 1, B 0 and B 1 are free parameters, and B n = 0 for n ≥ 2; 2. ADLER3, in which z 2 , B 0 , and B 1 are free parameters, and B n = 0 for n ≥ 2.
We use superscripts and subscripts to denote the corresponding channel, e.g. z 2 πK and B KK 0 . The relation of these parameters to the scattering length a 0 and effective range r 0 is given by Here we are using the convention in which a 0 is positive for repulsive interactions. The ERE form is simply an expansion in powers of q 2 , where the B n are dimensionless. We use the same names for the coefficients as in the ADLER forms, but it will always be clear from the context which fits are being used. The mass M sets the scale of q and may be chosen to be either of the masses of the particles in the scattering pair. In Ref. [49], we found that the Adler zero form was preferred for the ππ and πππ channels, while for kaons the ERE form was slightly preferred. This was not unexpected, as ChPT should work better for pions than kaons. Here we only use the ERE form for KK and KKK channels, where we compare it to the ADLER forms. Specifically, we use 3. ERE3, in which B 0 , B 1 , and B 2 are free parameters, and B n = 0 for n ≥ 3.
The relation of the ERE3 parameters to the scattering length and effective range is For = 1, which is present only for πK scattering, we use a one-parameter form, (3.10) Here P πK 0 is the p-wave scattering length, and the same notation for masses is used as above. The factor of E * 2 is adopted from standard continuum analyses [83,84]. We find that the signal for nonzero p-wave scattering is sufficiently weak that we cannot include higher-order terms and obtain stable fits.
We now turn to the three-particle K matrix. As noted above, K df,3 is a real, analytic function of Lorentz invariants. We use the analog of the ERE, which is an expansion about threshold, and has been worked out for 2+1 systems in Ref. [33]. We denote M 1 as the mass of the particle that appears twice, while M 2 is the mass of the singleton. Implementing the relevant particle interchange symmetries, as well as parity and time reversal, the resulting form is Here, K 0 , K 1 , K B , and K E are real, dimensionless constants, 7 to be determined from fits to the three-particle spectrum, while ∆, ∆ S 2 andt 22 are functions of the Mandelstam variables. The initial (incoming) momenta are {k 1 , k 1 , k 2 }, while the (outgoing) final momenta are {p 1 , p 1 , p 2 }, where k 2 and p 2 are the momenta of the singleton. The kinematic quantities appearing in eq. (3.11) are ∆, ∆ S 2 andt 22 are all of the same order in the threshold expansion, but only ∆ S 2 andt 22 have to a nontrivial angular dependence and lead to contributions with both = 0 and 1.
The expansion in eq. (3.11) can be continued to higher order, with terms of O(∆ 2 ) including d-wave ( = 2) contributions. However, this leads to a proliferation of parameters, and so we restrict ourselves here to max = 1. This is in contrast to Ref. [49], where the analysis used max = 2 for the 3π + and 3K + spectra, which was possible because there are many fewer allowed forms in K df,3 for identical particles. Here, where needed, we have redone the identical-particle fits with max = 1, implying that K df,3 contains only the K 0 and K 1 contributions, since the K B and K E terms in eq. (3.11) vanish for a fully symmetric system.
The explicit forms for the K B and K E terms in the {k m} basis have been worked out in Ref. [54], and we simply use the results.

Results from chiral perturbation theory
Analyzing the two-and three-particle spectra for two different pion masses allows us to investigate the mass dependence of the scattering parameters we derive from fits. In this subsection we collect results from ChPT that will be useful when we fit the dependence of these scattering parameters versus M 2 π /F 2 π . Since our focus is on systems including both pions and kaons, we consider results from SU(3) ChPT. We note that a criterion for the utility of these results is that M 2 K /(4πF K ) 2 1, where F K is the kaon decay constant in the convention that F π 92MeV. For ensembles D200 and N203 the ratio M 2 K /(4πF K ) 2 is 0.11 and 0.13 respectively.
The next-to-leading order (NLO) ChPT expressions for the ππ and KK scattering lengths can be found in Refs. [85,86]. Both depend on the mass of the η meson as well as the pion and kaon, but it is consistent in the NLO contribution to use the LO result 3M 2 η = 4M 2 K − M 2 π , as well as to treat F 2 π and F 2 K as interchangeable. In this way we can express the pion and kaon scattering lengths as functions of just M π , M K , F π , and F K : Here the chiral logarithms are given by while L ππ is an LEC, which is evaluated at the scale 4πF π . We also need the expression for the s-wave πK scattering length, which may be found in Ref. [86]: is the reduced mass, L 5 is another LEC, and (3.21) For the effective ranges, we compare only to the LO prediction, since NLO results from SU(3) ChPT are not given explicitly in the literature (see Refs. [87,88] for the full πK scattering amplitude at NLO). For identical particles, the LO predictions are We turn now to the p-wave scattering length, which is nonzero only for πK scattering. The NLO ChPT prediction for the πK scattering amplitude has been worked out in Refs. [87,88]. However, no closed form for the corresponding scattering length has been provided in the literature, and so we have worked it out and present the results in appendix C-see eq. (C.6). The expression, which is rather lengthy, depends on two different combinations of LECs.
Since we have only two data points, we opt to fit to the expected leading order chiral behavior. From eq. (C.6), one can determine that the p-wave scattering length is finite in the chiral limit (proportional to M K /F 4 π ), so that with higher-order terms suppressed by powers of M π /F π . Neglecting chiral logs, this is equivalent to considering only the effect of the LECs proportional to M 2 K in eq. (C.6). We conclude this section by presenting the LO chiral predictions for K df,3 for 2 + 1 systems, which were derived in Ref. [54]. At LO, only K 0 and K 1 in the expansion shown in eq. (3.11) are nonzero: To obtain these forms we have used the interchangeability of F π and F K in LO terms. We note that K df,3 is independent of the cutoff function at LO in ChPT, and is in this sense a physical quantity at this order. Cutoff dependence only enters at NLO. Predictions are not yet available for K B and K E , since these quantities are expected to appear first at NLO, and no calculation at this order has been done. Following Ref. [54], we choose generic forms for these quantities We stress that the full dependence on x π and x K will likely be much more complicated, but this form satisfies the correct chiral power counting and is sufficient given that we have only two values of pion masses. A potentially important, and so far unquantified, source of systematic errors in our results comes from working at a single lattice spacing. In this regard, it is useful to know the form of the prediction for discretization effects that is given by Wilson ChPT (WChPT), i.e. ChPT including contributions proportional to powers of the lattice spacing a [57,58]. In appendix B we have worked out the leading, O(a 2 ), terms for the two-particle scattering amplitudes and K df,3 , both for the degenerate and nondegenerate cases. The results for scattering lengths are where w 6 and w 8 are dimensionless LECs proportional to a 2 . Note that we make a πK 0 dimensionless by multiplying by the average mass M πK = (M π + M K )/2, rather than the reduced mass µ πK used in eq. (3.17). With this choice we see that all three quantities have the same offset.
-18 - The corresponding results for K df,3 are Thus the predicted shifts are proportional to the same combination of LECs as for the scattering lengths. In section 5.3, we will use these results to estimate the magnitude of discretization effects contributing to K df,3 .
We close with a comment on the appropriate power counting in WChPT. In the standard power counting one takes a 2 Λ 2 QCD ∼ M 2 π /(4πF ) 2 , and, using this, it would be inconsistent to include NLO terms proportional M 4 π /F 4 , while not including discretization contributions proportional to a 3 , a 2 M 2 π , etc. Since such terms have not been calculated, however, we simply attempt fits using the available information, namely using eqs. (3.29) to (3.33), in which we take the continuum NLO expressions given above for the a = 0 part. In effect, we are assuming that a 2 Λ 2 QCD ∼ M 4 π /(4πF ) 4 . We find, in section 5.3, that the a 2 contributions are in fact considerably smaller than the estimate from standard power counting, providing a posteriori justification for this approach.

Extraction of scattering parameters
In this section we discuss the extraction of scattering quantities from the energy levels obtained in lattice QCD. We start by discussing the different strategies that one can use to fit the spectrum using the quantization conditions. We then turn to the results of the fits using different methods. Finally, we compare the different approaches and discuss what seems to be the optimal one for this dataset.

Fitting strategies
In this work, we will use variations of the so-called spectrum method [90]. The main idea is to obtain the best fit parameters by minimizing a χ 2 function that depends only on some spectral quantity X: where X i is the lattice QCD result for that spectral quantity in the i-th energy level, X QC i ( p) represents the prediction from the finite-volume formalism assuming a certain parametrization of the K matrices with parameters p, and C is the covariance matrix of all the X i quantities, such that C ij = cov(X i , X j ).
Several choices for the spectral quantity X are possible. While, in the limit of infinite statistics, all should lead to the same answer, in practice some may be more advantageous. Specifically, some choices can lead to covariance matrices with larger condition numbers, such that the calculation of the inverse matrix appearing in eq. (4.1) may be more unstable. Some examples for X are listed below.
1. The original energy levels obtained from lattice QCD, which are, in general, in a moving (or "laboratory") frame. These are denoted E lab .
2. The energy levels boosted to the center-of-mass frame (CMF), assuming the continuum dispersion relation, i.e. ignoring possible lattice artifacts. The resulting energies are denoted E cm . 8 This is the choice of Ref. [49].
3. The shift with respect to the non-interacting finite-volume energy, with the latter calculated assuming a continuum dispersion relation. This can be done either in the lab or cm frame, yielding ∆E lab or ∆E cm , respectively.
4. In the two-particle sector, it is also possible to use the CM momentum, q 2 . This is the choice used, for example, in Ref. [64].
As we discuss in detail below, we use ∆E lab for our preferred fits. Once the best fit parameters have been obtained by minimizing eq. (4.1), the errors of the best fit parameters (and their covariance) need to be estimated. One possibility to do so is to perform a separate fit on each jackknife sample, and use the results to estimate the covariance of the parameters. While an option in the two-particle sector, present evaluations of the predictions of the three-particle quantization condition are too slow for this approach to be practical in general. Instead, we perform a fit only on the mean, while using the jackknife samples to estimate the covariance matrix C of the data, and then apply the derivative method. This method, discussed in Ref. [49], estimates the covariance between the fit parameters p n and p m as where the derivatives are evaluated numerically at the minimum of the χ 2 function, χ 2 min . A very similar approach is to find the 1σ interval by finding the contour such that χ 2 = χ 2 min + 1, and assuming that the dependence of χ 2 on p is quadratic. In practice, we find that these two methods give essentially identical results, and we use them interchangeably to quote errors in the following.
Another issue to address is how to combine information from the different channels, e.g. those with different numbers of particles. In particular, when analyzing three-particle energies, the two-particle interaction parameters are also needed. In previous work, e.g. in Refs. [42,49], the approach has been to perform a combined fit to both two-and threeparticle energy levels. For example, in the pion sector, one defines a χ 2 function that combines ππ and πππ levels, and uses the covariance matrix of all levels, including crosscorrelations between two-and three-particle energies. Fits using this approach will be referred as "fully correlated fits". In this work, we take this approach one step further by considering nondegenerate systems, where two different two-particle spectra are needed.
For example, for the ππK sector, we must consider also the ππ and πK levels. This leads to a larger number of total energy levels, and one may be concerned about the reliability of the calculation of the covariance matrix between levels. We find, in practice, that this is not an issue in the fits done here (either for the ππK or KKπ cases), but it certainly will become a problem eventually, as the number of channels and levels increases further. For example, we have not attempted fully correlated simultaneous fits to the ππ, πK, KK, ππK and KKπ levels.
Because the issue of fitting to multiple channels is a generic one in multiparticle physics, it is worthwhile investigating alternative approaches. We consider two approaches in which varying amounts of information about the correlations between two-and three-particle levels is dropped. We can imagine these being relevant in situations where one has limited information on correlations, because, for example, different ensembles have been used to calculate two-and three-particle quantities, or the determination of the full covariance matrix is unstable. The underlying idea here is that two-particle spectra might be able to pin down two-particle scattering quantities sufficiently well that the three-particle spectra can be used primarily to determine three-particle scattering quantities.
Our first alternative is to use what we refer to as a "chained fit". Here, we first perform a fit to the two-particle sector. Minimizing the two-body χ 2 function will lead to the twoparticle best fit parameters, p fit 2P and their covariance. We will use these values to construct the chained χ 2 function as: where p = ( p 2P , p 3P ) is a vector that contains the two and three-particle parameters, ∆ p 2P = p fit 2P − p 2P , and X 3P represents a specific spectral quantity for all three-particle energy levels. Note that the method requires also covariance between the three-particle energy levels and the two-particle best-fit parameters, cov( X 3P , p fit 2P ). These can be estimated using a resampling technique, e.g., jackknife.
This approach can be further simplified by neglecting completely the off-diagonal terms in the covariance matrix of eq. (4.3). In this case, the χ 2 functions becomes: This can be seen as the augmented χ 2 of a "Bayesian fit", where p fit 2P is the prior of those parameters. In the context of nondegenerate spectra, a further simplification of such Bayesian fits is possible. For example, for the ππK case, the Bayesian augmentation can include, or not, the correlations between the fit parameters in the ππ and πK channels.

Fit results
In this section we present results from fitting the spectra using the two-and three-particle quantization conditions. We begin with examples showing the impact of using the different fitting methods described above, and then present our core new results for the parameters -21 -describing the ππK and KKπ systems. Finally, we compare the results for the two-particle scattering parameters (scattering length and effective range) obtained using different fits.
For every fit we need to choose a maximum value of E cm for the levels to be included. The quantization conditions that we use formally break down above the lowest inelastic threshold, which for systems of pions occurs when two additional pions can be created (single-pion production being forbidden by G-parity), while for systems involving kaons the breakdown occurs when only a single additional pion can be produced. The issue is discussed in detail in Ref. [49], where it is noted that, in practice, the quantization conditions are likely to remain applicable some distance above the nominal maximal E cm , a conclusion supported by the numerical results of that work. Thus, here we also work with values of E cm that lie above the nominal maximal, making the same choices for the 2π, 3π, 2K, and 3K channels as in Ref. [49], and similar choices for the πK, ππK and KKπ channels.
As discussed above, we include only s-and p-wave terms in the quantization conditions, but not d-wave terms. For two identical particles, i.e. for ππ and KK, this implies that we must exclude from the fits levels that lie in nontrivial irreps, since such levels are only shifted by d-wave terms (p waves being absent). For three identical particles, levels in nontrivial irreps are shifted by s-wave terms because there can be relative angular momentum between the dimer pair and the spectator, but these shifts are incomplete due to the absence of dwave terms in the dimer. Thus we also keep only πππ and KKK levels in trivial irreps in the fits. This is different from the fits in Ref. [49], where were able to include d-wave terms, and thus also levels in nontrivial irreps.
By contrast, for the nondegenerate channels that are of central interest here, the inclusion of p waves implies that fits to levels in all available irreps are possible, and we include such levels in the fits.

Comparison of fitting strategies
We begin by showing examples of the differences between results obtained by fitting to E cm and ∆E lab . As discussed above, the latter fits have the advantage of being to quantities that are closer to those that are actually obtained from the lattice simulations, and are in this sense preferable. We used E cm fits in Ref. [49], where we studied the 3π and 3K systems, and our aim here is to study the impact of changing to ∆E lab fits. We expect that the latter fits will lead, in general, to larger values of χ 2 , but that this provides a more accurate reflection of the goodness of fit. We recall that our fits to correlator ratios yield ∆E lab a. This "primary" quantity is then converted to ∆E lab /M , where M is either M π or M K depending on the quantity being studied, and we use the value of M from the rest frame fit in the corresponding jackknife sample. The "∆E lab fits" are to ∆E lab /M . When fitting to E cm /M , a further conversion is needed, first from ∆E lab to E lab , and then by a boost to the rest frame. The key point is that this conversion depends on M L, and that fluctuations in this quantity between jackknife samples leads to an increase in the errors in E cm /M . Thus we expect, in general, that fitting to E cm /M will lead to a smaller χ 2 , but that this reduction is not due to having a better fit, but rather due to the errors being -22 -overestimated. What we do not know a priori is how the results for the fit parameters will change, and that is the focus of our investigation here.
In a few cases, the difference between E cm and ∆E lab fits is minimal. An example is provided by fits to the πK spectrum on the D200 ensemble, which are shown in table 3. We find good fits in both cases, leading to consistent fit parameters, although there is a small decrease in the errors in the fit parameters when fitting to ∆E lab . Most striking is the change in the condition number of the correlation matrix for the values of E cm or ∆E lab . The correlation matrix is closely related to the covariance matrix, differing by a normalization procedure that guarantees each element of the principal diagonal is 1, and each off-diagonal element lies in the range [−1, 1]: (4.5) By providing a more natural normalization, the correlation matrix is better suited to estimate the condition number, as it avoids issues of significantly different errors among the extracted energies which can lead to artificially large condition numbers. The reduction in the condition number of the correlation matrix by an order of magnitude for the ∆E lab fits implies that the fits will be more stable. We observe a substantial reduction in the condition number in fits to all quantities.  Table 3. Comparison of fitting approaches for the πK spectrum on the D200 ensemble. All quantities are in units in which M π = 1. The fits are to the 26 levels lying below the cutoff E cm = 5M π = (M π + M K ) + 1.62M π , and use the fit form eq. (3.5). The position of the Adler zero is fixed to its leading order value (i.e. z πK = 1).
In most cases, however, the value of χ 2 increases substantially when fitting to ∆E lab . As a first illustration of this behavior, we show, in table 4, fits to the KK spectrum on D200. First, we note that values of χ 2 ref = χ 2 /DOF are large in all cases, where DOF stands for degrees of freedom. One reason for this is that we are fitting without including the d-wave interaction, a choice we make, as explained above, in order to have practical fits when we consider nondegenerate three-particle systems. As an example of the impact of this omission, we note that, were we to include a d-wave scattering length in the E cm fit, χ 2 would be reduced by about 20 [49]. Our focus here, however, is on the differences between E cm and ∆E lab fits, and we see that χ 2 increases significantly, consistent with our general expectation discussed above. However, we note that the fit parameters themselves change little (B KK 1 is reduced by ∼ 2σ), and the errors are essentially unchanged.  -1700(330) -1400(340) -1100(320) Table 5. Fits to the ππ + πππ spectrum on ensemble N203. All quantities are in units in which M π = 1. We use cutoffs of E cm = 3.46M π and 4.46M π respectively from the ππ and πππ spectra, and fit only to levels in trivial irreps, leading to 27 levels for each channel. In the two-particle channel, the fit model used is the Adler zero form given in eq. (3.5), with the Adler zero fixed to its leading order position (z ππ = 1). The fit model used in the three-particle channel includes only the K 0 and K 1 terms of eq. (3.11).
Both tables show the same pattern as for the KK results above: there is an increase in χ 2 when using ∆E lab fits, while fit parameters and errors are largely consistent. The largest change is for B 1 , which decreases by more than 2σ. The conclusions drawn in Ref. [49], namely that K df,3 is significantly different from zero, and that K iso,1 df,3 is negative, remain valid for the ∆E lab fits. Since the fit parameters are highly correlated, one cannot judge the significance of a nonzero K df,3 from the tables alone; using the full covariance matrices we find this to be 7.0σ and 4.4σ for the E cm and ∆E lab fits to ππ + πππ, while the corresponding results for the KK + KKK fits are 6.4σ and 4.5σ, respectively. Thus the significance of the nonzero K df,3 is somewhat reduced, but remains high.
In the remainder of this section we only consider fits to ∆E lab . Our first task is to compare the results of fits using the different choices for χ 2 described in section 4.1 above. -10000(3500) -8300(3600) Table 6. Fits to the KK + KKK spectrum on ensemble D200. All quantities are in units in which M K = 1. We use cutoffs of E cm = 2.53M K = 2M K + 1.26M π and 3.53M K = 3M K + 1.26M π for the KK and KKK spectra, respectively, leading to 28 and 26 levels in the trivial irreps. In the two-particle channel, the fit model used is the Adler zero form given in eq. (3.5), with the Adler zero fixed to its leading order position (z KK = 1). The fit model used in the three-particle channel includes only the K 0 and K 1 terms of eq. (3.11).
We do so only for the ππ + πK + ππK fits, on ensemble D200, as the pattern we find is the same in all other fits. In table 7 we compare the results when using (a) the standard choice of χ 2 , given in eq. (4.1), (b) chained fits, using the χ 2 from eq. (4.3), and (c) Bayesian fits using the χ 2 given in eq. (4.4). (The final column will be discussed in section 4.2.2 below.) We recall that there are two types of Bayesian fits, depending on whether one keeps the correlations between the results of the parameters for the ππ and πK fits, or not. We find little difference in the results of these two approaches, and present results only for the latter.
The standard fit is to a total of 69 levels, 9 using 9 parameters. This is a challenging fit, but we see no signs of numerical instability in the calculation of χ 2 . Indeed, the main challenge is finding the minimal χ 2 with a large number of parameters, and our minimizer takes ∼ 10 3 iterations to converge. In all cases we have checked the fits by repeating them with different initial conditions, and by using three independent codes. The final χ 2 is high, but a large part of this arises from the fit to the ππ sector, where a fit to the 22 levels alone leads to χ 2 ≈ 52 due to the absence of d-wave terms (as discussed above in the context of KK fits). In addition, we have seen above that moving from E cm fits to ∆E lab fits leads to increased χ 2 .
The large number of levels that must be fit motivates investigating the alternatives provided by using chained and Bayesian fits. We stress that the values of χ 2 in the three fits cannot be compared. A rough comparison can be obtained by adding to the χ 2 for fits (b) and (c) the values of the χ 2 from the individual ππ and πK fits, which are 52 and 31, respectively, leading to total values of ∼ 115 for fits (b) and (c). However, this ignores the impact on the χ 2 in fit (a) of including the full correlations between the levels.
Comparing the results, we see that the central values for all two-particle scattering 9 We note here that one level (the second level in the A2(8) irrep) shows ∼ 3σ fluctuations due to the choice of GEVP parameters. We opt to keep it in the fit, but have checked that removing this level from the fit barely impacts the best fit parameters in  Table 7. Fits to the ππ + πK + ππK spectrum on ensemble D200, using ∆E lab . All quantities are in units in which M π = 1. For the first three columns, we use cutoffs E cm = 3.74M π , 5M π = (M π + M K ) + 1.62M π , and 5.4M π = (2M π + M K ) + 1.02M π for the ππ, πK and ππK channels, respectively, leading to 22, 26 and 21 levels. For fit (d), the cutoff for the πK spectra is reduced to 4.64M π = (M π + M K ) + 1.26M π , so that there are 16 πK levels. In the two-particle channel, the fits use the functions of eq. (3.5) and eq. (3.10), with the ππ and πK Adler zeros fixed to their lowest-order values. The fit model used in the three-particle channel is given by eq. (3.11). In fits (b) and (c), the number of DOF is given by the number of levels (21) plus the number of fit parameters in the ππ + πK fits (2+3), minus the number of parameters fit (9).
parameters are very similar, and have essentially the same errors, in all three fits. By contrast, while the values of the three-particle parameters are consistent within errors, those errors are significantly larger for the chained fit than the standard fit, and larger still for the Bayesian fit. In other words, the information on correlations between levels that is lost when using "sequential" fits makes it harder to pin down the (already challenging) three-particle parameters. It is thus no surprise that the significance of the nonzero K df,3 is reduced when moving from the standard fit to the chained and Bayesian fits: it is 3.4σ, 2.4σ and 2.8σ, respectively for the three fits.
A similar pattern is observed for all other channels, and thus we conclude that standard fits are clearly preferable if they are possible, as is the case here. We use only standard fits for our central results to be presented shortly.
The final issue that we address in this subsection is whether to rebin the results on the D200 ensemble. As discussed above, the results on the D200 ensemble are rebinned by N rebin = 3, leading to 771 jackknife samples. As discussed in Ref. [49], this rebinning is useful to account for autocorrelations, and indeed the errors in the energy levels increase as one increases the rebinning factor. However, for N203, where we have only 666 configurations, rebinning can lead to unstable fits when the number of samples becomes too close to the number of levels being considered. We have tested this in several cases, two examples being given in table 5 above and table 8 below. What we find is that the fits are stable for -26 -both N rebin = 2 or 3, that they lead to essentially the same results for all fit parameters, with no change in errors, but that χ 2 increases significantly (as do the condition numbers). We interpret these results as indicating that any autocorrelations have minimal effect on the scattering parameters, while the reduction in the number of samples leads to less well conditioned fits. Thus we choose to use no rebinning for our central fits N203. The choice of using N rebin = 3 for D200 is, therefore, conservative.

Fits to nondegenerate channels
As noted above, our central values are obtained from fits to ∆E lab , using the standard, fully-correlated χ 2 , and without rebinning on ensemble N203. Results for ππ + πK + ππK are shown in table 7 (above) and   Table 8. ∆E lab fits to the ππ + πK + ππK spectrum on ensemble N203. All quantities are in units in which M π = 1. In all fits, the cutoffs for ππ and πK are E cm = 3.464M π = 2M π + 1.464M π and 3.4M π = (M π + M K ) + 1.12M π , corresponding to 27 and 19 levels, respectively. For the first column, the cutoff for ππK is 4.1M π = (2M π +M K )+0.82M π , while for the remaining two columns it is 4.3M π = (2M π + M K ) + 1.02M π , corresponding to 26 and 36 levels, respectively. In the twoparticle channel, the fits are to eq. (3.5) and eq. (3.10) with ππ and πK Adler zeros fixed to their lowest-order values. The fit model used in the three-particle channel is given by eq. (3.11).
For the ππ +πK +ππK fits we show examples of the results of fitting to different levels. Table 7 compares our "standard" fit, in the first column, with an alternative, in the final column, in which the cutoff on the πK levels is reduced from 0.62M π above the inelastic threshold (2M π + M K ) to 0.26M π above this threshold. The fits are of similar quality and give consistent results, and errors, for all parameters, indicating that our more aggressive cutoff is acceptable.   table 7 by powers of (M π /M K ), and similarly for other parameters. Cutoffs are given by E cm = 2.53M K , 1.95M K = (M K + M π ) + 1.26M π , and 2.832M K = (2M K + M π ) + 0.98M π for the KK, πK and KKπ channels, respectively, corresponding to 28, 16, and 29 levels, respectively. In the KK channel, the ADLER2 and ADLER3 fits use eq. (3.5), with the Adler zero fixed to its leading-order position for the former fit and allowed to vary for the latter, while the ERE3 fit uses eq. (3.8). In all three fits the πK channel is fit using the ADLER2 form for the s-wave K matrix and eq. (3.10) for the p wave, while the fit model used in the three-particle channel is given by eq. (3.11). from the two fits are consistent within errors, and we use the results from the 82-level fit henceforth.
For the KK + πK + KKπ fits we show comparisons between different fit choices. 10 In Ref. [49] it was observed that the KK interactions were slightly better described on the D200 ensemble (for which the kaon is heavier) if, compared to our standard two-parameter Adler-zero fit, one either allowed the position of the Adler zero to float, or used a threeparameter ERE fit. Thus we have done all three fits for both ensembles and compare the results in the tables. Freeing the position of the Adler zero leads to no significant improvement in χ 2 ref on D200, but a mild improvement on N203. Switching to an ERE fit, leads to a slight improvement on D200, but a worse fit on N203. Results for the πK and KKπ fit parameters are essentially unchanged.
One of our major aims is to study how well three-particle interactions can be determined by using a large collection of energy levels. As can be seen from the tables, the significance of the nonzero values for individual terms in K df,3 varies, with the most significant parameter being K 1 on the D200 ensemble. The significance of the entire K df,3 10 Note that in the KK + πK + KKπ case we fit to ∆E lab /MK rather than to ∆E lab /Mπ, as in the ππ + πK + ππK fits. This makes a small difference as there are fluctuations of the single-particle masses between jackknife samples.  in table 8 by powers of (M π /M K ), and similarly for other parameters. Cutoffs are given by E cm = 2.9M K , 2.66M K = (M K + M π ) + 1.12M π , and 3.521M K = (2M K + M π ) + 0.95M π for the KK, πK and KKπ channels, respectively, corresponding to 23, 19, and 32 levels, respectively. For the KK channel, only the ground state in the frame with d 2 = 9 is kept. In the KK channel, the ADLER2 and ADLER3 fits use eq. (3.5), with the Adler zero fixed to its leading-order position for the former fit and allowed to vary for the latter, while the ERE3 fit uses eq. (3.8). In all three fits the πK channel is fit using the ADLER2 form for the s-wave K matrix and eq. (3.10) for the p wave, while the fit model used in the three-particle channel is given by eq. (3.11).
being nonzero, including the correlations between the parameters, is 3.4σ and 3.6σ for the ππ + πK + ππK standard fits on the D200 and N203 ensembles, respectively, and 5.0σ and 2.4σ for these ensembles in the KK + πK + KKπ fits.

Visualization of fits
We have shown several global views of the fits in figures 2 and 11 to 13, which illustrate that the fits match the spectrum well both in the fit range and also above the nominal inelastic threshold. To investigate this more carefully, we need to zoom in and show results for the quantities ∆E lab to which we actually fit. This is the purpose of this section.
We focus first on the KK + πK + KKπ fit on the D200 ensemble, as this is the fit for which both K df,3 and P 0 are determined to be nonzero with the greatest significance. Specifically, we consider the ADLER3 fit to 73 levels in table 9. In the upper panel of figure 3 we compare the results for ∆E lab from our simulations to those obtained using the quantization condition with the best fit parameters. The latter is shown by the red dots just above each of the data points. We see that the shifts are small (of order 1% of M K ), but are determined with errors ranging from a few percent to 10 − 15%. We note that all levels, whether in trivial or nontrivial irreps, are shifted by the same order of magnitude. This differs from the situation with the two particle spectrum, as will be seen below, and is due to the fact that the two-particle K matrix, K 2 , contributes to all irreps. Indeed, as is well known, and was clearly seen in Ref. [49], the dominant physical effect leading to these energy shifts is the two-particle interaction, and it is nontrivial to determine the subdominant contribution of K df,3 . To illustrate the impact of K df,3 , we also show, as blue squares, the results predicted by the quantization condition if we set K df,3 = 0 while keeping all other parameters unchanged. The shifts in the levels are small, reaching the size of the error bars in the data only for the higher energy levels. From this figure alone, it would appear that there is little chance of determining K df,3 , but this is misleading because all levels are correlated, and the fit includes not only these 29 levels, but also the 28 K + K + and 16 K + π + levels. One illustration of the claim that the nonzero value of K df,3 leads to a significant improvement in the fit is that χ 2 increases from 162 to 188 when K df,3 is turned off.
We also include, using orange triangles, the result of turning off the p-wave πK scattering amplitude by setting P 0 = 0, with all other parameters unchanged from the ADLER3 fit. This changes the levels by an amount that is typically smaller than that caused by setting K df,3 to zero.
We next investigate how the fit works for the levels that are not included in the fit, because they lie above the maximal CMF energy, E cm = (2M K + M π ) + 0.98M π . These levels thus lie at or above the inelastic threshold. We have determined eight of them, and the comparison of their lab shifts to the ADLER3 fit is shown in the lower panel of figure 3, along with the predictions if K df,3 or P 0 are set to zero. The values of the shifts are comparable to those for the fitted levels, and we see that the fit continues to work at a similar level of accuracy even in the inelastic regime.
We next display the corresponding results for the K + π + channel in the same fit. These are shown in figure 4: the upper panel shows 16 levels that are included the fit, and the lower panel shows 10 that are not. We recall from table 9 that the cutoff for the fit lies at E cm = (M K + M π ) + 1.26M π in this channel, and thus lies slightly above the inelastic threshold. Since two-particle levels do not depend on K df,3 , we show only the impact of setting P 0 to zero.
For the trivial irreps, the situation is similar to that for the KKπ levels: the fit works well, and continues to do so in the inelastic regime. The shifts due to the nonzero P 0 are small, although they increase for the higher-lying levels. The levels in nontrivial irreps, however, are shifted only by the p-wave interactions (as can be seen by the fact that the yellow triangles lie at ∆E lab = 0. This is as expected based on group-theoretical considerations. For some of these levels, there is evidence that the energy shift is negative, and it is this that leads to the result that the scattering length is slightly attractive.
The situation for the other fits and ensembles is qualitatively similar, and we do not show the corresponding plots for these other cases.
(not included in fit) -31 - (not included in fit)

Derived results for two-particle scattering quantities
In this section we collect the results for the scattering lengths and effective ranges for ππ, πK and KK scattering. These quantities are obtained from the fits presented above, and from additional fits to the following systems: ππ alone, ππ + πππ, KK alone, and KK + KKK. We stress that these additional fits are different from those presented in Ref. [49], as here we fit to ∆E lab , rather than E cm . Our overall aim in this section is to compare the results, and in particular the errors, obtained by using fits to different sets of levels. We discuss the chiral behavior of these results in the following section. One of the methods used here is to determine the scattering length from the 1/L expansion of the energy shift of the ground state (i.e. for which the noninteracting state has all particles at rest). For two nondegenerate scalars the result is [91] ∆E (2) g.s. = 2πa 0 where µ 12 is the reduced mass of the pair, a 0 the corresponding scattering length, and the constants are [11] For three degenerate particles of mass M the result is [11] We observe that the first two orders are a factor of 3 larger for three particles than for two particles, reflecting the number of two-particle pairs. This implies an increased sensitivity to a 0 in the three-particle channel. The above-described truncated 1/L expansions of the ground state energy shifts are sometimes used to determine scattering lengths from two-particle energy shifts (see, e.g., Ref. [92]). Truncation at the 1/L 5 term is needed to have a one-to-one relation between the energy shift and a 0 , since the 1/L 6 terms include the effective range and, for three particles, also a subtracted version of the three-particle amplitude at threshold [93]. Thus, in this approach, one must proceed by assuming the 1/L 6 term is numerically small, and then estimate the resulting systematic error due to truncation [92]. Here our interest is less in the central value obtained in this fashion, but rather in the size of the error obtained in this method compared to those from global fits.
We begin with results for ππ scattering, which are shown in table 11 and table 12, respectively, for the D200 and N203 ensembles. These include the results of the core fits presented above, as well as those to different numbers of ππ levels, and the results of fitting the truncated 1/L expansions to the ground-state energy shifts. We see that all central values are consistent within 1 − 2σ, including results obtained from in Ref. [49] using E cm fits. Our focus here is on a comparison of the errors. In particular, we want to know if anything is gained by increasing the number of levels in the fits, i.e. moving from a single level (in the ground-state-only fits), to the ππ levels alone, and finally to the fits involving   [49]). The ππ + πK + ππK fits are from table 7 (59 level fit). The next two rows give the results of fits to the ππ data alone, with cutoff for the two fits being E cm = 3.74M π and 3.0M π , respectively. The next block of two rows give the results obtained using the ground-state energy shifts alone, as described in the text. The final row gives the results from Ref. [49] from E cm fits, with the first error being statistical and the second a systematic error due to the variation between fits.
Fit   Table 5, while those for ππK fits are from Table 8. The next two rows give results from fits to the ππ data alone, using cutoffs E cm = 3.436M π (which is the same as that used in the ππ + πππ and ππ + πK + ππK fits) and 3.0M π . The next block of two rows give the results obtained using the ground-state energy shifts alone, as described in the text. The final row gives the results from Ref. [49] from E cm fits, with the first error being statistical and the second a systematic error due to the variation between fits.
two-and three-particle levels. Our results indicate that, for the scattering lengths, the errors decrease slightly as the number of levels in the fit increases. This trend is what we would have naively expected, but the size of the change is relatively small. We attribute -34 -this smallness to the fact that the higher levels constrain the phase shift at values away from threshold, and because of the strong correlations between two-and three-particle levels. A similar pattern is observed for the errors in a 0 r 0 . In light of this discussion, in the chiral fits below we will take the values from the fits involving two-and three-particle levels. Specifically, we use the ππ + πK + ππK fits, since the χ 2 ref of this fit is significantly smaller on the N203 ensemble, while the values and errors are essentially the same as those from the ππ + πππ fits on D200.
We next turn to the results for KK scattering parameters, which are collected in table 13 and table 14. Here we also show results with the ADLER3 and ERE3 choices for the KK phase shift, as these can lead to improved fits, as noted above. For a given choice of fit type, the sizes of the errors follow a similar pattern to those for ππ scattering, with those from fits involving two and three particles being smallest. As expected, errors increase for the fits in which the KK phase shift is described with three parameters.
All fits using ∆E lab lead to consistent results for a KK 0 . Additionally, these fits are consistent with those using E cm at the 2σ level. There is, however, much more variation in the results for a KK 0 r KK 0 . Finally, in table 15 and table 16, we present results for πK scattering parameters. In both cases we show results with two choices of E cm cutoff for the πK fits: the 16-level fit with excellent χ 2 ref and the more aggressive 32-level fit with slightly poorer fit quality. The results for P πK 0 , a πK 0 and r πK 0 are consistent within 1 − 2σ. The pattern of the sizes of errors is consistent with that described above.   in the ERE3 block. The KK and KKK fits use cutoffs of 2.9E K and 3.9M K , respectively, keeping only levels in trivial irreps, and discarding all but the lowest d 2 = 9 level in the KK channel. This leads to 23 levels in both channels.  Table 15. Comparison of πK scattering parameters from different fits on the D200 ensemble, using units in which M π = 1. All fits use the ADLER2 + ERE1 form for the πK phase shift. The first row gives the results from the 59-level fit in table 7, while the second gives the result from the  ADLER3-KK fit in table 9, with the latter converted to units in which M π = 1. The next two give the results of fits to the πK data alone, The cutoff energies for the these fits are, respectively, E cm = 3.64M π and 5.0M π , with the former also used for the πK channel in the ππ + πK + ππK and KK + πK + KKπ fits. The final row shows the result of fitting the ground-state energy shift to the 1/L expansion.   fit in table 10, with the latter converted to units in which M π = 1. The next two rows give the results of fits to the πK data alone, using cutoff energies 3.4M π and 4.3M π , respectively, with the former also used for the πK channel in the ππ + πK + ππK and KK + πK + KKπ fits. The final row shows the result of fitting the ground-state energy shift to the 1/L expansion.
- 38 -In this section, we present and discuss our final results for different two-and three-meson scattering quantities. First, in section 5.1, we provide our final numbers for these quantities by combining results from different fits. Then, in section 5.2, we compare these numbers to expectations and predictions from ChPT. Finally, in section 5.3, we discuss the size of discretization errors, based on the LO Wilson-ChPT results presented above.

Final results for scattering parameters
In section 4.2, we have presented results for the two-and three-meson scattering parameters using different fit forms and strategies. While we find overall consistency, it is useful to have a set of final results that contain both statistical uncertainties as well an estimate of the systematic spread due to different fit forms. Note that we always use lab-frame shift fits for this set of final results.
We first consider two-particle parameters. Results for scattering lengths and for the combination M X a XY 0 r XY 0 are given in table 17 and table 18, respectively. In each case, the central values are obtained by averaging the results from a pair of three-particle fits (which pair will be explained shortly), and include a systematic error that is obtained from the spread of the fit results (and which we call the "fit systematic"). For the ππ case, the central values are the average of those from fits to ππ + πππ and ππ + πK + ππK-see tables 11 and 12-while the fit systematic is the standard deviation obtained from the two results. For the πK case, the same procedure is used but now taking the ππ + πK + ππK and KK + πK + KKπ fit results from tables 15 and 16. For the KK case, the central value and statistical error are obtained using the same procedure but taking the KK + KKK and KK + πK + KKπ fits obtained using the ADLER3 parametrization for the KK phase shift from tables 13 and 14. The fit systematic is given by the standard deviation of the results from these two fits together with those obtained with the ADLER2 and ERE3 parametrizations (six fits in total). In all cases, we take the largest of the statistical errors when combining results.  For the three-particle parameters, we quote the final values in table 19. In this case we take the results from a single fit-the one we view as most reliable, based on the discussion in the previous section-and do not quote a fit systematic as any such error would be dwarfed by the statistical errors.

Comparison with chiral perturbation theory
We now turn to the comparison of the final results of tables 17 to 19 to ChPT. We perform fits to ChPT expressions wherever they are available, and compare with the form of the expected dependence on M 2 π in other cases. We also compare to previous results in the literature.
We start with the scattering lengths. Figure 5 shows the results for each of the dimensionless s-wave scattering lengths, 11 along with the LO chiral prediction, and a fit to the NLO expressions in eqs. (3.14), (3.15) and (3.17). We have not determined the correlations between the three quantities on a given ensemble, and thus use an uncorrelated fit. In order to plot the result as a function of (M π /F π ) 2 , we use an interpolating function for M K /F K as a function of M π /F π extracted from the results of Ref. [49]. We find a good fit with χ 2 /DOF = 1.5/4 (6 data points and 2 parameters). We determine the two LECs (evaluated at a renormalization scale µ = 4πF π ) to be L ππ = −8.77(36) × 10 −4 , L 5 = 0.0(1.5) × 10 −3 . (5.1) We can compare these values to previous determinations. In Ref. [49], a larger value of L ππ = −1.13(3) · 10 −3 was found, which is many standard deviations away from our new result. However, an important difference that may explain the discrepancy is that Ref. [49] included d-wave interactions, which is not possible here. For L 5 , we can compare to results from lattice QCD, which are summarized in the FLAG report [94], and based on Refs. [95,96]), or from phenomenological determinations [97]. These are, however, quoted at a different renormalization scale, µ = 770 MeV. Changing the renormalization scale with Γ 5 = 3/8, we find that the result from our fit yields L 5 (770 MeV) = 1.0(1.5) · 10 −3 .
Varying the choice of 4πF π to take for the initial scale (using the physical value of F π , or the value on either of the ensembles) leads to changes in L 5 that are significantly smaller than the error. Our result for L 5 is in agreement with all values in the literature, although we note that our error is much larger than that in the other values. Next, we discuss our results for the effective range parameters, which are presented in table 18 in the combination M 2 X r XY a XY 0 . For the case of identical particles (X = Y = π or K), the LO ChPT prediction from section 3.3 is that this quantity equals 3. For two pions, the results lie 15% and 25% below this prediction on the D200 and N203 ensembles, respectively, which is consistent with being due to an NLO correction. For two kaons, the results lie very far away from the LO prediction. Both findings are qualitatively similar to those obtained in Ref. [49].
For the πK channel, which is a novel result of this work, the LO ChPT prediction-given in eq. Our results in table 18 lie ∼ 25% and ∼ 30%, respectively, below the LO ChPT prediction. Again we view this as reasonable consistency, given the absence of NLO corrections.
-41 - We now turn to the p-wave π + K + scattering length, reported in the rightmost column of table 17 through the dimensionless combination P πK 0 = −M 3 π a πK 1 . Note that, in contrast to all the s wave results, the value of P πK 0 corresponds to slightly attractive interactions. We plot the results for the two ensembles in figure 6, including a fit to the leading chiral behavior given by eq. (3.23), which shows reasonable consistency.
We also plot the NLO ChPT prediction given in appendix C. To do so we use values for the requisite LECs determined in Ref. [89] from experimental data (specifically, fit 10 to O(p 4 ) from that work). As can be seen, the NLO ChPT result has the same sign as our results, but its magnitude is significantly smaller. The failure of NLO ChPT for this quantity was, in fact, expected, based on the observation of Ref. [88] that the NNLO contribution is two orders of magnitude larger than the NLO one at the physical point (see table 2 of that work). π /F 2 π . A fit to the leading chiral scaling of P πK 0 ∝ (M π /F π ) 3 is shown with the corresponding error band, as well as the NLO ChPT prediction as described in appendix C. Also included are the result at the physical point from the dispersive analysis of Ref. [98] (see Table 29 of that work), and the lattice QCD determination of the HadSpec collaboration at a heavier pion mass [56].
We can also compare to the expectations and results in the literature from experiment and dispersive analyses. The current understanding is summarized in figure 10 of Ref. [98]. Experimental results [99] for the p-wave phase shift point to a negative (repulsive) value at high energies. By contrast, the dispersive analysis indicates a change of sign for the phase at around √ s M K + 3M π (physical values of the masses), resulting in an attractive scattering length. The value from analysis of Ref. [98] is also shown in figure 6. As can be seen, our results for the two ensembles of this work are in qualitative agreement with the low-momentum behavior found by the dispersive analysis. We note that our fits only -42 -include levels in the region where the phase shift is expected to stay positive.
We are aware of two other LQCD results concerning p-wave π + K + scattering. First, Ref. [100], reports a single energy level far from threshold (at much higher energy than our levels, and in the inelastic regime) that is dominated by p-wave interactions. There, the p-wave πK interactions seems repulsive, which is consistent with what experiments find at those high energies. This result therefore gives no information concerning the scattering length.
Second, Ref. [56] computed the p-wave scattering length at heavy meson masses, M π 391 MeV and M K 549 MeV, and its sign and magnitude are consistent with our results at lighter pion masses. We include this result with the label "HadSpec" in the plot, although it is not strictly speaking comparable as Ref. [56] does not follow the same chiral trajectory. We conclude that, overall, the results from LQCD are in qualitative agreement with dispersive and experimental results. Finally, we compare our results for K df,3 for 2+1 systems to ChPT. In figures 7 and 8 we plot the results for K 0 and K 1 for ππK and KKπ scattering, respectively. We compare to the LO ChPT predictions of eqs. (3.24) and (3.25), and find substantial disagreement, most notably in the sign of K 1 , while the magnitudes are better matched. Similar disagreement has been observed for 3π and 3K systems [49]. There are two possible interpretations for this disagreement. First, it may be that we have underestimated the errors in the determinations of K 0 and K 1 . One possibility is that discretization errors might be large, although we present evidence against this option in section 5.3. Second, NLO terms in ChPT may be substantial, and invalidate the LO result, such as in the case of M 2 K a KK 0 r KK 0 discussed -43 -above. To address the latter possibility, a NLO ChPT calculation would be needed, but, while NLO results are available for the three-particle scattering amplitude [101,102], the relation to K df,3 has yet to be worked out. LO ChPT In figures 9 and 10 we plot the results for K B and K E for ππK and KKπ scattering, respectively. These quantities vanish at LO in ChPT; their first nontrivial contribution is expected to appear at NLO in ChPT. Since a NLO calculation has yet to be done, we have fit to the expected chiral scaling given in eqs. We find a reasonable description of the data based on these fit forms.

Discretization errors
Up to this point we have neglected the effects of discretization errors in our two-and threeparticle fits. Since the ensembles used in this work are O(a) improved, these errors are of O(a 2 ). Here we extend the fits by including the leading a 2 terms predicted by WChPT. As explained in section 3.3, this is only consistent with chiral power counting if we assume a 2 Λ 2 QCD ∼ M 4 π /(4πF π ) 4 . We begin with the two-particle scattering lengths. The WChPT results of eqs. (3.29) to (3.31) predict that each of these quantities receive a common offset proportional to a 2 . Repeating the global fit to the six s-wave scattering lengths shown in table 17, allows us -44 - π K E (ππK) Figure 9. Results for K B and K E for ππK scattering as a function of M 2 π /F 2 π . Fits to the expected leading chiral behavior given in eq. (3.27) are plotted alongside the data. For better visibility, the x-coordinates of the left-most datapoints have been shifted.  Figure 10. Results for K B and K E for KKπ scattering as a function of M 2 π /F 2 π . Fits to the expected leading chiral behavior given in eq. (3.28) are plotted alongside the data.
to find the value of this offset, which we denote as follows, The results for L ππ and L 5 are consistent with those found in the fit without discretization errors, given in eq. (5.1), and the offset is found to be consistent with zero. The result for δ a (M a 0 ) can be converted into a determination of the O(a 2 ) LECs, 2w 6 + w 8 = −0.14(23) . (5.7) As noted above, the magnitude of this effect is very small. For example, at the physical pion mass M π a ππ 0 ∼ 0.04 in our fits, which is about an order of magnitude larger than the central value of, or error in, δ a (M a 0 ).
Using the WChPT results given in eqs. (3.32) and (3.33), we can use eq. (5.7) to predict the size of the O(a 2 ) terms in the leading isotropic contribution to K df,3 . This is because they are proportional to the same combination of LECs. The resulting predictions are presented in table 20. We observe that the size of the discretization errors is an order of magnitude smaller than the central values and errors we obtain from fits to the spectrum This indicates that, given the precision with which we are able to calculate the parameters in K df,3 , discretization effects can be viewed as negligible. Similar results hold also for the 3π and 3K channels, for which the WChPT results are presented in appendix B.  Table 20. Final results for K 0 compared to the estimated discretization effects based on WChPT, δ a (K 0 ) = −6(2w 6 + w 8 )M 2 X /F 2 X , where X = π for the ππ + πK + ππK fits, and X = K for the KK + πK + KKπ fits.

Conclusion
This works represents the first determination, theoretical or experimental, of quantities related to ππK and KKπ scattering at maximal isospin. In particular, we have studied these systems using LQCD for two different ensembles along a trajectory of approximately constant trace of the quark mass matrix, 2m + m s const.
We have extracted more than 200 finite-volume energies across the ππK and KKπ systems, leading to roughly 35 -45 and 20 -30 levels below the inelastic thresholds on N203 and D200, respectively. As with our previous work in Ref. [49], the ability to determine such a large number of energies was enabled by advanced algorithms, like stochastic LapH and -46 -common subexpression elimination. In order to facilitate the use of our extracted spectrum by other collaborations, we have made the jackknife samplings for all energies determined in this work available in HDF5 format as ancillary files with the arXiv submission.
As explained in section 4.1, we have tested several strategies to extract scattering quantities from the finite-volume spectra. Our findings can be summarized as follows. First, in the systems of this work, it is better to perform fits to energy shifts with respect to the noninteracting energies, rather than to the energy levels themselves. This leads to less correlation between energy levels (better-conditioned correlation matrices), and to more constrained best-fit parameters. Second, it seems that the best approach when dealing with multiple processes is to fit all the available information at the same time (which here means a simultaneous fit to the two-and three-hadron spectra). Other fitting strategies, such as first fitting the two-particle spectra and then using the results as Bayesian priors for a fit to the three-particle levels, produce consistent results, but with significantly larger uncertainties. We expect that these conclusions will be relevant for other finite-volume systems.
The main results of this work are summarized in tables 17 to 19. First, we are able to extract two-particle threshold parameters, including the p-wave π + K + scattering length. Second, for each of the 2+1 systems, we extract the four parameters of K df,3 corresponding to first and second order in the threshold expansion. We observe the expected hierarchy, that is, interactions in the KKπ channel are stronger than in ππK. Moreover, the values for the different terms of the K matrix have the correct order of magnitude, at least based on chiral expectations. We are also able to determine their values to be nonzero with greater than 2σ significance in some cases, and we find the entire K df,3 to be nonzero with significance between 2.4 and 5.0σ in all channels on each ensemble.
Another novel feature of this work is the extraction of discretization errors in K df,3 , which at leading order in the appropriate power counting enter only in K 0 . By comparing the two-particle scattering lengths computed on the lattice to those calculated in WChPT, we are able to determine a numerical value for the magnitude of discretization effects in a WChPT calculation of K 0 . We find the discretization errors to be small compared to both the central values and uncertainties of K 0 for each scattering channel and on each ensemble.
Natural extensions of this work are to three pseudo-Goldstone bosons at nonmaximal isospin, for which two-and three-meson resonances are present, or to the doubly charmed tetraquark T cc . While there is a long way to go to determine properties of particles coupling to both two-and three-particle channels, and involving particles with spin, such as the Roper resonance, the strategies and tools of this work bring us one step closer.

A Energy levels used in fits
In this appendix we give further details of the fits described in the main text. First, in table 21 we list the energy levels used in the fits. Second, in figures 11 to 13 we show the CMF energy levels for the N203 KKπ, D200 ππK and D200 KKπ spectra, respectively, along with the predictions from the quantization condition using our standard fit parameters. The corresponding result for the N203 ππK spectrum is shown in the main text in figure 2.  Table 21. Energy levels used in the fits of this work. Notation is as follows: "A 1u + E u " means the lowest level in the A 1u irrep, and the lowest in the E u irrep. When two different energy ranges have been chosen, the levels in brackets contribute only to the fit with the higher energy cutoff. For instance, in D200 a fit with 16 πK levels is shown in the last column of table 7 ("(d) Standard *"), whereas another with 26 levels is shown in the first column of the same table ("(a) Standard"). In this case, "A 1g + T 1u + [A 1g ]" means that the lowest A 1g and T 1u levels contribute to the fit with 16 levels, and the second A 1g level is additionally included in the fit with 26 levels. The energy levels included in ππ, KK, πππ and KKK systems can be read off from 3.8

B ChPT result for K df,3 including O(a 2 ) errors
In this appendix we extend previous SU (3) ChPT results for two-and three-particle K matrices by including discretization errors. We work at LO, and use SU(3), rather than SU(2), ChPT due to the presence of both kaons and pions. The methodology for including the effects of nonzero a into ChPT was developed in Refs. [57,58], and is generally referred to as Wilson ChPT, or WChPT for short. Previous work has considered ππ scattering in SU(2) WChPT [103] (working at NLO), but, to our knowledge, no previous calculations of two-or three-particle scattering using SU (3) WChPT have been performed. For present-day lattice calculations, with light quark masses close to their physical values, and with nonperturbative O(a) improvement, the appropriate power counting to use is m q ∼ p 2 ∼ a 2 (leaving factors of Λ QCD implicit). In other words, discretization errors, which begin at O(a 2 ), are comparable to the squared masses of the psuedo-Goldstone bosons (PGBs), M 2 PGB = O(m q ). For an example of the appropriateness of this power counting, we point to Ref. [104], where it is noted, in the related context of simulations with twisted-mass fermions, that an O(a 2 ) contribution to PGB mass-squareds (w 8 f 2 in the notation below) is indeed of the same order as the physical M 2 π . The LO chiral Lagrangian (in Euclidean space), including discretization terms, is given by [58] where Σ ∈ SU(3) is the field that contains the PGBs, χ = 2B 0 M is proportional to M , the renormalized mass matrix, and F ≈ F π 92 MeV and B 0 are the usual continuum LO LECs. Additionally, we haveâ = 2W 0 a, where W 0 , W 6 , and W 8 are LECs associated with discretization errors. We assume exact isospin symmetry, so that M = diag(m , m , m s ). With the Lagrangian in hand, we can now compute the LO pion and kaon masses, where s is the usual two-particle Mandelstam variable. We thus find that, when we express the results in terms of the PGB masses, only a single linear combination of w 6 and w 8 appears in all three expressions. Indeed, if one carries out the calculation in SU(N ) WChPT, the same result is obtained, with all the N dependence absorbed into the PGB masses [3 becoming N in eq. (B.2) and eq. (B. 3)]. Thus we can compare the result for M ππ with that obtained in SU(2) WChPT in Ref. [58], and find complete agreement.
The O(a 2 ) terms in eq. (B.5), eq. (B.6), and eq. (B.7) again display a violation of chiral symmetry, for the scattering amplitudes do not vanish at threshold in the chiral limit. To see this explicitly, we give the result for the scattering lengths These results agree with the LO parts of eq. (3.14), eq. (3.17), and eq. (3.15), in the limit that a 2 → 0. By taking appropriate linear combinations of the scattering lengths, one can cancel the O(a 2 ) term. This prediction could be tested in simulations with several lattice spacings. We note that there are no O(a 2 ) corrections to the effective range at LO. Finally, we turn to K df,3 . This has previously been calculated without discretization effects for three pions in Ref. [42], three kaons in Ref. [49], and for ππK and KKπ scattering in Ref. [54]. We refer to those works for the methodology, and simply quote the final results. We find that [see eq. (3.12) for notation]  Table 22. π + K + operators used in this work. For each set of momenta that are equivalent up to allowed rotations, one representative choice of d ref is given, where the total momentum is P = (2π/L)d ref . The momentum squared (in units of (2π/L) 2 ) of the individual single particles in each operator is listed as d 2 sh , where "sh" indicates the particular single hadron. The lab-frame free energies in units of M π are listed for those operators that are used on each ensemble together with irreps that are included. Daggered and starred irreps were only used for the N203 and D200 ensembles, respectively.