Estimating the two-particle $K$-matrix for multiple partial waves and decay channels from finite-volume energies

An implementation of estimating the two-to-two $K$-matrix from finite-volume energies based on the L\"uscher formalism and involving a Hermitian matrix known as the"box matrix"is described. The method includes higher partial waves and multiple decay channels. Two fitting procedures for estimating the $K$-matrix parameters, which properly incorporate all statistical covariances, are discussed. Formulas and software for handling total spins up to $S=2$ and orbital angular momenta up to $L=6$ are obtained for total momenta in several directions. First tests involving $\rho$-meson decay to two pions include the $L=3$ and $L=5$ partial waves, and the contributions from these higher waves are found to be negligible in the elastic energy range.


Introduction
Currently, the best method of computing the masses and other properties of hadrons from quantum chromodynamics (QCD) involves estimating the QCD path integrals using Markov-chain Monte Carlo methods, which requires formulating the theory on a spacetime lattice. Such calculations are necessarily carried out in finite volume. However, most of the excited hadrons we seek to study are unstable resonances. Fortunately, it is possible to deduce the masses and widths of resonances from the spectrum determined in finite volume.
The idea that finite-volume energies can be related to infinite-volume scattering processes is actually rather old, dating back to Refs. [1,2] in the mid-1950s. First suggestions of applying such techniques for gauge field theories appeared in Ref. [3]. In Ref. [4] in 1986, Lüscher studied the volume dependence of the energy spectrum of stable particle states in massive quantum field theories, then examined the volume dependence of scattering states in Ref. [5] soon thereafter. In Ref. [6] in 1991, Lüscher then found relationships between finite-volume energies and infinite-volume scattering phase shifts in the case of two identical spinless particles having zero total momentum and interacting via a central potential. The advantages of using sectors with non-zero total momenta were then described by Rummukainen and Gottlieb in Ref. [7] in 1995. These calculations were later revisited in 2005 using a field theoretic approach by Kim, Sachrajda, and Sharpe in Ref. [8]. This work focused on the case of a single channel of identical spinless particles, but the total momentum could be any value allowed by the boundary conditions. The results of Ref. [8] were eventually generalized in Refs. [9][10][11][12][13][14], among others, to treat multi-channels with different particle masses and nonzero spins.
As new techniques, such as those introduced in Refs. [15] and [16], have made possible accurate determinations of the energies of states created by two-hadron operators, an increasing number of lattice QCD studies of scattering phase shifts have appeared .
The purpose of this paper is to provide explicit formulas, software, and new fitting implementations for carrying out two-particle scattering studies using energies obtained from lattice QCD. We first express the quantization condition that relates lattice QCD energies to the scattering matrix in terms of the more convenient K-matrix and a Hermitian matrix B which we refer to as the "box matrix" since it essentially describes how the spherical partial waves manage to fit into a cubic finite volume. We restrict our attention to periodic boundary conditions. We obtain explicit results for several spins and center-of-momentum orbital angular momenta up to L = 6 with total momentum of the form P = (0, 0, 0), (0, 0, p), (0, p, p), (p, p, p) for p > 0, many of which have never before appeared in the literature. More importantly, the software for evaluating all of these box matrix elements is made available. For multiple partial waves and channels, the quantization condition at an allowed energy is insufficient to determine the entire Kmatrix at that energy. Approximating the K-matrix elements using physically-motivated functions of the energy involving a handful of parameters is needed, with the hopes that these parameters can be estimated by appropriate fits using a sufficiently large number of different energies. We discuss two fitting strategies with such a goal in mind. In our first tests of the fitting methods, we incorporate the L = 3 and L = 5 partial waves in the decay of the ρ-meson to two pions and find their contributions to be negligible in the elastic energy region.

Quantization condition
We begin by summarizing the quantization condition that relates each finite-volume energy to the scattering S-matrix. We work in an L 3 spatial volume with periodic boundary conditions. For a total momentum P = (2π/L)d, where d is a vector of integers, we determine the total energy E in the lab frame for a particular two-particle interacting state in our lattice QCD simulations. We boost to the center-of-momentum frame and calculate 2 Let N d denote the number of two-particle channels that are open, and denote the masses and spins of the two scattering particles in channel a by m ja and s ja , respectively, for j = 1, 2. If |E 2 cm | ≥ |m 2 1a − m 2 2a |, then we can calculate the following quantities in each channel: The total energy E is then related to the dimensionless unitary scattering S-matrix through the quantization condition [6][7][8]14]: We use an orthonormal basis of states, each labelled by |Jm J LSa , where J is the total angular momentum of the two particles in the center-of-momentum frame, m J is the projection of the total angular momentum onto the z-axis, L is the orbital angular momentum of the two particles in the center-of-momentum frame (not to be confused with the lattice length here), S in the basis vector is the total spin of the two particles (not the scattering matrix). The angular momentum, orbital angular momentum, and spin operators are defined as usual in terms of the Pauli-Lubanski tensor W µ = − 1 J m J L S a |F (P ) |Jm J LSa = δ a a δ S S 1 2 where j 1 m 1 j 2 m 2 |JM are the familiar Clebsch-Gordan coefficients, and the W (P a) matrix elements are given by Z lm (s a , γ, u 2 a ) π 3/2 γu l+1 a (2L + 1)(2l + 1) (2L + 1) The Rummukainen-Gottlieb-Lüscher (RGL) shifted zeta functions Z lm , introduced in Refs. [6,7], are evaluated using [9,13] Z lm (s, γ, with Y lm (x) = |x| l Y lm ( x) and We choose Λ ≈ 1 for fast convergence of the summation, and the integral in Eq. (7) is efficiently done using Gauss-Legendre quadrature. The function F 0 (x) is given in terms of the Dawson or erf function: Many of these RGL shifted zeta functions are related to one another. Some of these relations and the symmetry properties of these functions are summarized in Appendix A. The above relations apply for both distinguishable and indistinguishable particles, since the associated symmetry factors cancel in the quantization condition in general [12]. The only difference that occurs with indistinguishable particles is that certain combinations of L and S cannot occur. In the absence of isospin, L + S must be even for identical particles. For identical particles of isospin I 1 , L + S + I − 2I 1 must be even, where I is the total isospin.

The K-matrix and box matrix
The quantization condition in Eq. (4) is a single relation between an energy E determined in finite-volume and the entire S-matrix. In the very limited case of a single channel with a single partial wave, this relationship can be used to directly extract the scattering phase shift at energy E. When multiple partial waves or multiple channels are involved, this single relation is clearly not sufficient to extract all of the S-matrix elements at energy E. The best way to proceed is to approximate the S-matrix elements using physically motivated functions of the energy E involving a handful of parameters. Values of these parameters can then hopefully be estimated by appropriate fits using a sufficiently large number of different energies. Such a procedure has long been standard in the partial wave analyses of particle experiments.
The S-matrix in Eq. (4) is dimensionless and unitary. Since it is easier to parametrize a real symmetric matrix than a unitary matrix, one usually employs the real and symmetric K-matrix. The K-matrix was first introduced by Wigner [49] and Wigner and Eisenbud [50] studying resonances in nuclear reactions. Its first use in particle physics was an analysis of resonance production in Kp scattering by Dalitz and Tuan [51]. A comprehensive review of the K-matrix formalism is given in Ref. [52], and a more recent review appears in Ref. [53]. Recent applications of the K-matrix to study resonances in lattice QCD can be found in Refs. [32,34].
Defining the transition operator T via S = 1 + iT , the K-matrix can be defined by 4 and hence, its relationship to the S-matrix is Using the above definition, it is straightforward to show that the unitarity of S implies the Hermiticity of K. Time reversal invariance of the S-matrix implies that the K-matrix must be symmetric. Given its Hermiticity, this means that K is real and symmetric. Rotational invariance implies that the K-matrix must have the form where a , a denote other defining quantum numbers, such as channel, and K (J) is a real, symmetric matrix that is independent of m J . Invariance under parity also gives us that where η P ja denotes the intrinsic parity of particle j in the channel associated with a. For a single channel of spinless particles, the S-matrix is diagonal with elements S L = e 2iδ L , so the K-matrix diagonal elements are K L = tan δ L , where δ L is the scattering phase shift of the L-th partial wave. For a short-ranged potential in the single-channel case, one can derive the effective range expansion, which states that where q cm = q 2 cm , the constants a L are called scattering lengths (although only the S-wave constant a 0 has the dimensions of a length), and r L are known as the effective ranges. The multichannel generalization [54][55][56] of the effective range expansion is where K −1 L S a ; LSa (E cm ) is a real, symmetric, and analytic function of the center-ofmomentum energy E cm .
The effective range expansion given in Eq. (17) suggests the convenience of writing since K −1 L S a ; LSa (E cm ) is real and symmetric and expected to behave smoothly with energy E cm . It is then straightforward to show that the quantization condition of Eq. (4) can be written where we define the box matrix by 5 This box matrix B (P ) is Hermitian for u 2 a real. Whenever det K = 0, which is usually true in the presence of interactions, the quantization condition can also be written The Hermiticity of B (P ) and the fact that K is real and symmetric for real u 2 a ensures that the determinants in the quantization conditions of Eqs. (19) and (21) are real. Note that K and B (P ) do not commute in general, which means 1 − B (P ) K and 1 − KB (P ) are not Hermitian. However, it is easy to show that their determinants must be real.
Again, rotational invariance of the K-matrix implies that K has the form When S = S = 0, then J = L and J = L yielding Given that K −1 is expected to be analytic in E cm , an obvious parametrization of the inverse of the K-matrix over a small range of energies is using a symmetric matrix of polynomials in E cm : where α, β are compound indices referring to orbital momentum L, total spin S, and channel a, and the c (Jk) αβ form a real symmetric matrix for each k. Another common parametrization (see, for example, Ref. [53]) expresses the K-matrix as a sum of poles with a background described by a symmetric matrix of polynomials: where the couplings g form a real symmetric matrix for each k. These can be written in Lorentz invariant form using E cm = √ s, where the Mandelstam variable s = (p 1 + p 2 ) 2 , with p j being the fourmomentum of particle j.

Block diagonalization
In the previous sections, we expressed the matrices F (P ) and B (P ) in terms of the orthonormal center-of-momentum frame basis states labelled by |Jm J LSa . In this basis, the quantization condition in each of Eqs. (4), (19) and (21) is problematic due to the need to evaluate the determinant of an infinite matrix. If we can transform to a basis in which both B (P ) and K are block diagonal, then we only need to examine the determinant separately in each block. Each block has infinite dimension, but if we truncate in the orbital angular momentum, keeping only states with L ≤ L max , then each truncated block has a finite and reasonably small size. 6 Under a symmetry transformation G which is either an ordinary spatial rotation R or spatial inversion I s , the total momentum P changes to GP , and if we define a unitary matrix Q (G) by where D (J) m m (R) are the familiar Wigner rotation matrices, one can show that the box matrix satisfies The result in Eq. (27) is very important since it allows us to block diagonalize the B (P ) matrix. If G is an element of the little group of P , then GP = P and Gs a = s a , and we have Since Q (G) is unitary, this implies that the B (P ) matrix commutes with the matrix Q (G) for all G in the little group of P . This means that we can simultaneously diagonalize B (P ) and Q (G) . By rotating into a basis formed by the eigenvectors of Q (G) , we can reduce the B (P ) matrix into a block diagonal form since the matrix elements of B (P ) between different eigenvectors of Q (G) must vanish. Rotations, reflections, and spatial inversion do not change J, L, S, a when acting on basis state |Jm J LSa . These symmetry operations only mix states of different m J . A partial diagonalization of B (P ) can be achieved by diagonalizing Q (G) for each J, L, S, a, or equivalently, by projecting onto the irreducible representations of the little group of P . These eigenvectors or projections can be labelled by the irreducible representation (irrep) Λ and irrep row λ of the little group, and an integer n identifying each occurrence of the irrep Λ in the |Jm J LSa reducible representation. In other words, to block diagonalize B (P ) , we need to apply a particular unitary change of basis: where η = (−1) L . From Eq. (26), one sees that the transformation coefficients depend on J and (−1) L , but not on L itself, and are independent of S, a. Our procedure for computing the transformation coefficients is as follows. To simplify notation, we suppress the S, a indices, define η = (−1) L , and abbreviate the state |Jm J LSa by |J η m J . For a given J, L, S, a, we apply the standard group theoretical projections onto the (2J + 1) basis vectors |J η , J , |J η , J − 1 , . . . , |J η , −J for the first row λ = 1 of each irrep of the little group, producing the vectors (one for each m J ) where G denotes the little group, g G is the number of elements in the little group, d Λ is the dimension of the irrep Λ, Γ (Λ) (G) is the unitary matrix representing G in the Λ irrep, and for rotations R and spatial inversion I s , we have If the irrep Λ does not occur in the J representation of SU (2) subduced to G, all of the resulting vectors will be zero. If the irrep Λ occurs once in the subduction, then only one nonzero vector will occur, which can then be suitably normalized. If the Λ irrep occurs more than once, then there is some freedom in choosing the basis vectors. We first look for linear combinations of the projected vectors that are simpler in form with less terms, order the resulting vectors by increasing complexity, then apply a Gram-Schmidt procedure to obtain orthonormal basis states. Once we have basis vectors |Λλn for the first row λ = 1 of all Λ irreps, where positive integer n is the occurrence label, we then obtain the partner basis vectors for the other rows µ using the transfer operation Our choices of irreducible representation matrices Γ (Λ) µν (G) are presented in Ref. [57], and the irrep labels for the various little groups are listed in Ref. [57] as well. Given these choices, we have applied the above procedure using software written in Maple. We have computed basis vectors for all J ≤ 8 for momenta P = (0, 0, 0), (0, 0, p), (0, p, p), (p, p, p), with p > 0. Given the very large number of such basis vectors, it is not possible to present them all here. However, explicit tables containing our basis vectors have been generated and are available from Ref. [58].
Expressing the box matrix in this basis, one can show that B (P ) is diagonal in Λ, λ, but not in the occurrence index n. Given Eq. (20), we find that we can write The box matrix depends on a only through u a and s a . Notice that in Eq. (33) we use the irrep label Λ B instead of Λ to label the matrix elements of B (P ) . We wish to reserve the irrep Λ to describe the symmetry of the block in question for the full system, which includes the intrinsic parities of the constituent particles. The B (P ) matrices transform independently of these intrinsic parities. If η P 1a η P 2a = 1, then Λ B = Λ, but for a channel in which η P 1a η P 2a = −1, Λ and Λ B may not be the same. When the total momentum P = 0, Λ and Λ B have opposite parity labels: g ↔ u in the O h label subscripts. For P = 0 when η P 1a η P 2a = −1, one must carefully deduce the relationship between Λ and Λ B . In each case, one consults the conjugacy classes of the little group and identifies the classes that involve improper transformations (those involving parity); for a given irrep Λ, one takes the character vector and flips the signs of the characters in the classes involving improper symmetry operations, and the resulting character vector identifies Λ B . The relationships of Λ B to Λ when η P 1a η P 2a = −1 for various momenta P are summarized in Table 1.
With these transformation coefficients, the box matrix elements can be obtained in this block diagonal basis by straightforward matrix multiplication. Maple code was written to perform these evaluations. In this software, the symmetries of the RGL zeta functions are used to express the results in terms of the minimum number of independent functions (those listed in Tables A.1 to A.4). To check our results, we explicitly verified that (a) elements of the box matrix between states in different blocks are all zero, (b) the box matrix in each block is Hermitian, and (c) the box matrix blocks are exactly the same for all rows in each irrep. 8 Table 1: Relationship of B-matrix irrep Λ B to full symmetry irrep Λ. When the product of intrinsic parities η P 1a η P 2a = 1, then Λ B = Λ. For a channel in which η P 1a η P 2a = −1, the relationship of Λ B to Λ is shown in the table below. LG denotes "little group", and n > 0 in P = (2π/L)d below. For details about the little group irreps (single-valued and double-valued) listed below, see Ref. [57].

d
LG Expressions in terms of the RGL shifted zeta functions for a small selection of the box matrix elements B (P Λ B Sa) J L n ; JLn (E) which we have determined are presented in Appendix B. We have obtained explicit expressions for all box matrix elements with L ≤ 6, total spin S ≤ 2, and total momentum P = (0, 0, 0), (0, 0, p), as well as all box matrix elements with L ≤ 6, S ≤ 3 2 , and P = (0, p, p), (p, p, p), with p > 0. We have developed and tested software, written in C++, to evaluate these box matrix elements. This software is freely available [58] and is described below.
Lastly, we need to express K in the new basis. Given Eq. (22) and the orthonormality of the states in both the |Jm J LSa basis and the block diagonal |ΛλnJLSa basis, one can show that (34) where η = (−1) L and η = (−1) L . If η = −η , the situation is much more complicated. However, in QCD, we should never need such matrix elements. Recall that invariance under parity and rotations implies that the K-matrix elements are nonzero only when All of the mesons we would ever want to consider in the 2-particle states needed to relate finite-volume energies to the K-matrix have negative parity: π, η, φ, K, D, B, η b , η c , J/ψ, Υ, and so on. All of the baryons we would ever want to consider in 2-particle states have positive parity: N , ∆, Σ, Ξ, Ω, Λ, Λ c , Λ b . Thus, all meson-meson states of interest have η P 1a η P 2a = 1, and all meson-baryon states of interest have η P 1a η P 2a = −1, and in QCD, we should never need 2-to-2 K-matrix elements between two-hadron states in which the products of the intrinsic parities are different. In other field theories, it may occur that The box matrix is diagonal in total spin S and in the compound index a. However, the K-matrix allows mixings between different spins and channels, which means that the block structure of the box matrix is not the same as that of 1 − B (P ) K and K −1 − B (P ) . For a given P , S, a, the box matrix blocks can be labelled by Λ B , λ, but the box matrix elements are independent of the irrep row λ. The K-matrix is diagonal in total angular momentum J and is also independent of the irrep row λ, as long as η = η . Given the independence of both matrices on the irrep row, λ is superfluous as a block identifier. Thus, for a given P , we can label the quantization blocks of 1 − B (P ) K and K −1 − B (P ) in the |ΛλnJLSa basis solely by the irrep label Λ, where Λ is the irrep associated with the K-matrix. 9

Software overview
As the tables in Appendix B show, the formulas for evaluating the box matrix elements are cumbersome, and the evaluations of the RGL shifted zeta functions are complicated and require care. One goal of this work is to make software available [58] that evaluates these quantities for use by practitioners. A brief overview of the software is given here.
To evaluate 1 − B (P ) K or K −1 − B (P ) in a given block for a particular P , the basis states involved in the block must first be determined. To accomplish this, the individual spins s 1a and s 2a of particles 1 and 2 in each channel a are needed, and a maximum orbital angular momentum L (a) max must be specified in each channel of particle species. Given the particle spins, the allowed total spins S (a) in each channel can be determined, and these are then combined with all allowed L values up to L (a) max to produce the allowed J values. The software also checks to see if there are any basis states α such that K αβ are zero for all states β, in which case, the state α is removed from the basis. Once the basis of states is determined, Eqs. (33) and (34) can then be used to evaluate 1 − B (P ) K or K −1 − B (P ) in this basis. Again, the irrep Λ B must be used for the box matrices contained in each Λ quantization block for a given P .
In our software, all masses and lengths are specified in terms of some reference energy m ref . The user must specify the vector of integers d specifying the total momentum, the little group irrep Λ, and the dimensionless quantity m ref L, as well as the number of channels. For an anisotropic lattice of temporal spacing a t and spatial spacing a s , the aspect ratio ξ = a s /a t must also be given. In each channel, the user must give the dimensionless quantities m 1a /m ref , m 2a /m ref the particle spins s 1a and s 2a , the product of intrinsic parities η P 1a η P 2a , and a maximum orbital angular momentum L (a) max . With this information, the software sets up the basis of states. When given a lab-frame energy E/m ref or a center-of-momentum energy E cm /m ref , the software evaluates the entire box matrix and the K-matrix, then returns the determinant of either 1 − B (P ) K or K −1 − B (P ) . To evaluate these quantities for different resamplings, the software can reset the masses and box lengths without rebuilding the basis. The fit forms of Eqs. (24) and (25) are used in the software, but these can be easily modified if other forms are desired.

Fitting
Let κ j , for j = 1, . . . , N K , denote the parameters that appear in the matrix elements of either the K-matrix or its inverse K −1 . Once a set of energies for a variety of twoparticle interacting states is determined, the primary goal is then to determine the bestfit estimates of the κ j parameters using the quantization determinant, as well as to determine the uncertainties in these estimates. In this section, we describe two methods to achieve this.
To set the stage, we first summarize the fitting procedure commonly used in lattice QCD. For an observable O, let E(O) denote a Monte Carlo estimate of the observable obtained using the entire Markov-chain ensemble of gauge configurations, and let E where N r is the number of resamplings and the factor N (r) depends on the resampling scheme. For the jackknife and bootstrap methods, it is given by It often occurs that a set of observables is believed to be reasonably well described by a set of model functions containing unknown parameters. In such cases, the goal is usually to find best fit estimates of these parameters. Arrange the observables into the components of a vector R and the fit parameters into a vector α. Denote the set of model functions by the vector M (α, R) which depend on the parameters and which might depend on the observables themselves. The i-th component of M (α, R) gives the model prediction for the observable corresponding to the i-th component of R. In lattice QCD, we generally determine the best fit estimates of the α parameters as the values which minimize a correlated−χ 2 of residuals given by where the vector of residuals is defined by r(R, α) = R − M (α, R) and σ ij = cov(r i , r j ) is the covariance matrix of the residuals. Since the observables are usually obtained using the same ensemble of gauge field configurations, the residuals in Eq. (38) are not statistically independent so the presence of the covariance matrix in the likelihood function is very important. Usually, the minimization of χ 2 with respect to the parameters α is accomplished using computer software, such as Minuit2 [59]. If the model estimates of any of the observables depend on any of the other observables, then the covariance matrix must be recomputed using Eq. (35) and inverted each time the parameters α are changed during the minimization process, making for a rather laborious minimization. A significant simplification occurs if the model estimates are completely independent of the observables. In this case, cov(r i , r j ) = cov(R i , R j ), which needs to be computed and inverted only once at the start of the minimization. The statistical uncertainties in the best-fit parameter estimates can usually be obtained from the minimization software, which typically provides the covariances of the parameter estimates under certain assumptions. An alternative approach to determining these covariances is to perform the minimization of for each resampling k and obtain the covariance of the fit parameter estimates cov(α k , α l ) using Eq. (35). 11 Best fit estimates of the K-matrix parameters can be improved by utilizing results from multiple Markov-chain ensembles and lattices. One approach to performing such fits is to minimize the χ 2 of Eq. (38) taking the elements of σ ij to be zero between the estimates from different ensembles, then obtain the covariances of the best fit parameter estimates from the minimization software. An alternative approach is to ensure that N r is the same for all ensembles, then use the resamplings of all ensembles in the χ 2 of Eq. (38) with the covariance matrix estimated using Eq. (35). Given the statistical independence of the different ensembles, Eq. (35) naturally yields covariances between observable estimates from different ensembles which are very nearly zero.
Again, the primary goal is to determine the best-fit estimates of the κ j parameters appearing in K or K −1 from the quantization condition, as well as to determine the uncertainties in these estimates. Having made the above introductory comments, we now describe two methods to achieve this.

Spectrum method
For each P and irrep Λ, one obtains as many lab frame two-interacting-particle energies E k as possible, staying below the thresholds for three or more particles. For the observations R i in the fit, an obvious choice would be to include the lab-frame energies E k or the center-of-momentum energies E cm,k . Here, we choose the E cm,k energies. The quantization condition with the chosen functional forms of K or K −1 then provides the model predictions of the observations. This involves scanning the quantization determinant in E cm to find the values that result in the determinant having zero value. Evaluating the determinant requires evaluating the box matrix elements, which requires knowing s a , u 2 a for each channel a. To determine s a , u 2 a , one needs to know the masses m 1a , m 2a in each decay channel, the spatial lattice volume L 3 , and the lattice aspect ratio ξ = a s /a t if an anisotropic lattice is used. Unfortunately, these quantities must be obtained from the Monte Carlo simulations, and hence, are observations. This poses the problem that the predictions cannot be obtained solely from the parameters of the model, independent of the observations.
A simple way around this problem is to include the masses m 1a , m 2a in each decay channel, the spatial lattice volume L 3 , and the lattice aspect ratio ξ = a s /a t as both observations and model parameters. In addition to the energy observations E In summary, the observations in the χ 2 minimization in this first method are for k = 1, . . . , N E and j = 1, . . . , N p . If there are N p particle species in all of the decay channels and N E energies found, then there are N obs = 2 + N p + N E observations on an 12 anisotropic lattice. Improved results can be obtained by increasing N E by using several different P , Λ blocks. The model parameters are for i = 1, . . . , N K and j = 1, . . . , N p , where N K is the total number of parameters in the K-matrix elements. The total number of fit parameters is N param = 2 + N p + N K . For an isotropic lattice, the anisotropy observable and fit parameter can be omitted. Evaluating the predictions M i (α) of the model for the N obs observations is done as follows. The parameters m  observations are not so easily done. One needs to scan the quantization determinant in E cm to find which values yield a zero value. For a given E cm , one uses the parameters κ j to evaluate the K-matrix or its inverse, and determines the box matrix elements in terms of the RGL zeta functions using the parameters m a . This is a rather onerous task. Computing the determinant for a given E cm is quite complicated, and this must be done many times in order to bracket and then numerically find all of the needed zeros of the determinant using bisection or a Newton-Raphson type algorithm. One must then match each root found with the appropriate observed E cm . Let E (model) cm,k denote each energy root found. In summary, the residuals in this method are We emphasize that, in this method, the E (model) cm,k are very difficult quantities to compute using the model parameters in Eq. (41). This method (without the model-parameter trick and without properly treating all covariances) has been used in Refs. [32][33][34]45].

Determinant residual method
The difficulty in calculating the model predictions in the first method leads us to seek other simpler methods. In this second method, we introduce the quantization determinant itself as a residual. In the determinant, we use the observed box matrix elements, which requires the observed energies and the observed values for the particle masses, lattice size, and anisotropy.
Expressing the quantization condition in terms of a vanishing determinant is just a convenient way of stating that one eigenvalue becomes zero. The determinant itself is not a good quantity to use as an observable since it can become very large in magnitude for larger matrices. Instead of the determinant, we express the quantization condition using the following function of matrix A, having real determinant, and scalar µ = 0: When one of the eigenvalues of A is zero, this function is also zero. This function can be evaluated as a product of terms, one for each eigenvalue of A. For eigenvalues of A 13 which are much smaller in magnitude than |µ|, the associated term in the product tends towards the eigenvalue itself, divided by |µ|. However, the key feature of this function is that for eigenvalues which are much larger than |µ|, the associated term in the product goes to e iθ for real θ. This function replaces the large unimportant eigenvalues with unimodular quantities so that the function does not grow with increasing matrix size. This is a much better behaved function, bounded between -1 and 1 when the determinant is real, which still reproduces the quantization condition. The constant µ can be chosen to optimize ease of numerical root finding or χ 2 minimization. In this method, the model fit parameters are just the κ i parameters, and the residuals are chosen to be cm,k ) could be used in the Ω function. Clearly, the model predictions in this method are dependent on the observations themselves, so the covariance of the residual estimates must be recomputed and inverted by Cholesky decomposition throughout the minimization as the κ j parameters are adjusted. However, this is still much simpler than the root finding required in the spectrum method.
An advantage of this method is that the complicated RGL zeta functions only need to be computed for the box matrix elements as observables; they do not need to be recomputed as model parameters are changed. Since we cannot completely remove the dependence of the model predictions on the observables in this method, there is no advantage in introducing model parameters for the energies, particle masses, and the lattice anisotropy. Hence, we do not need to recompute the box matrix elements as the model parameters are adjusted in the χ 2 minimization. The model predictions involve only the κ j parameters and the observed energies, particle masses, and anisotropy.
Treating the box matrix elements as observables enables a natural interface between the lattice calculation and phenomenology. If the box matrix elements and center-ofmomentum energies are calculated, then together with the covariances, they contain all the information required to extract the scattering amplitudes. Non-lattice practitioners can use them without, for example, implementing the RGL zeta functions. These quantities can act as a bridge between lattice QCD computations and phenomenological applications.

Tests of fitting procedures
As first tests, we applied the determinant residual method to the interacting ππ energies in the I = 1 nonstrange sector in irreps relevant for extracting the P -wave amplitude. The operators used and the energies obtained are described in Ref. [17]. These energies were obtained on a 32 3 × 256 anisotropic lattice with m π ≈ 240 MeV. Defining the fit forms we used are Table 2: First tests of the determinant residual method applied to the interacting ππ energies in the I = 1 nonstrange channel described in Ref. [17]. These energies were obtained on a 32 3 × 256 anisotropic lattice with mπ ≈ 240 MeV. In Ref. [17], the number of energy levels used was N E = 19. The fits below using N E = 20 include an additional energy from a B + 1 d 2 = 1 irrep in which the leading partial wave is L = 3. Fits which used det( K −1 − B) as residuals are indicated by − in the first column. Fits with a value of µ in the first column used Ω(µ, K −1 − B) as the residuals. Fits in which higher partial waves were excluded are indicated by − in the fifth and sixths columns corresponding to the scattering lengths of such waves. (24) Our results are summarized in Table 2. In Ref. [17], the number of energy levels used was N E = 19. In one half of our fits, we used just these N E = 19 levels. In the other half of our fits, we also included an elastic energy from an additional B + 1 d 2 = 1 irrep in which the leading partial wave is L = 3, as was done in Ref. [33]. These fits are presented in the rows in Table 2 in which N E = 20. Fits in which higher partial waves were excluded are indicated by − in the fifth and sixths columns of Table 2, corresponding to the scattering lengths of such waves.
In the first three rows of Table 2, results are shown using det( K −1 − B) as the residuals. The Ω function was not used in these three fits. Using the determinant alone, we were not able to perform the minimization including all L = 1, 3, 5 partial waves. Also, the errors on the fit parameters using the determinant alone and using just the N E = 19 levels were found to be dramatically larger. In all subsequent fits, we utilized Ω(µ, K −1 − B) as the residuals. Using the Ω function, we were able to find the minimum of the χ 2 function much more easily. For µ = 1, we found that the minimum χ 2 /dof values were uncomfortably large. This was remedied by increasing µ to a value around µ = 8 or larger.
The most important thing that the test fits in Table 2 demonstrate is that the effects of higher partial waves can be taken into account using the determinant residual method. Also, our results show that the phase shifts from the L = 3 and L = 5 waves are negligible in this energy range, justifying our neglect of these waves in Ref. [17]. This is consistent with a phenomenological determination of m 7 π a 3 = 5.65(21) × 10 −5 taken from Ref. [60]. The table also shows that the contribution from the L = 3 wave can be reliably estimated without using the additional B + 1 d 2 = 1 irrep. To further illustrate the effect of the higher partial waves, we define a quantity ∆ for each energy where L = 1 contributes by This quantity is shown in Fig. 1 for the µ = 10, N E = 20 fit including the three L = 1, 3, 5 partial waves. One sees that ∆ is consistent with zero throughout this energy range, again justifying our neglect of the L = 3 and L = 5 waves in Ref. [17].
In the future, we plan to utilize both the spectrum and residual determinant methods in the analysis of meson-meson and meson-baryon systems involving multiple channels. Studies involving the K * (892) and a 0 (980) should appear soon. Various baryon resonances are also being investigated.

Conclusion
The purpose of this paper was to provide explicit formulas, software, and new fitting implementations for carrying out two-particle scattering studies in lattice QCD. We introduced a so-called "box matrix" B which describes how the partial waves fit into the cubic finite volumes of lattice QCD simulations. The quantization condition was expressed in terms of this Hermitian matrix B and the real, symmetric scattering Kmatrix. The effective range expansion was used to introduce an analytic matrix K −1 . 16 We obtained explicit expressions for the box matrix elements for several spins and centerof-momentum orbital angular momenta up to L = 6 with total momentum of the form P = (0, 0, 0), (0, 0, p), (0, p, p), (p, p, p) for p > 0. More importantly, the software for evaluating all of these box matrix elements was announced and overviewed. Lastly, we discussed two fitting strategies for estimating the parameters used to approximate the K-matrix. First tests involving ρ-meson decay to two pions included the L = 3 and L = 5 partial waves, and the contributions from these higher waves were found to be negligible in the elastic energy range.
Acknowledgements: This work was supported by the U.S. National Science Foundation under award PHY-1613449. Computing resources were provided by the Extreme Science and Engineering Discovery Environment (XSEDE) under grant number TG-MCA07S017. XSEDE is supported by National Science Foundation grant number ACI-1548562. We acknowledge helpful conversations with Raul Briceno and Steve Sharpe. We dedicate this work to the late Keisuke Jimmy Juge, our longtime collaborator and good friend.

Appendix A. Symmetries of the RGL shifted zeta functions
In this appendix, we summarize the various symmetries of the RGL shifted zeta functions, and some of the relations between them. Using Eq. (7), it is straightforward to show the following relations and symmetries. For any rotation R, one finds When m 1 = m 2 , then the RGL shifted zeta functions must be zero for odd l. When m 1 = m 2 , then the RGL shifted zeta functions are no longer zero for odd values of l.
One last property of the RGL shifted zeta functions which follows from Eq. (7) is For s = (0, 0, 0), the little group is the orthogonal group O h which includes parity. Applying these relations and solving, one finds that the only nonzero values for angular momentum quantum numbers l ≤ 12, 0 ≤ m ≤ l are those listed in Table A.1. When d = (0, 0, 1), the little group is C 4v . Applying all group elements in the little group, and solving the resulting relations, we find that the only nonzero values for l ≤ 12, 0 ≤ m ≤ l are those listed in Table A.2. When d = (0, 1, 1), the little group is C 2v . We have applied all group elements in the little group and solved the resulting relations for l ≤ 12. Results for l ≤ 9, 0 ≤ m ≤ l are listed in Table A.3. When d = (1, 1, 1), the little group is C 3v . Again, we have determined all relations for l ≤ 12, but only results for l ≤ 8, 0 ≤ m ≤ l are presented in Table A.4.

Appendix B. Box matrix elements in the block diagonal basis
Expressions in terms of the RGL shifted zeta functions for a small selection of the box matrix elements B (P Λ B Sa) J L n ; JLn (E) which we have determined are presented in Tables B.1-B.8. These quantities depend on a only through s a and u a . We have obtained explicit expressions for all box matrix elements with L ≤ 6, total spin S ≤ 2, and total momentum P = (0, 0, 0), (0, 0, p), as well as all box matrix elements with L ≤ 6, S ≤ 3 2 , and P = (0, p, p), (p, p, p), with p > 0. These are available from Ref. [58]. The software for evaluating these box matrix elements was described in Sec. 5 and is freely available [58].       Table A.3: Nonzero elements of Z lm (s, γ, u 2 ) with s = (0, s, s), s > 0, for l ≤ 9, 0 ≤ m ≤ l. Negative values of m can be obtained using Z l,−m (s, γ, u 2 ) = (−1) m Z * lm (s, γ, u 2 ). If the two particle masses are equal m 1 = m 2 , then the elements below for odd l are zero. We have also determined these relations for 9 < l ≤ 12, but these are not listed here.
l Nonzero elements Dependent nonzero elements 0 ReZ 00 1 ReZ 10  Table A.4: Nonzero elements of Z lm (s, γ, u 2 ) with s = (s, s, s), s > 0, for l ≤ 8, 0 ≤ m ≤ l. Negative values of m can be obtained using Z l,−m (s, γ, u 2 ) = (−1) m Z * lm (s, γ, u 2 ). If the two particle masses are equal m 1 = m 2 , then the elements below for odd l are zero. We have also determined these relations for 8 < l ≤ 12, but these are not shown here.     J L n ; JLn (E) for various irreps with P = 0 and total spin S = 0.
These quantities depend on a only through sa and ua. R lm is short hand for (γπ 3/2 u l+1 a ) −1 Re Z lm (sa, γ, u 2 a ). The Hermiticity of B can be used to obtain other elements that are not shown.   J L n ; JLn (E) for various irreps with P = 0 and total spin S = 1.
These quantities depend on a only through sa and ua. R lm is short hand for (γπ 3/2 u l+1 a ) −1 Re Z lm (sa, γ, u 2 a ). The Hermiticity of B can be used to obtain other elements that are not shown. J L n ; JLn (E) for P = (2π/L)(0, 0, n) and total spin S = 0. These quantities depend on a only through sa and ua. R lm is short hand for (γπ 3/2 u l+1 a ) −1 Re Z lm (sa, γ, u 2 a ). The Hermiticity of B can be used to obtain other elements that are not shown.

J
L n J L n u  J L n ; JLn (E) for P = (2π/L)(0, 0, n) and total spin S = 1 2 . These quantities depend on a only through sa and ua. R lm is short hand for (γπ 3/2 u l+1 a ) −1 Re Z lm (sa, γ, u 2 a ). The Hermiticity of B can be used to obtain other elements that are not shown.

J
L n J L n u −(L +L+1) a B Λ B = G 2 (partial) J L n ; JLn (E) for P = (2π/L)(0, n, n) and total spin S = 0. These quantities depend on a only through sa and ua. R lm is short hand for (γπ 3/2 u l+1 a ) −1 Re Z lm (sa, γ, u 2 a ), and I lm is used to abbreviate (γπ 3/2 u l+1 a ) −1 Im Z lm (sa, γ, u 2 a ). The Hermiticity of B can be used to obtain other elements that are not shown.

J
L n J L n u