Three-boson spectrum in the presence of 1D spin-orbit coupling: Efimov's generalized radial scaling law

Spin-orbit coupled cold atom systems, governed by Hamiltonians that contain quadratic kinetic energy terms typical for a particle's motion in the usual Schr\"odinger equation and linear kinetic energy terms typical for a particle's motion in the usual Dirac equation, have attracted a great deal of attention recently since they provide an alternative route for realizing fractional quantum Hall physics, topological insulators, and spintronics physics. The present work focuses on the three-boson system in the presence of 1D spin-orbit coupling, which is most relevant to ongoing cold atom experiments. In the absence of spin-orbit coupling terms, the three-boson system exibits the Efimov effect: the entire energy spectrum is uniquely determined by the $s$-wave scattering length and a single three-body parameter, i.e., using one of the energy levels as input, the other energy levels can be obtained via Efimov's radial scaling law, which is intimately tied to a discrete scaling symmetry. It is demonstrated that the discrete scaling symmetry persists in the presence of 1D spin-orbit coupling, implying the validity of a generalized radial scaling law in five-dimensional space. The dependence of the energy levels on the scattering length, spin-orbit coupling parameters, and center-of-mass momentum is discussed. It is conjectured that three-body systems with other types of spin-orbit coupling terms are also governed by generalized radial scaling laws, provided the system exhibits the Efimov effect in the absence of spin-orbit coupling.


I. INTRODUCTION
Under which conditions do two, three, or more particles form weakly-bound states, i.e., bound states that are larger than the range of the two-, three-, and higherbody forces that bind the particles together? And under which conditions are the characteristics of these few-body bound states governed by underlying symmetries? These questions are of utmost importance across physics. For example, the existence of bound tetra-quark systems [1], first proposed in 1964 by Gell-Mann [2], has been challenging our understanding of QCD. The existence of the extremely weakly-bound triton has a profound effect on the nuclear chart, including the existence of larger exotic halo nuclei [3,4]. Historically, the triton has played an important role in the context of the Thomas collapse [5] and the Efimov effect [6,7], which is intimately tied to a discrete scaling symmetry of the three-body Schrödinger equation.
The three-boson system with two-body short-range interactions is considered the holy grail of few-body physics. It has captured physicists' attention since Efimov's bizarre and counterintuitive predictions in the early 70ies [6,7] and has spurred a flurry of theoretical and experimental works from nuclear to atomic to condensed matter to particle physics [8][9][10][11][12][13][14][15][16][17][18][19][20][21]. The unique scaling laws exhibited by Efimov trimers can be traced back to the existence of just one large length scale in the problem, namely the two-body s-wave scattering length. The main focus of the present work is on investigating what happens to the three-boson Efimov states in the presence of 1D spin-orbit coupling. Similar to few-body systems on the lattice [22], the 1D spin-orbit coupling introduces a parametric dependence of the relative Hamiltonian on the center-of-mass momentum. This center-ofmass momentum dependence leads, as we will show, to a modification of the lowest break-up threshold and has a profound effect on the binding energy. Despite this dependence on the center-of-mass momentum and despite the fact that the spin-orbit coupling terms depend on three additional parameters (namely, k so , Ω and δ; see below), it is argued that the three-boson system in the presence of 1D spin-orbit coupling possesses, in the zerorange limit, a discrete scaling symmetry and it is shown that the energy spectrum is described by a generalized radial scaling law.
The 1D spin-orbit coupling terms, which break the rotational symmetry, introduce an unusual single-particle dispersion. The HamiltonianĤ j of the j-th particle with mass m and momentum operatorˆ p j (with componentŝ p j,x ,p j,y , andp j,z ) is not simply given byˆ p 2 j /(2m) but includes a term that emulates a spin-1/2 particle interacting with a momentum-dependent "magnetic field" of infinite range [23][24][25][26][27], Here, I j denotes the 2x2 identity matrix that spans the spin degrees of freedom of the j-th particle, the vectorˆ σ j contains the three Pauli matricesσ j,x ,σ j,y , andσ j,z of the j-th particle, andˆ B represents the effective magnetic field,ˆ B = (Ω/2, 0, k sopj,z /m+δ/2), felt by the j-th particle. The Raman coupling Ω, detuning δ, and spin-orbit coupling strength k so , which characterize the two-photon Raman transition that couples (effectively) two hyperfine states of an ultracold atom, describe the deviations from the "normal" quadratic single-particle dispersion curves, where p j and p j,z (both without "hat") are expectation values of the corresponding operators. For large | p j |, the dispersion curves E j,± approach p 2 j /(2m). For small | p j |, in contrast, the E j,± curves deviate appreciably from p 2 j /(2m). The momenta p j are generalized momenta (sometimes also referred to as quasi-momenta) and not mechanical momenta (sometimes also referred to as kinetic momenta) [28]. Throughout this article, we frequently drop the prefix "generalized" and refer to p j as momentum vector of the j-th atom. The Hamiltonian given in Eq. (1) can also be realized by lattice shaking techniques as well as in photonic crystals and mechanical setups [27,[29][30][31].
If two-body short-range interactions are added, the modified single-particle dispersion curves can significantly alter the properties of weakly-bound two-and three-body states. This has been demonstrated extensively for two identical fermions for 1D, 2D, and 3D spin-orbit coupling [32][33][34][35][36][37][38][39][40][41] and for two identical bosons for 2D and 3D spin-orbit coupling [41][42][43][44][45] but not for the 1D spin-orbit coupling considered in this work. The present work presents the first study of how the experimentally most frequently realized 1D spin-orbit coupling terms modify the three-boson energy spectrum. We note, however, that several three-body studies for bosonic and fermionic systems with other types of spin-orbit coupling exist [46][47][48][49]. All of these earlier studies limited themselves to vanishing center-of-mass momentum. Our work, in contrast, allows for finite center-of-mass momenta.
The key objective of the present work is to show that the three-boson system in the presence of 1D spin-orbit coupling obeys a generalized radial scaling law, which reflects the existence of a discrete scaling symmetry in the limit of zero-range interactions. The scaling parameter λ 0 , λ 0 ≈ 22.694, is the same as in the absence of the spin-orbit coupling terms. The generalized radial scaling law relates the energy for a given 1/a s , k so , Ω, andδ [δ is a generalized detuning that is defined in terms of the detuning δ and the z-component of the center-of-mass momentum, see Eq. (21)] to the energy for a scaled set of parameters, namely λ 0 /a s , λ 0 k so , (λ 0 ) 2 Ω, and (λ 0 ) 2δ . Correspondingly, the term "radial" does not refer to the radius in a two-dimensional space as in the usual Efimov scenario but to the radius in a five-dimensional space. The fact that the discrete scaling symmetry "survives" when the spin-orbit coupling terms are added to the three-boson Hamiltonian with zero-range interactions can be intuitively understood from the observation that k so , Ω andδ can be thought of as introducing finite length scales into the system. In the standard Efimov scenario, a s introduces a finite length scale and the radial scaling law holds regardless of whether |a s | is larger or smaller than the size of the trimer, provided |a s | is much larger than the intrinsic scales of the underlying two-and threebody interactions. In the generalized Efimov scenario considered here, the parameters a s , k so , Ω, andδ each introduce a finite length scale. Correspondingly, the generalized radial scaling law holds regardless of whether these length scales are larger or smaller than the size of the trimer, provided the length scales are much larger than the intrinsic scales of the underlying two-and threebody interactions.
Our findings for the experimentally most frequently realized 1D spin-orbit coupling are consistent with Ref. [48]. References [46,48] considered an impurity with 3D spin-orbit coupling that interacts with two identical fermions that do not feel any spin-orbit coupling terms and interact with the impurity through short-range two-body potentials. Restricting themselves to vanishing center-of-mass momenta, Ref. [46] stated that the trimers for mass ratio 13.6 "no longer obey the discrete scaling symmetry even at resonance" because the spinorbit coupling "introduces an additional length scale". In Ref. [48], the same authors arrive at a seemingly different conclusion, namely "in the presence of SO [spin-orbit] coupling, the system exhibits a discrete scaling behavior" and "the scaling ratio is identical to that without SO [spin-orbit] coupling". The two statements can be reconciled by noting that the discrete scaling symmetry requires an enlarged parameter space, an aspect that was recognized in Ref. [48] but not in Ref. [46]. We conjecture that the discrete scaling symmetry holds for any type of spin-obit coupling and all center-of-mass momenta. Depending on the type of the spin-orbit coupling, the generalized Efimov plot is four-or five-dimensional and the generalized radial scaling law applies to the entire lowenergy spectrum. The dependence of the energy levels on the system parameters has to be calculated explicitly once for each type of spin-obit coupling.
The remainder of this article is organized as follows. To set the stage, Sec. II reviews the standard Efimov scenario for three identical bosons. Section III introduces the system Hamiltonian in the presence of 1D spin-orbit coupling and discusses the associated continuous and discrete scaling symmetries. The generalized radial scaling law for the three-boson system in the presence of 1D spin-orbit coupling is confirmed numerically in Sec. IV. Section V highlights the role of the center-of-mass momentum and discusses possible experimental signatures of this dependence. Finally, Sec. VI presents an outlook. Technical details are relegated to several appendices.

II. REVIEW OF STANDARD EFIMOV SCENARIO
The relative Hamiltonian for two identical bosons of mass m interacting through the zero-range contact inter-action V 2b,zr ( r), where a s denotes the two-body s-wave scattering length and r the internuclear distance vector (r = | r|), possesses a continuous scaling symmetry [8]. Performing the transformation a s → λa s , r → λ r, and t → λ 2 t, where t denotes the time and λ a real number (scaling parameter), the relative two-body time-dependent Schrödinger equation remains unchanged. Importantly, the continuous scaling symmetry extends to three identical mass m bosons with position vectors r j that interact through pairwise s-wave zero-range interactions V 2b,zr ( r jk ) [8]. To see this, we consider the timedependent Schrödinger equation for the relative threebody HamiltonianĤ rel , where ρ j denotes the j-th relative Jacobi vector and µ j the associated Jacobi mass. We use a "K-tree" (see Appendix A) in which µ 1 for the two-body system is given by m/2 and µ 1 and µ 2 for the three-body system are given by m/2 and 2m/3. The zero-range three-body potential V 3b,zr (R), is written in terms of a six-dimensional delta-function in the three-body hyperradius R, R 2 = r 2 12 + r 2 13 + r 2 23 . Since the coupling constant g 3 has units of length 4 , it can be rewritten as g 3 = Cκ −4 * , where C is a real constant and κ * the three-body binding momentum of one of the three-boson bound states at unitarity (infinite a s ). While V 3b,zr (R) has to be regularized in practice, the explicit regularization is irrelevant for our purpose. Performing the transformation a s → λa s , r jk → λ r jk , t → λ 2 t, and κ * → λ −1 κ * , (7) the Schrödinger equation for the relative Hamiltonian given in Eq. (5) remains unchanged, i.e., the three-body system possesses a continuous scaling symmetry.
Intriguingly, the three-body system with zero-range interactions additionally exhibits an exact discrete scaling symmetry [8]. The discrete transformation is given by and κ * → κ * , (8) where n = ±1, ±2, · · · , ±∞ and λ 0 ≈ 22.694. The discrete scaling transformation, which underlies the threebody Efimov effect, is illustrated in Fig. 1(a). Fixing the three-body parameter κ * [see Eq. (8)], the Efimov plot depicts K as a function of 1/a s , where and E denotes the eigen energy of the HamiltonianĤ rel given in Eq. (5). The thick solid line in Fig. 1(a) shows K for one of the three-body eigen energies. The thick solid line merges with the three-atom threshold on the negative a s -side and with the atom-dimer threshold (dashed line) on the positive a s -side. The thick solid line is obtained by solving the time-independent Schrödinger equation for the three-body HamiltonianĤ rel . Provided the thick solid line is known (a parametrization can be found in Refs. [8,11]), the thin solid lines-which correspond to other three-body eigen energies-can be obtained using the discrete scaling symmetry without having to explicitly solve the Schrödinger equation again. For the construction, it is convenient to switch from the vector y = (1/a s , K) T to a radius y = | y| and an angle ξ, and (a s ) −1 = y cos ξ, where ξ goes from π/4 to π. The limits π/4 and π are set by the atom-dimer and three-atom thresholds, respectively. To obtain the thin solid lines in Fig. 1(a) from the thick solid line, one fixes the angle ξ and reads off the values of the pair (1/a s , K) corresponding to the solid line. Using it can be seen that the discrete scaling transformation a s → (λ 0 ) n a s and E → (λ 0 ) −2n E implies y → (λ 0 ) −n y. Thus, dividing the radius y of the thick solid line by (λ 0 ) ±1 , (λ 0 ) ±2 , · · · and using the scaled value of y in Eqs. (10) and (11), one obtains the values of the vectors y = (1/a s , K) T corresponding to the thin solid lines. This construction, referred to as Efimov's radial scaling law, is a direct consequence of the discrete scaling symmetry. If the three-boson system is characterized by κ new * instead of κ * , the entire energy spectrum is scaled, i.e., if y = (1/a s , K) T describes a point on the Efimov plot for κ * , then (κ new * /κ * ) y describes a point on the Efimov plot for κ new * .

III. SYMMETRIES IN THE PRESENCE OF 1D SPIN-ORBIT COUPLING
This section generalizes the symmetry discussion presented in the previous section to the two-and three-boson systems in the presence of 1D spin-obit coupling. As a first step, we derive the relative two-and three-body Hamiltonian with zero-range interactions in the presence To make this plot, λ0 has been artificially set to 2 instead of 22.694. The dashed line shows the atom-dimer threshold. The thin radially outgoing solid lines and arrows illustrate the scaling law. Circles and squares mark the critical scattering lengths a− at which the trimer energy is degenerate with the three-atom threshold and the critical scattering lengths a * at which the trimer energy is degenerate with the atom-dimer threshold, respectively. (b) Collapse of neighboring energy levels for the finite-range interaction model [Ĥ rel in Eq. (5) with V 2b,zr and V 3b,zr replaced by V 2b,G and V 3b,G , respectively; R0 = √ 8r0 and (κ * ) −1 ≈ 66.05r0]. The solid line shows the fourth-root of the energy of the lowest three-boson state as a function of the square-root of the inverse of the s-wave scattering length. The dashed line shows the associated atom-dimer threshold. The dots show the energy of the second-lowest three-boson state, with the radial scaling law applied in reverse so as to collapse the second-lowest level (dots) onto the lowest level (solid line). For clarity, the scaled atom-dimer threshold for the second-lowest three-boson state is not shown.
of 1D spin-orbit coupling. In a second step, it is shown that these systems possess a continuous scaling symmetry. In a third step, it is argued that the three-boson system additionally exhibits a discrete scaling symmetry, suggesting the existence of a generalized radial scaling law. Numerical evidence that supports our claim that the three-boson system with 1D spin-orbit coupling is governed by a generalized radial scaling law is presented in Sec. IV.
We start with the first step. The N -boson Hamiltonian in the presence of 1D spin-orbit coupling readŝ where the non-interacting and interacting pieces are given bŷ and Here, I j,··· ,k with j < k spans the spin degrees of freedom of particles j through k, I j,··· ,k = I j ⊗ · · · ⊗ I k . For N = 2, only V 2b contributes. For N = 3, r jkl is equal to the three-body hyperradius R. The interaction model considered throughout this work assumes that the interactions are the same in all spin channels. It is, just as in the case without spin-orbit coupling, convenient to use Jacobi coordinates ρ j and associated momentum operatorsˆ q j instead of the single-particle quantities r j andˆ p j . Importantly, the N -th Jacobi "quantities" ρ N andˆ q N correspond to the center-ofmass vector and center-of-mass momentum operator. It can be shown straightforwardly that the HamiltonianĤ commutes with the center-of-mass momentum operator q N [50], i.e., the Schrödinger equationĤΨ = EΨ can be solved for each fixed q N . Using this and Jacobi coordinates, the non-interacting fixed-q N Hamiltonian, denoted byĤ ni , readsĤ whereĤ The explicit form of the operatorsΣ j,z with j = 1, · · · , N − 1 is given in Appendix A. Note that the first and second lines of Eq. (17) contain momentum operators while the fourth line of Eq. (17) and the second term on the right hand side of Eq. (16) contain expectation values of the center-of-mass momentum operators (and not operators). As "usual", the interactionV int depends on ρ 1 , · · · , ρ N −1 but not on the center-of-mass vector ρ N . This implies that the eigen states Ψ can be written as where [51] and where the Φ rel , which are eigen states ofĤ rel , depend on the Jacobi vectors ρ 1 , · · · , ρ N −1 and the spin degrees of freedom. Equation (17) shows that the eigen energies ofĤ depend on the generalized detuningδ, i.e., q N,z and δ enter as a combination and not as independent parameters. This observation suggests that the center-of-mass momentum may play a decisive role in determining the characteristics of the weakly-bound two-and three-body states (see also Refs. [35,37]). The parametric dependence of the HamiltonianĤ rel on the z-component q N,z of the center-of-mass momentum is a direct consequence of the fact that the presence of the spin-orbit coupling breaks the Galilean invariance [26]. One immediate consequence of the broken Galilean invariance is that knowing the energy of an eigen state with q N,z = 0 does not, in general, suffice for predicting the energy of an eigen state with q N,z = 0. Importantly, the eigen states Ψ depend, in general, explicitly on q N,z and δ and not just onδ.
We are now ready to address the second step. Parametrizing the two-body interactions V 2b by the zerorange potential V 2b,zr , the N = 2 relative Hamiltonian depends on four parameters, namely a s , k so , Ω, andδ. It can be readily checked that the corresponding timedependent Schrödinger equation is invariant under the transformation r → λ r, and t → λ 2 t, (22) i.e., the N = 2 system with zero-range interactions possesses a continuous scaling symmetry. The continuous scaling symmetry extends to the three-boson system with zero-range interactions [V 2b = V 2b,zr and V 3b = V 3b,zr in Eq. (15)] in the presence of spin-orbit coupling since the corresponding time-dependent N = 3 Schrödinger equation is invariant under the transformation Equations (22) and (23) generalize Eqs. (4) and (7) from Sec. II.
Paralleling the discussion of Sec. II, step three poses the question whether or not the three-boson system in the presence of spin-orbit coupling additionally possesses a discrete scaling symmetry in the zero-range interaction limit. Our claim is that it does and that the discrete transformation is given by where λ 0 is identical to the scaling factor of the standard Efimov scenario, i.e., λ 0 ≈ 22.694. Since no general analytical solutions exist to the three-boson Schrödinger equation in the presence of spin-orbit coupling, we rely on numerics to support our claim. The claim that the discrete scaling symmetry survives in the presence of the spin-orbit coupling terms can be understood intuitively by realizing that the spin-orbit coupling terms modify the low-but not the high-energy portions of the singleparticle dispersion curves. To set the stage for the numerical calculations presented in the next section, we discuss a number of consequences of the discrete scaling symmetry. The discrete scaling symmetry suggests a generalized radial scaling law for the three-boson system in the presence of 1D spin-orbit coupling in which the Efimov plot for y = (1/a s , K) T discussed in the previous section is replaced by a generalized Efimov plot for In the limit that the second, third, and fourth parameters vanish, each of the usual Efimov energies is fourfold degenerate due to the fact that the spin degrees of freedom enlarge the three-boson Hilbert space by a factor of four (from the 2 3 = 8 independent spin configurations, one can construct four fully symmetric spin functions). For non-vanishing k so , Ω, andδ, we expect that the three-boson system supports four "unique" energy levels. Each of the four energies, collectively referred to as a manifold, is characterized by a vector y.
Knowing the dependence of each of these energy curves on 1/a s , k so , Sign(Ω) m|Ω|/ 2 , and Sign(δ) m|δ|/ 2 , there should exist other energy manifolds for the same κ * that can be obtained from the manifold that has been mapped out without explicitly solving the three-boson Schrödinger equation again. To see how, we switch from the five parameters given in Eq. (25) to the length y = | y| and four angles ξ 1 , · · · , ξ 4 for each of the four energy levels in the "reference manifold", and 1/a s =y cos ξ 4 .
The full range of possible a s , k so , Ω, andδ is covered if The range of the angles is further constrained by the energy surfaces of the three-atom and atom-dimer thresholds (see Secs. IV and V). To obtain the K for other manifolds, one chooses a direction of the vector y by fixing the angles ξ 1 to ξ 4 and reads off the values of the components of y for each of the four known energy levels. Using it can be seen that the discrete transformation a s → Thus, dividing the "hyperradius" y corresponding to the j-th energy in the reference manifold by (λ 0 ) ±1 , (λ 0 ) ±2 , · · · and using the scaled value of y in Eqs. (26)- (30), one obtains the values of the components of y for the j-th energy level in the other manifolds. The generalized scaling law is tested in the next section by considering two neighboring energy manifolds and confirming that the energy manifolds collapse onto each other if the discrete scaling transformation is applied to the energy levels in the more weaklybound manifold.

IV. NUMERICAL TEST OF THE GENERALIZED RADIAL SCALING LAW
To facilitate the numerical calculations, we replace the two-body zero-range potential V 2b,zr by an attractive Gaussian V 2b,G with range r 0 and depth v 0 , where v 0 is negative and adjusted such that V 2b,G (r jk ) supports at most one two-body s-wave bound state. To reduce finite-range effects, we aim to work in the regime where the absolute value of the free-space s-wave scattering length a s is notably larger than r 0 . Parameter combinations where the absolute value of the free-space pwave scattering volume is large are excluded. The threebody zero-range potential V 3b,zr is replaced by a repulsive Gaussian V 3b,G with range R 0 and height V 0 , In our numerical calculations, R 0 is fixed at √ 8r 0 and V 0 (V 0 ≥ 0) is varied to dial in the desired three-body parameter κ * . Specifically, we define κ * to be the binding momentum of the energetically lowest-lying universal three-body state at unitarity (infinite a s ) for k so = Ω = δ = 0. Without the repulsive three-body potential, the lowest three-body state is not universal [52]. The repulsive three-body potential pushes the lowest three-body energy up and we adjust V 0 , for fixed v 0 (infinite a s ), such that the energy of the lowest three-body state for finite V 0 is identical to the energy of the first excited three-body state for V 0 = 0. This corresponds to (κ * ) −1 ≈ 66.05r 0 , i.e., the trimer is much larger than the intrinsic scales of the two-and three-body interactions. With the repulsive three-body potential turned on, the radial scaling law can be tested using the two lowest-lying energy manifolds.
We start the discussion of our numerical results by looking at the three-body spectrum for the standard Efimov scenario (k so = Ω =δ = 0). The reason for discussing this "reference system" is two-fold: it illustrates how to check the validity of the radial scaling law for a case where it is known to hold and it gives us a sense for the finite-range corrections expected in the presence of spin-orbit coupling. The solid line in Fig. 1(b) shows the relative three-body energy of the energetically lowestlying state as a function of the inverse of the s-wave scattering length for the Hamiltonian given in Eq. (5) with V 2b,zr and V 3b,zr replaced by V 2b,G and V 3b,G , respectively. To "compress" the data, the horizontal and vertical axis employ a square-root and fourth-root representation. The scattering length is scaled by r 0 and the energy by E sr , The trimer energy merges with the three-atom threshold on the negative scattering length side at r 0 /|a s | ≈ 0.01. To get a feeling for the finite-range effects, we assume that the radial scaling law holds and apply it "in reverse". Specifically, using numerically determined pairs (1/a s , K) corresponding to the excited state, the dots in Fig. 1(b) show the points (λ 0 r 0 /a s , λ 0 r 0 K), using-as for the lowest state-the square-root and fourth-root depiction. In the zero-range limit (r 0 → 0 and R 0 → 0), the dots would lie on top of the solid line. The nearly perfect agreement between the solid line and the dots in Fig. 1(b) indicates that the finite-range effects are negligibly small for the parameter combinations considered.
To test the generalized radial scaling law proposed in Sec. III, we calculate the eigen energies of states in the lowest and second-lowest manifolds ofĤ rel (there are at most four states in each manifold) and scale the energies in the second-lowest manifold assuming that the generalized radial scaling law holds. If the energy curves collapse, the generalized radial scaling law is validated.
In the presence of the 1D spin-orbit coupling, the generalized Efimov plot has five axes. Clearly, visualizing energy surfaces that depend on four parameters is impossible and fully mapping out these high-dimensional dependences is computationally demanding. Thus, we consider selected cuts in the five-dimensional space. Our first cut uses (k so ) −1 = 50r 0 , Ω = 2E so = 0.04E sr , where andδ = 0. For these parameters, we calculate the relative energy E of the states in the lowest energy manifold. The solid lines in Fig. 2(a) show the quantity −|(E − E aaa th )/E sr | 1/4 as a function of Sign(a s )|r 0 /a s | 1/2 , where E aaa th denotes the energy of the lowest three-atom threshold whose wave function has the same total momentum q 3,z along the z-axis as the three-body system. The determination of E aaa th is discussed in Appendix C. The energy E aaa th is independent of κ * and referencing E relative to the lowest three-atom threshold does not alter the generalized radial scaling law. Figure 2(a) shows that the lowest energy manifold consists of, depending on the value of r 0 /a s , zero, one, or two energy levels [the second and third excited states of the lowest manifold exist at larger r 0 /a s than those shown in Fig. 2(a)].
The lowest three-body energy merges with the threeatom threshold on the negative s-wave scattering length side and with the atom-dimer threshold [dashed line in Fig. 2(a)] on the positive scattering length side. The determination of the atom-dimer threshold energy E ad th is discussed in Appendix D. Just as the three-atom threshold, the atom-dimer threshold is independent of κ * and determined such that the momentum q 3,z of the atomdimer system is the same as that of the three-body system. The second lowest state does not merge with the three-atom threshold on the negative a s side but with the atom-dimer threshold on the positive a s side.
Having determined the energies of the states in the lowest energy manifold, the next step is to calculate the energies of the states in the second-lowest energy manifold. To map the energies of the states in the secondlowest manifold onto the energies of the states in the lowest energy manifold, we use the same r 0 , R 0 , and κ * and calculate the energies of the states in the secondlowest energy manifold for a k so that is λ 0 times smaller than the k so used to calculate the energies of the states in the lowest energy manifold [i.e., for (k so ) −1 ≈ 1, 135r 0 ], for a Ω that is (λ 0 ) 2 times smaller than the Ω used to calculate the energies in the lowest energy manifold (i.e., for Ω ≈ 7.77 × 10 −7 E sr ), and forδ = 0 (the scaling does not change zero) as a function of r 0 /a s . Hav- In all three panels, the dashed lines show the atom-dimer threshold. The dots show the energies of states in the secondlowest manifold, with the generalized radial scaling law applied in reverse so as to collapse the three-body energies of states in the second-lowest manifold (dots) onto the threebody energies of the states in the lowest manifold (solid lines). For clarity, the scaled atom-dimer thresholds for the secondlowest energy manifold are not shown in any of the panels.
ing calculated the energies of the states in the secondlowest manifold for the scaled k so , Ω, andδ, the pairs (1/a s , E − E aaa th ) are scaled (note that E aaa th for the excited state manifold is calculated using the scaled k so , Ω, andδ values). The dots in Fig. 2(a) show the scaled pairs (λ 0 r 0 /a s , −(λ 0 ) 2 |E − E aaa th |/E sr ), using-as for the lowest energy manifold-the square-root and fourth-root depiction. It can be seen that the solid lines and dots agree very well. Note that the atom-dimer threshold for the second-lowest energy manifold also needs to be recalculated using the scaled k so , Ω, andδ [the resulting energies lie essentially on top of the dashed line and are not shown in Fig. 2(a)]. The deviation between the solid line and dots is 0.025 % for (r 0 /a s ) 1/2 = 0 and 0.80 % for (r 0 /a s ) 1/2 = 0.19. These deviations are comparable to those between the corresponding atom-dimer thresholds [in this case, the deviations are 0.023 % for (r 0 /a s ) 1/2 = 0 and 0.82 % for (r 0 /a s ) 1/2 = 0.19]. We conclude that our numerical results are consistent with the generalized radial scaling law.
We emphasize that the scaling law has to be applied to all five axes of the generalized Efimov plot, i.e., to obtain the dots in Fig. 2(a) it is imperative to not only scale the two axes depicted but also the parameters corresponding to the three axes that are not depicted. The ratio of the lowest energy in neighboring manifolds at unitarity, e.g., is only equal to 22.694 2 if the direction ofŷ is the same for the two energy levels under consideration.
The energy scales E so , |Ω|, and |δ| are much smaller than |E − E aaa th | for a large portion of Fig. 2(a). The region close to the three-atom threshold is an exception. As such it might be argued that the spin-orbit coupling terms are too weak to notably influence the energy spectrum, possibly suggesting that the applicability of the generalized radial scaling law is trivial. One fact that speaks against this argumentation is that the shape of the energy levels is notably influenced by the spin-orbit coupling terms. This is, e.g., reflected by the fact that the energy levels in a given energy manifold are not degenerate. To more explicitly demonstrate that the generalized radial scaling law holds when one or more of the energy scales associated with the spin-orbit coupling terms is/are larger than the binding energy, we repeat the calculations for larger k so than those used in Fig. 2(a). Specifically, to determine the energy of the lowest state in the lowest energy manifold [solid line in Fig. 2(b)], we use (k so ) −1 = 25r 0 while keeping r 0 , R 0 , κ * , andδ unchanged. The Raman coupling strength Ω is set to be equal to 2E so . To demonstrate the collapse of the energies of the lowest states in the second-lowest and lowest manifolds, we apply the generalized radial scaling law in the same way as in Fig. 2(a). The energy of the lowest state in the second-lowest manifold is shown by dots in Fig. 2(b). The agreement with the solid line is excellent, supporting our claim that the generalized radial scaling law is not limited to the case where the energy scales associated with the spin-orbit coupling are smaller than the binding energy of the trimer, provided these energies are notably smaller than E sr .
To show the characteristics of the excited states in the lowest manifold in more detail, we consider a smaller k so , (k so ) −1 = 100r 0 , and as beforeδ = 0 and Ω = 2E so . The use of a smaller k so [solid lines in Fig. 2(c)] moves the merging points of the three-body energies corresponding to the excited states with the atom-dimer threshold to the left compared to Fig. 2(a). Again, scaling the parameters appropriately, the dots in Fig. 2(c) show the energies of the states in the second-lowest energy manifold. The dots agree nearly perfectly with the solid lines not only for the lowest state in the two manifolds but also for the excited states in the two manifolds, lending strong numerical support for the validity of the generalized ra-dial scaling law and hence for the existence of the discrete scaling symmetry in the presence of 1D spin-orbit coupling terms in the zero-range limit.
As already mentioned, the three-body parameter κ * for V 0 = 0, defined using the energy of the first excited state at unitarity in the absence of spin-orbit coupling, is identical to the κ * for the three-body interaction with finite V 0 used throughout this section. Turning on the spinorbit coupling, we checked that the energies of states in the second-lowest manifold for V 0 = 0 agree well with the energies of states in the lowest manifold for the finite V 0 . This provides evidence that the generalized radial scaling law is, just as the standard radial scaling law, independent of the details of the underlying microscopic interaction model. To confirm the continuous scaling symmetry of the three-body Hamiltonian in the presence of 1D spinorbit coupling, we checked that the energies for different κ * can be mapped onto each other: If y describes a point on the Efimov plot for κ * , then (κ new * /κ * ) y describes a point on the Efimov plot for the new κ new * .

V. EXPERIMENTAL IMPLICATIONS: ROLE OF CENTER-OF-MASS MOMENTUM
Measuring signatures associated with two consecutive trimer energy levels is challenging, especially for equalmass bosons, due to the relatively large discrete scaling factor of 22.694. The reason is that the absolute value of the scattering length should, on the one hand, be notably larger than the van der Waals length r vdW and, on the other hand, be smaller than the de Broglie wave length λ dB [8,11]. Despite these challenges, the discrete scaling symmetry underlying the standard Efimov scenario has been confirmed experimentally by monitoring the atom losses of an ultracold thermal gas of Cs atoms as a function of the s-wave scattering length [20] (for unequal mass mixtures, see Refs. [57,58]). When the trimer energy is degenerate with the three-atom threshold [dots in Fig. 1(a); the corresponding critical scattering lengths are denoted by a (n) − ] or with the atom-dimer threshold [squares in Fig. 1(a); the corresponding critical scattering lengths are denoted by a (n) * ], the losses are enhanced. Since the critical scattering lengths for consecutive trimer states are related to the scaling factor λ 0 , these atom-loss measurements provide a direct confirmation of the discrete scaling symmetry. In addition, other characteristics of the standard Efimov scenario have been measured [8,[11][12][13]. For example, the critical scattering lengths a can be viewed as a test of the functional form of the energy levels shown in Fig. 1(a). Other experimental tests of the standard Efimov scenario include the determination of the binding energy of an Efimov trimer via radiofrequency spectroscopy [17,18], the imaging of the quan- tum mechanical density of the helium Efimov trimer via Coulomb explosion [21], and the observation of four-and five-body loss features that are universally linked to the critical scattering lengths of the Efimov trimer [53][54][55][56]. Directly measuring the discrete scaling symmetry in the presence of spin-orbit coupling requires varying the inverse of the s-wave scattering length by the scaling factor λ 0 , as in the standard Efimov scenario, as well as varying the spin-orbit coupling parameters Ω, E so , and δ by (λ 0 ) 2 . Covering such a wide range of parameters is expected to be very challenging experimentally. In what follows we instead focus on the situation where the spinorbit coupling parameters k so and Ω are held fixed while the s-wave scattering length a s and generalized detuning δ are varied. An analogous study for the standard Efimov scenario would look at the three-boson system for a fixed finite s-wave scattering length. In this case, the energy spacing would not be (λ 0 ) 2 ; however, the energies of neighboring states would still be uniquely related to each other.
For concreteness, we consider the 133 Cs system [59], for which the three-atom resonances in the absence of spin-orbit coupling occur at the critical scattering lengths a  [20]. Here, a 0 denotes the Bohr radius and the superscripts "(0)" and "(1)" indicate that these critical scattering lengths are for the ground and first excited Efimov trimers, respectively. Applying our numerical result κ * a − = −1.505 to the first excited state, the Cs system is characterized by (κ * ) −1 ≈ 13, 416a 0 . Figure 3 shows the negative of the binding energy of the lowest state in the first excited manifold for k so /κ * ≈ 1.32 and Ω = 2E so as functions of the inverse of the swave scattering length and the generalized detuningδ using, as in the previous sections, that the scattering lengths are the same for all spin channels. Using Cs's a (1) − , these parameters correspond to (k so ) −1 ≈ 10, 156a 0 , E so /h ≈ 0.132kHz, and Ω/h ≈ 0.264kHz. Comparison with the 87 Rb experiment at NIST [28], which uses (k so ) −1 ≈ 3, 410a 0 (corresponding to E so /h ≈ 1.786kHz) and Ω/h values ranging from zero to about 10kHz, suggests that the parameter regime covered in Fig. 3 is reasonable. Figure 3 shows that the three-boson binding energy for a fixed scattering length is largest forδ = 0 (this is where the three-atom threshold has a degeneracy of six; see Appendix C). In addition, there exists an enhancement of the binding forδ/h ≈ 0.301kHz (this is where the three-atom threshold has a degeneracy of four; see Appendix C). Asδ goes to infinity, the trimer in the presence of the 1D spin-orbit coupling becomes unbound at the same scattering length as the corresponding trimer in the absence of spin-orbit coupling (i.e., at a s ≈ −20, 190a 0 ). Fig. 3 is calculated by enforcing that the three-boson threshold has the same center-of-mass momentum as the trimers (see Appendices C-E). If the detuning δ is equal to zero, the generalized detuningδ is directly proportional to the z-component q 3,z of the center-of-mass momentum [see Eq. (21)]. In this case, the trimer is bound maximally for q 3,z = 0. However, for finite detuning δ, the most strongly bound trimer has a finite center-of-mass momentum. A similar dependence on the center-of-mass momentum was pointed out in Refs. [35,37,60] for the two-fermion system.

The three-boson binding energy shown in
The dependence ofĤ rel on q N =3,z is a key characteristic of systems with 1D spin-orbit coupling. A similar dependence exists for three-body systems in the presence of 2D or 3D spin-orbit coupling (in these cases, the relative Hamiltonian depends on two or all three components of q N =3 ) and for three-body systems on a lattice (in this case, q N =3 is a lattice or quasi-momentum vector). In all works known to us [16,[46][47][48], the assumption q N =3 = 0 is made prior to obtaining concrete results. Table I contrasts studies for systems, which possess a center-of-mass momentum dependence, with the "standard" three-boson Efimov system (first row), for which the relative Hamiltonian is independent of q N =3 . In the standard Efimov case, the lowest atom-dimer threshold of the relative Hamiltonian is given by the energy E a of an atom with vanishing atom momentum vector q a (E a is equal to zero) plus the energy E d ( q d = 0) of a dimer with vanishing dimer momentum vector q d . When the relative Hamiltonian depends on q N =3 , the atom-dimer threshold needs to be determined carefully, since the trimer with fixed q N =3 can break up into an atom with finite momentum and into a dimer with finite momentum in such a way that the generalized three-body center-of-mass momentum is conserved. Of the many break-up configu-systemĤ rel =Ĥ rel ( qN=3)? restriction? atom-dimer threshold (rel. Ham.) 3-spinless bosons; "standard" Efimov scenario [6] no no E d ( q d = 0) + Ea( qa = 0) FFX; X feels 2D SOC; Borromean binding [47] yes qN=3 = 0 min q d E d ( q d ) + min qa Ea( qa) FFX; X feels 3D SOC; universal/Efimov trimers [46,48] yes qN=3 = 0 min q d + qa= q N =3 [E d ( q d ) + Ea( qa)] BBB quasi-particles on lattice; Efimov trimers [16] yes TABLE I: Summary of three-particle studies. The standard Efimov scenario (first row) is contrasted with three-particle systems for which the relative HamiltonianĤ rel depends parametrically on the generalized three-body center-of-mass momentum vector qN=3. The last column lists the atom-dimer threshold of the relative Hamiltonian. The symbols E d , Ea, q d , and qa denote the energy of the dimer, energy of the atom, generalized momentum of the dimer, and generalized momentum of the atom, respectively. "SOC" stands for "spin-orbit coupling", "F" for "fermion", "X" for a particle different from "F", and "B" for "boson".
rations that conserve the three-body center-of-mass momentum, the one with the lowest energy defines the atomdimer threshold. Table I shows that the definition of the lowest atom-dimer threshold of the relative Hamiltonian varies in the literature. The definitions employed in Refs. [16,47] disagree with the definition used in the present work (last row of Table I). While the definition of Ref. [47] may be meaningful in a many-body context (see also comment [60]), we fail to see how the definition of Ref. [16] can, in general, be correct.
It is proposed that the center-of-mass momentum dependence can be observed experimentally by performing atom-loss measurements for fixed k so , Ω, and δ on a cold thermal atomic gas. Tuning the s-wave scattering length, one expects-just as in the case where the spin-orbit coupling is absent-enhanced losses when the trimer energy is degenerate with the three-atom threshold. However, in contrast to the standard Efimov scenario, such a degeneracy exists for a range of scattering lengths provided the trimers embedded in the thermal gas have different threebody center-of-mass momenta (the exact distribution of center-of-mass momenta is set by the temperature of the gas sample). Figure 3 shows that the critical scattering length a (1) − changes, for the Cs example, from −20, 190a 0 for largeδ to −7, 791a 0 forδ = 0. Provided the threebody center-of-mass momenta are spread over the range covered on the vertical axis in Fig. 3, one expects enhanced losses over the entire scattering length window. The difference between the losses in the presence and absence of the spin-orbit coupling terms can be interpreted as a few-body probe of the breaking of the Galilean invariance in the presence of spin-orbit coupling.
An important question is whether the changes of the loss features related to the lowest state in the secondlowest manifold will be washed out by finite temperature effects. A comprehensive answer to this question will require performing three-body recombination calculations, which include thermal averaging, in the presence of spin-orbit coupling. Such calculations are beyond the scope of this work. Given that the energy scales associated with the spin-orbit coupling are, for the example considered in Fig. 3, comparable to 2 κ 2 * /m and that the binding energy of the trimer near the three-atom threshold is much smaller than 2 κ 2 * /m, we are hopeful that the temperatures realized in previous Cs experiments (T ≥ 7.7nK) [20] are low enough to observe the impact of the spin-orbit coupling terms on the loss features. For example, without spin-orbit coupling, the loss coefficient L 3 is maximal around −20, 000a 0 and reaches about half of its maximum value at around −10, 000a 0 [see Fig. 1(a) of Ref. [20]]. In the presence of spin-orbit coupling, the loss feature is expected to be centered over the range −20, 190a 0 to −7, 790a 0 , leading to an observable modification of the shoulder on the less negative scattering length side. The shape of the shoulder is expected to carry a signature of the non-monotonic dependence of the critical scattering length a − , associated with the lowest state in the lowest manifold, displays essentially no dependence onδ, i.e., the associated three-atom loss feature should only be minimally affected by the spin-orbit coupling terms. Intuitively, this can be understood by realizing that the energy scales associated with the spin-orbit coupling parameters are much smaller than the binding energy of the lowest-lying trimer state. The fact that the loss features for the lowest state in the lowest and second-lowest manifolds are expected to be very different can also be understood from the generalized radial scaling law. Fixing the spin-orbit coupling parameters corresponds to looking at particular cuts in the fivedimensional parameter space as opposed to looking along a specific radial direction. As a consequence, the loss features for the two manifolds can be very different even if the scattering lengths at which the loss features occur are, roughly, spaced by λ 0 .

VI. CONCLUDING REMARKS
This work analyzed what happens to the three-boson Efimov spectrum if 1D spin-orbit coupling terms, realizable in cold atoms as well as in photonic crystals and mechanical setups, are added to the Hamiltonian. The spinorbit coupling terms introduce a parametric dependence of the relative Hamiltonian on the center-of-mass momentum vector. A similar center-of-mass momentum vector dependence exists for few-body systems with short-range interactions on a lattice. The present work mapped out, for the first time, the three-boson spectrum as a function of the center-of-mass momentum vector. It was found that the three-boson system in the presence of 1D spinorbit coupling obeys a generalized radial scaling law in a five-dimensional parameter space, which is associated with a discrete scaling symmetry. Within the framework of effective field theory, the existence of the discrete scaling symmetry can be rationalized by scale separation: The discrete scaling symmetry of the standard Efimov scenario "survives" provided the additional length scales are much larger than the ranges of the intrinsic interactions. While our work focused on 1D spin-orbit coupling, the discrete scaling symmetry should persist for other types of spin-obit coupling schemes as well.
The spin degrees of freedom lead-for the type of spinorbit coupling considered in this work-to a quadrupling of each Efimov trimer (manifold of four states). The three-body states in a given manifold are tied to one of the three two-boson states [61]. The point (1/a s , k so , Ω,δ) = (0, 0, 0, 0) serves as an accumulation point for all four states of the manifold, i.e., in its vicinity, there exist infinitely many three-body bound states. The rich structure of two-and three-boson states should be amenable to experimental verification. Due to the dependence of the trimers on the center-of-mass momentum, the scattering length at which the lowest trimer in the second-lowest manifold merges with the atom-atomatom threshold is, in fact, a scattering length window. Similar scattering length windows exist for the excited states in the second-lowest manifold. It was argued that these scattering length windows should be observable in cold atom loss experiments, providing a direct few-body signature of the breaking of the Galilean invariance of systems with spin-orbit coupling.
If one considers a cut in the generalized fivedimensional Efimov plot, energy levels are not spaced by the scaling factor (λ 0 ) 2 . Let us consider the situation where a s is infinitely large and where k so , Ω, andδ are finite. In this case, the low-energy scales associated with the 1D spin-orbit coupling terms lead to a cut-off of the hyperradial −1/R 2 Efimov potential curve (R denotes the three-body hyperradius). As a consequence, the number of three-body bound states at unitarity is not infinitely large. Albeit due to a different mechanism, this is similar in spirit to the disappearance of Efimov states if an Efimov trimer is placed into a gas of bosons or fermions [62,63]. This is also similar in spirit to a rather different system, namely the H − ion. Taking only Coulomb interactions into account, one obtains a −1/r 2 attraction [64], where r is the distance between the extra electron and the atom. Relativistic effects introduce an additional length scale, which renders the number of bound states finite [64].
The calculations in this work were performed assuming that the interactions between the different spin-channels are all equal. If one of the scattering lengths is large and tunable while the others are close to zero, the discrete scaling symmetry should still hold (approximately). To find the functional form of the energies for this scenario, the spectrum has to be recalculated.
The study presented should be viewed as a first step toward uncovering the rich three-and higher-body physics that emerges as a consequence of the unique coupling between the relative and center-of-mass degrees of freedom in cold atom systems in the presence of artificial Gauge fields. While somewhat different in nature, the coupling of these degrees of freedom in the relativistic Klein Gordon and Dirac equations and quantum field theories has captured physicists' imagination for many decades. for the equal-mass two-particle system and for the equal-mass three-particle system. The transformation matrix U also defines the matricesΣ j,z (j = 1, · · · , N − 1). For the two-body system, we havê For the three-body system, we havê Appendix B: Explicitly correlated basis set expansion approach This appendix discusses our approach to solving the few-particle time-independent Schrödinger equation using a basis set expansion in terms of explicitly correlated Gaussian basis functions, which contain non-linear variational parameters that are optimized semi-stochastically.
As discussed in the main text, we are interested in bound states, i.e., eigen states that approach zero at large interparticle distances. To solve the time-independent Schrödinger equation, we expand the relative portion Φ rel of the eigen state Ψ sought in terms of a set of nonorthogonal eigen functions ψ j [65,66], The basis functions ψ j depend on the relative spatial and the spin degrees of freedom, where χ j denotes an N -particle spin function that is chosen from the complete set of 2 N possible spin functions. The spatial parts φ j are written in terms of a total of (N − 1)(N/2 + 3) non-linear variational parameters d Here, the superscript "(j)" serves to remind us that each basis function is characterized by a set of non-linear variational parameters. The non-linear variational parameters d k determine the spatial oscillations due to the k so -dependent one-body terms. We do find that the values of s (j) k can depend quite strongly on the spin basis function considered. This shows that the parameters s (j) k govern, to leading order, the interplay between the spatial and spin degrees of freedom. Of course, strictly speaking, the influence of the singleparticle and two-particle interaction terms in the Hamiltonian on the eigen states cannot be separated; rather, the eigen states are the result of the relative importance of each of these terms and the kinetic energy terms. The N (N − 1)/2 non-linear parameters d are optimized semi-stochastically, i.e., the basis set is constructed so as to minimize the energy of the eigen state under study. The symmetrizerŜ in Eq. (B2) ensures that the basis functions, and hence the eigen state sought, are fully symmetrized. For the two-boson system, e.g., the symmetrizer readsŜ = (1 +P 12 )/ √ 2, whereP 12 exchanges the spatial and spin degrees of freedom of particles 1 and 2. For N identical particles, the symmetrizer contains N ! terms.
The linear parameters c j are determined by solving the generalized eigen value problem that depends on the Hamiltonian matrix and, due to the non-orthogonality of the basis functions, the overlap matrix [65]. The results presented in this paper utilize basis sets consisting of up to N b = 1, 800 basis functions. Key strengths of the numerical approach employed are that the Hamiltonian and overlap matrix elements have compact analytical expressions and that the non-linear variational parameters can be adjusted so as to capture correlations that occur on scales smaller than r 0 and larger than (k so ) −1 .
The energies for the N = 2 system depend on the parameters a s , k so , Ω, andδ. For N = 3, they additionally depend on κ * . In practice, we set q N = 0 and scan Ω/E so , δ/E so , and a s k so (for N = 3, we fix κ * ). To obtain the relative eigen energies E and eigen states Φ rel for finite q N,z , we do not need to redo the numerical calculations. Instead, we use the "conversion" implied by Eq. (21), i.e., the relative eigen energy and eigen states are obtained by changing from δ toδ. The full eigen states Ψ are obtained by multiplying the q N = 0 eigen states by the center-of-mass piece Φ cm .
As in the case without spin-orbit coupling [67], it tends to be more efficient to describe each relative eigen state Φ rel by its own basis set as opposed to constructing one basis set that describes multiple eigen states well. We refer to the eigen state sought-this could be the ground state or one of the excited states-as target state. To construct the basis, we add one basis function at a time. Let us assume that we have a basis set of size N initial and that we want to enlarge the basis set by one basis function. To do this, we generate N trial trial basis functions (N trial is of the order of a few hundred to a few thou-sand) and calculate the energy E k (k = 1, · · · , N trial ) of the target state for each of these enlarged basis sets. Assuming that the generalized eigen value problem for the basis set of size N initial has been solved, the trial energies can be calculated via a root-finding procedure [65], which is, generally, computationally significantly faster than solving the generalized eigen value problem. The basis function to be added is determined by which of the N trial trial energies is lowest, i.e., by looking for the trial basis function that lowers the energy of the target state the most. After the "best" trial function has been added, the generalized eigen value problem is solved for the basis set of size N initial + 1 and the procedure is repeated to generate a basis set of size N initial + 2.
The convergence of the energy with increasing basis set size can be improved significantly by choosing "good" basis functions, i.e., by generating basis functions that efficiently cover the entire Hilbert space. Conversely, if the basis set is not constructed carefully, the energy may not even converge. In our implementation, the parameters d kl are chosen so as to cover length scales ranging from less than r 0 to a few times the length set by the binding energy of the target state. Since the target energy is, in general, not known a priori, the parameter windows are typically refined based on results obtained in preliminary calculations. To choose parameter windows for the x-, y-and z-components of s (j) k , we are guided by the non-interacting two-and three-particle dispersion curves. For example, if the non-interacting two-particle dispersion along the z-coordinate exhibits two minima, we choose s 1,x and s (j) 1,y are selected from windows that include zero. The widths of the parameter windows are adjusted through an educated trial and error procedure. Among other things, we check if the parameters selected by the code are clumped in a particular region of the parameter space.

Appendix C: Three-atom threshold
The three-atom system with center-of-mass momentum q N , N = 3, is bound if its energy is lower than that of three infinitely far separated atoms with the same centerof-mass momentum and, if a two-body bound state exists, lower than that of an infinitely far separated dimer and atom with the same center-of-mass momentum. To determine the lowest three-body scattering threshold, one thus needs to know the lowest dimer binding energy for all two-body center-of-mass momenta. As a consequence, the three-body scattering threshold depends oñ δ/E so and Ω/E so as well as on a s k so . We define the scattering threshold E th usingĤ rel and denote the eigen energies ofĤ rel by E.
We start by determining the lowest relative three-body scattering threshold in the absence of two-body bound states. In this case, the lowest relative scattering threshold is determined by the minimum energy of the noninteracting relative dispersion curves for fixedδ/E so and Ω/E so . The dispersion curves depend on two relative momenta (namely q 1,z and q 2,z ) and the total number of dispersion curves is eight. Figures 4(a)  ferent spin channels. As an example, Figs. 4(d)-4(f) show the lowest relative non-interacting dispersion curves for Ω/E so = 2 andδ/E so = 0, 2.278, and 4, respectively. As for vanishing Raman coupling, the number of global minima changes from six to three (not shown) to four to one with increasingδ. However, the critical generalized detuningδ at which these changes occur differs for Ω/E so = 2 and Ω = 0.
The minimum of the non-interacting relative threeatom dispersion curves defines, assuming two-body bound states are absent, the lowest three-atom scattering threshold. Figure 5 shows the lowest three-atom scattering threshold energy E aaa th as functions ofδ/E so and Ω/E so . The thick open circles and thick dashed line indicate the parameter combinations at which the number of global minima is six and four, respectively. For parameter combinations above the thick open circles and below the thick dashed line the number of global minima of the lowest non-interacting relative dispersion curve is equal to three. Above the thick dashed line the number of global minima is equal to one. If we assume that the three-body binding energy is, approximately, largest when the degeneracy of the lowest noninteracting relative dispersion curve is largest, then Fig. 5 suggests that the three-body system on the negative scattering length side, provided two-body bound states are absent, is enhanced the most compared to the energy of the system without spin-orbit coupling whenδ = 0 or q 3,z = −3mδ/(2 k so ). The main text shows that this reasoning provides an intuitive understanding for the behavior of the lowest three-boson state in each manifold. However, the situation for the excited states in a manifold is more intricate [61]. the z-component q 3,z of the three-body center-of-mass momentum, which is a conserved quantity, and on q 2,z , which is treated as a parameter. Assuming that the distance between the center-of-mass of the dimer and the atom is large compared to the size of the dimer and compared to the ranges r 0 and R 0 of the two-and three-body interactions, the coupling term V coupling , V coupling = [V 2b (r 13 ) + V 2b (r 23 ) + V 3b (r 123 )]I 1 ⊗ I 2 ⊗ I 3 , can be set to zero. Thus, the q 2,z -dependent relative atom-dimer dispersion curves are obtained by adding the eigen energies ofĤ 12 andĤ 3 , which depend parametrically on q 2,z . Equations (D1)-(D6) assumed that the dimer is formed by atoms 1 and 2. Alternatively, the dimer could be formed by atoms 1 and 3 or by atoms 2 and 3. These alternative divisions yield atom-dimer dispersion curves that depend on the z-component of the momentum that is associated with the distance vector between particle 2 and the center-of-mass of the 13-dimer and the zcomponent of the momentum that is associated with the distance vector between particle 1 and the center-of-mass of the 23-dimer, respectively. Since we are considering three identical bosons, the three divisions are equivalent.
In what follows, we use q ad,z to reflect that we could single out any of the three atoms. The corresponding atom-dimer energy is denoted by E ad th . Since there exist up to three two-boson bound states [61], the three-boson system supports up to six atom-dimer dispersion curves (there could be four or two). As an example, Fig. 6 shows the energy E ad th of the lowest relative atom-dimer dispersion curve as a function of q ad,z for (a s k so ) −1 = 0.01128, Ω/E so = 2, and variousδ, i.e.,δ = 0, 2.287, and 3.5. The system supports, for this Raman coupling strength and scattering length, one weakly-bound two-boson state for all twobody center-of-mass momenta. Forδ = 0 (solid line in Fig. 6), the atom-dimer dispersion is symmetric with respect to q ad,z = 0 and supports two global minima at finite q ad,z . The break-up into a dimer and an atom is energetically most favorable when q ad,z /( k so ) is equal to ±0.76. This translates, using Eqs. (D3) and (D5), intoδ 12,eff /E so = ±1.52 andδ 3,eff /E so = ±3.04. For δ > 0 (the dashed and dotted lines in Fig. 6 are for δ/E so = 2.287 and 3.5, respectively), the atom-dimer dispersions are asymmetric with respect to q ad,z = 0 and exhibit a global minimum at negative q ad,z , which approaches q ad,z = 0 in theδ → ∞ limit. Intuitively, the asymmetry can be understood by realizing that the atom and the dimer already see a detuning. Thus, moving in the positive momentum direction is not equivalent to moving in the negative momentum direction. The minimum of the lowest relative atom-dimer dispersion curve decreases with increasingδ. The contours show the energy E th , in units of Eso, of the lowest three-boson scattering threshold as functions of (askso) −1 andδ/Eso. To the left of the thick dotted line, bound dimers are not supported and the lowest threshold is given by E aaa th . In this regime, the lowest scattering threshold is independent of (askso) −1 . To the right of the thick dotted line, a weakly-bound bosonic dimer exists and the lowest threshold is given by E ad th . In this regime, the lowest scattering threshold depends on (askso) −1 . The thick open circles and thick dashed line mark the ((askso) −1 ,δ/Eso) combinations for which the lowest three-atom threshold has a degeneracy of six and four, respectively. The thick dash-dotted line marks the ((askso) −1 ,δ/Eso) combinations for which the lowest atom-dimer threshold has a degeneracy of two. In the region encircled by the thick open circles, the thick dotted line, the thick dashed line, and the left edge of the figure the degeneracy of the three-atom threshold is equal to three. In the region encircled by the thick dashed line, the thick dotted line, the upper edge of the figure, and the left edge of the figure, the degeneracy of the three-atom threshold is equal to one. In the region encircled by the thick dash-dotted line, the right edge of the figure, the upper edge of the figure, and the thick dotted line the degeneracy of the atom-dimer threshold is equal to one. The range of scattering lengths shown on the horizontal axis is, even though different units are used, the same as in Fig. 3.