# Component separation in harmonically trapped boson-fermion mixtures

## I Introduction

Since the recent experimental realization of Bose-Einstein condesation in dilute gases of rubidium , sodium , lithium , and hydrogen  a great deal of interest in Bose condensed systems has concentrated on the topic of multi-component condensates. This field was stimulated by the succesful demonstration of overlapping condensates in different spin states of rubidium in a magnetic trap  and of sodium in an optical trap , the (binary) mixtures being produced either by sympathetic cooling, which involves one species being cooled to below the transition temperature only through thermal contact with an already condensed Bose gas, or by radiative transitions out of a single component condensate. Since then a host of experiments has been conducted on systems with two condensates, exploring both the dynamics of component separation , and measuring the relative quantum phase of the two Bose-Einstein condensates . Most of the theoretical work concerning multi-component condensates  has been devoted to systems of two Bose condensates. However, other systems are of fundamental interest, one of these being a Bose condensate with fermionic impurities, a system reminiscent of superfluid <sup>3</sup>He-<sup>4</sup>He mixtures. In particular the possibilty of sympathetic cooling of fermionic isotopes has been predicted in both $`{}_{}{}^{6}\mathrm{Li}`$-$`{}_{}{}^{7}\mathrm{Li}`$ , $`{}_{}{}^{39}\mathrm{K}`$-$`{}_{}{}^{40}\mathrm{K}`$, and $`{}_{}{}^{41}\mathrm{K}`$-$`{}_{}{}^{40}\mathrm{K}`$ . Magneto-optical trapping of the fermionic potassium isotope $`{}_{}{}^{40}\mathrm{K}`$ has been reported .

The boson-fermion mixture was discussed in a previous paper  within the Thomas-Fermi approximation, which amounts to neglecting the kinetic energy of the bosons, and to apply a semi-classical filling of phase space of the fermions. For the bosons, this is a valid approximation in the limit of strong interactions or large particle numbers, see . In this paper we present a numerical analysis of the system, incorporating the correct operator form of the kinetic energy of the particles.

The paper is structured as follows. In Sec. II we study in detail the case of an isotropic external potential and we develop both the Thomas-Fermi approximation and the full quantum mechanical description of the fermions. The numerical procedure is briefly introduced. In Sec. III the case of the anisotropic harmonic oscillator trap is outlined within the Thomas-Fermi approximation for the fermions. In Sec. IV we present our quantitative results for the isotropic and anisotropic trapping potentials, demonstrating the accuracy of the predictions made in , and addressing the issue of symmetry breaking in elongated traps. Sec. V summarizes the main results.

Throughout, we assume that the bosons and fermions have the same mass, $`M`$, and that the atoms are all trapped in the same external harmonic oscillator potential. This choise is of course only a convenience; all our calculations are readily generalized to differing experimental parameters.

## II Isotropic traps

### A Gross-Pitaevskii equation for the bosons

In the mean field description the behavior of the single particle wavefunction $`\psi (\stackrel{}{r})`$, assumed to describe all $`N_B`$ bosons in the gas, is governed by the Gross-Pitaevskii equation . In the presence of fermions, this equation is modified by the addition of an interaction term proportional to the fermion density, $`n_F(\stackrel{}{r})`$

$$\left[\frac{\mathrm{}^2}{2M}^2+V_{ext}(\stackrel{}{r})+gN_B|\psi (\stackrel{}{r})|^2+hn_F(\stackrel{}{r})\right]\psi (\stackrel{}{r})=\mu \psi (\stackrel{}{r}),$$
(1)   
where $`V_{ext}(\stackrel{}{r})=\frac{1}{2}M(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)`$ is the external confining potential, and $`\mu `$ is the boson chemical potential (energy per particle). The value of $`\mu `$ is fixed by the normalisation condition, $`d^3rn_B(\stackrel{}{r})=N_B`$ on the boson density $`n_B(\stackrel{}{r})=N_B|\psi (\stackrel{}{r})|^2`$.

The low kinetic energies of the atoms permit the replacement of their short range interaction potential by a delta function potential of strength g or h. This is known as the pseudopotential method . There is no fermion-fermion interaction in this description, see below. In (1) g and h thus represent the boson-boson and the boson-fermion interaction strengths proportional to the respective $`s`$-wave scattering lengths .

In isotropic traps we have $`V_{ext}(\stackrel{}{r})=\frac{1}{2}M\omega ^2r^2`$, r being the distance from the trap center. By the substitution $`\chi =r\psi `$ in (1) we obtain the radial equation

$$\frac{\mathrm{}^2}{2M}\frac{d^2\chi }{dr^2}+\left[V_{ext}(r)+gN_B\left|\frac{\chi (r)}{r}\right|^2+hn_F(r)\right]\chi (r)=\mu \chi (r).$$
(2)   
In order to simplify the formalism, we rescale (2) in terms of harmonic oscillator units, that is

$`\stackrel{}{r}`$ $`=`$ $`a_0\stackrel{~}{\stackrel{}{r}},`$ (3)   
$`\mu `$ $`=`$ $`\mathrm{}\omega \stackrel{~}{\mu },`$ (4)   
$`\psi (\stackrel{}{r})`$ $`=`$ $`a_0^{3/2}\stackrel{~}{\psi }(\stackrel{~}{\stackrel{}{r}}),`$ (5)   
$`\chi (r)`$ $`=`$ $`a_0^2\stackrel{~}{\chi }(\stackrel{~}{r}),`$ (6)   
where $`a_0=\sqrt{\mathrm{}/M\omega }`$ is the width of the oscillator ground state. Defining

$$\stackrel{~}{g}=\frac{gM}{a_0\mathrm{}^2},\stackrel{~}{h}=\frac{hM}{a_0\mathrm{}^2},$$
(7)   
we arrive at the simplified equation for the radial function

$$\frac{1}{2}\frac{d^2\stackrel{~}{\chi }}{d\stackrel{~}{r}^2}+\left[\frac{1}{2}\stackrel{~}{r}^2+\stackrel{~}{g}N_B\left|\frac{\stackrel{~}{\chi }(\stackrel{~}{r})}{\stackrel{~}{r}}\right|^2+\stackrel{~}{h}\stackrel{~}{n_F}(\stackrel{~}{r})\right]\stackrel{~}{\chi }(\stackrel{~}{r})=\stackrel{~}{\mu }\stackrel{~}{\chi }(\stackrel{~}{r}).$$
(8)   
In the remaining parts of the paper we shall omit the tilde from this equation.

To solve the Gross-Pitaevskii equation for the bosons, we must find $`n_F(\stackrel{}{r})`$. To this end we invoke two methods: A semi-classical (Thomas-Fermi) approximation and a quantum mechanical treatment.

### B Thomas-Fermi approximation for the fermions

In the semi-classical (Thomas-Fermi) approximation the particles are assigned classical positions and momenta, but the effects of quantum statistics are taken into account. That is: The density in the occupied part of phase space is simply $`(2\pi )^3`$, and sums over states can be replaced by the corresponding integrals over $`\stackrel{}{r}`$ or $`\stackrel{}{k}`$. The fermions experience a potential $`V(\stackrel{}{r})=V_{ext}(\stackrel{}{r})+hn_B(\stackrel{}{r})`$ and for particle motion in such a potential it is posible to define a local Fermi vector $`\stackrel{}{k}_F(\stackrel{}{r})`$ by

$$E_F=\frac{\mathrm{}^2k_F(\stackrel{}{r})^2}{2M}+V(\stackrel{}{r}),$$
(9)   
so that the volume of the local Fermi sea in $`k`$ space is simply

$$\frac{4}{3}\pi k_F(\stackrel{}{r})^3=(2\pi )^3n_F(\stackrel{}{r}).$$
(10)   
In the low temperature limit, where $`p`$-wave (and higher multipole) scattering can be neglected, the supression of the $`s`$-wave scattering amplitude due to the antisymmetry of the many-body wavefunction implies that the spin polarized fermions constitute a noninteracting gas (for the case of an interacting Fermi gas, see ). Hence the density of the fermionic component is given by

$$n_F(\stackrel{}{r})=\left\{\frac{2M}{\mathrm{}^2}\left[E_FV_{ext}(\stackrel{}{r})hn_B(\stackrel{}{r})\right]\right\}^{3/2}/(6\pi ^2).$$
(11)   
As in the case of the bosons, where the chemical potential must be adjusted for the integral of the density over space to give the correct number of particles, the Fermi energy determines the proper normalisation; $`d^3rn_F(\stackrel{}{r})=N_F`$. For a thorough discussion of trapped fermions (also at $`T>0`$), and comments on the range of validity of the Thomas-Fermi approximation see .

### C Slater determinant description

The many-body wavefunction, $`\mathrm{\Psi }(\stackrel{}{r}_1\mathrm{}\stackrel{}{r}_{N_F})`$, may be represented by a Slater determinant

$$\mathrm{\Psi }(\stackrel{}{r}_1\mathrm{}\stackrel{}{r}_{N_F})=\frac{1}{\sqrt{N!}}𝒜\underset{i=1}{\overset{N_F}{}}\phi _i(\stackrel{}{r}_i),$$
(12)   
where $`𝒜`$ is the antisymmetrization operator. This Slater determinant solves a stationary Schrödinger equation

$$\widehat{H}\mathrm{\Psi }(\stackrel{}{r})=E\mathrm{\Psi }(\stackrel{}{r}),$$
(13)   
with a Hamiltonian that is the sum of $`N_F`$ independent single-particle operators

$`\widehat{H}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_F}{}}}\widehat{H}_i,`$ (14)   
$`\widehat{H}_i`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2M}}_{r_i}^2+{\displaystyle \frac{1}{2}}M\omega ^2r_i^2+hn_B(\stackrel{}{r}_i).`$ (15)   
The orbitals $`\phi _i(\stackrel{}{r}_i)`$ solve the eigenvalue equation

$$\widehat{H}_i\phi _i(\stackrel{}{r}_i)=E_i\phi _i(\stackrel{}{r}_i).$$
(16)   
We make the substitution $`\phi (\stackrel{}{r})=\frac{u_n\mathrm{}(r)}{r}Y_\mathrm{}m(\theta ,\varphi )`$, where $`Y_{lm}(\theta ,\varphi )`$ are the usual spherical harmonics, and we thus obtain a radial equation for the functions $`u_n\mathrm{}`$ in harmonic oscillator units

$$\frac{1}{2}\frac{d^2u_n\mathrm{}}{dr^2}+\left[\frac{1}{2}r^2+\frac{\mathrm{}(\mathrm{}+1)}{2r^2}+hn_B(r)\right]u_n\mathrm{}(r)=E_n\mathrm{}u_n\mathrm{}(r).$$
(17)   
It is important to keep in mind that the radial functions must satisfy the boundary condition $`u_n\mathrm{}(0)=0`$, to ensure a finite particle density at the center of the trap.

Equation (17) can be solved once for every $`\mathrm{}`$-value, thus producing the energy spectrum. The centrifugal term in the radial equation, implies that the fermions can be considered to move in an isotropic effective potential, $`V_{eff}(r)=V_{ext}(r)+hn_B(r)+\frac{\mathrm{}(\mathrm{}+1)}{2r^2}`$. The energy levels $`E_n\mathrm{}`$ are $`2\mathrm{}+1`$ times degenerate, and we have the fermion density given by

$$n_F(\stackrel{}{r})=\underset{\stackrel{occupied}{states}}{}\left|\frac{u_n\mathrm{}(r)}{r}Y_\mathrm{}m(\theta ,\varphi )\right|^2=\underset{n\mathrm{}}{}(2\mathrm{}+1)\frac{|u_n\mathrm{}(r)|^2}{4\pi r^2},$$
(18)   
since $`_{m=\mathrm{}}^{m=\mathrm{}}|Y_\mathrm{}m(\theta ,\varphi )|^2=(2\mathrm{}+1)/4\pi `$. Once found the eigenstates are sorted after energy and the energy levels are filled from below with $`N_F`$ particles. The Fermi energy is the energy of the highest occupied orbital. The maximum value of the angular momentum, may be estimated from the Thomas-Fermi expressions (10,11) for $`h=0`$ by maximizing $`rk_F(\stackrel{}{r})`$, the maximal length of $`\mathrm{}`$ at the point $`\stackrel{}{r}`$. This yields the simple result, $`\mathrm{}_{max}E_F`$, where the Fermi energy in the noninteracting limit is $`E_F=(6N_F)^{1/3}`$ in harmonic oscillator units. To test our numerical calculations for fermions not interacting with the bosons ($`h=0`$), we have compared our spatial density distributions with those of Schneider and Wallis  and found excellent agreement.

### D Numerical Procedure

We note that the solution of both Eqs. (11) and (17) require prior knowledge of the boson density $`n_B(\stackrel{}{r})`$. To obtain the density profiles of the two components, we insert iteratively the density of one component into the equation for the other until a desired convergence is reached.

To solve (8) for the boson density, we use the method of steepest descend, that is we propagate a trial function (which can be chosen initially almost arbitrarily) in imaginary time $`\tau `$, replacing $`\mu \chi (r)`$ by $`(/\tau )\chi (r,\tau )`$. In the long time limit the propagation “filters” the trial function to the condensate ground state Alternative methods for solving numerically the Gross-Pitaevskii equation are presented in .

The evaluation of the fermion density profile is done by two methods, as described above. In the case of the Thomas-Fermi approximation, $`n_F(\stackrel{}{r})`$ is found by direct insertion of $`n_B(\stackrel{}{r})`$ into (11), searching numerically for the energy $`E_F`$ giving the right number of particles. Within the Slater determinant method one obtains the density profile directly from (18), once the diagonalization of (17) has been done.

## III Anisotropic traps

In this section we treat the case of an anisotropic trapping potential with a cylindrical geometry ($`\omega _x=\omega _y=\omega _{}\omega _z`$) as this corresponds to current experimental setups We thus have 

$$V_{ext}=\frac{1}{2}M\omega _{}^2r^2+\frac{1}{2}M\omega _z^2z^2,$$
(19)   
where $`r=\sqrt{x^2+y^2}`$ and $`z`$ are the radial and axial coordinate respectively. We define the asymmetry parameter $`\lambda =\omega _z/\omega _{}`$.

As in the case of the isotropic potential we have for the bosons a non-linear Schrödinger equation corresponding to (1). By the substitution $`\chi =\sqrt{r}\psi `$ we obtain the equation

$`\mu \chi (r,z)=`$ $``$ $`{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{^2\chi }{r^2}}+{\displaystyle \frac{^2\chi }{z^2}}\right]{\displaystyle \frac{\chi (r,z)}{8r^2}}+{\displaystyle \frac{1}{2}}(r^2+\lambda ^2z^2)\chi (r,z)`$ (20)   
$`+`$ $`gN_B{\displaystyle \frac{\left|\chi (r)\right|^2}{r}}\chi (r,z)+hn_F(\stackrel{}{r})\chi (r,z),`$ (21)   
in harmonic oscillator units. Again the radial function has to vanish on the symmetry axis to remove potential problems of divergences near the origin. This boundary condition is implemented in our numerical procedure by imposing on the radial function a $`\sqrt{r}`$ dependence for small values of $`r`$, fitting to the value of $`\chi `$ at larger distances from the axis. As in the case of the isotropic trapping potential a Split-Step-Fourier technique is used to propagate the boson wavefunction to the ground state. An alternative method for solving the Gross-Pitaevskii equation in a cylindrical configuration, applying an alternating-direction implicit method to compute the derivatives is discussed in .

We shall limit ourselves to the Thomas-Fermi approximation for the fermions, both out of necessity and convenience. Already in the spherically symmetric, effectively 1-dimensional case, the full quantum mechanical analysis is very time consuming and as we shall demonstrate in the next section, the Thomas-Fermi approximation offers the same qualitative features as the exact description. The fermion density is thus evaluated using equation (11) with the external potential (19) and with the boson density obtained from (21).

## IV Results

The main conclusion of  is the prediction of a component separation under variation of the strength of the boson-boson and boson-fermion interaction. In the Thomas-Fermi approximation for both components, the density distributions solve the coupled equations

$`V_{ext}(\stackrel{}{r})+gn_B(\stackrel{}{r})+hn_F(\stackrel{}{r})`$ $`=`$ $`\mu `$ (22)   
$`{\displaystyle \frac{\mathrm{}^2}{2M}}(6\pi ^2n_F(\stackrel{}{r}))^{2/3}+V_{ext}(\stackrel{}{r})+hn_B(\stackrel{}{r})`$ $`=`$ $`E_F.`$ (23)   
In the case of $`N_F/N_B1`$ the fermions may be neglected in the equation for the bosons. For the fermions we then obtain the simple equation

$`{\displaystyle \frac{\mathrm{}^2}{2M}}(6\pi ^2n_F(\stackrel{}{r}))^{2/3}+(1{\displaystyle \frac{h}{g}})V_{ext}(\stackrel{}{r})+{\displaystyle \frac{h}{g}}\mu =E_F,`$ (24)   
where the terms proportional with $`h`$ are absent in regions with vanishing $`n_B(\stackrel{}{r})`$. We may distinguish between 3 different types of solutions; if $`h<g`$ the potential minimum of the fermions is located at the center of the trap, and if their number is small enough, they will constitute a ’core’ entirely enclosed within the Bose condensate. The two quantum gases are truly interpenetrating. If $`h=g`$ the fermions have a constant density throughout the Bose condensate, falling towards zero outside. If $`h>g`$ the effective potential for the fermions is that of an inverted harmonic oscillator having a minimum at the edge of the Bose condensate, where the fermions localize as a ’shell’ wrapped around the condensate.

### A Isotropic trap, quantum treatment

When we replace the Thomas-Fermi approximation by an exact description including the kinetic energy operator for the bosons and treating the fermions quantum mechanically, we expect to observe the same overall behaviour, but with minor corrections. The boson kinetic energy is expected to cause penetration into the fermionic component and a rounding off of the atomic distributions at the boundaries. Fig. 1 shows the spatial distribution of 1000 fermions in a condensate of $`10^6`$ bosons for different values of the boson-fermion interactions strength, $`h`$. The strength of the boson-boson interaction, $`g`$, is chosen to give maximal overlap between the two atomic clouds. In order to have clouds of comparable size, we equate the Thomas-Fermi expressions for the radius of the Bose condensate $`(15N_Bg/4\pi M\omega ^2)^{1/5}`$, and the radius of the zero temperature Fermi gas $`(48N_F)^{1/6}\sqrt{\mathrm{}/M\omega }`$. This gives the condition:

$$g/(\mathrm{}\omega a_{0}^{}{}_{}{}^{3})21.1N_{F}^{}{}_{}{}^{5/6}/N_B,$$
(25)   
which for the parameters of Fig. 1 requires $`g=0.015`$. The coupling $`g`$ differs for different atomic species and this value is in approximate agreement with the coupling stregth in the MIT Na setup , and we recall the possibility to achieve couplings of arbitrary strength by the recently demonstrated modification of the atomic scattering length by external fields . This allows a ’tuning’ of the scattering length through both positive and negative values. Finally we recall that we have insisted on equal masses and trapping potentials for the two components. If these constraints are relaxed, we may more easily vary the values of the scaled interaction strengths.

The oscillations in the fermion density distribution near the trap center reflect the matter wave modulation of the particles in the outermost shell. Their de Broglie wavelength can be estimated in the Thomas-Fermi approximation from (9): In the center of the trap the particles in the $`\mathrm{}=0`$ states experience a vanishing potential for $`h=0`$. As the Thomas-Fermi expression for the Fermi energy of $`N_F`$ fermions in a harmonic potential is $`E_F=(6N_F)^{1/3}\mathrm{}\omega `$ we find for the de Broglie wavelength

$$\lambda _{DB}\frac{2\pi }{k_F(0)}\frac{\sqrt{2}\pi a_0}{(6N_F)^{1/6}}1a_0,$$
(26)   
an estimate that is reproduced by the data, see inset.

We now turn to the case of equal numbers of bosons and fermions. The influence of the inter-species interaction grows as the number of fermions is increased with dramatical effects on the atomic distibutions, as we shall demonstrate. We study the case of $`10^6`$ fermions, and the same number of bosons with an interaction strength of $`g=2.11\mathrm{}\omega a_{0}^{}{}_{}{}^{3}`$. We again expect that for certain critical parameters, the components find it energetically favorable to separate into two distinct phases, but this time bosons are expelled from the trap center, minimizing their internal interaction energy by spreading in a ’shell’ around a fermionic bubble. Figs. 2 and 3 present our results. The essential features are again the spatial separation of the two components, this time manifesting itself by the exclusion of the bosonic component from the trap center and the existence of a constant fermion density through the boson distribution for $`h=g`$. For a different choise of parameters, for example by letting the fermions be trapped by a weaker potential, we are also capable of producing a multi-layered structure with fermions residing on both sides of the bosons.

We notice that as the bosons are expelled from the center of the trap, forming a ’mantle’ around the fermions, the fermionic component is compressed, having a higher peak density and covering a smaller portion of the trapping volume. A similar behavior has been noted for bi-condensate systems . One of the essential features predicted in the Thomas-Fermi approximation is the existence of a ’plateau’ of constant fermion density through the boson distribution for $`h=g`$. As illustrated by Fig. 3, which is just a magnification of the central parts of Fig. 2e, this phenomenon also appear in our quantum mechanical treatment, although with the parameters chosen it does not involve quite as many particles as obtained from the semi-classical calculations in .

It is interesting to compare the above mentioned results with those obtained by treating the fermions in the Thomas-Fermi approximation. This is done in Fig. 4 for $`N_B=N_F=10^6`$, $`h=g=2.11\mathrm{}\omega a_0^3`$, and we note that the semi-classical description gives a qualitatively correct description, in that it reliably predicts the phase separation. Thus it is reasonable to use this approximative treatment of the fermions in the anisotropic case, where the exact description is too cumbersome.

### B Anisotropic trap

We now turn our attention to the anisotropic potentials, where we will use only the semi-classical Thomas-Fermi approximation for the fermion density. We aim to reveal similar variations in the ground state density profiles as for the isotropic trap, but going to higher dimensions we now have the opportunity to investigate the phenomenon of spatial symmetry breaking. Intuitively, we assume that for critical parameters it may be preferable for the two components to break mirror symmetry ($`zz`$), thereby minimizing their mutual interaction, especially in elongated traps. Such a behavior has been predicted by Öhberg and Stenholm for bi-condensates in two dimensions .

It remains to be demonstrated though, that the features described in the case of the isotropic trap are still essential, when we consider the anisotropic scenario relevant in comparison with currently experimentally feasible setups. We present in Figs. 5 and 6 the analog of Fig. 1 with the same choise of parameters and $`\lambda =1/\sqrt{8}`$, *i.e.* the inverse of the value for the current traps which have the strongest confinement along the $`z`$-axis. We notice the appearance of the same qualitative features as in the isotropic trap, that is component separation for $`h>g`$ and a plateau of constant fermion density for $`h=g`$.

The $`10^6`$ bosons are in the condensate which is unaffected in form and location by the presence of the relatively few fermions. Not shown in Figs. 5 and 6 is the distribution of fermions for $`h`$ smaller than $`g`$. In this case the fermionic component overlaps the boson cloud at the center of the trap.

To address the issue of symmetry breaking we adopt the same procedure as Öhberg and Stenholm . This offers only suggestive evidence that symmetry breaking may occur. To investigate this behavior correctly one must use an altogether different approach, minimizing the energy functional to find the ground-state density profile . The point is that the solutions of the Gross-Pitaevskii equation are stationary points of the energy functional, not necessarily corresponding to minima. They may therefore be unstable in certain parameter regions. It is possible though to single out the more stable of two configurations by comparing their total energy as this is minimum in equilibrium.

The total energy functional of the two-component system is a sum of four terms

$$E=T_B+T_F+V_{ext}+V_{int}.$$
(27)   
The first term is the boson kinetic energy

$$T_B=d^3r\frac{\mathrm{}^2}{2M}|\psi (\stackrel{}{r})|^2.$$
(28)   
As a fermion with wave number $`\stackrel{}{k}(\stackrel{}{r})`$ has a kinetic energy of $`\mathrm{}^2k^2/2M`$, the total fermionic contribution to the kinetic energy is found by integrating this local term over all of phase-space, weighted by the phase-space density, $`1/(2\pi )^3`$,

$`T_F`$ $`=`$ $`{\displaystyle \frac{d^3r}{(2\pi )^3}_0^{k_F(\stackrel{}{r})}d^3k\frac{\mathrm{}^2k^2}{2M}}`$ (29)   
$`=`$ $`{\displaystyle \frac{d^3r}{2\pi ^2}\frac{\mathrm{}^2}{10M}\left[6\pi ^2n_F(\stackrel{}{r})\right]^{5/3}}.`$ (30)   
Calculating the potential energy terms is easy, as they involve only integrals over the atomic densities

$`V_{ext}`$ $`=`$ $`{\displaystyle d^3r\frac{1}{2}M\omega ^2r^2\left[n_B(\stackrel{}{r})+n_F(\stackrel{}{r})\right]}`$ (31)   
$`V_{int}`$ $`=`$ $`{\displaystyle d^3rn_B(\stackrel{}{r})\left[gn_B(\stackrel{}{r})+hn_F(\stackrel{}{r})\right]}.`$ (32)   
We have chosen the number of atoms to be $`N_B=N_F=10^3`$, while the asymmetry parameter is still set to $`\lambda =1/\sqrt{8}`$. The interaction parameters are $`g=6.67\mathrm{}\omega a_0^3`$ and $`h=5g`$.

Starting the iteration with two well separated clouds displaced along the cylinder axis, *i.e.* along the direction of the weaker trapping potential, the calculation converges to a situation where the fermions localize on both sides of a central concentration of the Bose condensate: a ’boson-burger’, see Fig. 7. Initiating the calculation with two overlapping clouds in the center of the trap results in just the reversed situation: a ’fermion-burger’, consisting of a central fermionic part surrounded on two sides by bosons, but this configuration has a larger energy. The ’boson-burger’ seems to be the stable solution.

In Fig. 8 we show the spatial distribution of 5000 fermions and 1000 bosons. The particles feel the same trapping potential as in Fig. 7, and the interaction strengths are kept unchanged. This configuration is the result when the starting point of the calculation is two separated clouds. When we start by placing both species at the center of the trap we achieve again the ’fermion-burger’, but at a higher energy. Thus we conclude that in this region of parameter space the system is unstable against breaking of the reflection symmetry.

We note that our approach provides two degenerate symmetry broken states, the one in Fig. 8 and its mirror image in the $`xy`$-plane. Going beyond our theoretical treatment (Hartree), we may construct superpositions of these two macroscopically states which do not break the spatial symmetry. One of these states will have a lower energy, but such a ’Schrödinger-cat’ state is exceedingly complicated to prepare, *c.f.* the discussion in . Thus the symmetry broken solution is most likely to be observed in an experiment.

## V Conclusion

In this paper we have investigated the zero temperature ground state of a mixture of boson and fermion gases in both isotropic and anisotropic trapping potentials. We have addressed the issue of component separation using nummerical techniques to solve coupled equations for the spatial density of the two species. Our calculations have confirmed and expanded upon the results of a previous paper , which treated the problem only within the Thomas-Fermi approximation for both components and which analyzed only the case of an isotropic trap. We have confirmed the existence of three distinct states of the system under variation of the ratio of the interaction strengths $`h/g`$: For small values of this parameter the gases are interpenetrable, overlapping throughout the occupied volume of the trap, as their mutual repulsion is not strong enough to cause separation. When the coupling strength $`h`$ exceeds the strength of the boson-boson interaction one of the species is expelled from the center of the trap. The spatial configuration in this case depends on the symmetry of the trapping potential. In an isotropic trap the separated phase is rotationally symmetric, the excluded component constitutes a spherical shell wrapped around a centrally compressed bulk. The anisotropic trap however has a parameter region where a breaking of symmetry ($`zz`$) may occur, and we have demonstrated such forms. In the limiting case $`h=g`$ there exists the possibility for the fermions to have a constant spatial density where the bosons are localized.

An aspect of this work is the availability of an almost isolated degenerate Fermi gas through the complete separation of the two species. The trapped, degenerate Fermi gas is interesting in view of the possibility of a BCS transition when two spin states are trapped simultaneously  and because of the analogies between this system and both atomic nuclei and the interior of neutron stars.

The details of sympathetically cooling the Fermi gas to the degeneracy level through thermal contact with the Bose condensate are of course of great importance in further research . In general the investigation of the cooling ability of the condensate should not be restricted to fermionic impurities. In view of the recent trapping of simple molecules in both optical  and magnetic  potentials, also more complex solutes with several internal degrees of freedom pose an interesting challenge for future research.

Another direction worth noticing is the prospect of trapping a boson-fermion mixture in the periodic potential of an optical lattice , both in its own right and as a study of solid state phenomena. With quantum gases well beyond the degeneracy level a complete filling of the potential wells may well be expected .

Finally it should be mentioned that in this work we have concentrated solely on systems with a positive coupling strength $`h`$. Allowing the interaction between the species to become attractive is known to induce a dramatic change in the macroscopic behavior of the system as it becomes unstable against collapse for large negative values of $`h`$ . We are currently setting up calculations to investigate this phenomenon in detail using the numerical procedure developed in this work.