Scattering from finite-volume energies including higher partial waves and multiple decay channels

A new implementation of estimating the two-to-two $K$-matrix from finite-volume energies based on the Luescher formalism is described. The method includes higher partial waves and multiple decay channels, and the fitting procedure properly includes all covariances and statistical uncertainties. The method is also simpler than previously used procedures. Formulas and software for handling total spins up to $S=2$ and orbital angular momenta up to $L=6$ are presented.


Introduction
A key goal in lattice QCD is to determine the spectrum of hadrons. Lattice simulations can only calculate the energies of stationary states in finite volume. Since excited hadrons are unstable resonances, their masses and decay widths must be deduced from the finite-volume energies obtained in lattice QCD using rather complicated formulas, developed over several decades in Refs. [1][2][3][4], among others. This talk reports on work completed in Ref. [5] to provide explicit formulas, software, and new fitting implementations for carrying out two-particle scattering studies using energies obtained from lattice QCD. In our first tests, 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
Let P = (2π/L)d, where d is a vector of integers, denote a total momentum in an L 3 spatial volume with periodic boundary conditions. The center-of-momentum energy E cm is related to the lab frame energy E by Speaker, e-mail: bulava@imada.sdu.dk 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. In each channel, we can define a quantity q 2 cm,a by q 2 cm,a + m 2 1a + q 2 cm,a + m 2 2a = E cm .
Since E 2 cm = E 2 − P 2 must be real, a real solution for q 2 cm,a exists 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 [1][2][3][4]: In an orthonormal basis of states, each labelled by |Jm J LS a , 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-ofmomentum frame (not to be confused with the lattice length here), and S is the total spin of the two particles (not the scattering matrix). The index a is generalized to refer to species, the spins s 1 , s 2 , intrinsic parities η P 1 , η P 2 , isospins I 1 , I 2 , isospin projections I z1 , I z2 , and possibly G-parities η G 1 , η G 2 of particle 1 and 2. The F (P) matrix in this basis is given by where j 1 m 1 j 2 m 2 |JM are the familiar Clebsch-Gordan coefficients, and the W (Pa) matrix elements are given by The Rummukainen-Gottlieb-Lüscher (RGL) shifted zeta functions Z lm , introduced in Refs. [1,2], are evaluated as detailed in Ref. [5].

The K-matrix and box matrix
Eq. (5) is a single relation between a lab-frame finite-volume energy E and the entire S -matrix. This single relation is insufficient to extract all of the S -matrix elements at energy E. We proceed by the usual method of approximating 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. The S -matrix in Eq. (5) 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 [6,7], defined by 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. The multichannel generalization [8][9][10] of the effective range expansion is where K −1 L S a ; LS a (E cm ) is a real, symmetric, and analytic function of the center-of-momentum energy E cm . The effective range expansion given in Eq. (11) suggests the convenience of writing since K −1 L S a ; LS a (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. (5) can be written where we define the box matrix by 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. (13) and (15) 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 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. [11]) expresses the K-matrix as a sum of poles with a background described by a symmetric matrix of polynomials: where the couplings g (J p) α are real and the background coefficients d (Jk) αβ 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 four-momentum of particle j.

Block diagonalization
So far, we have expressed the matrices F (P) and B (P) in terms of the basis states labelled by |Jm J LS a . In this basis, the quantization condition in each of Eqs. (5), (13) and (15) is difficult to use since the determinant of an infinite matrix must be evaluated. By transforming to a basis in which both B (P) and K are block diagonal, we can focus on the determinant separately in each block. Each block has infinite dimension, but by truncating in the orbital angular momentum, keeping only states with L ≤ L max , each truncated block has a finite and reasonably small size.
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 If G is an element of the little group of P, then G P = 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 LS a . These symmetry operations only mix states of different m J . To block diagonalize B (P) , we apply a particular unitary change of basis: 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. (14), we find that we can write Λ λ n J L S a | B (P) |ΛλnJLS a = δ Λ Λ δ λ λ δ S S δ a a B (PΛ B S a) J L n ; JLn (E).
Notice that in Eq. (24) 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. The relationships of Λ B to Λ when η P 1a η P 2a = −1 are summarized in Table 1 of Ref. [5].
We have determined expressions in terms of the RGL shifted zeta functions 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 [12].
Lastly, we need to express K in the new basis. Given Eq. (16) and the orthonormality of the states in both the |Jm J LS a basis and the block diagonal |ΛλnJLS a basis, one can show that where η = (−1) L and η = (−1) L . If η = −η , the situation is much more complicated. However, in QCD, we should never need such matrix elements. 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. Thus, for a given P, we can label the quantization blocks of 1 − B (P) K and K −1 − B (P) in the |ΛλnJLS a basis solely by the irrep label Λ, where Λ is the irrep associated with the K-matrix.

Fitting
Let κ j , for j = 1, . . . , N K , denote the parameters that appear in the matrix elements of either the Kmatrix or its inverse K −1 . Once a set of energies for a variety of two-particle interacting states is determined, the primary goal is then to determine the best-fit estimates of the κ j parameters using the quantization determinant, as well as to determine the uncertainties in these estimates.
One method, which we call the spectrum method, is described in Ref. [5], but it is very difficult to implement. A much simpler method, known as the determinant residual method, is advocated in Ref. [5]. In this method, we introduce the quantization determinant itself as a residual in the correlated χ 2 to be minimized. 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. Determinants are susceptible to round off errors, which can make numerical minimization difficult. 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 which are much smaller in 5 EPJ Web of Conferences 175, 05005 (2018) https://doi.org/10.1051/epjconf/201817505005 Lattice 2017 Table 1. First tests of the determinant residual method applied to the interacting ππ energies in the I = 1 nonstrange channel described in Ref. [13]. These energies were obtained on a 32 3 × 256 anisotropic lattice with m π ≈ 240 MeV. In Ref. [13], the number of energy levels used was N E = 19. The fits below use N E = 20 by including an additional energy from a B + 1 d 2 = 1 irrep in which the leading partial wave is L = 3. The fits shown used Ω(µ, K −1 − B) as the residuals.
µ N E m ρ /m π g m 7 π a 3 m 11 π a 5 χ 2 /dof 1 20 3.338 (13)  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 study, we chose µ by starting with µ = 1, then increasing µ until the χ 2 value at the minimum did not change very much as µ was further increased.
In this method, the model fit parameters are just the κ i parameters, and the residuals are chosen to be r k = Ω µ, 1 − B (P) (E (obs) cm,k ) K(E (obs) cm,k ) , or the matrix K(E (obs) cm,k ) −1 − B (P) (E (obs) 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.

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. [13]. These energies were obtained on a 32 3 × 256 anisotropic lattice with m π ≈ 240 MeV. Defining k 0 = 2π/(m π L), the fit forms we used are Some of our results are presented in Table 1. In Ref. [13], the number of energy levels used was N E = 19. In the fits shown in Table 1 irrep in which the leading partial wave is L = 3. In the fits listed, 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 1 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. [13]. This is consistent with a phenomenological determination of m 7 π a 3 = 5.65(21) × 10 −5 taken from Ref. [14].
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
This talk reported on work completed in Ref. [5] to provide explicit formulas, software, and new fitting implementations for carrying out two-particle scattering studies using energies obtained from 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 K-matrix. The effective range expansion was used to introduce a smooth, well-behaved matrix K −1 . We obtained explicit expressions for the box matrix elements 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. More importantly, the software for evaluating all of these box matrix elements was made available. Lastly, we described a fitting strategy 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.