One-dimensional Fermi gas with a single impurity in a harmonic trap: Perturbative description of the upper branch Seyed Ebrahim Gharashi, X. Y. Yin, Yangqian Yan, and D. Blume1 1Department of Physics and Astronomy, Washington State University, Pullman, Washington 99164-2814, USA (Dated: January 20, 2015) The transition from “few to many” has recently been probed experimentally in an ultracold 5 harmonically confined one-dimensional lithium gas, in which a single impurity atom interacts with 1 a background gas consisting of one, two, or more identical fermions [A. N. Wenz et al., Science 0 342, 457 (2013)]. For repulsive interactions between the background or majority atoms and the 2 impurity, the interaction energy for relatively moderate system sizes was analyzed and found to n convergetoward thecorresponding expression for an infinitelylarge Fermigas. Motivated bythese a experimentalresults,weapplyperturbativetechniquestodeterminetheinteractionenergyforweak J and strong coupling strengths and derive approximate descriptions for the interaction energy for repulsiveinteractions with varyingstrength between theimpurityand themajority atoms and any 9 1 numberof majority atoms. PACSnumbers: ] s a g I. INTRODUCTION perturbation theory to harmonically confined systems - t and derive approximate solutions whose accuracy can n One-dimensionalBoseandFermisystemswithcontact be improved systematically by considering successively a u interactions have been studied for many decades now, higher orders in the expansion in the small parame- q especially in the regime where the systems obey peri- ter. We focus on one-dimensional Fermi gases with a . odic boundary conditions [1–5]. A large fraction of the single impurity under external harmonic confinement. t a eigenstates can be thought of as corresponding to gas- This system is of particular interest since it has been m like states. A second subset of eigenstates corresponds realized experimentally in Jochim’s cold atom labora- - to self-bound droplet-like states. These states maintain tory [18, 19]. In the experiments, the impurity is a d their bound state character in the absence of periodic lithium atom that occupies a hyperfine state different n o boundary conditions, i.e., in free space. In many cases, from the hyperfine state that the majority atoms oc- c both the gas-like and droplet-like states can be obtained cupy. The trapping geometry is highly-elongatedand ef- [ analyticallyviatheBetheansatz. TheBetheansatztakes fectively one-dimensional. We willshowthatour pertur- advantage of the fact that the zero-range nature of the bativeresultsenableustocalculatetheenergyoftheup- 1 v interactions,combinedwiththefactthatparticlesinone per branch,which has beenstudied experimentally, with 9 dimension have to pass through each other to exchange fairly good accuracy for all N over a wide range of cou- 6 positions,allowsonetoidentifyconstantsofmotion. The pling strengths. In addition, our results provide bounds 3 solutionscanthenbe derivedinterms oftheseconstants on the energies in the weakly- and strongly-interacting 4 of motion. A closely related aspect is that a variety of regimes. These bounds can, e.g., be used to assess the 0 one-dimensional systems with two-body contact interac- accuracy of numerical solutions. . 1 tions are integrable [1, 4]. The remainder of this paper is organized as follows. 0 The solution of the homogeneous system can be ap- Section II introduces the system Hamiltonian and nota- 5 plied to one-dimensional systems under spatially vary- tion. Section III summarizes our perturbative results. 1 ing external confinement via the local density approxi- The perturbative results areanalyzedinSecs. IV andV. : v mation [6–10]. This approximation typically provides a Finally, Sec. VI concludes. i X highly accurate description for a large number of par- ticles but not necessarily for a small number of parti- r II. SYSTEM HAMILTONIAN a cles. It is thus desirable to derive moreaccuratedescrip- tions for small one-dimensional systems with two-body delta-function interactions under external confinement. We consider a single impurity immersed in a one- Unfortunately, extensions of the Bethe ansatz to inho- dimensional Fermi gas that consists of N identical mass mogeneous systems are, in general, not known. This m fermions. The mass of the impurity is equal to that can be understood intuitively by realizing that the rel- of the majority or background particles. The impurity, ative two-body momentum in inhomogeneous systems is with position coordinate z , interacts with the majority 0 not conserved due to the presence of the spatially vary- particles, with position coordinates z (j = 1, ,N), j ··· ingconfinement. Correspondingly,harmonicallytrapped through a zero-range two-body potential with strength one-dimensionalfew-bodysystemshavebeentreatednu- g, merically by various techniques [7, 11–17]. In this work, we apply standard Raleigh-Schr¨odinger V (z )=gδ(z z ), (1) 2b j0 j 0 − 2 wherez =z z . TheHamiltonianH fortheharmon- ically cojn0finedj−(N0,1) system then reads 10 (a) ) ω N N _h ( H = Hho(zj)+Hho(z0)+ V2b(zj0), (2) ) / 5 j=1 j=1 N X X ( E wherethesingleparticleharmonicoscillatorHamiltonian H (z) is given by ho 0 ~2 ∂2 1 2 H (z)= + mω2z2; (3) ho −2m∂z2 2 (b) ) N here, ω denotes the angular trapping frequency. The ( F delta-function interactions in Eq. (2) can be replaced by E 1 a set of boundary conditions on the many-body wave ) / N function Ψ(z ,z , ,z ), 0 1 N ( ··· ε ∂Ψ ∂Ψ gm = Ψ , (4) 0 ∂zj0(cid:12)zj0→0+ − ∂zj0(cid:12)zj0→0−! ~2 |zj0→0 -4 - _h ω 0a / g 4 (cid:12) (cid:12) ho where the(cid:12)(cid:12)limits zj0 (cid:12)(cid:12)0+, zj0 0−, and zj0 0 → → → are taken while keeping the other N coordinates, i.e., z , ,z ,z , ,z and (z +z )/2, fixed. 1 j−1 j+1 N j 0 ··· ··· Inthe following,we determine the eigenenergiesE(N) of the Hamiltonian H for various N. Throughout, we FIG. 1: (Color online) (a) Solid, dotted, and dashed lines restrict ourselves to the so-called upper branch. This show the energy of the upper branch for N = 1, 2, and 3 branch can be populated by preparing the system in the as a function of −1/g. The energies of the (1,1) system are non-interacting limit (g 0+) and by then adiabati- obtained by solving the transcendental equation derived in → cally first increasing g to large positive values, then con- Ref. [21]. The energies of the (2,1) and (3,1) system are tinuing acrossthe confinement-induced resonance[20] to taken from Refs. [14, 15, 22]. (b) Solid, dotted, and dashed infinitely negative g values and finally increasing g to lines show the interaction energy ǫ(N), normalized by the Fermi energy E (N), for N = 1, 2, and 3, respectively, as small negative values. Solid, dotted, and dashed lines F a function of −1/g. The harmonic oscillator length a is in Fig. 1(a) show the energy of the upper branch for ho definedin Eq. (10). N = 1 [21], 2, and 3 [14, 15, 22], respectively, as a func- tion of 1/g. For all N, the energy increases monotoni- − callyasafunctionofincreasing 1/g. Theupperbranch − forsystemswithN =1,2,and3majorityparticles. The corresponds to the ground state of the model Hamilto- energy E (N) is directly proportional to the number of nian when g is positive but not when g is negative. For F majority particles, negative g, the model Hamiltonian supports molecular- like bound states. In real cold atom systems, energet- E (N)=N~ω. (6) F ically lower-lying molecular states exist even for posi- tive g. However, it has been demonstrated experimen- Figure 1(b) shows that the normalized interaction en- tally [19] that the upper branch can be populated with ergy depends relatively weakly on the number of parti- reasonablyhighfidelity for positive g, motivating us—as cles. IndependentofN,wehaveǫ(N)=0forg =0+ and well as others [7, 12, 13, 16, 17, 23–26]—to investigate ǫ(N)=E (N)for g = . As canbe readoffFigs.1(a) F the properties of the upper branch within a stationary | | ∞ and 1(b), the energy increase of the upper branch is the zero-temperature quantum mechanics framework. Since sameonthepositiveg sideasitisonthenegativeg side, decay to states with molecular character can lead to sig- indicating that ǫ(N) approaches 2E (N) in the g = 0− F nificantdepopulationofthe upperbranchfornegativeg, limit for N = 1 3. We refer to E (N) as the Fermi F our primary focus in the following lies on the positive g − energy of the majority particles. It should be noted, portion of the upper branch. however, that the “exact” Fermi energy of the majority For g = 0+, the energy of the upper branch is equal particlesisE (N) ~ω/2,i.e.,E (N)correspondstothe to E (N) = (N2+1)~ω/2. We write the energy E(N) F − F ni leading order term of the Fermi energy of the majority of the upper branch in terms of the interaction energy particles in the large N limit. ǫ(N), One of the main goals of this paper is to derive ex- pansions for the interaction energy of the upper branch E(N)=E (N)+ǫ(N). (5) ni around g = 0+ and g = using standard Raleigh- | | ∞ Solid, dotted, and dashed lines in Fig. 1(b) show the Schr¨odinger perturbation theory for any N, i.e., for interaction energies, normalized by the energy E (N), N = 1, , . To this end, we express the interaction F ··· ∞ 3 energies ǫ(0+) and ǫ(|∞|) in the vicinity of g = 0+ and g = , respectively,inapowerseriesofthe dimension- TABLE I: Coefficients B(k)(N) for various (N,k) combina- | | ∞ tions. The numbers in parenthesis denote the uncertainty less interaction parameter γ (for g small) or in a power series of γ−1 (for g large) [27], | | thatarisesfromevaluatingtheperturbationtheorysumswith | | a finite energy cutoff. The numbers without errorbars have been rounded. kmax ǫ(0+)(N)= B(k)(N)γk E (N)+ (γkmax+1) (7) k=1 k=2 k=3 " # F O N=1 0.179587 −0.0223551 0.00179230 kX=1 N=2 0.190481 −0.0239838 0.00179523 N=3 0.194409 −0.0244852 0.00177603(1) and N=4 0.196423 −0.0247210 0.0017627(1) N=5 0.197647 −0.0248563 0.0017535(1) ǫ(|∞|)(N)= N=6 0.198469 −0.0249435 0.0017470(1) kmax N=7 0.199059 −0.0250042 1+ C(k)(N)γ−k EF(N)+ (γ−(kmax+1)), (8) N=8 0.199503 −0.0250488 " k=1 # O N=9 0.199849 −0.0250828 X N=10 0.200126 −0.0251096 wherethe dimensionlessinteractionparameterγ isgiven N=11 0.200353 −0.0251313 by N=12 0.200543 −0.0251491 N=∞ 0.202642 −0.0253303 0.00171100 π g γ = , (9) √2N ~ωaho with TABLE II: Coefficients C(k)(N) for various (N,k) combina- tions. Thenumbersinparenthesisdenotetheuncertaintythat ~ arisesfromevaluatingtheperturbationtheorysumswithafi- a = (10) niteenergycutoff. Thenumberswithouterrorbarshavebeen ho mω r rounded. k=1 k=2 k=3 denotingthe harmonicoscillatorlength. Aswewillshow N=1 −3.54491 3.85603 34.3007 below, the scaling of the interaction energy by E (N) F N=2 −3.17245 2.41904(1) 25.38(2) ensures a smooth connection between the energy shifts N=3 −3.02854 1.8142(2) 23.78(8) for finite and infinite N. In Eqs. (7)-(8), the dimension- N=4 −2.95040 less kth-order perturbation theory coefficients B(k)(N) N=5 −2.90081 andC(k)(N)dependonN andwillbe determinedinthe N=6 −2.86634 next section. N=7 −2.84091 N=8 −2.82133 N=9 −2.80578 III. PERTURBATIVE RESULTS N=∞ −2.66667 0 21.0552 N limit: The impurity problem for the homo- geneo→us ∞system with positive γ was solved by McGuire (1,1) system: The eigenenergies of the harmonically in 1965[28]. Within the localdensity approximationthe trapped (1,1) system can be obtained for any γ by solv- Fermi wave vector is replaced by the wave vector at the ing a simple transcendental equation [21]. Expanding trapcentersuchthattheinteractionenergyoftheground the transcendental equation around the known eigenen- state for the harmonically trapped system with N ergies for small and large γ, one obtains power series in becomes [19] →∞ theinteractionenergy. Invertingtheseseries,oneobtains analytical expressions for the B(k)(1) and C(k)(1) coeffi- ǫ( ) γ γ γ 2π γ cients. We find B(1)(1) = π−3/2, B(2)(1) = ln(2)/π3, E ∞( ) = π2 1− 4 + 2π + γ arctan 2π .(11) B(3)(1) = [π2 9ln(4)2]/(24π9/2), C(1)(1)−= 2π1/2, F ∞ (cid:20) (cid:18) (cid:19) (cid:16) (cid:17)(cid:21) C(2)(1) = −4π[ln−(2) 1], and C(3)(1) = π3/2[π2− 24 − − − − ExpandingEq.(11)aroundγ =0+ and γ = ,respec- 9(ln(4) 4)ln(4)]/3. The numericalvalues of these coef- tively, B(k)( ) and C(k)( ) can be ob|ta|ine∞d for k = ficients−are summarizedin Tables I-II. As in the N 1,2, . We∞findB(1)( )=∞2/π2, B(2)( )= 1/(4π2), case,thesmallandlargeγ seriesforN =1,Eqs.(7)→an∞d B(3)(···) = 1/(6π4), C∞(1)( ) = 8/3∞, C(2)(− ) = 0, (8), have a finite radius of convergence. Employing the and C∞(3)( ) = 32π2/15. T∞he num−erical values∞of these techniquesofRef.[29],wefind—usingupto50expansion coefficient∞s are summarized in Tables I-II. It is readily coefficients—that the small and large γ series converge shown that the small and large γ series, Eqs. (7) and for γ <1.0745(2) 2π and γ −1 <[1.0745(2) 2π]−1, | | × | | × (8), convergefor γ <2π and γ−1 <(2π)−1, respectively. respectively. Our result for the convergence of the small Table II shows that C(2)( ) vanishes. We will return γ series is consistent with what is reported in the litera- to this finding when we dis∞cuss the N-dependence of the ture [30]. C(k)(N) coefficients. Weakly-repulsive(N,1)system,N =2,3, : Totreat ··· 4 the weakly-interacting system with finite N, N > 1, we These matrix elements are closely related to the bound- rewrite the system Hamiltonian in second quantization aryconditionrepresentationofthe one-dimensionalodd- and expand the field operators for the majorityparticles parity pseudo-potential [35, 36]. We evaluate the inte- and the impurity in terms of single particle harmonic grals in Eq. (13) analytically for N =1 4. The analyt- − oscillator states (see, e.g., Ref. [31]). The interaction ical results for N = 1 and 2 read C(1)(1) = 2√π and − matrix elements can be evaluated analytically and the C(1)(2) = π/2(81/32). The analytical expressions − first-order perturbation theory treatment for positive g for N = 3 and 4 are lengthy and not reported here [37]. p yields ForlargerN,weperformallbutoneintegrationforeach oftheperturbationmatrixelementsanalytically. There- 2√NΓ(1/2+N) sultingnumericallydeterminedenergyshiftsareaccurate B(1)(N)= . (12) π2N! to morethan 10digits. Table II summarizesthe numeri- cal values for the coefficient C(1)(N) for N 9 obtained ≤ The first-order energy shift may be interpreted by us. The extension to larger N is, although tedious, as the leading-order mean-field shift. We find possible in principle. limN→∞B(1)(N) = 2/π2, which agrees with the coef- To determine the energy shift proportional to γ−2, ficient obtained by expanding Eq. (11). The evalua- we use second-order perturbation theory. Reference [38] tion of the second-order energy shift involves the eval- pointed out that the second-order perturbation theory uation of infinite sums. We find, as expected, that these energy shift of the (1,1) system diverges, thus requiring sums converge. The reason is that the one-dimensional regularization. Analogous divergencies arise in the per- delta-function interaction does not, unlike two- or three- turbativetreatmentofone-dimensionalsingle-component dimensional delta-function interactions [32, 33], require Fermi gases with generalized delta-function interactions any regularization if used in standard perturbation the- (see, e.g., Ref. [35]) and that of one-dimensional Bose oryapproaches. Wedidnotfindacompactanalyticalex- gases with effective range dependent zero-range interac- pression applicable to all N for the second-order energy tions. In the following, we discuss the impurity problem shift. For N = 1 and 2, we have B(2)(1) = ln(2)/π3 with N = 2 and 3. To evaluate the second-order energy − andB(2)(2)=[ 9+6√3+3ln(2+√3) 12ln(2)]/(4π3). shifts, we need to know the complete set of eigenstates − − For larger N, the expressions are lengthy. The numer- of the (2,1) and (3,1) systems with g = . For the ical values for N 12 are listed in Table I. Table I (2,1) system, we use the analytical wa|v|e fun∞ctions from ≤ also summarizes the numerically determined values for Ref. [39] and evaluate the integrals analytically. For the the third-order coefficients B(3)(N) for N = 2 6. The (3,1)system,wederivecompactformsfortheeigenstates B(3)(N) coefficient increases slightly as N chan−ges from using spherical coordinates and evaluate the relevant in- 1 to 2, and then decreases monotonically as N increases tegrals analytically. We then evaluate the second-order further. ThenumericallydeterminedB(3)(N)coefficients perturbation theory sums numerically, imposing an en- for N =2 6 approach the N = coefficient smoothly ergycutoff onthe relative energyof the intermediate (or − ∞ if plotted as a function of 1/N. virtual)states that are being summed over. The second- Strongly-interacting (N,1) system, N = 2,3, : The orderenergyshiftisfoundtocontainpowerlawdivergen- ··· strongly-interacting regime has been treated perturba- ciesintheenergycutoff. Thesedivergenciesarecanceled tively at leading order, i.e., at order 1/γ, for N 8 [26] through the introduction of a counterterm and the con- ≤ (note, though, that only the coefficients for N 4 were stant (and physically meaningful) part is extracted with ≤ reported explicitly, i.e., in equation or numerical form). high precisionby a regularizationscheme similar to that Toderivetheseresults,thetwo-bodyinteractionforlarge developedforharmonicallytrappedbosons[40]. TableII g ismodeledbyimposingthetwo-bodyboundarycondi- reports the resulting second-order perturbation theory | | tion on the many-body wave function when the distance coefficients with errorbars. Our perturbative coefficients betweenthe unlike particlesapproacheszero[23,26,34]. are consistent with the coefficients obtained by fitting Since the groundstate eigenenergy for g = is degen- the (2,1) and (3,1) energies reported in Refs. [15, 22] to | | ∞ erate,the perturbationshiftis obtainedbydiagonalizing a polynomial in γ−1. For the N = 2 and 3 systems, we the perturbation matrix constructed using the degener- extend the above treatment to the third order (see Ta- ate states for g = . For the many-body states Ψ and bleII).Thesethird-ordercalculationsrequiretheevalua- α ∞ Ψ , the perturbation matrix element V reads [23] tionofmatrixelementsV betweenexcitedstates. Since β αβ αβ the third-order perturbation expression is more involved ~4 thanthesecond-orderperturbationexpression,ourthird- V = αβ −m2g × order result has a larger errorbar than our second-order result [41]. N ∂Ψ∗ ∂Ψ∗ α α δ(z ) Thecalculationsofthe second-andthird-orderenergy j0 j=1Z ···Z ∂zj0(cid:12)zj0→0+ − ∂zj0(cid:12)zj0→0−! × shifts can, in principle, be extended to larger N. To do X (cid:12) (cid:12) so,twochallengesneedtobeovercome. First,anefficient (cid:12) (cid:12) ∂Ψβ (cid:12) ∂Ψβ (cid:12) dz dz dz (.13) method to generate the complete set of eigenstates at 0 1 N ∂zj0(cid:12)(cid:12)zj0→0+ − ∂zj0(cid:12)(cid:12)zj0→0−! ··· |g|=∞hastobedevised. Second,anefficientschemefor (cid:12) (cid:12) (cid:12) (cid:12) 5 TABLEIII:Fittingcoefficientsb(k)fork=2and3. Fork=2 -2.7 j and 3, we used j =6 and 4, respectively. max k=2 k=3 N) -3 j=0 −0.0253304 1.71100×10−3 (1)C( j=1 0.0019591 2.23905×10−4 -3.3 j=2 0.0033477 6.59881×10−6 (a) j=3 −0.0116972 −3.67025×10−4 -3.6 j=4 0.0379414 2.64073×10−4 4 j=5 −0.0684440 (b) j=6 0.0486346 3 ) N (2 2) evaluating the matrix elements andinfinite perturbation (C1 theory sums has to be developed. This is not pursued 0 here. 35 (c) IV. FITTING THE B(k)(N) AND C(k)(N) COEFFICIENTS )30 N ( 3) Tables I-II suggest that the coefficients B(k)(N) and (C25 C(k)(N) change,for fixedk, smoothly with N. This mo- tivates us to write 20 0 0.2 0.4 0.6 0.8 1 jmax 1 j 1/N B(k)(N)= b(k) (14) j N j=0 (cid:18) (cid:19) X and FIG. 2: (Color online) Symbols show the coefficients (a) C(k)(N)=jmaxc(k) 1 j. (15) C(1)(N), (b) C(2)(N), and (c) C(3)(N) as a function of 1/N. j N The solid lines show our fits with j =6, 3, and 3, respec- j=0 (cid:18) (cid:19) max X tively. It should be noted that the expressions (14)-(15) reduce tob(k) andc(k),respectively,inthe N limit. Inthe follo0wing, th0e parameters b(jk) and c(jk→) a∞re obtained by TkA=B1L,E2,IaVn:dF3i,ttwinegucsoedeffijcient=s6c,(jk3),faonrdk3,=re1s,p2e,catinvdely3.. For fitting the coefficients B(k)(N) and C(k)(N) for fixed k. max k=1 k=2 k=3 WestartwithB(2)(N). WefitEq.(14)totheB(2)(N) j=0 −2.66667 0.00000 21.05520 values for N = 1 80 (the values for N = 1 12 are − − j=1 −1.40749 7.06739 8.80915 reported in Table II), varying j from 2 20. We find max j=2 1.78704 −5.70589 −5.07455 − that the most reliable fit is obtained for jmax = 12 j=3 −4.21746 2.49453 9.51090 − 13. In this case, the fitting parameter b(2) differs from j=4 7.33094 0 1/(4π2) [the resultobtainedby expandingEq.(11)] by j=5 −7.13604 l−ess than 10−8. This suggests that not only the k = j=6 2.76476 1 coefficient (see discussion above) but also the k = 2 coefficient connects smoothly with the infinite N result. Table III reports the results of our fit to the B(2)(N) Symbols in Figs. 2(a)-2(c) show the C(k)(N) coeffi- coefficientswithN =1 80and byapolynomialwith cients with k =1,2, and 3, respectively, as a function of − ∞ j =6. 1/N. Our fits to these data (see Table II) using poly- max As mentioned earlier,the B(3)(1) coefficient is slightly nomials with j =6,3, and 3 are shown by solid lines max smaller than the B(3)(2) coefficient. The B(3)(N) coeffi- (see Table IV for the coefficients). cients for N 2, however, decrease monotonically. This The discussion so far has focused on the coefficients motivatesus≥tofittheB(3)(N)coefficientswithN =2 6 B(k)(N)andC(k)(N)withk=1 3. Itis,ingeneral,not and byapolynomialwithjmax =4. Thefitcoefficien−ts feasibletoextendthe perturbativ−ecalculationstohigher ∞ are reported in Table III. It can be seen that the coef- k for arbitrary N. However, for N = 1 and , the co- ficient b(3) agrees with the coefficient B(3)( ) reported efficients with larger k can be obtained readil∞y. We find in Table0I. We believe that our fit provides∞an accurate that B(k)(1) decreases monotonically with increasing k description of the 6<N < coefficients. (we c|hecked |this for k 50). The C(k)(1) coefficient ∞ ≤ | | 6 increases monotonically with increasing k for k <37; for (a) k 37,weobservesmallnon-monotonicoscillations. For 1.1 ≥ N = ,wefindthattheB(k)( )withk evenandk 4 vanish∞while the B(k)( ) wi∞th k odd decrease mo≥no- 2) | ∞ | ρ( tonically with increasing k (again, we checked this for k 50). Similarly, the C(k)( ) with k even and k 2 1 ≤ ∞ ≥ vanishwhilethe C(k)( ) withkoddincreasemonoton- | ∞ | ically with increasing k. Assuming a linear change with 1/N,interpolatingbetweenB(k)(1)andB(k)( )andbe- tweenC(k)(1)andC(k)( )fork >3yieldses∞timatesfor 1.1 (b) ∞ the finite N, N > 1, coefficients. While rough, these es- ) timatesmightprovideareasonablemeanstoconnectthe 3 ( ρ weakandstrongperturbationtheorylimitsforquantities such as those shown in Figs. 3 and 4. 1 We cannot accurately estimate the radius of conver- genceofthesmallandlargeγexpansionsfor1<N < . ∞ However,the fact thatthe radiusofconvergenceis given (c) by γ <1.0745(2) 2πforN =1andγ <2π forN = 1.1 for|th|e small γ ser×ies and by γ −1 < [1.0745(2) 2π]−∞1 for N = 1 and γ−1 < (2π)−1| f|or N = for th×e large ∞) ∞ ρ( γ series suggests two speculations: First, a convergent 1 series can be found for any γ and N. Second, the radius of convergence of the small γ series is approximately 2π for all N. Figures 3 and 4, which are discussed in the -1 -0.5 0 0.5 -1/γ next section, are consistent with these speculations. V. DISCUSSION FIG. 3: (Color online) Solid lines show the quantity ρ(N) This sectioncomparesthe perturbativeenergyexpres- as a function of −1/γ for (a) N = 2, (b) N = 3, and (c) sionswiththenumericallydeterminedenergiesoftheup- N = ∞, respectively. For comparison, dotted, dashed, and perbranch. Figure1(b)showsthatthescaledinteraction dash-dotted lines show the perturbative results for ρ(N) ac- energy ǫ(N)/EF(N) depends weakly on N if plotted as counting for terms up to order γ0, γ1 and γ2, respectively, a function of ~ωaho/g. The dependence on N is even in the weakly-interacting regime and accounting for terms weakerwhent−heinteractionstrengthisparameterizedby up to order γ−1, γ−2 and γ−3, respectively, in the strongly- γ asopposedtog. Tobenchmarktheapplicabilityofthe interacting regime. perturbative expressions we analyze the interaction en- ergyofthe systemwithN majorityatomsby comparing where the coefficients ρ(w)(N) are determined by the with that of the (1,1) system. Specifically, we consider k the quantities ρ(N), B(l)(N) and B(l)(1) with l k+1. Expanding Eq. (16) ≤ in the strongly-interacting (large γ ) regime, we find | | ǫ(N)/E (N) ρ(N)= F , (16) ρ(N)=1+ρ(s)(N)γ−1+ρ(s)(N)γ−2+ ǫ(1)/E (1) 1 2 F ρ(s)(N)γ−3+ (γ−4), (19) and δ(N,N′), 3 O where the coefficients ρ(s)(N) are determined by the k ǫ(N)/E (N) ǫ(1)/E (1) C(l)(N) and C(l)(1) with l k. Solid lines in Figs. 3(a)- δ(N,N′)= F − F . (17) ≤ ǫ(N′)/E (N′) ǫ(1)/E (1) 3(c) show the quantity ρ(N) for N = 2, 3, and , re- F − F spectively. The solid lines are obtained using th∞e nu- For finite N, ρ(N) reduces to ǫ(N)/[Nǫ(1)], i.e., ρ(N) merical (2,1) and (3,1) energies and the semi-analytical tells one the interaction energy per particle, normalized (1,1) and ( ,1) energies. For γ 0+ (i.e., for ∞ → bythe interactionenergyofthe(1,1)system. Thequan- 1/γ ), the quantity ρ(N) approaches the con- − → −∞ tity δ(N,N′) can alternatively be written as [ρ(N) stantρ(w)(N)=B(1)(N)/B(1)(1)[seethehorizontaldot- 1]/[ρ(N′) 1]. − ted line0s in Figs. 3(a)-3(c)], which increases monotoni- − ExpandingEq.(16)intheweakly-interacting(smallγ) callyfrom1.0607to 1.1284asN goes from2 to . This ∞ regime, we find portionoftheinteractionenergycanbeinterpretedasthe mean-field contribution. Inclusion of the next order cor- ρ(N)=ρ(w)(N)+ρ(w)(N)γ+ρ(w)(N)γ2+ (γ3),(18) rection[theρ(w)(N)γ term]andthe nexttwocorrections 0 1 2 O 1 7 [the ρ(w)(N)γ and ρ(w)(N)γ2 terms] yields the dashed 0.48 1 2 anddash-dottedlinesinFigs.3(a)-3(c). Thedash-dotted linesprovideafairlyaccuratedescriptionofthe quantity ) ∞ ρ(N) for −1/γ . −0.4. For |γ| → ∞, the leading-order 2, γ-dependent term [see the (non-horizontal) dotted lines δ(0.44 inFigs.3(a)-3(c)]increasesmonotonicallyfrom0.3725to (a) 0.8782as N changes from2 to . Inclusion of the next- ∞ order correction and the next two corrections yields the dashed and dash-dotted lines in Figs. 3(a)-3(c). It can 0.64 be seen that the dash-dotted lines provide a fairly accu- rate description of the quantity ρ(N) for 1/γ & 0.15. ) − − ∞ This value is close to the expected radius of convergence 3, oftheinteractionenergy[recall,theradiusofconvergence δ( 0.6 is 1/γ = (1.0745 2π)−1 0.148 for the (1,1) system]. (b) × ≈ Combining the perturbative descriptions for small and large γ , the expansions provide a fairly accurate de- | | scription of the interaction energy for the system with N majority particles, normalized by that for the (1,1) system, over a wide range of interaction strengths γ. )0.73 3 Expanding Eq. (17) in the weakly-interacting regime, 2, ( we find δ (c) δ(N,N′)=δ(w)(N,N′)+δ(w)(N,N′)γ+ 0.72 0 1 δ2(w)(N,N′)γ2+O(γ3), (20) -1.5 -1 -1/γ -0.5 0 where the coefficients δ(w)(N,N′) are determined by the k B(l)(N),B(l)(N′),andB(l)(1)withl k+1. Expanding ≤ Eq. (17) in the strongly-interacting regime, we find ′ δ(N,N′)=δ0(s)(N,N′)+δ1(s)(N,N′)γ−1+ F−I1G/γ..4:T(hCeolsoorliodnlliinnee)isThfoerq(uaa)nNtity=δ2(Na,nNd N)a′s=a∞funacntidonthoef δ(s)(N,N′)γ−2+ (γ−3), (21) circlesarefor(b)N =3andN′ =∞and(c)N =2andN′ = 2 O 3. In the large γ regime, the uncertainty of the numerically where the coefficients δ(s)(N,N′) are determined by the determined (3,1) energies leads to appreciable uncertainties k inδ(2,3)andδ(3,∞)[seetheerrorbarsinFigs.4(b)and4(c)]. C(l)(N), C(l)(N′),andC(l)(1)withl k+1. Thequan- ≤ For comparison, dotted, dashed, and dash-dotted lines show tity δ(N,N′) is shown by the solid line in Fig. 4(a) for ′ theperturbativeresultsforδ(N,N )accountingfortermsup (N,N′)=(2, )andbydotsinFigs.4(b)-4(c)for(3, ) toorderγ0,γ1,andγ2,respectively,intheweakly-interacting, ∞ ∞ and (2,3), respectively. We observe that the quantity small |γ| regime and accounting for terms up to order γ−1, δ(N,N′) changes only slightly as 1/γ goes from γ−2, and γ−3, respectively, in the strongly-interacting, large − −∞ to 0; this is particularly true for δ(2,3) [see Fig. 4(c)]. |γ| regime. The limiting values [see the dotted lines in Figs. 4(a)- 4(c)] are given by δ(w)(N,N′) and δ(s)(N,N′), respec- 0 0 tively. Dashed lines include the next order correction in thelargeγ regimetomeaningfullycomparewiththeper- the weakly- and strongly-interacting regimes, and dash- turbative results. dotted lines include the next two corrections. In the weakly-interacting regime, the dash-dotted lines provide afairlygooddescriptionofthe quantityδ(N,N′). Inthe VI. CONCLUSION strongly-interactingregime,however,the validity regime oftheperturbativeexpressionsisquitesmall. Forδ(2,3), This paper considered the upper branch of a non- e.g., the expansion coefficients are δ(s)(2,3) = 0.7213, 0 interacting harmonically trapped one-dimensional Fermi δ1(s)(2,3) = 0.0694(3), and δ2(s)(2,3) = −2.31(12), where gas with a single impurity. Zero-rangetwo-body contact the numbers in brackets denote the errorbars due to the interactions with strength g were assumed between the uncertaintiesofthesecond-andthird-orderperturbation majority atoms and the impurity. This system consti- theorycoefficients. Thefactthat δ(s)(2,3) δ(s)(2,3) tutes one of the simplest mesoscopic systems accessible | 2 |≫| 1 | isresponsibleforthe turn-aroundofthe dash-dottedline to experiment and theory. On the experimental side, it for large positive γ. We note that the errorbar of the hasbeendemonstratedbythe Heidelberggroupthatthe quantities δ(3, ) and δ(2,3), obtained from the numer- upperbranchofthemodelHamiltoniancanbeemulated ∞ ical energies [see dots in Figs. 4(b)-4(c)], is too large in reliably using ultracold atoms [18, 19]. On the theory 8 side, various numerical and analytical techniques have a program can be carried through explicitly beyond the beenapplied[7,11–17,23–26]. Thispaperpursuedaper- leading-order correction. (iii) The radii of convergence turbative approach, which determined expansions of the of the series were reported for N = 1 and . (iv) The energy of the upper branch in the weakly- and strongly- behavior of the expansion coefficients in the∞series in γk interacting regimes for various N. In the cases where and γ−k with k > 3 was discussed. (v) The perturba- we were not able to obtain general N expressions for a tive expressionswerebenchmarkedandfound to provide fixed order in the perturbative expansion, approximate a reliable description over a wide range of interaction expressions applicable to all N were obtained through strengths. fits. Through comparison with accurate numerical few- The results presented in this work can be used to cal- body energies, the perturbative expressions were shown culateperturbativeexpressionsforthecontactandother to provide a satisfactory description for a wide range of observables. Moreover, the second- and third-order re- interaction strengths. sultsintheγ−1seriesallowonetoassesstheapplicability The main results of this work are: (i) We determined regime of effective spin models [26, 34, 42]. an expansion for the energy of the upper branch of a one-dimensional harmonically trapped Fermi gas with a Acknowledgement: We thank G. Zu¨rn for email cor- single impurity in the weakly-repulsive regime up to or- respondence that motivated this work and suggested to der γ3, applicable to any system size. (ii) We deter- analyze the quantity δ(N,N′) shown in Fig. 4. We also mined an expansion for the energy of the upper branch thank N. Zinner for bringing Ref. [38] to our attention, in the strongly-interacting regime up to order γ−3, ap- and P. Johnson and E. Tiesinga for helpful discussions plicable to any system size. While the idea to treat the on the numerical regularization of diverging sums. Sup- coupling strength 1/γ as a small parameter is not new, port by the National Science Foundation through grant our work provides an explicit demonstration that such number PHY-1205443is gratefully acknowledged. [1] D. C. Mattis, The Many-Body Problem. An Encyclope- [19] A.N.Wenz,G.Zu¨rn,S.Murmann,I.Brouzos,T.Lompe, dia of Exactly Solved Models in One Dimension,World and S.Jochim, Science 342, 457 (2013). Scientific,1993. [20] M. Olshanii, Phys. Rev.Lett. 81, 938 (1998). [2] T.Giamarchi,Quantum Physics inOne Dimension,Ox- [21] T.Busch,B.-G.Englert,K.Rza¸z˙ewski,andM.Wilkens, ford UniversityPress, Oxford, 2004. Found.Phys. 28, 549 (1998). [3] A. Imambekov, T. L. Schmidt, and L. I. Glazman, Rev. [22] In the vicinity of |g| → ∞, we improved the numerical Mod. Phys. 84, 1253 (2012). energies reported in Ref. [15]. [4] X.-W. Guan, M. T. Batchelor, and C. Lee, Rev. Mod. [23] A. G. Volosniev, D. V. Fedorov, A. S. Jensen, M. Va- Phys.85, 1633 (2013). liente, and N. T. Zinner, Nat.Commun. 5, 5300 (2014). [5] Z. N. C. Ha. Quantum Many-Body Systems in One Di- [24] E.J.Lindgren,J.Rotureau,C.Forss´en,A.G.Volosniev, mension, Series on Advances in Statistical Mechanics - and N.T. Zinner,New J. Phys.16, 063003 (2014). Volume12, World Scientific, 1996. [25] X. Cui and T.-L. Ho, Phys.Rev.A 89, 023611 (2014). [6] C. Menotti and S. Stringari, Phys. Rev. A 66, 043610 [26] J. Levinsen, P. Massignan, G. M. Bruun, and M. M. (2002). Parish, arXiv:1408.7096. [7] G. E. Astrakharchik and I. Brouzos, Phys. Rev. A 88, [27] Throughout, we refer to the expansion around γ = 0+ 021602(R) (2013). and |γ−1| = 0 as small and large γ series. Sections III [8] G. E. Astrakharchik, D. Blume, S. Giorgini, and L. P. and IV discuss the applicability regime, i.e., the radius Pitaevskii, Phys.Rev.Lett. 93, 050402 (2004). of convergence, of these series. [9] I.V. Tokatly,Phys. Rev.Lett. 93, 090405 (2004). [28] J. B. McGuire, J. Math. Phys.6, 432 (1965). [10] G. Orso, Phys. Rev.Lett. 98, 070402 (2007). [29] J. Zamastil and F. Vinette, J. Phys. A: Math. Gen. 38, [11] M.Casula,D.M.Ceperley,andE.J.Mueller,Phys.Rev. 4009 (2005). A 78, 033607 (2008). [30] S. Kvaal, E. Jarlebring, and W. Michiels, Phys. Rev. A [12] I. Brouzos and P. Schmelcher, Phys. Rev. A 87, 023605 83,032505 (2011).Convertedtoourunits,thisreference (2013). reports that the small γ series converges for γ < 2π. [13] N.L. Harshman, Phys. Rev.A 86, 052122 (2012). InspectionoftheinsetofFig.2ofthisreferencesuggests, [14] S. E. Gharashi, K. M. Daily, and D. Blume, Phys. Rev. however, thattheradiusofconvergenceisslightly larger A 86, 042702 (2012). than 2π. [15] S. E. Gharashi and D. Blume, Phys. Rev. Lett. 111, [31] A.J.Leggett,QuantumLiquids: BoseCondensationand 045302 (2013). Cooper Pairing in Condensed-Matter Systems, Oxford [16] P. O. Bugnion and G. J. Conduit, Phys. Rev. A 87, University Press Inc., NewYork,2006. 060502(R) (2013). [32] K. Huangand C. N. Yang,Phys. Rev.105, 767 (1957). [17] T. Sowin´ski, T. Grass, O. Dutta, and M. Lewenstein, [33] K. W´odkiewicz, Phys. Rev.A 43, 68 (1991). Phys.Rev.A 88, 033607 (2013). [34] F. Deuretzbacher,D. Becker, J. Bjerlin, S.M. Reimann, [18] F. Serwane, G. Zu¨rn, T. Lompe, T. B. Ottenstein, and L. Santos. Phys.Rev.A 90, 013611 (2014). A.N. Wenz,and S. Jochim, Science 332, 336 (2011). [35] M. D. Girardeau and M. Olshanii, Phys. Rev. A 70, 9 023608 (2004). state with smaller slope. Lastly, the perturbation analy- [36] T. Cheon and T. Shigehara, Phys. Rev. Lett. 82, 2536 sis allows one to improve the accuracy of the coefficient (1999). c2 = 0.865(7), which is reported below Eq. (10) of the [37] Our C(1)(3) coefficient agrees with that reported in supplementalmaterialofRef.[15];theperturbativeanal- Ref. [26] but differs by 0.0012 % from the slope that ysis yields c2 ≈0.8626961. one obtains using the adiabatic wave function reported [38] D. Sen,Int.J. Mod. Phys. A 14 1789 (1999). in Ref. [15] [see Eq. (2) of the supplemental material]. [39] A.G. Volosniev, D.V.Fedorov,A.S. Jensen,and N.T. This implies (see also Ref. [26]) that the true adiabatic Zinner, arXiv:1408.6058. wave functions are linear combinations of Eqs. (2) and [40] P. R. Johnson, D. Blume, X. Y. Yin, W. F. Flynn, and (3) of the supplemental material of Ref. [15] with coeffi- E. Tiesinga, NewJ. Phys.14, 053037 (2012). cients α≈0.9999929 and −(1−α2)1/2 ≈−0.003764883 [41] For theN =2 system, weused energy cutoffs of 1920~ω for state 1 and coefficients (1−α2)1/2 and α for state and 708~ω for the second- and third-order perturbation theory calculations. For the N = 3 system, we used en- 2. An analogous analysis shows that the true adiabatic ergycutoffsof190~ω and80~ωforthesecond-andthird- wave functions of the (2,2) system are linear combina- order perturbation theory calculations. tions of Eqs. (11) and (12) of the supplemental material of Ref. [15] with coefficients −(1−α2)1/2 and α for the [42] A. G. Volosniev, D. Petrosyan, M. Valiente, D. V. Fe- dorov, A. S.Jensen, and N. T. Zinner. arXiv:1408.3414. state with larger slope and α and (1−α2)1/2 for the