Born-Oppenheimer study of two-component few-particle systems under one-dimensional confinement
BBorn-Oppenheimer study of two-component few-particle systems underone-dimensional confinement
N. P. Mehta
Department of Physics and Astronomy, Trinity University, San Antonio, TX 78212-7200 ∗ (Dated: September 23, 2018)The energy spectrum, atom-dimer scattering length, and atom-trimer scattering length for systemsof three and four ultracold atoms with δ -function interactions in one dimension are presented asa function of the relative mass ratio of the interacting atoms. The Born-Oppenheimer approachis used to treat three-body (“HHL”) systems of one light and two heavy atoms, as well as four-body (“HHHL”) systems of one light and three heavy atoms. Zero-range interactions of arbitrarystrength are assumed between different atoms, but the heavy atoms are assumed to be noninteractingamong themselves. Fermionic and bosonic heavy atoms with both positive and negative parity areconsidered. I. INTRODUCTION
Cold-atom experiments now have the ability to simul-taneously control the atom-atom scattering length andthe trapping geometry. Quantum gases with essentiallyzero-range interactions in one-dimensional (1D) trap ge-ometries have been realized [1–6]. At the same time,the variety of atomic species that have been trappedand cooled, including all of the alkali metals, contin-ues to grow, ranging in mass from Hydrogen [7] to Ra-dium [8]. Moreover, quantum degenerate mixtures ofatoms have been the subject of several experiments re-lated to, for example, the creation of a gas of degeneratepolar molecules. [9], the observation of heteronuclear Efi-mov states [10], and the realization of mixtures of alkaliatoms with alkaline-earth-like atoms [11].These recent experimental advances were preceded bya large body of literature on the few and many-bodyphysics of strongly interacting 1D systems [12–17]. Morerecent theory work includes the calculation of the 3-bosonhyperradial potential curves [18–20], three-body recom-bination rates and threshold laws [21], and benchmarkquality hyperspherical calculations of three-boson bind-ing energies and scattering amplitudes [22]. The three-body problem for unequal masses has been studied infree-space [23] and in an optical lattice[24].Of particular relevance to the present study is themass-dependent calculation of atom-dimer (2+1) scat-tering lengths and three-body binding energies performedin [23]. The calculations of [23] incorporate all of the adi-abatic hyperspherical potential curves necessary for nu-merical convergence. Here, instead of the (in principle)exact adiabatic hyperspherical representation [25], we usethe Born-Oppenheimer approach. For the three-bodycalculations presented here, the accuracy of the Born-Oppenheimer factorization is studied by direct compari-son to the results of [23], and that comparison gives somequantitative insight to the accuracy of the four-body cal-culations that follow. ∗ [email protected] It should also be noted that the HHHL system for spin-polarized heavy fermions in three-dimensions (3D) hasbeen studied by Castin et al. [26]. They found thatfor heavy fermions with J Π = 1 + symmetry, an infi-nite set of four-body states appears in the mass range13 . < m H /m L < . et al. argue thatthese states have Efimov character, however there seemsto be some debate in the literature. Other authors [27]have argued these are truly new states with propertiesdistinct from Efimov states. The authors of [27] con-sider particles interacting with attractive 1 /r interac-tions, basing their model on a Born-Oppenheimer cal-culation of the potential energy surface governing theheavy-particle dynamics. Better establishing the accu-racy of the Born-Oppenheimer approximation for short-range potentials could potentially play a role in the in-terpretation of these calculations.The Born-Oppenheimer approach has been success-fully applied to cold-atom systems in optical lattices tostudy novel crystalline phases in Fermi mixtures [28].The authors of [28] note that the large mass ratios neededto observe these crystalline phases can be achieved withsmall filling factors by tuning the effective mass for theheavy particles. We note that such a scheme could po-tentially be used to observe the tetramer states predictedin this work.In this paper, we consider 1D systems of three andfour particles in which one particle is “light” (of mass m L = βm H with 0 < β <
1) in comparison to theremaining “heavy” (mass m H ) particles. We restrictour attention to cases of noninteracting heavy particles( a HH → ∞ ). Here, a HH is the 1D heavy-heavy scat-tering length. We denote the 1D heavy-light scatteringlength simply by a . For cylindrical harmonic traps inwhich only the lowest transverse mode is significantlypopulated, the 1D scattering length may be expressedin terms of the 3D s-wave scattering length a D andthe transverse oscillator length a ⊥ by the Olshanii for-mula [29, 30]: a = − a ⊥ a D (cid:18) − C a D a ⊥ (cid:19) , (1) a r X i v : . [ c ond - m a t . qu a n t - g a s ] A p r where C ≈ . a ⊥ = Ca D , Eq. (1) predicts a “confinement inducedresonance” (CIR), and the 1D scattering length vanishes.The degree to which the renormalization of the 1Datom-atom scattering length by Eq. (1) accounts for thequasi-1D nature of the confinement in few-body calcu-lations is not a trivial question [31–33]. Fully quasi-1Dfew-body calculations are complicated by the fact thatcylindrical confinement breaks spherical symmetry, andthe total angular momentum of the three or four-bodysystem is not a good quantum number. In this paper, weproceed under the assumption that meaningful few-bodyobservables may be calculated with purely 1D δ -functioninteractions, renormalized according to Eq. (1).This paper is organized as follows. In Section II, wecalculate the Born-Oppenheimer potential curve describ-ing the effective heavy-heavy interaction as mediated bythe light particle. The HHL bound-state spectrum andthe H-HL scattering length is calculated as a functionof the heavy-light mass ratio. The accuracy of the Born-Oppenheimer approximation is studied by comparison tothe high-accuracy calculation of [23].In Section III, we calculate the two-dimensional po-tential energy surface describing the heavy-particle dy-namics in the HHHL system. The adiabatic wavefunc-tion describing the light particle is governed by a one-dimensional Schr¨odinger equation with three δ -functions.We choose coordinates such that for a given permuta-tion of heavy particles, the ordering of the δ -functionsalong the light-particle coordinate is fixed. The result-ing energy surface is then used in a calculation of the three-body adiabatic hyperradial potential curves for theheavy particles. From those potential curves, the HHHLbinding energies and H-HHL scattering lengths are cal-culated. II. THREE-BODY (HHL) PROBLEM
Let particles 1 and 2 have mass m = m = m H andparticle 3 have mass m L = βm H . Throughout this paper,we set ¯ h = 1. For a zero-range heavy-light interaction ofthe form V ij = gδ ( x i − x j ), the 1D H-L scattering lengthis a = − / ( µ HL g ), and assuming a >
0, the heavy-lightbinding energy is B = µ HL g / / (2 µ HL a ). Forparticle positions { x , x , x } , we introduce the followingunitless mass-scaled Jacobi coordinates (See Fig. 3): x = 1 a (cid:114) µ µ b ( x − x ) (2) y = 1 a (cid:114) µ , µ b (cid:18) m x + m x m + m − x (cid:19) (3)Here, µ = m H / µ , = m H [2 β/ (2 + β )] and µ b = √ µ µ , are reduced masses. The heavy-light reducedmass is µ HL = m H [ β/ (1 + β )]. It is convenient to scalethe Hamiltonian by the heavy-light binding energy: B = 1 m H a β + 12 β , (4)so that all energies are measured in units of B . TheSchr¨odinger equation then reads: − µ (cid:18) ∂ ∂x + ∂ ∂y (cid:19) Ψ( x, y ) + g ( λδ (2 x ) + δ ( y + x ) + δ ( y − x )) Ψ( x, y ) = E Ψ( x, y ) . (5)The parameter λ is the ratio of the heavy-heavy couplingto the heavy-light coupling. In this work, only λ → λ → ∞ are considered. The notational cost of scaling by B is contained in the definition of the following unitlessparameters: µ = 1 + β (cid:112) β (2 + β ) (6) g = − √ (cid:18) β β (cid:19) / (7) x = x (cid:115) β β (8)We now assume the wavefunction may be approximatedby the Born-Oppenheimer product:Ψ( x, y ) = Φ( x ; y ) ψ ( x ) (9) where Φ( x ; y ) is a solution to the fixed- x equation, (cid:20) − µ ∂ ∂y + g ( δ ( y + x ) + δ ( y − x )) (cid:21) Φ( x ; y )= u ( x )Φ( x ; y ) , (10)and u ( x ) < x ; y ) and the potential curve u ( x ) are independent of λ .Inserting Eq. (9) into Eq. (5) and making use of Eq. (10)yields, (cid:32) − µ ∂ ∂x + g λδ (2 x ) + u ( x ) + ˜ Q ( x )2 µ (cid:33) ψ ( x ) = Eψ ( x )(11)where, ˜ Q ( x ) = (cid:28) ∂ Φ ∂x (cid:12)(cid:12)(cid:12)(cid:12) ∂ Φ ∂x (cid:29) y (12)It is understood that the integration in the matrix ele-ment ˜ Q ( x ) is carried out over the y coordinate only, whilethe adiabatic coordinate x is held fixed. A. Solution to The Adiabatic Equation
Equation (10) is symmetric with respect to the oper-ation y → − y , and so the eigenstates Φ( x ; y ) must beeven or odd under that operation. The elementary solu-tions that vanish as | y | → ∞ are conveniently written forpositive y as:Φ( x ; y ) = (cid:40) A sinh( κy ) + B cosh( κy ) , if 0 ≤ y ≤ x ; De − κy , if x ≤ y ; (13)where κ ( x ) = (cid:112) − µ u ( x ). For the even solution, A = 0,while for the odd solution, B = 0. Matching thewavefunctions, and imposing the derivative discontinuityacross the delta-function at y = x leads to the followingtranscendental equation for the eigenvalue κ : κg µ + 1 = ( − P +1 e − κx . (14) Here, P = 0 corresponds to the (even) solution for which ∂ Φ ∂y (cid:12)(cid:12)(cid:12) y =0 = 0, and P = 1 corresponds to the (odd) solutionfor which Φ | y =0 = 0. Borrowing language from molecularphysics, one can view the P = 0 solution as belonging tothe “bonding” orbital, and the P = 1 solution to the“anti-bonding” orbital.The potential curves resulting from the x -dependentsolution to Eq. (14) for β − = 22 .
08 (for Li-Cs mixtures)are shown in Fig. 1(b) and 1(d). The potential curvesshown in these two graphs are identical because Eq. (5)is independent of the heavy-particle symmetry. Any ap-parent differences are due to the energy scales on thegraph. The bound-state structure, however, is dependenton the heavy symmetry through the boundary conditionplaced on ψ ( x ) at x = 0. Note that the λ = ∞ solu-tions to Eq. (11) for heavy bosons are identical to thosefor noninteracting heavy fermions. The boundary con-dition ψ (0) = 0 is applied for fermionic heavy atoms aswell as fermionized bosonic atoms, leading to the corre-spondence first recognized in [34]. For small heavy-atomseparations, there is no P = 1 negative energy solutionto Eq. (14), and the light particle is lost to the contin-uum where the excited state potential terminates at thezero-energy threshold.The atom-dimer scattering length and the HHL spec-trum are to a very good approximation determined solelyby the potential curve corresponding to the P = 0 solu-tion to Eq. (14), so we shall restrict our immediate focusto that solution. When the heavy atoms are far apart, | y − x | (cid:29) | x | , Eq. (5) behaves as though there is asingle δ -function of modified strength at the origin. Thesolutions to Eq. (14) when x → ∞ underestimate thecorrect threshold energy by β/ (2 + β ):lim | x |→∞ u ( x ) = − − β β (15)This result is not unexpected, since we have so far ne-glected the positive-definite contribution ˜ Q ( x ) / (2 µ ) tothe heavy-particle kinetic energy. It is known that ne-glecting this “diagonal correction” — which we call the“Extreme Adiabatic Approximation” (EAA) — yieldsa lower bound E EAA to the N -body bound-state en-ergy. Including the diagonal correction, but neglectingany couplings between Born-Oppenheimer curves — anapproximation we call the “uncoupled adiabatic approx-imation” (UAA) — yields an upper bound E UAA to thecorrect energy [35–38]. We find for this problem that thetrend in these inequalities E EAA < E < E
UAA is alreadypresent in the threshold values of the adiabatic potentialitself. In other words, we find that in the limit | x | → ∞ , u ( x ) < − < u ( x ) + ˜ Q ( x ) / (2 µ ). In the next section, weexplicitly calculate ˜ Q ( x ). B. The Diagonal Correction ˜ Q/ (2 µ ) Using the solutions Eq. (13) (with A = 0) along withthe normalization Eq. (16), we explicitly calculate the in-tegral involved in the nonadiabatic correction Eq. (12).Taking A = 0 in Eq. (13), continuity of the wavefunc-tion immediately yields D = B cosh( κx ) /e − κx . Theremaining normalization constant, B , depends on theH-H separation distance both explicitly, and implicitlythrough the eigenvalue κ : B ( x ) = 2 √ κ √ x κ + e x κ + 1 (16)Derivatives of the eigenvalue κ ( x ) are replaced by expres-sions involving κ itself by differentiating Eq. (14) withrespect to x and solving for κ (cid:48) . We find that the nona-diabatic correction ˜ Q ( x ) can be expressed as a rationalpolynomial in the separation distance x :˜ Q ( x ) = β β ) ( − hx + 2 κx + 1) × (cid:2) h + x (cid:0) − h + 24 h κ − hκ + 24 κ (cid:1) + x (cid:0) − h κ + 48 h κ − hκ + 16 κ (cid:1) + x (cid:0) h − h κ + 108 h κ − hκ + 48 κ (cid:1) + x (cid:0) − h κ + 64 h κ − h κ + 64 hκ − κ (cid:1) (cid:3) , (17) FIG. 1. (color online) Graph (a) shows the bosonic hyperradial potential curves describing the heavy-particle dynamics inthe HHHL system, while graph (b) shows the corresponding Born-Oppenheimer potential curves for the HHL system. Graphs(c) and (d) show similar curves, except for fermions. All graphs are for β − = 22 .
08, appropriate for Li-Cs mixtures. Thedashed-red lines indicate bound states. where we have defined the constant h = µ g . EvaluatingEq. (17) at the asymptotic value of the potential Eq. (15)gives ˜ Q/ µ → [ β/ (2 + β )] + [ β/ (2 + β )] , and including˜ Q/ (2 µ ) in Eq. (11) yields the correct threshold energyto order ( β β ) :lim | x |→∞ (cid:34) u ( x ) + ˜ Q µ (cid:35) = − (cid:18) β β (cid:19) . (18)In other words, for small β , the error in the thresholdenergy vanishes linearly without the diagonal correction,but quadratically when it is included. Interestingly, forthe equal mass case ( β = 1), the UAA gives the cor-rect threshold energy to within 11%. This may seem asomewhat surprising result since the Born-Oppenheimerfactorization is typically expected to fail catastrophicallyin this limit, however other authors [39] have found theBorn-Oppenheimer approach to work surprisingly wellfor short-range s-wave interactions in 3D for a wide vari-ety of mass ratios. It seems that the present 1D calcula-tion shares similar good-fortune. C. Numerical results for the HHL system
Here, we compare the present Born-Oppenheimer cal-culation for the HHL system to the high-accuracy calcu-lations of [23]. Binding energies and scattering solutionsare calculated in the UAA.For the scattering calculation, u ( x ) and ˜ Q ( x ) are cal-culated to 15 digits on a uniform grid, and the Numerovmethod is used to propagate the solution out from x = 0to some x max . The attractive well in u ( x ) widens as themass ratio β − increases. An x max ∼
40 is sufficientfor β − < ∼
10, but must be increased to x max ∼
120 for β − ∼ s , each integrationstep in the Numerov method can introduce an error of or-der s . For N s total steps, an upper bound to the asymp-totic values of the wavefunction of order N s s = x max s is maintained less than 10 − . The asymptotic wavefunc-tion is matched to: ψ ( x ) → (cid:40) cos( kx ) − tan( δ ) sin( kx ) Bosonssin( kx ) + tan( δ ) cos( kx ) Fermions (19)The atom-dimer scattering length is then extracted fromthe effective range expansion as:1 a AD = lim k AD → (cid:40) k AD tan δ Bosons − k AD cot δ Fermions (20)Here, k AD = (cid:112) µ , E rel , while k = √ µ E rel . Themass ratios in Table I (discussed below) are obtained bya bisection root-finding algorithm (either on 1 /a AD or on a AD ) to 6-digit precision. The number of digits reportedhere represents the precision of our calculation. The ac-curacy is best estimated by comparing to the calculationsof [23].Bound-state calculations are performed variationallyby expanding ψ ( x ) in a basis of b-splines, and solving theresulting generalized eigenvalue problem. We have veri-fied that the results are well converged with respect to thenumber of grid points used to interpolate the potential u ( x ), as well as the number and placement of b-splines.In Fig. 2(a), we show the three-body spectrum as afunction of β − / (recall β = m L /m H ). The mass ratiosat which a new state appears, marked by the red crossesfor λ = 0 and red dots for λ = ∞ , trace out a curvegoverned by the β -dependence of the threshold Eq. (18).In the hyperspherical calculation of [23], the threshold isreproduced exactly, and all dots and crosses appear at E /B = − − ( a AD /a ) /π as a function of β − / . Again, the red dots and crosses denote the massratios at which a new state appears and the atom-dimerscattering length a AD → ∞ . The blue stars indicate a AD →
0. In a manner similar to [23], we tabulate theseparticular values of the mass ratio in Table I. Note thatthe Born-Oppenheimer calculation consistently overes-timates the critical mass ratios β − by approximately0 . − .
4. The overestimate is understood, at least qual-itatively, by noting that the
U AA produces an upperbound to the binding energy, and the trend in the spec-trum is for deeper binding as β − increases. The percent-age error in the critical values of β − decreases monoton-ically, as one might expect.The β = 1 HHL ground state for λ = 0 bosons wasfound in [23] to be (in units of B ) E = − . E = − . E ,EAA = − . E ,UAA = − . E ,UAA /E thresh = − . . β − > ∼ .
17. Pricoupenko and Pedri [42] found similarstates in 2D for β − > ∼ .
33. Levinsen and Parish [43]established that these states are continuously connectedas confinement is increased. It is interesting to specu-
FIG. 2. Graph (a) shows the energy spectrum for the three-body (HHL) system as a function β − / . Graph (b) showstan − ( a AD /a ) /π , where a AD is the atom-dimer scatteringlength. For both graphs, the solid-black curves denote nonin-teracting bosonic H-particles. The dashed-blue curves denotenoninteracting fermionic particles, or equivalently, fermion-ized H-particles with λ → ∞ . Red crosses and dots, forbosons and fermions respectively, indicate the appearanceof a new bound state and | a AD | → ∞ . Blue stars indicate a AD → late whether the fermionic state that appears in 1D at β − = 1 ( β − ≈ .
170 in our calculation) is continuouslyconnected to these universal trimer states in higher di-mensions.
III. FOUR-BODY (HHHL) PROBLEM
Let us now turn to the calculation of four-body observ-ables. The basic three-step recipe for this calculation isas follows. First, the Born-Oppenheimer method is usedto calculate the 2D potential energy surface for the heavyparticles in the extreme adiabatic approximation (EAA).Next, this potential energy surface is inserted into a cal-culation of the hyperradial adiabatic potential curves andcouplings. Finally, the resulting set of coupled hyperra-dial equations is solved for the bound-states and atom-trimer scattering length. The entire procedure is thenrepeated for different values of β . If a sufficiently largenumber of hyperradial curves and couplings are includedin the final step, then the accuracy of the calculation islimited almost entirely by the EAA made in the first step. TABLE I. The values of the mass ratio β − = m H /m L forwhich the atom-dimer scattering length is infinite ( a AD → ∞ ,corresponding to the appearance of the n th trimer state), orzero ( a AD → λ →
0) and fermionic H atoms ( λ → ∞ ).Results are compared to Ref. [23]. An asterisk (*) denotes anexact result. λ = 0 λ = 0 n β − | a AD → β − | a AD →∞ This work Ref. [23] This work Ref. [23]1 − − − − .
357 0 .
971 3 .
255 2 . .
747 9 .
365 12 .
336 11 . .
333 22 .
951 26 .
602 26 . .
142 41 .
762 46 .
055 45 . .
168 65 .
791 70 .
695 70 . .
404 95 .
032 100 .
523 100 . .
845 129 .
477 135 .
539 135 . .
488 169 .
120 175 .
742 175 . .
331 213 .
964 221 .
133 220 . λ = ∞ λ = ∞ n β − | a AD → β − | a AD →∞ This work Ref. [23] This work Ref. [23]1 − ∗ .
170 1 ∗ .
499 5 . .
694 7 . .
456 16 . .
373 19 . .
650 32 .
298 36 .
235 35 . .
067 53 .
709 58 .
283 57 . .
697 80 .
339 85 .
518 85 . .
535 112 .
179 117 .
940 117 . .
577 149 .
222 155 .
550 155 . .
820 191 .
463 198 .
347 197 . .
262 238 .
904 246 .
331 245 . A. The adiabatic equations
For all four-body (HHHL) calculations that follow, wechoose particles 1, 2 and 3 to have mass m = m = m = m H and particle 4 to have mass m = βm H .The solution to the adiabatic equation is most easily car-ried out using the “K-type” Jacobi coordinates shown inFig. 3(b), with unitless mass-scaled coordinates definedas: x = 1 a (cid:114) µ µ b ( x − x ) y = 1 a (cid:114) µ , µ b (cid:18) m x + m x m + m − x (cid:19) z = 1 a (cid:114) µ , µ b (cid:18) m x + m x + m x m + m + m − x (cid:19) (21) x y x y z FIG. 3. A schematic diagram of the Jacobi coordinates for(a) the three-body problem and (b) the four-body problemare shown. The heavy particles are contained in the shadedregions. (a) (b)
Here, µ b = ( µ µ , µ , ) / is the four-body reducedmass. Again, we rescale the Schr¨odinger equation bythe heavy-light binding energy B . The full four-bodySchr¨odinger equation then reads: − µ ∇ Ψ( ρ, φ, z ) + (cid:104) g (cid:88) i =1 δ ( z − z i )+ λg (cid:88) i 2, and φ = − φ = π/ 6. Again, rescaling by B introduces the fol-lowing unitless parameters: µ = β + 12 β / (3 + β ) / (23) g = − √ (cid:18) β β (cid:19) / (24) z i = (cid:115) β β ρ sin ( φ − φ i ) (25)where φ = − π/ φ = 0, and φ = − π/ 3. The par-ticular choice of Jacobi coordinates Eq. (21) has the ad-vantage that the separation distances x , x and x areall independent of the z -coordinate. The heavy-particledynamics is restricted to the x - y plane, and the lightparticle can be integrated out by solving an equation inthe z -coordinate only, with fixed x and y . The transfor-mation to hyperspherical coordinates is accomplished byexpressing x , y , and z in terms of the usual spherical po-lar coordinates R , θ and φ . The heavy-particle subsectoris then described by x = ρ cos φ and y = ρ sin φ , where ρ = (cid:112) x + y is the projection of R onto the x - y plane: ρ = R sin θ , and z = R cos θ .Clearly, fixing x and y is equivalent to fixing ρ and φ .We make the Born-Oppenheimer factorization:Ψ = Φ( ρ, φ ; z ) ψ ( ρ, φ ) , (26)where the adiabatic equation for the Born-Oppenheimersurface is: (cid:32) − µ ∂ ∂z + g (cid:88) i =1 δ ( z − z i ) (cid:33) Φ( ρ, φ ; z ) = U ( ρ, φ )Φ( ρ, φ ; z )(27)The heavy-particle eigenstates now live on the potentialenergy surface U ( ρ, φ ), and satisfy (in the EAA): − µ (cid:18) ρ ∂∂ρ ρ ∂∂ρ + 1 ρ ∂ ∂φ (cid:19) ψ ( ρ, φ )+ U ( ρ, φ ) + λg (cid:88) i 6. For arbitrary λ , one would need to accountfor the λ -dependent derivative discontinuity at φ = π/ ρ , which are convenientlywritten in matrix form as: − µ (cid:18) ∂ ∂ρ + Q ( ρ ) + 2 P ( ρ ) ∂∂ρ (cid:19) (cid:126)f ( ρ )+ U eff ( ρ ) (cid:126)f ( ρ ) = E EAA (cid:126)f ( ρ ) (31)Here, U eff is a diagonal matrix with elements U n ( ρ ) − / µ ρ , P mn ( ρ ) = (cid:68) χ m (cid:12)(cid:12)(cid:12) ∂χ n ∂ρ (cid:69) φ and Q mn ( ρ ) = (cid:68) χ m (cid:12)(cid:12)(cid:12) ∂ χ n ∂ρ (cid:69) φ . When P and Q are included in the so-lution to Eq. (31), and enough channels are retained fornumerical convergence, the accuracy of the four-bodyenergy is (in principle) limited only by the omission offirst and second derivative couplings, (cid:68) Φ m (cid:12)(cid:12)(cid:12) (cid:126) ∇ Φ n (cid:69) z and (cid:10) Φ m (cid:12)(cid:12) ∇ Φ n (cid:11) z , that arise from generalizing Eq. (26) toinclude a sum: Ψ = (cid:80) n Φ n ψ n . Such a generalization isnot possible for our model without the introduction ofa confining potential because Eq. (27) admits only onesolution that vanishes as | z | → ∞ . Identical particle symmetry of the heavy particles al-lows one to restrict the domain of the four-body wave-function to the region 0 < φ < π/ 6. Thus, for a givenpermutation of heavy particles, the locus of points de-scribing the coalescence of a heavy particle and a lightparticle — i.e. when z is equal to z i — remain ordered z < z < z along the z -coordinate. Because the order-ing is independent of ρ and φ , the solution to Eq. (27)for all ρ and all φ ∈ [0 , π/ 6] is straightforward.The boundary condition on χ ( ρ ; φ ) at φ = 0 is de-termined by a combination of the parity operator, ˆΠ φ → φ + π , and the 1-2 permutation operator, ˆ P φ → π − φ , bythe rule: ˆ P ˆΠ φ → − φ . Considering positive parity, theboundary conditions on χ ( ρ ; φ ) for noninteracting bosonsare: ∂χ∂φ (cid:12)(cid:12)(cid:12) φ =0 = ∂χ∂φ (cid:12)(cid:12)(cid:12) φ = π/ = 0, while for noninteractingfermions, χ | φ =0 = χ | φ = π/ = 0. Note that the boundaryconditions for noninteracting fermions of positive parityare equivalent to those for bosons of negative parity, but λ → ∞ . B. Numerical Solutions for the HHHL System In Appendix A, we calculate the transcendental equa-tion for the eigenvalue of a 1D Schr¨odinger equationwith three δ -functions of arbitrary strength and arbi-trary placement. The resulting Eq. (A3) is applied tothe Eq. (27) by letting g a = g b = g c = 2 µ g , κ = − µ U ( ρ, φ ) > 0, and a = z , b = z and c = z .Figure. 4 shows the potential energy surface for theparticular mass ratio β − = 22 . 08 appropriate for anatomic mixture of Li-Cs. Potential surfaces like this oneare calculated by solving Eq. (A3) on a nonlinear gridwith typically 200 × 400 points in the ( ρ, φ ) plane. Thepoints are distributed so that more grid-points are con-centrated in the vicinity of the well at ρ = 0 and thevalley near φ = π/ 6. Particular care must be taken todescribe the valley near φ = π/ ρ , orelse the numerical solution to the fixed- ρ adiabatic equa-tion does not reproduce the correct threshold behaviorin any of the atom-trimer channels. This is because thefixed- ρ solutions as ρ → ∞ should approach the HHLbound-state energies from the spectrum in Fig. 2 withthe correct ρ -dependence. In particular, at large ρ wefind that Q ( ρ ) / µ → − / µ ρ , and exactly cancelsthe +1 / µ ρ in the U , eff . That is, the effective po-tential with the diagonal correction approaches a con-stant, and describes a 2-body channel to which Eqs. (19)and (20) may be applied with the replacements ψ → f , k AD → k AT = (cid:112) µ , E rel , a AD → a AT , and x → ρ ,along with the boundary condition f (0) = 0.Figures 1(a) and 1(c) show the hyperradial potentialcurves U n ( ρ ) obtained by solving the fixed- ρ Eq. (30) us-ing the potential energy surface shown in Fig. 4. Notethat at large ρ , the lowest potential curves converge tothe appropriate HHL bound state energy shown as reddashed lines in Fig. 1(b) and 1(d), as appropriate for an Ρ (cid:144) a 0.00.51.0 Φ (cid:144)(cid:72) Π (cid:144) (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) U (cid:72) Ρ , Φ (cid:76) B FIG. 4. The Born-Oppenheimer surface for the HHHL systemis shown. U ( ρ, φ ) is in units of B , and ρ is in units of a . atom-trimer channel. The red dashed lines in Fig. 1(a)and (c) indicate HHHL bound states obtained by solvingEq. (31) with 10 coupled channels. Typically, calcula-tions with only the lowest channel (but including the di-agonal correction) give bound-state energies converged tofour or five digits. The error incurred by ignoring excitedhyperradial potential curves is expected to be negligiblecompared to making the EAA in the calculation of thesurface U ( ρ, φ ).In Figures 5 and 6, we show the spectrum and atom-trimer scattering lengths of the HHHL system with non-interacting bosonic and fermionic heavy atoms, respec-tively. We show both positive (B+, F+), and negative(B-, F-) parity cases for each identical particle symmetry.The HHL ground-state energies for each symmetry fromFig. 2(a) are replotted here as dashed-red curves. Again,mass ratios at which a new tetramer state appears (and | a AT | → ∞ ) are marked by red crosses, while zeroes of a AT are indicated by blue stars. The particular numericalvalues for the coordinates ( β − / , E/B ) are also marked.As the mass-ratio β − increases, four-body bound statesenter at lower energies than one would expect from thethree-body calculation (i.e. the dashed-red curve). Thisdiscrepancy in the threshold energy is attributed to thefact that the EAA underestimates the potential surface inFig. 4 by neglecting the positive nonadiabatic correction − (cid:104) Φ | ∇ Φ (cid:11) z / µ , while the corresponding correction isincluded at the three-body level in Fig. 2.While we perform a multichannel calculation of thespectrum, we find that it is sufficient to use a simplesingle-channel calculation for a AT . Indeed, comparingthe critical mass ratios for which a AT → ∞ , and atetramer state lies at threshold, we find good agreementbetween the two calculations. This can be readily ob-served by comparing the positions of the red crosses ingraphs (a) and (b) of Fig. 5, and similarly in Fig. 6. FIG. 5. Graphs (B+)(a) and (B-)(a) show the HHHL spec-trum for bosonic heavy particles of positive and negative par-ity, respectively. The red crosses indicate the appearance ofa new HHHL bound state and a AT → ∞ . The red dashedcurve is the lowest (solid line) bosonic HHL bound state fromFig. 2(a). Graphs (B+)(b) and (B-)(b) show the arctangentof the atom-trimer scattering length for bosonic heavy par-ticles of positive and negative parity, respectively. The bluestars indicate a AT → 0. The specific ordinates of the crossesand stars are labeled. At β = 1, we find that noninteracting bosons of posi-tive parity admit an HHHL bound state with E , EAA ≈− . 55. The second tetramer state appears at β − ≈ . β − ≈ . 3. For Li-Cs mixtures,one might expect two universal tetramer states. Neg-ative parity bosons are less likely to bind than thosewith positive parity. The first tetramer state appearsat β − ≈ . 4, and the second at β − ≈ . FIG. 6. Graphs (F+)(a) and (F-)(a) show the HHHL spec-trum for fermionic heavy particles of positive and negativeparity, respectively. The red crosses indicate the appearanceof a new HHHL bound state and a AT → ∞ . The red dashedcurve is the lowest (solid line) fermionic HHL bound statefrom Fig. 2(a). Graphs (F+)(b) and (F-)(b) show the arctan-gent of the atom-trimer scattering length for fermionic heavyparticles of positive and negative parity, respectively. Theblue stars indicate a AT → 0. The specific ordinates of thecrosses and stars are labeled. down within the Born-Oppenheimer approximation atthese small mass ratios. The difficulty is magnified be-cause even in the UAA three-body calculations of Sec-tion II, the trimer state doesn’t appear until β − ≈ . β − / ,and the tetramer energy tracks along with it until about β − ≈ . 6. The increasing threshold energy is undoubt-edly an artifact of the approximation at the four-bodylevel since it is absent in the more accurate three-bodycalculations. We can nonetheless estimate the criticalmass ratio as β − ≈ . ± . 6, indicated by a red circle inFig. 6(F-)(a). The second negative parity fermionic stateappears at β − = 19 . β − ≈ . 5. In 2D, Levinsen and Parish [43] found a criticalmass ratio of β − ≈ . 0. It is interesting to speculatewhether these states are continuously connected to eachother, and to the universal state that appears in thesecalculations at β − ≈ . IV. SUMMARY AND OUTLOOK We have calculated three-body and four-body spectra,as well as the atom-dimer and atom-trimer scatteringlengths for two-component systems with one light par-ticle, as a function of the mass ratio. Heavy particlesare assumed to be noninteracting, and the four-particlesystem is assumed to be in free-space. Both bosonic andfermionic heavy particles are treated. For the HHL sys-tem, the Born-Oppenheimer method gives good quan-titative agreement with the hyperspherical calculationsof [23]. For the HHHL system, the potential energy sur-face governing the heavy-particle dynamics is calculatedin the “extreme adiabatic approximation”. That sur-face is then used to calculate hyperradial potential curvesand couplings. The values for the resulting atom-trimerthresholds converge to the appropriate three-body boundstate energies, lending some confidence to the HHHL cal-culation.Let us now discuss possible extensions of this work.Note that we have scaled away the only length scale, a ,that appears in our model. There are two immediategeneralizations that expand the parameter space consid-erably.First, there is the generalization to arbitrary H-H in-teractions, which introduces the H-H scattering length a HH . Such an extension was already treated at thethree-body level in [23], but no such four-body calcu-lations have appeared in the literature. In a hyperspher-ical calculation, the additional derivative discontinuity inthe angular wavefunction is treated analytically, and thehyperradial potential curves are calculated as the solu-tion to a single transcendental equation [21, 23]. TheHHHL Born-Oppenheimer calculation for bosons can beextended to arbitrary λ by choosing a b-spline basis setthat satisfies the boundary condition,lim (cid:15) → χ ∂χ∂φ (cid:12)(cid:12)(cid:12)(cid:12) π/ − (cid:15) = ρλ (1 + β ) √ β (cid:18) β β (cid:19) / (32)With this generalization, one can smoothly transition be-tween the energies shown in Fig. 5 and Fig. 6, passingfrom noninteracting bosons to the fermionized limit.0The bound-state calculation can be extended by theintroduction of a harmonic trapping potential, whichseparates into relative and center-of-mass parts underthe transformation to Jacobi coordinates. This exten-sion would establish a connection with several papersthat have appeared recently, treating equal-mass two-component systems [45–47]. The addition of a trappingpotential would introduce excited potential energy sur-faces and the possibility of interesting physics beyondthe Born-Oppenheimer approximation.Here, we have only considered the “3+1” branch (i.e.the HHHL system) of the few-component problem. Itmay be that other branches can be treated by similarmethods. For example, for the 2+2 (HHLL) problem, in-tegrating out the light atoms would result in a 1D heavy-heavy potential, but the adiabatic equation is a 2D par-tial differential equation, instead of a 1D equation likeEq. (27).Finally, it is worth emphasizing that ultimately a fully3D solution to the few-body problem with finite-range in-teractions is needed in order to understand the physics ofquasi-1D systems. A hyperspherical solution to the few-body problem in quasi-1D for finite-range interactionsremains a significant challenge, although recent advancesin the Correlated Gaussian Hyperspherical method [48]may make these calculations possible. ACKNOWLEDGMENTS I would like to thank Brett Esry, Chris Greene andJose D’Incao for early discussions related to this topic.Thanks also to Jesper Levinsen for a series of helpful correspondences. Appendix A: Triple δ -function problem Here, we give the solution for the eigenvalue to thefollowing Schr¨odinger equation: (cid:20) − ∂ ∂z + g a δ ( z − a ) + g b δ ( z − b ) + g c δ ( z − c ) (cid:21) Φ( z )= − κ Φ( z ) (A1)We assume that the positions of the δ -functions are or-dered as a < b < c , but no other assumptions regardingtheir placement are made. In particular, the Hamiltonianis not assumed to commute with the parity operator. The(unnormalized) solution satisfying asymptotic boundarycondition, Φ( | z | → ∞ ) = 0, is elementary:Φ( z ) = Φ I = Ae κz z < a Φ II = Be − κz + Ce κz a < z < b Φ III = De − κz + Ee κz b < z < c Φ IV = F e − κz c < z (A2)Matching the solutions and enforcing the derivative dis-continuities at z = a , z = b , and z = c yields, afterconsiderable algebra: g a g c ( g b − κ ) e κ ( a + b ) − g a g b ( g c + 2 κ ) e κ ( a + c ) + ( g a + 2 κ )( g b + 2 κ )( g c + 2 κ ) e κ ( b + c ) − g b g c e bκ ( g a + 2 κ ) = 0 (A3)For the special case of a quadrupolar potential ( a = − c , b = 0 and g a = g c = − g b / [1] M. Greiner, I. Bloch, O. Mandel, T. W. H¨ansch, andT. Esslinger, Phys. Rev. Lett. , 160405 (2001).[2] A. G¨orlitz, J. M. Vogels, A. E. Leanhardt, C. Raman,T. L. Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur,S. Gupta, S. Inouye, T. Rosenband, and W. Ketterle,Phys. Rev. Lett. , 130402 (2001).[3] B. Laburthe Tolra and K. M. O’Hara and J. H. Huckansand W. D. Phillips and S. L. Rolston and J. V. Porto,Phys. Rev. Lett. , 190401 (2004).[4] T. Kinoshita, T. Wenger, and D. S. Weiss, Science ,1125 (2004).[5] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature ,900 (2006).[6] A. E. Leanhardt, A. P. Chikkatur, D. Kielpinski, Y. Shin,T. L. Gustavson, W. Ketterle, and D. E. Pritchard,Phys. Rev. Lett. , 040401 (2002).[7] D. G. Fried, T. C. Killian, L. Willmann, D. Landhuis,S. C. Moss, D. Kleppner, and T. J. Greytak, Phys. Rev.Lett. , 3811 (1998). [8] R. H. Parker, M. R. Dietrich, K. Bailey, J. P. Greene,R. J. Holt, M. R. Kalita, W. Korsch, Z.-T. Lu, P. Mueller,T. P. O’Connor, J. Singh, I. A. Sulai, and W. L. Trimble,Phys. Rev. C , 065503 (2012).[9] K.-K. Ni, S. Ospelkaus, M. De Miranda, A. Pe’er,B. Neyenhuis, J. Zirbel, S. Kotochigova, P. Julienne,D. Jin, and J. Ye, Science , 231 (2008).[10] G. Barontini, C. Weber, F. Rabatti, J. Catani, G. Thal-hammer, M. Inguscio, and F. Minardi, Phys. Rev. Lett. , 043201 (2009).[11] H. Hara, Y. Takasu, Y. Yamaoka, J. M. Doyle, andY. Takahashi, Phys. Rev. Lett. , 205304 (2011).[12] E. Lieb and W. Liniger, Phys. Rev. , 1605 (1963).[13] J. B. McGuire, J. Math. Phys. , 622 (1964).[14] C. N. Yang, Phys. Rev. Lett. , 1312 (1967).[15] C. Yang and C. Yang, J. Math. Phys. , 1115 (1969).[16] L. Dodd, J. Math. Phys. , 207 (1970).[17] H. B. Thacker, Phys. Rev. D , 838 (1974). [18] W. G. Gibson, S. Y. Larsen, and J. Popiel, Phys. Rev.A , 4919 (1987).[19] A. Amaya-Tapia and S.Y. Larsen and J. Popiel, Few-Body Syst. , 87 (1997).[20] N. P. Mehta and J. R. Shepard, Phys. Rev. A , 032728(2005).[21] N. P. Mehta, B. D. Esry, and C. H. Greene, Phys. Rev.A , 022711 (2007).[22] O. Chuluunbaatar, A. Gusev, M. Kaschiev, V. Kaschieva,A. Amaya-Tapia, S. Larsen, and S. Vinitsky, J. Phys. B:At. Mol. Opt. Phys. , 243 (2006).[23] O. I. Kartavtsev, A. V. Malykh, and S. A. Sofianos,ZhETF , 419 (2009).[24] G. Orso, E. Burovski, and T. Jolicoeur, Phys. Rev. Lett. , 065301 (2010).[25] J. H. Macek, J. Phys. B: At. Mol. Opt. Phys. , 831(1968).[26] Y. Castin, C. Mora, and L. Pricoupenko, Phys. Rev.Lett. , 223201 (2010).[27] N. L. Guevara, Y. Wang, and B. D. Esry, Phys. Rev.Lett. , 213202 (2012).[28] D. S. Petrov, G. E. Astrakharchik, D. J. Papoular, C. Sa-lomon, and G. V. Shlyapnikov, Phys. Rev. Lett. ,130407 (2007).[29] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[30] T. Bergeman, M. G. Moore, and M. Olshanii, Phys. Rev.Lett. , 163201 (2003).[31] C. Mora, R. Egger, and A. O. Gogolin, Phys. Rev. A , 052705 (2005). [32] C. Mora, A. Komnik, R. Egger, and A. O. Gogolin, Phys.Rev. Lett. , 080403 (2005).[33] J. Levinsen, P. Massignan, and M. Parish, (2014),arXiv:1402.1859 [cond-mat.quant-gas].[34] M. D. Girardeau, J. Math Phys. , 516 (1960).[35] V. F. Brattsev, Sov. Phys. Dokl. (1965).[36] S. T. Epstein, J. Chem. Phys. , 836 (1966).[37] A. F. Starace and G. L. Webster, Phys. Rev. A , 1629(1979).[38] H. T. Coelho and J. E. Hornos, Phys. Rev. A , 6379(1991).[39] S. K. Adhikari, V. Brito, H. Coelho, and T. Das, Il NuovoCimento B Series 11 , 77 (1992).[40] M. Gaudin and B. Derrida, J. de Phys. , 1183 (1975).[41] O. Kartavstev and A. Malykh, J. Phys. B: At. Mol. Opt.Phys. , 1429 (2007).[42] L. Pricoupenko and P. Pedri, Phys. Rev. A , 033625(2010).[43] J. Levinsen and M. M. Parish, Phys. Rev. Lett. ,055304 (2013).[44] D. Blume, Phys. Rev. Lett. , 230404 (2012).[45] S. E. Gharashi and D. Blume, Phys. Rev. Lett. ,045302 (2013).[46] T. Sowi´nski, T. Grass, O. Dutta, and M. Lewenstein,Phys. Rev. A , 033607 (2013).[47] E. J. Lindgren, J. Rotureau, C. Forssn, A. G. Volosniev,and N. T. Zinner, (2013), arXiv:1304.2992v1 [cond-mat].[48] K. M. Daily and C. H. Greene, Phys. Rev. A , 012503(2014).[49] S. Patil, Eur. J. Phys.30