Implementing the three-particle quantization condition including higher partial waves

We present an implementation of the relativistic three-particle quantization condition including both $s$- and $d$-wave two-particle channels. For this, we develop a systematic expansion about threshold of the three-particle divergence-free K matrix, $\mathcal{K}_{\mathrm{df,3}}$, which is a generalization of the effective range expansion of the two-particle K matrix, $\mathcal{K}_2$. Relativistic invariance plays an important role in this expansion. We find that $d$-wave two-particle channels enter first at quadratic order. We explain how to implement the resulting multichannel quantization condition, and present several examples of its application. We derive the leading dependence of the threshold three-particle state on the two-particle $d$-wave scattering amplitude, and use this to test our implementation. We show how strong two-particle $d$-wave interactions can lead to significant effects on the finite-volume three-particle spectrum, including the possibility of a generalized three-particle Efimov-like bound state. We also explore the application to the $3\pi^+$ system, which is accessible to lattice QCD simulations, where we study the sensitivity of the spectrum to the components of $\mathcal{K}_{\mathrm{df,3}}$. Finally, we investigate the circumstances under which the quantization condition has unphysical solutions.


Introduction
There has been considerable recent progress developing the formalism necessary to extract the properties of resonances coupling to three-particle channels from simulations of lattice QCD, with three different approaches being followed [1][2][3][4][5][6][7]. For a recent review, see ref. [8]. The outputs of this work are quantization conditions, which relate the finite-volume spectrum with given quantum numbers to the infinite-volume two-and three-particle interactions. This development is timely since simulations now have extensive results for the finite-volume spectrum above the three-particle threshold; see, e.g., refs. [9][10][11] and the recent review in ref. [12]. Turning the formalism into a practical tool remains, however, a significant challenge. To date, this has been done only for the simplest case, in which all particles are spinless and identical, the total momentum vanishes, the two-particle interaction is purely s-wave, and three particles interact only via a momentum-independent contact interaction [4,6,[13][14][15]. 1 This is the analog in the three-particle system of the initial implementations of the two-particle quantization condition of Lüscher [16,17], which assumed only s-wave interactions and vanishing total momentum. In the two-particle case, such an approximation makes sense for levels close to the two-particle threshold, since higher partial waves are suppressed by powers of the relative momentum. In the meson sector it begins to fail for energies around 1 GeV. Indeed, recent applications of the two-particle quantization condition use multiple partial waves (see, e.g., refs. [18,19]). Similar considerations apply for three particles, and we expect that for many resonances of interest one will need to include higher partial waves.
The aim of this paper is to take the first step in this direction by including the first higher partial wave that enters in the case of identical, spinless particles, namely the d wave. 2 In the language of refs. [3,4,6], we include dimers (two-particle channels) with both = 0 and = 2. At the same time, for consistency, we make a corresponding extension of the three-particle interaction beyond its local (pure s-wave) form. We will explain how to implement the formalism in this generalized setting, and show examples for which the higher-order terms have a significant impact on the finite-volume spectrum.
Three-particle quantization conditions have been developed with three different approaches. These use, respectively, generic relativistic effective field theory analyzed diagrammatically to all orders in perturbation theory (the RFT approach) [1,5,7], nonrelativistic effective field theory (NREFT) [3,4], and unitarity constraints on the two-and three-particle S-matrix elements applied to finite-volume amplitudes (the finite-volume unitarity or FVU approach) [6]. To date, only in the RFT approach has the formalism been worked out explicitly with no limitations on the two-particle partial waves, whereas in the other two approaches the quantization condition has been written down only for s-wave dimers. 3 Therefore we adopt the RFT approach in this work. Specifically, we use the formalism of ref. [1], which applies to identical, spinless particles, with a G-parity-JHEP03(2019)106 2 Threshold expansion of the three-particle quantization condition As noted above, we consider a theory of identical, scalar particles, with interactions constrained only by the imposition of a Z 2 global symmetry that prevents odd-legged vertices. In such a theory, the spectrum of odd-particle-number states in a cubic box of length L, with periodic boundary conditions, is determined by solutions to the quantization condition [1] det F 3 (E, L) −1 + K df,3 (E) = 0 . (2.1) This holds up to finite-volume corrections that are exponentially suppressed, i.e., which fall as exp(−mL) up to powers of L, where m is the mass of the particle. In eq. (2.1), F 3 and K df,3 are matrices with index space { k, , m}, where k ∈ (2π/L)Z 3 is the finite-volume momentum assigned to one of the particles (the "spectator"), while and m specify the angular momentum of the other two (the "dimer"). 4 This matrix space will be truncated, as explained in section 3 below, so that the quantization condition (2.1) becomes tractable. The matrix F 3 is a complicated object given in eq. (3.1) below; all we need to know for now is that it depends on the two-particle K matrix, K 2 . Thus the infinite-volume quantities that enter into the quantization condition are K 2 and the three-particle quasilocal interaction K df,3 . 5 The quantization condition (2.1) is valid only when the CM (center of momentum) energy lies in the range m < E * < 5m, within which the only odd-particle-number states that can go on shell involve three particles (rather than one, five, seven, etc.). Here E * = E 2 − P 2 , with (E, P ) the total four-momentum of the state. As in the previous numerical studies [3,6,13,14], we further restrict our considerations to the overall rest frame, with P = 0, implying E * = E henceforth. We also recall that eq. (2.1) assumes that there are no poles in K 2 in the kinematic regime of interest. We discuss the constraints that this places on the two-particle scattering parameters in section 3.
The aim of this section is to develop a systematic expansion of K df,3 about the threeparticle threshold at E = 3m. To that end, we make use of the fact that, unlike the matrix F 3 , K df,3 is an infinite-volume quantity, and so is defined for arbitrary choices of the three incoming and three outgoing on-shell momenta in the scattering process, and not just for finite-volume momenta. It is also important that it can be chosen to be relativistically invariant, if an appropriate choice of the kinematic function G entering F 3 is made [5] [see eq. (A. 3)].
In the remainder of this section, we first recall the threshold expansion of K 2 and its relation to the partial wave decomposition, and then describe the generalization of the threshold expansion to K df, 3 , extending an analysis given in ref. [13]. Finally, we show how the terms in this expansion are decomposed into the matrix form needed for eq. (2.1).

Warm up: expanding K 2 about threshold
To illustrate the method that we employ for K df,3 , we first consider the simpler, and wellunderstood, case of the two-particle K matrix, K 2 . Since K 2 is relativistically invariant, it depends only on the standard Mandelstam variables s 2 , t 2 and u 2 = 4m 2 − s 2 − t 2 . It is convenient to use dimensionless variables that vanish at threshold, where q * 2 is the magnitude of the momentum of each particle in the CM frame, and c θ is the cosine of the scattering angle. For physical scattering, ∆ 2 , − t 2 and − u 2 are all non-negative, and satisfy implying that − t 2 and − u 2 are both bounded by ∆ 2 .
Since K 2 is known to be analytic near threshold, we can expand it in powers of ∆ 2 , t 2 and u 2 . The previous considerations imply that, for generic kinematics (i.e., θ = 0 or π), all three quantities are of the same order. Bose symmetry implies that the expression must be symmetric under t 2 ↔ u 2 . Thus, through quadratic order we have where the c i are constants (which are real since K 2 is real), and we have used the constraint (2.3) to reduce the number of independent terms. We now decompose this result into partial waves, using (2 + 1)K ( ) 2 ( ∆ 2 )P (cos θ) . (2.5) All odd partial waves vanish by Bose symmetry, while eq. (2.4) leads to (2.7) The first equation gives the first three terms in the effective range expansion for K 2 , while from the second equation we recover the well-known result that K 2 ∝ q * 4 2 near threshold. By extending this analysis, one can show that K ( ) 2 only enters when we include terms of O( ∆ 2 ) in the threshold expansion [13].
The threshold expansion has a finite radius of convergence. In particular, we know that K 2 has a left-hand cut at ∆ 2 = −1, so that the radius of convergence cannot be greater than | ∆ 2 | = 1. In practice, we truncate the expansion at the order shown in eqs. (2.6) and (2.7) (and set K ( ) 2 = 0 for ≥ 3), use a cutoff function such that ∆ 2 > −1, and restrict E < 5m implying that ∆ 2 < 3. We are thus assuming that the deviations from the truncated threshold expansion are small over this kinematic range.

Invariants for three-particle scattering
To extend the analysis to the three-particle amplitude K df,3 , we begin by listing the generalized Mandelstam variables, where p i (p i ), i = 1-3, are the incoming (outgoing) momenta. As in the two-particle case, it is convenient to use dimensionless quantities that vanish at threshold, where in the definitions of ∆ i and ∆ i , (i, j, k) form a cylic permutation of (1, 2, 3). These sixteen quantities are constrained by the following eight independent relations, Thus only eight are independent: the overall CM energy (parametrized here by ∆) and seven "angular" degrees of freedom. 6 This counting is as expected: six on-shell momenta with total incoming and outgoing 4-momentum fixed have 3 · 6 − 4 · 2 = 10 degrees of freedom, which is reduced to 7 by overall rotation invariance. For physical scattering, it is straightforward to show that ∆ i , ∆ i , − t ij are all nonnegative, and the constraint equations then lead to the inequality (2.12) Thus all the variables {∆, ∆ i , ∆ i , t ij } can be treated as being of the same order in an expansion about threshold.

Expanding K df,3 about threshold
By construction, K df,3 is a smooth function for some region around threshold. 7 Thus it can be expanded in a Taylor series in the variables {∆, ∆ i , ∆ i , t ij }, which are all treated as being of O(∆). Since K df,3 is real, the coefficients in this expansion must also be real. The expansion must also respect the symmetries of K df,3 , which is invariant under [5]: 8 • Interchange of any two incoming particles: p i ↔ p j ⇒ ∆ i ↔ ∆ j and t ik ↔ t jk • Interchange of any two outgoing particles: 6 We call these variables angular since they span a compact space. 7 More precisely, what is shown in ref. [1] is that K df,3 has no kinematic singularities at threshold, a result that is checked by the explicit perturbative calculations of refs. [21,22]. There can be dynamical singularities due to a three-particle resonance, but, generically, this will lie away from threshold. 8 The first two symmetries hold because we are considering identical bosons. They would not hold in the more general case of nonidentical particles, allowing additional terms to be present in K df,3 .

JHEP03(2019)106
It is then a tedious but straightforward exercise to write down the allowed terms at each order in ∆, and simplify them using the constraints (2.10)- (2.11). Through quadratic order we find df,3 are real, dimensionless constants. We thus see that there is a single term both at leading (zeroth) order and at first order, while there are three independent terms at quadratic order. The particular linear combinations of the quadratic terms that appear in eqs. (2.15) and (2.16) (and in particular the subtraction of ∆ 2 in ∆ (2) A and ∆ (2) B ) are chosen based on our numerical experiments described below in order to ensure that their contributions to the finite-volume spectrum are distinct.
As noted in ref. [13], the leading order contribution to K df,3 in eq. (2.13) is independent of momenta p i and p j . This shows that the isotropic approximation to K df,3 , defined as independence of the seven angular variables, arises naturally in the same way as the s-wave approximation to K 2 . What we add here is the result that K df,3 remains isotropic at O(∆), having only an overall linear dependence on s. Furthermore, at quadratic order, we find only two terms that depend on angular variables (∆ (2) A and ∆ (2) B ), compared to the seven angular variables that are needed to fully characterize three-particle scattering. Thus, if it is a good approximation to truncate the threshold expansion at O(∆ 2 ), the number of parameters needed to describe K df,3 is smaller than one might naively have expected.
For most of our numerical investigations, we have restricted ourselves to quadratic order in the expansion of K df,3 . It is interesting, however, to push the classification to higher order for at least three reasons. First, in order to know how rapidly the number of parameters grows; second, to see which dimer partial waves enter; and, third, to investigate the issue of solutions to the quantization condition with energies given by those of three noninteracting particles (see section 4.4.3). Thus we have classified all terms of cubic order. We find eight independent terms: three that are just ∆ times each of the terms of quadratic order, plus five new angular terms, where σ ∈ S 3 is a permutation of the indices (1,2,3). Thus the number of terms is growing fairly rapidly with order. 9

Decomposing K df,3
In order to use K df,3 in the quantization condition, we need to decompose it into the variables used in its matrix form. This is the analog of the partial wave decomposition of K 2 , described in section 2.1 above. The steps in this decomposition were presented in ref. [1] and we recall them here. The total four-momentum P µ is fixed, in our case to (E, 0). One each of the initial and final particles is designated as the spectator, with three-momenta denoted k and p, respectively. Since K df,3 is symmetric separately under initial and final particle interchange, it does not matter which particles are chosen as the spectators, and we take k = p 3 and p = p 3 . The remaining two particles form the (initial and final) dimers. The total momenta of both dimers are fixed, e.g. to P − p 3 in the initial state. For each dimer, we can boost to its CM frame, and the only remaining degree of freedom is the direction of one of the particles in the dimer in this frame. We take this particle to be p 1 in the initial state, and denote its direction in the dimer CM frame byâ * . Similarly, the direction of p 1 in the final-state-dimer CM frame is calledâ * . Using these variables we can write 10 (2.20) The next step is to set each spectator momentum to one of the allowed finite-volume values, e.g. k = n(2π/L), with n a vector of integers. The final step is then to decompose the dependence onâ * andâ * into spherical harmonics where there is an implicit sum over all angular-momentum indices. This defines the entries in the matrix form of K df,3 . 11 In practice, we use the real version of spherical harmonics, so the complex conjugation in eq. (2.21) has no impact. The simplest example of this decomposition is for the isotropic terms in K df,3 , namely K iso in eq. (2.14). Recalling that E, and thus ∆, is fixed, K iso is simply a constant. This implies that the matrix form of K iso vanishes unless = = 0, and is independent of p, k: The approximation K df,3 = K iso is studied in ref. [13]. 9 We do not think that there is any significance to the fact that the number of terms depending on angular variables through cubic order, i.e. 2 + 5 = 7, equals the number of independent angles in threeparticle scattering. The dependence on these angles can be arbitrarily complicated, so there is not a one-to-one correspondence between variables and functions. 10 As above, the 2·(3+2) = 10 momentum components are reduced to seven independent angular variables by rotation invariance. 11 Note that we follow ref. [1] and drop the vector symbol on the momenta in the matrix indices, in order not to overly clutter the notation.

JHEP03(2019)106
We next work out the decomposition of ∆ (2) A , eq. (2.15), which is conveniently written as ∆ (2) (2.23) The first term depends on k 2 and p 2 , but not onâ * orâ * . This can be seen from with ω k = k 2 + m 2 , and the corresponding result for ∆ 3 . Thus the first term in eq. (2.23) leads to a purely s-wave ( = = 0) contribution to K df,3 , although now with nontrivial dependence on k and p, so this differs from an isotropic contribution. The second term in eq. (2.23) can be rewritten using To obtain the second form one must explicitly boost to the dimer CM frame, in which p − equals 2 a * , with a * 2 = 9m 2 ∆ 3 /4. The first term on the right-hand side of eq. (2.25) is independent ofâ * , and thus again contributes only an s-wave component. The second term, however, depends quadratically onâ * , and thus, through the addition theorem for spherical harmonics, 12 leads to both s-and d-wave contributions. In other words, both K df,3;p00;k00 and K df,3;p00;k2m are nonvanishing. These contributions are straightforward to work out from the above equations, and we do not display them explicitly. The final term in eq. (2.23) differs from the second term only by changing unprimed quantities to their primed correspondents. Thus one finds contributions both to K df,3;p00;k00 and K df,3;p2m ;k00 . Overall, we conclude that the angular dependence in ∆ (2) A leads to both s-and d-wave dimer interactions, although there are no terms with both = 2 and = 2. The latter result arises from the fact that there are no terms in ∆ (2) A that depend on both incoming and outgoing momenta.
Finally, we consider ∆ (2) B , given in eq. (2.16). This is more complicated to decompose because t ij contains both incoming and outgoing momenta, but this same property leads to contributions with = = 2. We provide only a sketch of the decomposition, as the details are tedious, lengthy, and straightforward to automate. Expanding ∆ (2) B , one finds terms that are similar to those dealt with in ∆ (2) A , which lead to additional contributions to K df,3;p00;k00 , K df,3;p00,k2m , and K df,3;p2m ;k00 , and a term proportional to where p ± = p 1 ± p 2 , i, j, r, and s are now spatial vector indices, and S is a tensor that depends on k and p and is symmetric separately under i ↔ j and r ↔ s. By decomposing 12 Again, in practice, we use real spherical harmonics, so the complex conjugation is not needed.

JHEP03(2019)106
S into the spherical tensor basis one finds contributions to the = = 2 part of K df,3 , K df,3;p2m ;k2m , as well as to the other three components. In summary, because the terms of O(∆ 2 ) in K df,3 are at most quadratic in a * and/or a * , they give rise to dimer interactions that are either s-or d-wave. This is the analog of the result derived in section 2.1 that, at the same order, only K (0) 2 and K (2) 2 are present. The generalization to higher order is straightforward. Terms of O(∆ 3 ), can, in principle be cubic in a * and/or a * , but Bose symmetry forbids odd powers. Thus O(∆ 3 ) terms lead only to s-and d-wave contributions to K df,3 , as we have checked explicitly. In order to obtain contributions with = 4 or = 4 one must work at O(∆ 4 ) in the threshold expansion. The pattern continues similarly at higher order.

Implementing the quantization condition
In this section we describe how we numerically implement the quantization condition, eq. (2.1), when working to quadratic order in the threshold expansion. The expression for F 3 is 13 where all quantities are matrices with indices {k, , m}. K 2 is a diagonal matrix 1 2ωK 2 p m ;k m = δ pk δ δ m m 1 where the only nonzero elements are the s-and d-wave terms Here E * 2 2,k = (P − k) 2 is the invariant mass of the dimer, while q * k = E * 2 2,k /4 − m 2 is the momentum of each particle composing the dimer in its CM frame. 14 The expression (3.4) is the standard form for the effective range expansion through quadratic order, with a 0 the s-wave scattering length, r 0 the effective range, and P 0 the shape parameter. Expanding the overall factor of E * 2,k about threshold, and for now ignoring the 1 − H( k) term, one recovers the form given in eq. (2.6). Similarly, aside from the 1 − H term, the expression JHEP03(2019)106 for K (2) 2;k , eq. (3.5), is equivalent to the earlier result, eq. (2.7). Here the leading order term is parametrized in terms of the d-wave scattering length a 2 . 15 The 1 − H terms in the expressions (3.4) and (3.5) arise from the need to introduce a smooth cutoff function H( k) that vanishes for E * 2 2,k ≤ 0. We refer the reader to refs. [1,23] for further explanation of both the need for this cutoff and the manner in which it enters these expressions. It is sufficient to note here that the 1 − H term turns on smoothly only well below the dimer threshold at E * 2,k = 2m. The explicit form of H( k) that we use is given in appendix A.
As noted above, the quantization condition holds only if there are no poles in K 2 in the kinematic regime under study. The kinematic range of q * 2,k is given by −m 2 < q * 2 2,k < 3m 2 (corresponding to 0 < E * 2 2,k < 16m 2 ). The parameters in eqs. (3.4) and (3.5) are thus constrained so that neither right-hand side vanishes for this range of q * 2 2,k . In our numerical investigations, we always use values of the scattering parameters that satisfy these constraints. For a 2 the constraint is that ma 2 < 1, with arbitrarily negative values allowed.
The other two quantities appearing in F 3 are the finite-volume kinematic functions F and G. The former is essentially a two-particle quantity, and thus is diagonal in spectator momenta, though not in the angular-momentum indices: 16 G is a kinematic function that arises from one-particle exchange between dimers, and is thus a quantity that involves all three particles. In particular, it is not diagonal in any of the indices. We give the explicit forms of F and G in appendix A, and provide some details of their numerical evaluation of F in appendix B. An important property is that G p m ;k m is proportional to H( p)H( k), and is thus truncated to the finite number of values of spectator momenta for which H( k) = 0. We call this number N spect (E, L). The same truncation applies to F , due to the factor of H( k) in eq. (3.6). Both matrices are, however, infinite-dimensional in angular-momentum space. This is to be contrasted to K 2 and K df,3 , which are (by approximation) truncated in angular momenta but not in spectator-momentum space. In angular momentum space the dimension is 1 + 5 = 6 when keeping both s and d waves.
Nevertheless, it turns out that these two truncations are sufficient to reduce the quantization condition, eq. (2.1), to a determinant of a 6N spect -dimensional matrix. To show this, we first write the quantization condition as It appears from this rewriting that there will be solutions to the quantization condition when det[F 3 ] → ∞, i.e., when F 3 has a diverging eigenvalue. However, in that case, the second determinant will, for a general K df,3 , also diverge, leading to a finite product. Thus

JHEP03(2019)106
we expect that the only solutions of the quantization condition (2.1) for general K df,3 will be those that also satisfy This also makes sense intuitively, since we expect all finite-volume energies to depend upon the three-particle interaction. The advantage of the form (3.8) is that it has been shown in ref. [1] that it effectively truncates all matrices that appear (i.e., F , G, K 2 and K df,3 ) to N spect entries in spectator-momentum space and to s and d waves in angular-momentum space. By "effectively" we mean that elements of the matrices that lie outside the truncated space do not contribute to the determinant. In the following, we also consider at times the further truncation to only s-wave dimer interactions. This is effected by setting to zero all entries in the matrices having = 2, so that their dimension becomes N spect .
We have now explained how all the matrices contained in the quantization condition eq. (2.1) are constructed, for given values of E and L. We combine these matrices to form F −1 3 + K df,3 , and calculate its eigenvalues. For a given choice of L, the finite-volume spectrum is then given by those values of E for which an eigenvalue vanishes.
The practical calculation of this spectrum is facilitated by decomposing into irreducible representations (irreps) of the symmetry group of finite-volume scattering. For a cubic box with P = 0, this is the cubic group, O h . For the case of pure s-wave dimers, this decomposition has been worked out for the NREFT and FVU quantization conditions in ref. [14]. It has also been used implicitly in the numerical study of the isotropic approximation to the RFT quantization condition in ref. [13], since the isotropic approximation automatically involves a projection onto the trivial (A + 1 ) irrep. 17 The new result that we now present is the generalization of the decomposition to the case in which one has both s-and d-wave dimers.

Projecting onto cubic group irreps
We begin by recalling some useful properties of the cubic group, O h . It has dimension [O h ] = 48, and ten irreps. Its character table can be found, e.g. in ref. [24]. The labels for, and dimensions of, the irreps can be seen in table 1 below. Each finite-volume momentum, k = (2π/L) n k , lies in a "shell" (also known as an orbit) composed of all momenta related to k by cubic group transformations. We refer to this shell as o k . There are seven types of shell, differing by the symmetry properties of the individual elements. We label these by the form of n k : (000), (00a), (aa0), (aaa), (ab0), (aab) and (abc), where a, b and c are all different, nonzero components. They have dimensions N o = 1, 6, 12, 8, 24, 24 and 48, respectively. For example, n k =x lies in the (001) shell of type (00a), and n k =x + 2ẑ lies in the (120) shell of type (ab0). Each element in a shell is invariant under rotations in a subgroup of O h , called its little group, L k . The little groups for all elements in a shell are isomorphic, with dimension The four matrices that enter the quantization condition eq. (2.1), namely 2ωK 2 , K df,3 , F and G, are all invariant under a set of orthogonal transformations U (R), where R ∈ O h . 17 For a more detailed discussion of the isotropic approximation, see appendix F. Specifically, if M is one of these matrices, then Here the Wigner D-matrix is defined in eq. (A.7), while S(R) permutes the spectator momenta within shells: For 2ωK 2 and K df,3 this result follows because they are invariant under rotations, while for F and G it follows from the fact that they are form-invariant under cubic-group rotations if the quantization axis that defines the spherical harmonics is rotated along with the spectator momenta.
One may decompose this reducible representation into irreps I of the cubic group using projection matrices (see, e.g., ref. [25]) where d I is the dimension of I and χ I (R) its character. 18 An important simplifying property of U (R), which carries over to P I , is that it is block-diagonal. For the spectator-momentum indices, this follows because which implies that each U (R) is block diagonal in shells, o. We label the resulting "shell blocks" of P I as P I,o . These shell blocks inherit from D(R) the property of being block diagonal in , and we label the corresponding sub-blocks as P I,o( ) , with = 0 or 2. The result is that we can write P I in the form  This simplified structure allows for more efficient computation of the P I matrices, as explained in appendix C.1.

JHEP03(2019)106
Using these projectors, we can decompose the quantization condition into separate conditions for each irrep. From eq. (3.9) we know that [P I , M ] = 0, for each of the four matrices M , from which it follows that (3.17) Using I P I = 1, and the orthogonality of the projectors onto different irreps, one can then show that the determinant factorizes into that for each irrep where the subscript indicates that the determinant is taken only over the subspace onto which P I projects. Thus the quantization condition for irrep I becomes det sub,I If desired, one can also apply the projectors to all the matrices contained in F 3 , eq. (3.1), so that the entire evaluation of the quantization condition involves matrices of reduced dimensionality.
The number of eigenvalues in a given irrep is given by the dimension of the projected subspaces, d(P I ). This is obtained by summing the dimensions of the sub-blocks, where the sum over o runs over all shells that are "active", i.e., that lie below the cutoff. We explain how the d(P I,o( ) ) are calculated in appendix C.2, and list the results in table 1. From this we learn, for example, that the k = 0 shell contains one A + 1 irrep for = 0, and one each of the E + and T + 2 irreps for = 2. Note that shells can contain multiple versions of a given irrep, e.g., the (00a) shell-type with = 2 contains two versions each of the E + , T + 2 , T − 1 and T − 2 irreps. At this stage it is useful to give an example of how shells become active as E and L are increased. With our cutoff, described in appendix A, the maximum value of | n k |, n k,max , is determined by the vanishing of E * 2 2,k : Although each P I is block diagonal in o and , F −1 3 + K df,3 is generally not. Thus even though each eigenvector of F −1 3 + K df,3 lies in a single irrep, it will generally be a nontrivial linear combination of vectors lying in the subspaces projected onto by P I,o( ) . However, we can still use table 1 to determine how many eigenvalues will be present in a given irrep for (0, 0) (0, 6) (3, 12) (0, 6) (6, 24) (3, 21) (9, 45) Table 1. Dimension of irrep projection sub-blocks for each shell-type and angular momentum, (d(P I,o(0)) , d(P I,o(2) )). Each row corresponds to an irrep of the cubic group O h , whose dimension is also listed for completeness.  000) and (001), are active, so that N spect = 1 + 6 = 7. Then the table tells us that F −1 Table 2. Irreps appearing in the lowest energy levels of three identical noninteracting particles. The first column gives the level number (for values of mL ∼ 5), starting at zero. The states are labeled by the squares of the three vectors n i that determine the momenta of the particles -see eq. (3.23) -and these are given in the second column. The third column gives the degeneracy, and the final column the irreps that appear.
Looking at the other irreps, we see that in this regime there is 1 eigenvalue in in T − 1 , and 6 in T − 2 giving the correct total of 6N spect = 42 eigenvalues. We stress that eigenvalues lying in a given irrep always come in degenerate multiplets corresponding to the dimension of the irrep. Thus, for example, the eight eigenvalues in the E + irrep in the two-shell regime consist of four two-fold-degenerate pairs.
A point that may lead to confusion when we present results in the following section is that the number of eigenvalues of F −1 3 + K df,3 bears no direct relation to the number of solutions to the quantization condition. For there to be a solution, an eigenvalue must vanish, and this occurs only for a subset of the eigenvalues in the energy range of interest. This point can be seen explicitly if the interactions K 2 and K df,3 are weak, for then we expect the number of states to be the same as for noninteracting particles. We quote in table 2 the irreps that appear in the first few three-particle levels for noninteracting particles. These states have energies where n i are integer vectors. As an example of the difference between the dimensions of F −1 3 + K df,3 and the number of solutions, we consider mL = 5 and the A + 1 irrep, and focus on the energy range E/m = 3-5. From figure 1 we see that the number of active momentum shells begins at 2 for E = 3m, increases to 3 at some point, and then reaches 4 below E = 5m. From table 1 we deduce that the corresponding number of eigenvectors in the A + 1 irrep are 3, 6 and 8. By contrast, the free levels in this irrep occur at E = 3m, E = 4.21m, E = 5.08m, . . . . For weak interactions, we expect solutions to the quantization condition only near these three values, and thus we find that, in all cases, the number of eigenvalues of F −1 3 + K df,3 significantly exceeds the number of solutions at, or below, the given energy.
We close this section by noting that the components of K df,3 , given in eq. (2.13), can themselves be decomposed into different irreps. While it is clear that K iso df,3 , eq. (2.14), lies purely in the A + 1 irrep, we also find that the same is true for the K (2,A) df,3 term. The K (2,B) df,3 term, however, has components that lie in the A + 1 , E + , T + 2 and T − 1 irreps. For components lying in the remaining irreps one must go to cubic or higher order in the threshold expansion.

Results
The goal of this section is to illustrate the impact of including d-wave interactions in the quantization condition. In particular, we aim to determine which energy levels and which irreps are particularly sensitive to such interactions. We begin, however, with a case where the impact of d-wave interactions is small, namely the ground state energy with a weak two-particle interaction. This allows us to test of our implementation of the quantization condition in a regime where we can make an analytic prediction. We then consider the impact of a strong d-wave interaction, m|a 2 | ∼ 1, comparing its effect on the ground and excited states, and for different irreps. Next we study the sensitivity of the finite-volume spectrum of the physical 3π + state, with K 2 taken from experiment, to the various terms in K df,3 . And, finally, we discuss the different types of unphysical solutions to the quantization condition that appear.

Threshold expansion including a 2
In this section we consider the energy of the lightest two-and three-particle states in the case of weak two-particle interactions, and with the three-particle interaction K df,3 set to zero. The energy of these states (called E . Here c L , C 3 , and the c i and d i , are numerical constant available in the aforementioned references, and M 3,thr is a subtracted three-particle threshold scattering amplitude, whose definition will be discussed in appendix D. What we observe from these results is that they depend, through O(L −5 ), only on the s-wave scattering length, a 0 , with the effective range r 0 first entering at O(L −6 ).

JHEP03(2019)106
There is no explicit dependence on the d-wave scattering amplitude at this order. We can understand this pattern qualitatively as follows. 20 The typical relative momentum, q, satisfies ∆E ∼ q 2 /m, and thus, since ∆E ∼ a 0 /L 3 , we learn that q 2 ∼ a 0 /L 3 . Using the effective range expansion, eq. (3.4), we then expect that the relative contribution from the r 0 term will be r 0 a 0 q 2 ∼ r 0 a 2 0 /L 3 , and this is indeed what is seen in eqs. (4.2) and (4.3). By the same argument, we expect the q 4 terms proportional to both P 0 and a 5 2 to appear first at relative order O(L −6 ), and thus contribute to ∆E n at O(L −9 ). If this were the case, it would be very challenging to see the dependence of the threshold energies on a 2 .
However, it turns out that there is an additional contribution of O(L −6 ) to ∆E 3 that depends on a 2 , and indeed on all higher partial waves, hidden in M 3,thr . In appendix D we calculate the leading dependence on a 2 in a perturbative expansion in the scattering amplitudes, finding The appearance of a 5 2 , rather than a 2 , follows from our parametrization of the d-wave K matrix, eq. (3.5). In order to isolate the a 2 dependence of ∆E 3 , we consider the difference Substituting eq. (4.4) into the expression for ∆E 3 , eq. (4.3), we obtain the theoretical prediction We have checked that the results from numerically solving the quantization condition are consistent with eq. (4.6). In particular, we have verified that the leading dependence on a 0 , a 2 and 1/L is as predicted. An example of the comparison, showing the L dependence, is given in figure 2. Agreement at the 10% level holds over many orders of magnitude. Based on our tests, we find that the major source of this small discrepancy arises from terms of higher order in a 0 .
This comparison provides a strong cross-check of our numerical implementation. However, for weakly interacting system, such as mesons in QCD, one cannot achieve, using lattice calculations, results for the spectrum with the precision shown in the figure, nor can one work at such large values of mL. We now turn to situations in which a 2 has a numerically more significant effect.

Effects of a 2 on the three-particle spectrum
We begin by studying the strongly interacting regime, where m|a 2 | ∼ 1. This regime, although hardly conceivable in particle physics, represents an interesting academic problem that is relevant in the physics of cold atoms [29,30].
In figure 3, we show the three particle spectrum for E < 4m in two irreps, A + 1 and E + , as a function of negative ma 2 . Here we have fixed the volume to mL = 8.1, and chosen 20 See also appendix C in ref. [28].
2 , eq. (3.5), for which our formalism breaks down. Note that negative a 2 corresponds, at least for small magnitudes, to an attractive interaction, as seen from the result for δE d , eq. (4.6). Since we use a small value of m|a 0 |, the energy levels at the right-hand edges of both plots (where a 2 = 0) lie close to the energies of three noninteracting particles (which are E/m = 3, 3.53, 3.97, 4.02, . . . for mL = 8.1). As m|a 2 | increases, the energies are almost flat, until at a value m|a 2 | ∼ 1, the levels shift rapidly downwards. This shift begins at smaller values of m|a 2 | for excited states. As the magnitude of a 2 increases, the excited states approach lower-lying states until an avoided level crossing occurs. We also observe that states in the E + irrep are more sensitive to d-wave interactions, which seems to be a general feature, as will be seen in the following section.
Another interesting observation from figure 3 is the presence of a deep subthreshold state for m|a 2 | > 1. This resembles the Efimov effect, which describes a three-particle bound state arising from an attractive two-particle interaction m|a 0 | 1 [31]. The Efimov bound state has been reproduced numerically with only s-wave interactions present, both in the NREFT approach [4,14] and in the isotropic approximation of the RFT formalism [13]. Moreover, there is some theoretical work regarding the existence of this generalized Efimov scenario in the presence of d-wave interactions [30], although there is no result concerning the asymptotic volume dependence, unlike in the s-wave case [32]. We have been able to solve the quantization condition numerically up to mL = 37.5 and the bound state energy barely changes, which strongly suggests that it is indeed an infinite volume bound state. 3.8   Results for ma 2 = −1.3 are shown in figure 4. The volume dependence of the energy is dominated by effects of the UV cut-off, which manifest themselves as small oscillations when a new shells become active. These are similar to oscillations observed in several quantities in ref. [13].
We close by commenting on the impact of using a relativistic formalism, as opposed to a NR approach, on the results of this section. We expect that the qualitative features of the results will be unchanged, but that the quantitative energy levels will be changed once they differ significantly from 3m. Thus, for example, we expect that the energy of the subthreshold state will be only slightly changed, since it lies at the border of the NR regime.

Application: spectrum of 3π + on the lattice
The simplest application in QCD for the three-particle quantization condition is the 3π + system, not only from the theoretical point of view -no resonant subchannels -but also from the technical side -no quark-disconnected diagrams and a good signal/noise -19 -JHEP03(2019)106 ratio. Here we use our formalism to predict the 3π + spectrum, using values for the twobody scattering parameters determined from experiment, and a range of choices for the parameters in K df,3 . 21 Our focus will be on how to differentiate effects arising from the different components of K df,3 , listed in eq. (2.13).
An important point in the following is that there is no natural size for the parameters in K df,3 : the magnitudes of the dimensionless coefficients K iso df,3 , df,3 , and K (2,B) df,3 are not constrained. Strictly speaking, we know this only for K iso df,3 , because, in the nonrelativistic limit, it is related to the three-particle contact interaction in NREFT (a relation given explicitly in ref. [8]), and it is well known that the latter interaction varies in a log-periodic manner from −∞ to ∞ as the cutoff varies [33]. But we see no reason why this should not also apply to the other coefficients. In particular, we note that the physical three-particle scattering amplitude, M 3 , does not diverge when K df,3 does [2,13].
We take the parameters describing isospin-2 ππ scattering from ref. [34]: In a lattice simulation, these parameters would be extracted from the two-pion spectrum, using the two-particle quantization condition. Indeed, there is considerable recent work on the 2π + system using lattice methods, in some cases incorporating d-wave interactions [10,[35][36][37][38][39]. We emphasize that one must determine these parameters with high precision in order to disentangle the two-and three-body effects in the three-particle spectrum.
For the relatively weak two-particle interactions of eq. (4.7), the energy levels lie close to the noninteracting energies of eq. (3.23). For the regime of box sizes available in current lattice simulations, 4 m π L 6, there are at most three such levels below the fiveparticle threshold, E = 5m π (above which the quantization condition breaks down). For these levels, the solutions lie in three irreps: Γ = A + 1 , E + , T + 2 (see table 2). We denote the difference between the actual energy and its noninteracting value as where n = 0, 1, . . . labels the levels following the numbering scheme of table 2. It is known that, asymptotically, [40] We stress, however, that the asymptotic result is not numerically accurate for the range of mL that we consider. Let us start from the ground state, which lies in the A + 1 irrep. Here our expectations are guided by the threshold expansion, eq. (4.3). In addition to explicit dependence on a 0 and r 0 , and the implicit dependence on a 2 worked out in section 4.1, the energy depends on K df,3 through the M 3,thr /L 6 term. Following the arguments given in section 4.1, we expect that only K iso df,3 will enter at this order, with dependence on K iso,1 df,3 suppressed by 1/L 3 and that on K iso,2 df,3 , K df,3 by 1/L 6 . This is borne out by our numerical results, shown 21 We ignore QED effects, which are numerically small, and, in any case, cannot be incorporated into the present formalism. s-and d-wave K iso df,3 = 300  = 3m. The two-particle scattering parameters are those in eq. (4.7), aside from the orange curve in the left panel, where only a 0 is nonzero. The three particle scattering parameters are as indicated in the legend, and explained further in the text. We use the convention that a parameter value not given explicitly is set to the value given earlier. For example, the blue line in the left panel has the parameters set to K iso df,3 = 300 and K iso, in figure 5. The left panel compares results with several choices of parameters: (i) those of eq. (4.7) plus K df,3 = 0 (labeled "s-and d-wave" -black, dotted line); (ii) the same as (i) but with K iso df,3 = 300 and all other parameters in K df,3 vanishing (magenta); (iii) the same as (ii) but with K iso,1 df,3 also turned on, taking the three values 135 (blue), 270 (cyan) and 810 (grey); and (iv) the isotropic approximation, i.e., with only s-wave interactions, and a 0 the only nonzero scattering parameter (orange). We see that adding d-wave two-particle interactions has a similar impact to adding K iso df,3 = 300, but that adding K iso,1 df,3 with a similar magnitude has almost no impact.
The right panel shows the dependence on K iso df,3 , with other parameters fixed at the values in eq. (4.7). The range we consider is K iso df,3 = [−1000, +1000]. In order to have sensitivity to K iso df,3 in this range, a determination of ∆E 0 /m with an error of ≈ 0.01 is needed. Such an error can be achieved with present methods. Thus, as noted in ref. [13], if one has a sufficiently accurate knowledge of the two-particle scattering parameters, one can use the ground state energy to determine the leading three-particle parameter K iso df,3 . Indeed, this approach has been carried out successfully in refs. [11,41].
In figure 6, we investigate the sensitivity of the energy of the first excited state to the various two-particle scattering parameters, comparing the two irreps that are present. The magnitude of the energy shifts are comparable to those for the ground state, but the dependence on the scattering parameters differs markedly. This can be understood because the relative momenta between the particles is nonvanishing for the excited state. Denoting generically the relative momenta by q, this satisfies q/m ≈ 2π/(mL) ∼ O(1). Because of this we expect that the higher-order terms in the effective range expansion, i.e. r 0 and P 0 , should play a much more significant role. This is borne out by the results in the figure, particularly for the A + 1 irrep. We observe that the effect of these additional terms is opposite in the two irreps, which is consistent with the prediction of the threshold s-and d-wave Figure 6. Energy shift of the first excited state in the A + 1 irrep (left) and E + irrep (right). In the range of mL shown, E free 1 /m = 4.7-3.9. The quantization condition is solved with only two-particle scattering parameters being nonzero, while K df,3 = 0. When a parameter is nonzero, its value is given by eq. (4.7). The solid orange and red curves include only s-wave dimers, the former having only a 0 turned on ("only a 0 "), with the latter having all three s-wave parameters in K 2 nonzero ("a 0 , r 0 , P 0 "). The dotted black line shows the impact of adding d-wave dimers, with a 2 nonzero ("s-and d-wave").   Figure 7. Energy shift of the first excited state in the A + 1 irrep (left) and E + irrep (right) with various choices of the parameters in K df,3 . The two-particle scattering parameters are given by eq. (4.7) for all curves. The choices of K df,3 parameters is explained by the legend, with the convention that a parameter value not given explicitly is set to the value given earlier. For example, the black line has the parameters set to K iso df,3 = 100, K iso,1 df,3 = 90, and K iso,2 df,3 = 40, while K expansion generalized to excited states [40]. We also see that adding d-wave dimers has almost no impact on the A + 1 irrep (indeed, the effect is smaller than for the ground state) while the impact is comparable to that of r 0 and P 0 for the E + irrep. Qualitatively, this is as expected, since the averaging over orientations in the A + 1 irrep suppresses the overlap with d-wave dimers.
In figure 7 we illustrate the dependence of the same two excited states on the five parameters in K df,3 , eq. (2.13). Because q/m ∼ O(1) we expect that, unlike for the ground state, the energy should be sensitive to all five parameters, and not just to K iso df,3 . s-and d-wave This is borne out for the A + 1 irrep, where there is strong sensitivity to all three isotropic parameters, and a somewhat weaker dependence on K  figure 8. We show results only for those volumes for which the states lie below the five-particle threshold, which requires mL 5.2. The A + 1 energy-shift depends on all parameters in K df,3 , while the E + and T + 2 irreps depend only on K (2,B) df,3 . The results show a similar dependence on parameters as for the first excited states. We also find that the E + and T + 2 irreps show the greatest sensitivity to a 2 of all the states considered.
To sum up, a possible program for determining the coefficients in K df,3 up to quadratic order in the threshold expansion is as follows: 1. Determine a 0 , r 0 , P 0 , and a 2 from the two-body sector using standard two-particle methods.
2. Extract K iso df,3 from the threshold state.
3. Use states in the E + and T + 2 irreps to calculate K

JHEP03(2019)106
Further information could be obtained using moving frames, as has been done very successfully in the two-particle case. The formalism of ref. [1] is still valid, but the detailed implementation along the lines of this paper has yet to be worked out.
We close by commenting on the importance of using a relativistic formalism for the results that we have presented in this section. We note that the excited states whose energies we consider lie in the relativistic regime. For example, at mL = 5.5, the relativistic noninteracting energy of the second excited state is E free 2 = 4.80m, to be compared to the nonrelativistic energy 3m + 2m(2π/(mL)) 2 = 5.61m. Nevertheless, it may be that the energy splittings ∆E Γ n are much less sensitive to relativistic effects, and it would be interesting to implement the NREFT approach including d waves in order to study this. We do expect, however, that the parametrization of the three-particle interaction will require additional terms once the constraints of relativistic invariance are removed.

Unphysical solutions
In this section we describe solutions to the quantization condition that are, for various reasons, unphysical. These fall roughly into two classes (although there is some overlap): solutions that occur at the energies of three noninteracting particles (which we refer to as "free solutions", occurring at "free energies"), and solutions that correspond to poles in the finite-volume correlator that have the wrong sign of the residue. The latter were first observed in ref. [13] within the isotropic approximation. In the following, we begin with a general discussion of the properties of physical solutions, and then discuss the two classes of unphysical solutions in turn.

General properties of physical solutions
We recall here the properties that physical solutions to the quantization condition, eq. (2.1), must obey. This extends the analysis presented in ref. [13] for the isotropic approximation.
The key quantity is the two-point correlation function in Euclidean time, where the operator O † has the correct quantum numbers to create three particles (and here also has P = 0). We stress that its hermitian conjugate is used to destroy the states. Inserting a complete set of finite-volume states with appropriate quantum numbers, we find the standard result where E j > 0 are the energies relative to the vacuum, and the c j are real and positive. Fourier transforming to Euclidean energy and Wick rotating yields where E is the Minkowski energy that appears in the quantization condition. Thus C L (E) is composed of single poles whose residues, for E > 0, are given by i times real, positive coefficients.

JHEP03(2019)106
Next we recall from the analysis of ref. [1] that the correlator can also be written as where A is a column vector, and to obtain the second form we have decomposed F −1 3 +K df,3 in terms of its eigenvalues λ j (E) and eigenvectors v j (E). 22 Since F −1 3 + K df,3 is real and symmetric, the eigenvalues are real.
It follows from comparing eqs. (4.12) and (4.13) that (a) λ j (E) cannot have double zeros. This is because, in the vicinity of a double zero at The same prohibition applies to higher-order zeros.
(b) Eigenvalues of F −1 3 + K df,3 that pass through zero (and thus lead to solutions to the quantization condition) must do so from below as E increases. To understand this, note that, if λ j (E) has a single zero at E = E j , then (4.14) Comparing to eq. (4.12) we learn that This is the generalization of a condition found in ref. [13] for the isotropic approximation (where there is only a single relevant eigenvalue).
Any solutions to the quantization condition that do not satisfy both of these conditions we refer to as unphysical.
We are aware of only three possible sources for unphysical solutions. First, they can arise from the truncation of the quantization condition to a finite-number of partial waves. Second, they could be the result of an unphysical parametrization of K 2 and K df,3 ; for example, the truncation of the threshold expansion for K df,3 could be unphysical. And, finally, the exponentially-suppressed terms that we have dropped could be large in some regions of parameter space, particularly for small mL. We now present examples of unphysical solutions that we have found in our numerical investigation.

Solutions with the wrong residue
In this section we give examples of unphysical solutions to the quantization condition that do not satisfy eq. (4.15), i.e. which lead to single poles whose residues have the wrong sign. These were observed in the isotropic approximation in ref. [13], where it was found that they occurred only when |K iso df,3 | was very large. Here we investigate how this result generalizes in the presence of d-wave dimers. 22 For the sake of brevity, we do not show explicitly that the quantities also depend on L.

JHEP03(2019)106
We first investigate whether unphysical solutions can be induced by adding d-wave interactions alone, with K df,3 = 0. We do not find such solutions for large negative values of ma 2 -the results obtained in section 4.2 all correspond to zero crossings in the correct direction. However, as ma 2 approaches unity (which, as we saw in section 3, is the upper bound allowed for the formalism), we do find examples of unphysical solutions. Since we have seen in sections 4.2 and 4.3 that the impact of d-wave interactions is greater for irreps other than A + 1 , we focus on the E + irrep, and work in the vicinity of the energy of the first noninteracting excited state, E free 1 . In figure 9, we plot the smallest eigenvalue in magnitude in the E + irrep as a function of energy, for two different values of mL and a range of positive values of ma 2 approaching unity. The only other nonvanishing scattering parameter is ma 0 = −0.1. Consider first the left panel, with mL = 8.1. When a 2 = 0, there is a solution at E ≈ E free 1 = 3.53m, as shown by the lowest level in figure 3b. As a 2 is increased, the energy shifts upwards, as expected since positive a 2 corresponds to a repulsive interaction. When ma 2 = 0.9, the level is at E 1 ≈ 3.6m, and moves to yet higher energies as ma 2 increases. These solutions are physical, as shown in the bottomleft inset. For ma 2 = 0.9 and 0.91, however, there is also a single unphysical solution near E = 3.85m, which displays the additional unphysical behavior of having a decreasing energy with increasingly repulsive a 2 . Furthermore, for ma 2 = 0.92, there is a triplet of solutions -two unphysical and one physical. Since they are clearly related, we consider all three to be unphysical. For even larger ma 2 , there are no solutions in the energy range shown.
The right panel, figure 9b, displays a similar pattern, with an additional twist. Here mL = 10, so that E free 1 = 3.36m. The energy of the physical solution lies above this, and increases with increasing ma 2 . There is also an unphysical solution at higher energy, whose energy decreases with increasing ma 2 . The new feature is the presence of a double zero at E free 1 . As discussed above, this is manifestly unphysical since it leads to a double pole in C L (E). It is also unexpected, as its energy lies at that of noninteracting particles. We discuss such solutions in detail in the following section.
Another example of unphysical solutions in shown in figure 10, this time induced by a large, negative value of K (2,B) df,3 . Recall that, out of the parameters in K df,3 , the E + irrep is only sensitive to K  We do not yet understand the source of these unphysical solutions, i.e. which of the three possible sources mentioned at the end of the previous section are most important. This is a topic for future study. Our attitude is that, if a physical solution is well separated from an unphysical one, and its behavior as interactions are made more attractive or repulsive is reasonable, then we accept the physical solution and reject the unphysical one. The examples we have shown occur when the interactions are strong and repulsive, in which limit the two solutions come close together, and at some point become unreliable. For attractive interactions, the two solutions are far apart, often with the unphysical one lying outside the range in which the quantization condition is valid. In this regime, which includes that discussed in section 4.2, we trust the physical solutions.  We conclude by stressing that, in the case of three pions in QCD, the interactions are relatively weak, and we do not expect unphysical solutions to be relevant.

Solutions at free particle energies
This section concerns "free solutions": solutions to the quantization condition that, even in the presence of interactions, lie at one of the energies given in eq. (3.23). We expect that, in general, there will be no such solutions. Exceptions can occur only if the symmetry of the finite-volume three-particle state is such that the chosen interactions do not couple to it. An example in the two-particle sector is that, if P = 0, a finite-volume state lying in the E + irrep would not be shifted from its noninteracting value if only s-and p-wave    interactions were included, since the lowest wave contributing to E + has = 2. One question we address here is where such examples occur in the three-particle sector.
We were prompted to study this issue by finding examples of free solutions in our numerical study. One example has already been seen above, in figure 9b, and further examples are shown in figure 11. The first two plots show solutions with only s-wave channels included. In figure 11a, which shows results for the A + 1 irrep, we see a double zero at the first excited free energy, E free 1 , as well as a solution shifted to slightly higher energies. The latter is expected, since the repulsive interactions should raise the energy of the free state. In the E + irrep, by contrast, there is a single zero at E free 1 , with the unphysical sign -28 -

JHEP03(2019)106
Level Irreps with zeros Zeros removed by ≥ quartic for each Table 3. Irreps in which free zeros appear for the first two excited levels when K df,3 = 0. The "(1)" in the first row denotes that the E free 1 , = 0 free zeros in the E + irrep are single roots with unphysical residue; all other free zeros in the table are (unphysical) double roots. Also noted are the lowest-order terms in the threshold expansion of K df,3 that remove the free zeros. The notation "≥ quartic" indicates that a term of at least quartic order is needed. Note that cubic-order terms are needed to remove the E free 2 , = 0 free zeros in the T − 2 irrep, as neither of the quadratic terms K  figure 11d, have a double-zero at E free 1 . We find similar results for higher excited free energy levels, in which case they appear in an increasing number of irreps. We list these irreps for the first two excited free energies in table 3. There are, however, no free solutions for the lowest free energy E free 0 = 3m. 23 In all the examples we have found, the free solutions are also unphysical -they are either double zeros or single zeros with the wrong residue. We do not know if this is a general result. Also, although the examples shown above are for K df,3 = 0, free solutions also occur when some components of K df,3 are turned on. Indeed, one of the questions we address in the following is which components of K df,3 are required to either remove the free solutions or move them away from E free n . Our first task, however, is to understand in more detail when and why free solutions occur. All such solutions originate from the fact that F and G have single poles at all the free energies. These can lead to poles in F 3 and thus zeros in F −1 3 . We analyze in detail only the lowest two free energies, i.e. those with level number n = 0 and 1 in the notation of table 2, and then draw some general conclusions.
For E ≈ E free 0 = 3m, the only elements of F and G that have poles at E free 0 have vanishing spectator momenta and = 0, 24 specifically .
Here we are using the symbol ∼ to indicate "up to nonpole parts". All other elements of these matrices, and of K 2 , either vanish or are of O(1). From table 1 it now follows that poles in F and G only appear in the A + 1 irrep, and the issue is whether these lead to a pole in F 3 . 23 Strictly speaking, this is only true when one uses the improved form of the quantization condition given in eq. (A.13), and described in appendix A, which removes spurious solutions to eq. (2.1). 24 Pole contributions with = 2 and/or = 2 vanish because, at the pole, a * = a * = 0.

JHEP03(2019)106
To address this we consider the simplest case in which the volume is chosen such that only the lowest two momentum shells are active, which is the case for mL ≈ 5. From We will use a 1 + 2 block notation for the matrices, since this conveys all the necessary information. Close to E free 0 the matrices have the form 25 where O(1) elements are constrained only by the fact that F and G are symmetric. K 2 is a diagonal matrix with O(1) elements. From this it follows that

.19) and thus in turn that
We thus find that free poles at E free 0 cancel in F 3 . This argument generalizes to any number of active shells, since there are no additional poles, and the only change is that the second block in the above analysis is enlarged. The result agrees with our numerical finding that there are no free poles at E free 0 . Next we consider poles at the second free energy, E free 1 . For mL ≈ 4-6 there are then three active shells, so the matrices to consider become larger, e.g. six-dimensional in the A + 1 irrep, and the analysis correspondingly more complicated. We work out the case of the A + 1 irrep in appendix E, both with = 0 channels only and with = 0 and 2 channels included. In both cases we find that is spanned by The factors in eqs. (4.22) and (4.24) result from the form of the spherical harmonics and the size of the first two shells. They are thus kinematical. These analytic results confirm what we find numerically. For example, the double zero at E free 1 shown in figure 11a exactly matches that expected from the analysis of appendix E, and we have checked numerically that it lies in the predicted subspace.
We now discuss how the single zeros at free energies arise. There is a particularly simple case in which we can easily understand these analytically: the E + irrep when we keep only s-wave channels and choose mL such that only the first two shells are active. We must also choose mL such that E free 1 < 5m (so that the formalism applies); one example is mL = 3.8, for which E free 1 = 4.86m. In fact, as shown in table 1, the first shell has no E + component for = 0, so this simple case actually involves only the second shell, for which the E + irrep appears once. Although the E + irrep is two-dimensional, within this space all matrices are proportional to the identity. Thus the matrices are effectively one-dimensional.
The second shell consists of six elements, which we label by the direction of the spectator momentum k in the following order where .

(4.28)
It immediately follows that Thus F 3 indeed has a single pole at E = E free 1 , and F −1 3 a single (doubly degenerate) zero. Increasing L so that there are more active shells does not change the pole structure or the presence of the single zero. We also see that the zero in F −1 3 has a negative coefficient, implying that it decreases through zero, consistent with the behavior seen in figure 11b.
Thus we have understood in a few simple cases why the free zeros listed in table 3 appear. It is interesting to contrast this to the results of ref. [13], where the quantization -31 -

JHEP03(2019)106
condition was studied numerically in the isotropic approximation. In that work no free zeros in F −1 3 were found. At first this may seem puzzling, because the isotropic approximation is a subset of our analysis when we restrict to = 0 channels. The resolution is that the additional isotropic projection that is used is orthogonal to the subspace in which the zeros live. This is demonstrated in appendix F, along with a derivation of the precise relation between the isotropic approximation and the analysis carried out here.
The final stage of our analysis is to study whether the inclusion of components of K df,3 removes the free zeros. Here by "remove" we mean that there is no longer a solution to the quantization condition at a free energy. This can be accomplished either by removing the solution altogether (which is possible for a double zero, which only touches the axis) or by moving it away from the free energy (the likely solution for a single zero). We expect that if K df,3 were not truncated then there would be no free zeros, since there would be some overlap between the state and the three-particle interaction. This is indeed consistent with what we find. What turns out to be surprising, however, is which components of K df,3 that are needed to remove the free zeros.
We first consider the = 0, A + 1 case. To remove the double zero, it must be that the projection of K df,3 into the space of zeros is nonvanishing: where |x 1 is defined in eq. (4.22). Here the square brackets indicate the matrix that results when K df,3 is decomposed into the k m basis and projected into an irrep. Note that this equation need only hold for E = E free 1 , i.e. at the energy of the free zero. The isotropic parts of K df,3 , eq. (2.14), do not solve the problem. These terms have the matrix form Since this vector is orthogonal to |x 1 , it follows that, for all energies, so that eq. (4.30) is not satisfied. The form of |1 K follows from the fact that K iso is independent of the spectator momentum, so that the A + 1 projection simply gives factors of the square root of the multiplicity of the shells. We thus expect that the inclusion of any dependence on the spectator momentum will lead to a [K df,3 ] satisfying eq. (4.30). This is what we find in practice with both of the quadratic terms, i.e. those with coefficients K . This result is an example of a general pattern: the part of K df,3 that "removes" the free zeros comes from terms that involve higher values of than those being included in F −1 3 . Here, we need quadratic terms, which have both = 0 and 2 components, in order to remove the free zeros from the = 0 part of F −1 3 . To be clear, the = 2 components of the quadratic terms play no role; it is simply that by going to higher order one obtains -32 -

JHEP03(2019)106
a more complicated form of the = 0 parts, and this is sufficient to remove the unwanted free zeros. Further examples of this are shown in the last column of table 3, where we list, for all irreps that enter in a given free momentum shell, the terms in K df,3 that remove the free zero.
The second example we consider is the combined = 0 and 2 part of F −1 irrep. In this case, we need [K df,3 (E free 1 )]|x 1 = 0 (4.34) [with |x 1 given in eq. (4.24)] in order to remove the free zeros. We find numerically that this equation is not satisfied by any of the quadratic or cubic terms contributing to K df,3 , but that quartic terms do satisfy it. 26 This exemplifies the general pattern discussed above: quadratic and cubic terms contain only = 0 and 2, while quartic terms include also = 4 parts. We were initially surprised by this result, because K df,3 is an infinitevolume quantity, while |x 1 arises from finite-volume considerations. However, we show analytically in appendix G that orthogonality follows solely from the rotation invariance and particle-interchange symmetry of K df,3 , together with the fact that quadratic and cubic terms contain only = 0 and 2 parts. Thus it is an example of the phenomenon described at the beginning of this section, in which symmetries make the finite-volume state transparent to certain interactions. It is also clear from the arguments in appendix G that all that is required for eq. (4.34) to be satisfied is to use contributions to K df,3 that involve ≥ 4, i.e. terms of quartic or higher order in the threshold expansion. Finally, we consider the case of the single zero in the E + irrep for = 0 channels only, shown in figure 11b. Here we aim to shift the zero away from the free energy. This is accomplished by including a contribution from K df,3 that lives in the E + irrep. As noted in the final paragraph of section 3, the lowest-order term in the threshold expansion for which this is the case is the K (2,B) df,3 term. Thus, once again, we have to use a term in K df,3 that contains higher values of (here = 2) than are included in F 3 .
These theoretical arguments are supported by our numerical results. We show two examples in figure 12. These correspond to the two cases shown in figures 11a and 11b, except that we have turned on K (2,A) df,3 and K (2,B) df,3 , respectively. We expect the double-zero in the former case (A + 1 irrep) to removed by the addition of any quadratic term in K df,3 , and the figure shows that K (2,A) df,3 does the job. In figure 12b, corresponding to the E + irrep, we need to use the K (2,B) df,3 term, since K (2,A) df,3 does not contain an E + component. Since this is a single zero, it is not removed, but is rather shifted to a non-free energy. Note, however, that it remains unphysical because it decreases through zero. In fact, for higher values of K (2,B) df,3 , the zeros coalesce and then disappear.
We close this section with two general comments on the nature of the resolution that we have presented to the problem of unwanted free solutions. The first concerns the result that we need higher-order terms in the threshold expansion of K df,3 in order to remove the free zeros of a given order in F −1 3 . On its face, this invalidates the threshold expansion, for we are evaluating distinct terms in the quantization condition at different orders. We do not think this is the case, however, because we know that, above threshold, all terms in 26 In this case it is crucial to set the energy to E free 1 ; for other energies eq. (4.34) is satisfied.   the expansion of K df,3 are present at some level, and it only takes an infinitesimal value for the coefficient of the requisite higher-order term to remove the unwanted solution. Thus we conclude that we can proceed, in practice, by truncating the expansion of all quantities at the same order in the threshold expansion, and simply ignore the free solutions.
The second comment concerns the fact that our resolution fails if the coefficient of the required parts of K df,3 vanish. In fact, this would require the simultaneous vanishing of an infinite number of terms in the threshold expansion, since higher-order terms in the correct irrep can remove the free solutions. Thus it would require an enormous fine-tuning, which seems highly implausible, especially because there is no enhancement of the symmetry of K df,3 at the tuned point.

Conclusions
The work presented in this paper is the first step towards the systematic inclusion of higher partial waves in the three-particle quantization condition. We have used the generic relativistic field theory (RFT) approach, formulated so that the three-particle scattering quantity, K df,3 , is Lorentz invariant. This invariance proves very important in simplifying the threshold expansion of K df,3 . Indeed, we find that, at quadratic order and for identical particles, only five parameters control the contribution from the three-particle sector, of which only two describe dependence on angular degrees of freedom. This provides a simple starting point for studying the impact of K df,3 . Working at quadratic order implies keeping both s-and d-wave two-particle channels (dimers). We have numerically implemented the quantization condition at this order, and obtained several new results that we now highlight.

JHEP03(2019)106
The first of these is to determine the projection onto irreps of the cubic group including higher partial waves. This has previously been done only for the case of s-wave dimers [14]. The generalization is nontrivial, since both the spectator momentum and the parameters of the dimer transform. While we have worked this out explicitly only for coupled s-and d-wave dimers, the formalism holds for dimers with any angular momentum.
Second, we have understood how the two-particle scattering amplitudes in higher partial waves enter in the 1/L expansion of the energy of the three-particle ground state. We find that all even partial waves enter at O(1/L 6 ), and have calculated analytically the dependence on the d-wave amplitude in the weak-coupling limit and for K df,3 = 0. Although this contribution itself is likely too small to be seen in present simulations of three-particle systems, we have used it as a nontrivial check of our implementation.
Third, we have shown that d-wave interactions, if they are moderately strong, can have a sizable effect on the finite-volume three-particle spectrum. For example, we have presented evidence for a generalized Efimov-like three-particle bound state induced by a strongly attractive d-wave two-particle interaction.
Fourth, we have shown how the five parameters describing K df,3 lead to distinguishable effects on the spectrum of the 3π + system, suggesting that they can be separately determined in a dedicated lattice study. Indeed, this is the system within QCD to which our truncated formalism is most applicable.
Finally, we have characterized solutions to the quantization condition that are unphysical. These presumably arise because of the truncation to a small number of partial waves, and the fact that we have dropped terms that are exponentially suppressed in mL. One class of solutions generally appears when either the two-or the three-particle interactions are strong and repulsive. Our approach is to use parameters such that there are no unphysical solutions near to the physical solutions of interest. The second class of solutions are those that occur at the energies of three noninteracting particles. We have presented numerical evidence and analytical arguments that these are removed if sufficiently highorder terms in K df,3 are included. We expect that other approaches to the three-particle quantization condition will face similar issues, for which our observations may be relevant.
There remain many directions for future study. In order to make our implementation more useful, it is important to generalize it to moving frames. The underlying formalism of ref. [1] applies in all finite-volume frames, but the projectors onto irreps will need to be generalized to account for the reduced symmetry. Another important generalization is to include subchannel resonances, i.e., dynamical poles in K 2 . For this one must implement the formalism of ref. [7], and go beyond the threshold expansion. Finally, we recall that K df,3 is an intermediate quantity, related to the physical three-particle scattering amplitude, M 3 , by integral equations. Since it is only by looking for complex poles in M 3 that one can study three-particle resonances, it is crucial to develop methods to solve the necessary integral equations.
To conclude, we would like to restate that, as it is a relativistic approach, our implementation can simultaneously be useful to both the lattice QCD community and the field of cold atom physics. We are very grateful to MIAPP for its hospitality and stimulating environment.

A Definitions
Here we collect definitions of quantities appearing in F 3 , eq. (3.1), that are not given in the main text.
We begin with the cutoff function: with α H ∈ [−1, 3) a constant. We choose α H = −1, corresponding to the highest cutoff, in all our numerical investigations. For G we use the relativistic form suggested in ref. [5], where b = P − p − k is the momentum of the exchanged particle, p * is the result of boosting p to the CM frame of the dimer for which k is the spectator momentum, and vice versa. Explicitly, we have with k * given by p ↔ k. Finally, Y m ( k) are harmonic polynomials, where Y m are the real spherical harmonics. The elements of G are clearly straightforward to evaluate numerically.

JHEP03(2019)106
For completeness, we quote the real d-wave harmonic polynomials √ The associated Wigner D-matrices are where R is a rotation matrix. They are orthogonal matrices, and implement rotations of the spherical harmonics: Finally, F ( k) is a sum-integral difference that is proportional to the zeta functions that appear in the two-particle quantization condition [16,17]. It requires ultraviolet (UV) regularization, and can be written in various forms that are equivalent up to exponentiallysuppressed corrections. The form that follows from that presented in ref. [1] is where b = P − k − a here, and a * is the result of boosting a to the dimer rest frame, with k the spectator. Here the UV regularization is provided by the product of H functions, and the integral over the pole is defined by the principle value prescription (leading to a real result). Instead, we use a different form that is simpler to evaluate numerically. Following the steps similar to those used in ref. [42], we change variables and introduce a new regularization, finding that, up to exponentially-suppressed corrections, F can be rewritten as where a = n a (2π/L), x = q * 2,k L/(2π), and r( n k , n a ) = n a + n k n a · n k n 2 with k = n k (2π/L). The UV regularization is now provided by the exponential in the integrand, and is parametrized by α > 0. What is shown in ref. [42] is that the α dependence -37 -

JHEP03(2019)106
is exponentially suppressed in L, and that, in practice, one should choose a value that is small enough that the dependence on α lies below the accuracy required. We find that α ≈ 0.5 is usually small enough. An important technical point is that, as seen from eq. (3.6), in the full matrix form F p m ;k m , F ( k) is always multiplied by H( k), from which it follows that γ k is always finite and real whenever F p m ;k m is nonvanishing.
We close this appendix by commenting on the factors of q * (which we use generically for q * 2,k or q * 2,p ) in the denominators of G and F . These lead to poles for particular kinematic configurations, which in turn can lead to solutions to the quantization condition. These solutions appear to be similar to free solutions discussed in section 4.4.3, but are in fact spurious. To understand this we need an argument given in appendix A of ref. [1], which shows that the factors of q * in the denominators are always canceled by corresponding factors in the numerators of K 2 , K df,3 , and the end cap factors A † and A in the finitevolume correlation function C L (E) [see eq. (4. 13)]. The presence of the necessary factors of q * in K 2 can be seen from eq. (3.5), while those in K df,3 arise from the quadratic dependence on a * and a * described in section 3.1. Indeed, one can derive a version of the quantization condition in which all such factors are absent. To do so, we define the matrix Then, from the arguments of ref. [1] we know that we can write the end caps as A = Q A and A † = A † Q, with A and A † nonsingular. Thus an alternative, improved form of the quantization condition is Now we observe that, by simple algebraic manipulations, we can rewrite this form of the quantization condition in terms of Q F Q, Q GQ, Q −1 K 2 Q −1 and Q −1 K df,3 Q −1 , in all of which the factors of q * cancel. Since the difference between the two quantization conditions is a factor of det(Q 2 ), it follows that the solutions to the new form, eq. (A.13), are the same as those to eq. (2.1), except that spurious solutions to the latter, arising from the factors of q * , are removed. In conclusion, we can use the original form of the quantization condition, eq. (2.1), as long as we ignore the spurious solutions.

B Numerical evaluation of F
In this appendix we describe some technical details concerning the evaluation of F ( k).

B.1 Evaluating the integrals
An advantage of the form eq. (A.10) is that the integrals can be evaluated analytically. Dropping overall factors, the integral that is needed is

JHEP03(2019)106
Changing variables to r, we find The remaining integral can be evaluated analytically for all . The explicit result for = 0 was worked out in ref. [13], and we have extended this to the = 2 case. For convenience, we quote both results

B.2 Cutting off the sum
The sum in eq. (A.10) is convergent, but in practice we must introduce a cutoff in order to evaluate it numerically. We use a spherical cutoff, | n a | < n max , and in this section explain how we choose n max . The basic idea is to split the sum S as where S < is the contribution from below the cutoff, and S > the remainder. Assuming that the pole in the summand lies well below the cutoff, then S > can be well-approximated by a remainder integral, R > . We evaluate this integral, and then choose n max such that R > lies below our desired accuracy. The resulting n max depends on E, L and the orbit of k. Dropping overall factors, and changing the overall sign, the sum of interest from eq. (A.10) is Here we have included the cutoff function H( k) that enters in the expression for F p m ;k m , eq. (3.6). Although this is an overall factor, it will play an important role in the determination of n max . The integral R > that results when replacing the sum over n a with an integral is more easily evaluated by changing variables to r. The relation between n a and r, eq. (A.11), can be rewritten as γ k r = n a, − n k 2 , r ⊥ = n a,⊥ , (B.8) with and ⊥ defined relative to k. The cutoff is chosen such that n max n k , implying that the n k /2 term in the expression for r is subleading. Dropping this term, we find -39 -JHEP03(2019)106 that a spherical cutoff on n a corresponds to an ellipsoidal cutoff on r. This makes the integral difficult to evaluate, so we replace this with a spherical cutoff, | r| < Λ, choosing Λ = n max /γ k . We call the resulting integral R Λ . The resulting spherical region is a superset of the original ellipsiodal region, so that we overestimate the remainder, R Λ > R > , since the integrand is positive.
To evaluate R Λ we make two further approximations. First, we drop the x 2 term in the denominator, which is subleading since r 2 x 2 within the region of integration. Second, we make the replacement 4πY m (r)Y m (r) → 1, which leads to an overestimate of the integral. Then we find (B.10) The overall factor of γ k is the Jacobian from changing the integration variable from n a to r. We choose the Λ by specifying a tolerance (we use = 10 −9 ) and numerically solving R Λ = . 27 Given Λ, we then obtain the cutoff for the sum using n max = γ k Λ.
We can now explain why we include the factor of H( k) in S. As | k| approaches the value where H( k) vanishes, γ k diverges. This leads to an increase in n max , both from the factor of γ k in R Λ , and because n max /Λ = γ k . However, this increase is more than compensated by the very rapid drop in H( k) near the end point, so that n max is always finite.

B.3 Using cubic symmetries
Symmetries can be exploited to optimize the computation of F . It follows from eq. (3.6) that F (R k) can be obtained from F ( k) via an orthogonal transformation for any cubic- Here D(R) is the Wigner D-matrix defined in eq. (A.7). Thus once one has computed F ( k) for some finite-volume momentum k, one can use eq. (B.11) to obtain F ( k ) for all k in the same momentum shell. Furthermore, for each initial F ( k) that one computes directly, any symmetries of k can be used to simplify the construction of F ( k). In particular, if R is in the little group of k (so that R k = k), then eq. (B.11) says that F ( k) is invariant under the transformation. This often leads to linear relationships between several matrix elements F m , m ( k), in which case one need only compute the linearly-independent elements in order to construct the full matrix.

JHEP03(2019)106
C Further details of the projection onto cubic group irreps We collect here some results that we have found useful in the computation of the projection matrices and the determination of their properties.

C.1 Computing P I efficiently
The projector P I is defined in eq. (3.14). As explained in the main text, it is block diagonal in momentum shells and in angular momentum, with blocks P I,o( ) . Here we explain how to simplify the computation of P I,o( ) by reducing the sum in eq. (3.14), which runs over all 48 elements of O h , to a sum over the elements of the little group of an element of the shell under consideration. Let k and k be two elements of the orbit. Then, from eqs. (3.14) and (3.15), we have where δ k R k is unity if R k = k and zero otherwise. Thus the sum is restricted to those elements of O h that rotate k into k . A convenient representation of these elements makes use of an (arbitrarily chosen) canonical element of the orbit, denoted k. Let R L k be an element of the little group L k of k. Then all the elements of O h that rotate k to k can be written as R k k R L k R kk , where R kk is any choice of transformation from k to k, and R k k is any choice of transformation from k to k . Thus the number of elements contributing to the sum in eq. (C.1) is [L k ], the dimension of L k . This allows us to rewrite the projector as is the number of elements in the orbit. Once we have constructed the block projectors, we combine them into P I using eq. (3.16). In practice, we want to reduce our original matrices (M = F etc.) down to the part that lives in the projected subspace, which has dimension d(P I ). To do so, we evaluate the eigenvalues and eigenvectors of P I . Since P I is a projector, its eigenvalues λ i are either zero or unity. We keep only the eigenvectors with unit eigenvalues, for these span the projection subspace. We orthonormalize the eigenvectors, and label them . The reduced matrix is then given by . Since, as already noted, M df,3 vanishes when K df,3 = 0, we see that it is the IR finite terms that must contain the contribution to M 3,thr from higher partial waves. What we have learned so far is that, for K df,3 = 0, the a 2 dependence of M 3,thr is given by that of D evaluated at threshold. Here we are interested in determining the leading dependence, which, as discussed in the main text, is proportional to a 5 2 . This is given by .

(D.4)
Here D thr is D( p,â * ; k,â * ) evaluated at E = 3m and p = k = 0, so that there is no dependence onâ * andâ * . In fact, D itself diverges in this limit, but the derivative in eq. (D.4) does not.
To proceed, we need the explicit expression for D, given in ref. [2]. It is obtained by symmetrizing over initial and final momenta the quantity D (u,u) , which is given by (D.5) Here s ≡ d 3 s/(2π) 3 , and theâ * andâ * dependence has been decomposed into partial waves, so that all quantities are implicitly matrices in angular momentum space. The spectator-momentum dependence is, however, kept explicit. M 2 ( p) is the two-particle scattering amplitude for the dimer when the spectator-momentum is p. As for K 2 [see eq. It contains all (even) partial waves, including, in particular, the d-wave amplitude. Finally, G ∞ is given by where the kinematic quantities are the same as those appearing in eq. (A.3). Equation (D.7) is the relativistically-invariant version of the definition given in eq. (81) of ref. [2]. At threshold, only the s-wave part of D (u,u) is nonzero, and symmetrization simply leads to an overall factor of 9: dD thr d(a 5 2 ) given by the leading terms in eqs. (3.4) and (3.5). Inserting these results, we find that I d is IR and UV convergent, so we do not need to actually take the derivative in eq. (D.8). By numerical evaluation we find This gives the leading term in the result (4.4) quoted in the main text. The corrections in (D.10) arise from the subleading terms in the expressions for K ( ) 2 . We close with two further observations. First, a similar calculation with M (2) 2 in I d replaced by any (even) higher-order amplitude leads to a nonzero contribution to M 3,thr . Thus all higher partial waves contribute to ∆E 3 at O(L −6 ). Second, higher-order terms in D (u,u) will also contribute to M 3,thr , although suppressed by powers of a . For example, the first term not shown in eq. (D.5), which has four factors of M 2 , leads to contributions to ∆E 3 proportional to a 3 0 a 5 2 /L 6 and a 2 0 a 10 2 /L 6 . These are of the same order as the corrections in eq. (D.10).

E Free solutions at the first excited energy
In this appendix we analyze free solutions to the quantization condition in the A + 1 irrep at the energy of the first excited noninteracting state, E free 1 = m + 2ω 1 (with ω 1 = m 2 + k 2 L and k L = 2π/L). Our aim is to understand when F −1 3 has zeros at this energy, and to determine their properties. We work with box lengths 4 mL 6 such that there are three active shells, although the final result generalizes straightforwardly to any number of shells.
E.1 A + 1 irrep with s and d waves We first consider the case in which both = 0 and = 2 channels are included. The matrices that enter into the quantization condition are then six dimensional: the first three indices as in eq. (4.17), and the remaining three from the third shell (one with = 0, and two with = 2; see table 1). The free poles enter only in the first two shells, and are proportional to It will be useful to introduce the vectors As in the example discussed in section 4.4.3, all we need to know about the O(1) contributions are that they are real and symmetric. The relative factor of √ 5 between the two terms in v 2 | arises from Y 20 (ẑ) = √ 5Y 00 (ẑ). Combining the results for F and G we find Thus, while the pole parts of F and G are both of rank 2, that of H is of rank 1, due to a partial cancelation.
In the following, we determine the pole structure of F 3 , aiming to find a basis in which this structure is simple. We begin by changing to a more convenient basis, namely |w 1 combined with and any choice of four other vectors filling out the orthonormal set. We use a 1 + 1 + 4 block notation, in which The inverse of H has the form where the quantities α 12 , α 22 , β 22 etc. are given in terms of the O(1) parts of H in a way that is not pertinent. At this stage we can see that F H −1 F will contain a double pole proportional to α 22 that will have the form of an outer product, as well as a complicated single-pole term. Performing the algebra we find where a, b and z are given in terms of the α ij and β 22 . Thus we have learned that F 3 contains a free double pole that can be written -45 -

JHEP03(2019)106
The form of |x 1 is determined entirely by the pole structure of F and H, although the overall coefficient is determined by the O(1) parts. Qualitatively we can say that although F contains two independent poles in this irrep, the H −1 factor cancels one of them, leading to a left-over double pole. We conclude by discussing the impact of the single pole contribution to F 3 . First we note that the coefficient of p can be written as where the new normalized basis vector is Thus in the basis consisting of |x 1 , |x 2 and four other orthonormal vectors, F 3 has the 1 + 1 + 4 block form where f , g and h are known constants. This matrix can be diagonalized using a final, fourth basis. All we need to know here is that, close to the pole, when |p| 1, the shift in the eigenvalues due to the off-diagonal hp term is ±(hp) 2 /(f p 2 + gp) ∼ O(1). Thus in the final basis we have Note that the size of the change to this final basis is proportional to 1/p, and thus vanishes at the zero of F −1 3 , so that the double zero lies in the subspace spanned by |x 1 . In summary, we find that the single pole in F 3 is hidden beneath the double pole, such that in the inverse there is simply a double zero. As L is increased, there are more active shells, but the only change to the result of this section is that the number of vanishing components of |x 1 increases [see eq. (E. 10)]. The nonvanishing components are unchanged.

E.2 A +
1 irrep with only s waves We have repeated the previous analysis for the case of only = 0 contributions and three active shells. 28 The matrices are now three dimensional, with one entry per shell. We do not present the details, except to note that we follow the same steps as in the previous subsection, and find very similar conclusions aside from some changes in factors. In particular, F −1 3 still has a double zero, but this now lives in the space spanned by the vector where entries are ordered as in eq. (4.21). 28 This builds upon, and corrects, the analysis given in appendix C of ref. [1].

F Properties of the isotropic approximation
This appendix recalls the definition of the isotropic approximation, describes its relation to the work of this paper, and explains why the free solutions discussed in section 4.4.3 are absent in this approximation. The isotropic approximation was introduced in ref. [1] and used in the numerical investigation of ref. [13]. It involves three components: (1) Only = 0 dimer channels are included in F , G, K 2 and K df,3 ; (2) The resulting K df,3 is taken to be independent of the spectator momentum, although dependence on the total energy E is allowed; (3) F 3 is projected onto the isotropic vector |1 K , which has a unit entry for every available choice of spectator momentum. Note that the third step automatically picks out solutions in the A + 1 irrep. The isotropic approximation is thus a subset of an approach we use several times in this paper, namely restricting dimers to = 0, keeping only the isotropic part of K df, 3 in the expansion about threshold, and projecting onto the A + 1 irrep. We refer to this as the "low-energy A + 1 approximation". The major difference is the absence of the third stepwe do not project onto |1 K . A minor difference is that, for K df,3 to be purely isotropic, we must work only at linear order in the threshold expansion. Thus we can have at most a linear dependence of K df,3 on E 2 , as opposed to the arbitrary dependence allowed in the isotropic approximation.
To explain the relationship between the two approximations, we begin in the low-energy A + 1 approximation. All matrices, including F 3 , are labeled by an index denoting the shell of the spectator momentum, as shown in eq. (4.21). All matrices have the same finite dimension given by the number of shells lying below our cutoff. Since K df,3 is isotropic, the quantization condition is 29 where the square braces indicate the A + 1 , = 0 matrix, and in this basis. The entries here are the square roots of the sizes of the shells. We can rewrite the determinant in the quantization condition as where we have used det(1+M ) = exp tr ln(1+M ), expanded in M , used the cyclicity of the trace, and resummed. The isotropic approximation consists of keeping only the solutions arising from the numerator on the right-hand side of eq. (F.3), i.e. those satisfying -47 -

JHEP03(2019)106
It follows from eq. (F.3) that any solution in the isotropic approximation is also a solution in the low-energy A + 1 approximation, barring an accidental, and unexpected, juxtaposition with a zero of det([F 3 ]). 30 Thus, aside from this caveat, which appears to be irrelevant in practice, all solutions to the low-energy A + 1 approximation that require a nonzero K iso are also obtained in the isotropic approximation.
What are lost in the isotropic approximation are solutions to the quantization condition (F.1) that arise when an eigenvector of F 3 diverges (so that det([F 3 ]) → ∞) while F iso 3 remains finite. This requires that the corresponding eigenvector of F 3 is orthogonal to |1 K . In our experience, this only happens for solutions that occur at free energies (which, we recall, means one of the energies of three noninteracting particles in the given volume), although we do not know of a fundamental reason why this should be so. Furthermore, it was found numerically in ref. [13] that there are no free solutions in the isotropic approximation. Taken together, these observations suggest that the isotropic approximation picks out all the non-free solutions to the quantization condition obtained in the low-energy A + 1 approximation.
In the remainder of this appendix we explain analytically the result found numerically in ref. [13], namely that there are no free solutions in the isotropic approximation. As discussed in section 4.4.3, such solutions occur first at E = E free 1 , and there yield a double pole in det(F 3 ) lying in the space spanned by |x 1 , eq. (4.22). This pole is, however, absent in the isotropic approximation because 1 K |x 1 = 0, so the pole is removed from F iso 3 . Our aim is to generalize this argument to any excited free energy. We will do so for P = 0, and for an excited state in which the three momenta, labeled k, p and b = − k − a, lie in different shells, e.g. k = k L (0, 0, 1), p = k L (1, 1, 0) and b = k L (−1, −1, −1), with k L = 2π/L. We denote the degeneracies of these shells by N 1 , N 2 , and N 3 , respectively (6, 12 and 8 in our example). For each choice of k from shell 1, we define N 12 as the number of choices of p from shell 2 that can lead to a free solution, and define N 13 analogously. By cubic symmetry N 12 and N 13 do not depend on the choice of k from shell 1. Clearly we have N 12 = N 13 , since each solution contains both a p and b. We define N 23 = N 21 and N 31 = N 32 analogously. The total degeneracy of free-particle solutions is then N sol = N 1 N 12 = N 2 N 23 = N 3 N 31 . (F.5) As above, we denote the = 0, A + 1 parts of F and G by [ F ] and [ G], which are indexed by the shell number. The poles in these matrices occur only when both indices lie in one of the three shells discussed above, and thus we can focus on this three-dimensional subspace. The matrices in this subspace have the form . (F.7) The coefficients in [ F ] count the number of choices of a in eq. (A.9) that lead to the pole. For example, for the (1, 1) element, there are N 12 + N 13 = 2N 12 choices, which combines with the overall factor of 1/2 in F to give the quoted result N 12 . To understand the form of [ G] consider the (1, 2) element of the pole part. This arises from each of the N sol solutions, multiplied by the normalization factors for the A + 1 projections, 1/ √ N 1 N 2 . Then we use Here we are assuming that K 2 does not have a zero at E = E free 1 . It follows from eq. (F.9) that [H] −1 has the form (see, e.g., eq. (C14) of ref. [1]): (F.11) Here |W 2 and |W 3 are any choice for the other two members of an orthonormal basis of which |W 1 is a member. Note that only the coefficient of the first term is known; for all other terms only the power of p is known. We can now calculate the pole part of F iso 3 , which requires projection with 1 K |. As claimed, all poles have canceled from F iso 3 . It is straightforward to generalize this result to the case that two or more shells are the same, and also to moving frames, i.e. P = 0, although we do not present the details here.
-49 -JHEP03(2019)106 G Failure of eq. (4.34) for quadratic and cubic terms in the threshold expansion As noted in the main text, we find numerically that the following results hold, [K (2) df,3 (E free 1 )]|x 1 = [K df,3 (E free 1 )]|x 1 = 0 , (G.1) where the superscript on K df,3 indicates the order in the threshold expansion of K df,3 . The vector |x 1 is given in eq. (4.24), and the square brackets indicate the A + 1 projection of K df,3 expressed in the k m basis. Our aim here is to give an analytic explanation for these results.
We can rewrite eq. (G.1), using the symmetry of K df,3 and the form of |x 1 , as The ordering of the indices is given in eq. (4.23). We recall that the √ 6 here arises because the first shell has 6 elements, while the √ 5 arises because Y 20 (ẑ) = √ 5Y 00 . The superscript on K df,3 indicates that the equation should hold for both the quadratic and cubic terms in the threshold expansion.
We wish to demonstrate eq. (G.2) for any choice of i. To do so we first change notation, recalling from section 2.4 that the k, , m indices can be replaced by dependence on k,â * . Here we are also replacing the spectator-momentum index k with k, both in order to be more explicit, and because K df,3 is an infinite-volume quantity that is defined for all k. At first, we make this change only for the initial-state indices, leading to the hybrid notation K df,3 (E; p, , m ; k,â * ). 31 In terms of this new quantity, we claim that eq. (G.2) holds for any choice of the index i if K (2,3) df,3 (E free 1 ; 0, 0, 0; k,â * ) + c K (2,3) df,3 (E free 1 ; 0, 2, 0; k,â * ) = K (2,3) df,3 (E free 1 ; k Lẑ , 0, 0; k,â * ) + √ 5K (2,3) df,3 (E free 1 ; k Lẑ , 2, 0; k,â * ) , (G.3) is valid for all k andâ * , and for one choice of c. To understand this, first note that (G.3) applies for an arbitrary initial state, and this subsumes all possible values of the finitevolume index i. As for the final state, to obtain eq. (G.2) we need to project onto the A + 1 irrep. Doing so, the second term on the left-hand side of eq. (G.3) vanishes, as can be seen from the absence of an = 2 entry in the A + 1 row of the (000) shell column in table 1. This is why it is sufficient if eq. (G.3) holds for one value of c. The A + 1 projections of the remaining three terms in eq. (G.3) leads to the three terms in eq. (G.2). The averaging over the first shell leads to the factors of √ 6 in the latter result. Note that to perform this averaging one must also use the rotation invariance of K df,3 . It is also important that m = 0 in the last term in eq. (G.3), since this is the component that lives in the A + 1 irrep when the spectator momentum lies in theẑ direction.

JHEP03(2019)106
In the following, we demonstrate that eq. (G.3) holds if c = √ 5. There are three inputs needed for this demonstration. The first is the observation that the same configuration of final-state particles can contribute to both sides of eq. (G.3). To explain this we need to write both initial and final states in the form used prior to their decomposition into harmonics, so that we have K df,3 (E; p,â * ; k,â * ). Then one can show, using permutation symmetry alone, that K df,3 (E free 1 ; 0,ẑ; k,â * ) = K df,3 (E free 1 ; k Lẑ ,ẑ; k,â * ) .

(G.4)
This result holds for any term in the threshold expansion of K df,3 (or, indeed, for the entire quantity), and thus we do not include a superscript. To understand eq. (G.4), note that the three particles in the final state have momenta 0, k Lẑ and −k Lẑ . Calling 0 the spectator momentum yields the left-hand side of eq. (G.4), while calling k Lẑ the spectator yields the right-hand side. Since both choices describe the same momentum configuration, they must be equivalent due to the permutation symmetry of K df,3 . The second input is that K (2,3) df,3 is either independent of, or quadratic in,â * . This is explained in section 2.4, and is in one-to-one correspondence with the fact that only s-and d-waves contribute.
The final key input concerns angular averaging of a quadratic form: where V ij is an arbitrary tensor. In other words, the combination appearing on the lefthand side can be evaluated by settingn =ẑ. The same is trivially true for a quantity that is independent ofn. Combining the second and third key inputs, we deduce that K (2,3) df,3 (E; p, 0, 0; k,â * ) + √ 5K (2,3) df,3 (E; p, 2, 0; k,â * ) = K (2,3) df,3 (E; p,â * =ẑ; k,â * ) (G. 6) holds for any choice of E and p. Applying this to both sides of eq. (G.3), with E = E free 1 , and p = 0 for the left-hand side and p = k 1ẑ for the right-hand side, we find that eq. (G.3) with c = √ 5 is equivalent to the first key identity eq. (G.4). This establishes the desired result. This derivation will fail for terms of quartic and higher order in K df,3 , since the combination of = 0 and 2 parts that appears in eq. (G.6) will no longer allow the replacement ofâ * withẑ, implying that eq. (G.4) cannot be used. For example, considering one of the terms that arises in quartic terms, we find We have checked this numerically by decomposing the simplest of the quartic terms and finding that eq. (G.1) does not hold.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.