A new effective interaction for the trapped fermi gas: the BEC-BCS crossover

# A new effective interaction for the trapped fermi gas: the BEC-BCS crossover

C. N. Gilbreth and Y. Alhassid Center for Theoretical Physics, Sloane Physics Laboratory, Yale University, New Haven, CT 06520
###### Abstract

We extend a recently introduced separable interaction for the unitary trapped Fermi gas to all values of the scattering length. We derive closed expressions for the interaction matrix elements and the two-particle eigenvectors and analytically demonstrate the convergence of this interaction to the zero-range two-body pseudopotential for -wave scattering. We apply this effective interaction to the three- and four-particle systems along the BEC-BCS crossover, and find that their low-lying energies exhibit convergence in the regularization parameter that is much faster than for the conventional renormalized contact interaction. We find similar convergence properties of the three-particle free energy at unitarity.

###### pacs:
03.75.Ss, 05.30.Fk, 21.60.Cs, 31.15.-p

## I Introduction

Atomic Fermi gas systems have generated much interest in the last several years giorgini:1215 (); bloch:885 (). While these systems can be well-described by a relatively simple model Hamiltonian, at low temperatures they exhibit a rich phenomenology whose theoretical description has proven to be challenging. At low momenta, or when the range of the interaction is sufficiently short, the dominant scattering process occurs in the -wave channel, so that a single parameter, the -wave scattering length , suffices to characterize the inter-particle interaction. Depending on the value of (where is the Fermi momentum), widely different behavior is observed. This ranges from the formation of tightly bound dimers in the Bose-Einstein condensate (BEC) regime when is small and positive, to Cooper pairing in the Bardeen-Cooper-Schrieffer (BCS) regime when is small and negative. These behaviors are connected by a strongly interacting, nonperturbative regime, known as the unitary regime, where is much larger than any other length scale in the system. Remarkably, each of these regimes is accessible experimentally, and all exhibit superfluid behavior below a certain -dependent critical temperature.

While accurate mean-field theories exist for the BEC and BCS regimes, there is no simple approximation that can accurately describe the transition between these regimes through the unitary limit. As a result much interest has been taken in applying numerical methods, such as quantum Monte Carlo and numerical diagonalization. For systems varying from to hundreds of particles, quantum Monte Carlo calculations have been carried out to calculate ground state carlson:2003 (); astrakharchik:2004 (); chang:2004 (); chang:2005 (); stecher:2007 (); blume:2007 (); jauregui:2007 (); greene:2008 (); gezerlis:2008 (); zinner:2009 () and thermodynamic properties bulgac:2006-9 (); burovski:2006 (); chang:2007 (); akkineni:2007 (); goulko:2010 (). For smaller numbers of trapped atoms, energy spectra have been calculated using a basis set expansion method with correlated Gaussians in coordinate space (for up to six atoms) stecher:2007 (); blume:2007 () and a stochastic variational approach (for four atoms) daily:2010 ().

The three-body problem has been solved analytically in the unitary regime werner:2006 (); tan3body () and numerically to high accuracy along the BEC-BCS crossover Kestner:2007 (); Liu:2009 (). In fact the few-body trapped cold atom problem has become more interesting recently, since it was pointed out that the virial expansion for the partition function in a harmonic trap works well at unexpectedly low temperatures into the quantum degenerate regime Liu:2009 (). Moreover, the scaling relations which exist for the dilute gas can allow larger systems to be addressed from the study of smaller ones adhikari:2009 (); zhang:2009 ().

Here we discuss the cold atom problem in the context of the configuration-interaction (CI) approach to interacting many-particle fermionic systems. In this approach, widely used in atomic, molecular and nuclear physics, a single-particle basis in the laboratory frame is used to construct a many-particle basis of Slater determinants for a fixed number of fermions. The Hamiltonian matrix is then calculated in this many-particle basis and diagonalized. For a harmonically trapped system, a natural choice for the single-particle basis is that of the three-dimensional harmonic oscillator.

The short-range interaction of the cold atom problem is often approximated by a zero-range interaction. However, the obvious form of such a potential, a contact interaction , is ill-defined in three dimensions and must be regularized. This is usually accomplished by introducing a momentum or energy cutoff in relative motion and renormalizing the strength of the interaction so it reproduces the physical scattering length for the uniform gas, or the lowest bound state energy in relative motion for the trapped system. The many-particle energies are then calculated as a function of the regularization cutoff parameter.

A contact interaction is non-vanishing only for relative angular momenta of the two particles. For a basis of harmonic oscillator wave functions, a natural regularization parameter is the number of oscillator waves in relative motion. However, the convergence of the many-particle energies as a function of this regularization parameter for the renormalized contact interaction is slow (as a low negative power) Stetcu:2008 (); Alhassid:2008 ().

In Ref. Alhassid:2008, a new effective interaction was introduced in the unitary limit for which the low-lying energies of the three- and four-particle systems converge substantially faster than for the conventional renormalized contact interaction as a function of the regularization parameter. While this effective interaction is no longer a contact interaction, it reproduces the same many-particle energies in the limit when the regularization parameter is sent to infinity. The faster convergence of the many-particle energies enables their calculation to higher accuracy when using this new effective interaction. Another approach to improve the convergence relative to the renormalized contact interaction was investigated within the framework of effective field theory by including perturbatively next-to-leading order and next-to-next-to-leading order interactions Stetcu:2010 ().

Here we generalize the construction of such an effective interaction away from unitarity. We study the low-energy spectra of the three- and four-particle systems and find convergence properties versus the regularization parameter along the complete BEC-BCS crossover that are similar to the unitary regime. We also use this effective interaction to study the thermodynamics of the three-particle system at unitarity. In particular, we demonstrate the exponential convergence of the free energy at finite temperature as a function of the regularization parameter. The converged free energy is compared with the exact free energy constructed from the known analytical spectrum of the three-particle system werner:2006 (); tan3body ().

This paper is organized as follows. In Sec. \@slowromancapii@, we briefly review the two-particle cold atom problem in a harmonic trap. In Sec. \@slowromancapiii@, we derive closed expressions for the matrix elements of the new effective interaction for any scattering length and demonstrate convergence (as a function of the regularization parameter) of its two-particle eigentates to the exact eigenstates in the unitary limit. In Sec. \@slowromancapiv@, we use this interaction to study the convergence (vs. the regularization parameter) of the spectra of the three- and four-body systems along the BEC-BCS crossover and the three-particle free energy at unitarity. The latter is compared with exact results. Finally, in Sec. \@slowromancapv@ we present our conclusion.

## Ii Two-particle problem

The trapped two species cold atom system is modeled by the Hamiltonian

 H=−N∑i=1ℏ22m∇2i+N∑i=112mω2r2i+∑i

where , and are the number of atoms for each species (spin-up and spin-down), is the trap frequency, and is a short-range interaction. A natural length scale for this Hamiltonian is the oscillator length .

The interaction is modeled by a zero-range, pure -wave interaction, which due to Pauli exclusion acts only between particles of differing species. This can be expressed by a regularized -function,

 V(r)=gδ3(r)(∂/∂r)r, (2)

where relates the scattering length to the strength of the interaction giorgini:1215 (). Equivalently, one can impose the Bethe-Peierls contact conditions on the wavefunctions; see, e.g., in Ref. werner:2006, . Several interesting consequences of the form (2) of the interaction potential for the many-body system have been explored in Ref. tan:2008, .

For a two-particle system, the Hamiltonian is separable in the center-of-mass and relative coordinates and , so that , with

 Hrel=−ℏ22μ∇2r+12μω2r2+V(r)≡H0+V(r),

where is the harmonic oscillator Hamiltonian in relative motion with reduced mass .

The eigenstates of the non-interacting two-particle system in a harmonic trap may be labeled as , where are the radial, angular momentum, and magnetic quantum numbers for the center of mass motion, and are the corresponding quantum numbers for the relative motion. The associated energies are . A pure -wave interaction leaves the states and energies unchanged while mixing the states into eigenstates we denote by , with . The energies in relative motion for scattering length are given by the solutions to the transcendental equation Busch:1998 (); Jonsell:2002 ()

 Γ(−ε/(2ℏω)+3/4)Γ(−ε/(2ℏω)+1/4)=1√2a/aosc. (3)

The exact wave functions in relative motion are also known. They are given by Busch:1998 ()

 |u(i)⟩=∞∑n=0u(i)n|n00⟩,

where

 u(i)n=Aiφ∗n(0)αn−εi.

Here are the non-interacting relative energies, is a normalization factor, and is a harmonic oscillator wave function, with

 φn(0)=π−3/4[(2n+1)!!(2n)!!]1/2.

In the unitary limit, where , the normalization factor has the simple form

 A−2i = π−3/2∞∑n=0(2n+1)!!(2n)!!1[2(n−i)+1]2 = π−1/22(2i)!!(2i−1)!!.

## Iii Effective interaction

We first discuss the conventional contact interaction. It can be regularized by introducing an energy cutoff in relative motion. This defines a sequence of interactions , , with

 ⟨nlm|^V(q)c|n′lm⟩=χqψn(0)ψn′(0)δl,0δm,0, (4)

for , and where . All other matrix elements of are taken to be zero. The coefficient is determined for each so as to reproduce the lowest relative-motion energy . We refer to this interaction as a renormalized contact interaction.

The contact interaction in Eq. (4) is separable in the relative oscillator basis. A general separable interaction for has the form

 ⟨nlm|^V(q)eff|n′lm⟩=f∗nfn′δl,0δm,0(n,n′⩽q). (5)

An effective interaction for the trapped cold atom system can be determined by choosing the coefficients to reproduce the lowest relative energies of the two-particle system. This idea was used in Ref. Alhassid:2008, to construct a new effective interaction in the unitary limit of infinite scattering length. Here we consider the more general problem for any value of the scattering length .

### iii.1 General problem and solution

The matrix elements of the relative-coordinate Hamiltonian take the form

 ⟨n00|^H(q)rel|n′00⟩=δn,n′αn−f∗nfn′, (6)

where are the non-interacting eigenvalues of . To determine the coefficients we derive the following result.

Theorem. Let and be real numbers such that , and for all , and let be complex. Then the -dimensional matrix has the as its eigenvalues if and only if

 |fn|= ⎷∏k(αn−εk)∏k≠n(αn−αk), (7)

in which case the th eigenvector has components

 b(i)n=Cif∗nαn−εi,Ci= ⎷∏k(αk−εi)∏k≠i(εk−εi). (8)

A solution satisfying (7) exists if and only if .

Proof. An eigenvector with eigenvalue must satisfy

 q∑n=0(δmnαn−f∗mfn)b(i)n=εib(i)m (9)

which implies

 b(i)n=Cif∗nαn−εi (10)

where is a normalization constant. Inserting (10) into (9), we obtain

 q∑n=0Min|fn|2=1,

where the matrix is defined by

 Min=1αn−εi. (11)

Thus the are determined from the eigenvalues by the matrix equation , where and .

The Matrix is of a well-known form, called a Cauchy matrix, whose inverse is given by knuthv1 ()

 (M−1)ni=∏k(αn−εk)(αk−εi)(αn−εi)∏k≠i(εk−εi)∏k≠n(αn−αk).

We therefore have

 |fn|2 = (M−11)n (12) = ∑i∏k(αn−εk)(αk−εi)(αn−εi)∏k≠i(εk−εi)∏k≠n(αn−αk) = ∏k(αn−εk)∏k≠n(αn−αk)∑i∏k≠n(εi−αk)∏k≠i(εi−εk).

The sum in the last expression can be evaluated using the identity knuthv1 ()

 (13)

where are all distinct. This is done by observing that , where is a polynomial of degree smaller than . There are of the , which we take to be in Eqs. (13). According to the first case in Eqs. (13), the contribution to the last sum in Eq. (12) from vanishes, while according to the second case of (13), the contribution of to the sum is exactly 1. Thus, the components satisfy (7). On the other hand, it is easy to verify that if the satisfy (7) then given by (10) is an eigenvector with eigenvalue .

If then it can easily be seen that the argument to the square root in Eq. (7) is positive, and therefore a solution for the exists. Conversely, we can show that if the argument to the square root is positive for all , the eigenvalues must have this ordering. For each let be the number of eigenvalues for which . Then

 sign∏qk=0(αn−εk)∏k≠n(αn−αk) = (−1)q+1−kn(−1)q−n.

In order for the sign to be positive for all , we therefore require that is always even. But since , the only possibility is , i.e., .

The eigenvectors are given by Eq. (10). To determine the normalization constant , we require and use Eq. (7) for to find

 C−2i = q∑n=0∏0⩽k⩽q(αn−εk)∏k≠n,0⩽k⩽q(αn−αk)1(αn−εi)2.

This sum can also evaluated using the identity (13). To do this, we define and rewrite

 C−2i = q∑n=0∏k≠i,0⩽k⩽q(αn−εk)∏k≠n,0⩽k⩽q+1(αn−αk).

We now add and subtract the term to obtain

 C−2i = q+1∑n=0∏k≠i,0⩽k⩽q(αn−εk)∏k≠n,0⩽k⩽q+1(αn−αk) −∏k≠i,0⩽k⩽q(αq+1−εk)∏0⩽k⩽q(αq+1−αk).

Taking , the first case of (13) may be applied to find that the sum vanishes and thus

 C−2i = ∏k≠i(εk−εi)∏k(αk−εi).

This completes the proof.

Since the eigenvalues , for the cold atom problem are indeed ordered as the theorem requires, we choose to be positive real numbers according to Eq. (7).

The complete set of the two-particle eigenvectors of the effective interaction is therefore given by (i) the noninteracting states for or ( and ) with eigenvalues , and (ii) the interacting states

 |NLMb(i)⟩=[∏k(αk−εi)∏k≠i(εk−εi)]1/2×q∑n=0fnαn−εi|NLMn00⟩ (14)

with eigenvalues and where the are given by the solutions of Eq. (3).

We also derive closed expressions for the trace and norm of the -dimensional matrix defined by . Its trace has a simple form which can be obtained via Eqs. (7) and (13). Taking and applying the second and third cases of Eq. (13), we find

 tr V(q)eff = q∑n=0∏k(αn−εk)∏k≠n(αn−αk) = q∑n=0(αn−εn).

This result implies that the Frobenius norm of is

 ∑n,n′|(V(q)eff)nn′|2=(q∑n=0(αn−εn))2.

### iii.2 Properties in the unitary limit

In this section we consider the properties of the separable effective interaction in the unitary limit. The simple form of the relative-motion unitary eigenvalues allows us to verify Eq. (9) in Ref. Alhassid:2008, for the interaction parameters at unitarity and prove analytically that converges to the pseudopotential (2) in the limit .

#### iii.2.1 Interaction parameters fn

In the unitary limit are the harmonic oscillator energies shifted down by one oscillator quantum. In this case the expression for reduces to

 f2n = ∏qk=0(2(n−k)+1)∏k≠n2(n−k) (15) = ∏nk=0(2k+1)∏nk=12k∏q−nk=1(−2k+1)∏q−nk=1(−2k) = (2n+1)!!(2n)!!(2(q−n)−1)!!(2(q−n))!!,

which is exactly Eq. (9) of Ref. Alhassid:2008, .

#### iii.2.2 Convergence of two-particle eigenvectors

To prove the convergence of the effective interaction to the pseudopotential in Eq. (2), we study the convergence of the two-particle eigenstates of to the corresponding eigenstates of as . For any scattering length, the square of the -th component of is given by

 (b(i)n)2 = ∏k(αk−εi)∏k≠i(εk−εi)∏k(αn−εk)∏k≠n(αn−αk)1(αn−εi)2.

In the unitary case

 (b(i)n)2 = ∏k(2(k−i)+1)∏k≠i2(k−i)∏k(2(n−k)+1)∏k≠n2(n−k)1(αn−εi)2 = 1(αn−εi)2(2i−1)!!(2i)!!(2(q−i)+1)!!(2(q−i))!! ×(2n+1)!!(2n)!!(2(q−n)−1)!!(2(q−n))!!.

The asymptotic behavior for large is

 b(i)n(q)∼ √2/π(αn−εi)[(2i−1)!!(2i)!!]1/2[(2n+1)!!(2n)!!]1/2 ×[q−i+3/4q−n+1/2]1/4.

We compare to the exact wave function by computing the relative difference for each component,

 |b(i)n−u(i)n|u(i)n ∼ ∣∣ ∣∣(1+n−i+1/4q−n+1/2)1/4−1∣∣ ∣∣ ∼ n−i+1/44(q−n+1/2).

This diminishes like for large , proving convergence to the exact two-body eigenstates in relative motion. Importantly, this is faster than the convergence of the renormalized contact interaction Alhassid:2008 ().

## Iv Applications

To demonstrate the advantage of the new effective interaction as compared with the renormalized contact interaction, we study the convergence versus of low-lying energies of the three and four-particles systems at various values of the scattering length. We also study the convergence versus of the three-particle free energy at unitarity (for which an exact solution exists). The latter is important for thermodynamical studies of the trapped gas. The three-particle system we study consists of two spin-up particles and one spin-down particle (), while the four-particle system is unpolarized ().

### iv.1 Spectroscopy

The CI method works in the laboratory frame and requires a certain truncation scheme. As in Ref. Alhassid:2008, , we use as our many-particle basis the set of all Slater determinants which can be constructed from the single-particle states in the harmonic oscillator shells . For a given regularization parameter , the two-body interaction matrix elements in the laboratory frame are calculated by transforming to the relative and center of mass coordinates via the Talmi-Moshinsky brackets Talmi:1952 () and using the interaction (5) in relative motion. We carry out direct diagonalization of the CI Hamiltonian using a new code for three and four particles which works in a two-species formalism with a basis of good total orbital angular momentum and parity . In order to determine the many-particle energies, we take the following limits. For fixed we calculate the energies as a function of and extrapolate to following the method of Ref. Alhassid:2008, . The resulting energies are then studied as a function of to obtain an estimate of the limit.

In the CI calculations we used PARPACK for parallel diagonalization on a Linux cluster using between 1 and 320 CPU cores. The largest system studied, the four-particle system with and total (involving million configurations), required a total of about 2300 GB to store the matrix elements of the Hamiltonian and 24 hours to run on 320 cores, while more typical calculations, e.g., four particles at and ( configurations) required about 64 GB of storage and four hours of computation time on 32 cores.

#### iv.1.1 Three particles

We study a selected set of the low-lying energies of the three-particle system as a function of the inverse scattering length. Here the main topic of interest is the convergence of the eigenvalues (i) as the size of the model space (i.e., the maximal number of oscillator shells ) increases, and (ii) as the regularization parameter increases.

(i) Overall, convergence in is fast for both the renormalized contact interaction and the new effective interaction, and similar to what was already discussed in Ref. Alhassid:2008, in the unitary limit. Convergence slows somewhat as we follow the crossover from the BCS to the BEC regime and tends to require more shells as increases. At , , we find that the lowest four eigenvalues for each of , , , , and are converged in to about in the BCS regime (), in the unitary limit and in the BEC regime (), on average.

(ii) As a function of (for ), we identify two convergence patterns for low-lying three-body eigenvalues. The first occurs for most eigenvalues, and is characterized by fast, monotonic convergence for or at an exponential or faster-than-exponential rate. The values at larger then give small, often nonmonotonic corrections. As an example, the lowest state at unitarity has the exact value . With the effective interaction, we find the energy of this state decreases monotonically to at (slightly overshooting), and at . The convergence then switches direction at and approaches the exact value with a low power-law behavior, so that the values provide only a correction to the value. For such eigenvalues we report the energy at the highest calculated, usually .

The second convergence pattern we encounter is completely monotonic for , but with a noticeably slower rate of convergence. For these eigenvalues, the convergence is smooth, and we estimate the value by fitting , where , to a quadratic polynomial. This provides an interpolating function (with ) which we then use to numerically sum the residual terms in the expansion . An example of such a state is the state at . (We label energy levels as , where indexes the levels with total angular momentum and parity in increasing order starting from zero, so that is the fourth-largest eigenvalue of those with and positive parity.) In this case, the value is and the extrapolated value (from ) is , indicating the state with the exact energy of , a error.

Comparing with the renormalized contact interaction, we find significantly faster convergence with the new interaction for both convergence patterns. For example, the state at unitarity attains the value at with the contact interaction, a difference of from the correct value, compared to at for the new interaction. The difference for faster-converging eigenvalues is not as dramatic, but the new interaction is still clearly preferred. This is shown in Fig. 1 for the state on either side of the crossover, at and . Importantly, we see that the advantage of the new interaction is maintained away from unitarity. This holds true for all the eigenvalues we studied, and for the four-particle system as well.

Because convergence is only truly attained for large , and it is difficult to assess the effects of non-monotonicity, we estimate the systematic error in our three-particle energies along the crossover based on comparison with the unitary calculations, where we compare the convergence in with the known exact results. From this we estimate an upper bound for the systematic error due to non-convergence of approximately in our three-particle energies.

Using the new effective interaction, we also present in Fig. 2 results for a selected set of the low-lying energies of the three-particle system as a function of the inverse scattering length. Although the three-particle system is now well-studied, this provides a clear verification of the capability of the method. Our results agree with those published previously Stetcu:2008 (); Kestner:2007 (). As in the the earlier studies, we find crossings between levels with different quantum numbers, for instance the state and the state at approximately , and the state and the state at approximately , which changes the parity and angular momentum of the ground state.

#### iv.1.2 Four particles

The four-particle system is less well-studied and provides a more interesting test of the new effective interaction. As with the three-particle system, we examined a selected set of the low-lying energies as a function of the inverse scattering length.

(i) Overall we find fast convergence in , very similar to the three-particle case. More precisely, at , , the level of convergence ranges from an average of about at to at for the lowest four states of each of the , , and configurations.

(ii) As a function of , the two convergence patterns observed for the three particle case – the more common fast, slightly non-monotonic convergence and the less common slower, monotonic convergence – are still prevalent, but with the addition of a few eigenvalues which are non-monotonic for smaller . An example is the state at unitarity, which becomes non-monotonic at . For eigenvalues with the first two convergence patterns, we estimate the values in the same manner as the three-particle case. For the new non-monotonic eigenvalues, we report the result of the calculation at largest , usually . For instance, for the state at unitarity, the value is , and appears to be an upper bound based on the direction of convergence. This compares favorably with the result (a difference) of Ref. daily:2010, , in which the center of mass coordinate is separated out and a stochastic variational approach is employed. Note we add a center of mass excitation of to the results of Ref. daily:2010, to allow comparison.

An example of a fast, slightly non-monotonic energy is the state, which was calculated in Ref. daily:2010, to have a value . Our estimate of the new interaction is , well within error.

Compared with the renormalized contact interaction, the new interaction provides much improved convergence for four particles, as it did for the three particle system. For example, Fig. 3 shows the convergence of the state using the new and the renormalized contact interactions on either side of the crossover. As with all the energies we have studied, this displays a marked improvement in convergence. On the BEC side, in particular, we find a difference between the two interactions at .

In Fig. 4 we present results for a selected set of low-lying energies of the four-particle system as a function of inverse scattering length. As for the three-particle case, we find crossings between levels with different quantum numbers, for instance, the state and the state at approximately . In contrast to the three-particle system, the angular momentum and parity of the ground state (here ) remain fixed through the full range of scattering lengths.

### iv.2 Thermodynamics

As thermodynamics is one of the central topics of interest for cold atoms, we test the usefulness of the new separable effective interaction for thermodynamical studies by studying the convergence of the free energy (at a given temperature) versus for the three-particle system. Our calculations were done at temperatures between and as a function of and via a direct diagonalization of the Hamiltonian through . We compare our calculations with exact results for the three-particle system.

An example of the convergence in for and is shown in Fig. 5, where we observe smooth behavior for both interactions. For higher temperatures, the smooth character of the vs. is preserved, while the curvature flattens out to a straight line (on a logarithmic scale), allowing for accurate extrapolations in .

The convergence as a function of is demonstrated in the inset to Fig. 6 for . The convergence using the new effective interaction is significantly faster than for the renormalized contact interaction. In general, we found that the convergence in for the new interaction was quite uniform and fast for temperatures up to approximately . In fact, at higher temperatures, the abundance of interactionless states (see below) diminishes the error for even small values so that at , the difference between and the exact value is only about , compared to at . By comparison, the difference in the free energy between the non-interacting and interacting systems is at and at .

Our thermodynamics studies are summarized in Figs. 6 and 7, where we compare our calculations of the free energy and entropy as a function of temperature with exact calculations for the interacting and non-interacting systems (described below). We find excellent agreement. Although this is only a three-body system, the entropy curve indicates that our results apply in the entirety of the interesting quantum degenerate regime.

Since the three-body problem has been solved exactly at unitarity werner:2006 (); tan3body (), we can compare the thermodynamic functions calculated using the new effective interaction with exact results.

We summarize the exact calculation as follows. There are two classes of eigenstates for the trapped Fermi gas, interacting states and interactionless states. The interacting states satisfy the Bethe-Peierls contact condition werner:2006 ()

 φ(r1,r2,r3)=(1rij−1a)A(Rij,rk)+O(rij)

(as ) via a function which is non-zero. The interactionless states satisfy a similar condition but with identically zero. They are also eigenstates of the non-interacting Hamiltonian and are identical for all values of the scattering length.

To calculate the energy levels of the interacting states we follow the procedure of Ref. werner:2006, . However, as noted there, at large temperatures the number of interactionless states grows much larger than that of the interacting states, so that energies for both are necessary to calculate thermodynamic quantities for the system. The wave functions and energies of the interactionless states for the three-particle system have been characterized completely in Ref. tan3body, . Below we summarize the main results.

The energy of any state of the system (interacting or interactionless) is given by

 E=(γ+nc+2ν+9/2)ℏω

where is the number of oscillator quanta in the center-of-mass excitation, is a scaling exponent, and is the hyperradial quantum number. The corresponding intrinsic wave functions are of the form , where is the hyperradius. The function () satisfies the Bethe-Peierls contact conditions as well as additional equations. States with total internal angular momentum have exponents where . For the interacting states where are an infinite sequence of solutions to an -dependent transcendental equation obtained from the three-particle Schrodinger equation werner:2006 (). For the interactionless states, are integers and have degeneracies; these degeneracies are non-trivial and must be determined in order to calculate thermodynamic quantities for the system. We followed the method of Ref. tan3body, to determine these degeneracies (see, e.g., Table II of Ref. tan3body, ).

Using the above classification of interacting and interactionless states, we computed the entire three-body spectrum to the extent needed to ensure convergence of the free energy at any temperature of interest (see Figs. 6 and 7). Table 1 provides a few values of the free energy and entropy at various temperatures for reference. For the non-interacting calculations we use simple closed formulas that are derived in the Appendix.

## V Conclusion

We have extended a new separable effective interaction for the unitary trapped gas in the configuration-interaction framework to regions away from unitarity, i.e., the BEC-BCS crossover. The main advantage of this interaction, when compared with the conventional renormalized contact interaction, is its much improved convergence of the many-particle energies in the regularization parameter. This allows accurate calculations in the configuration interaction approach. In particular, we calculated the low-lying energies for three and four particles as a function of the inverse scattering length.

We also studied the three-particle system at finite temperature in the unitary limit and found similar convergence properties of the free energy as a function of the regularization parameter. This separable effective interaction may therefore facilitate accurate studies of the thermodynamics of the trapped cold atom system.

###### Acknowledgements.
We acknowledge S. Tan for providing us with Ref. tan3body, and G.F. Bertsch for useful discussions. We also thank I. Stetcu for sharing with us certain three-body calculations using the new interaction. This work was supported in part by the U.S. DOE grant No. DE-FG02-91ER40608, by facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center, and by the NSF grant No. CNS 08-21132 that partially funded acquisition of the facilities. Computational cycles were also provided by the NERSC high performance computing facility at LBL.

## Appendix A Non-interacting thermodynamics

In Figs. 6 and 7 we compared our calculations with the non-interacting three-body trapped systems (dashed lines). For this purpose we derived simple formulas for the canonical partition functions of small systems of trapped non-interacting spin- particles. For three particles (), we proceed as follows. Label the state of each particle as , where are the -D harmonic oscillator quantum numbers and is the -component of the spin. Then the canonical partition function is

 Z=13!∑s1,s2,s3σ1+σ2+σ3=1/2e−β(Es1+Es2+Es3)S(s1,s2,s3)

where and is the selection function

 S(s1,s2,s3) ={0if any two si % are equal1otherwise (16) =∏i

To perform the sum, one expands the selection function to obtain eight terms involving sums over the Maxwell-Boltzmann factor times a product (say ) of zero to three Kronecker ’s. By collapsing the sum over the ’s one obtains a geometric series which then easily yields a closed-form expression. As an example, for , and with the notation , the corresponding sum is

 ∑s1,s2,s3σ1+σ2+σ3=1/2e−β(Es1+Es2+Es3)δs1,s2=(11−e−2β)3(11−e−β)3.

Note that in this case, the sum over and is exactly 1 (and is zero for the two and three- terms for this system).

The final result is

 \omit\span\omit\span\omit\@ADDCLASSltxeqnlefteqn$Z=12e−9β/2[(11−e−β)9−$%(11−e−β)3(11−e−2β)3].

## References

• (1) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
• (2) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
• (3) Q. Chen, J. Stajic, S. Tan, and K. Levin, Phys. Rep. 412, 1-88 (2005).
• (4) J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003).
• (5) G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett. 93, 200404 (2004).
• (6) S. Y. Chang, V. R. Pandharipande, J. Carlson, K. E. Schmidt, Phys. Rev. A 70, 043602 (2004).
• (7) S. Y. Chang and V. R. Pandharipande, Phys. Rev. Lett. 95, 080402 (2005).
• (8) J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. A 76, 053613 (2007).
• (9) D. Blume, J. von Stecher, and C. H. Greene, Phys. Rev. Lett. 99, 233201 (2007).
• (10) R. Jáuregui, R. Paredes, and G. T. Sánchez, Phys. Rev. A 76, 011604(R) (2007)
• (11) J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. A. 77, 043619 (2008)
• (12) A. Gezerlis and J. Carlson, Phys. Rev. C 77, 032801(R) (2008)
• (13) N. T. Zinner, K. Mølmer, C. Özen, D. J. Dean, and K. Langanke, Phys. Rev. A 80, 013613 (2009).
• (14) A. Bulgac, J. E. Drut, and P. Magierski, Phys. Rev. Lett. 96, 090404 (2006); Phys. Rev. Lett. 99, 120401 (2007); Phys. Rev. A 78, 023625 (2008); P. Magierski, G. Wlazlowski, A. Bulgaac, and J. E. Drut, Phys. Rev. Lett. 103, 210403 (2009).
• (15) E. Burovski, N. Prokof’ev, B. Svistunov, and M. Troyer, Phys. Rev. Lett. 96, 160402 (2006).
• (16) S. Y. Chang and G. F. Bertsch, Phys. Rev. A 76, 021603(R) (2007).
• (17) V. K. Akkineni, D. M. Ceperly, and N. Trivedi, Phys. Rev. B 76, 165116 (2007).
• (18) O. Goulko and M. Wingate, Phys. Rev. A 82, 053621 (2010).
• (19) K. M. Daily and D. Blume, Phys. Rev. A 81, 053615 (2010); D. Blume and K. M. Daily, arXiv:1008.3191.
• (20) F. Werner and Y. Castin, Phys. Rev. Lett. 97, 150401 (2006).
• (21) S. Tan, preprint (2006) and private communication.
• (22) J. P. Kestner and L.-M. Duan, Phys. Rev. A 76, 033611 (2007).
• (23) X.-J. Liu, H. Hu, and P. D. Drummond, Phys. Rev. Lett. 102, 160401 (2009); Phys. Rev. A 82, 023619 (2010).
• (24) S. K. Adhikari, Phys. Rev. A 79, 023611 (2009)
• (25) S. Zhang and A. J. Legget, Phys. Rev. A 79, 023601 (2009)
• (26) I. Stetcu, B. R. Barrett, U. van Kolck, and J. P. Vary, Phys.Rev. A 76, 063613 (2007).
• (27) Y. Alhassid, G. F. Bertsch, and L. Fang, Phys. Rev. Lett. 100, 230401 (2008).
• (28) I. Stetcu, J. Rotureau, B.R. Barrett, U. van Kolck, arXiv:1001.5071; J. Rotureau, I. Stetcu, B. R. Barrett, M. C. Birse, and U. van Kolck, Phys. Rev. A 82, 032711 (2010).
• (29) S. Tan, Ann. Phys. (N.Y.) 323, 2952 (2008); Ann. Phys. (N.Y.) 323, 2971 (2008); Ann. Phys. (N.Y.) 323, 2987 (2008).
• (30) S. Jonsell, Few-Body Systems 31, 255 (2002).
• (31) T. Busch, B.-G. Englert, K. Rza̧żewski, and M. Wilkens, Found. Phys. 28, 549 (1998).
• (32) D. E. Knuth, The art of computer programming, vol. 1, fundamental algorithms, (Addison-Wesley, 1997) pp. 36-38, 3rd ed.
• (33) I. Talmi, Helv. Phys. Acta 25, 185 (1952); M. Moshinsky, Nucl. Phys. 13, 104 (1959).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters