Numerical study of the relativistic three-body quantization condition in the isotropic approximation

We present numerical results showing how our recently proposed relativistic three-particle quantization condition can be used in practice. Using the isotropic (generalized $s$-wave) approximation, and keeping only the leading terms in the effective range expansion, we show how the quantization condition can be solved numerically in a straightforward manner. In addition, we show how the integral equations that relate the intermediate three-particle infinite-volume scattering quantity, $\mathcal K_{\text{df},3}$, to the physical scattering amplitude can be solved at and below threshold. We test our methods by reproducing known analytic results for the $1/L$ expansion of the threshold state, the volume dependence of three-particle bound-state energies, and the Bethe-Salpeter wavefunctions for these bound states. We also find that certain values of $\mathcal K_{\text{df},3}$ lead to unphysical finite-volume energies, and give a preliminary analysis of these artifacts.


I. INTRODUCTION
Studies of few-hadron systems based on lattice quantum chromodynamics (LQCD) are advancing rapidly. Recent results highlighting this progress include the first study of multiple, strongly coupled scattering channels [1,2], the first determination of resonant electroweak amplitudes [3,4], and the first study of a meson-baryon scattering amplitude in a resonant channel [5]. Each of these calculations has been made possible by a series of theoretical developments, stemming from seminal work by Lüscher [6,7]. This formalism and its subsequent generalizations explain how the desired infinite-volume observables, namely scattering and transition amplitudes, can be obtained from the finite-volume correlation functions evaluated using numerical LQCD. We point the reader to Ref. [8] for a recent review on the topic.
Current theoretical work is focused on extending the finite-volume relations to extract observables with initial or final states composed of three or more hadrons. To this end, in a series of papers published in the past few years, we have derived a quantization condition that relates the finitevolume energies of states containing a three-particle component to infinite-volume, two-and three-particle scattering amplitudes [9][10][11]. 1 This quantization condition accounts for all power-law volume dependence while dropping dependence that falls exponentially with the box length, L. The formalism is relativistic and encompasses arbitrary interactions aside from two restrictions: (i) the particles must be spinless and identical, and (ii) the two-particle K matrix cannot have poles in the kinematical regime of interest. From our past experience in the two-body sector [22,23], we expect the former restriction to be straightforward to remove, and now understand how to remove the latter [24]. The relation to physical scattering amplitudes involves two steps. In the first, the quantization condition is used to determine an infinitevolume K matrix like quantity, K df;3 [10]. In the second, K df;3 is related to the physical scattering amplitudes via integral equations. 2,3 The formalism has been tested in several ways, most notably by reproducing the known finite-volume dependence of a weakly interacting threshold state and of an Efimov-like bound state [26][27][28][29].
A crucial issue yet to be considered, however, is whether the formalism is usable in practice. Indeed, in recent papers introducing an alternative approach based on nonrelativistic effective field theory (NREFT), Refs. [17,18] have suggested that our formalism may be too complicated to use in the analysis of real lattice data. It is the purpose of this work to investigate this issue. We find, in fact, that the status with regard to applicability is more or less identical to that for the NREFT approach: the steps are the same, the number of parameters are the same (when using analogous approximations), and the numerical implementation seems to be of comparable difficulty. There are, however, technical differences that we discuss briefly here and return to in the conclusions. 4 In particular, we note here four advantages of using a relativistic formalism for three-particle physics. First, one aims to constrain the physical observables over the widest energy range possible, and our formalism applies for threeparticle center-of-momentum (c.m.) energies reaching up to 4m (5m if there is a Z 2 symmetry forbidding odd-legged vertices), clearly in the regime of relativistic momenta. Second, in Ref. [11] we describe how to determine 2 → 3 scattering amplitudes from finite-volume energies. Such processes are intrinsically relativistic since the incoming particles must have enough kinetic energy to produce a new particle. Third, it is known in the 1=L threshold expansion that, for weakly interacting systems, three-body effects and relativistic effects enter at the same order in 1=L. Thus it is natural to pursue a formalism that includes both. Fourth, as we describe below, for three noninteracting particles the second and third excited states (as well many higher groups of states) become degenerate in the nonrelativistic limit. Thus the basic counting and locations of noninteracting states, as well as their deformations due to interactions, is very different between the relativistic and nonrelativistic theories. This final point is discussed further in Sec. III A.
In this work, to address the issue of applicability, we primarily use a dynamical approximation similar to that used in the numerical example worked out in Ref. [18], referred to here as the low-energy isotropic approximation. However, in all calculations presented here, we make no kinematical approximations, i.e., we keep the relativistic form throughout. In addition, we restrict attention to theories in which there is a Z 2 symmetry forbidding transitions between even-and odd-particle-number sectors. This is a simplifying approximation that we know, at least formally, how to remove [11]. In short, we conclude that the three-body formalism we have previously derived [9][10][11] is indeed in a form that is suitable for the analysis of some realistic lattice systems.
The remainder of this paper is organized as follows. In Sec. II we present a brief summary of the three-body formalism, and explain the justification for the isotropic approximation, in which the matrix quantity, K df;3 , is replaced by a single function of the total three-particle energy, K iso df;3 . In Sec. III we present several results concerning the three-particle spectrum obtained using the quantization condition, starting with the simplest case of vanishing K iso df;3 and then turning on nonzero values. In cases where this leads to a three-particle bound state, we compare the volume dependence of the bound-state energy to an analytic prediction. We close Sec. III by studying the volume dependence of the threshold state and comparing it to analytic predictions. In Sec. IV we implement the relation between K iso df;3 and the physical scattering amplitude, beginning below threshold and then working directly at threshold. This illustrates how our complete, two-stage formalism can be implemented. In Sec. V we describe how, in certain regimes of parameters, unphysical solutions to the quantization condition can appear, and we discuss their possible origin. We conclude and describe directions for future work in Sec. VI. Two appendices describe some technical details of our numerical implementation of the quantization condition and our methods for solving the integral equations.

II. SUMMARY OF FORMALISM IN THE ISOTROPIC LOW-ENERGY APPROXIMATION
In this section we recall the essential results for the Z 2symmetric case; further details can be found in Refs. [9,10]. The spatial volume is a cube of length L with periodic boundary conditions, so that finite-volume momenta have the form ⃗ k ¼ 2π ⃗ n=L, with ⃗ n a three-vector of integers. The total momentum, ⃗ P, can take any value in this finitevolume set.
Within this setup, the result of Ref. [9] is that, for any fixed values of L and ⃗ P, the finite-volume energy spectrum, fE n ðLÞg, is given by solutions to the quantization condition 5 det½F 3 ðE; ⃗ P; LÞ þ K −1 df;3 ðE Ã Þ ¼ 0: Here the finite-volume-frame energy, E, is related to the c.m.frame energy, E Ã , by the standard dispersion relation, 4 The steps in our approach are also similar to those in the recent relativistic proposal of Ref. [19]. This parametrizes threeparticle interactions using an isobar (dimer) formalism that maintains unitarity. This parametrization is then used both in finite volume to predict the spectrum, and, in a separate calculation, in infinite volume to give the scattering amplitude. We suspect that this formalism will yield similar results to ours. 5 The ultimate aim is for this result to be used to interpret results from lattice QCD simulations. These results inevitably involve errors due to working at nonvanishing lattice spacing. Such effects are not incorporated into the quantization condition, which is a continuum quantum field theory result. Thus, strictly speaking, lattice results should be extrapolated to the continuum limit before they can be used in the quantization condition.
In Eq. (1) the quantities F 3 and K df;3 are matrices in a space labeled by the finite-volume momentum, ⃗ k, of one of the particles (denoted the "spectator") and the angular momentum of the other two in their two-particle c.m. frame. The determinant above acts on this space. K df;3 is an infinite-volume quantity characterizing the underlying local three-particle interaction. It is analogous to the three-body contact terms in the NREFT approach of Ref. [18]. F 3 incorporates both the effects of two-particle scattering and of the finite volume. More specifically, it depends on the two-particle K matrix, K 2 , and on known kinematic finite-volume functions. Its explicit form is given in Eq. (8) below in the approximation we use.
Just as in the two-body sector [6,7,[30][31][32], to use the quantization condition in practice one must truncate the partial waves that contribute, thus reducing the matrices to finite size [9]. This applies here not only to K 2 , as in the twoparticle case, but also to K df;3 . One then proceeds as follows 6 : (1) Perform a two-particle finite-volume analysis to determine K 2 as a function of the two-particle c.m. energy, E Ã 2 , using Lüscher's quantization condition [6,7] and its generalizations.
(2) Use the quantization condition, Eq. (1), and the three-particle spectrum to constrain K df;3 . (3) Determine the relativistic three-to-three scattering amplitude, M 3 , from K df;3 and K 2 by solving the integral equations given in Ref. [10]. Our aim here is to show how this procedure works when we truncate to a single partial wave and make a few further simplifying approximations.
An important technical point is that our formalism includes a smooth cutoff function, Hð ⃗ kÞ, that depends on the spectator momentum ⃗ k. For fixed E, as ⃗ k is increased the c.m. energy in the remaining two-particle subsystem, E Ã 2;k , decreases, dropping first below the two-particle threshold and eventually becoming complex. Our formalism requires that E Ã 2;k is real and positive, E Ã 2;k > 0, and the cutoff function ensures that this condition is satisfied. This means that the sum over ⃗ k is truncated to a finite number of terms. There are two reasons for requiring E Ã 2;k > 0. First, K 2 has a singularity (the left-hand cut) at this point, and this can lead to additional power-law finite-volume effects that are not accounted for in the formalism. Second, the boost to the two-particle c.m. frame becomes unphysical if the condition is not satisfied. There remains, however, considerable latitude in the choice of cutoff function. In particular, the lower limit on E Ã 2;k can lie anywhere in the range from 0 to ð2 − δÞm, with δ a positive constant of order 1. The final results for physical quantities should be independent of this cutoff (up to terms suppressed by e −δmL ). We stress that, if δ is order 1, then the cutoff occurs for spectator momenta satisfying j ⃗ kj ∼ m and thus lying in the relativistic regime. 7 In this work we set δ ¼ 2 throughout.

A. Definition and motivation of the isotropic approximation
The approximation we consider here consists of three parts. First, we restrict K 2 and K df;3 to contain only s-wave interactions between the nonspectator pair. This implies that all matrices appearing in the quantization condition have only the spectator-momentum indices, e.g., K df;3 ¼ K df;3 ðE Ã ; ⃗ k; ⃗pÞ. As noted above, these indices are truncated by the cutoff function. Second, we assume that K df;3 depends only on E Ã and not on the spectator momenta, so that K df;3 ðE Ã ; ⃗ k; ⃗pÞ ≡ K iso df;3 ðE Ã Þ, independent of ⃗ k and ⃗p. Together these give the "isotropic approximation" introduced in Ref. [9]. Finally, we neglect the energy dependence of q Ã 2 cot δðq Ã 2 Þ appearing within K 2 . This corresponds to taking only the leading order (scattering-length-dependent) term in the effective range expansion.
In the remainder of this section we explain why the isotropic approximation is the natural generalization of the s-wave approximation in the two-body case. We begin by recalling the argument for the latter case. We make use of the two independent Mandelstam variables, which we denote by The key input is that, at fixed s 2 and away from isolated poles, K 2 is a finite and thus square-integrable function of cos θ. This means that it admits a convergent decomposition in the Legendre polynomials, P l ðcos θÞ. Alternatively, at fixed s 2 , K 2 is an analytic function of t 2 near threshold so that one can perform a Taylor expansion about t 2 ¼ 0. Combining these two expansions, we deduce that the coefficient of the lth polynomial, call it K 2;l ðq Ã 2 Þ, must scale as q Ã 2 2l as q Ã 2 → 0. This holds because the lth polynomial contains a term proportional to cos l θ and this must correspond to the ðt 2 Þ l term in the Taylor expansion. Thus the s-wave contribution dominates close to threshold.
To justify the isotropic approximation in the three-body case, it is convenient to work with the full divergence-free K matrix, without the decomposition into interacting-pair partial waves. This quantity is a function of generalized Mandelstam variables, which we label 6 This description applies to theories with a Z 2 symmetry. For the general case there are more quantities to determine but the overall approach is the same [11]. 7 In the NREFT approach of Refs. [17,18] there is no corresponding constraint on the sum over spectator momentum, nor is there a need for the cutoff to be smooth. While this simplifies practical calculations, it comes at the price that physical singularities such as the left-hand cut have to be dealt with in some fashion.
where i, j ¼ 1-3, while p i are the initial and p 0 j the final fourmomenta. Note that, at threshold, s ¼ 9m 2 , s ij ¼ 4m 2 ¼ s 0 ij , and t ij ¼ 0. There are many relations between these variables, so that, in addition to s, there are only seven independent kinematic variables. 8 For fixed s, the remaining variables are all "angular," in the sense that they span a compact sevendimensional space [33]. In particular, for fixed s ¼ 9m 2 þ Δ, the quantities that measure the distance from threshold, namely δ ij ≡ s ij − 4m 2 , δ 0 ij ≡ s 0 ij − 4m 2 and t ij , are all bounded in magnitude by cΔ, where c ¼ Oð1Þ. This follows because of the relations together with the fact that δ ij , δ 0 ij and −t ij are all positive. The key input now is that, at fixed s, K df;3 should be an analytic function of the kinematic variables in the vicinity of the threshold. Performing a Taylor expansion about threshold, the leading term is independent of δ ij , δ 0 ij and t ij , with the leading dependence on these variables proportional to Δ. Thus, close to threshold, the dominant contribution is independent of the angular variables. One choice of these variables is given by those introduced in Ref. [9], namely the initial and final spectator momenta introduced above, ⃗ k and ⃗p, together with the initial and final directions of the nonspectator pairs in their respective c.m. frames,â Ã andâ 0Ã . These ten variables are reduced to seven by overall rotation invariance. Thus we conclude that the dominant near-threshold contribution is not only independent ofâ Ã andâ 0Ã (which is the s-wave approximation for K df;3 already introduced above), but also of ⃗ k and ⃗p, yielding the isotropic approximation. 9 We close by commenting that, in the two-particle sector, the s-wave approximation holds both for the K matrix, K 2 , and the scattering amplitude, M 2 . Indeed the harmonic components of these two objects have the same lowmomentum scaling, the usual ðq Ã 2 Þ 2l . This differs from the situation in the three-particle sector, where the argument holds for K df;3 but fails for the scattering amplitude, M 3 . The reason is that the latter exhibits kinematic singularities, discussed at length in Refs. [9,10]. In particular, M 3 is not smooth (indeed it diverges) at threshold and one cannot expect its harmonic coefficients to show low-energy suppression. This is a key advantage of K df;3 over M 3 .

B. Quantization condition in the isotropic approximation
We now return to the main argument. As shown in Ref. [9], the isotropic approximation reduces the quantization condition to an algebraic equation, To reach this form we first note that the determinant over angular momentum appearing in Eq. (1) is trivial given that only the l ¼ 0 contribution to the K matrix is nonzero. Second, in the isotropic approximation, the K matrix is independent of the spectator momentum. Therefore, the only eigenvector of K df;3 in the space of spectator momenta with nonzero eigenvalue is that in which every entry is unity, i.e., j1i ¼ ð1; 1; …; 1Þ. 10 In this way only a onedimensional block of the matrices contributes, leading to Eq. (6). As noted above, this form is analogous to the s-wave approximation of the two-particle formalism. In Fig. 1 we give an example of how this condition is used and compare to the s-wave two-particle case. For any fixed L, ⃗ P and any given finite-volume energy, E n ðL; ⃗ PÞ, Eq. (6) directly gives the value of K iso df;3 ðE Ã Þ at as this is needed to determine F iso 3 , defined below. Given K iso df;3 ðE Ã Þ, one can determine the corresponding M 3 ðE Ã ; Ω 3 0 ; Ω 3 Þ at the same energy. Note that, although we are considering K df;3 only in the isotropic approximation, the three-to-three scattering amplitude still depends on the incoming and outgoing three-particle phase space, indicated here with the shorthand Ω 3 ≡ ð ⃗ k;â Ã Þ. The primary motivation of this work is to demonstrate the practical utility of our result. Thus, for the sake simplicity, we consider only the ð ⃗ P ¼ 0Þ frame. This allow us to use E rather than E Ã to denote the simultaneous finite-volume and c.m.-frame energy. In the same spirit, and following Ref. [18], we take K s 2 to be given by the leading-order term in the threshold expansion, i.e., the term involving the scattering length a.
The expression for F iso where j1i has been defined above, and the sum over the momenta k, p is of finite range because F s 3 is truncated by the cutoff function, Hð ⃗ kÞ. Here and below we keep dependence on E and L implicit. The matrix F s 3 is given by 8 One choice is s 12 , s 13 , s 0 12 , s 0 13 , t 11 , t 22 , and t 33 . 9 It would be interesting to extend this argument to determine the form of the OðΔÞ corrections in terms of ⃗ k, ⃗p,â Ã andâ 0Ã , but this is beyond the scope of the present work. 10 The other eigenvectors, which have vanishing eigenvalues of K df;3 , lead to free three-particle states, as discussed in Ref. [9].
Here ω k and ω ka are the on-shell energies for particles with momenta ⃗ k and ⃗ k þ ⃗ a, respectively, i.e., Other ωs are defined analogously. The two-particle c.m. energy and relative momentum are given by The sum over ⃗ a in Eq. (11) runs over all finite-volume momenta, while the integral is defined as 3 . The principal value (PV) prescription is defined such that the integral is an analytic function of ⃗ k 2 (and is referred to in Ref. [9] as the f PV prescription). Finally, the cutoff function is 11 JðzÞ ¼ This is the form introduced in Ref. [11], chosen to smoothly interpolate between 0 and 1 as E Ã 2;k =m ranges from ffiffiffiffiffiffiffiffiffiffiffi 1 þ α p to the threshold value of 2. In the following we consider α ¼ −1, which gives the maximum allowed range. 12 In Eq. (11) we have labeled both the sum and the integral with a superscript "UV" indicating that an ultraviolet cutoff is required to separately evaluate the sum and integral. In Refs. [9,10] a specific choice of cutoff is used, namely the product of two of the smooth cutoff functions, Hð⃗ aÞHð−⃗ a − ⃗ kÞ. We primarily use this definition in this work as well, but we also make use of the definition given in Ref. [30] for some quantities. These two definitions are described in more detail in Appendix B, where we also explain our method of numerical evaluation. In places where we use both definitions, we refer to that using H-functions for FIG. 1. Examples of solving the quantization conditions in the two-particle (left) and three-particle (right) sectors for ⃗ P ¼ 0 and mL ¼ 6. The two-particle condition in the left panel can be written asF where ω is the energy of the spectating third particle. This is satisfied when the two curves intersect, as indicated by the gray circles. (Here we take the spectator to have ⃗ k ¼ 0 and thus ω ¼ m.) This is closely analogous to the isotropic three-particle quantization condition given by Eq. (6), again satisfied at the indicated intersection points in the right panel. For this example, we take K s 2 from the leading order effective-range expansion with ma ¼ −10, corresponding to an attractive two-particle interaction that pulls the lowest level below E Ã 2 =m ¼ 2. We take 1=K iso df;3 to be a simple polynomial in E=m.

11
Note that, for ⃗ P ¼ 0, Hð ⃗ kÞ ¼ HðkÞ. Nevertheless we keep the more general notation for consistency with Refs. [9][10][11] and because H does depend on ⃗ k when ⃗ P ≠ 0. 12 The relationship between α and the parameter δ used earlier in this section is the UV cutoff asF s HS , and that using the approach of Ref. [30] asF s KSS . The subscripts abbreviate the authors of the article where each cutoff was first introduced.
It is important to note that the freedom to adjust the ultraviolet cutoff here is logically separate from the freedom in the choice of Hð ⃗ kÞ in Eqs. (9) and (12). Varying the UV regulator in Eq. (11) changes the value ofF s only by the exponentially suppressed corrections that we are ignoring throughout. 13 Thus we can choose the regulator that is most convenient for numerical evaluation. By contrast, varying factors of Hð ⃗ kÞ outside a sum-integral difference, such as in Eqs. (9) and (12), leads to changes F s 3 that are, in general, not exponentially suppressed. These are such, however, that K df;3 can in principal be adjusted to keep the low-energy physics unchanged. In other words, an adjustment in the external H functions corresponds to a change in the renormalization scheme. 14 The form of the result for 1=K s 2 in Eq. (9) deserves further explication. Above threshold, where Hð ⃗ kÞ ¼ 1, this form arises from the standard s-wave K matrix, 16πE Ã 2;k tan δ 0 ðq Ã 2;k Þ=q Ã 2;k , keeping only the leading order in the threshold expansion. Below threshold, the result interpolates smoothly to the subthreshold s-wave scattering amplitude, M s 2 , reaching this amplitude when Hð ⃗ kÞ → 0. As explained in Ref. [9], this behavior follows from the choice of pole prescription inF s .
From these definitions we see that F iso 3 depends on E, L and a. For fixed L and a, the spectrum is determined by those values of E for which F iso 3 ðEÞK iso df;3 ðEÞ¼−1. In Appendix A, we describe how we implement this numerically. Here we note two caveats. First, the formalism breaks down as E approaches 5m, where the five-particle channel becomes important. Second, the formalism does not hold if K s 2 has a pole in the region of E Ã 2;k that enters into the calculation, namely ffiffiffiffiffiffiffiffiffi ffi 1þα p <E Ã 2;k =m<ðE Ã −mÞ=m. Note that this restriction includes poles below as well as above threshold.
With the form of K s 2 that we use, Eq. (9), we see that there are no poles above threshold, but there is a pole below threshold if One can show that the right-hand side lies between 0 and m for all allowed values of E, ⃗ k and α. Thus to avoid the poles in general the scattering length must satisfy 15 We stress that negative values of a having arbitrarily large magnitude are allowed, so we can investigate the unitary limit. Indeed, as can be seen from Eqs. (9) and (A1), we can work directly at 1=a ¼ 0, although we do not make use of this possibility in our numerical studies.
C. Relation between K iso df;3 and M 3 We close this section by recalling from Ref. [10] the relation between the infinite-volume quantities K df;3 and M 3 . In the isotropic approximation, this requires solving only one integral equation. This is for the quantity D ðu;uÞ ð ⃗ k; ⃗pÞ that sums up repeated two-particle scattering in which the two particles involved can switch any number of times. It satisfies where, as usual, ⃗ k and ⃗p are spectator momenta, which are now continuous variables. M s 2 ð ⃗ kÞ is the physical s-wave two-particle scattering amplitude with two-particle c.m. energy E Ã 2;k , which in the low-energy approximation is given by and G ∞ is an infinite-volume quantity related toG s , The cutoff functions imply that the integral is of finite range. Note that we are using the iϵ pole prescription here. This is correlated with the appearance of the scattering amplitude M s 2 , rather than K s 2 , in the integral equation. Above threshold, M 3 has singularities at particular, physical kinematic points, and so in Ref. [9] we introduced a divergence-free version of the amplitude The notation here is that S is a symmetrization operator that sums over the three choices of spectator momentum for both initial and final states. The need for such symmetrization implies that M df;3 and M 3 depend not only on the spectator momenta, but also on the directions of the other two particles in their relative c.m. frame, which are given bŷ a Ã andâ 0Ã respectively for the initial and final states. M df;3 13 Strictly speaking, this holds only if the regulator only modifies the terms satisfying jE − ω k − ω a − ω ka j ≫ m, and equals unity when E − ω k − ω a − ω ka ¼ 0.
14 Despite this expectation, we discuss below examples where exponentially suppressed finite-volume artifacts can lead to significant effects, e.g., the unphysical solutions discussed in Sec. V. 15 In fact, for α > −1, a somewhat higher, α-dependent upper limit applies. has the advantage compared to M 3 of being a smooth function of momenta and E, so that, in particular, it is well defined at threshold. It has the disadvantage of depending on the cutoff function H.
In the isotropic approximation, M df;3 is related to K iso df;3 by where These relations involve only integrals over a finite range of momenta. In Appendix A, we discuss how these quantities can be readily determined at or below threshold.

III. NUMERICAL RESULTS
In this section we present a sampling of numerical results, aiming both to provide checks by comparing with several known analytic results, and to give examples of the finite-volume spectrum that emerges for various choices of the scattering parameters. Throughout this section we use units in which m ¼ 1, with the exception of the figures, where for clarity we add back in appropriate factors of m. Most of the details regarding the numerical evaluation of the finite-volume functions are described in the appendices.
A. Energy spectrum with K iso df;3 = M df;3 = 0 We begin by studying the finite-volume spectrum for the special case of K iso df;3 ¼ 0. This in turn implies that M df;3 ¼ 0 and thus that M 3 ¼ SfD ðu;uÞ ð ⃗ k; ⃗pÞg [see Eqs. (25) and (26)]. In words this says that the three-to-three scattering amplitude is given by the sum over all pairwise scattering diagrams in which the two-particle subprocesses are mediated by on-shell two-to-two scattering amplitudes. 16 In this case the quantization condition simply becomes This is a useful starting point because it provides a benchmark for three-particle lattice calculations. If threeparticle energies were found to be consistent with the K iso df;3 ¼ 0 predictions, then it would only be possible to place upper limits on K iso df;3 . By contrast, resolving a shift from these values would give a direct indication of the strength of this local three-body interaction. The solutions to Eq. (30) occur at the poles in F iso 3 . The numerical determination of the positions of these poles is straightforward, as described in Appendix A 1. Examples of the form of F iso 3 are shown in Figs. 1 and 21. In Fig. 2 we plot the low-lying finite-volume spectrum for K iso df;3 ¼ 0 and a ¼ −10, together with the noninteracting three-particle levels. The latter are given by where ⃗m 1 and ⃗m 2 are integer vectors determining the momentum of two of the particles, while ⃗m 12 ¼ − ⃗m 1 − ⃗m 2 determines the momentum of the third. In Table I we collect some information about the low lying noninteracting levels. The values of L that are shown correspond to those used in present lattice QCD simulations (4 ≲ m π L ≲ 6) as well as somewhat larger values that may be accessible in the future. Our interpretation of the a ¼ −10 levels is that they correspond to the first four noninteracting levels, but pushed to significantly lower energies by the strongly attractive two-particle interaction. In particular, the lowest state is not a bound state. We can see this by extending the calculation to larger values of L and observing that it FIG. 2. The four lowest-energy solutions to the quantization condition for K iso df;3 ¼ 0 and ma ¼ −10 (thick, orange) together with the lowest eleven noninteracting levels, i.e., solutions for K iso df;3 ¼ a ¼ 0, (thin, black). We only show results in the energy range for which our formalism is valid, namely E n < 5m. Noninteracting levels are clustered according to momenta that are degenerate in the nonrelativistic theory, as discussed in the text. 16 We stress that both K df;3 and M df;3 are scheme dependent in the sense that physical predictions for a given K df;3 can only be made once a particular form of H has been specified. Thus, when we say that K iso df; approaches the threshold energy E ¼ 3. We refer to this state below as the threshold state.
We have shown several additional noninteracting levels in Fig. 2 in order to illustrate the clustering of excited states. This clustering can be understood by doing a nonrelativistic expansion of the energies. In particular, keeping only the leading term-that present in nonrelativistic quantum mechanics (NRQM)-the noninteracting energies are E n:r: As a result, all states for which the sum of squared momenta are equal become degenerate. This increased degeneracy is indicated by the groupings in Table I. The gaps within the clusters scale as E n ðLÞ − E n:r: whereas the gaps between different clusters scale as 1=L 2 . 17 We note that the splittings within clusters become significant for the values of mL used in present simulations, i.e., those at the lower end of the range displayed. This indicates the importance of including relativistic kinematics in order to gain sufficient precision in the spectrum. One issue that is potentially confusing concerns the degeneracies of the levels shown in Fig. 2. Solutions to the quantization condition in the isotropic approximation are nondegenerate, whereas the noninteracting levels are highly degenerate, as can be seen from Table I. As explained in Ref. [9], the resolution is that, even in the presence of interactions, all but one of the degenerate levels remain at the noninteracting energy when working in the isotropic approximation. These remaining levels will be shifted and split upon inclusion of nonisotropic interactions. We also note that one can project onto the states shown in the figure in practice by using three-particle operators living in A þ 1 irreducible representation (irrep) of the cubic group. This irrep has overlap with the state j1i, and also picks out a single state from each of the noninteracting levels.
In Fig. 3 we show the result of varying the scattering length. The upper left panel shows a ¼ −8, which is very similar to a ¼ −10, while subsequent panels halve the value of a, with the exception of the final panel, which shows the result for a small, positive a. (We recall that the maximum value for which our formalism holds is a ¼ 1.) In these figures we extend the range of mL up to 20, which allows one to see clearly the approach of the levels to the noninteracting curves as a → 0. The larger range also allows us to show additional noninteracting levels, and thus further emphasize the clustering discussed above. Finally, we add to the plot the prediction for the energy of the threshold state in an expansion in powers of a=L, Eq. (36). We observe that this expansion works well for the smallest values of jaj. We investigate this expansion in more detail in Sec. III D below.
B. Energy spectrum for nonzero K iso df;3 We now consider solutions to the quantization condition with nonzero K iso df;3 . We first take energy-independent, negative values of K iso df;3 . As with the two-particle K matrix, small negative values of K iso df;3 correspond to repulsive interactions, and thus push the levels up. We illustrate this in Fig. 4 for the case of a ¼ −10 shown previously for K iso df;3 ¼ 0 in Fig. 2. The levels increase monotonically as K iso df;3 becomes more negative. Large magnitudes of K iso df;3 are required to see a noticeable shift because, as we discuss in more detail below, for small values of K iso df;3 and a, the TABLE I. Summary of noninteracting three-particle energies, in units where the particle mass is m ¼ 1. The first column gives the index of the level, ordered by increasing energy for large L (which means L ≳ 9.5 for these levels). The second column gives the three squared integers describing the individual momenta. The third column gives the degeneracy for identical particles. The final two columns give the energies for L ¼ 4, 6, and 10. Horizontal lines group levels having the same value of the sum ⃗m 2 1 þ ⃗m 2 2 þ ⃗m 2 12 , which are thus degenerate in the nonrelativistic limit. We show all the levels having values of this sum up to 12. Interestingly, for the threshold state, the effect of interactions scales with the power between these two, i.e., as 1=L 3 (as discussed in detail below). effect of the three-body contact interaction on the energy is suppressed by 1=L 6 . In this regard, we stress that such large values of jK iso df;3 j are not unphysical. Indeed, as can be seen from Eq. (26), the three-particle scattering amplitude is finite in the jK iso df;3 j → ∞ limit. This is analogous to the twoparticle sector where K 2 → ∞ corresponds to the unitary limit, Fig. 4 is the appearance of a "bump" in the curves around L ¼ 5.5. If K iso df;3 is made even more negative the spectral lines double back, which is an unphysical result. We discuss this issue further in Sec. V.
What we want to stress here is that, for most values of K iso df;3 , a and L, the quantization condition in the isotropic approximation gives reasonable results, with energy levels that are sensitive to the three-particle interaction.
A more striking example of this sensitivity is shown in Fig. 5, where we use the freedom to allow K iso df;3 to depend on energy to model a three-particle resonance. The ansatz we use is with a "resonance mass" of M R ¼ 3.5. This form is inspired by the standard Breit-Wigner parametrization of the two-particle K matrix, although further investigation is needed to understand if this gives a physical description of three-particle resonances. At the very least, however, it gives a unitary description of three-to-three scattering that, as c → 0, smoothly deforms to a decoupled system of a stable state with mass M R together with three-particle scattering states. For nonzero values of c the two sectors couple and the avoided-level crossings characteristic of a resonance are observed, with the gap increasing with c. For a physical system described by this ansatz, fitting lattice-determined finite-volume levels would give constraints on c, M R and the scattering length a. Consideration of how this ansatz for K iso df;3 converts to M 3 , and whether this gives a useful three-particle resonance description, is a topic for future study.

C. Volume dependence of the energy of a bound state
In this section we provide a quantitative test of our numerical results by studying the volume dependence of the energy of a bound state E B ðLÞ in the unitary regime, jaj ≫ 1. This can be compared with the analytic result of Ref. [34], where κ is defined in terms of the infinite-volume value of the bound state energy, E B ð∞Þ ¼ 3 − κ 2 , jAj 2 is a normalization factor that is expected to be close to unity, 18 while α is of Oð1Þ. This result is valid if κ ≪ 1 (nonrelativistic regime), jaj ≫ 1 (two-particle unitary regime), and κL ≫ 1. In addition, two-particle interactions are assumed to be s-wave dominated, and the bound state is assumed to have J ¼ 0.
We note that this result can also be derived analytically from our quantization condition, under the same assumptions [29]. This derivation applies also within the lowenergy isotropic approximation. However, this derivation requires crucial external input beyond the quantization condition itself, namely the long-distance part of the Schrödinger wave function in the three-body system. Thus agreement with Eq. (35) tests not only our numerical methods, but also that the quantization condition itself correctly reproduces the physics of the bound state. We can also learn where the formula breaks down, i.e., where subleading volume dependence enters.
With this in mind, we have numerically determined the bound-state energy for the parameters a ¼ −10 4 (assuring that we are in the unitary regime) and K iso df;3 ¼ 2500. Note that, in contrast to the previous section, here we choose K iso df;3 positive, as we find that this generically produces a bound state. 19 The results are shown in Fig. 6. We find that for κL > 4 (L > 37) E B ðLÞ is well described by the asymptotic form given in Eq. (35). To be conservative we do our final fit only to data for L > 59 (corresponding to κL > 6.3), as shown in Fig. 6(a). The fit gives κ ¼ 0.106844, corresponding to a binding energy of E B ¼ 2.98858. In addition we find jAj 2 ¼ 0.948, and the  18 For a detailed discussion of the significance of A see Ref. [17]. 19 Why this is the case will become clear in the following section. This result is analogous to the fact that, in the twoparticle case, a bound state occurs when a is large and positive. fact that this result lies close to unity is a strong check on the applicability of the asymptotic form. Figure. 6(b) compares the spectrum to the fit for smaller volumes, 20 < L < 40. None of the data shown in this figure are used in the fit, so the good qualitative agreement for L > 35 provides a strong check that the result of Eq. (35) is consistent with our quantization condition over a wide range of volumes. The deviation as one drops below L ≈ 30 is also expected since κL then becomes too small and the asymptotic form no longer holds. We stress, however, that the solution to the quantization condition continues to be valid for all volumes shown, including the lowest range, 4 < L < 10, shown in Fig. 6(c). For smaller volumes the exponentially suppressed corrections that we are ignoring would start to become sizable.
These results illustrate the potential utility of the quantization condition for analyzing three-particle bound states. Given the value of a from two-particle scattering, one can constrain K iso df;3 near threshold using multiple three-particle scattering states. Extrapolating the results for K iso df;3 to subthreshold energies, one can use the quantization condition to predict the volume dependence of the bound state. We see from Fig. 6(c) that, in the regime of mL accessible to simulations, the finite-volume energy shifts are large, and the asymptotic formula does not hold. Thus the full quantization condition is needed to remove the finitevolume shift and determine the infinite-volume binding energy. We also stress that, in this regime, the bound-state energy is pushed so far below threshold that relativistic momenta are sampled. Thus a relativistic formalism is required to reliably describe even the near threshold state.

D. Volume dependence of the threshold-state energy
In this section we investigate in detail the energy of the threshold state. We have already shown examples of this energy for various values of a in Fig. 3, and our aim here is to provide a detailed comparison with the predicted largevolume behavior. The analytic prediction is where the coefficients are (using the fact that, in our approximation, the effective range, r, vanishes) The numerical values of the constants entering these expressions are 20 The terms through Oð1=L 5 Þ were derived in NRQM in Refs. [35,36]. Relativistic effects first enter at Oð1=L 6 Þ, and the relativistic form ofc 6 was determined in Ref. [26] from our three-particle quantization condition. The derivation was done including all partial waves in K 2 and K df;3 , but holds also in the isotropic limit. Considering only terms through Oð1=L 5 Þ, we see from Eq. (36) that the expansion parameter is a=L. Because of this, for a fixed range of L, we expect the expansion to break down as jaj increases. This is borne out by the results shown in Fig. 3, where only for jaj ≲ 1 does the threshold expansion-shown only through Oð1=L 5 Þ in the plotsprovide a good description over most of the range of L.
Three-particle interactions enter Eq. (36) only at Oð1=L 6 Þ, through the quantity M 3;thr , which is a particular definition of the three-particle divergence-free scattering amplitude at threshold, and is discussed in Sec. IV C and Appendix A 2. As noted earlier, the appearance only at high order implies that the spectrum is only sensitive to threeparticle interactions at smaller values of L, which is the region where simulations are done. But in this small L region, the finite-order threshold expansion might not apply, and one must then use the full quantization condition. By contrast, in this section we are aiming to test our numerical methods by working in a regime where the threshold expansion does hold, namely small jaj and large L. Specifically, we consider a single value of the scattering length, a ¼ 0.41315, and determine the threshold energy to very high accuracy for the range L ¼ 5-60, and with 1=K iso df;3 ¼ 0.04-0. 16. By doing so we are able to extract a value for M 3;thr -which is the only undetermined parameter in Eq. (36). Our results for M 3;thr as a function of K iso df;3 can then be checked against the predictions from the infinite-volume integral equations, as will be discussed in Sec. IV C below.
In Fig. 7 we show that the numerical results from the quantization condition are very well described by the threshold expansion for our choice of scattering parameters. The top curve shows the results from the quantization condition for EðLÞ − 3. Here we suppress the comparison to c 3 =L 3 þ c 4 =L 4 þ c 5 =L 5 as the curves are indistinguishable at this scale (indeed, c 3 =L 3 þ c 4 =L 4 is already indistinguishable from the top curve). The plot also shows the residuals as successively more terms are subtracted from the threshold expansion, as labeled. We see nice convergence for L ≳ 10, with each successive term improving the agreement, and the residuals decreasing in the expected way with L. Note that we subtract the log-dependent piece ofc 6 together with c 5 in the FIG. 7. Log-log plot of EðLÞ=m − 3 vs mL (top curve) determined from the quantization condition for ma ¼ 0.41315 and m 2 K iso df;3 ¼ 10, together with various subtracted curves as labeled. The results are indistinguishable for the two definitions ofF s , with the exception of the lowest (maximally subtracted) curve, where the H-function regulator is shown in orange and that based on Ref. [30] in blue. The oscillations in the former are discussed in the text. 20 Note that we need I to greatest accuracy, followed by J , while K, C 3 etc., are needed to lower accuracy. second-to-last residual, as these terms are of similar numerical magnitude. We stress that we must solve the quantization condition with a numerical accuracy of better than 1 part in 10 8 in order to pick out the maximally subtracted result. This turns out to be straightforward.
The maximally subtracted residual shows oscillatory behavior. To investigate this, we have repeated the calculation replacing the sum-integral difference regulated using H-functions,F s HS , with that regulated following Ref. [30], F s KSS . The residues are indistinguishable for all but the lowest curve, in which we find that the results obtained usingF s KSS do not oscillate. Since the difference between the two choices ofF s is exponentially suppressed, we conclude that the oscillations represent a class of neglected exponentially suppressed finite-volume effects. They are visible here presumably because we are investigating tiny contributions to the energy. Other examples of such effects will be seen below.
As noted above, we can determine M 3;thr from the maximally subtracted results. To do so, we scale up the residual by L 6 and define This quantity is shown in Fig. 8 as a function of 1=L for L ≳ 20. Here we again show the results using the two regulators forF s . The oscillations withF s HS are more pronounced with the new scale, and it is easier to use theF s KSS results to extrapolate to the infinite-volume limit. Averaging quadratic and cubic fits in 1=L to the latter yields M 3;thr =48 ¼ 60.0 AE 0.8, with the uncertainty determined by half the difference between the two fits.
We close this subsection by considering one additional infinite-volume quantity that can be extracted from the threshold energy. With little additional effort we can determine the dependence of the extracted M 3;thr on K iso df;3 , using 21 We determine the derivative numerically by varying K iso df;3 close to 10. 22 The extrapolation to L ¼ ∞ is done either linearly or quadratically in 1=L. An example is shown in Fig. 9. Comparing to the results for R 6 ðLÞ, we see that the derivative removes much of the oscillatory volume dependence, in addition to the first three orders in 1=L. We show the resulting K iso df;3 dependence of the extrapolated derivative in Fig. 10. We take the average of linear and quadratic extrapolations as the central value and half the difference as the uncertainty. The solid line shows the infinite-volume prediction found by solving the integral equation relating M 3;thr to K df;3 , discussed in Sec. IV C below. We stress that this is not a fit to the data, but rather the result of an independent calculation. The agreement between the two results provides a strong check of our numerical implementation of the quantization condition, as well as of the analytic derivation of the threshold expansion in Ref. [26]. The solid curves show quadratic and cubic fits in 1=ðmLÞ to theF s KSS data up to 1=ðmLÞ ¼ 0.05. We take the average of these curves at 1=ðmLÞ ¼ 0 as the central value for the infinite-volume limit, and half the difference as the uncertainty.
FIG. 9. Extrapolation in 1=ðmLÞ of the left-hand side of Eq. (44) evaluated at 1=ðm 2 K iso df;3 Þ ¼ 0.1 and ma ¼ 0.41315. Linear and quadratic fits are done to the region of points indicated by the curves. We stress that this data was generated usingF s HS , but in this case there are only weak oscillations, unlike in Fig. 8. 21 We take the derivative with respect to 1=K iso df;3 because this, rather than K iso df;3 itself, is the more natural quantity entering the quantization condition in the form we use. 22 Given the weak dependence of E on M 3;thr we need to vary E over a very small range. For example, for L ¼ 20, the range E ¼ 3.002067695-3.002067697 leads to a variation in K iso df;3 from ≈6-13 when a ¼ 0.41315.

IV. RELATING K iso df;3 TO THE SCATTERING AMPLITUDE
As explained in Sec. II, to obtain physical infinite-volume quantities given knowledge of K iso df;3 requires solving an integral equation and performing several integrals. In this section we show how this can be done by straightforward extensions of the numerical methods used to solve the quantization condition, as long as we work below or at threshold. We divide this section into three parts. In the first we show results for the quantities needed to relate K iso df;3 to M df;3 below threshold. In the second, we show how, in the case of a three-particle bound state, we can determine a quantity related to the infinite-volume Bethe-Salpeter amplitude. This quantity can then be compared to the predictions of NRQM. Finally, we work directly at threshold and calculate the relation between K iso df;3 and the quantity that enters into the threshold expansion, M 3;thr .
A. Relating K iso df;3 to M df;3 below threshold The relationship between K iso df;3 and M df;3 is given in Eq. (26), which we reproduce here for clarity, making use of the results that Lð ⃗ kÞ ¼ Rð ⃗ kÞ and that L depends only on the magnitude of k in the isotropic approximation, In this subsection we illustrate how to calculate the quantities on the right-hand side of this equation when working below threshold. The methods for doing so are explained in Appendix A 2. The infinite-volume quantities F ∞ 3 and LðkÞ can be obtained simply by taking the L → ∞ limit of appropriate finite-volume quantities. In the case of F ∞ 3 one choice of finite-volume quantity is simply F iso 3 . In Figs. 11 and 12 we show the approach to the L ¼ ∞ limit for F iso 3 and Lð0Þ, respectively, taking E ¼ 2.99 as an example. For fixed a, the approach to the limit is exponential, allowing a controlled extrapolation to L ¼ ∞, although larger values of L are needed as jaj increases. Figure 11 illustrates why, generically, there are bound states for a range of values of K iso df;3 . We recall that, for any finite L, there is a solution to the quantization condition if F iso 3 ¼ −1=K iso df;3 . Since F iso 3 approaches a limiting function of a as L → ∞, which we observe to be monotonically increasing, there will be a bound state with E ¼ 2.99 at some value of a for all values of K iso df;3 in the range −1=F iso Since the limiting function is negative for almost all values of a, most bound states occur with positive values of K iso df;3 . One example (for a different value of E) is the bound state discussed in Sec. III C. Figure 13 shows examples of the k-dependence of LðkÞ for various choices of a. This quantity describes the effect of multiple two-to-two scattering with the scattering pair changing each time to include the spectator of the previous event. As k increases the scattered pair lies increasingly far below threshold. For a bound state, LðkÞ is related to the Bethe-Salpeter amplitude, as discussed in the following subsection. FIG. 10. Dependence of −ðm 4 =48Þ∂M 3;thr =∂ð1=K iso df;3 Þ vs 1=ðm 2 K iso df;3 Þ for ma ¼ 0.41315. The points are obtained from the threshold energy, using extrapolations such as that in Fig. 9, while the solid curve is obtained from solving the integral equation relating M 3;thr and K iso df;3 .
FIG. 11. F iso 3 =m 2 vs ma for E ¼ 2.99m and mL ¼ 40-65, together with an extrapolation to L → ∞ using mL ¼ 50-65. The inset shows the small ma region, in which F iso 3 changes sign.
Here Lð ⃗ kÞ for finite L is given by Eq. (A15). The inset shows the small ma region, within which Lð0Þ changes sign. Note that Lð0Þ ¼ 1=3 when ma ¼ 0.
The results for F ∞ 3 and LðkÞ can be combined to determine results for M df;3 , using Eq. (45). We choose not to quote results here since the symmetrization that is needed is complicated, and the results produced are not transparent. We will, however, quote the corresponding results below when working at threshold.

B. Determining the wave function of the bound state
A specific application of the subthreshold relation between K iso df;3 and M df;3 is provided by the bound state studied in Sec. III C. For the fixed values of K iso df;3 ¼ 2500 and a ¼ −10 4 , one can calculate F ∞ 3 and identify the infinitevolume bound state pole in M df;3 , as described in the previous subsection. Since this is equivalent to solving the quantization condition K iso df;3 ¼ −1=F iso 3 for asymptotically large volumes, one finds the same result for the infinitevolume bound-state energy as from the fit in Sec. III C, namely E B ¼ 2.98858 (corresponding to κ ¼ 0.106844).
The residues of the pole in M df;3 contain information about the Bethe-Salpeter amplitudes of the bound state. Specifically, as discussed in Ref. [29], the unsymmetrized version of M df;3 takes the following factorized form near the bound state: This assumes that pairwise scattering occurs only in the s-wave, as is the case in the isotropic approximation. The quantity Γ ðuÞ ðkÞ is related to the Bethe-Salpeter amplitude by amputating and going on shell, as explained in detail in Appendix B of Ref. [29]. We call Γ ðuÞ ðkÞ the residue function. Combining this expression with Eq. (45) we find that Γ ðuÞ ðkÞ is proportional to LðkÞ, In our approach both F ∞ 3 ðEÞ and LðkÞ are determined by taking infinite-volume limits of appropriate finite-volume quantities. For the purposes of extracting jΓ ðuÞ ðkÞj 2 it turns out to be convenient to define a finite-volume version as jΓ ðuÞ ðkÞj 2 ðLÞ ¼ lim where L L ðE; k; LÞ is defined as the argument of the limit in Eq. (A15). Using this quantity, the infinite-volume limit, is approached more rapidly. Figure 14 shows numerical results for jΓ ðuÞ ðkÞj 2 ðLÞ, calculated by setting E ¼ E B ðLÞ þ δE (with δE ¼ −0.001) and using mL ¼ 60, 65, 70. The results fall on a common curve giving confidence that we have reached the infinite-volume limit. In Ref. [29] we showed that, in NRQM in the unitary limit, the residue function is given by 23  23 It is interesting to note that the leading finite-volume dependence of the bound state energy, given in Eq. (35), is obtained using the leading term in the expansion of the result presented here for Γ ðuÞ ðkÞ about the singularity at k 2 ¼ −κ 2 . This leading term is given in Eq. (100) of Ref. [29]. When evaluated on the real axis, however, it differs substantially from the full result. Thus it is essential to use the full form given here when studying the function for real k jΓ ðuÞ ðkÞ NR j 2 ¼ jcjjAj 2 256π 5=2 3 1=4 with s 0 ¼ 1.00624 and jcj ¼ 96.351, and jAj the quantity entering into Eq. (35). This prediction is also plotted in Fig. 14, and is in excellent agreement with our numerical results. We stress that this curve is a parameter-free prediction and not a fit. However, we do expect there to be relativistic corrections to the relationship between Γ ðuÞ ðkÞ and Γ ðuÞ ðkÞ NR . These should vary in magnitude between of Oðκ 2 =m 2 Þ ¼ Oð1%Þ at k ¼ 0 to of Oðk=mÞ ¼ Oð1Þ for k ≈ m. These expectations are consistent with the small differences we find. What do we learn from this agreement? The derivation of Eq. (50) in Ref. [29] does not use the quantization condition in any way. Instead, it relies only on the definition of the relativistic scattering amplitude and the standard NRQM determination of the bound-state wave function. Thus the agreement is not a consistency check, but rather shows that the relation (45) reproduces the physics leading to the Efimov bound-state solution of the NRQM problem. This is also true for the predicted volume dependence of the bound-state energy, discussed in Sec. III C, but here the test is even more stringent because we are predicting a function and not just a number.
Finally, we note that the curves in Fig. 13 are proportional to the residue functions for bound states that are not in the unitary regime. This is because, for all values of a < 1, one can tune K iso df;3 to give a bound state at E ¼ 2.99, and then use Eq. (47). Since the k dependence comes only from LðkÞ, it follows that jΓ ðuÞ ðkÞj ∝ jLðkÞj. We observe that, away from the unitary regime, the dependence on k varies substantially with a. It would be interesting to compare these results to predictions from NRQM.
C. Relating K iso df;3 at threshold to M 3;df;thr and M 3;thr As discussed above, M df;3 is finite for all energies and choices of momenta (aside from bound-state poles). In particular it is finite at threshold, and we denote its value there by M df;3;thr . This divergence-free scattering amplitude is defined by subtracting an infinite series of terms from the usual three-to-three scattering amplitude, M 3 . An alternate definition of a finite, three-particle threshold amplitude was introduced in Ref. [26], based on subtracting from M 3 only those parts of D ðu;uÞ that contain IR divergences. It is this new quantity, called M 3;thr , that appears in the threshold expansion, Eq. (36) above.
In Sec. III D above, we studied the threshold state predicted by the quantization condition for a ¼ 0.41315 and K df;3 ¼ 10 and found that the volume dependence is very well described by the threshold expansion with M 3;thr =48 ¼ 60.0 AE 0.8. In this section we aim to test this result by directly applying the relation between K iso df;3 and M df;3;thr as well as that between M df;3;thr and M 3;thr [26].
We begin by solving the integral equations relating K iso df;3 to M df;3;thr . At threshold, the general relationship of Eq. (45) simplifies to with the factor of 9 arising from symmetrization. Thus we need only to determine Lð0Þ and F ∞ 3 at E ¼ 3, for the chosen value of a. This is slightly more complicated than the subthreshold determinations discussed above because the finite-volume analogs of L and F ∞ 3 both diverge for E ¼ 3. We have two methods to circumvent this issue. One option is to take the L → ∞ limit for a set of subthreshold values of E (using the method described in previous subsections) and then extrapolate E → 3. An alternative, direct approach is to define modified versions of the finitevolume objects in which the singularity at E ¼ 3 is removed. As explained in Ref. [26], this removal does not affect the L → ∞ limit. The direct approach has the advantage that only one limit need be considered. We have confirmed that the two methods give consistent results and in this subsection only show results using the direct approach. The details of its numerical implementation are summarized in Appendix A 3.
We show the extrapolations for F ∞ 3 and Lð0Þ in Figs. 15 and 16, respectively. Note that here the use of finite volume is simply a tool to discretize the equations, and is not related to the volume of any simulation. We have worked up to L ¼ 100, which, as the figures show, is enough to provide reasonable control over the extrapolation. For F ∞ 3 , the results show oscillations for L ≲ 40, but for larger L settle to a constant value. We estimate the infinite volume value by fitting the large L results to a constant. For Lð0Þ we take the average of the linear and quadratic fits as the central value, and half the difference as the error. We find FIG. 15. F ∞ 3 =m 2 vs 1=ðmLÞ for ma ¼ 0.41315 at threshold. The line indicates a fit of the points shown in blue to a constant, which we use as our estimate of the L → ∞ value.
Inserting these results into Eq. (51) we obtain We now turn to the relation between M 3;thr and M 3;df;thr . The latter is given by whereĨ 1 ,Ĩ 2 and S I are defined in Appendix A 3. In all cases the quantities are obtained by taking infinite-volume limits of appropriate finite-volume quantities. As above, this is only a tool for discretizing integral equations and the parameter L used here does not correspond to the finite volume of the system. Results forĨ 1 , obtained using Eq. (A20) are shown in Fig. 17. Values of L up to 100 are easily attained, and the extrapolation to L ¼ ∞ is well controlled. We find The corresponding extrapolation forĨ 2 , based on Eq. (A21), is shown in Fig. 18. Our result for the infinite-volume limit is Finally, the extrapolation leading to S I is shown in Fig. 19, based on Eq. (A22), yielding Combining these results we find that This is in good agreement with M 3;thr =48 ¼ 60.0 AE 0.8, the indirect value found using the threshold expansion Sec. III D. This provides an important cross-check on our calculations and formalism. Finally, we calculate the dependence of M 3;thr on K iso df;3 , in order to compare to the results obtained using the threshold expansion. Noting that the relation between M 3;thr and M 3;df;thr is independent of K iso df;3 we find, using Eq. (51), that − 1 48 Since we have determined Lð0Þ and F ∞ 3 , we can immediately calculate this quantity. We plot the result in Fig. 10 above as the solid line. The uncertainty in this line from the volume extrapolation, which comes dominantly from the uncertainty in Lð0Þ, is less than the width of the curve. In the figure we compare this result to that obtained above using the threshold expansion and find good agreement.

V. UNPHYSICAL SOLUTIONS
In this section we briefly discuss unphysical solutions that we have identified numerically when solving the quantization condition for certain values of a and K iso df;3 . As we explain below, to guarantee physical solutions for all choices of a and K iso df;3 , F iso 3 ðE; L; aÞ must be a monotonically decreasing function of E. However, we find that this condition is violated for certain parameter choices. An example is shown in Fig. 20. The monotonicity condition on F iso 3 is derived as follows. We consider the finite-volume correlator, where O † ðxÞ is any operator that creates three-particle states. This is the correlator that is used in the derivation of the quantization condition in Ref. [9]. From Eq. (42) of that work, we find that, when restricted to the isotropic approximation, the correlator in the vicinity of a pole has the form where AðEÞ is an infinite-volume matrix element connecting the vacuum to a three-particle state by the action of O.
From the spectral decomposition of C L ðEÞ, we know that the pole has the form with E n ðLÞ the pole position and c a positive, real constant. For this to be true for the pole in Eq. (62), the following condition must be satisfied: FIG. 20. Plot showing an example of the unphysical solutions that arise for certain choices of parameters, here for ma ¼ −10 and large, negative values of K iso df;3 . The left panel shows the lowest finite-volume level, with −10 −4 m 2 K iso df;3 ranging from 1 to 19 in unit steps. This level also appears in Fig. 4, but here we extend the results to more negative values of K iso df;3 . The unphysical behavior is the doubling back of the spectrum, so that there are three, rather than one, levels for a range of mL around 5.4 (the value shown by the vertical dashed line). To understand how this doubling back arises, we show in the right two panels plots of F iso 3 =m 2 vs E=m for mL ¼ 5.4, together with the values of −1=ðm 2 K iso df;3 Þ whose intersections give the solutions. The middle panel looks reasonable, but the enlargement shown in the right panel reveals that F iso 3 is not decreasing monotonically, leading to the triplet of solutions for the largest three values of jK iso df;3 j. As explained in the text, the middle solution corresponds to a pole with an unphysical residue. For a constant K iso df;3 , as used in most of this work and, specifically, in Fig. 20, this implies that F iso 3 must be decreasing at the crossing point. Assuming that there are physical theories with all values of K iso df;3 , the crossing point can occur anywhere along the curve, and thus F iso 3 must be monotonically decreasing. Assuming that all values of a are physical, this monotonicity property must hold in general.
We now return to Fig. 20. This shows an example where F iso 3 does not decrease monotonically with E, but instead, as shown in the right panel, has a small upward excursion. This implies that, in a small range of K iso df;3 there are three solutions to the quantization condition, the middle of which violates the condition Eq. (64). To obtain three states, the spectral curves must double back, as shown in the left panel. Thus this doubling back is an alternative criterion for unphysicality.
Clearly the appearance of such solutions is problematic and needs to be understood. This is work in progress, but based on our tests so far we can offer some remarks. We find that F iso 3 only develops a positive slope in regions where its magnitude is small and, as the volume is increased, these regions always go away and the function becomes "healthy." This leads us to suspect that we are seeing a form of neglected finite-volume effects that are formally exponentially suppressed but with oscillatory energy dependence. These cause problems when small values of F iso 3 are sampled by largemagnitude values of K iso df;3 . In addition, the oscillations in F iso 3 share similarities with the oscillations observed in the threshold state of Fig. 7, and seem to be connected. In that case we found that using a different definition ofF s removed the oscillations and this points to the fact that the smooth cutoff function, Hð ⃗ kÞ, may be the source of the issue. This is plausible because, although smooth, the cutoff function does have vanishing support above a certain value of k. It is well known that sharp cutoffs lead to oscillatory behavior, and the oscillations here might be a related phenomenon.
It is also possible that the unphysical solutions are an indication that the isotropic, low-energy approximation is breaking down for certain choices of parameters. Our approach for now is to avoid the (relatively small) regions in which unphysical solutions occur, while at the same time actively investigating their source.
We close with some more general remarks about the condition Eq. (64). First, we stress that, since K iso df;3 can depend on E, one should in general use the full condition including the derivative of 1=K iso df;3 ðEÞ. Second, we note that a similar condition holds for the two-particle quantization condition, in the s-wave approximation, namely where F s ¼ 2ω kFs ð ⃗ kÞ, E is the total two-particle energy, and E n ðLÞ is here a solution to the two-particle quantization condition, F s ¼ −1=K s 2 . Note that here we are considering also a moving frame, with − ⃗ k being the total momentum. The result (65) can be derived using the result of Ref. [30] for the two-particle correlation function, following similar steps to those outlined above. To our knowledge it has not been presented before. The solutions to the two-particle quantization condition shown in the left panel of Fig. 1 all satisfy this condition.
We can use Eq. (65) to learn about the way in which such consistency conditions can fail. In all examples that we are aware of, F s is a monotonically decreasing function of E. Thus a violation of this condition requires 1=K s 2 to rise sufficiently rapidly at the crossing point. If this is the case, there will be spectral lines that double back, similar to those shown in Fig. 20. If this occurs for some L, the problem will go away as the box size increases, because F s becomes an increasingly steep function of E as L is taken larger. In this case it seems that there are two possible causes for the problem. One is that it is caused by neglected exponentially suppressed corrections, the other that the choice of rapidly increasing 1=K s 2 is simply unphysical within the s-wave approximation. The problem cannot, however, lie with the H functions, since F s can be regulated in other ways.
Our third observation is that the consistency condition, Eq. (65), turns into the condition introduced in Ref. [37] if a subthreshold solution persists as L → ∞. This is because F s → ρ in this limit, and one can then show algebraically that the conditions are equivalent. This is as expected since the persistence of a subthreshold solution is equivalent to the existence of a bound state, and Eq. (65) is just the requirement that this pole, which remains isolated as L → ∞, has a residue with the proper sign.
Finally, we can relate Eq. (65) to a result from Ref. [31]. In that work it was noted that Lellouch-Lüscher factors were physical only if the condition ∂ðδ s þ ϕ P Þ=∂E Ã > 0 was satisfied, where δ s is the s-wave phase shift and ϕ P is the Lüscher pseudophase, related toF s . It is straightforward to show that this condition is equivalent to Eq. (65). We note also that, from the perspective of Refs. [23,38], this equivalence is clear since the Lellouch-Lüscher-like relation is derived there via the same matching condition that leads to Eq. (65) here.

VI. CONCLUSIONS AND OUTLOOK
In this work we have numerically explored the relativistic three-particle quantization condition derived in Refs. [9,10]. In order to capture the key features of the formalism, and to compare the work flow to that described in Ref. [18], we have made a number of simplifications.
Specifically we have restricted attention to vanishing total momentum in the box, truncated the infinite-volume scattering observables to the s-wave, isotropic approximation, and taken the two-particle sector to be dominated by the scattering length in the effective-range expansion.
Within this reduced setup, we find that the quantization condition is numerically straightforward to implement and that a great deal of interesting physics is buried in the simple formula. For example, as summarized in Fig. 3, the condition provides a useful benchmark, by predicting the part of the volume dependence of three-particle energies that is due only to two-particle scattering-i.e., the case where the three-particle contact interaction is neglected, K iso df;3 ¼ 0. This is a useful starting point in lattice calculations since infinite-volume, three-particle scattering information can only be recovered by measuring deviations from these benchmark curves.
Going beyond this, we show how turning on nonzero values of K iso df;3 predicts a rich set of phenomena, including resonancelike avoided level crossings (Fig. 5) and the finite-volume energy shift for an Efimov-like three-particle bound state (Fig. 6). The latter case is particularly clean as we can study the state over a vast range of volumes, mL ¼ 4 to 70, and show that the predicted level matches the asymptotic predictions for κL ≫ 1 (in our case implying mL ≫ 10), but also show how the level deviates from the asymptotic form and thus that the full formalism is needed to make reliable predictions for realistic volumes. Finally, our result also describes the regime of weak interactions and, as we discuss in Sec. III D, we can numerically resolve all known terms in the 1=L expansion of the threshold state, including the logðmLÞ=L 6 dependence.
Beyond predicting detailed finite-volume behavior, our formalism also provides a powerful tool in understanding the infinite-volume scattering of three-particle states. This is because a simple form of K iso df;3 corresponds to a complicated three-particle scattering amplitude, M 3 , with nontrivial phase space dependence generated dynamically by the integral equations relating K iso df;3 to M 3 . The most dramatic example of this is summarized in Fig. 14, where we take two inputs designed to produce a shallow bound state (K iso df;3 ¼ 2500 and a ¼ −10 4 ) and from this predict the wave function with no free parameters. The numerical reproduction of this complicated functional form, which spans many orders of magnitude, gives us confidence that the approach of relating K df;3 to M 3 should be a useful tool in describing three-particle physics for a variety of systems.
Despite these successes, future work is needed to bring this formalism to maturity for use in realistic numerical LQCD calculations. In this direction it is instructive to first compare our approach to that using NREFT, described in Refs. [17,18]. One key difference between our formalism and the NREFT proposal is that the latter uses a hard cutoff in place of our smooth cutoff function Hð ⃗ kÞ, and places this cutoff at much higher spectator momenta.
Recalling that ðE; ⃗ PÞ is fixed, the spectator momentum ðω k ; ⃗ kÞ determines the invariant mass squared of the nonspectator pair to be E Ã2 taking ⃗ k very large takes E Ã2 2;k not only below 4m 2 but in fact to negative values with large magnitude. From the perspective of our approach, dependence on the deep subthreshold values of the two-to-two scattering amplitude is undesirable, leading us to introduce Hð ⃗ kÞ. Among other things, this avoids the region of the left-hand cut in the twoto-two amplitude. This region is accessed in the approach of Refs. [17,18], but the left-hand cut is avoided by restricting the nonrelativistic expansion to a few terms.
We emphasize the role of Hð ⃗ kÞ once more here, because we suspect this to be related to unphysical finite-volume energies that we find for certain values of K iso df;3 . As we describe in Sec. V, the finite-volume function F iso 3 is generally monotonically decreasing with energy, but can exhibit small upward oscillations for volumes up to mL ≈ 6, i.e., including nearly all present-day lattice calculations. These oscillations lead to unphysical solutions when jK iso df;3 j is large enough to intersect them. Understanding the exact nature of these artifacts and modifying the formalism to remove them is clearly crucial. As a first step we note that varying the width of the cutoff function, Hð ⃗ kÞ, can show which regions suffer from these effects. Thus one can identify values of F iso 3 where the artifacts do not arise and restrict attention to values of K iso df;3 that only intersect these regions. This is only a first step as our ultimate goal is a formalism that works for all possible scattering parameters, with no need to identify safe regions numerically.
In addition to addressing the issues mentioned above, future projects include going beyond the isotropic approximation in a systematic way, including the role of the mixing of different angular-momentum states in finite volume. Such mixing is already built into our full quantization condition so the task is one of block diagonalization, or subduction, onto the irreps of the finite-volume symmetry group. Additional formal steps include incorporating K 2 poles, multiple two-and three-particle channels and nonidentical and nondegenerate particles into our formalism, as well as particles with intrinsic spin. Finally, on the side of infinite-volume physics, we are in the process of developing tools to numerically relate K df;3 and M 3 also above threshold. Here it will be crucial to develop realistic parametrizations of three-particle scattering amplitudes, especially in resonant channels. Jefferson Science Associates, LLC, manages and operates Jefferson Lab. We thank Akaki Rusetsky, as well as all other participants of the INT workshop "Multi-Hadron Systems from Lattice QCD", for useful discussions, and Tyler Blanton for comments on the manuscript.

APPENDIX A: NUMERICAL IMPLEMENTATION
The numerical implementation naturally falls into two parts: (i) applying the finite-volume quantization condition, Eq. (6), and (ii) solving the infinite-volume integral equation, Eq. (21), and doing the integrals in Eqs. (27)- (29). We consider these parts in turn. As in the rest of this work, it is convenient to use units in which m ¼ 1. Factors of m can be added back using dimensional analysis.

Implementing the quantization condition
The matricesF s ,G s , etc., entering the quantization condition (6)  To simplify the numerical computation, and, in particular, to allow a straightforward determination of the dependence on a, we rewrite F s 3 as ðA1Þ H FG is a real, symmetric matrix that can be diagonalized as where λ n and jni are its eigenvalues and eigenvectors respectively. To implement the sums over ⃗ k and ⃗p in the definition of F iso 3 , Eq. (7), we use the vector j1i introduced following Eq. (6) in the main text. Putting this together we find Thus, in order to determine F iso 3 , it is convenient to construct H FG , then diagonalize it, and finally calculate the (real) matrix elements h1jF s j1i and h1jF s =ζjni. Given the eigenvalues and these matrix elements we know F iso 3 (at the chosen values of E and L) for all values of a. An example of the a dependence is shown in Fig. 21 where the energy and volume have been fixed to E ¼ 4 and L ¼ 10 respectively. Because of the overall factor of 1=L 3 , F iso 3 has a small magnitude except near the poles at a ¼ 1=λ n . The figure shows a typical example in that most of the poles are in the region a ≥ 1 where our formalism does not hold.
We can substantially reduce the size of the matrices needed in the calculation using group theory [9]. The finitevolume momenta fall into "shells," the members of which are related by elements of the octahedral group O h (which we define to include parity). For example, the first four shells have 1, 6, 12 and 8 members, respectively, with representative elements being 2π ⃗ n=L with ⃗ n ¼ ð0; 0; 0Þ, (0,0,1), (0,1,1) and (1,1,1). The matricesρ,F s andG s are invariant under cubic group transformations, implying that the eigenvectors of H FG lie in irreducible representations (irreps) of the group. The state j1i projects each shell onto the fully symmetric A þ 1 irrep, and the invariance ofF s and ξ implies that only eigenvectors jni lying in the A þ 1 irrep contribute to F iso 3 . The net result is that we need only invert the A þ 1 block of H FG , which is an N sh × N sh matrix, with N sh being the number of shells for which Hð ⃗ kÞ ≠ 0. This drastically reduces the matrix size and concomitantly speeds up the numerical evaluation. For the examples given above, where E ¼ 4, and L ¼ 5, 10 and 20, the values of N sh are 3, 8 and 40 (compared to N ¼ 19, 93 and 895). Note also that the eight poles in Fig. 21 directly correspond to the eight shells for E ¼ 4, L ¼ 10. Although we focus on systems at rest here, one may speed up calculations for systems with nonzero total momenta by consider the irreps of the little groups of O h . F s is our version of the zeta function introduced by Lüscher in Refs. [6,7]. We give a brief description of the methods we use to calculate it in Appendix B, below.

Implementing the integral equations and integrals below threshold
As discussed in Sec. II C, in order to relate K iso df;3 to physical quantities, it is general necessary to solve an integral equation, Eq. (21). In this study we have restricted our attention to energies below and at threshold, where this procedure is relatively straightforward. Here we explain how the infinite-volume limit can be achieved for these kinematics. The below-threshold case is simplest and we discuss it first. In particular, we are interested in determining the quantities appearing in the equation for M df;3 , Eq. (26). In particular, there will be a bound state whenever F ∞ 3 ¼ −1=K iso df;3 , with Lð ⃗ kÞ ¼ Rð ⃗pÞ being proportional to the on-shell Bethe-Salpeter wave function of this bound state.
The determination of F ∞ 3 below threshold is straightforward: we can simply take the L → ∞ limit of our numerical evaluations of F iso 3 , i.e., This is an example of the limit employed in Ref. [10] in order to relate the finite-and infinite-volume scattering amplitudes.
To explain why Eq. (A6) is valid, we first note that, as explained in Ref. [10], This holds for all values of E, but we are interested only in E < 3, whereρ is real. Naively one would have expected the right-hand side to vanish, but the nonzero result arises because we use the PV prescription in the integral inF s . 24 The expression for F s 3 , Eq. (8) Here D ðu;uÞ kp ¼ D ðu;uÞ ð ⃗ k; ⃗pÞ for finite-volume momenta, and we have used the definitions ofG s and G ∞ in Eqs. (12) and (24). 25 Since Eq. (A11) is now a finite matrix equation, its solution is D ðu;uÞ ¼ −L 3 1 1 þð2ωM s 2 ÞG s ð2ωM s 2 ÞG s ð2ωM s 2 Þ: ðA12Þ Using this we can rewrite Eq. (A9) as Comparing to Eqs. (27) and (28), and using the fact that 1=L 3 P k → R ⃗ k as L → ∞, we find the claimed result, Eq. (A6).
By similar arguments, one can evaluate the Bethe-Salpeter amplitudes using either of the forms We stress that, when using the results Eqs. (A6)-(A15), L is no longer playing the role of the spatial volume in the lattice calculation. Instead, it allows for a convenient discretization of integral equations. In particular, we are here interested only in the limiting values as L → ∞, and not to the form of the finite-L corrections.

Implementing the integral equations and integrals at threshold
We now explain how we solve the integral equations and perform the integrals when working directly at threshold. The only change compared to the subthreshold case is that G s has a pole when ⃗ k ¼ ⃗p ¼ 0, which leads to an IR divergence in D ðu;uÞ . However, this IR divergence is absent for all the quantities of interest, either because it is multiplied byρð ⃗ 0Þ, which vanishes at threshold, or because it appears in an IR finite integral. Thus we can regularize in the IR simply by removing the single divergent entry inG s : 24 This result also serves as a useful check of our numerical evaluation ofF s . 25 In general, to take the L → ∞ limit ofG s , we have to introduce a pole prescription [9], but this is not the case below threshold, because there are no poles. Note thatG s does not come with a built-in pole prescription, unlikeF s .
The finite-volume version of D ðu;uÞ is then given by which is simply Eq. (A12) withG s replaced by = G. Now, sinceF s →ρ even at threshold, we can still use Eqs. (A8), (A9) and (A13), as long asG s → = G and the quantity being calculated is IR finite. This leads to the results (all at threshold) Lð The quantities on the right-hand sides of these equations can be evaluated numerically by a slight extension of the work needed to solve the quantization condition, and taking L → ∞ gives the left-hand sides. These are then combined to determine M df;3;thr using Eq. (51) in the main text. Similar methods allow the numerical determination of the relation between M 3;thr and M 3;df;thr that is given in Eq. (123) of Ref. [26]. The basic relation is given in Eq. (54), and we give here the definitions of the quantities I 1 ,Ĩ 2 and S I that appear in that equation.
The final quantity is S I ¼ P ∞ n¼3 I n , where I n is defined in Eq. (124) of Ref. [26]. Given this definition, it is straightforward to evaluate the geometric sum to arrive at S I ¼ 9L 3 1 1 þ 2ωM s 2 = G ð2ωM s 2 = GÞ 4 2ωM s 2 00 þ Oð1=LÞ: In the numerical evaluation of F ∞ 3 , Lð ⃗ kÞ,Ĩ 1 ,Ĩ 2 and S I , we can use the same group-theoretical simplifications as described in the numerical solution of the quantization conditions. Thus the matrices involved have dimensions given by the number of momentum shells. We also note that I 1 is the simplest of the quantities to calculate, since it involves only a column of = G rather than the full matrix.

APPENDIX B: EVALUATION OFF s
In this Appendix we describe how we numerically evaluate the two versions ofF s that we use. Since both differ from the more standard choice of Refs. [7,39], based of zeta-function regulation, we think it is useful to present a short description.
We begin withF s HS . We rewrite Eq. (11) where δF s is exponentially suppressed.F s HS is the form introduced in Eqs. (43) and (44) of Ref. [26], and contains a sum and integral over the vector of integers ⃗ n a , with ⃗ a ¼ ð2π=LÞ⃗ n a and ⃗ b ¼ −⃗ a − ⃗ k, while x 2 ≡ q Ã2 2;k L 2 =ð4π 2 Þ is a quantity that can have either sign. Finally, r is magnitude of a vector whose parts parallel and perpendicular to − ⃗ k are with γ ¼ ðE − ω k Þ=E Ã 2;k . The reason for this rewriting is that it is now easier to numerically implement the PV prescription, following a method introduced in Appendix A of Ref. [26]. Using d 3 n a ¼ γd 3 r, the integral can be rewritten as PV Z ⃗n a Hð ⃗aÞHð ⃗ bÞ x 2 − r 2 ¼ γPV Z ⃗r Hð ⃗aÞHð ⃗ bÞ