K-Rb Fermi-Bose mixtures: vortical states and sag

We study a confined mixture of bosons and fermions in the quantal degeneracy regime with attractive boson-fermion interaction. We discuss the effect that the presence of vortical states and the displacement of the trapping potentials may have on mixtures near collapse, and investigate the phase stability diagram of the K-Rb mixture in the mean field approximation supposing in one case that the trapping potentials felt by bosons and fermions are shifted from each other, as it happens in the presence of a gravitational sag, and in another case, assuming that the Bose condensate sustains a vortex state. In both cases, we have obtained an analytical expression for the fermion effective potential when the Bose condensate is in the Thomas-Fermi regime, that can be used to determine the maxima of the fermionic density. We have numerically checked that the values one obtains for the location of these maxima using the analytical formulas remain valid up to the critical boson and fermion numbers, above which the mixture collapses.


I. INTRODUCTION
Recent experiments on degenerate Fermi-Bose mixtures have opened the possibility of studying in a direct way the effect of quantum statistics in Bose-Einstein condensates (BECs). In practice, the direct evaporative cooling techniques used to obtain BEC are not applicable in a Fermi gas, as the Pauli exclusion principle forbids s-wave collisions between fermions. This difficulty has been overcome using a gas of bosons as a coolant, which has given further relevance to the study of these mixtures.
One of the systems studied in recent experiments is the 6 Li-7 Li mixture [1,2]. It is characterized by having a positive boson-fermion scattering length an order of magnitude larger than the boson-boson one [2]. Due to this repulsive boson-fermion interaction, the system does not exhibit a large overlap between the two species [3]. Moreover, this mixture may undergo a two-component separation [3,4]. Both facts conspire against having the fermions well inside the boson cloud, which is desirable to obtain an efficient sympathetic cooling. In contrast, this desirable overlap can be achieved in 40 K-87 Rb mixtures due to the large attractive boson-fermion interaction. This mixture has been recently obtained [5,6,7], and it has been shown that if the particle numbers are above some critical values the system collapses [7].
From a theoretical point of view, a fairly amount of work has been devoted to the study of boson-fermion mixtures [3,8,9]. In particular, a systematic study of the structure of binary mixtures has been performed in Ref. [3], where all possible sign combinations of scattering lengths between boson-boson and boson-fermion s-wave interactions have been discussed. The purpose of our work is, firstly, to present an analysis of the location of the minima of the effective potential felt by fermions submitted to a large number of condensate bosons when the boson-boson interaction is repulsive, and secondly, to extend this study to systems with a large number of fermions with an attractive boson-fermion interaction.
An essential characteristic of superfluid systems is the occurrence of quantized vortices [10]. Quantized vortices in BECs were first produced experimentally by Matthews et al. [11], and have been studied since in detail (see Refs. [12,13,14] and Refs. therein). It may thus be interesting to see what is the physical appearance of such vortices, arising in the bosonic component, in the presence of a fermionic cloud which is in the normal (non superfluid) but quantum degenerate phase. This is another aim of our work, which bears some similarities with the description of quantized vortices in 3 He-4 He nanodroplets recently addressed [15]. This work is organized as follows. In Sec. II we describe the mean-field model we have used. In Sec. III we present an analysis of the effective potential felt by the fermions when they are inside a large condensate that exhibits a boson-boson repulsive interaction, and the number of fermions (N F ) is much smaller than the number of bosons (N B ). We shall consider the case in which the boson and fermion external confining potentials are displaced from each other, which may be caused by gravity for instance, and also when bosons are in a vortical state. In Sec. IV we present the results obtained by solving the mean-field coupled equations for different N B and N F values up to the critical values where collapse of the mixture occurs. A summary is given in Sec. V. Finally, in the Appendix we derive some expressions using a scaling transformation which, on the one hand, constitute a generalization of the virial theorem, and on the other hand are especially useful for testing the accuracy of the numerical procedure.

II. THE MEAN-FIELD MODEL
We consider a mixture of a Bose condensate (B) and a degenerate Fermi gas (F) at zero temperature confined in an axially symmetric harmonic trap. Assuming that the minimum of the trapping potential felt by each species may be displaced in the z direction in a value d i (i=B,F), the confining potentials in cylindrical coordinates are and where ω rB , ω zB and ω rF , ω zF are the trapping radial (axial) angular frequencies for bosons and fermions, respectively. M B and M F are the corresponding masses.
Since the number of fermions we consider in the numerical calculations is fairly large, the fermionic kinetic energy density can be written in the Thomas-Fermi-Weizsäcker (TFW) approximation as a function of the local fermion density n F and its gradients. For fully polarized spin 1/2 fermions, it reads: and the fermionic kinetic energy is The value of the β coefficient in the Weizsäcker term is fixed to 1/18. This term contributes little to the kinetic energy, and it is usually neglected [3]. However, its inclusion [16] has some advantages, see below. We refer the interested reader to Ref. [17] for a discussion on the accuracy of the TFW approximation (see also Refs. therein).
Neglecting all p-wave interactions, the energy density functional for the boson-fermion mixture at zero temperature has the form where n B = |Ψ| 2 is the local boson atomic density. The boson-boson and boson-fermion coupling constants g BB and g BF are written in terms of the s-wave scattering lengths a B and a BF as g BB = 4π a Bh 2 /M B and g BF = 4π a BFh 2 /M BF , respectively. We have defined When bosons sustain a quantized vortex line along the z-axis, the condensate wave function can be written as Ψ = |Ψ|e i m φ , where m = 1, 2, 3 . . . is the circulation number, yielding for the kinetic energy densitȳ We have considered only singly quantized vortices, that is m = 1. If a vortex is present in the condensate, bosons flow around the vortex core with quantized circulation, which yields the centrifugal term in the kinetic energy. We assume that the Fermi component is not superfluid, and consider that it is in a stationary state. This situation could be achieved experimentally waiting enough time after the generation of the vortex in the condensate, to let the drag force between bosons and fermions to dissipate.
Variation of E with respect to Ψ and n F keeping N B and N F fixed yields the following coupled Euler-Lagrange (EL) equations The inclusion of the Weizsäcker term in the energy density has the major advantage that it yields a EL equation for n F that is well behaved everywhere, avoiding the classical turning point problem when this term is neglected (Thomas-Fermi approximation). Moreover, solving Eq. (8) is not more complicated than solving the GP equation. This can be readily seen writing the later in terms of n B : which is formally equivalent to Eq. (8). We have discretized these equations using 9point formulas and have solved them on a 2D (r, z) mesh using a sufficiently large box.
The results have been tested for different sizes of the spatial steps (we have mostly used ∆r = ∆z ∼ 0.1µm). We have employed the imaginary time method to find the solution of these coupled equations written as imaginary time diffusion equations [18]. After every imaginary time step, the densities are normalized to the corresponding particle numbers. To start the iteration procedure we have used positive random numbers to build both normalized densities. This avoids to introduce any bias in the final results. We have also checked that the solutions fulfill the generalized virial theorem deduced in the Appendix.

III. EFFECTIVE FERMION POTENTIAL
Before we present the numerical results obtained by solving Eqs. (8) and (9), it is useful to analyze the features of the effective potential felt by a small number of fermions in the presence of a boson condensate, paying special attention to the location of its minima. We are interested in studying the cases in which either the condensate hosts a vortex line or there exists a displacement between the minima of the boson and fermion external potentials. The later situation is routinelly met in the experiments -a gravitational sag in the z direction-, in which case the displacement is g(1/ω 2 zF − 1/ω 2 zB ), being g the acceleration of gravity. An interesting issue is that the analytic expressions we obtain in this Sec. for the location of the minimum of the effective potentials remain valid, as we will show in the next Sec., when the number of fermions and bosons are similar. These locations coincide with the positions of the maxima of the fermion density and are relevant because the collapse of the mixture originates precisely around them. This will be discussed in the next Sec.
The effective fermion potential has two main terms, the external potential and the meanfield term arising from the interaction with the boson condensate, which is proportional to n B . We consider a large condensate in the Thomas-Fermi (TF) regime with positive scattering length and we assume that its density profile is not affected by the fermion presence. To obtain n B we may thus use Eq. (7) with n F = 0 and neglect the first kinetic energy term.

A. Effect of a shift in the minimum of the external potentials
For simplicity, we restrict ourselves to the vortex-free case m = 0. Without loss of generality, the origin of coordinates can be fixed at the minimum of the boson trapping potential, and we will assume that the displacement is only in the z direction. Thus, the effective trapping potential experienced by the fermions is where d z = d F − d B is the z displacement between the trapping potential centers. Assuming that the density profile of the condensate is not affected by the interaction with fermions, the number of bosons is large, and that the interaction between bosons is repulsive (g BB > 0), one obtains from Eq. (7) setting m = 0 the boson density profile that corresponds to the uncoupled TF density with Θ(x) = 1 if x > 0 and zero otherwise. The TF condensate boundary is given by the ellipsoid arg(Θ) = 0. For an axially symmetric trap it yields the well-known TF radius of the condensate R i = 2µ B /(M B ω 2 iB ) with i = r, z in the radial and axial direction, respectively. Replacing n B into Eq. (10) and defining the dimensionless parameter with i = r, z, it follows that the effective fermion potential inside the Bose condensate -where whereas outside the condensate This effective potential can be viewed as having a renormalized frequency inside the boson condensate, a feature already discussed in Ref. [9] for a mixture with concentric external potentials. This model has been called the double-parabola model by Vichi et al. [9], and it has also been used by Capuzzi and Hernández [19]. However, in these works the possibility of a displacement between the minima of the potentials has not been considered.
The extremum of V ef f d inside the Bose condensate is attained at the point Note that if the shift between the external potentials is d = (d x , d y , d z ), the effective potential has its extremum at d ′ = (γ −1 . Depending on the signs of γ i , d ′ may correspond to a maximum, a minimum or a saddle point.

B. Vortices
We consider now the effect of a quantized boson vortex line in the fermion distribution, without gravitational sag. In this case, the effective fermion potential is where n v B is the boson density hosting a vortex line along the z axis. This density is zero at r = 0, and reaches its maximum value on a circle of radius r 0 around that axis. In the TF approximation, n v B can be derived using Eq. (7) keeping the centrifugal term proportional to m, and may be approximately written as [12] n with r 0 = h m/(M B ω rB ). Using n v B we calculate the mimimum of V ef f v inside the condensate, which is located at z ′ = 0 and Thus, the minimum is attained at a circle of radius r ′ . We will see that for the 40 K-87 Rb mixture r ′ is very close to the radius at which the boson density has its maximum.

IV. NUMERICAL RESULTS
The system under consideration is a confined 40 K-87 Rb mixture. We have assumed spher- We display in Fig. 1 several boson and fermion density profiles as a function of z, considering an arbitrary displacement d z = 10 µm. They all correspond to N B = 10 5 , but for three different fermion numbers: N F = 0 , N F = 10 3 , and N F = 2.5 ×10 4 . The attractive boson-fermion interaction produces an enhancement of the density of both species in the overlap volume. However, this effect is reduced with respect to the concentric case due to z displacement [7].
Using the TF approximation, the radius of the Bose condensate can be estimated as R B = (15N B a B /a HO ) 1/5 a HO , with a HO = h/M B ω B being the oscillator length of the boson trap; this yields R B ∼ 6µm. Following the analysis performed in Sec. III.A, and using γ −1 z = 0.136, we find that d ′ z = 1.36µm < R B ; thus, within that model there are two minima in the effective fermion potential. One is at d = 10 µmẑ, i.e., beyond the TF radius, and the other is inside the boson cloud, at d ′ = 1.36 µmẑ. It is interesting to note from Fig.   1 that in all cases, the maxima of the fermion density are precisely at z = d ′ z and z = d z , which we have displayed in the inset as vertical lines.
We thus see that for N F /N B equal to 1% (solid line) or even 25% (dashed line), the maxima of the fermion density appear at the values obtained using the model of Sec. III.A. This can be understood by solving Eqs. (7) and (8) 7), we obtain the following expression for n B : Replacing it into Eq. (8) and deriving the resulting expression with respect to z we get The extremum of the fermion density is found by setting ∂n F /∂z = 0, with n F = 0. After some algebra we find that the following equation has to be fulfilled: whose solution coincides precisely with the value z = d ′ z we have found in the previous Sec. Whereas the maximum of the fermion density is located at z = d ′ z irrespective of the value of the dilution (N F /N B ), it can be seen from the numerical results that the maximum of the boson density moves from z = 0 towards z = d ′ z as N F increases. In fact, Fig.  1 shows that even for a rather small N F value (N F = 10 3 , N B = 10 5 ), the shape of the condensate differs from the parabolic-type profile yielded by the TF approximation for a fermion-free condensate, Eq. (11). It can be also appreciated in the inset that for a dilution of N F /N B = 25%, a fairly large number of fermions remains without mixing (long tail in the fermion density profile outside the Bose condensate).
We plot in Fig. 2  Following a procedure analogous to that used before, we can justify the position of the maximum of the fermion density in the presence of a vortex in a TF condensate. From Eq.
(7) we obtain Replacing it into Eq. (8) and deriving the resulting expression with respect to r we get Once again, the extremum in the fermionic density is found by setting ∂n F /∂r = 0 with n F = 0, which means that the following equation has to be fulfilled: whose solution again coincides with r = r ′ as found in the previous Sec. [Eq. (18)].
Due to the attractive boson-fermion interaction, stable trapped 40 K-87 Rb mixtures may only have a limited number of fermions and bosons. If the atom numbers increase above some critical values N c B and N c F , an instability occurs [3]. It has been shown that the meanfield framework is able to reproduce the critical numbers for collapse [7]. We have calculated the stability diagram of the 40 K-87 Rb mixture by solving the coupled mean-field equations (8) and (9) for different values of N B and N F . In our study, the instability signature is the failure of the numerical iterative process. In particular, the instability onset appears as an indefinite growth of the maximum of the densities, which triggers the collapse.
We display in Fig. 3  From Fig. 3 one can conclude that, for the 40 K-87 Rb mixture, the presence of a vortex in the condensate, or any other mechanism that increases the distance between the maxima of the -still-overlapping densities, as for example a sag displacement between the two clouds, allows to have a stable mixture for larger particle numbers. The reason is that these mecha- In this Appendix we use a scaling transformation [14,20,21] to derive an expression that generalizes the virial theorem and is useful to test the accuracy of the numerical procedure.
Performing a scaling transformation of the vector position r → λ r, it is easy to check that to keep the normatization of the order parameter Ψ( r ) and of the boson and fermion atomic densities, they have to transform as: Splitting the energy of the system into a kinetic T , a harmonic oscillator U H , and an interaction part U g , the above transformations induce the following ones: where the last two expressions hold for fermions and bosons as well. The expression for U H supposes a spherically symmetric harmonic trap, the generatization to deformed harmonic traps is obvious.
It is easy to check that in the case of fermions, the kinetic energy in the TFW approximation also scales as T → T λ = λ 2 T , and that in the case of boson-fermion mixtures, the interaction term U g BF scales as U g BF → U g BF λ = λ 3 U g BF . Hence, the total energy of the mixture scales as If the scaling is carried out from the equilibrium configuration corresponding to the value λ = 1, one has One sees that at equilibrium, for non-interacting fermions and bosons one recovers the virial theorem, namely When a vortex is present, its superfluid kinetic energy scales as This is quite a natural result in view of the first Eq. (26) and, as a consequence, T V may be incorporated in the definition of T B which then represents the total kinetic energy of the bosonic component of the mixture.
If a gravitational sag is considered that displaces the atomic clouds in the z direction, its effect can be taken into account changing the confining potential in that direction into Apart from a trivial constant term proportional to the number of fermions and bosons in the trap that is invariant under the scaling transformation, there appear new contributions to the total energy of the kind U s = −Mz 0 ω 2 d r z n( r ) which scales as U s → U s λ = 1 λ U s . Eq. (28) then becomes The above expression is, in a way, a generalization of the virial theorem. We have used it to routinely check the accuracy in the numerical solution of the GP and TFW coupled equations. Values ∼ O(10) are found for that expression when U H is ∼ O(10 6 ). This makes us confident on the numerical method we have used to find the mean-field structure of the atomic mixture.