Efficient solution of the multi-channel Lüscher determinant condition through eigenvalue decomposition
JJLAB-THY-20-3133
Efficient solution of the multi-channel L¨uscher determinant condition througheigenvalue decomposition
Antoni J. Woss, ∗ David J. Wilson, † and Jozef J. Dudek
2, 3, ‡ (for the Hadron Spectrum Collaboration) DAMTP, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge, CB3 0WA, UK Thomas Jefferson National Accelerator Facility, 12000 Jefferson Avenue, Newport News, VA 23606, USA Department of Physics, College of William and Mary, Williamsburg, VA 23187, USA (Dated: January 24, 2020)We present a method for efficiently finding solutions of L¨uscher’s quantisation condition, theequation which relates two-particle scattering amplitudes to the discrete spectrum of states in aperiodic spatial volume of finite extent such as that present in lattice QCD. The approach proposedis based on an eigenvalue decomposition in the space of coupled-channels and partial-waves, whichproves to have several desirable and simplifying features that are of great benefit when consideringproblems beyond simple elastic scattering of spinless particles. We illustrate the method with atoy model of vector-vector scattering featuring a high density of solutions, and with an applicationto explicit lattice QCD energy level data describing J P = 1 − and 1 + scattering in several coupledchannels. I. INTRODUCTION
Hadron spectroscopy is enjoying a renaissance drivenby a global experimental program producing new andunexpected excited states. Excited states are observedas resonant enhancements in hadron-hadron partial-wavescattering amplitudes , a common quantity appearing inboth experiment and theory that allows these resonancesto be understood from the fundamental equations gov-erning hadron interactions, Quantum Chromodynamics(QCD).One approach for making predictions from QCD is lattice QCD , in which a finite Euclidean spacetime is dis-cretized and the path integrals of the quantum field theoryare sampled numerically using Monte-Carlo techniques.This approach makes only controlled and improvable ap-proximations, and enables correlation functions and thusthe finite-volume spectrum to be computed. Providedsufficiently large volumes and sufficiently small latticespacings are used, reliable predictions can be made.The past decade has seen the state of the art in lat-tice QCD calculations move to consider scattering sys-tems of ever increasing complexity, from elastic cases like ππ [1–16], to the first calculations of coupled-channelscattering [17–21], to recent studies of multiple coupledchannels such as ππ, K ¯ K, ηη considering both J P = 0 + and 2 + which resonate at low energies [22], and systemswhich feature scattering hadrons of non-zero spin, suchas πω [23, 24] or πN [25–29].In these calculations, hadron-hadron scattering am-plitudes are inferred from their effect on the discretespectrum of states in a periodic spatial volume, usuallycubic, defined by the lattice. The precise relationship is ∗ [email protected] † [email protected] ‡ [email protected] encoded in the L¨uscher quantization condition [30–43],which can be written as det (cid:2) D ( E cm ) (cid:3) = 0 , where D ( E cm ) = + i ρ ( E cm ) · t ( E cm ) · (cid:0) + i M ( E cm , L ) (cid:1) . (1)In this expression, M is a matrix of known kinematic fac-tors 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 = + 2 i √ ρ · t ·√ ρ . The matrix indexing runs overhadron-hadron channels and partial-waves, with a ver-sion of Eqn. 1 existing for each irreducible representation( irrep ) of the relevant symmetry group. Because the sym-metry of the cubic boundary is reduced with respect to thefull rotation group, each irrep features multiple partial-waves . The space of hadron-hadron scattering channelsis determined by those channels which are kinematicallyaccessible in the energy region of interest.The complex scalar function det (cid:2) D ( E cm ) (cid:3) features notonly 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 paircould go on-shell in a theory without interactions, i.e.where the individual hadrons have allowed finite-volumemomenta, π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 This is valid for any two-hadron to two-hadron scattering process– significant recent steps towards a general quantization conditionrelevant for three-hadron channels have been made [44–53]. Formally an infinite number of partial-waves are present, buthigher angular momenta are suppressed at low energies by barrierfactors and may be neglected in practice. Explicit relations forthe subduction in the case of systems with net momentum canbe found in [54]. a r X i v : . [ h e p - l a t ] J a n by finding all zeroes in a desired energy region, and inpractice this is done numerically using root-finding algo-rithms applied to the scalar function det (cid:2) D ( E cm ) (cid:3) . Thisrequires computation of the elements of the matrix M which comes with some non-negligible computational over-head, so approaches are generally preferred that evaluatethe function less often. It is common to encounter caseswhere hadron-hadron interactions are relatively weak insome energy regions, such that zeroes lie very close to non-interacting energies where det (cid:2) D ( E cm ) (cid:3) diverges. Thiscan pose a challenge to conventional approaches to root-finding which rely upon bracketing a zero by locatingfunction values of opposite sign either side of the zero.The hadspec collaboration has pioneered an approachfor determining coupled-channel scattering amplitudesin which the energy dependence of the t -matrix is pa-rameterized, and the best description of the lattice QCDcomputed finite-volume spectrum is found by varying thefree parameters in the scattering amplitude. In practicethis requires finding the solutions of Eqn. 1 for each consid-ered volume and irrep, at every iteration of the scatteringamplitude parameters which are being varied under a χ minimization, where the χ measures the differencebetween the found ‘model’ spectrum (the zeroes) and thelattice spectrum. In this approach it is vitally importantto find all the relevant solutions at every iteration, as ifany are missed in any iteration, the sampling of the χ surface will feature discontinuous jumps, which are likelyto cause failure in the minimization routine.A more fundamental challenge arises when stable par-ticles with nonzero spin feature in the scattering process.In these cases it becomes possible for there to be zeroes ofdet (cid:2) D ( E cm ) (cid:3) with multiplicity larger than one [24, 55],which can appear as the function touching zero ratherthan 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 multiplicityis. As interactions are typically of different strengthsin different partial-waves, this degeneracy of zeroes isbroken, but this leads to multiple zeroes appearing in avery small energy range. This remains a challenge forbracketing approaches: for example, suppose two zeroesare separated by a small energy, and the bracketing ‘grid’samples points such that the two closest neighbouringpoints lie below the first zero and above the second zero,no change of sign will be observed, and no hint providedthat any zeroes are present.In this paper, we will present a pragmatic approach tofinding zeroes of the determinant condition that is bothcomputationally fast and reliable. The approach makesuse of the eigenvalue decomposition of a suitably trans-formed version of the matrix, D , appearing under thedeterminant in Eqn. 1, using the well-known result thata zero of det D corresponds to at least one eigenvalue of D taking a zero value. In practice we find that consid-ering the energy-dependence of the eigenvalues reducesan n -channel problem to one which is no more difficult tosolve than n independent one-channel problems, where bracketing approaches typically work well.Several considerations are made in search of a practicalimplementation. Firstly, whether there are transforma-tions on D defined in Eqn. 1 that preserve the zeroesof det D while also simplifying the numerical problem offinding the zeroes. Secondly, how one should match theeigenvalues between neighbouring energy steps to ensurethe λ p ( E cm ) are continuous functions, and finally how tohandle isolated special points like kinematic thresholds,and energies at which the eigenvalues diverge.In the remainder of the paper, after describing ourimplementation, we will illustrate the method using a toymodel of scattering of two vector mesons featuring verydense zeroes, before presenting results from an explicitlattice 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 spaceof partial-waves and hadron-hadron scattering channelsto some finite number, rendering D ( E cm ) an n × n matrixwhose determinant can be decomposed as the product ofits n eigenvalues λ p ( E cm ),det [ D ( E cm )] = n (cid:89) p =1 λ p ( E cm ) , where D ( E cm ) v ( p ) ( E cm ) = λ p ( E cm ) v ( p ) ( E cm ) , defines the corresponding eigenvectors . Solutions ofdet [ D ( E cm )] = 0 occur when at least one eigenvaluetakes a zero value, so our approach will be to find eacheigenvalue as a function of energy and then performroot-finding on each independently. In order to do this,the eigenvalues must be reliably matched at each step inenergy over the entire range considered, otherwise ‘jumps’due to mismatching could be mistaken for true zeroes andvice-versa, true zeroes could be missed.For the form of the quantization condition given inEqn. 1, at least one of the eigenvalues diverges at everynon-interacting energy owing to a pole in M , originat-ing in the L¨uscher zeta-functions . Such a divergencecan present a problem when attempting to match theeigenvalues or eigenvectors at energies close to these di-vergences – the assumption that a small change in energyleads to a small perturbation of D ( E cm ) is invalid nearthe non-interacting energies as the matrix is unbounded. A related approach is proposed in the context of the three-bodyquantization condition in Ref. [57]. See Eq. A2 in Appendix A.
In those scattering processes where interactions are rel-atively weak, it is very common that solutions of thequantization condition lie in a neighbourhood of thesenon-interacting energies, so removing these divergences isa priority.Fortunately it is quite straightforward to transform D ( E cm ) in such a way as to remove the non-interactingdivergences while leaving the zeroes of det (cid:2) D ( E cm ) (cid:3) un-changed. For example, D = + i ρ · t · ( + i M )= ρ · (cid:0) + S · V (cid:1) · (cid:0) − i M (cid:1) · ρ − , (2)where V ( E cm , L ) = ( + i M )( − i M ) − . The diver-gences at non-interacting energies have been explicitlyfactorised off in ( − i M ), leaving V which is free of thesedivergences. It follows that the zeroes of det (cid:2) D ( E cm ) (cid:3) are also zeroes of det (cid:2) D V ( E cm ) (cid:3) , where D V ( E cm ) = + S ( E cm ) · V ( E cm , L ) . (3)A closely related alternative isdet [ D U ( E cm )] = det [ S ( E cm ) + U ( E cm , L )] = 0 , (4)where U = ( − i M )( + i M ) − = V − , which is alsofinite at non-interacting energies. Note that like M , both V and U are diagonal in hadron-hadron channel butdense in partial-wave space.The properties of the matrices S , U and V dependupon energy, and in particular whether any channel iskinematically closed at any given energy. We begin byconsidering the case of energies high enough that allincluded channels are kinematically open. A. With all channels open
Consider a scattering system with m hadron-hadron channels which open at thresholds, E (1)thr. < E (2)thr. < · · · < E ( m )thr. . For each hadron-hadronchannel a , there may be multiple partial-waves underconsideration. Suppose channel a has n a partial waves,then (cid:80) m a =1 n a = n , the dimension of the matrix D ( E cm ). For energies E cm (cid:62) E ( m )thr. , the finite-volume matrix M ( E cm , L ) is hermitian , and then since i M is skew-hermitian, it has purely imaginary eigenvalues meaningthat ( − i M ) − exists for all such energies. It immedi-ately follows that V is a unitary matrix for all E cm (cid:62) E ( m )thr. .As S is also unitary in this energy region by conserva-tion of probability, commonly referred to as unitarity , D V ( E cm ) is diagonalisable with n linearly independent Which appeared originally in the elastic case in Ref. [30]. Some partial waves may appear embedded more than once. See Appendix A. orthonormal eigenvectors and bounded eigenvalues of theform λ p ( E cm ) = 2 cos (cid:0) θ p ( E cm ) (cid:1) exp (cid:0) i θ p ( E cm ) (cid:1) , (5)where the impact of unitarity is to ensure that the complexeigenvalues are each controlled by a single real parameter, θ p .Since D V ( E cm ) is continuous in E cm , a small pertur-bation in energy, E (cid:48) cm = E cm + ∆ E cm , which we will referto as an iteration in energy, results in a small perturba-tion of the matrix, D V ( E (cid:48) cm ) = D V ( E cm ) + ∆ D V ( E cm ),where || ∆ D V ( E cm ) || = O (∆ E cm ). It immediately fol-lows from standard matrix perturbation theory thatthe eigenvalues and eigenvectors at the next iteration, λ (cid:48) p ≡ λ p ( E (cid:48) cm ), v ( p ) (cid:48) ≡ v ( p ) ( E (cid:48) cm ) are approximated interms of the eigenvalues and eigenvectors at the currentiteration, λ p ≡ λ p ( E cm ), v ( p ) ≡ v ( p ) ( E cm ), by λ (cid:48) p = λ p + v ( p ) † · ∆ D V · v ( p ) + O (cid:0) || ∆ D V || (cid:1) v ( p ) (cid:48) = v ( p ) + (cid:88) k (cid:54) = p v ( k ) † · ∆ D V · v ( p ) λ p − λ k v ( k ) + O (cid:0) || ∆ D V || (cid:1) , so that for the case of non-degenerate eigenvalues, theeigenvectors are indeed small perturbations for ∆ E cm sufficiently small. The inner product between eigenvectorson consecutive energy iterations is given by v ( p ) (cid:48)† · v ( q ) = δ pq + ∆ pq , (6)where ∆ pq = v ( q ) † · ∆ D V · v ( p ) λ p − λ q + O (cid:0) || ∆ D V || (cid:1) , evaluated at E cm will be (cid:28) E cm . It follows that eigenvectors (and their corre-sponding eigenvalues) can be matched between energyiterations by examining Eqn. 6, with the order of theeigenvectors at E (cid:48) cm = E cm + ∆ E cm taken so as to max-imise the sum of the norm-squared of the inner productswith the eigenvectors at E cm .In the case that a number of eigenvaluesare degenerate or very nearly degenerate, i.e.when λ p − λ k = O (cid:0) || ∆ D V || (cid:1) for at least one k (cid:54) = p ,matching the corresponding eigenvectors between energyiterations is potentially ambiguous. As discussed inAppendix B, in these cases it is common to observean avoided crossing , where the eigenvalues repel eachother rather than cross, but crossings are possible whenchannels or partial-waves are decoupled. For an energy E cm where the eigenvalues are degenerate and cross, theeigenvectors span the corresponding eigenspace and thereis an arbitrariness in their direction. At this energy, it isnot necessary to match the eigenvectors to the previousiteration, as we are interested in the eigenvalues, whichare equal. At the next energy iteration, E cm + ∆ E cm ,where the eigenvalues are now distinct, rather thanmatching to the previous set of eigenvectors at E cm , wecan 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 ascattering system with a single hadron-hadron channel in an arbitrary number of partial-waves. For a particular energ,y ( E (1) cm ) ,below threshold, det[ − i M (( E (1) cm ) , L )] = 0, so an energy region is defined in which det (cid:2) S + U (cid:3) is considered rather thandet (cid:2) + S · V (cid:3) . B. With at least one closed channel
One might expect that the finite-volume spectrum atenergies below the kinematic threshold for some channel a to be insensitive to amplitudes describing channel a –after all, there are an infinite number of channels whosethresholds lie higher than any given energy, and we cannotever possibly know the amplitude in all of them. Infact, the finite-volume spectrum obtained from Eqn. 1 isonly sensitive to closed channels in a very small energyregion below the threshold, with the effect of the channelreducing exponentially (as e − κL where κ is the magnitudeof the imaginary scattering momentum). Nevertheless itis important to account for these closed-channel effectsin practical calculation, particularly when a resonanceappears near a threshold.At energies below the highest threshold under consid-eration, E cm (cid:54) E ( m )thr. , the finite-volume matrix M is ingeneral no longer hermitian and subsequently V is nolonger unitary. ( − i M ) is not prevented from havingzero eigenvalues and correspondingly V may diverge atsome energies.Recalling that M is block-diagonal in channel space,we can solve for roots of det[ − i M ] = 0 in each channelseparately. At decreasing energies below the thresholdfor channel a , M aa tends to i exponentially [58] and sosufficiently far below threshold, det[ − i M aa ] → n a . Foreach channel a root can only be found below threshold,as above threshold we recall M aa is hermitian and nosolutions to det[ − i M aa ] = 0 exist. In order to ensurethat zeroes of multiplicity higher than one are found, itis good practice to make use of an eigen-decompositionof M aa .A practical approach to dealing with these isolated en-ergies at which V diverges, is to eigen-decompose insteadthe matrix D U in a small region around the divergentenergy, since U is typically well defined at those energies.If it should happen by coincidence that det[ + i M aa ]has a zero very close to where det[ − i M aa ] has a zero,then we could eigen-decompose the original matrix D .We stress that in all practical applications, this finelytuned scenario seldom occurs, and can be avoided by con-sidering sufficiently small energy regions. Fig. 1 providesa pictorial representation of the approach.Finally, there can be divergences in the determinantcondition below the lowest threshold, where no channel is kinematically accessible, when the scattering systemfeatures stable bound-states which manifest as simplepole singularities in S ( E cm ). The divergence appears atthe volume-independent bound-state energy, which differsfrom a zero corresponding to a finite-volume energy levelonly 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 thereis 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-diagonalwith respect to hadron-hadron scattering channels, beingdense in the partial-wave space only. The threshold di-vergence in M is such that at threshold V and U arediagonal and equal to − .In general S is block-diagonal in partial-wave space, butdense in hadron-hadron channels, however at the open-ing of a new threshold, say for channel a , the vanishingof the phase-space for this channel means that for anypartial-wave, S aa ( E ( a )thr. ) = 1 and S ab ( E ( a )thr. ) = 0 for b (cid:54) = a ,while elements S bb ( E ( a )thr. ) are unconstrained. It followsthat at each threshold there will be at least one zeroeigenvalue of D V or D U , but these zeroes should not bemistakenly considered to be actual finite-volume energylevels – as they appear at exactly the known thresholdenergies, there should be no risk of this misidentification. D. Summary of the approach
Standard eigen-decomposition routines can be appliedto D V , D U or D (whichever is appropriate given the dis-cussion in Section. II B) on a discrete sampling of energies,and the eigenvalues λ p ( E cm ) matched between energy it-erations using similarity of the corresponding eigenvectorsestablished through magnitude of inner products. Eacheigenvalue can then be considered separately, and typicallywe find that root-finding on any one of these eigenvaluesis no more difficult than in the case of elastic scattering.Energy regions in which root-finding is performed areseparated by hadron-hadron thresholds and by regionsaround divergences in V , but as these regions are disjoint,the complete set of solutions is the combination of thesets of solutions of each eigenvalue in each energy region.It is therefore not necessary to match eigenvalues acrossboundaries.The combination of the sets of solutions in each eigen-value gives the full set of roots of det (cid:2) D ( E cm ) (cid:3) = 0, andof course with the eigenvalue roots in hand it is trivial tocheck that they are in fact roots of det (cid:2) D ( E cm ) (cid:3) . III. APPLICATION: VECTOR-VECTORSCATTERING
Vector-vector scattering provides a useful demonstra-tion of the eigen-decomposition approach described in theprevious section, as even in the simple case of total mo-mentum zero we can have a situation in which solutionsof high multiplicity appear in the non-interacting limit,while in the interacting case the corresponding zeroes maybe very closely spaced.We will illustrate our approach using a simple toymodel of vector-vector scattering in which we have asingle hadron-hadron channel comprising two identicalvector mesons. Considering the finite-volume spectrumin the [000] E + irrep, we note that Bose symmetry con-straints limit the number of contributing partial-waves,and we summarize those that can be non-zero for J (cid:54) (cid:96) = 2. As indicated in Table I, this resultsin the contribution of four partial waves, S , D , D and D .The vector particles are taken to have mass m = 0 . L = 70. In this case the non-interacting energies and the corresponding multiplicitiesare presented in Table II, which indicates that we expectfor a weakly interacting scattering system to find one levelclose to threshold, three closely-spaced levels somewhathigher in energy and four more levels a little higher still.Since there may be finite-volume energy levels belowthe opening of our threshold, we should first establish ifthere are energies around which we should use D U rather J P S +1 (cid:96) J + 5 S , D , D , G + 5 D , G , G , I TABLE I. Subduction of partial-waves, S +1 (cid:96) J , for J (cid:54) E + irrep for two identical vector mesons. The rowsdenote the partial-wave content for a given J P , with multiple S +1 (cid:96) J entries indicating partial-waves which mix dynamically.Partial-waves with (cid:96) (cid:62) (cid:0) L π (cid:1) (cid:12)(cid:12) (cid:126)p (cid:12)(cid:12) mult. E n.i . . . E n.i , for the lowest allowedback-to-back momenta, along with the corresponding multi-plicities that appear in [000] E + for a pair of vector particlesof mass m = 0 . L = 70. than our preferred choice D V . This can be done withoutknowledge of the scattering amplitude, by looking forroots of det[ − i M ] = 0 below threshold. Fig. 2 showsthe eigenvalues of − i M and their roots. Note that eacheigenvalue diverges at the threshold, and asymptoticallyapproaches a value of 2 (as M → i ) as energies go furtherbelow threshold. Two roots are found: a double root at( E cm ) = 0 . E cm ) = 0 . around each of these values ischosen to be treated with D U in place of D V .To describe scattering in this system of four partial-waves we use a K -matrix parameterisation which automat-ically implements the required unitarity of the S -matrix.The K -matrix, K ( s ), where s = E cm , is a real symmetricmatrix, related to the t -matrix via (cid:2) t − ( s ) (cid:3) S(cid:96)J,S (cid:48) (cid:96) (cid:48) J = 1(2 k cm ) (cid:96) (cid:2) K − ( s ) (cid:3) S(cid:96)J,S (cid:48) (cid:96) (cid:48) J k cm ) (cid:96) (cid:48) + δ SS (cid:48) δ (cid:96)(cid:96) (cid:48) I ( s ) . (7)Unitarity is guaranteed if Im I ( s ) = − ρ ( s ) for energiesabove threshold and Im I ( s ) = 0 below. We make useof the Chew-Mandelstam prescription [60] which definesRe I ( s ) through a dispersive integral featuring ρ ( s ), andchoose to subtract so that Re I = 0 at threshold. FIG. 2. Eigenvalues of − i M ( E cm , L ) for E cm (cid:54) E (1)thr. . Theimaginary part of each eigenvalue is identically zero belowthreshold. For the m, L combination used here, an energy region 1 . × − either side of the root is quite reasonable. FIG. 3. Phase-shifts and mixing-angles for the toy modelamplitude given in Eq. 8 following the n -channel Stapp pa-rameterisation described in Ref. [24]. To implement weak scattering we use a constant K -matrix parameterisation featuring small values, K ( s ) = S D D D S −
10 10 D − D − D , (8)where the block diagonalization in J is made explicit.The somewhat larger values in the D -wave amplitudesare compensated by the threshold factors k (cid:96) cm in Eq. 7.The resulting phase-shifts and mixing-angles for this toyamplitude are presented in Fig. 3 using the n -channelgeneralization of the Stapp parameterisation described inRef. [24]. Note that all amplitudes are modest by design,with the S -wave an order of magnitude larger than anyof the D -wave amplitudes.Following the procedure set out in Sec. II, we seek tofind solutions of det[ D ] = 0. The top panel of Fig. 4 showsdet[ D ] as defined in Eqn. 1 as a function of energy, wheredivergences at each of the non-interacting energies can beclearly seen, along with the near-degenerate roots locatedin very close proximity – as anticipated from our weaklyinteracting amplitudes. Finding the zeroes directly fromthis determinant requires ultra-fine sampling in energy,and roots are easily missed in automated root-findingalgorithms. Furthermore, the toy amplitudes could beeasily modified to make the near degenerate roots, shownin the zoomed regions, exactly degenerate, in which caseno amount of improvement in sampling resolution wouldfind the appropriate number of roots.The second panel in Fig. 4 shows det[ D V ] as a functionof energy. As discussed in Sec. II, the divergences that appear at the non-interacting energies have been removedand above threshold the determinant is bounded. In aneighbourhood around the near degenerate roots, thedeterminant is more clearly seen to be locally quadraticand cubic respectively.The remaining lower panels in Fig. 4 show the foureigenvalues 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 , − , V is singular as expected– see Sec. II B. The multiplicities of these divergences alsoagree, for example as shown in Fig. 2 the determinantdet[ − i M ] has a double root at ( E cm ) , which manifestsin the eigenvalues λ and λ , and a single root at ( E cm ) ,which appears in λ . At threshold, each eigenvalue isidentically zero as proven in Sec. II C and these thresholdzeroes are removed from the set of solutions. It is clearfrom counting the number of roots that we have agreementwith 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 moredifficult than solving for roots in elastic scattering andthis is illustrated by the simplicity of the zero crossings inthe eigenvalues. The found roots are recorded in Table III. λ λ λ λ D ] = det[ D V ] = 0 and eigen-values, λ p , of D V in which they appear. Values groupedaccording to the non-interacting multiplicities provided inTable II. IV. APPLICATION: + , − SCATTERING
Previous lattice QCD calculations have determined thelowest lying isospin=1 J P C = 1 −− , + − vector [8] andaxial-vector [24] resonances, the ρ and b , in studies onthree lattice volumes with a pion mass ∼
391 MeV. Inthe determination of the b resonance, energy levels usedto constrain the scattering amplitudes were obtained inthe [000] T +1 and ( (cid:126)P (cid:54) = (cid:126) A irreps, as these irreps receiveno contribution from J P = 1 − scattering amplitudes.This avoids the complication of disentangling the low-lying, resonating 1 − channels that mix due to the reducedsymmetry of the finite-volume. The eigen-decompositionmethod we have introduced in this manuscript enables us Details of the lattices, and the correlator construction technologycan be found in Refs. [61, 62].
FIG. 4. Quantization conditions for the toy model of vector-vector scattering described in the text. Kinematic threshold at E cm = 1 . (cid:2) D ] as defined in Eqn. 1. Second panel: det (cid:2) D V ] as defined in Eqn. 3. Bottom four panels:Eigenvalues of det (cid:2) D V ] – divergences below threshold at energies presented in Fig. 2. to lift such practical restrictions and perform a combinedscattering analysis of the 1 −− and 1 + − sectors. Using amuch larger set of irreps, including now those in whichboth parities feature, allows us to significantly increasethe number of energy levels constraining the scatteringamplitudes, with a total of 144 energy levels being avail-able.Several hadron-hadron channels in various partial-wavesfeature in these scattering systems. For the J P = 1 − sec-tor, in addition to ππ { P } considered in Ref. [8] (whichrestricted itself to elastic energies), we include KK { P } and πω { P } which are kinematically open in the en-ergy region where where the b resonance appears. For J P = 1 + , we consider πω { S } , πω { D } and πφ { S } as in Ref. [24]. As discussed in great detail in Ref. [24],three-body thresholds, ππη and πKK , open below the4 π -threshold in the region where the b resonates and weensure that in each irrep, we keep below the lowest E non-interacting energy .Other partial-waves also feature due to the reducedsymmetry of the finite-volume. For ππ and KK withisospin=1, a F wave can appear (see Table III of Ref. [8])but will be suppressed by the high value of (cid:96) . For vector-pseudoscalar scattering, the nonzero intrinsic spin givesrise to a triplet of partial-waves for each (cid:96) (cid:62)
1, e.g. P , P and P for (cid:96) = 1 (see Table 1 of Ref. [55] for irrepsat rest, and Tables 5-7 for irreps in flight). Followingthe analysis in Ref. [24], we include only πω { P } , with J P = 1 − which has not been previously shown to becompatible with zero. The resulting set of partial wavesto be considered is presented in Table IV. − ππ { P } KK { P } πω { P } + πω { S } πω { D } πφ { S } TABLE IV. Partial-wave basis to describe finite-volume spec-trum.
As in the case of the toy amplitudes in Sec. III, singu-larities of V must be determined in order to solve theappropriate form of quantisation condition in each energyregion. For reference, tables VI – VIII in Appendix Creport the singularities of V on each volume and latticeirrep.In this case we are seeking to describe an explicit spec-trum of lattice QCD energy levels (with uncertaintiesand correlations) using a parameterised amplitude in a χ minimization . In the approach followed by hadspec ,typically a range of parameterisations would be consid-ered, but in this case, where we seek only to demonstratethe effectiveness of the eigen-decomposition approach, we As introduced in Ref. [24], E non-interacting energies arefinite-volume energies calculated in a two-meson subsystem thatare then combined with the third meson, where no interactionsare assumed between the two-meson subsystem and third meson. The χ is defined explicitly in Eqn. 9 of Ref. [8]. use just one , a K -matrix parameterisation of the ‘poleplus constant’ form for each J P , K ( s ) = (cid:20) K − ( s ) K + ( s ) (cid:21) , K − ( s ) = 1 m − − s g ππ { P } + γ ππ { P } γ KK { P }
00 0 γ πω { P } , K + ( s ) = 1 m − s g πω { S } g πω { S } g πω { D } g πω { S } g πω { D } g πω { D }
00 0 0 + γ πω { S } γ πφ { S } . (9)We use the Chew-Mandelstam prescription for I ( s ), choos-ing to subtract at the pole masses so that Re I aa ( m − ) = 0where a runs over channels with J P = 1 − and similarlyRe I aa ( m ) = 0 for channels with J P = 1 + . The use ofonly a single pole in each J P is motivated by the expecta-tion of a narrow ρ resonance, a narrow b resonance, andno other resonances until much higher energies [63, 64].The best description of the lattice spectra in ir-reps [000] T ± , [001] ( A , A , E ), [011] ( A , A , B , B ),[111] ( A , A , E ) and [002] ( A , A ), varying the ten pa-rameters in Eqn. 9, is found with a quite reasonable χ /N dof = 229 ./ (144 −
10) = 1 .
71, after ∼
400 iterationsof the MINUIT minimisation routine [65].Fig. 5 shows as orange curves the finite-volume spec-trum corresponding to the best-fit amplitude of Eqn. 9 forthose irreps in which both J P = 1 − and 1 + are subduced.Also shown are the relevant non-interacting energy curvesand the lattice QCD energy levels which were fitted toobtain the best-fit parameter values.In Fig. 6 we illustrate the solution of the quantiza-tion condition in the case of the [111] E irrep on the L/a s = 24 lattice using the central values of the best fitparameters. A single embedding of each partial-wave ineach hadron-hadron channel recorded in Table IV fea-tures in this irrep, and thus D V is a 6 × V (see Table VIII) are manifested inthe eigenvalue divergences, i.e. three singularities withmultiplicity one. The determinant of D V is indeed boundabove the πφ threshold and moreover, as no singularityof D V appears above πω threshold in this case, D V is and in addition we do not consider the effects of varying the stablehadron masses, or the anisotropies within their uncertainties. { }
16 20 24 16 20 24 { }
16 20 24 { }{ }{ } FIG. 5. For those irreps which feature both J P = 1 + and1 − , the lattice QCD energy spectrum (black points, or grey ifpoints not used in the fit), non-interacting curves (dashed lines)and in orange the finite-volume spectrum corresponding to theamplitude in Eqn. 9 using the best-fit parameters, includingthe propagation of their correlated errors. also bound between πω and πφ threshold. For this par-ticular parameterisation, the S -matrix is decoupled inhadron channel (see Eq. 9). As V is diagonal in hadronchannel, D V is block-diagonal and the eigenvectors aretherefore orthogonal in channel space. This is manifestin the corresponding eigenvalues and we can make thefollowing associations: λ corresponds to ππ { P } , λ to KK { P } , λ , λ , λ to πω { S , P , D } and λ to πφ { S } . Of course this is not a general feature of theeigen-decomposition method as the S -matrix is generallydense in channel space. The roots and the correspondingeigenvalues of D V in which they appear are recorded inTable V. λ λ λ λ λ λ D ] = det[ D V ] = 0 in the [111] E irrep on the L/a s = 24 lattice for the amplitude in Eqn. 9listed by the eigenvalue of D V in which they appear. Greyentries denote roots that were found at energies higher thanthose we are considering. For this parameterisation, we find the ππ { P } am-plitude appears to agree very well with the amplitudesdetermined in the previous calculation [8], with the energydependence above the resonance being quite modest. Theamplitude features a complex conjugate pair of resonance poles, which we interpret as the ρ , at, a t E cm = 0 . ± i . , (10)which is in quite reasonable agreement with the valuefound previously using a more limited set of irreps. The KK { P } and πω { P } amplitudes were found to beweak over the entire energy range considered, with phase-shifts not exceeding 1 ◦ and 17 ◦ respectively, as anticipatedowing 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 agreevery well with those found in the previous calculation [24].The amplitude features a complex conjugate pair of reso-nance poles, which we interpret as the b , at, a t E cm = 0 . ± i . . (11)in close agreement with the pole determined in Ref. [24]. V. SUMMARY
We have presented an approach for efficiently solvingthe L¨uscher quantization condition which determines thediscrete spectrum of states in a finite-volume given ascattering system. The technique makes use of the eigen-value decomposition of the matrix appearing under thedeterminant in the quantization condition, ultimatelyreducing the complexity of the numerical problem forcoupled-channels and coupled partial-waves down to oneof root-finding in a number of independent continuousfunctions which typically each feature far less rapid varia-tion 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 vectormesons shows the ability of the approach to determinethe very dense set of zeroes lying close to non-interactinglevels of high multiplicity. In the second, a set of latticeQCD calculated energy levels on three volumes were usedto constrain coupled-channel scattering with J P = 1 − and 1 + . This provides a stress-test of the approach, as thesolving of the quantization condition needs be successfulon multiple irreps in multiple volumes, at many hundredsof iterations of the scattering amplitude parameter values,varied under a χ minimization controlled by a minimiza-tion routine. The approach was found to be both reliableand computationally fast.Contemporary examples where application of this ap-proach may prove to make possible previously impracticalcalculations include the determination of the resonancespectrum in charmonium where a decay product vectormeson D ∗ is well approximated as a stable hadron, andwhere the D ¯ D , D ∗ ¯ D and D ∗ ¯ D ∗ thresholds lie very closeto each other. In-flight irreps will receive contributionsfrom many partial-waves, and there is the expectationof resonances in all of J P C = 1 −− , −− , −− with rathersimilar masses. Nucleon resonances provide another ap-plication where a relatively large number of scattering0 FIG. 6. Quantization condition for the [111] E irrep on the L/a s = 24 lattice described in the text. Kinematic thresholds for ππ , KK , πω and πφ at a t E cm = 0 . , . , . . V recorded in Table VIII in Appendix C. Red, blue curves indicate real,imaginary parts respectively. Top panel: det (cid:2) D V ] as defined in Eqn. 3. Bottom six panels: Eigenvalues of det (cid:2) D V ]. Blackpoints at the foot of the figure denote the lattice computed energies transcribed from Figure 5. Orange points are the roots ofdet (cid:2) D V ] with statistical uncertainties reflecting the uncertainties on the scattering parameters. ACKNOWLEDGMENTS
We thank our colleagues within the Hadron SpectrumCollaboration. AJW is supported by the U.K. Scienceand Technology Facilities Council (STFC) [grant num-ber ST/P000681/1]. JJD acknowledges support from theU.S. Department of Energy contract de-sc0018416, andcontract DE-AC05-06OR23177, under which JeffersonScience Associates, LLC, manages and operates JeffersonLab. DJW acknowledges support from a Royal Soci-ety University Research Fellowship and partial supportfrom the U.K. Science and Technology Facilities Council(STFC) [grant number ST/P000681/1].The software codes
Chroma [66] and
QUDA
Appendix A: Hermiticity of M above threshold Bounding of the eigenvalues and orthogonality of theeigenvectors of D V followed from hermiticity of M abovethreshold. Here we show this hermiticity, starting withthe general expression for matrix elements M in thebasis of intrinsic spin S , orbital angular momentum (cid:96) ,total angular momentum J , azimuthal component of totalangular momentum m and channel index a , M (cid:96)SJma, (cid:96) (cid:48) S (cid:48) J (cid:48) m (cid:48) b = δ ab δ SS (cid:48) (cid:88) m (cid:96) ,m (cid:48) (cid:96) m S (cid:104) (cid:96)m (cid:96) ; Sm S | Jm (cid:105) (cid:104) (cid:96) (cid:48) m (cid:48) (cid:96) ; Sm S | J (cid:48) m (cid:48) (cid:105)× (cid:88) ¯ (cid:96), ¯ m (cid:96) (4 π ) / k ¯ (cid:96) +1 cm c (cid:126)n ¯ (cid:96), ¯ m (cid:96) ( k cm ; L ) (cid:90) d Ω Y ∗ (cid:96)m (cid:96) Y ∗ ¯ (cid:96) ¯ m (cid:96) Y (cid:96) (cid:48) m (cid:48) (cid:96) . (A1)The spin-orbit coupling is expressed through the real SU(2) Clebsch-Gordan coefficients (cid:104) (cid:96)m (cid:96) ; Sm S | Jm (cid:105) and thevolume dependence is encoded in the functions c (cid:126)n(cid:96),m (cid:96) ( k cm ; L ), defined as c (cid:126)n(cid:96),m (cid:96) ( k cm ; L ) = √ πγL (cid:18) πL (cid:19) (cid:96) − Z (cid:126)n(cid:96),m (cid:96) (cid:20) (cid:18) k cm L π (cid:19) (cid:21) , Z (cid:126)n(cid:96),m (cid:96) [ s ; x ] = (cid:88) (cid:126)r ∈P (cid:126)n | (cid:126)r | (cid:96) Y (cid:96)m (cid:96) ( (cid:126)r )( | (cid:126)r | − x ) s , (A2)where (cid:126)n = L π (cid:126)P , and k cm is the cm -frame momentum corresponding to E cm . Definitions of the boost-factor, γ , and theconstruction of the set of real vectors to be summed over, (cid:126)r ∈ P (cid:126)n , can be found in Ref. [43]. Using Y ∗ (cid:96),m = ( − m Y (cid:96), − m ,it is easy to show that c (cid:126)n ∗ (cid:96),m (cid:96) = c (cid:126)n(cid:96), − m (cid:96) . The presence of the factor k − (¯ (cid:96) +1) cm in Eqn. A1 distinguishes the above and belowthreshold 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.2The integral over the product of three spherical harmonics can be expressed in terms of Clebsch-Gordan coefficients, (cid:90) d Ω Y ∗ (cid:96)m (cid:96) Y ∗ ¯ (cid:96) ¯ m (cid:96) Y (cid:96) (cid:48) m (cid:48) (cid:96) = (cid:115) (2 (cid:96) + 1)(2¯ (cid:96) + 1)4 π (2 (cid:96) (cid:48) + 1) (cid:104) (cid:96)m (cid:96) ; ¯ (cid:96) ¯ m (cid:96) | (cid:96) (cid:48) m (cid:48) (cid:96) (cid:105) (cid:104) (cid:96)
0; ¯ (cid:96) | (cid:96) (cid:48) (cid:105) , and thus the integral is manifestly real. We leave it in the integral form as it more directly exposes the relevantsymmetries.Under complex conjugation of Eqn. A1 we can replace the resulting c (cid:126)n ∗ (cid:96),m (cid:96) with c (cid:126)n(cid:96), − m (cid:96) and Y ¯ (cid:96), ¯ m (cid:96) with ( − ¯ m (cid:96) Y ∗ ¯ (cid:96), − ¯ m (cid:96) ,and then noting that the sum over ¯ m (cid:96) is equivalent to a sum over − ¯ m (cid:96) obtain in the above threshold case where k cm is real, M ∗ (cid:96)SJma, (cid:96) (cid:48) S (cid:48) J (cid:48) m (cid:48) b = δ ab δ SS (cid:48) (cid:88) m (cid:96) ,m (cid:48) (cid:96) m S (cid:104) (cid:96) (cid:48) m (cid:48) (cid:96) ; Sm S | J (cid:48) m (cid:48) (cid:105) (cid:104) (cid:96)m (cid:96) ; Sm S | Jm (cid:105)× (cid:88) ¯ (cid:96), ¯ m (cid:96) (4 π ) / k ¯ (cid:96) +1 cm c (cid:126)n ¯ (cid:96), ¯ m (cid:96) ( k cm ; L ) (cid:90) d Ω Y ∗ (cid:96) (cid:48) m (cid:48) (cid:96) Y ∗ ¯ (cid:96) ¯ m (cid:96) Y (cid:96)m (cid:96) , (A3)from which the symmetry under exchange of the totality of matrix indices is evident under a relabelling m (cid:96) ↔ m (cid:48) (cid:96) ,so that M ∗ (cid:96)SJma, (cid:96) (cid:48) S (cid:48) J (cid:48) m (cid:48) b = M (cid:96) (cid:48) S (cid:48) J (cid:48) m (cid:48) b, (cid:96)SJma , and the hermiticity of M above threshold is established. Subductioninto irreducible representations of the relevant cubic symmetry groups does not affect this result as the subductionmatrices are unitary. Appendix B: Avoided Eigenvalue Crossings
Eigenvalues of D V = + S · V are observed to exhibita behavior where they avoid crossing each other as E cm is varied. Examples from the toy model considered inSection III are shown in Figure 7. This behaviour is infact a general property that follows from the unitarity ofthe matrix S · V above threshold.Defining unitary matrices at neighbour-ing energy points, W = S ( E cm ) · V ( E cm ) and W (cid:48) = S ( E cm + ∆ E cm ) · V ( E cm + ∆ E cm ), we canfind hermitian matrices, H , H (cid:48) such that W = exp i HW (cid:48) = exp i H (cid:48) , and clearly the eigenvalues and eigenvectors of these ma-trices are trivially related, Hv ( p ) = θ p v ( p ) W v ( p ) = e iθ p v ( p ) , with θ p being real, and simply related to the complexeigenvalues λ p of D V as given in Eqn. 5. The avoided levelcrossings that manifest in both the real and imaginaryparts between some λ p ’s therefore appear also betweenthe corresponding θ p ’s.By choosing to examine H we reduce the problem toone familiar in quantum mechanics: as H and H (cid:48) − H are hermitian this is the case considered in (degener-ate) perturbation theory, where ‘avoided level cross-ings’ are a generic feature whenever the perturbation( H (cid:48) − H ) connects the (degenerate) eigenstates. A classicillustrative example is the two-state case where in theeigenbasis of H , H (cid:48) = (cid:18) θ θ (cid:19) + (cid:18) P P P ∗ P (cid:19) , and the corresponding exact eigenvalues of H (cid:48) are θ (cid:48) = a + b , θ (cid:48) = a − b , where a = (cid:0) θ + θ + P + P (cid:1) b = (cid:113)(cid:0) θ − θ + P − P (cid:1) + 4 (cid:12)(cid:12) P (cid:12)(cid:12) . It follows that unless both terms under the square rootare zero, (cid:12)(cid:12) θ (cid:48) − θ (cid:48) (cid:12)(cid:12) >
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 calculationpresented in Sec. IV. We record singularities for irrepsfeaturing subductions of J P = 1 − and not J P = 1 + inTable VI, irreps featuring J P = 1 + and not J P = 1 − inTable VII, and irreps featuring both in Table VIII. Forclarity, the singularities arising in different hadron-hadronchannels (‘ a ’) is made explicit in each of the tables. [1] S. R. Sharpe, R. Gupta, and G. W. Kilcup, Nucl. Phys. B383 , 309 (1992). [2] S. R. Beane, T. C. Luu, K. Orginos, A. Parreno, M. J.Savage, A. Torok, and A. Walker-Loud, Phys. Rev.
D77 , FIG. 7. Eigenvalues, λ , . . . , λ of D V ( E cm ). The real and imaginary components of each eigenvalue are shown in red and bluerespectively and the different dashed lines correspond to the different eigenvalues. Highlighted are examples of avoided levelcrossings between the real and imaginary parts of different eigenvalues as discussed in the text. L/a s ‘ a ’ [000] T − [001] A [011] A [111] A [002] A ππ KK πω ππ KK πω ππ KK πω E cm of the matrix V ( E cm ) as de-scribed in the text. L/a s ‘ a ’ [000] T +1 [001] A [011] A [111] A [002] A πω { } πφ πω { } πφ πω { } πφ D78 , 014511 (2008),arXiv:0804.2941 [hep-lat].[4] X. Feng, K. Jansen, and D. B. Renner, Phys. Lett.
B684 ,268 (2010), arXiv:0909.3255 [hep-lat].
L/a s ‘ a ’ [001] E [011] B [011] B [111] E ππ KK πω πφ ππ KK πω πφ ππ KK πω πφ D83 , 071504(2011), arXiv:1011.6352 [hep-ph].[6] S. R. Beane, E. Chang, W. Detmold, H. W. Lin, T. C.Luu, K. Orginos, A. Parreno, M. J. Savage, A. Torok,and A. Walker-Loud (NPLQCD), Phys. Rev.
D85 , 034505(2012), arXiv:1107.5023 [hep-lat].[7] J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys.Rev.
D86 , 034031 (2012), arXiv:1203.6041 [hep-ph].[8] J. J. Dudek, R. G. Edwards, and C. E. Thomas (HadronSpectrum), Phys. Rev.
D87 , 034505 (2013), [Erratum:Phys. Rev.D90,no.9,099902(2014)], arXiv:1212.0830 [hep-ph].[9] Z. Bai et al. (RBC, UKQCD), Phys. Rev. Lett. ,212001 (2015), arXiv:1505.07863 [hep-lat]. [10] T. Blum et al. , Phys. Rev. D91 , 074502 (2015),arXiv:1502.00263 [hep-lat].[11] C. Helmes, C. Jost, B. Knippschild, C. Liu, J. Liu, L. Liu,C. Urbach, M. Ueding, Z. Wang, and M. Werner (ETM),JHEP , 109 (2015), arXiv:1506.00408 [hep-lat].[12] L. Liu et al. , Phys. Rev. D96 , 054516 (2017),arXiv:1612.02061 [hep-lat].[13] R. A. Briceno, J. J. Dudek, R. G. Edwards, andD. J. Wilson, Phys. Rev. Lett. , 022002 (2017),arXiv:1607.05900 [hep-ph].[14] J. Bulava, B. Fahy, B. Horz, K. J. Juge, C. Morn-ingstar, and C. H. Wong, Nucl. Phys.
B910 , 842 (2016),arXiv:1604.05593 [hep-lat].[15] C. Alexandrou, L. Leskovec, S. Meinel, J. Negele, S. Paul,M. Petschlies, A. Pochinsky, G. Rendon, and S. Syritsyn,Phys. Rev.
D96 , 034525 (2017), arXiv:1704.05439 [hep-lat].[16] C. Andersen, J. Bulava, B. H¨orz, and C. Morningstar,Nucl. Phys.
B939 , 145 (2019), arXiv:1808.05007 [hep-lat].[17] J. J. Dudek, R. G. Edwards, C. E. Thomas, and D. J.Wilson (Hadron Spectrum), Phys. Rev. Lett. , 182001(2014), arXiv:1406.4158 [hep-ph].[18] D. J. Wilson, J. J. Dudek, R. G. Edwards, and C. E.Thomas, Phys. Rev.
D91 , 054008 (2015), arXiv:1411.2004[hep-ph].[19] D. J. Wilson, R. A. Briceno, J. J. Dudek, R. G. Ed-wards, and C. E. Thomas, Phys. Rev.
D92 , 094502(2015), arXiv:1507.02599 [hep-ph].[20] J. J. Dudek, R. G. Edwards, and D. J. Wilson(Hadron Spectrum), Phys. Rev.
D93 , 094506 (2016),arXiv:1602.05122 [hep-ph].[21] G. K. C. Cheung, C. O’Hara, G. Moir, M. Peardon, S. M.Ryan, C. E. Thomas, and D. Tims (Hadron Spectrum),JHEP , 089 (2016), arXiv:1610.01073 [hep-lat].[22] R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J.Wilson, Phys. Rev. D97 , 054513 (2018), arXiv:1708.06667[hep-lat].[23] C. B. Lang, L. Leskovec, D. Mohler, and S. Prelovsek,JHEP , 162 (2014), arXiv:1401.2088 [hep-lat].[24] A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards,and D. J. Wilson, (2019), arXiv:1904.04136 [hep-lat].[25] A. Torok, S. R. Beane, W. Detmold, T. C. Luu,K. Orginos, A. Parreno, M. J. Savage, and A. Walker-Loud, Phys. Rev. D81 , 074506 (2010), arXiv:0907.1913[hep-lat].[26] C. B. Lang and V. Verduci, Phys. Rev.
D87 , 054502(2013), arXiv:1212.5055 [hep-lat].[27] W. Detmold and A. Nicholson, Phys. Rev.
D93 , 114511(2016), arXiv:1511.02275 [hep-lat].[28] C. B. Lang, L. Leskovec, M. Padmanath, and S. Prelovsek,Phys. Rev.
D95 , 014510 (2017), arXiv:1610.01422 [hep-lat].[29] C. W. Andersen, J. Bulava, B. H¨orz, and C. Morningstar,Phys. Rev.
D97 , 014506 (2018), arXiv:1710.01557 [hep-lat].[30] M. Luscher, Nucl. Phys.
B354 , 531 (1991).[31] K. Rummukainen and S. A. Gottlieb, Nucl. Phys.
B450 ,397 (1995), arXiv:hep-lat/9503028 [hep-lat].[32] P. F. Bedaque, Phys. Lett.
B593 , 82 (2004), arXiv:nucl-th/0402051 [nucl-th].[33] S. He, X. Feng, and C. Liu, JHEP , 011 (2005),arXiv:hep-lat/0504019 [hep-lat].[34] C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl.Phys. B727 , 218 (2005), arXiv:hep-lat/0507006 [hep-lat]. [35] Z. Fu, Phys. Rev.
D85 , 014506 (2012), arXiv:1110.0319[hep-lat].[36] L. Leskovec and S. Prelovsek, Phys. Rev.
D85 , 114507(2012), arXiv:1202.2145 [hep-lat].[37] M. Gockeler, R. Horsley, M. Lage, U. G. Meissner, P. E. L.Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti,Phys. Rev.
D86 , 094513 (2012), arXiv:1206.4141 [hep-lat].[38] V. Bernard, M. Lage, U. G. Meissner, and A. Rusetsky,JHEP , 019 (2011), arXiv:1010.6018 [hep-lat].[39] M. Doring, U. G. Meissner, E. Oset, and A. Rusetsky,Eur. Phys. J. A48 , 114 (2012), arXiv:1205.4838 [hep-lat].[40] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D86 , 016007(2012), arXiv:1204.0826 [hep-lat].[41] R. A. Briceno and Z. Davoudi, Phys. Rev.
D88 , 094507(2013), arXiv:1204.1110 [hep-lat].[42] P. Guo, J. Dudek, R. Edwards, and A. P. Szczepaniak,Phys. Rev.
D88 , 014501 (2013), arXiv:1211.0929 [hep-lat].[43] R. A. Briceno, Phys. Rev.
D89 , 074507 (2014),arXiv:1401.3312 [hep-lat].[44] M. T. Hansen and S. R. Sharpe, Phys. Rev.
D90 , 116003(2014), arXiv:1408.5933 [hep-lat].[45] R. A. Brice˜no, M. T. Hansen, and S. R. Sharpe, Phys.Rev.
D95 , 074510 (2017), arXiv:1701.07465 [hep-lat].[46] H. W. Hammer, J. Y. Pang, and A. Rusetsky, JHEP ,115 (2017), arXiv:1707.02176 [hep-lat].[47] H.-W. Hammer, J.-Y. Pang, and A. Rusetsky, JHEP ,109 (2017), arXiv:1706.07700 [hep-lat].[48] M. Mai and M. D¨oring, Eur. Phys. J. A53 , 240 (2017),arXiv:1709.08222 [hep-lat].[49] M. D¨oring, H. W. Hammer, M. Mai, J. Y. Pang, A. Ruset-sky, and J. Wu, Phys. Rev.
D97 , 114508 (2018),arXiv:1802.03362 [hep-lat].[50] R. A. Brice˜no, M. T. Hansen, and S. R. Sharpe, Phys.Rev.
D99 , 014516 (2019), arXiv:1810.01429 [hep-lat].[51] M. T. Hansen and S. R. Sharpe, (2019), arXiv:1901.00483[hep-lat].[52] R. A. Brice˜no, M. T. Hansen, S. R. Sharpe, andA. P. Szczepaniak, Phys. Rev.
D100 , 054508 (2019),arXiv:1905.11188 [hep-lat].[53] F. Romero-L´opez, S. R. Sharpe, T. D. Blanton, R. A.Brice˜no, and M. T. Hansen, JHEP , 007 (2019),arXiv:1908.02411 [hep-lat].[54] C. E. Thomas, R. G. Edwards, and J. J. Dudek, Phys.Rev. D85 , 014507 (2012), arXiv:1107.1930 [hep-lat].[55] A. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards,and D. J. Wilson, JHEP , 043 (2018), arXiv:1802.05580[hep-lat].[56] C. Morningstar, J. Bulava, B. Singha, R. Brett, J. Fallica,A. Hanlon, and B. H¨orz, Nucl. Phys. B924 , 477 (2017),arXiv:1707.05817 [hep-lat].[57] T. D. Blanton, F. Romero-L´opez, and S. R. Sharpe,JHEP , 106 (2019), arXiv:1901.07095 [hep-lat].[58] R. A. Briceno, M. T. Hansen, and A. W. Jackura, (2019),arXiv:1909.10357 [hep-lat].[59] R. C. Johnson, Phys. Lett. , 147 (1982).[60] G. F. Chew and S. Mandelstam, Phys. Rev. , 467(1960).[61] H.-W. Lin et al. (Hadron Spectrum), Phys. Rev. D79 ,034502 (2009), arXiv:0810.3588 [hep-lat].[62] M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek,R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, andK. J. Juge (Hadron Spectrum), Phys. Rev.
D80 , 054506(2009), arXiv:0905.2160 [hep-lat]. [63] J. J. Dudek, R. G. Edwards, P. Guo, and C. E. Thomas(Hadron Spectrum), Phys. Rev. D88 , 094505 (2013),arXiv:1309.2608 [hep-lat].[64] J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G.Richards, and C. E. Thomas, Phys. Rev.
D82 , 034508(2010), arXiv:1004.4930 [hep-ph].[65] F. James and M. Roos, Comput. Phys. Commun. , 343(1975).[66] R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD), Lattice field theory. Proceedings, 22nd International Sym- posium, Lattice 2004, Batavia, USA, June 21-26, 2004 ,Nucl. Phys. Proc. Suppl. , 832 (2005), [,832(2004)],arXiv:hep-lat/0409003 [hep-lat].[67] M. A. Clark, R. Babich, K. Barros, R. C. Brower, andC. Rebbi, Comput. Phys. Commun. , 1517 (2010),arXiv:0911.3191 [hep-lat].[68] R. Babich, M. A. Clark, and B. Joo, in