Analytical solution for the spectrum of two ultracold atoms in a completely anisotropic confinement
AAnalytical solution for the spectrum of two ultracold atoms in a completelyanisotropic confinement
Yue Chen, Da-Wu Xiao, Ren Zhang, and Peng Zhang
1, 2, ∗ Department of Physics, Renmin University of China, Beijing, 100872, China Beijing Computational Science Research Center, Beijing, 100193, China School of Science, Xi’an Jiaotong University, Xi’an, 710049, China
We study the system of two ultracold atoms in a three-dimensional (3D) or two-dimensional(2D) completely anisotropic harmonic trap. We derive the algebraic equation J ( E ) = 1 /a ( J ( E ) = ln a ) for the eigen-energy E of this system in the 3D (2D) case, with a and a beingthe corresponding s -wave scattering lengths, and provide the analytical expressions of the functions J ( E ) and J ( E ). In previous researches this type of equation was obtained for spherically oraxially symmetric harmonic traps (T. Busch, et. al. , Found. Phys. , 549 (1998); Z. Idziaszek andT. Calarco, Phys. Rev. A , 022712 (2006)). However, for our cases with a completely anisotropictrap, only the equation for the ground-state energy of some cases has been derived (J. Liang andC. Zhang, Phys. Scr. , 025302 (2008)). Our results in this work are applicable for arbitraryeigen-energy of this system, and can be used for the studies of dynamics and thermal-dynamicsof interacting ultracold atoms in this trap, e.g. , the calculation of the 2nd virial coefficient or theevolution of two-body wave functions. In addition, our approach for the derivation of the aboveequations can also be used for other two-body problems of ultracold atoms. I. INTRODUCTION
The two-body problems of trapped interacting ultra-cold atoms are basic problems in cold atom physics [1–13]. They are of broad interest because of the followingreasons. First, the two-body systems are “minimum” in-teracting systems of trapped ultracold atoms, and onecan obtain a primary understanding for the interactionphysics of an ultracold gas from the analysis of such sys-tems [14]. Second, the solutions to these problems can bedirectly used to calculate some important few- or many-body quantities [15–26], e.g. , the 2nd virial coefficientwhich determines the high-temperature properties of theultracold gases. Third, these systems have been alreadyrealized in many experiments[14, 27–39] where the trapof the two atoms can be created via an optical lattice site[27, 28, 31], an optical tweezer [32–36], or nano-structure[37–39]. In these experiments, by measuring or control-ling the energy spectrum or dynamics of these two atoms,one can, e.g. , create a cold molecule in an optical tweezer[33, 34], derive the parameters of inter-atomic interactionpotential [27, 28, 31, 35] or study various dynamical ef-fects such as the interaction-induced density oscillation[32]. Theoretical results for the corresponding two-bodyproblems are very important for these experimental stud-ies.The most fundamental systems of two trapped ultra-cold atoms are the ones with a harmonic trap whichhas the same frequencies for each atom, and a s -waveshort-range inter-atomic interaction. For these systems,the center-of-mass (CoM) and relative motion of the twoatoms can be decoupled with each other, and the CoMmotion is just the same as a harmonic oscillator. Nev- ∗ [email protected] ertheless, the dynamics of the relative motion of thesetwo atoms is nontrivial. For the three-dimensional (3D)systems, T. Busch, et. al. [1] and Z. Idziaszek and T.Calarco [2, 3] studied the cases with a spherically and ax-ially symmetric harmonic trap, respectively. They showthat the eigen-energy E of the relative motion satisfiesan algebraic equation with the form J ( E ) = 1 a , (1)with a being the 3D s -wave scattering length, and pro-vided the analytical expressions of the function J ( E )for these two cases [1–3], as shown in Table I, respectively.Moreover, for the two-dimensional (2D) systems with anisotropic harmonic trap, T. Busch, et. al. obtained theequation J ( E ) = ln a , (2)for the relative-motion eigen-energy E , where a is the2D s -wave scattering length, and derive expression of thefunction J ( E ) [1]. With these results one can easilyobtain the complete energy spectrum of these two-atomsystems, as well as the corresponding eigen-states. Thus,these results have been widely used in the researches ofcold atom physics [5, 14, 27, 28, 31, 35, 40–44].However, for two atoms trapped in a 3D completelyanisotropic harmonic confinement, the algebraic equationfor arbitrary eigen-energy of relative-motion has not beenobtained so far. Only the equation of the ground-state energy of a system with positive scattering length hasbeen obtained (i.e., the eigen-energy E which satisfies E < E , with E being the relative ground-state energyfor the non-interacting case) in Ref. [4], by J. Liangand C. Zhang in 2008. This equation has the form ofEq. (1), and the corresponding function J ( E ) is also a r X i v : . [ phy s i c s . a t m - c l u s ] F e b type E J ( E ) Ref. ω x = ω y = ω z arbitrary E √ Γ ( − E + ) Γ ( − E + ) Ref. [1] ω x = ω y (cid:54) = ω z arbitrary E − η x √ ∞ (cid:80) n =0 (cid:20) Γ (cid:16) − E − E + nη x (cid:17) Γ (cid:16) − E − E + nη x (cid:17) − √ η x √ n +1 (cid:21) − (cid:112) η x ζ (cid:0) (cid:1) Ref. [3] ( ∗ ) ω x (cid:54) = ω y (cid:54) = ω z E < E − √ π (cid:82) ∞ dt (cid:26) √ η x η y exp[( E − E ) t/ √ − e − ηxt √ − e − ηyt √ − e − t − t / (cid:27) Ref. [4] ω x (cid:54) = ω y (cid:54) = ω z arbitrary E Eq. (39) or Eq. (40) this workTABLE I. Expressions of the function J ( E ) of Eq. (1) for various 3D traps. Here ω α ( α = x, y, z ) is the trapping frequency inthe α -direction, and we use the natural unit (cid:126) = 2 µ = ω z = 1 with µ being the reduced mass of the two atoms. The parameters η x and η y are defined as η x = ω x /ω z and η y = ω y /ω z , respectively, E = ( η x + η y + 1) / z ) and ζ ( z ) are the Gamma and Riemann zeta function, respectively. Thetwo expressions (39) and (40) of J ( E ) given by this work are mathematically equivalent with each other. ( ∗ ): Notice thatthere are typos in Eqs. (21, 23) of Ref. [3]. shown in Table I. On the other hand, since this kind ofconfinements are easy to be experimentally prepared, andare thus used in many experiments of ultracold gases [27–31, 35, 44], it would be helpful if we can derive a generalequation satisfied by all the eigen-energies.In this work we derive the equation for arbitrary eigen-energy of two atoms in a completely anisotropic harmonictrap. We show that this equation also has the form ofEq. (1) and Eq. (2) for the 3D and 2D cases, respec-tively, and provide the corresponding analytical expres-sions of the functions J ( E ) (Table I) and J ( E ) (Eq.(60) or Eq. (62)). As the aforementioned results [1–3] for the spherically and axially symmetric traps, ourequations are useful for the studies of various thermal-dynamcial or dynamical properties of ultracold gases inthe completely anisotropic harmonic confinements. Fur-thermore, our calculation approach used in this work canalso be generalized to other two-body problems of ultra-cold atoms in complicated confinements.The remainder of this paper is organized as follows. InSec. II and Sec. III, we derive the equations for eigen-energies of atoms in 3D and 2D completely anisotropicharmonic traps, respectively. In Sec. IV we discuss howto generalize our approach to other problems. A briefsummary and some discussions are given in Sec. V. Somedetails of our calculations are shown in the appendix. II. 3D SYSTEMS
We consider two ultracold atoms 1 and 2 in a 3D com-pletely anisotropic harmonic trap which has the same fre-quencies for each atom. Here we denote ω α ( α = x, y, z )as the trapping frequency in the α -direction, which sat-isfy ω x (cid:54) = ω y (cid:54) = ω z . (3) For convenience, in this work we use the natural unit (cid:126) = 2 µ = ω z = 1 , (4)with µ being the reduced mass of the two atoms. Wefurther introduce the aspect ratios η x = ω x ω z ; η y = ω y ω z ; η z = ω z ω z = 1 , (5)where η x and η y describe the anisotropy of the trappingpotential.As mentioned in Sec. I, for this system we can separateout the CoM degree of freedom of these two atoms andfocus on the inter-atomic relative motion. The Hamilto-nian operator of our problem is given byˆ H = ˆ H + ˆ V I . (6)Here ˆ H is the free Hamiltonian for the relative motionand can be expressed asˆ H = ˆ p + 14 (cid:2) η x ˆ x + η y ˆ y + ˆ z (cid:3) , (7)with ˆ p and ˆ r ≡ (ˆ x, ˆ y, ˆ z ) being the relative momentumand coordinate operators, respectively. In Eq. (6) ˆ V I isthe inter-atomic interaction operator, which is modeledas the s -wave Huang-Yang pseudo potential. Explicitly,for any state | ψ (cid:105) of the realtive motion we have (cid:104) r | ˆ V I | ψ (cid:105) = 4 πa δ ( r ) ∂∂r [ r · (cid:104) r | ψ (cid:105) ] , (8)with | r (cid:105) being the eigen-state of the relative-coordinateoperator ˆ r with eigen-value r , r = | r | , and a being the3D s -wave scattering length.For our system the parity with respect to the spatialinversion r → − r is conserved, and the contact pseudopotential ˆ V I only operates on the states with even par-ity. Therefore, in this work we only consider the eigen-energies and eigen-states of ˆ H in the even-parity sub-space.Now we deduce the algebraic equation for the eigen-energy E of the total Hamiltonian ˆ H . We begin fromthe schr¨oedinger equation (cid:104) ˆ H + ˆ V I (cid:105) | Ψ (cid:105) = E | Ψ (cid:105) , (9)satisfied by E and the corresponding eigen-state | Ψ (cid:105) ofˆ H . This equation can be re-expressed as | Ψ (cid:105) = 1 E − ˆ H ˆ V I | Ψ (cid:105) . (10)Using Eq. (8) we find that Eq. (10) yields (cid:104) r | Ψ (cid:105) = 4 πa G ( E, r ) (cid:20) ∂∂r [ r · (cid:104) r | Ψ (cid:105) ] (cid:12)(cid:12)(cid:12)(cid:12) r = (cid:21) , (11)where G ( E, r ) is the Green’s function of the free Hamil-tonian ˆ H , which is defined as G ( E, r ) = (cid:104) r | E − ˆ H | (cid:105) . (12)We can derive the the equation for the eigen-energy E by doing the the operation ∂∂r ( r · ) (cid:12)(cid:12) r =0 on both sides ofEq. (11). In this operation, without loss of generality,we choose r = z e z and e z being is the unit vector alongthe z -direction. Then we find that E satisfies J ( E ) = 1 a , (13)which is just Eq. (1) of Sec. I, with the function J ( E )being defined as J ( E ) = 4 π (cid:26) ∂∂ | z | [ | z | · G ( E, z e z )] (cid:12)(cid:12)(cid:12)(cid:12) z → (cid:27) . (14) A. Expression of J ( E ) Next we derive the expression of the function J ( E ).For our system, the ground-sate energy of the relativemotion of two non-interacting atoms is E = 12 ( η x + η y + 1) . (15)In the following we first consider the case with E < E ,which was also studied by Ref. [4], and then investigatethe general case with arbitrary E .
1. Special Case:
E < E When
E < E , the Green’s function G ( E, z e z ) ( z >
0) can be expressed as the Laplace transform of theimaginary-time propagator [45–49],i.e., G ( E, z e z ) = − (cid:90) + ∞ K ( z, E, β ) dβ, (16)with the function K ( z, E, β ) being defined as K ( z, E, β )= e βE (cid:104) z e z | e − β ˆ H | (cid:105) = exp (cid:20) βE − z β (cid:21) (cid:89) α = x,y,z (cid:114) η α π sinh ( η α β ) . (17)Here we emphasize that, when E < E the function K ( z, E, β ) exponentially decays to zero in the limit β →∞ , and thus the integration in Eq. (16) converges for anyfixed non-zero z . Nevertheless, this integration divergesin the limit z →
0. That is due to the behavior of theleading term e − z / (4 β ) / (4 πβ ) of the function K ( z, E, β )in the limit β → + [45–49].We can separate this diver-gence by re-expressing the integration as G ( E, z e z ) = − (cid:90) + ∞ dβ e − z β (4 πβ ) − (cid:90) + ∞ dβ ˜ K ( r, E, β )= − π | z | − (cid:90) + ∞ dβ ˜ K ( z, E, β ) , (18)where ˜ K ( z, E, β ) = K ( z, E, β ) − e − z β (4 πβ ) . (19)In Eq. (18) the integration (cid:82) + ∞ dβ ˜ K ( z, E, β ) uniformlyconverges in the limit z →
0. Using this result, we obtainthe expansion of G ( E, z e z ) in this limit:lim z → G ( E, z e z ) = − π | z | − (cid:90) + ∞ dβ ˜ K (0; E, β ) + O ( z ) . (20)Substituting this result into Eq. (14) and using Eqs. (17,19), we finally obtain the expression J ( E )= − (cid:90) + ∞ dβ (cid:40) e βE √ π (cid:89) α = x,y,z (cid:114) η α sinh ( η α β ) − √ π β (cid:41) (for E < E ) , (21)where the the integration converges for E < E , as theone in Eq. (16). This result was also derived by J. Liangand C. Zhang in Ref. [4].
2. General Case: Arbitrary E In the general case with arbitrary energy E we cannotdirectly use the above result in Eq. (21), because the inte-gration in this equation diverges for E > E . For the sys-tems with spherically or axially symmetric confinements,the authors of Refs. [1, 2] successfully found the analyt-ical continuation of this integration for all real E . How-ever, for the current system with completely anisotropictraps, to our knowledge, so far such analytical continua-tion has not been found.Now we introduce our approach to solve this problem.For convenience, we first define the eigen-energy of thefree Hamiltonian ˆ H , which is just a Hamiltonian of a 3Dharmonic oscillator, as E n . Here n = ( n x , n y , n z ) , (22)with n α = 0 , , , ... ( α = x, y, z ) being the quantum number of the α -direction. It is clear that we have E n ≡ (cid:15) n x + (cid:15) n y + (cid:15) n z , (23)with (cid:15) n α = (cid:18)
12 + n α (cid:19) η α , ( α = x, y, z ) . (24)In addition, we further denote the eigen-state of ˆ H cor-responding to E n as | n (cid:105) .Similar to above, the key step of our approach isto calculate free Green’s function G ( E, z e z ) defined inEq. (12). Since the result in Eq. (16) cannot be used forour general case because the integration in this equationdiverges for E > E , we need to find another expres-sion for G ( E, z e z ), which converges for any E . To thisend, we separate all the eigen-states {| n (cid:105)} of ˆ H into twogroups, i.e., the ones with n ∈ L E and n ∈ U E , respec-tively, with the sets L E and U E being defined as L E : (cid:26) ( n x , n y , n z ) | n x,y,z = 0 , , ..., (cid:15) n x + (cid:15) n y + 12 ≤ E (cid:27) , (25)and U E : (cid:26) ( n x , n y , n z ) | n x,y,z = 0 , , ..., (cid:15) n x + (cid:15) n y + 12 > E (cid:27) , (26)respectively. It is clear that we have E n > E, for all states with n ∈ U E . (27)Furthermore, we can re-express the free Green’s operator 1 / [ E − ˆ H ] as1 E − ˆ H = (cid:88) n ∈ U E | n (cid:105)(cid:104) n | E − E n + (cid:88) n ∈ L E | n (cid:105)(cid:104) n | E − E n (28)= − (cid:90) + ∞ dβ (cid:34) (cid:88) n ∈ U E | n (cid:105)(cid:104) n | e − β ( E n − E ) (cid:35) + (cid:88) n ∈ L E | n (cid:105)(cid:104) n | E − E n . (29)Due to the fact (27), the integration in of Eq. (29) converges . Moreover, the integrand in Eq. (29) can be re-writtenas (cid:88) n ∈ U E | n (cid:105)(cid:104) n | e − β ( E n − E ) = (cid:88) n | n (cid:105)(cid:104) n | e − β ( E n − E ) − (cid:88) n ∈ L E | n (cid:105)(cid:104) n | e − β ( E n − E ) = e βE e − β ˆ H (3D)0 − (cid:88) n ∈ L E | n (cid:105)(cid:104) n | e − β ( E n − E ) , (30)where we have used e − β ˆ H (3D)0 = (cid:80) n | n (cid:105)(cid:104) n | e − βE n . Substituting Eq. (30) into Eq. (29) and then into Eq. (12), weobtain G ( E, z e z ) = − (cid:90) + ∞ dβ [ K ( z ; E, β ) − F ( z ; E, β )] + Q ( z ; E ) , (31)where K ( z ; E, β ) is defined in Eq. (17), and the functions F ( z ; E, β ) and Q ( z ; E, β ) are defined as F ( z ; E, β ) = (cid:88) n ∈ L E (cid:104) z e z | n (cid:105)(cid:104) n | (cid:105) e − β ( E n − E ) ; (32) Q ( z ; E, β ) = (cid:88) n ∈ L E (cid:104) z e z | n (cid:105)(cid:104) n | (cid:105) E − E n . (33)Due to the convergence of the integration in Eq. (29), the integration in Eq. (31) also converges for non-zero z , nomatter if E > E or E < E . Therefore, Eq. (31) is the convergent expression of G ( E, z e z ) for arbitrary E .Furthermore, similar to in Sec. II. A, the integration in Eq. (31) diverges in the limit | z | →
0, due to the leadingterm e − z β / (4 πβ ) of the integrand in the limit β → + . We can separate this divergence via the technique used inEqs. (18-20), i.e., re-express this integration as − (cid:90) + ∞ dβ [ K ( z ; E, β ) − F ( z ; E, β )] = − (cid:90) + ∞ dβ e − z β (4 πβ ) − (cid:90) + ∞ dβ [ ˜ K ( z ; E, β ) − F ( z ; E, β )]= − π | z | − (cid:90) + ∞ dβ [ ˜ K ( z ; E, β ) − F ( z ; E, β )] , (34)with the function ˜ K ( z ; E, β ) being defined in Eq. (19). Using this technique we obtainlim | z |→ G (0) ( E, z e z ) = − π | z | + (cid:20) W ( E ) + (cid:90) + ∞ I ( E, β ) dβ (cid:21) + O ( z ) , (35)where the z -independent functions W ( E ) and I ( E, β ) are defined as W ( E ) ≡ Q (0 , E, β ) and I ( E, β ) ≡− ˜ K (0; E, β )+ F (0; E, β ), respectively. The straightforward calculations (Appendix A) show that they can be expressedas W ( E ) = − π (cid:114) η x η y (cid:88) ( n x ,n y ) ∈ C (3D) E n x + n y − Γ (cid:20) − ( E − (cid:15) nx − (cid:15) ny ) (cid:21) Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1)Γ (cid:20) − ( E − (cid:15) nx − (cid:15) ny ) (cid:21) , (36)and I ( E, β ) = − e βE (cid:89) α = x,y,z (cid:114) η α π sinh ( η α β ) + (cid:18) πβ (cid:19) + (cid:114) πη x η y β (cid:88) ( n x ,n y ) ∈ C (3D) E n x + n y − Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1) e β ( E − (cid:15) nx − (cid:15) ny ) , (37)where Γ( z ) is the Gamma function and C (3D) E is a set of two -dimensional number array ( n x , n y ), which is defined as C (3D) E : (cid:26) ( n x , n y ) | n x,y = 0 , , , , ... ; (cid:15) n x + (cid:15) n y + 12 ≤ E (cid:27) . (38)Finally, substituting Eq. (35) into Eq. (14), we obtain the expression of J ( E ). J ( E ) = 4 π (cid:20) W ( E ) + (cid:90) + ∞ I ( E, β ) dβ (cid:21) , (39)with W ( E ) and I ( E, β ) being given by Eqs. (36) and (37), respectively. Notice that the summations in Eqs. (36,37) are done for finite terms, and the integration in Eq. (39) converges. Therefore, using Eqs. (39, 36, 37) one canefficiently calculate J ( E ).It is clear that Eq. (13) for the eigen-energy and the expression (39) of the function J ( E ) are correct for not onlythe systems in completely anisotropic traps but also the ones in spherically or axially symmetric traps. For the lattertwo cases Eq. (39) is equivalent to the results derived by Ref. [1, 2], which are shown in Table I. B. Techniques for Fast Calculation of J ( E ) There are some techniques which may speed up the calculation for the function J ( E ). First, as shown in AppendixB, J ( E ) can be re-expressed as J ( E ) = 4 π (cid:34)(cid:90) Λ0 A ( E, β ) dβ + B (1)3D ( E, Λ) + B (2)3D ( E, Λ) + (cid:18) π (cid:19) / √ (cid:35) , (40)where Λ is an arbitrary finite positive number, and the functions A ( E, β ) and B (1 , ( E, Λ) are defined as A ( E, β ) = − e βE (cid:89) α = x,y,z (cid:114) η α π sinh ( η α β ) + (cid:18) πβ (cid:19) ; (41) B (1)3D ( E, Λ) = ( − × (cid:88) ( n x ,n y ) ∈ C (3D) E n x + n y − √ πη x η y Γ (cid:16) − E − (cid:15) nx − (cid:15) ny (cid:17) e ( E − (cid:15) nx − (cid:15) ny − )Λ Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1)Γ (cid:16) − E − (cid:15) nx − (cid:15) ny (cid:17) × (cid:112) e − × F (cid:20) , − E − (cid:15) n x − (cid:15) n y , − E − (cid:15) n x − (cid:15) n y , e − (cid:21) ; (42)and B (2)3D ( E, Λ) = (cid:112) πη x η y csch Λ × (cid:88) ( n x ,n y ) / ∈ C (3D) E n x + n y − e ( E − (cid:15) nx − (cid:15) ny − ( e − × F (cid:104) , − E − (cid:15) nx − (cid:15) ny , − E − (cid:15) nx − (cid:15) ny , e − (cid:105) Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1) (cid:2) (cid:0) E − (cid:15) n x − (cid:15) n y (cid:1) − (cid:3) , (43)respectively, with F being the Hypergeomtric function. In Appendix B we prove that (40) is exactly equivalent toEqs. (39) for any finite positive Λ. Thus, in the numerical calculation for a specific problem one can choose the valueof Λ by convenience.Due to the following two facts, the calculation of J ( E ) based Eq. (40) may be faster than the one basedEq. (39): (A) An important difference between the expressions(39) and (40) is that, the integrand I ( E, β ) of Eq. (39)includes a summation (cid:80) ( n x ,n y ) ∈ C (3D) E while the integrand A ( E, β ) of Eq. (40) does not include any summation.On the other hand, in the numerical calculations for theseintegrations for a fixed E , one needs to calculate thevalues of the integrands for many points in the β -axis.Thus, to calculate the term (cid:82) + ∞ I ( E, β ) dβ of Eq. (39)one need to do the the summation (cid:80) ( n x ,n y ) ∈ C (3D) E manytimes, while to calculate the term (cid:82) Λ0 A ( E, β ) dβ of Eq.(40) one does not require to do that. This advantage ismore significant for the large- E cases where this summa-tion includes a lot of terms. (B) Furthermore, the terms in the summation (cid:80) ( n x ,n y ) / ∈ C (3D) E of Eq. (43) decays to zero in the lim-its n x,y → ∞ and the decaying is faster than exponential (Appendix B). Therefore, although this summation in-cludes infinite terms, it can converges fast.Another technique which may be helpful for the fastcalculation of J ( E ) is based on direct conclusions ofEq. (39) or Eq. (40), i.e., for two energies E and E we have J ( E ) − J ( E )= 4 π (cid:20) δW + (cid:90) + ∞ δI ( β ) dβ (cid:21) (44)= 4 π (cid:34)(cid:90) Λ0 δA ( β ) dβ + δB (1)3D (Λ) + δB (2)3D (Λ) (cid:35) (45)where δW = W ( E ) − W ( E ), δI ( β ) = I ( E , β ) − I ( E , β ), δA ( β ) = A ( E , β ) − A ( E , β ) and δB (1 , (Λ) = B (1 , ( E , Λ) − B (1 , ( E , Λ). Therefore, if the value of J ( E ) isalready derived, one can calculate J ( E ) by eitherEq. (39) (Eq. (40)) or Eq. (44) (Eq. (45)), and it ispossible that the summations and integration in Eq. (44)(Eq. (45)) can converge faster. C. Energy Spectrum
By solving Eq. (13) with the expression (39) of J ( E ),we can derive the complete eigen-energy spectrum of thecomplete Hamiltonian ˆ H of the two-atom relative mo-tion. In Fig. 1(a-d) we show the energy spectrum of therelative motion of two atoms in 3D harmonic traps withvarious aspect ratios η x,y , which is derived via Eqs. (13,39). In Fig. 1(a) and (b) we consider the traps withspherical symmetry ( η x = η y = 1) and axial symmetry - 1 0 - 5 0 5 1 0- 5051 0 (cid:1) x = (cid:1) y = 1 E a T h i s w o r k
R e f . [ 1 ] ( a ) - 1 5 - 1 0 - 5 0 5 1 0 1 5- 202468 (cid:1) x = (cid:1) y = 1 . 1 ( b ) a T h i s w o r k
R e f . [ 2 ] E - 1 5 - 1 0 - 5 0 5 1 0 1 5- 20246 (cid:1) x = 1 . 2 , (cid:1) y = 1 . 1 a ( c ) E - 2 0 - 1 0 0 1 0 2 0- 5051 01 52 0 (cid:1) x = 5 , (cid:1) y = 3 ( d ) E a FIG. 1. (color online) The energy spectrum of the relative motion of two atoms in a 3D harmonic trap with aspect rations( η x = η y = 1) (a), ( η x = η y = 1 .
1) (b), ( η x = 1 . , η y = 1 .
1) (c) and ( η x = 5 , η y = 3) (d). We show the results given by Eq. (1)and the expression (39) of function J ( E ) with solid lines. In (a) and (b) we also show the results given by Refs. [1] and [2]with blue dots, respectively. Here we use the natural unit (cid:126) = 2 µ = ω z = 1. ( η x = η y = 1 . η x = 1 . , η y = 1 .
1) and quite different ( η x = 5 , η y = 3)with each other, respectively.As an example of the application of our method, wemake a simple investigation for the energy-level distri-bution for the systems with integer aspect ratios η x and η y . This property is crucial for many important problemssuch as thermalization and quantum chaotic behaviors.It is clear that for non-interacting cases (i.e., for a = 0)the energy levels of these systems are distributed withequal spacing. However, when a (cid:54) = 0 the energy-leveldistributions become uneven, i.e., the energy spacings be-come unequal with each other. Nevertheless, for systemswith different aspect ratios, the “degree of unevenness”of the energy-level distributions, or the fluctuation of theenergy spacings, are quite different. In Ref. [2] it wasfound that for the pancake-shape traps with η x (cid:29) η y = 1, the energies are almost distributed withequal spacing (i.e., the level spacings are almost equal),while for the cigar-shape traps with η x = η y (cid:29) a = ∞ , η x = 5and various values of η y ( η y = 1 , , , , l (cid:1) x = 5 , (cid:1) y = 1 (cid:1) x = 5 , (cid:1) y = 2 (cid:1) x = 5 , (cid:1) y = 3 (cid:1) x = 5 , (cid:1) y = 4 (cid:1) x = 5 , (cid:1) y = 5 (cid:1) l l FIG. 2. (color online) The normalized energy spacing δ l de-fined in Eq. (46) for the lowest 20 excited states of the systemswith a = ∞ , η x = 5 and η y ( η y = 1 , , , , (cid:126) = 2 µ = ω z = 1 ( l = 1 , , , ... ; ∆ l >
0) as the spacing between the l -thand ( l + 1)-th excited-state energy of our systems, andcalculate the “normalized” energy spacing δ l ≡ ∆ l (cid:80) l =1 ∆ l , ( l = 1 , ...,
20) (46)for the lowest 20 excited states for each system. As shownin Fig. 2, when the system is crossed from the pancake-shape limit to the cigar-shape, i.e., η y is increased from1 to 5, the fluctuation of δ l monotonically increases, i.e.,the energy-level distribution becomes more and more un-even. In future works we will perform more systematicinvestigations for the energy-level distribution of two ul-tracold atoms in a completely anisotropic harmonic trap. D. Eigen-State of ˆ H Furthermore, Eqs. (11) and (12) yields that for a giveneigen-energy E of the total Hamiltonian ˆ H for the two-atom relative motion, the corresponding the eigen-state | Ψ (cid:105) satisfies (cid:104) r | Ψ (cid:105) ∝ (cid:104) r | E − ˆ H | (cid:105) = (cid:88) n (cid:104) r | n (cid:105)(cid:104) n | (cid:105) E − E n . (47)Thus, considering the normalization condition (cid:104) Ψ | Ψ (cid:105) =1, we obtain the expression for wave function of the eigen-state: (cid:104) r | Ψ (cid:105) = 1 (cid:113)(cid:80) n | c n | (cid:88) n c n φ n x ( η x , x ) φ n y ( η y , y ) φ n z ( η z , z ) , (48)with c n = 1( E − E n ) φ n x ( η x , φ n y ( η y , φ n z ( η z , . (49)Here ( x, y, z ) are the components of r , the function φ n ( η, X ) is defined in Eq. (A3). It is just the eigenwave function of a one-dimensional harmonic oscillatorand satisfies (cid:104) r | n (cid:105) = φ n x ( η x , x ) φ n y ( η y , y ) φ n z ( η z , z ). Us-ing Eqs. (66, 67) one can easily derive the normalizedeigen-state | Ψ (cid:105) of ˆ H from the eigen-energy E [1, 3, 32]. III. 2D SYSTEMS
Now we consider two atoms confined in a 2Danisotropic harmonic trap in the x − z plane. As above,here we assume the trapping frequencies ω x ( ω z ) in the x ( z )-direction are same for each atom. Thus, by separat-ing out the CoM degree of freedom, we can obtain thefree Hamiltoian operator for the relative motion. In ournatural unit with (cid:126) = 2 µ = ω z = 1 this free Hamiltoniancan be expressed asˆ H (2D)0 = ˆ p + 14 (cid:0) η x ˆ x + ˆ z (cid:1) , (50) where η x = ω x /ω z is the aspect ratio, as in Sec. II.In Eq. (50), ˆ p and ˆ ρ ≡ (ˆ x, ˆ z ) are the relative mo-mentum and coordinate operators, respectively. We fur-ther assume that the two atoms experience a 2D s -wavezero-range interaction with 2D scattering length a , andmodel this potential with the 2D Bethe-Perierls bound-ary condition (BPC). Explicitly, in our calculation theeigen-energy E and corresponding eigen-state | Ψ (cid:105) forthese two interacting atoms should satisfy the equation (cid:104) ρ | ˆ H (2D)0 | Ψ (cid:105) = E (cid:104) ρ | Ψ (cid:105) ; for ρ > , (51)and the BPC [50]lim ρ → (cid:104) ρ | Ψ (cid:105) ∝ (ln ρ − ln a ) , (52)where | ρ (cid:105) is the eigen-state of the relative-coordinate op-erator ˆ ρ , with the corresponding eigen-energy ρ , and ρ = | ρ | . Namely, the wave function (cid:104) ρ | Ψ (cid:105) satisfies theeigen-equaiton of the free Hamiltonian ˆ H (2D)0 in the re-gion other than the origin (i.e., the region with ρ > ρ → ρ → − ρ , and thus in this work we only considerthe eigen-energies in this subspace. In addition, in thelong-range limit ρ → ∞ the wave function (cid:104) ρ | Ψ (cid:105) shouldsatisfy lim ρ →∞ (cid:104) ρ | Ψ (cid:105) = 0 . (53) A. Expression of J ( E ) The solution to Eq. (51) and the long-range condition(53) is proportional to the 2D Green’s function, i.e., (cid:104) ρ | Ψ (cid:105) ∝ G (2D)0 ( E ; ρ ) ≡ (cid:104) ρ | E − ˆ H (2D)0 | (cid:105) . (54)Therefore, we can derive the algebrac equation for theeigen-energy E by matching G (2D)0 ( E ; ρ ) with the BPC(52). To this end we need to expand G (2D)0 ( E ; ρ ) in thelimit ρ →
0. This expansion can be done with the sim-ilar approach as in Sec. II, and we show the detail inAppendix C. Finally we obtainlim ρ → G (2D)0 ( E ; ρ ) = 12 π ln ρ − (cid:20) W ( E ) + (cid:90) + ∞ I ( E, β ) dβ (cid:21) , (55)where the functions W ( E ) and I ( E, β ) being defined as W ( E ) = √ η x (cid:88) n x ∈ C (2D) E Γ (cid:104) − ( E − (cid:15) nx )2 (cid:105) n x − Γ (cid:0) − n x (cid:1) Γ( n x + 1)Γ (cid:104) − ( E − (cid:15) nx )2 (cid:105) − γ π − π ln (cid:16) κ (cid:17) , (56)and I ( E, β ) = e βE (cid:89) α = x,z (cid:114) η α π sinh ( η α β ) − πβ e − κβ − (cid:114) η x β (cid:88) n x ∈ C (2D) E (cid:40) n x − Γ (cid:0) − n x (cid:1) Γ( n x + 1) e β ( E − (cid:15) nx ) (cid:41) , (57)respectively, in our natural unit. Here γ = 0 . ... is the Euler’s constant, η z = 1, (cid:15) n x = ( n x + 1 / η x , the parameter κ could be any positive number and the result is independent of the value of κ . In Eq. (57) C (2D) E is a number setdefined as C (2D) E : (cid:26) n x | n x = 0 , , , , ... ; (cid:15) n x + 12 ≤ E (cid:27) . (58)Clearly, for E < E the set C (2D) E is empty and the summation (cid:80) n x ∈ C (2D) E in the expressions (56, 57) of W and I becomes zero.Combining Eq. (C8), Eq. (54) and Eq. (52), we obtain the equation for the eigen-energy E of our 2D system,which has the form of Eq. (2): J ( E ) = ln a . (59)Here the function J ( E ) is given by J ( E ) = 4 π (cid:20) W ( E ) + (cid:90) + ∞ I ( E, β ) dβ (cid:21) , (60)with W ( E ) and I ( E, β ) being defined in Eq. (56) and Eq. (57), respectively. One can derive the eigen-energiesof the relative motion of these two atoms by solving Eq. (59), and further obtain the corresponding eigen-states.As in the 3D cases, the summations in Eq. (56, 57) are done for finite terms, and the integration in Eq. (60)converges. Therefore, using Eqs. (60, 56, 57) one can efficiently calculate J ( E ). In addition, for the systems with a2D isotropic trap ( ω x = ω z or η x = 1), the equaiton (59) for the eigen-energy and the expression (60) for the function J ( E ) are equivalent to the results derived by Ref. [1]. Explicitly, we have [51] J ( E ) = ddz ln [Γ( z )] (cid:12)(cid:12)(cid:12) z = − E − γ + 12 ln 2 , (for ω x = ω z ) . (61) B. Techniques for Fast Calculation of J ( E ) The techniques shown in Sec. II.B for the fast calculation of J ( E ) can also be directly generalized to the 2D case(Appendix D). In particular, as proved in Appendix D, J ( E ) given by Eq. (60) can be re-expressed as J ( E ) = 4 π (cid:34)(cid:90) Λ0 A ( E, β ) dβ + B (1)2D ( E, Λ) + B (2)2D ( E, Λ) − Γ(0 , κ
Λ)8 π (cid:35) , (62)where Λ and κ are arbitrary finite positive numbers, as in Eqs. (56, 57, 40), Γ[ a, z ] is the incomplete Gamma functionand the functions A ( E, β ) and B (1 , ( E, Λ) are defined as A ( E, β ) = e βE (cid:89) α = x,z (cid:114) η α π sinh ( η α β ) − πβ e − κβ ; (63) B (1)2D ( E, Λ) = − γ π − π ln (cid:16) κ (cid:17) + (cid:88) n x ∈ C (2D) E √ η x n x − e ( E − (cid:15) nx − )Λ √ e − · Γ (cid:16) − E − (cid:15) nx (cid:17) Γ (cid:0) − n x (cid:1) Γ( n x + 1)Γ (cid:16) − E − (cid:15) nx (cid:17) × F (cid:20) , − E − (cid:15) n x , − E − (cid:15) n x , e − (cid:21) ; (64)0and B (2)2D ( E, Λ) = − (cid:88) n x / ∈ C (2D) E n x − · √ η x csch ΛΓ (cid:0) − n x (cid:1) Γ( n x + 1) · e ( E − (cid:15) nx − ( e − × F (cid:104) , − E − (cid:15) nx , − E − (cid:15) nx , e − (cid:105) E − (cid:15) n x ) − , (65)respectively, with F being the Hypergeomtric function. As in the 3D cases, the expression (62) of J ( E ) has theadvantages (A) and (B) shown in Sec. II. B. Therefore, the numerical calculations based on Eq. (62) is quite possiblyto be faster then the one based on Eq. (60), especially for the high-energy cases with large E . C. Energy Sepectrum and Eigen-States
In Fig. 3 (a) we illustrate the energy spectrum for thecases with an isotropic 2D trap ( η x = 1), and show thatour results are same as the ones from Ref. [1]. In Fig. 3(b) and (c) we show the results for anisotropic 2D trapswhich have similar ( η x = 1 .
1) and quite different ( η x = 5)frequencies in the x - and z -directions.In addition, similar to Sec. II.D, one can also derive thecorresponding eigen-state | Ψ (cid:105) with eigen-energy E , whichsatisfies Eq. (51) and the boundary conditions (52, 53),as well as the normalization condition (cid:104) Ψ | Ψ (cid:105) = 1. WithEq. (54) and the similar approach as in Sec. II. D weobtain (cid:104) ρ | Ψ (cid:105) = 1 (cid:113)(cid:80) m | d m | (cid:88) m d m φ m x ( η x , x ) φ m z (1 , z ) , (66)with d m = 1 (cid:16) E − E (2D) m (cid:17) φ m x ( η x , φ m z (1 , . (67)Here ( x, z ) are the components of ρ , m = ( m x , m y ) with m x,y = 0 , , , ... , E (2D) m = ( m x + 1 / η x + ( m z + 1 / φ n ( η, X ) is defined in Eq. (A3). IV. GENERALIZATION OF OUR APPROACHTO OTHER PROBLEMS
At the end of this work, we briefly summarize the mainideas of our approach used in the above calculations, anddiscuss how to generalize these ideas. Here we take the3D systems as an example.The key step for solving the two-body problem withzero-rang interaction is the calculation of the free Green’soperator ˆ G ( E ) ≡ / [ E − ˆ H ]. When E is less than theground-state energy E of ˆ H , ˆ G ( E ) can be directlyexpressed as the Laplace transform of of the imaginary-time evolution operator e − β ˆ H , i.e., we have ˆ G ( E ) = − (cid:82) + ∞ e β ( E − ˆ H ) dβ . However, when E > E this in-tegration diverges. To solve this problem, we separateˆ G ( E ) into two parts, i.e., the contributions from thehigh-energy states with n ∈ U E and the ones from the low-energy states with n ∈ L E , as shown in Eq. (28).Furthermore, as shown in Eqs. (29, 30), the first partcan still be expressed as the Laplace transform of the op-erator e − β ˆ H (3D)0 − (cid:80) n ∈ L E | n (cid:105)(cid:104) n | e − β ( E n ) , which convergesfor for any non-zero z , and the second part only includesfinite terms.It is pointed out that, the definitions of the sets L E and U E are not unique. The only requirements are:(i) : L E is the complement of U E .(ii) : E n > E for ∀ n ∈ U E . Namely, the set L E in-cludes (but is not limited to) all n which satisfies E n < E .For instance, an alternative definition of these two setsis: U E : { n | E n > E } and L E : { n | E n ≤ E } .Using this method we can derive the helpful expression(31) of Green’s function G (0) ( E, z e z ), which is just thematrix element (cid:104) z e z | ˆ G ( E ) | (cid:105) of the free Green’s oper-ator ˆ G ( E ). In this expression there is just a summa-tion for finite terms and a one-dimensional integration (cid:82) + ∞ dβ [ K ( z ; E, β ) − F ( z ; E, β )] which converges for anyfinite z .To complete the calculation we still require to removethe divergent term of the above integration in the limit z →
0. This divergent term is contributed by the lead-ing term e − z β / (4 πβ ) of the integrand K ( z ; E, β ) − F ( z ; E, β ) in the limit β → + . Therefore, we can re-move it via the technique used in Eqs. (18, 19, 20).Our above approach for the calculation of the two-body free Green’s function can be directly generalized toother few-body problems of ultracold atoms, especiallythe ones where the analytical expressions of the eigen-states and imaginary-time propagator of the free Hamil-tonian are known, e.g. , the few-body problems in mixed-dimensional systems [52]. Here we emphasize that, withthe help of the Laplace transformation for imaginary-time propagator, the free Green’s operator can alwaysbe expressed as a one-dimensional integration, no matterhow many degrees of freedoms are involved in the system.Thus, the free Green’s function given by our method al-ways includes a converged one-dimensional integrationand summations for finite terms.1 - 1 5 - 1 0 - 5 0 5 1 0 1 5- 5051 01 5 (cid:1) x = 1 ( a ) T h i s w o r k
R e f . [ 1 ] E l n a - 1 5 - 1 0 - 5 0 5 1 0 1 5- 4- 20246 (cid:1) x = 1 . 1 E l n a ( b ) - 1 5 - 1 0 - 5 0 5 1 0 1 5- 5051 01 52 0 (cid:1) x = 5 E l n a ( c ) FIG. 3. (color online) The energy spectrum of the relative motion of two atoms in a 2D harmonic trap with η x = 1 (a), η x = 1 . η x = 5 (c). We show the results given by Eq. (2) and the expression (60) of the function J ( E ) with solid lines. In(a) we also show the results given by Ref. [1] with blue dots. Here we use the natural unit (cid:126) = 2 µ = ω z = 1. V. SUMMARY
In this work we derive the algebraic equations for theeigen-energies of two atoms in a 2D or 3D harmonic trap.Our results is applicable for general cases, no matter ifthe trap is completely anisotropic or has spherical or ax-ial symmetry. Using our results one can easily derive thecomplete energy spectrum, which can be used for the fur-ther theoretical or experimental studies of dynamical orthermodynamical problems. Our approach can be usedin other few-body problems of confined ultracold atoms.
Note added:
When we finished this work, we realizedthat recently there is a related work [53]. The authors derived the expression of J ( E ) for E < E , and a re-currence relation of J ( E ) for arbitrary E . With this re-currence relation they also obtained the complete energyspectrum of two atoms in a 2D anisotropic confinement,as well as the eigen-states. ACKNOWLEDGMENTS
This work is supported by the National Key R&DProgram of China (Grant No. 2018YFA0306502 (PZ),2018YFA0307601 (RZ)), NSFC Grant No. 11804268(RZ), 11434011(PZ), 11674393(PZ), as well as the Re-search Funds of Renmin University of China under GrantNo. 16XNLQ03(PZ).
Appendix A: Proof of Eqs. (36) and (37)
In this appendix we prove Eqs. (36) and (37) in Sec. II.A. To this end, we first show some results on 1D harmonicoscillator, which will be used in our calculation.2
1. Some properties of 1D harmonic oscillator
Let us consider a 1D harmonic oscillator with frequency η and mass µ . The Hamiltonian of this oscillator is( (cid:126) = 2 µ = 1) ˆ H ho = ˆ P + η ˆ X , (A1)with ˆ X and ˆ P being the coordinate and momentum operator, respectively. The eigen-energy of ˆ H ho is E n = (cid:18) n + 12 (cid:19) η ; ( n = 0 , , , ... ) , (A2)and the wave function of the eigen-state | n (cid:105) corresponding to E n can be expressed as (cid:104) X | n (cid:105) ≡ φ n ( η, X ) = (cid:16) η π (cid:17) e − ηX (cid:112) n Γ( n + 1) H n (cid:18)(cid:114) η X (cid:19) , (A3)with | X (cid:105) being the eigen-state of the position operator ˆ X with eigen-value X , H n ( X ) and Γ( α ) being the Hermitianpolynomial and the Gamma function, respectively. The wave function φ n ( η, X ) also satisfies (cid:104) X | e − β ˆ H ho | (cid:105) = + ∞ (cid:88) n =0 φ n ( η, X ) φ ∗ n ( η, e − β ( n +1 / η = (cid:114) η π sinh( ηβ ) exp (cid:20) − ηX coth ( ηβ )4 (cid:21) (A4)for β > g ( ξ ; η ; X ) of the 1D harmonic oscillator, which is defined as g ( ξ ; η ; X ) ≡ (cid:104) X | ξ − ˆ H ho | (cid:105) = + ∞ (cid:88) n =0 φ n ( η, X ) φ ∗ n ( η, ξ − E n . (A5)This function satisfies the differential equation ξ · g ( ξ ; η ; X ) + d dX g ( ξ ; η ; X ) − η X g ( ξ ; η ; X ) = δ ( X ) (A6)and the boundary condition lim | X |→∞ g ( ξ ; η ; X ) = 0 . (A7)To derive g ( ξ ; η ; X ), we can first solve the equation (A6) in the regions X >
X < X = 0, which is given by the term δ ( X ) in Eq.(A6). With this approach we obtain g ( ξ ; η ; X ) = + ∞ (cid:88) n =0 φ n ( η, X ) φ ∗ n ( η, ξ − E n = − Γ( − ξ η )2 + ξ η √ πη D ξ − η/ η ( √ ηX ) , (A8)where D λ ( α ) is the parabolic cylinder function. Eq. (A8) and the property of the parabolic cylinder function furtheryields g ( ξ ; η ; 0) = + ∞ (cid:88) n =0 | φ n ( η, | ξ − ( n + 1 / η = − Γ( − ξ η )2 √ η Γ( − ξ η ) . (A9)
2. Proof of the two equations
Now we prove Eq. (36) and Eq. (37) in our main text, which are about the expressions of the functions W ( E )and I ( E, β ), respectively. As shown in Sec. II.A, these two functions are defined as W ( E ) ≡ Q (0 , E, β ) , (A10)3and I ( E, β ) ≡ − ˜ K (0; E, β ) + F (0; E, β ) , (A11)with the functions Q ( z, E, β ), ˜ K ( z ; E, β ) defined in Eq. (32) and Eq. (19), respectively. Thus, we have Q (0; E, β ) = (cid:88) n ∈ L E (cid:104) | n (cid:105)(cid:104) n | (cid:105) E − E n ; (A12) F (0; E, β ) = (cid:88) n ∈ L E (cid:104) | n (cid:105)(cid:104) n | (cid:105) e − β ( E n − E ) . (A13)It is clear that the eigen-state | n (cid:105) of the 3D free Hamiltonian ˆ H , which is defined in Sec. II.A, satisfies (cid:104) | n (cid:105) = 0when n x or n y is odd. Using this fact and the definitions of the sets L E and C (3D) E , which are given in Eqs. (25, 38)of our main text, we obtain Q (0; E, β ) = (cid:88) ( n x ,n y ) ∈ C (3D) E + ∞ (cid:88) n z =0 (cid:104) | n (cid:105)(cid:104) n | (cid:105) E − E n = (cid:88) ( n x ,n y ) ∈ C (3D) E (cid:40) | φ n x ( η x , | | φ n y ( η y , | (cid:34) + ∞ (cid:88) n z =0 | φ n z ( η z , | (cid:0) E − (cid:15) n x − (cid:15) n y (cid:1) − (cid:15) n z (cid:35)(cid:41) , (A14)and F (0; E, β ) = (cid:88) ( n x ,n y ) ∈ C (3D) E + ∞ (cid:88) n z =0 (cid:104) | n (cid:105)(cid:104) n | (cid:105) e − β ( E n − E ) = (cid:88) ( n x ,n y ) ∈ C (3D) E (cid:40) | φ n x ( η x , | | φ n y ( η y , | e β ( E − (cid:15) nx − (cid:15) ny ) (cid:34) + ∞ (cid:88) n z =0 | φ n z ( η z , | e − β(cid:15) nz (cid:35)(cid:41) , (A15)where η x,y,z and (cid:15) n x,y,z are defined in Sec. II, and the function φ n ( η, X ) is defined in Eq. (A3).Substituting Eq. (A9) into Eq. (A14) and then into Eq. (A10), and using the property H n (0) = 2 n √ π Γ( − n ) (A16)of the Hermitian polynomial, we can derive Eq. (36). Moreover, Substituting Eq. (A4) into Eq. (A15) and then intoEq. (A11), and using Eqs. (A16), (19) and (17), we can derive Eq. (37). Appendix B: Proof of Eq. (40)
In this appendix we prove Eq. (40) in Sec. II.B. To this end, we separate the integration in Eq. (39) into two parts,i.e., (cid:90) ∞ I ( E, β ) dβ = (cid:90) Λ0 I ( E, β ) dβ + (cid:90) ∞ Λ I ( E, β ) dβ, (B1)4with Λ being an arbitrary finite positive number. In addition, using the definition (37) of I ( E, β ), we immediatelyobtain (cid:90) Λ0 I ( E, β ) dβ = (cid:90) Λ0 A ( E, β ) dβ + (cid:88) ( n x ,n y ) ∈ C (3D) E (cid:90) Λ0 n x + n y − √ πη x η y e β ( E − (cid:15) nx − (cid:15) ny ) Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1) √ sinh β dβ = (cid:90) Λ0 A ( E, β ) dβ + (cid:88) ( n x ,n y ) ∈ C (3D) E n x + n y − √ πη x η y Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1) (cid:40) √ π Γ (cid:16) − E − (cid:15) nx − (cid:15) ny (cid:17) Γ (cid:16) − E − (cid:15) nx − (cid:15) ny (cid:17) − Γ (cid:16) − E − (cid:15) nx − (cid:15) ny (cid:17) Γ (cid:16) − E − (cid:15) nx − (cid:15) ny (cid:17) e ( E − (cid:15) nx − (cid:15) ny − )Λ (cid:112) e − × F (cid:20) , − E − (cid:15) n x − (cid:15) n y , − E − (cid:15) n x − (cid:15) n y , e − (cid:21)(cid:41) (B2)with A ( E, β ) being defined in Eq. (41).Now we calculate the term (cid:82) ∞ Λ I ( E, β ) dβ in Eq. (B1). We first notice that, according to Eqs. (A11,A15,19,17), I ( E, β ) can be re-expressed as I ( E, β ) = (cid:18) πβ (cid:19) − (cid:88) n x ,n y ,n z (cid:104) | n (cid:105)(cid:104) n | (cid:105) e − β ( E n − E ) + (cid:88) ( n x ,n y ) ∈ C (3D) E + ∞ (cid:88) n z =0 (cid:104) | n (cid:105)(cid:104) n | (cid:105) e − β ( E n − E ) = (cid:18) πβ (cid:19) − (cid:88) ( n x ,n y ) / ∈ C (3D) E + ∞ (cid:88) n z =0 (cid:104) | n (cid:105)(cid:104) n | (cid:105) e − β ( E n − E ) = (cid:18) πβ (cid:19) − (cid:88) ( n x ,n y ) / ∈ C (3D) E (cid:40) | φ n x ( η x , | | φ n y ( η y , | e β ( E − (cid:15) nx − (cid:15) ny ) (cid:34) + ∞ (cid:88) n z =0 | φ n z ( η z , | e − β(cid:15) nz (cid:35)(cid:41) , (B3)where (cid:15) n x,y,z is defined in Sec. II, and the function φ n ( η, X ) is defined in Eq. (A3). Moreover, Substituting Eq. (A4)and Eq. (A16) into Eq. (B3), we further obtain I ( E, β ) = (cid:18) πβ (cid:19) − (cid:114) πη x η y β (cid:88) ( n x ,n y ) / ∈ C (3D) E n x + n y − Γ (cid:0) − n x (cid:1) Γ (cid:16) − n y (cid:17) Γ( n x + 1)Γ( n y + 1) e β ( E − (cid:15) nx − (cid:15) ny ) , (B4)where η x,y,z is defined in Sec. II. Doing the integration (cid:82) ∞ Λ I ( E, β ) dβ in both sides of Eq. (B4), we further obtain (cid:90) ∞ Λ I ( E, β ) dβ = B (2)3D ( E, Λ) + (cid:18) π (cid:19) / √ , (B5)where B (2)3D ( E, Λ) is defined in Eq. (43).Substituting Eqs. (B2, B5) into Eq. (B1) and then into Eq. (39), and further using Eqs. (36), we can derive Eq.(40).
Appendix C: Proof of Eqs. (56) and (57)
In this appendix we prove Eqs. (56) and (57) in Sec. III, which is related to the behavior of the 2D Green’s function G (2D)0 ( E ; ρ ) in the limit | ρ | →
0. Our approach is similar to the method used in Sec. II.A.5As in Sec. II and Appendix A, by choosing ρ = ρ e z ( ρ >
0) and making direct calculations we can find that G (2D)0 ( E ; ρ e z ) = − (cid:90) ∞ [ K ( ρ ; E, β ) − Y ( ρ ; E, β )] dβ + Z ( ρ ; E ) , (C1)which is similar to the expression (31) of the 3D free Green’s function. Here the functions K ( ρ ; E, β ), Y ( ρ ; E, β )and Z ( ρ ; E ) are defined as K ( ρ ; E, β ) = e βE (cid:104) ρ e z | e − β ˆ H (2 D )0 | (cid:105) = (cid:32) (cid:89) α = x,z (cid:114) η α π sinh ( η α β ) (cid:33) · exp (cid:20) βE − coth ( β )4 ρ (cid:21) ; (C2) Y ( ρ ; E, β ) = (cid:88) n x ∈ C (2 D ) E + ∞ (cid:88) n z =0 e β ( E − E n ) (cid:104) ρ e z | n (cid:105)(cid:104) n | (cid:105) ; (C3) Z ( ρ ; E ) = (cid:88) n x ∈ C (2 D ) E + ∞ (cid:88) n z =0 (cid:104) ρ e z | n (cid:105)(cid:104) n | (cid:105) E − E n , (C4)where η x,z and (cid:15) n x,z have the same definition as in Sec. II and C (2 D ) E is defined in Eq. (58). Furthermore, theintegration in Eq. (C5) diverges in the limit z →
0, and we can remove this divergence via the similar approach as inSec. III. To this end, we re-express Eq. (C5) as G (2D)0 ( E ; ρ e z ) = − (cid:90) ∞ πβ · exp (cid:18) − βκ − β ρ (cid:19) dβ − (cid:90) ∞ (cid:104) ˜ K ( ρ ; E, β ) − Y ( ρ ; E, β ) (cid:105) dβ + Z ( ρ ; E ) , (C5)with κ being any positive number and˜ K ( ρ ; E, β ) = K ( ρ ; E, β ) − πβ · exp (cid:18) − βκ − β ρ (cid:19) . (C6)Furthermore, using the result (cid:90) ∞ πβ · exp (cid:18) − βκ − β ρ (cid:19) dβ = − π ln ρ − π γ − π ln κ O ( ρ ) (for ρ > , (C7)with γ = 0 . ... being the Euler’s constant, we derive a result with the same form of Eq. (C8):lim ρ → G (2D)0 ( E ; ρ ) = 12 π ln ρ − (cid:20) W ( E ) + (cid:90) + ∞ I ( E, β ) dβ (cid:21) . (C8)In this step the functions W ( E ) and I ( E, β ) are given by W ( E ) = − π γ − π ln κ − Z (0 , E ); (C9) I ( E, β ) = 12 (cid:104) ˜ K (0; E, β ) − Y (0; E, β ) (cid:105) . (C10)Moreover, with the method in Appendix A, we can derive the alternative expressions of W ( E ) and I ( E, β ), i.e.,Eqs. (56) and (57).
Appendix D: Techniques for fast calculation of J ( E ) In this appendix we generalize the techniques shown in Sec.II.B and Appendix B to the 2D case. We first generalizeEqs. (40-43) to the 2D cases and prove Eq. (62). This can be done via direct calculations with the method shown inAppendix B. We separate the integration in Eq. (60) into two parts, i.e., (cid:90) ∞ I ( E, β ) dβ = (cid:90) Λ0 I ( E, β ) dβ + (cid:90) ∞ Λ I ( E, β ) dβ, (D1)6with Λ being an arbitrary positive number. In addition, using the definition (57) of I ( E, β ), we immediately obtain (cid:90) Λ0 I ( E, β ) dβ = (cid:90) Λ0 A ( E, β ) dβ − (cid:88) n x ∈ C (2D) E (cid:90) Λ0 n x − √ η x e β ( E − (cid:15) nx ) Γ (cid:0) − n x (cid:1) Γ( n x + 1) √ sinh β dβ = (cid:90) Λ0 A ( E, β ) dβ + (cid:88) ( n x ) ∈ C (2D) E n x − √ η x Γ (cid:0) − n x (cid:1) Γ( n x + 1) (cid:40) − √ π Γ (cid:16) − E − (cid:15) nx (cid:17) Γ (cid:16) − E − (cid:15) nx (cid:17) +Γ (cid:16) − E − (cid:15) nx (cid:17) Γ (cid:16) − E − (cid:15) nx (cid:17) e ( E − (cid:15) nx − )Λ (cid:112) e − × F (cid:20) , − E − (cid:15) n x , − E − (cid:15) n x , e − (cid:21)(cid:41) (D2)with A ( E, β ) being defined in Eq. (63).Furthermore, using the method in Appendix B, we find that the function I ( E, β ) defined in Eq. (57) has analternative expression I ( E, β ) = − πβ e − κβ + (cid:114) η x β (cid:88) n x / ∈ C (2D) E (cid:40) n x − Γ (cid:0) − n x (cid:1) Γ( n x + 1) e β ( E − (cid:15) nx ) (cid:41) , (D3)which is similar to Eq. (B4). Thus, doing the integration (cid:82) ∞ Λ I ( E, β ) dβ in both sides of Eq. (D3), we further obtain (cid:90) ∞ Λ I ( E, β ) dβ = B (2)2D ( E, Λ) − Γ(0 , κ Λ )8 π (D4)where Γ[ a, z ] is the incomplete Gamma function and B (2)2D ( E, Λ) is defined in Eq. (65).As in Appendix B, substituting Eqs. (D2, D4) into Eq. (D1) and then into Eq. (60), and further using Eqs. (56),we can derive Eq. (62).In addition, the second technique shown in Sec.II.B is the one based on Eqs. (44, 45). It is clear that this techniquecan be directly generalized to the 2D case. [1] T. Busch, B.-G. Englert, K. Rza˙zewski, and M. Wilkens,Found. Phys. , 549 (1998).[2] Z. Idziaszek and T. Calarco, Phys. Rev. A , 050701(2005).[3] Z. Idziaszek and T. Calarco, Phys. Rev. A , 022712(2006).[4] J.-J. Liang and C. Zhang, Phys. Scr. , 025302 (2008).[5] D. Blume, Rep. Prog. Phys. , 046401 (2012).[6] S. Grishkevich and A. Saenz, Phys. Rev. A , 013403(2009).[7] S. Grishkevich, S. Sala, and A. Saenz, Phys. Rev. A ,062710 (2011).[8] S. Sala, G. Z¨urn, T. Lompe, A. N. Wenz, S. Murmann,F. Serwane, S. Jochim, and A. Saenz, Phys. Rev. Lett. , 203202 (2013).[9] S. Sala and A. Saenz, Phys. Rev. A , 022713 (2016).[10] B. Sun, W. X. Zhang, S. Yi, M. S. Chapman, and L. You,Phys. Rev. Lett. , 123201 (2006).[11] L. M. A. Kehrberger, V. J. Bolsinger, and P. Schmelcher,Phys. Rev. A , 013606 (2018).[12] G. Bougas, S. I. Mistakidis, and P. Schmelcher, Phys. Rev. A , 053602 (2019).[13] L. Budewig, S. I. Mistakidis, and P. Schmelcher,Molecular Physics , 2043 (2019),https://doi.org/10.1080/00268976.2019.1575995.[14] A. N. Wenz, G. Z¨urn, S. Murmann, I. Brouzos, T. Lompe,and S. Jochim, Science , 457 (2013).[15] D. Blume and C. H. Greene, Phys. Rev. A , 013601(2002).[16] X.-J. Liu, H. Hu, and P. D. Drummond, Phys. Rev. Lett. , 160401 (2009).[17] J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev.A , 043619 (2008).[18] K. M. Daily and D. Blume, Phys. Rev. A , 053615(2010).[19] S. E. Gharashi, K. M. Daily, and D. Blume, Phys. Rev.A , 042702 (2012).[20] X.-J. Liu, Phys. Rep. , 37 (2013).[21] X. Y. Yin, D. Blume, P. R. Johnson, and E. Tiesinga,Phys. Rev. A , 043631 (2014).[22] X. Y. Yin, D. Blume, P. R. Johnson, and E. Tiesinga,Phys. Rev. A , 043631 (2014). [23] S.-G. Peng, S.-H. Zhao, and K. Jiang, Phys. Rev. A ,013603 (2014).[24] X. Y. Yin and D. Blume, Phys. Rev. A , 013608 (2015).[25] D. Blume, M. W. C. Sze, and J. L. Bohn, Phys. Rev. A , 033621 (2018).[26] X. Y. Yin, H. Hu, and X.-J. Liu, Phys. Rev. Lett. ,013401 (2020).[27] F. Scazza, C. Hofrichter, M. H¨ofer, P. C. De Groot,I. Bloch, and S. F¨olling, Nat. Phys. , 779 (2014).[28] G. Cappellini, M. Mancini, G. Pagano, P. Lombardi,L. Livi, M. Siciliani de Cumis, P. Cancio, M. Pizzocaro,D. Calonico, F. Levi, C. Sias, J. Catani, M. Inguscio,and L. Fallani, Phys. Rev. Lett. , 120402 (2014).[29] M. A. Norcia, A. W. Young, and A. M. Kaufman, Phys.Rev. X , 041054 (2018).[30] A. Cooper, J. P. Covey, I. S. Madjarov, S. G. Porsev,M. S. Safronova, and M. Endres, Phys. Rev. X , 041055(2018).[31] G. Cappellini, L. F. Livi, L. Franchi, D. Tusi, D. Bene-dicto Orenes, M. Inguscio, J. Catani, and L. Fallani,Phys. Rev. X , 011028 (2019).[32] Q. Guan, V. Klinkhamer, R. Klemt, J. H. Becher,A. Bergschneider, P. M. Preiss, S. Jochim, and D. Blume,Phys. Rev. Lett. , 083401 (2019).[33] L. R. Liu, J. D. Hood, Y. Yu, J. T. Zhang, N. R. Hutzler,T. Rosenband, and K.-K. Ni, Science , 900 (2018).[34] L. Anderegg, L. W. Cheuk, Y. Bao, S. Burchesky,W. Ketterle, K.-K. Ni, and J. M. Doyle, Science ,1156 (2019).[35] J. D. Hood, Y. Yu, Y.-W. Lin, J. T. Zhang, K. Wang,L. R. Liu, B. Gao, and K.-K. Ni, “Multichannel in-teractions of two atoms in an optical tweezer,” (2019),arXiv:1907.11226.[36] K. Wang, X. He, R. Guo, P. Xu, C. Sheng, J. Zhuang,Z. Xiong, M. Liu, J. Wang, and M. Zhan, Phys. Rev. A , 063429 (2019).[37] D. E. Chang, J. S. Douglas, A. Gonz´alez-Tudela, C.-L. Hung, and H. J. Kimble, Rev. Mod. Phys. , 031002(2018).[38] Y. Meng, A. Dareau, P. Schneeweiss, and A. Rauschen-beutel, Phys. Rev. X , 031054 (2018).[39] A. Dareau, Y. Meng, P. Schneeweiss, and A. Rauschen-beutel, Phys. Rev. Lett. , 253603 (2018).[40] C. Ospelkaus, S. Ospelkaus, L. Humbert, P. Ernst,K. Sengstock, and K. Bongs, Phys. Rev. Lett. , 120402(2006).[41] T. St¨oferle, H. Moritz, K. G¨unter, M. K¨ohl, andT. Esslinger, Phys. Rev. Lett. , 030401 (2006).[42] M. J. Mark, E. Haller, K. Lauber, J. G. Danzl, A. J.Daley, and H.-C. N¨agerl, Phys. Rev. Lett. , 175301(2011).[43] L. Riegger, N. Darkwah Oppong, M. H¨ofer, D. R. Fer-nandes, I. Bloch, and S. F¨olling, Phys. Rev. Lett. ,143601 (2018).[44] R. Chapurin, X. Xie, M. J. Van de Graaff, J. S. Popowski,J. P. D’Incao, P. S. Julienne, J. Ye, and E. A. Cornell,Phys. Rev. Lett. , 233402 (2019).[45] P. Massignan and Y. Castin, Phys. Rev. A , 013616(2006).[46] R. Zhang and P. Zhang, Phys. Rev. A , 043627 (2018).[47] R. Zhang and P. Zhang, Phys. Rev. A , 063607(2019).[48] D. Xiao, R. Zhang, and P. Zhang, Few-Body Systems , 63 (2019).[49] R. Zhang and P. Zhang, Phys. Rev. A , 013636(2020).[50] B. J. Verhaar, J. P. H. W. van den Eijnde, M. A. J.Voermans, and M. M. J. Schaffrath, Journal of PhysicsA: Mathematical and General , 595 (1984).[51] Notice that the definitions of the 2D scattering length inRef. [1] and in this work are different.[52] Y. Nishida and S. Tan, Phys. Rev. Lett.101