Self-trapping of a Fermi superfluid in a double-well potential in the Bose-Einstein-condensate–unitarity crossover S. K. Adhikari,1,* Hong Lu,2 and Han Pu1,2,† 1Instituto de Física Teórica, São Paulo State University (UNESP), 01.140-070 São Paulo, SP, Brazil 2Department of Physics and Astronomy, and Rice Quantum Institute, Rice University, Houston, Texas 77251, USA �Received 27 July 2009; revised manuscript received 13 October 2009; published 3 December 2009� We derive a generalized Gross-Pitaevskii density-functional equation appropriate to study the Bose-Einstein condensate �BEC� of dimers formed of singlet spin-half Fermi pairs in the BEC-unitarity crossover while the dimer-dimer scattering length a changes from 0 to �. Using an effective one-dimensional form of this equation, we study the phenomenon of dynamical self-trapping of a cigar-shaped Fermi superfluid in the entire BEC- unitarity crossover in a double-well potential. A simple two-mode model is constructed to provide analytical insights. We also discuss the consequence of our study on the self-trapping of an atomic BEC in a double-well potential. DOI: 10.1103/PhysRevA.80.063607 PACS number�s�: 03.75.Ss, 71.10.Ay, 67.85.Lm, 03.75.Lm I. INTRODUCTION After the experimental realization of Bose-Einstein con- densate �BEC� and its controlled study under different trap- ping conditions �1�, there have been many interesting experi- ments with a cigar-shaped BEC in a quasi-one-dimensional �1D� trap with a tight transverse confinement �2�. Along the axial direction several different types of traps have been em- ployed: harmonic �1�, double-well �3�, periodic optical- lattice �4�, and bichromatic optical-lattice �5� traps. Many phenomena have been predicted and observed in such quasi-1D setting. Of these, the ones worth mentioning in- clude the formation of bright �6�, gap �7�, and dark �8� soli- tons, self-trapping �3,9� and Josephson oscillation �3,10,11�. Macroscopic dynamical self-trapping and Josephson os- cillation were predicted theoretically �9,12–15� and observed experimentally �3,10�. Josephson effect was observed in su- perfluid �SF� 3He �16� and 4He �17�. After the experimental observation of BEC in a optical-lattice trap �10�, controlled studies of Josephson oscillation and self-trapping in a cigar- shaped BEC seems well under control �3�. The studies of such phenomena in a cigar-shaped BEC usually employ a double-well potential �3�. In the simplest case of such a sym- metric 1D potential with the origin of the axial coordinate x set at the trap center, under certain initial conditions, when a BEC is released with a population imbalance between two sides of x=0, it executes undamped Josephson oscillation on both sides of the trap center maintaining a time-averaged population imbalance equal to zero. Under different initial conditions, the BEC exhibits self-trapping, occupying pref- erably one side of the trap, thus maintaining a definite non- zero value of time-averaged population imbalance. The un- derstanding of the transition from Josephson oscillation to self-trapping and vice versa has been the topic of many re- cent investigations �3�. A SF Fermi gas in a double-well potential is perhaps even more interesting, nevertheless much less studied �18�. �There have been studies of Josephson oscillation of a Fermi gas in a optical lattice �OL� potential �19�.� Such a trapped SF Fermi gas gives us the unique opportunity to study the Bardeen-Cooper-Schrieffer �BCS� to BEC crossover in a two-component Fermi gas under an entirely different setup. The BCS-BEC crossover can be realized by varying the at- traction between the spin-half fermions forming pairs using the Feshbach resonance technique. As the attraction is in- creased from zero, the simple BCS SF turns into a complex Cooper-pair-induced strongly interacting SF and at unitarity �when the Fermi-Fermi scattering length af tends to infinity�, it is possible for the Cooper pairs to turn spontaneously into Fermi dimers �two-body bound state of fermionic atoms� and the BCS SF turns into a BEC of dimers. After the experimental realization �20� of the BCS-BEC crossover in a trapped Fermi SF by varying the atomic inter- action near a Feshbach resonance, there have been renewed interests �21–23� in the study of a Fermi SF at unitarity and beyond in the BEC region where we have the BEC of dimers. One can thus recover the bosonic behavior in the BEC limit of the crossover �when the dimer-dimer scattering length a tends to zero� while expecting new and distinct behavior in the vicinity of the unitarity regime. Moreover, on the experimental front it is easier to realize a controlled BEC-unitarity crossover �BEC side of the BCS-BEC cross- over� of the Fermi SF than the BCS-unitarity crossover �BCS side of the BCS-BEC crossover� as the superfluid transition temperature in the BCS side of the crossover is very low and difficult to achieve. Here we present a unified Galilean-invariant dynamical equation for the study of the BEC-unitarity crossover of a cigar-shaped BEC of dimers formed of Fermi atoms. In the BEC limit of small dimer-dimer scattering length a, the present equation reduces to the usual Gross-Pitaevskii �GP� equation �24� for bosons, and in the unitarity limit it yields a density-functional �DF� equation �25� for fermions. Hence we call this equation a DF GP equation for a Fermi SF valid in the BEC-unitarity crossover. For the study of a cigar- shaped Fermi SF along the BEC-unitarity crossover, we re- duce the present DF GP equation to a quasi-1D form and use this reduced equation to the study of the self-trapping of *adhikari@ift.unesp.br; www.ift.unesp.br/users/adhikari †hpu@rice.edu PHYSICAL REVIEW A 80, 063607 �2009� 1050-2947/2009/80�6�/063607�10� ©2009 The American Physical Society063607-1 http://dx.doi.org/10.1103/PhysRevA.80.063607 Fermi SF in a double-well potential. This analytic develop- ment is presented in Sec. II. This reduced equation also de- scribes an atomic BEC with repulsive atomic interaction in the BEC-unitarity crossover with different numerical val- ue�s� of certain parameter�s�, and hence the results of the present investigation are also applicable to the self-trapping of a repulsive BEC in a double-well potential. The numerical simulation with the time-dependent quasi-1D equation for a cigar-shaped Fermi SF in a symmet- ric double-well potential with an initial population imbalance between two wells reveals interesting features of the Joseph- son oscillation and self-trapping across the BEC-unitarity crossover. In the limit of zero nonlinearity one has ac Joseph- son oscillation. As nonlinearity is increased by increasing the dimer-dimer scattering length or the number of particle, the Josephson oscillation stops and self-trapping emerges for a double-well potential with appropriate parameters. With fur- ther increase in nonlinearity self-trapping is destroyed and the population in the two wells executes irregular oscillation. For very large nonlinearity, however, the regular Josephson oscillation comes back. Nevertheless, for a small number of particles, the critical nonlinearity required for one of these phenomenon may not be attained and that particular phe- nomenon may not be realized. �The nonlinearity actually saturates for a large value of a and hence cannot be arbi- trarily increased by increasing a as one approaches unitarity for a small number of atoms.� These features are discussed in detail in Sec. IV where we present the numerical results. In Sec. III we present a simple analytic two-mode model to understand the essential features of the numerical results reported in Sec. IV and also point out the limitation of the two-mode model. Finally, in Sec. V we present a brief sum- mary and conclusion of the present investigation. II. DF GP EQUATION FOR A FERMI SF IN THE BEC- UNITARITY CROSSOVER At unitarity the following DF equation for trapped SF fermions �26,27� has produced results for energy in close agreement with independent Monte Carlo calculations �28� �− �2 2m �2 + U + ��n� − i � �t ���r,t� = 0, �1� where U is the trapping potential, m is the mass of a dimer �twice the atomic mass�, ��n�=��2n2/3 /m is the bulk chemi- cal potential of dimers with density �of dimers� n= ���2, and �=2�6�2�2/3�. The normalization condition of the DF wave function is ����2d3r=N, where N is the number of dimers. At unitarity the only length scale is n−1/3, and from dimensional argument the chemical potential—and all energies of the trapped SF fermions—have the above universal form �29�. There have been many theoretical �30,31� and experimen- tal �32� investigations which extracted the value of the con- stant � for fermions, and the most accurate value of this constant is given by independent Monte Carlo calculations by two groups �30�: � 0.44, consequently, � 13.37 and we shall use this value of � in the present study. For a trapped atomic BEC, the energy and chemical potential have the same universal form: ��2n2/3 /m �33�, where now mass m and density n refer to bosons. For fundamental bosonic at- oms, microscopic numerical calculation based on Jastrow variational wave function yielded for the constant � a slightly different value �=22.22 �33�. With this modification in the value of �, the present investigation could be applied to the study of self-trapping of an atomic BEC. At unitarity the Fermi pair can stay as a Cooper pair or a dimer and transform into each other without transfer of en- ergy and Eq. �1� can describe both the Cooper pair and dimer phases. Here we interpret Eq. �1� as the equation for dimers. At unitarity the scattering length a of two dimers goes to infinity: a→�. �Actually, at unitarity, the scattering length of two Fermi atoms af →�. Model studies have indicated that a�af �34�. Consequently, at unitarity we take a→�.� Although Eq. �1� describes both the dimer SF and the Cooper-pair induced BCS SF at unitarity, the bulk chemical potential ��n� appearing in this equation should be inter- preted differently. For the BCS SF it originates from the kinetic energy of Fermi atoms put in different quantum or- bitals consistent with the Pauli principle discounted for by the negative attractive energy due to atomic interaction. For the dimer SF it originates solely from the repulsive interaction energy among dimers. In the BCS limit as af →0−, we have the finite nonlinear term ��n�=2EF =2�6�2�2/3�2n2/3 /m �26� in Eq. �1� originating from the ki- netic energy of Fermi atoms with negligible contribution from interatomic attraction. On the other hand, in the BEC limit as a→0+ the nonlinear term for dimers reduces to zero and at unitarity the nonlinear term in Eq. �1� originates from the saturation of repulsive dimer-dimer interaction as a→�. For the Fermi SF of dimers �and also for an atomic BEC� in the BEC-unitarity crossover the following two leading terms of the bulk chemical potential of a dilute uniform gas can be obtained �35� from the expression for energy per par- ticle as obtained by Lee, Huang, and Yang �36� ��n,a� = �4��2an/m��1 + �n1/3a�3/2 + ¯� , �2� where =32 / �3��� and n1/3a is the dimensionless gas pa- rameter. In this expression the scattering length a must be positive �a 0� corresponding to a repulsive interaction. Higher-order terms of expansion �2� has also been consid- ered �37�; the lowest-order term was derived by Lenz �38�. Considering only the lowest-order term in expansion �2�, ap- propriate in the BEC limit as a ,af →0+, the dimers obey the usual GP equation �24� �− �2 2m �2 + U + 4��2a m ���2 − i � �t ���r,t� = 0. �3� Considering the second term in expansion �2�, in the BEC limit, the following modified GP equation for dimers can be written following the suggestion of Fabrocini and Polls �39�: �− �2 2m �2 + U + 4��2a m ���2�1 + a3/2���� − i � �t ���r,t� = 0. �4� Equation �4� provides an adequate correction to the GP Eq. �3� for small a. But as a increases and diverges at unitarity, the nonlinear term should saturate to the finite universal non- ADHIKARI, LU, AND PU PHYSICAL REVIEW A 80, 063607 �2009� 063607-2 linear term ��n� of Eq. �1� and not diverge like the nonlinear terms of the GP Eq. �3� and of the Fabrocini-Polls Eq. �4�. The chemical potential and energy should not diverge at uni- tarity as the interaction potential remains finite in this limit although the scattering length a diverges. In the weak- coupling GP limit, the scattering length serves as a faithful measure of interaction. But as a increases, it ceases to be a measure of interaction. For a general scattering length, an exact expression of the chemical potential is not available. However, a recent quantum Monte Carlo study maps out the equation of state in the entire BEC-BCS crossover regime �30�. Following a recent suggestion �40�, for the full BEC- unitarity crossover we consider the DF GP equation for the dimer SF providing a smooth interpolation between Eqs. �1� and �4�: �− �2�2 2m + U + ��n,a� − i � �t ���r,t� = 0, �5� ��n,a� = 4��2a m ���2 � � 1 + �1 + ��a3/2��� 1 + �a3/2��� + �1 + ��a5/2���5/3 , �6� where n= ���2 and � and are yet unknown constants satis- fying =4� /�. �These equations are also valid for an atomic BEC with m and a representing atomic mass and scattering length, respectively, and with �=22.22 �33�.� By construction, Eq. �6� yields limit �2� for small a; it also has the correct behavior at unitarity. Using a similar expression for ��n ,a� in the BCS-unitarity crossover �41�, the constant � was calculated �42� by requiring that the first derivative of ��n ,a� with respect to �an1/3� be continuous at unitarity. The condition for continuity yields a small value for �� 0.04�. However, we shall take �=0 in this study. This will make further analytical development easier while maintaining the first derivative of ��n ,a� with respect to �an1/3� approxi- mately continuous at unitarity. A set of equations, similar to Eqs. �5� and �6�, for fundamental bosons, and not for com- posite dimers, produced results for energy �26� of a trapped condensate in agreement with Monte Carlo calculations �43�. A similar equation for the BCS-unitarity crossover produced results for energy �26,41� of a trapped BCS SF in agreement with Monte Carlo calculations �44�. Different parametriza- tion of the chemical potential in the BEC-BCS crossover have been proposed in the literature �45�. We do not expect our results in this work will be sensitive to which specific form for ��n ,a� we choose to use here. Furthermore, it has been shown that the DF GP Eq. �5� is equivalent to the quan- tum hydrodynamic equations for dimers �26,45� �n �t + � · �nv� = 0, m �v �t + ��− �2 2m �2�n �n + mv2 2 + U + ��n,a�� = 0, if we identify ��r,t� = �n�r,t�eis�r,t�, v = � � s/m , where s is a phase and v is the velocity. For a cigar-shaped SF, where the transverse trapping is very strong, the interesting dynamics is confined in the axial direction and in the transverse direction the system is con- fined in its ground state. In such a quasi-1D geometry, the axial and transverse coordinates decouple and it is useful to write an effective 1D equation for the dynamics of a cigar- shaped SF and we perform the same in the following. For the cigar-shaped double-well trap, U�r� = V�x� + m�2�x2 + �2z2 + �2y2�/2, V�x� = A�� exp�− �m�x2/�� , where ��1 and A and � are two dimensionless parameters characterizing the strength and width of the barrier, respec- tively, it is appropriate to take ��r , t�=��x , t���y���z� with ��y�= �m�� / �����1/4exp�−m��y2 / �2��� representing the harmonic-oscillator ground state in the transverse direction and ��x , t� representing the essential dynamics in the x di- rection. The potential V�x� together with the harmonic trap m�2x2 /2 simulate a double well in the axial x direction. Multiplying Eqs. �1� and �4� by ��y���z� and integrating over y and z we get the following 1D equations �27,46�: �− 1 2 �2 �x2 + U�x� + � 3 5 ��n � 2/3 − i � �t ���x,t� = 0, �7� �− 1 2 �2 �x2 + U�x� + 2a�n�1 + a3/24 5 �n� �� − i � �t ���x,t� = 0, �8� where n= ���x , t��2 , U�x�=Ae−�x2 +x2 /2 represents the double-well potential and n= ���2, and we use harmonic os- cillator dimensionless units �=m=�=1. All lengths are now expressed in oscillator unit �� / �m�� and time in �−1 and � is normalized as �−� � dx���x , t��2=N. A simple DF GP equation interpolating between Eqs. �7� and �8� is �− 1 2 �2 �x2 + U�x� + ��n,a� − i � �t ���x,t� = 0, �9� ��n,a� = 2a�n�1 + 4 5 �� �� a3/2�n 1 + �a5/2n5/6 � , �10� where n= ���2 and �=8 �5/6�1/6 / �3��. Equation �9� reduces in the BEC a→0+ limit to Eq. �8� and in the unitarity a→� limit to Eq. �7�. We shall use Eq. �9� for the descrip- tion of self-trapping and Josephson oscillation in the BEC- unitarity crossover of fermions. SELF-TRAPPING OF A FERMI SUPERFLUID IN A … PHYSICAL REVIEW A 80, 063607 �2009� 063607-3 III. TWO-MODE MODEL OF THE FERMI SF A. Two-mode model Before we present the full numerical results, it is instruc- tive to consider the so-called two-mode model �9� which is widely used in the study of BEC in a double-well potential, and more recently, has been used in the investigation of Fermi SF across a weak link �47�. Due to its simplicity, the two-mode model can provide many useful insights. Here we construct the corresponding two-mode model of the Fermi SF based on Eq. �9�. To this end, we decompose ��x , t� as ��x,t� = �1�t��1�x� + �2�t��2�x� , �11� where the spatial mode functions �1,2�x� are assumed to be real, satisfy the orthonormal condition � �i�x�� j�x�dx = �ij , and are localized in each of the two wells, respectively. �1,2�t� are in general complex and satisfy the condition ��1,2�t��2=N1,2�t� so that ��1�t��2 + ��2�t��2 = N1�t� + N2�t� = N . Inserting decomposition �11� into Eq. �9�, integrating out the spatial degrees of freedom, we obtain the following equa- tions of motion for �1,2�t�: i�̇1 = E1�1 + E1���1���1 − K�2, �12� i�̇2 = E2�1 + E2���2���2 − K�1, �13� where Ei =� dx�i�x��− 1 2 d2 dx2 + U�x� �i�x� , �14� Ei���i�� =� dx�i�x���ni,a��i�x� , �15� K = −� dx�1�x��− 1 2 d2 dx2 + U�x� �2�x� , �16� with ni= ��i�i�2. Here we have neglected integrals involving spatial overlaps of �1�x� and �2�x�. For simplicity, we assume a symmetric double well with U�x�=U�−x� so that �1�x�=�2�−x� and consequently, we have E1=E2 and E1�����=E2������E�����. Let us write the waves �1,2 in terms of its amplitude �N1,2 and phase ��1,2� �1,2 = �N1,2ei�1,2 and define a pair of conjugate variables: S � �N1 − N2�/N, � � �2 − �1. Here the variable S denotes the population imbalance be- tween the two wells and � is the phase difference. After some straightforward algebra, the following equations of motion for S and � can be derived from Eqs. �12� and �13� by equat- ing the real and imaginary parts of both sides: Ṡ = − 2K�1 − S2 sin � , �̇ = �E��N1� − E��N2�� + 2K S �1 − S2 cos � . These are the two-mode equations for the Fermi SF. Note that the two-mode equations describing weakly interacting bosons in double-well potential �9� are recovered if we take ��n ,a�=2a�n and, correspondingly, E��Ni�=2a�Ni���i�4dx. The two-mode equations can be cast into the canonical form Ṡ = − �H �� , �̇ = �H �S , with the classical Hamiltonian defined as H =� �E��N1� − E��N2��dS − 2K�1 − S2 cos � . �17� By studying the properties of this Hamiltonian, we can tell whether the system should exhibit self-trapping or Josephson oscillation. B. Fermi SF at unitarity As a concrete example, let us consider the Fermi SF at unitarity where a→� and from Eq. �10� we find ��n,a� = ��n� = 3� 5 ��n � 2/3 . It follows that E��Ni� = UNi 2/3, with U= �3� /5��� /��2/3���i�10/3dx. The classical Hamil- tonian takes the form H 2K = 3� 5 ��1 + S�5/3 + �1 − S�5/3� − �1 − S2 cos � , where �= �N /2�2/3U / �2K� measures the ratio of the strength of the nonlinearity �N /2�2/3U and the tunneling energy 2K. As can be seen from Fig. 1, if we draw equienergy con- tours of H in the phase space of �S ,��, we can see two types of contours: those that form closed loop and those do not. The division of these two types of contours occurs at the critical energy �we take 2K as the units for energy� Ecrit = H�S = 0,� = �� = 6� 5 + 1. If the system has energy E�Ecrit, its dynamics will follow the closed contours, and both S and � will oscillate in time. In particular, the population imbalance S oscillates around 0 and the time averaged population imbalance vanishes. This corresponds to the ac Josephson regime. On the other hand, if the system has energy E Ecrit, it will follow the open contours where � will grow indefinitely and S will oscillate ADHIKARI, LU, AND PU PHYSICAL REVIEW A 80, 063607 �2009� 063607-4 around a nonzero value and will never cross the S=0 line. This corresponds to the self-trapping regime. Given the initial values S�0� and ��0�, the system moves on a contour of constant energy given by E0 = 3� 5 ��1 + S�0��5/3 + �1 − S�0��5/3� − �1 − S�0�2 cos ��0� . The condition for self-trapping E0 Ecrit = 6� 5 + 1 may be recast into the form � �c = 5 3 1 + �1 − S�0�2 cos ��0� �1 + S�0��5/3 + �1 − S�0��5/3 − 2 . �18� In other words, within the two-mode model, the ratio of the nonlinear strength and the tunneling energy, �, determines whether the system should be self-trapped or not. To determine the values of these quantities, we need to choose properly the spatial mode functions �1,2�x�. A reason- able choice is given by �9,48� �1,2�x� = �+�x� � �−�x� �2 , where the normalized wave functions ���x� are the lowest- energy symmetric and antisymmetric stationary solutions to the time-independent DF equation: �̄��� = �− 1 2 d2 dx2 + U�x� + 3� 5 �N� � 2/3 ����4/3���, with chemical potential ��. In Figs. 2�a�–2�c� we illustrate the nonlinear strength �N /2�2/3U, tunneling energy 2K and their ratio � as func- tions of N�, respectively. As N� increases, the strength of the �repulsive� nonlinearity increases. As a result, the �1 and �2 become widened and enjoy more overlap. This leads to an increased tunneling energy. However, the ratio of the nonlin- ear strength and the tunneling energy � does not have a monotonic behavior as N� is increased. As shown in Fig. 2�c�, � initially increases for small N�, reaches a peak and then decreases. In a recent study, Salasnich et al. �47� used the local-density approximation on top of quantum Monte Carlo data of Ref. �30� to explore the phase diagrams and find regimes of Josephson tunneling and of dynamical self- trapping of a three-dimensional �3D� Fermi superfluid. In the two-mode approach reported in Ref. �47�, a constant tunnel- ing energy is arbitrarily chosen for the whole crossover re- gime. This is an inappropriate oversimplification. In Fig. 2�d�, we show the critical value �c as a function of the initial population imbalance S�0�. One can see that as S�0� in- creases, �c decreases rapidly. The fact that � is bounded from above even though the interaction strength can increase without bound has impor- tant consequences. For example, for certain initial condi- tions, self-trapping may only occur within an intermediate range of nonlinearity. Both too small and too large a nonlin- earity will destroy self-trapping. This statement is also true away from the unitarity even in the BEC limit. Furthermore, for a sufficiently small S�0�, � may never exceed the corre- sponding �c. When this is the case, the system will always stay in the Josephson oscillation regime. For example, ac- cording to Fig. 2�d�, �c 300 for S�0�=0.1. Figure 2�c� shows that the system can therefore never reach the self- trapping regime if S�0� equals 0.1 or smaller. This is consis- tent with our numerical results. We want to remark that even though results obtained from the simple two-mode model may provide significant qualita- tive insights, they are not expected to be accurate quantita- tively. Particularly for large nonlinearity, predictions from -1 -0.5 0 0.5 1 (a) -0.8 -0.4 0 0.4 0.8 1.2Λ = 8 E =12 11.5 11 10.5 10 9.5 9 θ/π S -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S 6 -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S 6.5 -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S 7 -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S 7.5 -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S 8 -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S -1 -0.5 0 0.5 1 E = 10 Λ = 8.5 (b) -0.8 -0.4 0 0.4 0.8 1.2 θ/π S 9 FIG. 1. �Color online� Equienergy phase contour plot of the unitary Fermi SF. �a� Contours of different energies for �=8. �b� Contours of E=10 for different values of �. Energy is in units of 2K. -1 -0.5 0 0.5 1 1.5 2 2.5 3 0 20 40 60 (N /2 )2 /3 U -1 -0.5 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 2 K -1 -0.5 0 0.5 1 1.5 2 2.5 3 0 100 200 300 log (Nλ) Λ (a) (b) (c) lo g (Λ c) 0 0.2 0.4 0.6 0.8 1 4 3 2 1 0 S(0) (d) 1 0 10 FIG. 2. �a� Strength of nonlinearity �N /2�3/2U, �b� tunneling energy 2K, and �c� their ratio � �c� as functions of N� at unitarity. �d� The critical �c as function of S�0� for ��0�=0 �see Eq. �18��. SELF-TRAPPING OF A FERMI SUPERFLUID IN A … PHYSICAL REVIEW A 80, 063607 �2009� 063607-5 the two-mode model can deviate greatly from the numerical results �48,49�. The error mainly occurs in estimating the tunneling energy K. The two-mode Eqs. �12� and �13� are obtained by neglecting many terms involving overlap inte- grals of the mode functions �1 and �2 and hence in general greatly underestimates the tunneling energy, particularly for large nonlinearity when overlap between �1 and �2 can be significant. Furthermore, the two-mode approximation itself becomes questionable for large nonlinearity. When there is exchange of atoms between the two wells, the mode func- tions will change accordingly due to the modification of the nonlinear mean field. Indeed, in our numerical calculations to be presented below, we observe that the spatial wave func- tion of the system changes in time. In certain regimes, this change is significantly enough to invalidate the two-mode model. IV. NUMERICAL RESULTS In this section we present an account of the numerical study of self-trapping and Josephson oscillation in a double- well potential by solving the full quasi-1D DF GP Eq. �9� valid for a cigar-shaped SF. The double-well potential is taken as U�x� = x2/2 + Ae−�x2 . �19� We shall take the parameters A and � of this double well similar to the ones employed in Ref. �49� in a study of self- trapping with dipolar bosonic atoms. To create an initial state with desired population imbal- ance for a given set of parameters N� and a, we search for the ground state of an asymmetric well comprised of an off- centered harmonic potential and the Gaussian barrier poten- tial U��x� = �x − x0�2/2 + Ae−�x2 . �20� The ground state of this asymmetric well is obtained by solv- ing the time-independent version of the DF GP Eq. �9� using the imaginary time evolution method. The parameter x0 in Eq. �20� is chosen so that the population imbalance S�t� = �N1�t� − N2�t��/N , �21� has a fixed predetermined initial value S�0�. Here N1�t� and N2�t� are the number of dimers in the first and the second well of the double-well potential. Experimentally, this is in- deed the method used to generate the initial population im- balance �3�. We have also considered other forms of initial wave functions and found that the final results are qualita- tively insensitive to the specific forms provided the initial population imbalance S�0� is kept fixed at a small value. However, at a quantitative level the results could be sensitive to the form of the initial wave function. The sensitivity of the result to the initial wave form increases as S�0� is increased. Moreover, the results are quite sensitive to the initial S�0� employed. Once the initial wave function is chosen, Eq. �9� is solved numerically after discretization with the Crank-Nicolson scheme �50,51� employing space and time steps 0.025 and 0.0002, respectively, using real-time propagation with the FORTRAN programs provided in Ref. �50�. The results are also independently confirmed using a MATLAB code based on the split fast Fourier transform method. Now we present results of dynamical evolution of a Fermi SF, where we take �=13.37 in Eqs. �9� and �10�. The numeri- cal study of self-trapping and Josephson oscillation with the Fermi SF of dimers along the BEC-unitarity crossover re- veals interesting features. To start the investigation of self- trapping we fix the trap parameters A and � in Eqs. �19� and �20� at nontrivial values �a not too small value of A and a not too large ��, which permit smooth and free Josephson oscil- lation in the BEC limit �a=0�. Note that a very small value of A and a very large value of � tend to reduce the double well �Eq. �19�� to a single well where there cannot be any self-trapping and Josephson oscillation should appear for all values of a and N. In order to have self-trapping, A cannot be too small and � cannot be too large. This is illustrated in Fig. 3�a� for A=8, �=10, N�=100, S�0�=0.2 and for a=0.1, 0.01, and 0.001 where we plot S�t� vs t. There is no self- trapping for a very small value of S�0��=0.1�. The quantity S�t� is experimentally measurable and S�t� vs t dynamics provides information about self-trapping and Josephson os- cillation. From Fig. 3�a� we find that there is Josephson os- cillation for all values of a and there is no sign of self- trapping. �A nonzero time average �S�t�� ensures self- trapping.� But a completely new scenario emerges as A is increased to 16 from 8, as can be seen from Fig. 3�b� where we show the S�t� vs t dynamics for a=0.01, 0.1, 0.2, and 0.5. The plots for a=0.01 and 0.5 of Fig. 3�b� are quite similar to the plots for a=0.001 and 0.1 of Fig. 3�a� illustrating regular �periodic� Josephson oscillation with no sign of self- trapping. But, for intermediate values 0.1 and 0.2 of a, self- trapping and irregular �nonperiodic� oscillation can be seen in Fig. 3�b�. In Figs. 3�c� and 3�d� we illustrate two more cases of S�t� vs t dynamics with a different value of S�0��=0.3� and different N� and trap parameters, respec- tively, where one can clearly find self-trapping. In the following, we discuss in detail the results for three initial population imbalance S�0�=0.1, 0.2, and 0.3, which are representative for a general case. Population imbalance S�0�=0.1: for this relatively small initial population imbalance, we found that for any values of N� and a, the system is always in the Josephson regime: the population imbalance S�t� oscillates sinusoidally between S�0�=−0.1 and S�0�=0.1. The system never exhibits self- trapping. The frequency of oscillation increases as the strength of nonlinearity increases. Note that the strength of nonlinearity is increased by increasing either N� or a. How- ever, the nonlinear interaction among dimers saturates as scattering length a→� at unitarity, it increases indefinitely with N�. This result is consistent with our previous discus- sion of the two-mode model: for a sufficiently small S�0�, the required critical value of �c for self-trapping cannot be achieved by increasing the strength of the nonlinearity and the system stays in the Josephson regime for all values of N� and a. Population imbalance S�0�=0.2: the results for the S�t� vs t dynamics for this initial population imbalance is illustrated in Figs. 3�a� and 3�b� for fixed N�=100 and various values ADHIKARI, LU, AND PU PHYSICAL REVIEW A 80, 063607 �2009� 063607-6 of a and two different traps as described earlier. In Fig. 3�b�, for a small scattering length of a=0.01 �solid line�, Joseph- son oscillation is observed. When a is increased to 0.1 �dashed line�, self-trapping is clearly seen—S�t� does not de- viate from S�0� too much and never crosses zero. With fur- ther increase in a to a slightly larger value of 0.2 �dotted line�, self-trapping is destroyed and S�t� exhibits irregular oscillations around zero. Accompanied with this irregular population oscillation, the density profile ���x , t��2 also de- velops complex and irregular structures. Remarkably, upon further increase in a, as the dot-dashed curve for a=0.5 shows, regular oscillation returns and the population dynam- ics once again exhibits sinusoidal Josephson oscillations. Population imbalance S�0�=0.3: finally, let us discuss this relatively large initial population imbalance. If we use N�=100 as we did above for S�0�=0.2, when a is increased from zero to �, the system sequentially makes transitions from Josephson, to self-trapping and finally to irregular os- cillation regimes. The Josephson oscillation is never recov- ered for N�=100 for very large values of scattering length a �results not shown here�. In Fig. 3�c� we plot the results for S�t� vs t dynamics for N�=1000. In this case, in addition to the three regimes just mentioned, for a sufficiently large a�=10�, Josephson oscillation is restored, just as in the case of S�0�=0.2 and N�=100 discussed above. In Fig. 3�d� we show another example of S�t� vs t dynamics for a different trap, which is quite similar to that in Fig. 3�c�. We also did some calculation with larger S�0� where a similar panorama emerges and we do not report the details here. To summarize the general characteristics of the population dynamics, we find that for any given initial population im- balance and for either sufficiently small or sufficiently large nonlinear interaction strength, the system is in the Josephson oscillation regime. For intermediate interaction strength, the system can make transition to self-trappping and irregular oscillation regimes as schematically shown in Fig. 4. The critical interaction strength at which the system makes the transition to self-trapping is sensitive to the initial population imbalance and increases sharply as S�0� increases. �It is also sensitive to the parameters for the Gaussian barrier that cre- ates the double-well potential.� It is possible that for a suffi- ciently small S�0�, the system always stays in the Josephson regime. The restoration of the Josephson oscillation at large interaction strength may seem surprising at first sight. How- ever, one can understand it in the following intuitive way. For a sufficiently large interaction strength, the chemical po- tential is large and the effect of the Gaussian barrier becomes relatively unimportant. The wave functions on opposite sides of the barrier have sufficient overlap and hence the cloud tunnel back and forth without difficulty. The appearance of self-trapping is best illustrated through a study of the time-averaged population imbalance �S�t�� vs nonlinearity aN� and we do that next. In Fig. 5�a� we plot �S�t�� vs aN� by varying the scattering length from 0 to � -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t S(0) = 0.2 A = 8 κ = 10 (a) Nλ = 100 a = 0.001 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t S(0) = 0.2 A = 8 κ = 10 (a) Nλ = 100 a = 0.001 = 0.01 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t S(0) = 0.2 A = 8 κ = 10 (a) Nλ = 100 a = 0.001 = 0.01 = 0.1 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t A = 16 κ = 10 (b) Nλ = 100 S(0) = 0.2 a = 0.01 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t A = 16 κ = 10 (b) Nλ = 100 S(0) = 0.2 a = 0.01 = 0.1 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t A = 16 κ = 10 (b) Nλ = 100 S(0) = 0.2 a = 0.01 = 0.1 = 0.2 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 S( t) t A = 16 κ = 10 (b) Nλ = 100 S(0) = 0.2 a = 0.01 = 0.1 = 0.2 = 0.5 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3 A = 16 κ = 10Nλ = 1000(c) a = 0.001 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3 A = 16 κ = 10Nλ = 1000(c) a = 0.001 = 0.006 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3 A = 16 κ = 10Nλ = 1000(c) a = 0.001 = 0.006 = 0.02 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3 A = 16 κ = 10Nλ = 1000(c) a = 0.001 = 0.006 = 0.02 = 10 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3(d) A = 12 κ = 8Nλ = 100 a = 0.01 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3(d) A = 12 κ = 8Nλ = 100 a = 0.01 = 0.05 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3(d) A = 12 κ = 8Nλ = 100 a = 0.01 = 0.05 = 0.1 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 60 80 100 S( t) t S(0) = 0.3(d) A = 12 κ = 8Nλ = 100 a = 0.01 = 0.05 = 0.1 = 0.5 FIG. 3. �Color online� Population imbalance S�t� vs t dynamics with �a� N�=100, S�0�=0.2, �=10, A=8, �b� N�=100, S�0� =0.2, �=10, A=16, �c� N�=1000, S�0�=0.3, �=10, A=16, and �d� N�=100, S�0�=0.3, �=8, A=12 for dimer-dimer scat- tering length a varying over the BEC-unitarity crossover. Josephson self-trapping irregular oscillation Josephson increasing nonlinear interaction strength FIG. 4. Different dynamical regimes of Fermi superfluid in the BEC-unitary crossover. SELF-TRAPPING OF A FERMI SUPERFLUID IN A … PHYSICAL REVIEW A 80, 063607 �2009� 063607-7 for a fixed N� with trap parameters A=16 and �=10. The initial population imbalance is chosen as S�0�=0.3. In Fig. 5�a�, with the increase in aN�, self-trapping appears for aN� slightly greater than unity. With further increase in aN�, self- trapping increases with an increase in �S�t��. For N�=10, self-trapping never disappears and continues even at unitar- ity. However, beyond aN� 5 self-trapping decreases with the increase in aN� for N�=20, 100 and 1000. For larger N�, �S�t�� eventually goes to zero as the nonlinear repulsion becomes too large to maintain all dimers in a single trap, except for N�=10. In Fig. 5�b� we exhibit a similar plot for S�0�=0.6 with trap parameters A=16 and �=10. Finally, to check the validity of the quasi-1D approxima- tion, we performed full 3D numerical simulations based on Eqs. �5� and �6�. The quasi-1D approximation should be valid when �����. We have chosen different sets of pa- rameters, some of which satisfy and the rest violate the quasi-1D condition. For parameters such that the quasi-1D condition is satisfied, we indeed find that our 3D numerical results are nearly identical to the 1D results presented here. For parameters that the quasi-1D condition is violated, the 3D results show deviations from the 1D ones. However, the qualitative features �i.e., the dependence of different dynami- cal regimes on the initial population imbalance and the strength of nonlinearity� presented in Figs. 3 and 5 remain valid. Specifically, we verified that the results reported in Figs. 3 remain essentially valid in the 3D model. This assures us of the reliability of the present study in 1D. V. SUMMARY AND CONCLUSION To summarize, we have studied the dynamical properties of a Fermi SF confined in a double-well potential in the BEC-unitary crossover regime. To this purpose, we have de- veloped a nonlinear Schrödinger equation valid in the whole regime based on a density-functional approach and on the equations of state from quantum Monte Carlo calculations. This equation is equivalent to the hydrodynamic equations with the quantum pressure term included. In the BEC side of the crossover, it describes accurately the equilibrium and low-energy dynamical properties of the Fermi SF. In particu- lar, Josephson effect has been investigated using this method �18� and the results have been shown to agree with those obtained from the microscopic approach by solving the Bogoliubov–de Gennes equations �52�. Compared with the latter, the great advantage of the current approach is its math- ematical simplicity. We have identified three dynamical regimes of the sys- tem: the Josephson regime, the self-trapping regime and the irregular oscillation regime. For a given initial population imbalance, these regimes are accessed according to the strength of nonlinearity as schematically shown in Fig. 4. The Josephson regime is always reached at either sufficiently small or sufficiently large interaction strength. For a small initial population imbalance, Josephson regime may be the only regime that the system can have access to. Note that the strength of nonlinearity can be increased by either increasing the number of dimers N or the scattering length a. However, it saturates as a tends to infinity while no saturation occurs for large N. The quasi-1D model �Eqs. �9� and �10�� presented and used in the study of dynamical evolution of a Fermi SF in the BEC-unitarity crossover in this paper is also valid for an atomic BEC with a slightly modified value for the parameter �. Hence the present results for self-trapping of a Fermi SF in a double-well potential are also applicable for a repulsive atomic BEC when the atomic scattering length varies from 0 to �. However, in this case there could be practical difficulty with three-body loss in the experimental realization of the system for large scattering length. We have also developed a simple analytical two-mode model, analogous to the much studied system of a BEC in a double-well potential. We show that the properties of the system can be described by a classical Hamiltonian with population imbalance and relative phase as a pair of conju- gate variables. The great advantage of the two-mode model is its simplicity which makes analytical studies possible. The key parameters that characterize the two-mode model are the strength of nonlinearity and the tunneling energy. We calcu- lated these parameters using the spatial mode function ob- tained by numerically solving the full time-independent non- linear Schrödinger equation. From this calculation we show that the ratio of the interaction strength and the tunneling rate cannot increase indefinitely when the interaction strength in- 0 0.1 0.2 0.3 100 101 102 < S( t) > aNλ S(0) = 0.3 Nλ = 10 0 0.1 0.2 0.3 100 101 102 < S( t) > aNλ S(0) = 0.3 Nλ = 10 = 20 0 0.1 0.2 0.3 100 101 102 < S( t) > aNλ S(0) = 0.3 Nλ = 10 = 20 = 100 0 0.1 0.2 0.3 100 101 102 < S( t) > aNλ S(0) = 0.3 Nλ = 10 = 20 = 100 = 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 10-1 100 101 102 103 104 < S( t) > aNλ S(0) = 0.6 Nλ = 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 10-1 100 101 102 103 104 < S( t) > aNλ S(0) = 0.6 Nλ = 10 Nλ = 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 10-1 100 101 102 103 104 < S( t) > aNλ S(0) = 0.6 Nλ = 10 Nλ = 100 Nλ = 1000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 10-1 100 101 102 103 104 < S( t) > aNλ S(0) = 0.6 FIG. 5. �Color online� Time-averaged population imbalance �S�t�� vs nonlinearity aN� for different N� for trap parameters A=16, �=10 and �a� S�0�=0.3 and �b� S�0�=0.6. Different curves are generated by varying the scattering length a across the BEC- unitarity crossover for S�0�=0.3 and 0.6, respectively, based on initial wave form �20�. ADHIKARI, LU, AND PU PHYSICAL REVIEW A 80, 063607 �2009� 063607-8 creases. This explains the numerical observation that for suf- ficiently small initial population imbalance, the system may always stay in the Josephson regime. However, care must be taken when making quantitative comparisons with numerical results. In particular, for strong nonlinearity, the two-mode model can be even qualitative incorrect. For example, this model predicts the existence of the Josephson and the self- trapping regime but not the irregular oscillation regime found in the numerical calculation, which occurs at relatively large nonlinearity and lies in the regime where the two-mode model is no longer valid. ACKNOWLEDGMENTS FAPESP and CNPq �Brazil� provided partial support. H.P. acknowledges support from NSF and the Robert A. Welch Foundation �Grant No. C-1669�. �1� L. Pitaevskii and S. Stringari, Bose-Einstein Condensation �Oxford University Press, Oxford, 2003�; C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases �Cam- bridge University Press, Cambridge, England, 2002�. �2� A. Görlitz et al., Phys. Rev. Lett. 87, 130402 �2001�. �3� M. Albiez, R. Gati, J. Folling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 �2005�; R. Gati and M. K. Oberthaler, J. Phys. B 40, R61 �2007�. �4� A. Kastberg, W. D. Phillips, S. L. Rolston, R. J. C. Spreeuw, and P. S. Jessen, Phys. Rev. Lett. 74, 1542 �1995�. �5� G. Roati et al., Nature �London� 453, 895 �2008�; S. K. Adhikari and L. Salasnich, Phys. Rev. A 80, 023606 �2009�. �6� K. E. Strecker et al., Nature �London� 417, 150 �2002�; L. Khaykovich et al., Science 296, 1290 �2002�. �7� B. Eiermann, T. Anker, M. Albiez, M. Taglieber, P. Treutlein, K. P. Marzlin, and M. K. Oberthaler, Phys. Rev. Lett. 92, 230401 �2004�; O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 �2006�. �8� B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark, and E. A. Cornell, Phys. Rev. Lett. 86, 2926 �2001�. �9� A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Phys. Rev. Lett. 79, 4950 �1997�; S. Raghavan, A. Smerzi, S. Fan- toni, and S. R. Shenoy, Phys. Rev. A 59, 620 �1999�. �10� F. S. Cataliotti et al., Science 293, 843 �2001�. �11� S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Nature �London� 449, 579 �2007�. �12� L. Morales-Molina and J. B. Gong, Phys. Rev. A 78, 041403�R� �2008�. �13� E. A. Ostrovskaya, Y. S. Kivshar, M. Lisak, B. Hall, F. Cattani, and D. Anderson, Phys. Rev. A 61, 031601�R� �2000�. �14� A. Trombettoni and A. Smerzi, Phys. Rev. Lett. 86, 2353 �2001�. �15� J. Javanainen, Phys. Rev. Lett. 57, 3164 �1986�; J. E. Will- iams, Phys. Rev. A 64, 013610 �2001�; S. Giovanazzi, A. Smerzi, and S. Fantoni, Phys. Rev. Lett. 84, 4521 �2000�; I. Zapata, F. Sols, and A. J. Leggett, Phys. Rev. A 57, R28 �1998�; S. K. Adhikari, ibid. 72, 013619 �2005�; Eur. Phys. J. D 25, 161 �2003�. �16� S. Pereverzev et al., Nature �London� 388, 449 �1997�; S. Backhaus et al., Science 278, 1435 �1997�. �17� K. Sukhatme et al., Nature �London� 411, 280 �2001�. �18� F. Ancilotto, L. Salasnich, and F. Toigo, Phys. Rev. A 79, 033627 �2009�; L. Salasnich, F. Ancilotto, N. Manini, and F. Toigo, Laser Phys. 19, 636 �2009�. �19� L. Pezzè, L. Pitaevskii, A. Smerzi, S. Stringari, G. Modugno, E. De Mirandes, F. Ferlaino, H. Ott, G. Roati, and M. Inguscio, Phys. Rev. Lett. 93, 120401 �2004�; S. K. Adhikari, Eur. Phys. J. D 47, 413 �2008�. �20� M. Greiner et al., Nature �London� 426, 537 �2003�; C. Chin et al., Science 305, 1128 �2004�; C. A. Regal, M. Greiner, and D. S. Jin, Phys. Rev. Lett. 92, 040403 �2004�; J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and J. E. Thomas, ibid. 92, 150402 �2004�; M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Raupach, A. J. Kerman, and W. Ketterle, ibid. 92, 120403 �2004�; M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. H. Denschlag, and R. Grimm, ibid. 92, 203201 �2004�; M. W. Zwierlein et al., Nature �Lon- don� 435, 1047 �2005�. �21� J. Kinast et al., Science 307, 1296 �2005�; T. Bourdel, J. Cu- bizolles, L. Khaykovich, K. M. F. Magalhaes, S. J. J. M. F. Kokkelmans, G. V. Shlyapnikov, and C. Salomon, Phys. Rev. Lett. 91, 020402 �2003�. �22� S. Giorgini et al., Rev. Mod. Phys. 80, 1215 �2008�. �23� A. Bulgac and G. F. Bertsch, Phys. Rev. Lett. 94, 070401 �2005�; S. Stringari, ibid. 102, 110406 �2009�; G. M. Bruun, A. Recati, C. J. Pethick, H. Smith, and S. Stringari, ibid. 100, 240406 �2008�. �24� E. P. Gross, Nuovo Cimento 20, 454 �1961�; L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 �1961� �Sov. Phys. JETP 13, 451 �1961��. �25� R. M. Dreizler and E. K. U. Gross, Density Funtional Theory �Springer-Verlag, Berlin, 1990�. �26� S. K. Adhikari and L. Salasnich, Phys. Rev. A 78, 043616 �2008�. �27� S. K. Adhikari and L. Salasnich, New J. Phys. 11, 023011 �2009�. �28� 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�; S. Y. Chang and G. F. Bertsch, ibid. 76, 021603�R� �2007�. �29� H. Heiselberg, Phys. Rev. A 63, 043606 �2001�; G. A. Baker, Jr., Phys. Rev. C 60, 054311 �1999�; Int. J. Mod. Phys. B 15, 1314 �2001�. �30� G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett. 93, 200404 �2004�; J. Carlson, S. Y. Chang, V. R. Pandharipande, and K. E. Schmidt, ibid. 91, 050401 �2003�. �31� J. R. Engelbrecht, M. Randeria, and C. A. R. Sa de Melo, Phys. Rev. B 55, 15153 �1997�; S. Y. Chang, V. R. Pandharipande, J. Carlson, and K. E. Schmidt, Phys. Rev. A 70, 043602 �2004�. �32� M. E. Gehm, S. L. Hemmer, S. R. Granade, K. M. OHara, and SELF-TRAPPING OF A FERMI SUPERFLUID IN A … PHYSICAL REVIEW A 80, 063607 �2009� 063607-9 J. E. Thomas, Phys. Rev. A 68, 011401�R� �2003�; M. Barten- stein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. H. Den- schlag, and R. Grimm, Phys. Rev. Lett. 92, 120401 �2004�; G. B. Partridge, J. P. Gaebler, C. A. Regal, and D. S. Jin, Science 311, 503 �2006�; J. T. Stewart et al., Phys. Rev. Lett. 97, 220406 �2006�. �33� S. Cowell, H. Heiselberg, I. E. Mazets, J. Morales, V. R. Pan- dharipande, and C. J. Pethick, Phys. Rev. Lett. 88, 210403 �2002�. �34� D. S. Petrov, C. Salomon, and G. V. Shlyapnikov, Phys. Rev. Lett. 93, 090404 �2004�. �35� J. K. Nilsen, J. Mur-Petit, M. Guilleumas, M. Hjorth-Jensen, and A. Polls, Phys. Rev. A 71, 053610 �2005�. �36� T. D. Lee and C. N. Yang, Phys. Rev. 105, 1119 �1957�; T. D. Lee, K. Huang, and C. N. Yang, ibid. 106, 1135 �1957�. �37� T. T. Wu, Phys. Rev. 115, 1390 �1959�; E. Braaten and A. Nieto, Eur. Phys. J. B 11, 143 �1999�. �38� W. Lenz, Z. Phys. 56, 778 �1929�. �39� A. Fabrocini and A. Polls, Phys. Rev. A 60, 2319 �1999�; 64, 063610 �2001�. �40� S. K. Adhikari and L. Salasnich, Phys. Rev. A 77, 033618 �2008�. �41� S. K. Adhikari, Phys. Rev. A 79, 023611 �2009�; Laser Phys. Lett. 6, 901 �2009�. �42� S. K. Adhikari �unpublished�. �43� D. Blume and C. H. Greene, Phys. Rev. A 63, 063601 �2001�. �44� J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. A 76, 053613 �2007�. �45� N. Manini and L. Salasnich, Phys. Rev. A 71, 033625 �2005�; Y. E. Kim and A. L. Zubarev, ibid. 70, 033612 �2004�; S. K. Adhikari, ibid. 77, 045602 �2008�. �46� L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 043614 �2002�; A. Muñoz Mateo and V. Delgado, ibid. 77, 013617 �2008�; C. A. G. Buitrago and S. K. Adhikari, J. Phys. B 42, 215306 �2009�. �47� L. Salasnich, N. Manini, and F. Toigo, Phys. Rev. A 77, 043609 �2008�. �48� D. Ananikian and T. Bergeman, Phys. Rev. A 73, 013604 �2006�. �49� B. Xiong, J. Gong, Han Pu, W. Bao, and B. Li, Phys. Rev. A 79, 013626 �2009�. �50� P. Muruganandam and S. K. Adhikari, Comput. Phys. Com- mun. 180, 1888 �2009�. �51� P. Muruganandam and S. K. Adhikari, J. Phys. B 36, 2501 �2003�; S. K. Adhikari and P. Muruganandam, ibid. 35, 2831 �2002�. �52� A. Spuntarelli, P. Pieri, and G. C. Strinati, Phys. Rev. Lett. 99, 040401 �2007�. ADHIKARI, LU, AND PU PHYSICAL REVIEW A 80, 063607 �2009� 063607-10