Efficient solution of the multi-channel L\"uscher determinant condition through eigenvalue decomposition

We present a method for efficiently finding solutions of L\"uscher's quantisation condition, the equation which relates two-particle scattering amplitudes to the discrete spectrum of states in a periodic spatial volume of finite extent such as that present in lattice QCD. The approach proposed is based on an eigenvalue decomposition in the space of coupled-channels and partial-waves, which proves to have several desirable and simplifying features that are of great benefit when considering problems beyond simple elastic scattering of spinless particles. We illustrate the method with a toy model of vector-vector scattering featuring a high density of solutions, and with an application to explicit lattice QCD energy level data describing $J^P=1^-$ and $1^+$ scattering in several coupled channels.


I. INTRODUCTION
Hadron spectroscopy is enjoying a renaissance driven by a global experimental program producing new and unexpected excited states. Excited states are observed as resonant enhancements in hadron-hadron partial-wave scattering amplitudes, a common quantity appearing in both experiment and theory that allows these resonances to be understood from the fundamental equations governing hadron interactions, Quantum Chromodynamics (QCD).
One approach for making predictions from QCD is lattice QCD, in which a finite Euclidean spacetime is discretized and the path integrals of the quantum field theory are sampled numerically using Monte-Carlo techniques. This approach makes only controlled and improvable approximations, and enables correlation functions and thus the finite-volume spectrum to be computed. Provided sufficiently large volumes and sufficiently small lattice spacings are used, reliable predictions can be made.
In these calculations, hadron-hadron scattering amplitudes are inferred from their effect on the discrete spectrum of states in a periodic spatial volume, usually cubic, defined by the lattice. The precise relationship is encoded in the Lüscher quantization condition [30][31][32][33][34][35][36][37][38][39][40][41][42][43], which can be written as 1 In this expression, M is a matrix of known kinematic factors that encodes the dependence on the L×L×L volume, ρ is a diagonal matrix of phase-space factors, and t is the t-matrix describing scattering, related to the S-matrix by S = 1 + 2i √ ρ · t · √ ρ. The matrix indexing runs over hadron-hadron channels and partial-waves, with a version of Eqn. 1 existing for each irreducible representation (irrep) of the relevant symmetry group. Because the symmetry of the cubic boundary is reduced with respect to the full rotation group, each irrep features multiple partialwaves 2 . The space of hadron-hadron scattering channels is determined by those channels which are kinematically accessible in the energy region of interest. The complex scalar function det D(E cm ) features not only zeroes corresponding to the finite-volume spectrum, but also divergences which originate in the matrix M, that appear at each energy where a hadron-hadron pair could go on-shell in a theory without interactions, i.e. where the individual hadrons have allowed finite-volume momenta, 2π L (n x , n y , n z ) with n i ∈ Z. Considering Eqn. 1, for a given scattering matrix, t(E cm ), one can determine the finite-volume spectrum 1 This is valid for any two-hadron to two-hadron scattering process -significant recent steps towards a general quantization condition relevant for three-hadron channels have been made [44][45][46][47][48][49][50][51][52][53]. 2 Formally an infinite number of partial-waves are present, but higher angular momenta are suppressed at low energies by barrier factors and may be neglected in practice. Explicit relations for the subduction in the case of systems with net momentum can be found in [54].

arXiv:2001.08474v1 [hep-lat] 23 Jan 2020
by finding all zeroes in a desired energy region, and in practice this is done numerically using root-finding algorithms applied to the scalar function det D(E cm ) . This requires computation of the elements of the matrix M which comes with some non-negligible computational overhead, so approaches are generally preferred that evaluate the function less often. It is common to encounter cases where hadron-hadron interactions are relatively weak in some energy regions, such that zeroes lie very close to noninteracting energies where det D(E cm ) diverges. This can pose a challenge to conventional approaches to rootfinding which rely upon bracketing a zero by locating function values of opposite sign either side of the zero. The hadspec collaboration has pioneered an approach for determining coupled-channel scattering amplitudes in which the energy dependence of the t-matrix is parameterized, and the best description of the lattice QCD computed finite-volume spectrum is found by varying the free parameters in the scattering amplitude. In practice this requires finding the solutions of Eqn. 1 for each considered volume and irrep, at every iteration of the scattering amplitude parameters which are being varied under a χ 2 minimization, where the χ 2 measures the difference between the found 'model' spectrum (the zeroes) and the lattice spectrum. In this approach it is vitally important to find all the relevant solutions at every iteration, as if any are missed in any iteration, the sampling of the χ 2 surface will feature discontinuous jumps, which are likely to cause failure in the minimization routine.
A more fundamental challenge arises when stable particles with nonzero spin feature in the scattering process. In these cases it becomes possible for there to be zeroes of det D(E cm ) with multiplicity larger than one [24,55], which can appear as the function touching zero rather than crossing zero, as here bracketing will fail altogether. If the zero is found by some non-bracketing approach [56], it is generally not immediately clear what its multiplicity is. As interactions are typically of different strengths in different partial-waves, this degeneracy of zeroes is broken, but this leads to multiple zeroes appearing in a very small energy range. This remains a challenge for bracketing approaches: for example, suppose two zeroes are separated by a small energy, and the bracketing 'grid' samples points such that the two closest neighbouring points lie below the first zero and above the second zero, no change of sign will be observed, and no hint provided that any zeroes are present.
In this paper, we will present a pragmatic approach to finding zeroes of the determinant condition that is both computationally fast and reliable. The approach makes use of the eigenvalue decomposition of a suitably transformed version of the matrix, D, appearing under the determinant in Eqn. 1, using the well-known result that a zero of det D corresponds to at least one eigenvalue of D taking a zero value. In practice we find that considering the energy-dependence of the eigenvalues reduces an n-channel problem to one which is no more difficult to solve than n independent one-channel problems, where bracketing approaches typically work well.
Several considerations are made in search of a practical implementation. Firstly, whether there are transformations on D defined in Eqn. 1 that preserve the zeroes of det D while also simplifying the numerical problem of finding the zeroes. Secondly, how one should match the eigenvalues between neighbouring energy steps to ensure the λ p (E cm ) are continuous functions, and finally how to handle isolated special points like kinematic thresholds, and energies at which the eigenvalues diverge.
In the remainder of the paper, after describing our implementation, we will illustrate the method using a toy model of scattering of two vector mesons featuring very dense zeroes, before presenting results from an explicit lattice QCD calculation of I = 1, G = + scattering of ππ, KK, πω and πφ, where both J P = 1 + and J P = 1 − amplitudes appear.

II. EIGENVALUE DECOMPOSITION
We consider Eqn. 1 where we have truncated the space of partial-waves and hadron-hadron scattering channels to some finite number, rendering D(E cm ) an n × n matrix whose determinant can be decomposed as the product of its n eigenvalues λ p (E cm ), defines the corresponding eigenvectors 3 . Solutions of det [D(E cm )] = 0 occur when at least one eigenvalue takes a zero value, so our approach will be to find each eigenvalue as a function of energy and then perform root-finding on each independently. In order to do this, the eigenvalues must be reliably matched at each step in energy over the entire range considered, otherwise 'jumps' due to mismatching could be mistaken for true zeroes and vice-versa, true zeroes could be missed.
For the form of the quantization condition given in Eqn. 1, at least one of the eigenvalues diverges at every non-interacting energy owing to a pole in M, originating in the Lüscher zeta-functions 4 . Such a divergence can present a problem when attempting to match the eigenvalues or eigenvectors at energies close to these divergences -the assumption that a small change in energy leads to a small perturbation of D(E cm ) is invalid near the non-interacting energies as the matrix is unbounded.
In those scattering processes where interactions are relatively weak, it is very common that solutions of the quantization condition lie in a neighbourhood of these non-interacting energies, so removing these divergences is a priority.
Fortunately it is quite straightforward to transform D(E cm ) in such a way as to remove the non-interacting divergences while leaving the zeroes of det D(E cm ) unchanged. For example, where V (E cm , L) = (1 + iM)(1 − iM) −1 . The divergences at non-interacting energies have been explicitly factorised off in (1−iM), leaving V which is free of these divergences. It follows that the zeroes of det D(E cm ) are also zeroes of det D V (E cm ) , where A closely related alternative 5 is which is also finite at non-interacting energies. Note that like M, both V and U are diagonal in hadron-hadron channel but dense in partial-wave space. The properties of the matrices S, U and V depend upon energy, and in particular whether any channel is kinematically closed at any given energy. We begin by considering the case of energies high enough that all included channels are kinematically open.
thr. < E (2) thr. < · · · < E (m) thr. . For each hadron-hadron channel a, there may be multiple partial-waves under consideration. Suppose channel a has n a partial waves, then m a=1 n a = n, the dimension of the matrix D(E cm ). 6 For energies E cm E (m) thr. , the finite-volume matrix M(E cm , L) is hermitian 7 , and then since iM is skewhermitian, it has purely imaginary eigenvalues meaning that (1 − iM) −1 exists for all such energies. It immediately follows that V is a unitary matrix for all E cm E (m) thr. . As S is also unitary in this energy region by conservation of probability, commonly referred to as unitarity, D V (E cm ) is diagonalisable with n linearly independent 5 Which appeared originally in the elastic case in Ref. [30]. 6 Some partial waves may appear embedded more than once. 7 See Appendix A. orthonormal eigenvectors and bounded eigenvalues of the form where the impact of unitarity is to ensure that the complex eigenvalues are each controlled by a single real parameter, θ p . Since D V (E cm ) is continuous in E cm , a small perturbation in energy, E cm = E cm + ∆E cm , which we will refer to as an iteration in energy, results in a small perturbation of the matrix, . It immediately follows from standard matrix perturbation theory that the eigenvalues and eigenvectors at the next iteration, are approximated in terms of the eigenvalues and eigenvectors at the current iteration, so that for the case of non-degenerate eigenvalues, the eigenvectors are indeed small perturbations for ∆E cm sufficiently small. The inner product between eigenvectors on consecutive energy iterations is given by where evaluated at E cm will be 1 for a sufficiently small choice of ∆E cm . It follows that eigenvectors (and their corresponding eigenvalues) can be matched between energy iterations by examining Eqn. 6, with the order of the eigenvectors at E cm = E cm + ∆E cm taken so as to maximise the sum of the norm-squared of the inner products with the eigenvectors at E cm .
In the case that a number of eigenvalues are degenerate or very nearly degenerate, i.e. when λ p − λ k = O ||∆D V || for at least one k = p, matching the corresponding eigenvectors between energy iterations is potentially ambiguous. As discussed in Appendix B, in these cases it is common to observe an avoided crossing, where the eigenvalues repel each other rather than cross, but crossings are possible when channels or partial-waves are decoupled. For an energy E cm where the eigenvalues are degenerate and cross, the eigenvectors span the corresponding eigenspace and there is an arbitrariness in their direction. At this energy, it is not necessary to match the eigenvectors to the previous iteration, as we are interested in the eigenvalues, which are equal. At the next energy iteration, E cm + ∆E cm , where the eigenvalues are now distinct, rather than matching to the previous set of eigenvectors at E cm , we can simply match to those at E cm − ∆E cm .
FIG. 1. A schematic example of how the various determinant equations are used to find roots of the quantisation condition for a scattering system with a single hadron-hadron channel in an arbitrary number of partial-waves. For a particular energ,y (E cm )1, L)] = 0, so an energy region is defined in which det S + U is considered rather than det 1 + S · V .

B. With at least one closed channel
One might expect that the finite-volume spectrum at energies below the kinematic threshold for some channel a to be insensitive to amplitudes describing channel aafter all, there are an infinite number of channels whose thresholds lie higher than any given energy, and we cannot ever possibly know the amplitude in all of them. In fact, the finite-volume spectrum obtained from Eqn. 1 is only sensitive to closed channels in a very small energy region below the threshold, with the effect of the channel reducing exponentially (as e −κL where κ is the magnitude of the imaginary scattering momentum). Nevertheless it is important to account for these closed-channel effects in practical calculation, particularly when a resonance appears near a threshold.
At energies below the highest threshold under consideration, E cm E (m) thr. , the finite-volume matrix M is in general no longer hermitian and subsequently V is no longer unitary. (1 − iM) is not prevented from having zero eigenvalues and correspondingly V may diverge at some energies.
Recalling that M is block-diagonal in channel space, we can solve for roots of det[1 − iM] = 0 in each channel separately. At decreasing energies below the threshold for channel a, M aa tends to i1 exponentially [58] and so sufficiently far below threshold, det[1−iM aa ] → 2 na . For each channel a root can only be found below threshold, as above threshold we recall M aa is hermitian and no solutions to det[1 − iM aa ] = 0 exist. In order to ensure that zeroes of multiplicity higher than one are found, it is good practice to make use of an eigen-decomposition of M aa .
A practical approach to dealing with these isolated energies at which V diverges, is to eigen-decompose instead the matrix D U in a small region around the divergent energy, since U is typically well defined at those energies. If it should happen by coincidence that det[1 + iM aa ] has a zero very close to where det[1 − iM aa ] has a zero, then we could eigen-decompose the original matrix D. We stress that in all practical applications, this finely tuned scenario seldom occurs, and can be avoided by considering sufficiently small energy regions. Fig. 1 provides a pictorial representation of the approach.
Finally, there can be divergences in the determinant condition below the lowest threshold, where no channel is kinematically accessible, when the scattering system features stable bound-states which manifest as simple pole singularities in S(E cm ). The divergence appears at the volume-independent bound-state energy, which differs from a zero corresponding to a finite-volume energy level only by an exponentially small finite-volume correction.

C. At kinematic thresholds
The energy at which a new channel opens, its threshold, requires special consideration when D V or D U are used. If the threshold coincides with a non-interacting energy, such as in S-wave scattering in the rest frame, then there is a divergence due to a pole singularity of finite order in D when E cm = E thr. , which is removed in D V and D U .
Recall that the matrices U and V are block-diagonal with respect to hadron-hadron scattering channels, being dense in the partial-wave space only. The threshold divergence in M is such that at threshold V and U are diagonal and equal to −1.
In general S is block-diagonal in partial-wave space, but dense in hadron-hadron channels, however at the opening of a new threshold, say for channel a, the vanishing of the phase-space for this channel means that for any partial-wave, S aa (E thr. ) are unconstrained. It follows that at each threshold there will be at least one zero eigenvalue of D V or D U , but these zeroes should not be mistakenly considered to be actual finite-volume energy levels -as they appear at exactly the known threshold energies, there should be no risk of this misidentification.

D. Summary of the approach
Standard eigen-decomposition routines can be applied to D V , D U or D (whichever is appropriate given the discussion in Section. II B) on a discrete sampling of energies, and the eigenvalues λ p (E cm ) matched between energy iterations using similarity of the corresponding eigenvectors established through magnitude of inner products. Each eigenvalue can then be considered separately, and typically we find that root-finding on any one of these eigenvalues is no more difficult than in the case of elastic scattering.
Energy regions in which root-finding is performed are separated by hadron-hadron thresholds and by regions around divergences in V , but as these regions are disjoint, the complete set of solutions is the combination of the sets of solutions of each eigenvalue in each energy region. It is therefore not necessary to match eigenvalues across boundaries.
The combination of the sets of solutions in each eigenvalue gives the full set of roots of det D(E cm ) = 0, and of course with the eigenvalue roots in hand it is trivial to check that they are in fact roots of det D(E cm ) .

III. APPLICATION: VECTOR-VECTOR SCATTERING
Vector-vector scattering provides a useful demonstration of the eigen-decomposition approach described in the previous section, as even in the simple case of total momentum zero we can have a situation in which solutions of high multiplicity appear in the non-interacting limit, while in the interacting case the corresponding zeroes may be very closely spaced.
We will illustrate our approach using a simple toy model of vector-vector scattering in which we have a single hadron-hadron channel comprising two identical vector mesons. Considering the finite-volume spectrum in the [000] E + irrep, we note that Bose symmetry constraints limit the number of contributing partial-waves, and we summarize those that can be non-zero for J 4 in Table I.
We engineer a dense spectrum of energy levels by considering relatively weak interactions in all partial-waves, truncating at = 2. As indicated in Table I, this results in the contribution of four partial waves, 5 S 2 , 1 D 2 , 5 D 2 and 5 D 4 .
The vector particles are taken to have mass m = 0.5, and we consider a volume L = 70. In this case the noninteracting energies and the corresponding multiplicities are presented in Table II, which indicates that we expect for a weakly interacting scattering system to find one level close to threshold, three closely-spaced levels somewhat higher in energy and four more levels a little higher still.
Since there may be finite-volume energy levels below the opening of our threshold, we should first establish if there are energies around which we should use D U rather than our preferred choice D V . This can be done without knowledge of the scattering amplitude, by looking for roots of det[1 − iM] = 0 below threshold. Fig. 2 shows the eigenvalues of 1 − iM and their roots. Note that each eigenvalue diverges at the threshold, and asymptotically approaches a value of 2 (as M → i1) as energies go further below threshold. Two roots are found: a double root at (E cm ) 1 = 0.9962 and a single root at (E cm ) 2 = 0.9988. A small energy region 8 around each of these values is chosen to be treated with D U in place of D V .
To describe scattering in this system of four partialwaves we use a K-matrix parameterisation which automatically implements the required unitarity of the S-matrix. The K-matrix, K(s), where s = E 2 cm , is a real symmetric matrix, related to the t-matrix via Unitarity is guaranteed if Im I(s) = −ρ(s) for energies above threshold and Im I(s) = 0 below. We make use of the Chew-Mandelstam prescription [60] which defines Re I(s) through a dispersive integral featuring ρ(s), and choose to subtract so that Re I = 0 at threshold.
thr. . The imaginary part of each eigenvalue is identically zero below threshold.
FIG. 3. Phase-shifts and mixing-angles for the toy model amplitude given in Eq. 8 following the n-channel Stapp parameterisation described in Ref. [24].
To implement weak scattering we use a constant K-matrix parameterisation featuring small values, where the block diagonalization in J is made explicit. The somewhat larger values in the D-wave amplitudes are compensated by the threshold factors k cm in Eq. 7. The resulting phase-shifts and mixing-angles for this toy amplitude are presented in Fig. 3 using the n-channel generalization of the Stapp parameterisation described in Ref. [24]. Note that all amplitudes are modest by design, with the S-wave an order of magnitude larger than any of the D-wave amplitudes.
Following the procedure set out in Sec. II, we seek to find solutions of det[D] = 0. The top panel of Fig. 4 shows det[D] as defined in Eqn. 1 as a function of energy, where divergences at each of the non-interacting energies can be clearly seen, along with the near-degenerate roots located in very close proximity -as anticipated from our weakly interacting amplitudes. Finding the zeroes directly from this determinant requires ultra-fine sampling in energy, and roots are easily missed in automated root-finding algorithms. Furthermore, the toy amplitudes could be easily modified to make the near degenerate roots, shown in the zoomed regions, exactly degenerate, in which case no amount of improvement in sampling resolution would find the appropriate number of roots.
The second panel in Fig. 4 shows det[D V ] as a function of energy. As discussed in Sec. II, the divergences that appear at the non-interacting energies have been removed and above threshold the determinant is bounded. In a neighbourhood around the near degenerate roots, the determinant is more clearly seen to be locally quadratic and cubic respectively.
The remaining lower panels in Fig. 4 show the four eigenvalues of det[D V ] over the considered energy range, about which we make a number of observations. Firstly, the real part of each eigenvalue is in the interval [0, 2] above threshold and the imaginary component in [−1, 1] as proven in Sec. II A. Below threshold, the eigenvalues indeed diverge at energies where V is singular as expected -see Sec. II B. The multiplicities of these divergences also agree, for example as shown in Fig. 2 the determinant det[1 − iM] has a double root at (E cm ) 1 , which manifests in the eigenvalues λ 2 and λ 3 , and a single root at (E cm ) 2 , which appears in λ 1 . At threshold, each eigenvalue is identically zero as proven in Sec. II C and these threshold zeroes are removed from the set of solutions. It is clear from counting the number of roots that we have agreement with the multiplicities of the non-interacting energies, shown in Table II, as we expect for weak interactions. Most notably, finding roots in each eigenvalue is no more difficult than solving for roots in elastic scattering and this is illustrated by the simplicity of the zero crossings in the eigenvalues. The found roots are recorded in Table III.  Table II. IV. APPLICATION: 1 + , 1 − SCATTERING Previous lattice QCD calculations have determined the lowest lying isospin=1 J P C = 1 −− , 1 +− vector [8] and axial-vector [24] resonances, the ρ and b 1 , in studies on three lattice volumes 9 with a pion mass ∼ 391 MeV. In the determination of the b 1 resonance, energy levels used to constrain the scattering amplitudes were obtained in the [000] T + 1 and ( P = 0) A 2 irreps, as these irreps receive no contribution from J P = 1 − scattering amplitudes. This avoids the complication of disentangling the lowlying, resonating 1 − channels that mix due to the reduced symmetry of the finite-volume. The eigen-decomposition method we have introduced in this manuscript enables us to lift such practical restrictions and perform a combined scattering analysis of the 1 −− and 1 +− sectors. Using a much larger set of irreps, including now those in which both parities feature, allows us to significantly increase the number of energy levels constraining the scattering amplitudes, with a total of 144 energy levels being available.
Several hadron-hadron channels in various partial-waves feature in these scattering systems. For the J P = 1 − sector, in addition to ππ{ 1 P 1 } considered in Ref. [8] (which restricted itself to elastic energies), we include KK{ 1 P 1 } and πω{ 3 P 1 } which are kinematically open in the energy region where where the b 1 resonance appears. For J P = 1 + , we consider πω{ 3 S 1 }, πω{ 3 D 1 } and πφ{ 3 S 1 } as in Ref. [24]. As discussed in great detail in Ref. [24], three-body thresholds, ππη and πKK, open below the 4π-threshold in the region where the b 1 resonates and we ensure that in each irrep, we keep below the lowest E 2+1 n.i non-interacting energy 10 .
Other partial-waves also feature due to the reduced symmetry of the finite-volume. For ππ and KK with isospin=1, a 1 F 3 wave can appear (see Table III of Ref. [8]) but will be suppressed by the high value of . For vectorpseudoscalar scattering, the nonzero intrinsic spin gives rise to a triplet of partial-waves for each 1, e.g. 3 P 0 , 3 P 1 and 3 P 2 for = 1 (see Table 1 of Ref. [55] for irreps at rest, and Tables 5-7 for irreps in flight). Following the analysis in Ref. [24], we include only πω{ 3 P 1 }, with J P = 1 − which has not been previously shown to be compatible with zero. The resulting set of partial waves to be considered is presented in Table IV.
1 + πω{ 3 S1} πω{ 3 D1} πφ{ 3 S1} As in the case of the toy amplitudes in Sec. III, singularities of V must be determined in order to solve the appropriate form of quantisation condition in each energy region. For reference, tables VI -VIII in Appendix C report the singularities of V on each volume and lattice irrep.
In this case we are seeking to describe an explicit spectrum of lattice QCD energy levels (with uncertainties and correlations) using a parameterised amplitude in a χ 2 minimization 11 . In the approach followed by hadspec, typically a range of parameterisations would be considered, but in this case, where we seek only to demonstrate the effectiveness of the eigen-decomposition approach, we 10 As introduced in Ref. [24], E 2+1 n.i non-interacting energies are finite-volume energies calculated in a two-meson subsystem that are then combined with the third meson, where no interactions are assumed between the two-meson subsystem and third meson. 11 The χ 2 is defined explicitly in Eqn. 9 of Ref. [8].
use just one 12 , a K -matrix parameterisation of the 'pole plus constant' form for each J P , πω{ 3 S1} g πω{ 3 S1} g πω{ 3 D1} 0 g πω{ 3 S1} g πω{ 3 D1} g 2 We use the Chew-Mandelstam prescription for I(s), choosing to subtract at the pole masses so that Re I aa (m 2 − ) = 0 where a runs over channels with J P = 1 − and similarly Re I aa (m 2 + ) = 0 for channels with J P = 1 + . The use of only a single pole in each J P is motivated by the expectation of a narrow ρ resonance, a narrow b 1 resonance, and no other resonances until much higher energies [63,64].
The best description of the lattice spectra in ir- A 2 ), varying the ten parameters in Eqn. 9, is found with a quite reasonable χ 2 /N dof = 229./(144 − 10) = 1.71, after ∼ 400 iterations of the MINUIT minimisation routine [65]. Fig. 5 shows as orange curves the finite-volume spectrum corresponding to the best-fit amplitude of Eqn. 9 for those irreps in which both J P = 1 − and 1 + are subduced. Also shown are the relevant non-interacting energy curves and the lattice QCD energy levels which were fitted to obtain the best-fit parameter values.
In Fig. 6 we illustrate the solution of the quantization condition in the case of the [111] E 2 irrep on the L/a s = 24 lattice using the central values of the best fit parameters. A single embedding of each partial-wave in each hadron-hadron channel recorded in Table IV features in this irrep, and thus D V is a 6 × 6 matrix. We make a number of observations. First, as demonstrated in the toy model example in Sec III, the multiplicities of the singularities of V (see Table VIII) are manifested in the eigenvalue divergences, i.e. three singularities with multiplicity one. The determinant of D V is indeed bound above the πφ threshold and moreover, as no singularity of D V appears above πω threshold in this case, D V is also bound between πω and πφ threshold. For this particular parameterisation, the S-matrix is decoupled in hadron channel (see Eq. 9). As V is diagonal in hadron channel, D V is block-diagonal and the eigenvectors are therefore orthogonal in channel space. This is manifest in the corresponding eigenvalues and we can make the following associations: λ 1 corresponds to ππ{ 1 P 1 }, λ 2 to KK{ 1 P 1 }, λ 3 , λ 4 , λ 5 to πω{ 3 S 1 , 3 P 1 , 3 D 1 } and λ 6 to πφ{ 3 S 1 }. Of course this is not a general feature of the eigen-decomposition method as the S-matrix is generally dense in channel space. The roots and the corresponding eigenvalues of D V in which they appear are recorded in Table V. irrep on the L/as = 24 lattice for the amplitude in Eqn. 9 listed by the eigenvalue of D V in which they appear. Grey entries denote roots that were found at energies higher than those we are considering.
For this parameterisation, we find the ππ{ 1 P 1 } amplitude appears to agree very well with the amplitudes determined in the previous calculation [8], with the energy dependence above the resonance being quite modest. The amplitude features a complex conjugate pair of resonance poles, which we interpret as the ρ, at, which is in quite reasonable agreement with the value found previously using a more limited set of irreps. The KK{ 1 P 1 } and πω{ 3 P 1 } amplitudes were found to be weak over the entire energy range considered, with phaseshifts not exceeding 1 • and 17 • respectively, as anticipated owing to the absence of any 1 − resonances between the KK and 4π thresholds at this pion mass. Similarly, the J P = 1 + amplitudes were found to agree very well with those found in the previous calculation [24]. The amplitude features a complex conjugate pair of resonance poles, which we interpret as the b 1 , at, a t E cm = 0.2454(9) ± i 2 0.0199 (12). (11) in close agreement with the pole determined in Ref. [24].

V. SUMMARY
We have presented an approach for efficiently solving the Lüscher quantization condition which determines the discrete spectrum of states in a finite-volume given a scattering system. The technique makes use of the eigenvalue decomposition of the matrix appearing under the determinant in the quantization condition, ultimately reducing the complexity of the numerical problem for coupled-channels and coupled partial-waves down to one of root-finding in a number of independent continuous functions which typically each feature far less rapid variation than does the determinant itself.
Two applications were given for illustration. In the first, a toy model of weak scattering in a system of two vector mesons shows the ability of the approach to determine the very dense set of zeroes lying close to non-interacting levels of high multiplicity. In the second, a set of lattice QCD calculated energy levels on three volumes were used to constrain coupled-channel scattering with J P = 1 − and 1 + . This provides a stress-test of the approach, as the solving of the quantization condition needs be successful on multiple irreps in multiple volumes, at many hundreds of iterations of the scattering amplitude parameter values, varied under a χ 2 minimization controlled by a minimization routine. The approach was found to be both reliable and computationally fast.
Contemporary examples where application of this approach may prove to make possible previously impractical calculations include the determination of the resonance spectrum in charmonium where a decay product vector meson D * is well approximated as a stable hadron, and where the DD, D * D and D * D * thresholds lie very close to each other. In-flight irreps will receive contributions from many partial-waves, and there is the expectation of resonances in all of J P C = 1 −− , 2 −− , 3 −− with rather similar masses. Nucleon resonances provide another application where a relatively large number of scattering channels and partial-waves contribute in the irreps of interest. We expect our approach to significantly simplify the finding of finite-volume energy levels in systems like these, and hence to make possible the determination of larger parts of the resonance content of QCD.

ACKNOWLEDGMENTS
We thank our colleagues within the Hadron Spectrum Collaboration. AJW is supported by the U. K Bounding of the eigenvalues and orthogonality of the eigenvectors of D V followed from hermiticity of M above threshold. Here we show this hermiticity, starting with the general expression for matrix elements M in the basis of intrinsic spin S, orbital angular momentum , total angular momentum J, azimuthal component of total angular momentum m and channel index a, The spin-orbit coupling is expressed through the real SU(2) Clebsch-Gordan coefficients m ; Sm S |Jm and the volume dependence is encoded in the functions c n ,m (k 2 cm ; L), defined as where n = L 2π P , and k cm is the cm-frame momentum corresponding to E cm . Definitions of the boost-factor, γ, and the construction of the set of real vectors to be summed over, r ∈ P n , can be found in Ref. [43]. Using Y * ,m = (−1) m Y ,−m , it is easy to show that c n * ,m = c n ,−m . The presence of the factor k −(¯ +1) cm in Eqn. A1 distinguishes the above and below threshold cases: above threshold k cm is real and this factor is unchanged by complex conjugation, below threshold k cm is imaginary, and complex conjugation can have an effect.
The integral over the product of three spherical harmonics can be expressed in terms of Clebsch-Gordan coefficients, and thus the integral is manifestly real. We leave it in the integral form as it more directly exposes the relevant symmetries. Under complex conjugation of Eqn. A1 we can replace the resulting c n * ,m with c n ,−m and Y¯ ,m with (−1)m Y * ,−m , and then noting that the sum overm is equivalent to a sum over −m obtain in the above threshold case where k cm is real, Eigenvalues of D V = 1 + S · V are observed to exhibit a behavior where they avoid crossing each other as E cm is varied. Examples from the toy model considered in Section III are shown in Figure 7. This behaviour is in fact a general property that follows from the unitarity of the matrix S · V above threshold.
Defining unitary matrices at neighbouring energy points, W = S(E cm ) · V (E cm ) and W = S(E cm + ∆E cm ) · V (E cm + ∆E cm ), we can find hermitian matrices, H, H such that W = exp iH W = exp iH , and clearly the eigenvalues and eigenvectors of these matrices are trivially related, with θ p being real, and simply related to the complex eigenvalues λ p of D V as given in Eqn. 5. The avoided level crossings that manifest in both the real and imaginary parts between some λ p 's therefore appear also between the corresponding θ p 's.
By choosing to examine H we reduce the problem to one familiar in quantum mechanics: as H and H − H are hermitian this is the case considered in (degenerate) perturbation theory, where 'avoided level crossings' are a generic feature whenever the perturbation (H − H) connects the (degenerate) eigenstates. A classic illustrative example is the two-state case where in the eigenbasis of H, H = θ 1 0 0 θ 2 + P 11 P 12 P * 12 P 22 , and the corresponding exact eigenvalues of H are θ 1 = a + b, θ 2 = a − b, where a = 1 2 θ 1 + θ 2 + P 11 + P 22 b = 1 2 θ 1 − θ 2 + P 11 − P 22 2 + 4 P 12 2 .
It follows that unless both terms under the square root are zero, θ 2 − θ 1 > 0, and the eigenvalues do not cross.

Appendix C: Tables of singularities
In this appendix, we record the singularities of V (E cm ) in each irrep and volume relevant for the calculation presented in Sec. IV. We record singularities for irreps featuring subductions of J P = 1 − and not J P = 1 + in Table VI, irreps featuring J P = 1 + and not J P = 1 − in Table VII, and irreps featuring both in Table VIII. For clarity, the singularities arising in different hadron-hadron channels ('a') is made explicit in each of the tables.