Superfluid Bose-Fermi mixture from weak coupling to unitarity S. K. Adhikari1,* and Luca Salasnich2,† 1Instituto de Física Teórica, UNESP—São Paulo State University, 01405-900 São Paulo, SP, Brazil 2CNR-INFM and CNISM, Unit of Padua, Department of Physics “Galileo Galilei,” University of Padua, 35122 Padua, Italy �Received 23 August 2008; published 30 October 2008� We investigate the zero-temperature properties of a superfluid Bose-Fermi mixture by introducing a set of coupled Galilei-invariant nonlinear Schrödinger equations valid from weak coupling to unitarity. The Bose dynamics is described by a Gross-Pitaevskii-type equation including beyond-mean-field corrections possessing the correct weak-coupling and unitarity limits. The dynamics of the two-component Fermi superfluid is de- scribed by a density-functional equation including beyond-mean-field terms with correct weak-coupling and unitarity limits. The present set of equations is equivalent to the equations of generalized superfluid hydrody- namics, which take into account also surface effects. The equations describe the mixture properly as the Bose-Bose repulsive �positive� and Fermi-Fermi attractive �negative� scattering lengths are varied from zero to infinity in the presence of a Bose-Fermi interaction. The present model is tested numerically as the Bose-Bose and Fermi-Fermi scattering lengths are varied over wide ranges covering the weak-coupling to unitarity transition. DOI: 10.1103/PhysRevA.78.043616 PACS number�s�: 03.75.Ss, 03.75.Hh I. INTRODUCTION The macroscopic quantum phenomenon of superfluidity has been clearly demonstrated in recent experiments with ultracold and dilute gases made of alkali-metal atoms �1,2�. Both bosonic and fermionic ultracold and dilute superfluids can be accurately described by the zero-temperature hydro- dynamical equations of superfluids �1–5�. Trapped Bose-Fermi mixtures, with Fermi atoms in a single hyperfine state �normal Fermi gas�, were investigated by various authors both theoretically �6–14� and experimen- tally �15–17�. Recently, there was a study of a dilute mixture of superfluid bosons and fermions across a Feshbach reso- nance of the Fermi-Fermi scattering length af, obtaining the phase diagram of the mixture in a box �18�. In the strict one-dimensional case, we found that the superfluid Bose- Fermi mixture exhibits phase separation and solitons by changing the Bose-Fermi interaction �19,20�. In the present paper we analyze in detail a three- dimensional �3D� superfluid Bose-Fermi mixture under har- monic trapping confinement when the Bose-Bose repulsive �positive� and Fermi-Fermi attractive �negative� scattering lengths are varied from very small �weak coupling� to infi- nitely large �strong coupling� values. The variation of these scattering lengths from weak to strong coupling is achieved in laboratory �21� by varying a background magnetic field near a Feshbach resonance. In the first part of the paper we derive Galilei-invariant nonlinear Schrödinger equations for superfluid Bose and Fermi dynamics, valid from weak-coupling to unitarity in each case, which are equivalent to the hydrodynamical equa- tions. For bosons the equation is a generalized �22–24� zero- temperature Gross-Pitaevskii �GP� equation �25� including beyond-mean-field corrections to incorporate correctly the effect of bosonic interaction for large and positive �repulsive� scattering lengths ab. In the ordinary GP equation the non- linear term is linearly proportional to ab and hence highly overestimates the effect of bosonic interaction for large posi- tive values of ab. In fact there is a saturation of bosonic interaction in the so-called unitarity ab→ +� limit, properly taken care of in the present generalized GP equations. For superfluid fermions the present formulation is based on the density-functional �DF� theory �26,27� for fermion pairs �23,28–31� including beyond-mean-field corrections to in- corporate properly the effect of fermionic attraction between spin-up and -down fermions especially for large negative �at- tractive� values of Fermi-Fermi scattering length af. The DF theory, in different forms, has already been applied to the problem of ultracold fermions �32,33�. The usual DF equa- tion is valid for the weak-coupling Bardeen-Cooper- Schrieffer �BCS� limit: af →−0. The present generalized DF equation correctly accounts for the dynamics as the Fermi- Fermi attraction is varied from the weak-coupling BCS limit to unitarity: af →−�. In the unitarity limit there is a satura- tion of the effect of fermion interaction, properly taken care of in the present set of equations, which seems to be appro- priate to study the crossover �34� from the weak-coupling BCS limit to the molecular Bose limit. In the second part of the paper, by including the interac- tion between a boson and fermion pair, we obtain a set of coupled equations for superfluid Bose-Fermi dynamics in the presence of a Bose-Fermi interaction. The model correctly describes the dynamics as the Bose-Bose repulsion and Fermi-Fermi attraction are varied from weak coupling to uni- tarity. This model is tested numerically as the different scat- tering lengths are varied over wide ranges covering the weak-coupling to unitarity transition for a superfluid Bose- Fermi mixture. In Sec. II we present the hydrodynamical equations for bosons and fermions and present the appropriate bulk chemi- cal potentials valid in the weak-coupling limit as well as possessing the saturation in the strong-coupling unitarity limit consistent with the constraints of quantum mechanics. *adhikari@ift.unesp.br; URL www.ift.unesp.br/users/adhikari †salasnich@pd.infn.it; URL www.padova.infm.it/salasnich PHYSICAL REVIEW A 78, 043616 �2008� 1050-2947/2008/78�4�/043616�10� ©2008 The American Physical Society043616-1 http://dx.doi.org/10.1103/PhysRevA.78.043616 In Sec. III we derive mean-field equations for bosons and fermions consistent with the hydrodynamical equations for a large number of atoms. Next we derive the present model for interacting Bose-Fermi superfluid introducing a contact inter- action between bosons and fermions. In Sec. IV we perform numerical calculations for densities and chemical potentials of a trapped Bose-Fermi superfluid mixture to show the ad- vantage of the present model. Finally, in Sec. V we present a brief summary and conclusion. II. SUPERFLUID HYDRODYNAMICS FOR BOSONS AND FERMIONS At zero temperature, for a large number of atoms, statical and dynamical collective properties of bosonic and fermionic superfluids are expected to be properly described by the hy- drodynamical equations of superfluids �1,2�. For bosons the hydrodynamical equations are given by �1,3� � �t nb + � · �nbvb� = 0, �1� mb � �t vb + ��1 2 mbvb 2 + Ub + �b�nb,ab�� = 0, �2� where Ub�r� is the external potential acting on bosons, mb is the mass of a bosonic atom, nb�r , t� is the local density of bosons, and vb�r , t� is the local superfluid velocity �1–3�. The total number of bosons is given by Nb =� nb�r,t�d3r , �3� and the nonlinear term �b�nb ,ab� is the bulk chemical poten- tial of the bosonic system with ab the Bose-Bose scattering length. Equation �1� is the equation of continuity of hydro- dynamic flow, while Eq. �2� is the equation of conservation of the linear momentum. Equation �2� establishes the irrota- tional nature of the superfluid motion: ��vb=0, meaning that the velocity vb can be written as the gradient of a scalar field. Equations �1� and �2� differ from the corresponding equations holding in the collisional regime of a nonsuper- fluid system because of the irrotationality constraint. For superfluid fermions, the fundamental entity governing the superfluid flow is the Cooper pair. In the case of a two- component �spin up and down� superfluid Fermi system one can introduce the local density of pairs np�r , t�=nf�r , t� /2, where nf�r , t� is the total local density of fermions. The hy- drodynamical equations of a fermionic superfluid �2–5� can be written as � �t np + � · �npvp� = 0, �4� mp � �t vp + ��1 2 mpvp 2 + Up + �p�np,af�� = 0, �5� where the nonlinear term �p�np ,af� is the bulk chemical po- tential of pairs with af the attractive Fermi-Fermi scattering length. Here mp is the mass of a pair, that is twice the mass mf of a single fermion, i.e.; mp=2mf. In addition, the trap potential Up�r� is twice the trap potential Uf�r� acting on a single fermion, i.e., Up�r�=2Uf�r�. Note that the chemical potential �p�np ,af� is twice the total chemical potential � f�nf ,af� of the Fermi system, i.e., �p�np ,af�=2� f�2np ,af�. The total number Nf of fermions is Nf = 2Np = 2� np�r�d3r =� nf�r�d3r , �6� where Np is the number of pairs. It is important to stress that Eqs. �1� and �2� are formally similar to Eqs. �4� and �5�, but the quantities involved are strongly different. In particular, the bulk chemical potential of bosons is completely different from the bulk chemical potential of Fermi pairs �1,2,35�. In the full crossover from the small-gas-parameter regime to the large-gas-parameter regime we use the following ex- pression for the bulk chemical potential of the Bose super- fluid �22�: �b�nb,ab� = �2 mb nb 2/3f�nb 1/3ab� , �7� where f�x� = 4� x + �x5/2 1 + �x3/2 + �x5/2 . �8� In the small-gas-parameter �x→0� regime, Eq. �8� becomes f�x� = 4��x + �� − ��x5/2 + ¯ � , �9� which is the analytical result for bulk chemical potential found by Lee, Yang, and Huang �36,37� in this limit, pro- vided that ��−��=32 / �3���. In Eq. �8� we shall also choose �=4�� / , with =22.22; this will make the bulk chemical potential �7� satisfy the correct unitarity limit �b�nb ,ab� =22.22�2nb 2/3 /mb as established by Cowell et al. �38�. Hence Eq. �7� with Eq. �8� can be made to satisfy the correct weak- coupling and unitarity limits. Next we need to determine the constants � and � consistent with ��−��=32 / �3���. We consider the following three choices for � and � from many possible choices: �=32 / �3���, �=32� −1� / �3��� with �a� =1.05, �b� =1.1, and �c� =1.15. All three choices are consistent with the unitarity and the Lee-Yang-Huang limits �36,37�. We shall present a critical numerical study of the three choices in the next section. In the full crossover from BCS regime �af →−0� to uni- tarity �af →−�� we use the following expression for the bulk chemical potential of the Fermi superfluid �39�: �p�np,af� = 2 �2 mp �6�2np�2/3g�21/3np 1/3af� , �10� where g�x� = 1 + �x 1 − �x , �11� where � and � are fitting parameters. The close agreement of the interpolation function g�x� for negative x values with the results of Monte Carlo calculation of energy of a uniform gas S. K. ADHIKARI AND LUCA SALASNICH PHYSICAL REVIEW A 78, 043616 �2008� 043616-2 of Fermi superfluid with the parameters �=4� / �3�2�2/3, � =� / �1− �, and =0.44 �40,41� was established in Ref. �39�. The inclusion of the lowest-order term g�x�=1 in Eq. �10� leads to the bulk chemical potential in the absence of Fermi- Fermi interaction �af =0� of a uniform gas. The next-order term g�x�=1+�x includes a known analytical result in the small-gas-parameter regime as obtained by Lee and Yang �36,42� and by Galitskii �43�. The model �10� with Eq. �11� provides a smooth interpolation between the bulk chemical potential in the BCS limit �p�np ,af�=2�2�6�2np�2/3 /mp +16�np�2af /mp �36� for small np 1/3af values and that in the unitarity limit �39� �p�np ,af�=2�2�6�2np�2/3 /mp for large np 1/3af values. Manini and Salasnich �29� and Kim and Zubarev �28� also proposed similar interpolation formulas on the basis of Monte Carlo data of a uniform Fermi superfluid. Here we deal with a Fermi gas in a spherical harmonic trap rather than a uniform Fermi gas. By a direct comparison with the results of Monte Carlo calculation of a superfluid Fermi gas in a spherical harmonic trap, we shall show in Sec. III that Eqs. �10� and �11� still present a good approximation to the energy of the system, but now with a slightly different value of the parameters. The hydrodynamical equations are valid to describe equi- librium properties and dynamical properties of long wave- length for both bosons and fermions. In particular, one can introduce �1,2� a healing �or coherence� length such that the transport phenomena under investigation must be character- ized by a length scale much larger than the healing length. As suggested by Combescot, Kagan, and Stringari �44�, the healing length can be defined for bosons as �b = � mbvb cr , �12� where vb cr is the Landau critical velocity above which the system gives rise to energy dissipation. This critical velocity coincides with the first sound velocity and is given by vb cr =� nb mb ��b �nb . �13� In Fig. 1 we plot the critical sound velocity vb cr of a uniform Bose superfluid for the present �b given by Eqs. �7� and �8� with the parameter =1.1. The interesting feature of this plot is that vb cr=0, when either nb or ab is zero. However, vb cr increases monotonically with nb, whereas vb cr attains a con- stant value as ab increases. These features are present in the plot of Fig. 1. However, the critical sound velocity vb cr cal- culated from �b given by Eqs. �7� and �9� increases mono- tonically with both nb and ab. For superfluid fermions the healing length of Cooper pairs can be defined as �p = � mpvp cr , �14� where the critical velocity vp cr is related to the breaking of pairs through the formula vp cr =���p 2 + � 2 − �p mp , �15� where � is the energy gap �2,44�. We notice that in the deep BCS regime of weakly interacting attractive Fermi atoms �corresponding to � ��p� Eq. �15� approaches the expo- nentially small value vp cr= � /�mp�p /2. III. NONLINEAR SCHRÖDINGER EQUATIONS The bosonic superfluid can be described by a GP �1,25� complex order parameter �b�r , t� given by �b�r,t� = �nb�r,t�ei�b�r,t�. �16� The probability current density j��r , t� is given by j� nbvb = � 2imb �� b * � �b − �b � � b *� . �17� Equations �16� and �17� relate the phase �b�r , t� to the super- fluid velocity field vb by the formula �1� vb = � mb � �b. �18� The bosonic order parameter �b�r , t� is, apart from a nor- malization �4,5�, nothing but the condensate wave function �b�r,t� = ��̂�r,t�� , �19� that is, the expectation value of the bosonic field operator �̂�r , t� �4,23�. The order parameter �b�r , t� satisfies the fol- lowing generalized GP nonlinear Schrödinger equation �1� i� � �t �b = �− �2 2mb �2 + Ub + �b�nb,ab���b, �20� where the bulk chemical potential is given by Eq. �7�. The use of f�x�=4�x in Eq. �7� leads to the GP equation and the use of the two leading terms of Eq. �9�, e.g., f�x� = 4��x + 32 3�� x5/2� �21� in Eq. �7� leads to the modified GP �MGP� equation sug- gested by Fabrocini and Polls �24�, which includes correction to the GP equation for medium values of Bose-Bose scatter- ing length ab. The use of Eq. �8� in Eq. �7� leads to a non- 0 5 10 04812 nb 0 1 2 ab 0 5 10 vb cr(nb,ab) FIG. 1. �Color online� The critical sound velocity vb cr of a uni- form Bose superfluid as a function of Bose-Bose scattering length ab and density nb. SUPERFLUID BOSE-FERMI MIXTURE FROM WEAK … PHYSICAL REVIEW A 78, 043616 �2008� 043616-3 linear Schrödinger equation valid from the weak-coupling to the unitarity limit and is called the unitary Schrödinger �US� equation for bosons �22�. In Fig. 2 we plot the various possibilities for the interpo- lation function f�x�; e.g., those corresponding to the US equation �8� with �a� =1.05, �b� =1.1, and �c� =1.15 �as suggested after Eq. �9��, the MGP equation �21�, the GP equation f�x�=4�x. The asymptotic limit limx→�f�x� =22.22 is also shown. All three choices �a�, �b�, and �c� represent smooth interpolation of f�x� between small and large x values. For small x this choice agrees with the MGP result �21� and for large x with the asymptotic result. To test these three choices we actually solve the US Eq. �20�, by imaginary time propagation after a Crank-Nicholson discreti- zation �45–47� �detailed in the beginning of Sec. IV�, for the trap parameters of a possible experimental setup for 87Rb atoms in a spherical trap and compare the results for energy with those obtained by Blume and Greene �48� by the diffu- sion Monte Carlo �DMC� method. Actually, we solved Eqs. �20� and �7� with nb normalized to �Nb−1� in place of Eq. �3�, as appropriate for a small number of bosons, cf. Eq. �2� of �48�. The energy of the calculation is performed for a boson-boson scattering length ab=0.433l where l is the har- monic oscillator length for various Nb and the results for energy are displayed in Table I. From this table we find that the present US models always provide a much better ap- proximation to the energies than the GP and MGP models. More results of energy for larger number of bosons may help in fixing the parameters of the present model more accu- rately. In the present paper we shall use the US model �b� with =1.1. The relationship between Eq. �20� and the hydrodynami- cal equations �1� and �2� can be established by inserting Eq. �16� into Eq. �20�. After some straightforward algebra equat- ing the real and imaginary parts of both sides and taking into account Eq. �18�, one finds two generalized hydrodynamical equations �1�, � �t nb + � · �nbvb� = 0, �22� mb � �t vb + ��− �2 2mb �2�nb �nb + mbvb 2 2 + Ub + �b�nb,ab�� = 0, �23� which include the quantum pressure term Tb QP = − �2 2mb �2�nb �nb , �24� which depends explicitly on the reduced Planck constant �. Neglecting the quantum pressure term, one gets from Eqs. �22� and �23� the classical hydrodynamical equations �1� and �2�. Equation �22� is the continuity equation whereas Eq. �23� establishes the irrotational nature of hydrodynamic flow. Similarly, the fermionic superfluid can be described by a DF �32� complex order parameter �p�r , t� of pairs, with a modulus and a phase �p�r , t� such that �p�r,t� = �np�r,t�ei�p�r,t�. �25� The expression for the probability current density of fermi- ons leads to the following expression for the velocity field vp of pairs �2�: vp = � mp � �p. �26� Note that the mass of a pair mp, and not that of a fermion, appears in the denominator of the phase-velocity relation. The fundamental entity responsible for superfluid flow has the mass mp of a pair of fermions. For paired fermions in the superfluid state, the order parameter is, apart from a normal- 0 5 10 15 20 25 0 2 4 6 8 10 f( x) x US Eq. (8); ν=1.05 0 5 10 15 20 25 0 2 4 6 8 10 f( x) x US Eq. (8); ν=1.05 US Eq. (8); ν=1.1 0 5 10 15 20 25 0 2 4 6 8 10 f( x) x US Eq. (8); ν=1.05 US Eq. (8); ν=1.1 US Eq. (8); ν=1.15 0 5 10 15 20 25 0 2 4 6 8 10 f( x) x US Eq. (8); ν=1.05 US Eq. (8); ν=1.1 US Eq. (8); ν=1.15 MGP Eq. (21) 0 5 10 15 20 25 0 2 4 6 8 10 f( x) x US Eq. (8); ν=1.05 US Eq. (8); ν=1.1 US Eq. (8); ν=1.15 MGP Eq. (21) GP 0 5 10 15 20 25 0 2 4 6 8 10 f( x) x US Eq. (8); ν=1.05 US Eq. (8); ν=1.1 US Eq. (8); ν=1.15 MGP Eq. (21) GP asymptotic: 22.22 FIG. 2. �Color online� The Bose interpolation function f�x� given by Eq. �8� with �=32 / �3���, �=32� −1� / �3���, and � =4�� /22.22 with = �a� 1.05, �b� 1.1, and �c� 1.15. The results of MGP �Eq. �21�� and GP �f�x�=4�x� models together with the asymptotic value are also given. TABLE I. Ground state energies for different number Nb of bosonic atoms in a spherical trap and for ab=0.433 from a solution of GP equation, MGP equation, and US equation for = �a� 1.05, �b� 1.1, and �c� 1.15. The results for energies obtained with two potentials of DMC calculation �48� are quoted for compari- son. Length and energies are expressed in oscillator units. Nb DMC GP MGP US �a� US �b� US �c� 3 5.553�3�; 5.552�2� 5.329 5.611 5.570 5.564 5.558 5 10.577�2�; 10.574�4� 9.901 10.772 10.629 10.608 10.588 10 26.22�8�; 26.20�8� 23.61 26.84 26.24 26.16 26.07 20 66.9�4�; 66.9�1� 57.9 68.5 66.38 66.07 65.77 S. K. ADHIKARI AND LUCA SALASNICH PHYSICAL REVIEW A 78, 043616 �2008� 043616-4 ization �4,5,49�, the condensate wave function of the center of mass of the Cooper pairs: �p�r,t� = ��̂↑�r,t��̂↓�r,t�� , �27� that is the average of pair operators, with �̂��r , t� the fermi- onic field operator with spin component �= ↑ ,↓ �4,23�. The order parameter �p�r , t� satisfies the following zero- temperature nonlinear DF Schrödinger equation �23,28,29,31,39�: i� � �t �p = �− �2 2mp �2 + Up + �p�np,af���p, �28� where the bulk chemical potential �p�np ,af� is given by Eq. �10�. The use of g�x�=1 in �p�np ,af� of Eq. �10� leads to the lowest order DF equation. If we use the next order solution of Eq. �11� in Eq. �10�, e.g., g�x� = 1 + �x , �29� we obtain a modified density-functional �MDF� equation with corrections for medium value Fermi-Fermi scattering length af. The use of Eq. �11� in Eq. �10� leads to a DF equation valid in the weak-coupling BCS to unitarity limit and the resultant nonlinear equation �28� will be called the unitary Schrödinger �US� equation for fermion pairs. In Fig. 3 we plot the various possibilities for the interpo- lation function g�x�; e.g., that corresponding to Eq. �11�, the MDF equation �29� for two different values of �: �i� 4� / �3�2�2/3 and �ii� 20� / �3�2�2/3. The parameter � is al- ways taken as �=� / �1− �, =0.44 �40,41�. The asymptotic limit limx→−� g�x�=0.44 is also shown in Fig. 3. The expres- sion �11� represents a good approximation of g�x� for small and large x . For small x this choice agrees with the MDF result �29� and for large x with the asymptotic result for both choices of �. To test the above two choices �i� and �ii� of the parameter � in Eq. �11� for a Fermi superfluid in a spherical trap, we calculate the energy of the system for different values of the scattering length af and number of atoms Nf by solving di- rectly Eq. �28� by imaginary time propagation after a Crank- Nicholson discretization �45–47� �detailed in the beginning of Sec. IV�. We also compare the results with those obtained by the fixed-node Monte Carlo calculation �FNMC� �50–53�. The results of our investigation are shown in Figs. 4. In Fig. 4�a� we plot energy E vs af for Nf =4 and 8 obtained from a solution the US equation �28� with the choices �i� 4� / �3�2�2/3 and �ii� 20� / �3�2�2/3 for � and the FNMC data of Refs. �52,53�. The present hydrodynamical formulation is expected to be good for a large number of fermions �see, Fig. 4�b��, yet for a small number Nf =4 the result is reason- able. In Fig. 4�b� we plot energy E /Nf 2/3 vs Nf 2/3 at unitarity af →−� obtained from a solution of the US equation �28� with choice �ii� 20� / �3�2�2/3 for �. In this limit both choices of � lead to the same energy. The results of FNMC �50� and Green function Monte Carlo �GFMC� �54� calculations as well as of local density approximation �LDA� are also shown in Fig. 4�b�. The LDA result is analytically known in this case as E�N�= �3Nf�4/3�� /4, �=0.44 �33�. We plot E /Nf 2/3 vs Nf 2/3 in Fig. 4�b� because of the linear correlation between these two variables explicit in the LDA result. From Figs. 4 we find that the results with choice �ii� 20� / �3�2�2/3 of � 0.4 0.5 0.6 0.7 0.8 0.9 1 -10 -8 -6 -4 -2 0 g( x) x US Eq. (11) (i) 0.4 0.5 0.6 0.7 0.8 0.9 1 -10 -8 -6 -4 -2 0 g( x) x US Eq. (11) (i) MDF Eq. (29) (i) 0.4 0.5 0.6 0.7 0.8 0.9 1 -10 -8 -6 -4 -2 0 g( x) x US Eq. (11) (i) MDF Eq. (29) (i) US Eq. (11) (ii) 0.4 0.5 0.6 0.7 0.8 0.9 1 -10 -8 -6 -4 -2 0 g( x) x US Eq. (11) (i) MDF Eq. (29) (i) US Eq. (11) (ii) MDF Eq. (29) (ii) 0.4 0.5 0.6 0.7 0.8 0.9 1 -10 -8 -6 -4 -2 0 g( x) x US Eq. (11) (i) MDF Eq. (29) (i) US Eq. (11) (ii) MDF Eq. (29) (ii) asymptotic: 0.44 FIG. 3. �Color online� The Fermi interpolation function g�x� given by Eq. �11� with �i� �=4� / �3�2�2/3 and �ii� 20� / �3�2�2/3 and �=� /0.56. The MDF results given by Eq. �29� and the asymptotic value are also shown. 10 15 20 0.001 0.01 0.1 1 10 100 1000 E |af| X 2 US (Nf=8) (ii) 10 15 20 0.001 0.01 0.1 1 10 100 1000 E |af| X 2 US (Nf=8) (ii) US (Nf=8) (i) 10 15 20 0.001 0.01 0.1 1 10 100 1000 E |af| X 2 US (Nf=8) (ii) US (Nf=8) (i) FNMC (Nf=8) 10 15 20 0.001 0.01 0.1 1 10 100 1000 E |af| X 2 US (Nf=8) (ii) US (Nf=8) (i) FNMC (Nf=8) US (Nf=4) (ii) 10 15 20 0.001 0.01 0.1 1 10 100 1000 E |af| X 2 US (Nf=8) (ii) US (Nf=8) (i) FNMC (Nf=8) US (Nf=4) (ii) FNMC (Nf=4) 0 2 4 6 8 0 2 4 6 8 10 E /N f2/ 3 Nf 2/3 US 0 2 4 6 8 0 2 4 6 8 10 E /N f2/ 3 Nf 2/3 US FNMC 0 2 4 6 8 0 2 4 6 8 10 E /N f2/ 3 Nf 2/3 US FNMC GFMC 0 2 4 6 8 0 2 4 6 8 10 E /N f2/ 3 Nf 2/3 US FNMC GFMC LDA (b) (a) FIG. 4. �Color online� �a� Energy of a superfluid Fermi gas in a spherical trap in oscillator units vs Fermi-Fermi scattering length af also in oscillator units obtained from the solution of Eq. �28� for the choices �= �i� 4� / �3�2�2/3 and �ii� 20� / �3�2�2/3 and from the fixed node Monte Carlo �FNMC� calculation �52,53�. �b� Energy of a superfluid Fermi gas in a spherical trap in oscillator units vs number of atoms Nf in the unitarity limit af →−� obtained from a solution of Eq. �28�, from FNMC �50� and Green function Monte Carlo �GFMC� calculations �54�, and from local density approximation �LDA�. SUPERFLUID BOSE-FERMI MIXTURE FROM WEAK … PHYSICAL REVIEW A 78, 043616 �2008� 043616-5 agree well with the Monte Carlo data in both cases and this choice will be used in the present study. The Monte Carlo data clearly favor the present model over the LDA. Again, inserting Eq. �25� into Eq. �28�, after some straightforward algebra equating the real and imaginary parts and taking into account Eq. �26�, we find two generalized hydrodynamical equations for the Fermi superfluid, � �t np + � · �npvp� = 0, �30� mp � �t vp + � − �2 2mp �2�np �np + mpvp 2 2 + Up + �p�np,af�� = 0, �31� which include the quantum pressure term Tp QP = − �2 2mp �2�np �np , �32� involving �. Neglecting this term leads to the classical hy- drodynamical equations �4� and �5� of superfluid Fermi pairs. This specific quantum pressure term is a consequence of the proper phase-velocity relation �26� with the factor mp in the denominator �2�. Had we taken a different mass factor in the denominator, e.g., mf, a different quantum pressure term would have emerged. The present choice is physically moti- vated by Galilei invariance �2,23,31� and leads to energies of trapped Fermi superfluid in close agreement �55� with those obtained �50,54� from the Monte Carlo methods. We stress that, from the point of view of DF theory, the quantum pressure terms for superfluid bosons and fermions correspond to gradient correction �also called surface correc- tions� to the local density approximation �LDA�, where LDA gives exactly the classical hydrodynamical equations ob- tained by setting the gradient kinetic-energy term to zero. The gradient term, which can also be seen as a next-to- leading contribution in a momentum expansion to the classi- cal superfluid hydrodynamics �56�, takes into account correc- tions to the kinetic energy due to spatial variations in the density of the system, and it is essential to have regular be- havior at the surface of the cloud where the density goes to zero �57�. In the case of superfluid bosons, the quantum- pressure term, Eq. �24�, is exactly the one of the familiar cubic GP equation �1,25�. In the case of superfluid fermions the gradient term is consistent with the motion of paired fermions of mass mp=2mf. For normal fermions various au- thors have proposed and used different gradient terms �57–61�. For superfluid fermions we are using the familiar von Weizsäcker term �58� in the case of pairs. In fact, in our approach, the explicit form of the quantum-pressure term, Eq. �32�, is strictly determined by the Galilei-invariant rela- tion, Eq. �26�, between phase and superfluid velocity �2,23,31�. Finally, we observe that Eq. �32� is practically the same as that recently obtained in Ref. �62� in the case of a superfluid Fermi gas at unitarity using an � expansion around d=4−� spatial dimensions. The generalized hydrodynamical equations �22� and �23� for bosons and Eqs. �30� and �31� for fermions are thus com- pletely equivalent to the respective nonlinear Schrödinger equations �20� and �28�. In the limit of zero velocity fields, vb, vp→0, Eqs. �23� and �31� yield the stationary versions of Eqs. �20� and �28� describing the superfluid densities: �b0 �nb = �− �2 2mb �2 + Ub + �b�nb,ab���nb, �33� �p0 �np = �− �2 2mp �2 + Up + �p�np,af���np, �34� showing the agreement between hydrodynamical and nonlin- ear equations for superfluid bosons and fermions. Here �b0 and �p0 are the chemical potentials for bosons and fermions in the stationary state. Although Eq. �34� is written in terms of Fermi pair den- sity np �and we use this equation in the following�, an equivalent equation for Fermi density nf �=2np� can be writ- ten as � f0 �nf = �− �2 8mf �2 + Uf + � f�nf,af���nf , �35� � f�nf,af� = �2 2mf �3�2nf�2/3g�nf 1/3af� , �36� where we have used �p0=2� f0, Up=2Uf, and �p�np ,af� =2� f�nf ,af�. Equations �20� and �28� are the Euler-Lagrange equations of the following Lagrangian density: L j = i� 2 �� j * � �t � j − � j � �t � j *� − �2 2mj �� j 2 − Uj�r� � j 2 − E j�nj,aj� � j 2, �37� with j=b , p �here ap=af, the Fermi-Fermi scattering length� where E j�nj,aj� = 1 nj � 0 nj � j�nj�,aj�dnj�, �38� is the bulk energy per particle of the superfluid. Finally, we introduce a Lagrangian density for Bose- Fermi interaction to Eq. �37� of the form Lbp = Gbp �b 2 �p 2, �39� where Gbp=4��2abf /mbf, and mbf is the Bose-Fermi reduced mass and abf is the Bose-Fermi scattering length. We note that there are no coupling terms between the hydrodynamical equations �23� and �31� involving the veloc- ity of two superfluid components. The only coupling in the Lagrangian density �39� involves the density of the two su- perfluid components. However, the present hydrodynamical equations �23� and �31� are a special case of the three-fluid �two superfluids and a normal fluid� hydrodynamical formu- lation of Landau and Khalatnikov �63� �appropriate for the 4He-3He mixture�. The Galilei-invariant Landau-Khalatnikov formulation in its general form also has a coupling term in- S. K. ADHIKARI AND LUCA SALASNICH PHYSICAL REVIEW A 78, 043616 �2008� 043616-6 volving two superfluid velocities �not considered in this pa- per�, as discussed �64� by Andreev and Bashkin and also by Ho and Shenoy. However, the present static model defined by Eqs. �33� and �34� is obtained by setting all velocities equal to zero, and hence these terms depending on superfluid velocities would leave the present model unchanged. In this paper we consider abf to have not too large positive values only �moderately repulsive Bose-Fermi interaction�. There are other domains of Bose-Fermi interaction which require special attention. If this interaction is attractive and strong enough, then the ground state is very different, e.g., the composite fermions f̃ =bf are created at higher tempera- tures and at zero temperature we have Cooper pairs �or local pairs� of the type � f̃ f̃� consisting of two elementary fermions and two elementary bosons �65�. The complete Lagrangian density given by Eqs. �37� and �39� leads to the following set of coupled equations for the Bose-Fermi-Fermi mixture made of superfluid bosons and two-component superfluid fermions: i� � �t �b = �− �2�2 2mb + Ub + �b�nb,ab� + Gbp �p 2��b, �40� i� � �t �p = �− �2�2 2mp + Up + �p�np,af� + Gbp �b 2��p, �41� valid from weak-coupling to unitarity for both bosons and fermions. The coupled set of equations �40� and �41� with bulk chemical potentials �b�nb ,ab� and �p�np ,af� given by Eqs. �7� and �8� �with =1.1� and Eqs. �10� and �11� �with �=20� / �3�2�2/3, �=� / �1− �, =0.44�, respectively, is the US equation for the interacting superfluid Bose-Fermi mix- ture and is the required set. However, when the bulk chemi- cal potentials �b�nb ,ab� and �p�np ,af� are approximated us- ing Eqs. �21� and �29�, the set of equations �40� and �41� will be termed modified Gross-Pitaevskii density-functional �MGPDF� equation including some correction for medium values of scattering lengths ab and af. Finally, if we approxi- mate f�x�=4�x and g�x�=1 in Eqs. �40� and �41� with Eqs. �8� and �11�, the resultant model is termed the Gross- Pitaevskii density-functional �GPDF� model. We present a comparative numerical study of these three models—US, MGPDF, and GPDF—in the next section. IV. NUMERICAL RESULTS Here we solve numerically the coupled set of US Eqs. �40� and �41� with bulk chemical potentials �b and �p given by Eqs. �8� and �11� and compare the results so obtained with those of the MGPDF and GPDF models. We numerically solve the US, MGPDF, and GPDF partial differential equa- tions by discretizing them by the semi-implicit Crank- Nicholson algorithm with imaginary time propagation �45–47�. The space and time steps used in discretization were typically 0.05 and 0.001, respectively. Although we are not interested here in fitting experimen- tal data, we perform numerical calculation for variables for a possible experimental setup �16� of a 87Rb-40K mixture in spherical harmonic traps satisfying Ub=mb�b 2r2 /2 =mf� f 2r2 /2=Uf with an oscillator length l=�� / �mb�b� =1 �m for 87Rb atoms. The oscillator length for 40K is modified according to the trap. The fermion pair trap in Eq. �41� is Up=mf� f 2r2. In Figs. 5�a� and 5�b� we plot densities for the superfluid Bose-Fermi mixture as calculated from Eqs. �40� and �41� for Nb=10 000, Np=10 000, af =0, abf =100 nm with �a� ab =200 nm and �b� ab=2000 nm for the US, MGPDF, and GPDF models. For small ab both MGPDF, and GPDF mod- els are good approximations to the bosonic density of the US model �not shown in the figure�. But as ab increases, the MGPDF model is a better approximation to the bosonic den- sity of the US model than the GPDF model as can be seen from Fig. 5�a�. However, for very large ab, neither the MG- PDF nor the GPDF model provides the saturation of bulk chemical potential for bosons as in the US model. Conse- quently, the MGPDF and GPDF models acquire large bosonic nonlinearity compared to the US model and hence produce bosonic density extending to a larger region in space. As ab is increased, the MGPDF model yields rapidly divergent bosonic nonlinearity �compared to the GPDF model� and hence produces bosonic density inferior to the GPDF model as can be seen in Fig. 5�b�. In Figs. 6�a� and 6�b� we plot superfluid densities for Nb=100 000, Np=100, abf =100 nm, ab=5 nm with �a� af =−100 nm and �b� af =−300 nm for the US, MGPDF, and 0 0.0004 0.0008 0.0012 0 2 4 6 8 10 12 14 16 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=200 nm af=0 nm nb, US 0 0.0004 0.0008 0.0012 0 2 4 6 8 10 12 14 16 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=200 nm af=0 nm nb, US np, US 0 0.0004 0.0008 0.0012 0 2 4 6 8 10 12 14 16 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=200 nm af=0 nm nb, US np, US nb, MGPDF 0 0.0004 0.0008 0.0012 0 2 4 6 8 10 12 14 16 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=200 nm af=0 nm nb, US np, US nb, MGPDF np, MGPDF 0 0.0004 0.0008 0.0012 0 2 4 6 8 10 12 14 16 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=200 nm af=0 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF 0 0.0004 0.0008 0.0012 0 2 4 6 8 10 12 14 16 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=200 nm af=0 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF np, GPDF 0 0.0002 0.0004 0.0006 0.0008 0 2 4 6 8 10 12 14 16 18 20 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=2000 nm af=0 nm nb, US 0 0.0002 0.0004 0.0006 0.0008 0 2 4 6 8 10 12 14 16 18 20 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=2000 nm af=0 nm nb, US np, US 0 0.0002 0.0004 0.0006 0.0008 0 2 4 6 8 10 12 14 16 18 20 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=2000 nm af=0 nm nb, US np, US nb, MGPDF 0 0.0002 0.0004 0.0006 0.0008 0 2 4 6 8 10 12 14 16 18 20 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=2000 nm af=0 nm nb, US np, US nb, MGPDF np, MGPDF 0 0.0002 0.0004 0.0006 0.0008 0 2 4 6 8 10 12 14 16 18 20 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=2000 nm af=0 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF 0 0.0002 0.0004 0.0006 0.0008 0 2 4 6 8 10 12 14 16 18 20 n j (r ) r (µm) Nb=10000 Np=10000 abf=100 nm ab=2000 nm af=0 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF np, GPDF (b) (a) FIG. 5. �Color online� The normalized densities nj�r�= � j 2 with j=b , p from a solution of Eqs. �40� and �41� for the US, MG- PDF, and GPDF models with Nb=10 000, Np=10 000, abf =100 nm, af =0, and �a� ab=200 nm, and �b� ab=2000 nm. In this plot �nj�r�d3r=1. SUPERFLUID BOSE-FERMI MIXTURE FROM WEAK … PHYSICAL REVIEW A 78, 043616 �2008� 043616-7 GPDF models. For small af the fermionic densities of both MGPDF and GPDF models are a good approximation to the US model. But as af increases, the GPDF model continues with BCS nonlinearity for af =0, whereas the fermionic non- linearity of the MGPDF model acquires excessive attraction. This is clear from Figs. 6�a� and 6�b� showing a peaked fermionic density for the MGPDF model for large af corre- sponding to increased attraction. For very large af , the MG- PDF model does not provide the saturation of fermionic bulk chemical potential as in the US model. We note that com- pared to the study of Figs. 5 with comparable number of bosons and fermions, in the study of Figs. 6 we have a much reduced number of fermions compared to bosons and the system has undergone a mixing-demixing transition �20,66� expelling the fermions from the central to the peripheral re- gion. Next we calculate the chemical potentials �b0 and �p0 of the stationary superfluid Bose-Fermi mixture as described by Eqs. �33� and �34�. In Figs. 7�a� and 7�b� we plot these chemical potentials under the situations described in Figs. 4 and 5, e.g., in the �a� weak-coupling BCS limit when ab varies from a small to a large value and �b� weak-coupling GP limit when af varies from a small to a large value, respectively. In the first case the fermionic chemical potential �p0 remains practically constant over the full range of varia- tion of ab and we plot only the bosonic chemical potential �b0 in Fig. 7�a�. In the second case bosonic chemical poten- tial �b0 remains practically constant over the full range of variation of af and we plot only the fermionic chemical po- tential �p0. From Fig. 7�a� we find that, although for small values of ab, the MGPDF model produces results for �b0 in close agreement with the US model, for medium values of ab the GPDF model produces results closer to the US model. However, for larger values of ab the US model should show saturation, whereas the GPDF chemical potential �b0 should increase monotonically with ab and should be larger than the chemical potential of the US model �not shown in figure�. From Fig. 7�b� we find that none of the MGPDF or GPDF models produces a reasonable result for �p0 except for very small values of af . The GPDF model does not have a de- pendence on af and hence produces a constant result for �p0. V. CONCLUSION The dynamics of cold trapped atoms �both bosons and fermions� should be correctly handled by a proper treatment of many-body dynamics or through a field-theoretic formu- lation. However, both approaches could be much too com- plicated to have an advantage in phenomenological applica- tion. This is why mean-field approaches are widely used for phenomenological application. In the weak-coupling limit, for bosonic atoms the dynamics is handled through the mean-field GP equation �1�, whereas for fermionic atoms the 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-100 nm nb, US 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-100 nm nb, US np, US 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-100 nm nb, US np, US nb, MGPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-100 nm nb, US np, US nb, MGPDF np, MGPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-100 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-100 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF np, GPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-300 nm nb, US 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-300 nm nb, US np, US 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-300 nm nb, US np, US nb, MGPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-300 nm nb, US np, US nb, MGPDF np, MGPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-300 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF 0 0.001 0.002 0.003 0.004 0 2 4 6 8 n j (r ) r (µm) Nb=100000 Np=100 abf=100 nm ab=5 nm af=-300 nm nb, US np, US nb, MGPDF np, MGPDF nb, GPDF np, GPDF (b) (a) FIG. 6. �Color online� The normalized densities nj�r�= � j 2 with j=b , p from a solution of Eqs. �40� and �41� for the US, MG- PDF, and GPDF models with Nb=100 000, Np=100, abf =100 nm, ab=5 nm, and �a� af =−100 nm, and �b� abf =−300 nm. In this plot �nj�r�d3r=1. 0 20 40 60 80 0 200 400 600 800 1000 µ b 0 ab (nm) Nb=10000 Np=10000 abf=100 nm af=0 nm US 0 20 40 60 80 0 200 400 600 800 1000 µ b 0 ab (nm) Nb=10000 Np=10000 abf=100 nm af=0 nm US MGPDF 0 20 40 60 80 0 200 400 600 800 1000 µ b 0 ab (nm) Nb=10000 Np=10000 abf=100 nm af=0 nm US MGPDF GPDF 50 52 54 56 -1000 -800 -600 -400 -200 0 µ p 0 af (nm) Nb=100000 Np=100 abf=100 nm ab=5 nm US 50 52 54 56 -1000 -800 -600 -400 -200 0 µ p 0 af (nm) Nb=100000 Np=100 abf=100 nm ab=5 nm US MGPDF 50 52 54 56 -1000 -800 -600 -400 -200 0 µ p 0 af (nm) Nb=100000 Np=100 abf=100 nm ab=5 nm US MGPDF GPDF (b) (a) FIG. 7. �Color online� �a� Chemical potential �b0 vs ab of Eqs. �33� and �34� for the US, MGPDF, and GPDF models for Nb=Np =10 000, af =0, abf =100 nm and �b� chemical potential �p0 vs af for the US, MGPDF, and GPDF models for Nb=100 000, Np=100, ab=5 nm, abf =100 nm, respectively. The chemical potentials are expressed in oscillator unit ��b. S. K. ADHIKARI AND LUCA SALASNICH PHYSICAL REVIEW A 78, 043616 �2008� 043616-8 LDA �local density approximation� formulation is widely used �2�. The LDA approach �often called the Thomas-Fermi approximation� lacks the kinetic energy gradient term �quan- tum pressure term�. Here we introduce a quantum pressure term consistent with the correct phase-velocity relation of the superfluid fermions assuring the Galilei invariance of the model �23�. This modified LDA equation for fermions and the GP equation for bosons, both valid in the weak-coupling limit, are both shown to be equivalent to the hydrodynamical equation for superfluid flow. In the strong-coupling unitarity limit as the bosonic and fermionic scattering lengths ab and af tend to infinity, both the bosonic and fermionic interac- tions are to exhibit saturation due to constraints of quantum mechanics. We introduce proper unitarity saturation effects in these equations. In addition in our Bose-Fermi formulation we introduce a not too strong contact interaction between bosons and fermions. The resultant superfluid Bose-Fermi dynamics is described by a set of coupled partial differential equations, which is the principal result of this paper. In the time-dependent form with a Bose-Fermi interaction this set is described by Eqs. �40� and �41� with the stationary version without Bose-Fermi interaction given by Eqs. �33� and �34�. The present model for the superfluid Bose-Fermi mixture is termed the US model, valid in both weak- and strong- coupling unitarity limits for bosons and fermions. We also consider a model called the GPDF model where the boson dynamics is handled by the GP Lagrangian and the fermion dynamics is handled by the LDA formulation with a proper kinetic energy gradient term. Finally, we also consider an- other model called the MGPDF model including lowest or- der correction to the GPDF model due to the nonzero bosonic and fermionic scattering lengths ab and af . How- ever, the MGPDF model does not provide a saturation of the nonlinear interaction in the unitarity limit when ab and af tend to infinity. In our numerical studies of a 87Rb-40K mix- ture in a spherical harmonic trap, we find that for medium interactions the MGPDF model produces improved results. However, for stronger interactions the full US model should be used. Although slightly complicated for analytic applica- tions, the present US model is no more complicated than the usual GP and LDA equations for numerical treatment and is expected to find wide applications in future studies of Bose- Fermi superfluid mixture. ACKNOWLEDGMENTS We thank D. Blume for kindly providing additional data of FNMC calculation �53�. S.K.A. has been partially sup- ported by the FAPESP and CNPq of Brazil. L.S. has been partially supported by the Fondazione CARIPARO and GNFM-INdAM. �1� L. Pitaevskii and S. Stringari, Bose-Einstein Condensation �Oxford University Press, Oxford, 2003�; F. Dalfovo, S. Giorgini, L. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 �1999�; C. J. Pethick and H. Smith, Bose-Einstein Conden- sation in Dilute Gases �Cambridge University Press, Cam- bridge, England, 2002�. �2� S. Giorgini, L. P. Pitaevskii, and S. Stringari, e-print arXiv:0706.3360, Rev. Mod. Phys. �to be published�. �3� L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Course of Theoretical Physics No. 6 �Pergamon Press, London, 1987�. �4� L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 2: Theory of the Condensed State, Course of Theoretical Physics No. 9 �Pergamon Press, London, 1987�. �5� A. J. Leggett, Quantum Liquids �Oxford University Press, Ox- ford, 2006�. �6� K. Mölmer, Phys. Rev. Lett. 80, 1804 �1998�. �7� N. Nygaard and K. Mølmer, Phys. Rev. A 59, 2974 �1999�. �8� M. J. Bijlsma, B. A. Heringa, and H. T. C. Stoof, Phys. Rev. A 61, 053601 �2000�. �9� H. Heiselberg, C. J. Pethick, H. Smith, and L. Viverit, Phys. Rev. Lett. 85, 2418 �2000�. �10� L. Viverit, C. J. Pethick, and H. Smith, Phys. Rev. A 61, 053605 �2000�. �11� L. Viverit, Phys. Rev. A 66, 023605 �2002�. �12� K. K. Das, Phys. Rev. Lett. 90, 170403 �2003�. �13� X. J. Liu, M. Modugno, and H. Hu, Phys. Rev. A 68, 053605 �2003�; S. T. Chui and V. N. Ryzhov, ibid. 69, 043607 �2004�; S. K. Adhikari, ibid. 70, 043617 �2004�. �14� D. W. Wang, Phys. Rev. Lett. 96, 140404 �2006�. �15� G. Roati, F. Riboli, G. Modugno, and M. Inguscio, Phys. Rev. Lett. 89, 150403 �2002�. �16� M. Modugno, F. Ferlaino, F. Riboli, G. Roati, G. Modugno, and M. Inguscio, Phys. Rev. A 68, 043626 �2003�. �17� C. Ospelkaus, S. Ospelkaus, K. Sengstock, and K. Bongs, Phys. Rev. Lett. 96, 020401 �2006�. �18� L. Salasnich and F. Toigo, Phys. Rev. A 75, 013623 �2007�. �19� L. Salasnich, S. K. Adhikari, and F. Toigo, Phys. Rev. A 75, 023616 �2007�. �20� S. K. Adhikari and L. Salasnich, Phys. Rev. A 75, 053603 �2007�. �21� K. M. O’Hara, S. L. Hemmer, S. R. Granade, M. E. Gehm, J. E. Thomas, V. Venturi, E. Tiesinga, and C. J. Williams, Phys. Rev. A 66, 041401�R� �2002�; K. Dieckmann, C. A. Stan, S. Gupta, Z. Hadzibabic, C. H. Schunck, and W. Ketterle, Phys. Rev. Lett. 89, 203201 �2002�; C. A. Regal, M. Greiner, and D. S. Jin, ibid. 92, 083201 �2004�; S. Inouye, M. R. Andrews, J. Stenger, H. J. Miesner, D. M. Stamper-Kurn, and W. Ketterle, Nature �London� 392, 151 �1998�. �22� S. K. Adhikari and L. Salasnich, Phys. Rev. A 77, 033618 �2008�. �23� L. Salasnich, e-print arXiv:0804.1277. �24� A. Fabrocini and A. Polls, Phys. Rev. A 60, 2319 �1999�; 64, 063610 �2001�. �25� E. P. Gross, Nuovo Cimento 20, 454 �1961�; L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 �1961� �Sov. Phys. JETP 13, 451 �1961��. �26� P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 �1964�; W. Kohn, Rev. Mod. Phys. 71, 1253 �1999�; R. M. Dreizler and SUPERFLUID BOSE-FERMI MIXTURE FROM WEAK … PHYSICAL REVIEW A 78, 043616 �2008� 043616-9 E. K. U. Gross, Density Functional Theory: An Approach to the Quantum Many-Body Problem �Springer, Berlin, 1990�. �27� W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 �1965�. �28� Y. E. Kim and A. L. Zubarev, Phys. Lett. A 327, 397 �2004�; Phys. Rev. A 70, 033612 �2004�. �29� N. Manini and L. Salasnich, Phys. Rev. A 71, 033625 �2005�; G. Diana, N. Manini, and L. Salasnich, ibid. 73, 065601 �2006�; L. Salasnich and N. Manini, Laser Phys. 17, 169 �2007�. �30� S. K. Adhikari and B. A. Malomed, Europhys. Lett. 79, 50003 �2007�; Phys. Rev. A 76, 043626 �2007�. �31� L. Salasnich, N. Manini, and F. Toigo, Phys. Rev. A 77, 043609 �2008�. �32� L. N. Oliveira, E. K. U. Gross, and W. Kohn, Phys. Rev. Lett. 60, 2430 �1988�. �33� A. Bulgac, Phys. Rev. A 76, 040502�R� �2007�. �34� D. M. Eagles, Phys. Rev. 186, 456 �1969�; A. J. Leggett, J. Phys. �Paris�, Colloq. 41, C7-19 �1980�; P. Nozières and S. Schmitt-Rink, J. Low Temp. Phys. 59, 195 �1985�. �35� L. Salasnich, J. Math. Phys. 41, 8016 �2000�. �36� T. D. Lee and C. N. Yang, Phys. Rev. 105, 1119 �1957�. �37� T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106, 1135 �1957�. �38� S. Cowell, H. Heiselberg, I. E. Mazets, J. Morales, V. R. Pan- dharipande, and C. J. Pethick, Phys. Rev. Lett. 88, 210403 �2002�. �39� S. K. Adhikari, Phys. Rev. A 77, 045602 �2008�. �40� S. Y. Chang, V. R. Pandharipande, J. Carlson, and K. E. Schmidt, Phys. Rev. A 70, 043602 �2004�; J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 �2003�. �41� G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett. 93, 200404 �2004�; 95, 230405 �2005�. �42� H. Heiselberg, Phys. Rev. A 63, 043606 �2001�; G. A. Baker, Jr., Phys. Rev. C 60, 054311 �1999�. �43� V. M. Galitskii, Zh. Eksp. Teor. Fiz. 34, 151 �1958� �Sov. Phys. JETP 7, 104 �1958��. �44� R. Combescot, M. Yu. Kagan, and S. Stringari, Phys. Rev. A 74, 042717 �2006�. �45� S. E. Koonin and D. C. Meredith, Computational Physics For- tran Version �Addison-Wesley, Reading, MA, 1990�. �46� E. Cerboneschi, R. Mannella, E. Arimondo, and L. Salasnich, Phys. Lett. A 249, 495 �1998�; L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 64, 023601 �2001�. �47� S. K. Adhikari and P. Muruganandam, J. Phys. B 35, 2831 �2002�. �48� D. Blume and C. H. Greene, Phys. Rev. A 63, 063601 �2001�. �49� L. Salasnich, N. Manini, and A. Parola, Phys. Rev. A 72, 023621 �2005�; L. Salasnich, ibid. 76, 015601 �2007�. �50� D. Blume, J. von Stecher, and C. H. Greene, Phys. Rev. Lett. 99, 233201 �2007�; J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. A 77, 043619 �2008�. �51� D. Blume, Phys. Rev. A 78, 013613 �2008�; 78, 013635 �2008�. �52� J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. A 76, 053613 �2007�. �53� D. Blume �private communication� kindly provided the data of FNMC calculation of Fig. 4�a�. �54� S. Y. Chang and G. F. Bertsch, Phys. Rev. A 76, 021603�R� �2007�. �55� S. K. Adhikari and L. Salasnich �unpublished�. �56� D. T. Son and M. Wingate, Ann. Phys. �N.Y.� 321, 197 �2006�. �57� N. H. March and M. P. Tosi, Ann. Phys. �N.Y.� 81, 414 �1973�; P. Vignolo, A. Minguzzi, and M. P. Tosi, Phys. Rev. Lett. 85, 2850 �2000�. �58� C. F. von Weizsäcker, Z. Phys. 96, 431 �1935�. �59� E. Zaremba and H. C. Tso, Phys. Rev. B 49, 8147 �1994�. �60� L. Salasnich, J. Phys. A 40, 9987 �2007�. �61� S. K. Adhikari, Phys. Rev. A 72, 053608 �2005�; 70, 043617 �2004�. �62� G. Rupak and T. Schäfer, e-print arXiv:0804.26782v2. �63� L. D. Landau, J. Phys. �Moscow� 5, 71 �1941�; I. M. Khalat- nikov, Pis’ma Zh. Eksp. Teor. Fiz. 17, 534 �1973�; �JETP Lett. 17, 386 �1973��; I. M. Khalatnikov, An Introduction to the Theory of Superfluidity �W. A. Benjamin, New York, 1965�. �64� A. F. Andreev and E. P. Bashkin, Zh. Eksp. Teor. Fiz. 69, 319 �1975� �Sov. Phys. JETP 42, 164 �1975��; T. L. Ho and V. B. Shenoy, J. Low Temp. Phys. 111, 937 �1998�. �65� M. Yu. Kagan, I. V. Brodsky, D. V. Efremov, and A. V. Klaptsov, Phys. Rev. A 70, 023607 �2004�. �66� S. K. Adhikari and L. Salasnich, Phys. Rev. A 76, 023612 �2007�. S. K. ADHIKARI AND LUCA SALASNICH PHYSICAL REVIEW A 78, 043616 �2008� 043616-10