text:
# Effective s- and p-Wave Contact Interactions in Trapped Degenerate Fermi Gases

## I Introduction

The achievement of Bose-Einstein condensation in trapped dilute gases of bosonic atoms Anderson et al. (1995) triggered a wide spread interest in the field of ultracold atomic gases. Meanwhile these systems appear as a unique lab for the study of all kinds of fundamental quantum phenomena.

After a series of excellent experiments on bosonic systems the question arises, whether a dilute gas of fermionic atoms can also be cooled to temperatures where quantum effects dominate. In 1999 the group of Deborah Jin managed to cool a sample of typically $`10^6`$ fermionic 40K atoms to temperatures significantly below the Fermi energy of the system DeMarco and Jin (1999); in recent experiments they achieved a temperature corresponding to one fifth of the Fermi energy DeMarco et al. (2001). In this temperature regime the system can be described as a degenerate Fermi gas, where the majority of the atoms successively fills the lowest available one-body states according to the Pauli principle.

One of the goals of the investigations on dilute and ultracold Fermi gases is the observation of Cooper pairing and the transition to a superfluid state. The transition temperature depends on the density and on the strength of the attractive interaction that is necessary for the formation of Cooper pairs Houbiers et al. (1997); Stoof and Houbiers (1999). In order to increase the transition temperature one may increase the density or the interaction strength, where the latter seems to be more promising. An atomic species favored for the experimental observation of a BCS transition is 6Li due to its large s-wave scattering length of $`a_02160a_\text{B}`$ Abraham et al. (1997). But also 40K, which shows a rather small natural scattering length, is a possible candidate for Cooper pair formation DeMarco et al. (2001), because a simultaneous s- and p-wave Feshbach resonance is predicted Bohn (2000), which allows tuning of the s- and p-wave scattering lengths over a very wide range.

A serious constraint on the way towards a superfluid Fermi gas is the mechanical stability of the progenitor, i.e., the normal degenerate Fermi gas. As was experimentally demonstrated for the bosonic 85Rb system Cornish et al. (2000), an attractive interaction between the atoms leads to a mean-field instability of the trapped gas if the density exceeds a critical value. A similar collapse is expected in fermionic systems with attractive interactions. In contrast to the bosonic systems p-wave interactions contribute in a Fermi gas and may have strong influence on the stability of the system Roth and Feldmeier (2000).

In the following we address the question of the stability of degenerate one- and two-component Fermi gases in presence of s- and p-wave interactions within a simple and transparent model. In Section II we derive an Effective Contact Interaction (ECI) for all partial waves that is suited for a mean-field treatment of the many-body problem. Using the s- and p-wave part of this interaction we construct in Section III the energy functional of a trapped multi-component Fermi gas in Thomas-Fermi approximation. From that the ground state density profile can be determined by functional variation. In Section IV we discuss the structure and stability of single-component Fermi gases, where only the p-wave interaction contributes according to the Pauli principle. Section V deals with two-component Fermi gases where s- and p-wave interactions are present and lead to a subtle dependence of the stability on the two scattering lengths.

## II Effective Contact Interaction

### II.1 Concept

The approximate solution of the many-body problem in a restricted low-momentum sub-space of the full Hilbert space faces a fundamental problem: Many realistic two-body interactions — like van-der-Waals-type interactions between atoms or the interactions between two nucleons in an atomic nucleus — exhibit a strong short-range repulsion. This repulsion generates particular short-range correlations in the many-body state, which cannot be described within a low-momentum model-space Feldmeier et al. (1998). Therefore short-range correlations inhibit the use of a realistic two-body interactions in the framework of a naive mean-field model. One way to overcome this problem is to replace the original interaction by a suitable effective interaction. How this effective interaction has to be constructed from the original potential depends on the properties of the actual physical system under investigation.

Cold and very dilute quantum gases allow an effective interaction of simple structure. The typical wave length of the relative motion of two particles is always large compared to the range of the interaction. Therefore the particles experience only an “averaged” two-body potential and do not probe the detailed radial dependence. Moreover the gases are in a metastable not self-bound state that is kept together by the external trapping potential. Thus the bound states of the two-body potential are not populated and have only indirect influence. We make use of these facts and replace the original potential by a contact interaction. The strength of the contact terms is related to the properties of the original potential.

We establish this relation via the energy spectrum of the two-body system. Consider a two-body problem with an auxiliary boundary condition at some large radius such that the spectrum is discrete even for positive energies. For noninteracting particles this spectrum shows a sequence of levels at positive energies $`E_{nl}`$ which are labeled by an angular momentum quantum number $`l`$ and a radial quantum number $`n`$. A schematic sketch of this spectrum for some fixed $`l`$ is shown in Figure 1. If we switch on an interaction between the particles two things happen: A number $`n_l^\text{b}`$ of bound two-body states with relative angular momentum $`l`$ may appear at negative energies and the levels at positive energies $`\overline{E}_{nl}`$ are shifted compared to the noninteracting spectrum. For our purpose this energy-shift $`\mathrm{\Delta }E_{nl}=\overline{E}_{nl}E_{nl}`$ is defined between the positive energy states only. The lowest level of positive energy is labeled with the quantum number $`n=0`$.

The energy-shift is the relevant “averaged” property of the two-body potential that the effective interaction should reproduce. We require that the expectation value of the effective interaction calculated with eigenstates $`|nlm>`$ of the noninteracting two-body problem equal the energy shift $`\mathrm{\Delta }E_{nl}`$ induced by the original potential

$$<nlm\left|𝒗^{\text{eff}}\right|nlm>\stackrel{!}{=}\overline{E}_{nl}E_{nl}=\mathrm{\Delta }E_{nl}.$$
(1)
Due to the angular momentum dependence of the energy-shifts we have to construct the effective contact interaction as a sum of contact interaction for each partial wave. The strength of each contact term is fixed by this condition.

### II.2 Construction of the ECI

The construction of the Effective Contact Interaction (ECI) is organized in two steps: Firstly, the two-body energy-shift induced by the original potential is evaluated. Secondly, the operator of the Effective Contact Interaction is formulated for each partial wave and the interaction strengths are fixed by condition (1).

#### II.2.1 Energy Shift

Assume a system of two particles with reduced mass $`\mu =m_1m_2/(m_1+m_2)`$. Their wave function can be decomposed into a center of mass and a relative wave function, where the latter separates into a radial and a angular component because of rotational symmetry

$$<\stackrel{}{r}|nlm>=R_{nl}(r)Y_{lm}(\mathrm{\Omega }).$$
(2)
The radial wave functions $`R_{nl}(r)`$ of the noninteracting two-body system are solutions of the Schrödinger equation

$$\left[\frac{1}{r}\frac{^2}{r^2}r+\frac{l(l+1)}{r^2}2\mu E_{nl}\right]R_{nl}(r)=0$$
(3)
(we use units with $`\mathrm{}=1`$). We require the radial wave function to vanish at some arbitrary but large radius $`\mathrm{\Lambda }`$

$$R_{nl}(\mathrm{\Lambda })\stackrel{!}{=}0.$$
(4)
This auxiliary boundary condition leads to a discrete energy spectrum, which is needed to enumerate the energy levels and to evaluate the energy shift.

In the noninteracting case the solution of the radial Schrödinger equation (3) is given by a spherical Bessel function $`j_l(x)`$

$$R_{nl}(r)=A_{nl}j_l(q_{nl}r),$$
(5)
where $`q_{nl}`$ denotes the relative momentum of the two particles and $`A_{nl}`$ a normalization constant. The discrete momenta $`q_{nl}`$ are determined by the boundary condition (4) and thus are related the to the zeros of the Bessel function. Since the radius $`\mathrm{\Lambda }`$ can be chosen arbitrary large it is sufficient to use the asymptotic expansions of the spherical Bessel and Neumann functions

$$\begin{array}{cc}\hfill j_l(x)& \stackrel{xl}{=}\frac{1}{x}\mathrm{sin}\left(x\pi l/2\right),\hfill \\ \hfill n_l(x)& \stackrel{xl}{=}\frac{1}{x}\mathrm{cos}\left(x\pi l/2\right).\hfill \end{array}$$
(6)
Evaluating the boundary condition (4) with the asymptotic form of the Bessel function we obtain for the possible relative momenta

$$q_{nl}\mathrm{\Lambda }=\pi (n+l/2).$$
(7)
Accordingly the two-body energy spectrum in the noninteracting case is given by

$$E_{nl}=\frac{1}{2\mu }q_{nl}^2=\frac{\pi ^2}{2\mu \mathrm{\Lambda }^2}(n+l/2)^2.$$
(8)
The normalization constant can be determined explicitly

$$A_{nl}^2=_0^\mathrm{\Lambda }drr^2j_l^2(q_{nl}r)=\frac{\mathrm{\Lambda }^3}{2}j_{l+1}^2(q_{nl}\mathrm{\Lambda }).$$
(9)
Inserting the asymptotic form of the Bessel function we finally get

$$A_{nl}^2\stackrel{q_{nl}\mathrm{\Lambda }l}{=}\frac{\mathrm{\Lambda }}{2q_{nl}^2}.$$
(10)
In presence of a two-body potential $`v(r)`$ of finite range $`\lambda `$ the solution $`\overline{R}_{nl}(r)`$ of Schrödinger equation

$$\left[\frac{1}{r}\frac{^2}{r^2}r+\frac{l(l+1)}{r^2}+2\mu \left[v(r)\overline{E}_{nl}\right]\right]\overline{R}_{nl}(r)=0$$
(11)
outside the range of the potential ($`r>\lambda `$) is given by a superposition of spherical Bessel and Neumann functions

$$\overline{R}_{nl}(r)=\overline{A}_{nl}\left[j_l(\overline{q}_{nl}r)\mathrm{tan}\eta _l(\overline{q}_{nl})n_l(\overline{q}_{nl}r)\right],$$
(12)
where $`\eta _l(q)`$ denotes the phase shift of the potential for the $`l`$-th partial wave. The bar distinguishes quantities in presence of the interaction from those in the noninteracting case. Imposing the boundary condition (4) we get the following implicit equation for the momenta $`\overline{q}_{nl}`$ in the interacting case

$$\frac{j_l(\overline{q}_{nl}\mathrm{\Lambda })}{n_l(\overline{q}_{nl}\mathrm{\Lambda })}=\mathrm{tan}\eta _l(\overline{q}_{nl}).$$
(13)
Expressing the Bessel and Neumann functions by their asymptotic expansion (6) this reduces to

$$\mathrm{tan}(\overline{q}_{nl}\mathrm{\Lambda }\pi l/2)=\mathrm{tan}\eta _l(\overline{q}_{nl}).$$
(14)
In order to associate the energy levels with same quantum number $`n`$ in the interacting and noninteracting case as described above, the lowest positive-energy state should be labeled with the quantum number $`n=0`$. To achieve this for a potential with $`n_l^\text{b}`$ bound two-body states with angular momentum quantum number $`l`$ we add the phase $`\pi (n+n_l^\text{b})`$ to the argument on the r.h.s., i.e., the bound states contribute according to Levinson’s theorem. The momenta $`\overline{q}_{nl}`$ in presence of the interaction are thus determined by the equation

$$\overline{q}_{nl}\mathrm{\Lambda }=\eta _l(\overline{q}_{nl})+\pi (n+n_l^\text{b}+l/2).$$
(15)
The momentum shift $`\mathrm{\Delta }q_{nl}`$ induced by the interaction is obtained by subtracting the momenta of the noninteracting system (7) from equation (15)

$$\begin{array}{cc}\hfill \mathrm{\Delta }q_{nl}\mathrm{\Lambda }& =(\overline{q}_{nl}q_{nl})\mathrm{\Lambda }\hfill \\ & =[\eta _l(\overline{q}_{nl})\pi n_l^\text{b}]=\stackrel{~}{\eta }_l(\overline{q}_{nl}).\hfill \end{array}$$
(16)
Here $`\stackrel{~}{\eta }_l(q)=\eta _l(q)\pi n_l^\text{b}`$ denotes the phase shift reduced by the contribution of the bound states. In a final step we expand the phase shifts around the momenta $`q_{nl}`$ of the noninteracting system

$$\stackrel{~}{\eta }_l(\overline{q}_{nl})=\stackrel{~}{\eta }_l(q_{nl})+\stackrel{~}{\eta }_l^{}(q_{nl})\mathrm{\Delta }q_{nl}+\mathrm{}.$$
(17)
Already the term linear in $`\mathrm{\Delta }q_{nl}`$ can be neglected in good approximation because the momentum shift is of the order $`1/\mathrm{\Lambda }`$ according to equation (16). Thus we retain the following simple expression for the momentum shift

$$\mathrm{\Delta }q_{nl}\mathrm{\Lambda }=\stackrel{~}{\eta }_l(q_{nl}).$$
(18)
The shift of the energy levels of the interacting two-body system compared to the noninteracting spectrum is connected with the momentum shift by

$$\frac{\mathrm{\Delta }E_{nl}}{E_{nl}}=2\frac{\mathrm{\Delta }q_{nl}}{q_{nl}}+\left(\frac{\mathrm{\Delta }q_{nl}}{q_{nl}}\right)^2.$$
(19)
The term quadratic in the momentum shift can be neglected because $`\mathrm{\Delta }q_{nl}/q_{nl}`$ is small. The final expression for the energy shift reads

$$\frac{\mathrm{\Delta }E_{nl}}{E_{nl}}=\frac{2}{\mathrm{\Lambda }}\frac{\stackrel{~}{\eta }_l(q_{nl})}{q_{nl}}.$$
(20)
The proportionality between energy shift and phase shifts is well known Gottfried (1966) and was used in different applications before.

#### II.2.2 ECI in Scattering Length Approximation

In a second step we construct an operator form of the Effective Contact Interaction that obeys condition (1). For the application to ultracold dilute quantum gases the general form (20) of the energy shift $`\mathrm{\Delta }E_{nl}`$ can be simplified considerably. The relative momenta in these systems are extremely low, i.e., typical relative wave lengths are large compared to the range of the interaction. This allows an expansion of the phase shifts $`\stackrel{~}{\eta }_l(q)`$ in a power series in relative momentum. In lowest order approximation ($`qa_l1`$) the phase shifts of the $`l`$-th partial wave can be expressed in terms of the corresponding scattering length $`a_l`$ 111Some authors Bohn (2000) use a different definition of the p-wave scattering length without the $`l`$-dependent prefactors. We use this general form (see Huang and Yang (1957)) where the scattering length for a hard sphere equals the sphere radius for all $`l`$.

$$\frac{\stackrel{~}{\eta }_l(q)}{q^{2l+1}}\frac{\mathrm{tan}\eta _l(q)}{q^{2l+1}}\frac{(2l+1)}{[(2l+1)!!]^2}a_l^{2l+1}.$$
(21)
The energy shift in this scattering length approximation is given by

$$\frac{\mathrm{\Delta }E_{nl}}{E_{nl}}=\frac{2}{\mathrm{\Lambda }}\frac{(2l+1)}{[(2l+1)!!]^2}q_{nl}^{2l}a_l^{2l+1}.$$
(22)
Based on this form we can construct a simple operator for the ECI. For applications where the approximation (21) is not sufficient higher order terms of the power series can be included successively. We will come back to this point in the following section.

According to the dependence of the energy shifts on the angular momentum quantum number $`l`$ the operator of the effective interaction $`𝒗^{\text{eff}}`$ is formulated as a sum of independent operators $`𝒗_l^{\text{eff}}`$ for each partial wave

$$𝒗^{\text{eff}}=\underset{l=0}{\overset{\mathrm{}}{}}𝚷_l𝒗_l^{\text{eff}}𝚷_l,$$
(23)
where $`𝚷_l`$ denotes the projection operator on the subspace spanned by states of relative angular momentum $`l`$. We want to use a contact interaction for each partial wave. This requires nonlocal interaction operators beyond $`l=0`$, i.e., derivative couplings. The simplest ansatz for the operator of the effective contact interaction for the $`l`$-th partial wave is

$$\begin{array}{cc}\hfill 𝒗_l^{\text{eff}}& =(\stackrel{}{𝒒}\widehat{𝒓})^lg_l\delta ^{(3)}(\stackrel{}{𝒓})(\widehat{𝒓}\stackrel{}{𝒒})^l\hfill \\ & =\mathrm{d}^3r|\stackrel{}{r}>\frac{\stackrel{}{}^l}{r^l}g_l\delta ^{(3)}(\stackrel{}{r})\frac{\stackrel{}{}^l}{r^l}<\stackrel{}{r}|.\hfill \end{array}$$
(24)
Here $`\stackrel{}{𝒒}=\frac{1}{2}(\stackrel{}{𝒑}_1\stackrel{}{𝒑}_2)`$ denotes the operator of the relative momentum of two particles, $`\widehat{𝒓}=\stackrel{}{𝒓}/𝒓`$ the unit vector along the relative coordinate and $`\delta ^{(3)}(\stackrel{}{r})=\frac{1}{4\pi r^2}\delta (r)`$ the radial component of the 3-dimensional delta function. The arrows above the derivatives indicate to which side they act.

The interaction strength $`g_l`$ is a constant that contains the relevant information on the original two-body potential. The connection is provided by condition (1) via the shift of the energy levels with respect to the free spectrum. To evaluate (1) we calculate the expectation value of the contact interaction (24) for the $`l`$-th partial wave taken with the noninteracting two-body states $`|nlm>`$

$$\begin{array}{cc}& <nlm\left|𝒗_l^{\text{eff}}\right|nlm>\hfill \\ & =g_l\mathrm{d}^3r\delta ^{(3)}(\stackrel{}{r})\left|\frac{^l}{r^l}R_{nl}(r)Y_{lm}(\mathrm{\Omega })\right|^2\hfill \\ & =\frac{g_l}{4\pi }\left|\frac{^l}{r^l}R_{nl}(r)\right|_{r=0}^2\hfill \\ & =\frac{g_l}{4\pi }\left[\frac{l!}{(2l+1)!!}\right]^2A_{nl}^2q_{nl}^{2l},\hfill \end{array}$$
(25)
where we used the expansion of the radial wave-function around the origin

$$R_{nl}(r)=A_{nl}\frac{(q_{nl}r)^l}{(2l+1)!!}\left[1\frac{(q_{nl}r)^2}{2(2l+3)}+\mathrm{}\right].$$
(26)
Equating the expectation value (LABEL:eq:eci\_expectationvalue) with the energy shift results in the following equation for the interaction strengths

$$g_l=4\pi \left[\frac{(2l+1)!!}{l!}\right]^2A_{nl}^2\frac{\mathrm{\Delta }E_{nl}}{q_{nl}^{2l}}$$
(27)
Inserting the normalization constant (10) and the energy shift in scattering length formulation (22) gives the final expression for the interaction strengths

$$g_l=\frac{4\pi }{2\mu }\frac{(2l+1)}{(l!)^2}a_l^{2l+1}.$$
(28)
Together with equation (24) this defines the scattering length formulation of the Effective Contact Interaction. Notice that these expressions do not depend on the auxiliary boundary condition (4), which was introduced to obtain a discrete energy spectrum.

For $`l=1`$ sometimes a gradient operator with respect to the relative coordinate is used instead of the radial derivative (24). This alternative ansatz reads

$$\begin{array}{cc}\hfill 𝒗_1^{\text{eff}}& =\stackrel{}{𝒒}\stackrel{~}{g}_1\delta ^{(3)}(\stackrel{}{𝒓})\stackrel{}{𝒒}\hfill \\ & =\mathrm{d}^3r|\stackrel{}{r}>\stackrel{}{}\stackrel{~}{g}_1\delta ^{(3)}(\stackrel{}{r})\stackrel{}{}<\stackrel{}{r}|.\hfill \end{array}$$
(29)
It should be noted that the p-wave interaction strength $`\stackrel{~}{g}_1`$ in the gradient formulation is related to the interaction strength (28) by

$$\stackrel{~}{g}_1=g_1/3.$$
(30)
#### II.2.3 Beyond Scattering Length Approximation

The operator of the Effective Contact Interaction can be generalized systematically beyond the scattering length approximation shown in the preceding section. A formal scheme emerges from the expansion of the phase shifts in a power series in $`q`$

$$\frac{\stackrel{~}{\eta }_l(q)}{q^{2l+1}}=\underset{\nu =0}{\overset{\mathrm{}}{}}\frac{1}{\nu !}c_l^{(\nu )}q^{2\nu }.$$
(31)
The momentum independent lowest order term of this expansion matches the scattering length approximation. From equation (21) the connection between the coefficient $`c_l^{(0)}`$ and the the scattering length $`a_l`$ becomes obvious

$$c_l^{(0)}=\frac{(2l+1)}{[(2l+1)!!]^2}a_l^{2l+1}.$$
(32)
The coefficient $`c_l^{(1)}`$ of the quadratic term of the expansion is connected to the so called effective range or effective volume of the potential. We will discuss this contribution in more detail later. By inserting the expansion (31) into the general formula (20) for the energy shifts we obtain

$$\frac{\mathrm{\Delta }E_{nl}}{E_{nl}}=\frac{2}{\mathrm{\Lambda }}\underset{\nu =0}{\overset{\mathrm{}}{}}\frac{1}{\nu !}c_l^{(\nu )}q_{nl}^{2l+2\nu }.$$
(33)
Following the basic concept of the ECI these energy shifts have to be generated by the operator of the ECI according to condition (1). To include the momentum dependence of the energy shifts we have to generalize the ansatz for the ECI operator compared to the simple momentum-independent scattering-length formulation (24). The further calculation will show that

$$\begin{array}{cc}\hfill 𝒗_l^{\text{eff}}=\underset{\nu =0}{\overset{\mathrm{}}{}}\frac{1}{2}& g_l^{(\nu )}[(\stackrel{}{𝒒}\widehat{𝒓})^l\delta ^{(3)}(\stackrel{}{𝒓})(\widehat{𝒓}\stackrel{}{𝒒})^{l+2\nu }\hfill \\ & +(\stackrel{}{𝒒}\widehat{𝒓})^{l+2\nu }\delta ^{(3)}(\stackrel{}{𝒓})(\widehat{𝒓}\stackrel{}{𝒒})^l]\hfill \end{array}$$
(34)
is a proper ansatz for the effective interaction operator for the $`l`$-th partial wave. Besides the more complex nonlocal structure a set of interaction strengths $`g_l^{(\nu )}`$ ($`\nu =0,1,\mathrm{}`$) for each partial wave is included. They are related to the coefficients $`c_l^{(\nu )}`$ and thus correspond to the different powers of the momentum in equation (33). To employ condition (1) we calculate the expectation value of $`𝒗_l^{\text{eff}}`$ with noninteracting two-body states $`|nlm>`$

$$\begin{array}{cc}& <nlm\left|𝒗_l^{\text{eff}}\right|nlm>\hfill \\ & =\frac{1}{4\pi }\underset{\nu =0}{\overset{\mathrm{}}{}}g_l^{(\nu )}\left[\frac{^l}{r^l}R_{nl}(r)\right]_{r=0}\left[\frac{^{l+2\nu }}{r^{l+2\nu }}R_{nl}(r)\right]_{r=0}\hfill \\ & =\frac{A_{nl}^2}{4\pi }\underset{\nu =0}{\overset{\mathrm{}}{}}g_l^{(\nu )}\frac{l!}{(2l+1)!!}\frac{(1)^\nu (l+2\nu )!}{2^\nu \nu !(2l+2\nu +1)!!}q_{nl}^{2l+2\nu },\hfill \end{array}$$
(35)
where the full expansion of the noninteracting radial wave function (5) around $`r=0`$ was used Abramowitz and Stegun (1972)

$$R_{nl}(r)=A_{nl}\underset{\mu =0}{\overset{\mathrm{}}{}}\frac{(1)^\mu }{2^\mu \mu !(2l+2\mu +1)!!}(q_{nl}r)^{l+2\mu }.$$
(36)
By inserting the expansion of the energy shifts (33) and the expectation value into condition (1) and comparing the coefficients for the different powers of the momentum $`q_{nl}`$ we obtain

$$g_l^{(\nu )}=(1)^{\nu +1}\frac{4\pi }{2\mu }\frac{2^\nu (2l+1)!!(2l+2\nu +1)!!}{l!(l+2\nu )!}c_l^{(\nu )}.$$
(37)
Thus the interaction strengths $`g_l^{(\nu )}`$ of the general operator form of the ECI (34) are proportional to the coefficients $`c_l^{(\nu )}`$ of the expansion of the phase shifts (31). The equations (34) and (37) define the most general form of the Effective Contact Interaction.

For the application to dilute degenerate Fermi gases we will use the ECI up to quadratic terms in the momentum, i.e., we include the scattering length term of the s- and p-wave part as well as the s-wave effective range correction. At this point we have to discuss the connection between the quadratic term of the expansion (31) and the usual effective range theory. For the s-wave phase-shifts the effective range expansion reads

$$q\mathrm{cot}\eta _0(q)\frac{1}{a_0}+\frac{1}{2}r_0q^2,$$
(38)
where $`r_0`$ is the effective range of the potential. If we convert this into an expression for $`\stackrel{~}{\eta }_0(q)/q`$ and expand in $`q`$ we obtain

$$\frac{\stackrel{~}{\eta }_0(q)}{q}=a_0b_0q^2+\mathrm{}$$
(39)
with an effective volume $`b_0`$ that depends on the scattering length $`a_0`$ and the effective range $`r_0`$

$$b_0=\frac{1}{2}a_0^2r_0\frac{1}{3}a_0^3.$$
(40)
Rather than using this relation we will adjust $`b_0`$ in order to get the best representation of the phase shift $`\stackrel{~}{\eta }_0(q)/q`$ with the truncated expansion (39).

Finally, inserting $`c_0^{(1)}=b_0`$ into equation (37) gives an expression for the interaction strength of the s-wave effective range term

$$g_0^{(1)}=\frac{12\pi }{2\mu }b_0.$$
(41)
### II.3 Example: Square-Well Potential

To illustrate the concept of the Effective Contact Interaction we use the simple toy-problem of two particles interacting by an attractive square-well potential of radius $`\lambda `$ and depth $`V_0`$.

First we look at typical wave functions to which the idea of the ECI applies. The major condition is that the typical wavelength of the relative motion is large compared to the range of the interaction. This is ensured by choosing the radius $`\lambda `$ of the square-well much smaller than the radius $`\mathrm{\Lambda }`$ associated with the boundary condition (4); in the following we use $`\lambda =0.01\mathrm{\Lambda }`$. Figure 2 shows the radial wave function of the lowest $`l=0`$ state with positive energy for different potential depths $`V_0`$. Outside the potential the structure of the wave functions is very similar for the different interaction strengths. Only the wavelength is changed slightly due to the different matching to the wave function in the interior (see insert of Figure 2). This change of the relative momentum translates immediately into an energy shift. From that picture the connection between energy shifts and phase shifts is evident.

A second point becomes clear from this simple example: The detailed structure of the radial dependence of the potential or the number of bound states is irrelevant for the energy shift, only the phase shifts $`\stackrel{~}{\eta }_l(q)`$ matter. The insert in Figure 2 shows that the wave functions behave very different within the range $`\lambda `$ due to the different potential depths. Moreover the potentials have a different number of bound states, e.g., the thick solid curve is associated with a potential with one bound state but zero phase shift. Hence the behavior outside the potential is identical to the noninteracting case (thin solid curve) and the energy shift is zero.

Next we investigate the dependence of the $`l=0`$ energy shifts on the strength of the attractive potential. Figure 3 shows the relative energy shift $`\mathrm{\Delta }E^{\text{exact}}/E`$ versus $`\lambda \sqrt{V_0}`$, where the radius $`\lambda =0.01\mathrm{\Lambda }`$ of the square-well potential is kept fixed and the depth $`V_0`$ is increased. A characteristic pattern appears: In the vicinity of interaction strengths where the potential gains another bound state the relative energy shift assumes large positive and negative values. Large positive energy shifts occur for potentials that have a very weakly bound state, negative energy shifts for those that have an almost bound state. In between these interaction strengths extended plateaus of nearly constant energy shift appear. Within the plateaus the energy shift is independent of the radial quantum number $`n`$ of the level or the relative momentum. This is a special property of the s-wave channel; for higher partial waves the relative energy shift (22) is proportional to $`q^{2l}`$. Only at the edges of the plateaus a slight dependence on the relative momentum shows up (see Figure 3).

This structure is closely related to the behavior of the s-wave scattering length. The $`l=0`$ energy shift induced by the ECI in scattering length approximation (22) is proportional to the scattering length $`a_0`$. For the square-well potential we get

$$\frac{\mathrm{\Delta }E_{n0}}{E_{n0}}=2\frac{a_0}{\mathrm{\Lambda }}\text{with}\frac{a_0}{\lambda }=1\frac{\mathrm{tan}(\lambda \sqrt{V_0})}{\lambda \sqrt{V_0}}.$$
(42)
This ECI energy shift is right on top of the solid curve in Figure 3, i.e., it agrees very well with the exact energy shift for the lowest positive energy state. Even for higher momenta the agreement is very good provided that the magnitude of the scattering length is not too large. Significant deviations occur only if momentum *and* scattering length are large.

To obtain a quantitative measure of for the applicability of the ECI in scattering length approximation we investigate the relative deviation $`(\overline{E}\overline{E}^{\text{exact}})/E`$ of the ECI energy levels compared to the exact ones for the square-well potential. We expect that the agreement gets worse if either the relative momentum or if the scattering length is large. Therefore, Figure 4 shows the relative energy deviation for $`l=0`$ and $`l=1`$ states as function of the product of momentum and scattering length, $`qa_l`$, which was assumed to be small in order to introduce the scattering length (21). The different curves correspond to different values of the momentum $`q_{nl}`$ and were obtained by varying the depth of the square-well potential and thus the scattering length.

As expected the deviation increases with increasing value of $`qa_l`$. Nevertheless, the deviation of the energy calculated with the ECI in scattering length approximation is below 1% up to rather large values of $`qa_l1`$. If we tolerate a maximum deviation of 5% then the scattering length formulation can be used up to values $`qa_l1.5`$.

It should be noted that the relative deviation of the general form (20) of the ECI in the parameter range discussed above is below $`10^4`$. Thus all approximations made to obtain equation (20) are valid on a high level of accuracy. The restrictions on the validity of the scattering length formulation (22) originate for the replacement of the phase shifts by the scattering length alone, which is no inherent part of th ECI concept. If the simple scattering length formulation is not sufficient for a special application one can go beyond that.

For example, the inclusion of effective volume corrections (see eq. (39)) improves the agreement with the exact energy shifts. In this way we can reduce the maximum deviation to only 1% up to $`qa_01.5`$.

### II.4 ECI versus Pseudopotential

The idea to simulate the effect of a complicated finite-range two-body potential by a simple s-wave contact interaction dates back to E. Fermi Fermi (1936) and was used by several authors Blatt and Weisskopf (1979) in various physical contexts. K. Huang and C.N. Yang Huang and Yang (1957); Huang (1963) generalized this idea and constructed the so called *pseudopotential* that acts in all partial waves.

The aim of the pseudopotential is (a) to generate the phase shifts of the original potential by a boundary condition at $`r=0`$ and (b) to reformulate this by an additional inhomogenous term in the Schrödinger equation of the two-body scattering problem. This additional term is interpreted as the pseudopotential, which can be phrased in the following operator form Roth (2000)

$$𝒗_l^{\text{pseudo}}=\mathrm{d}^3r|\stackrel{}{r}>\frac{1}{r^l}g_l^{\text{pseudo}}\delta ^{(3)}(\stackrel{}{r})\frac{\stackrel{}{}^{2l+1}}{r^{2l+1}}r^{l+1}<\stackrel{}{r}|$$
(43)
with an interaction strength

$$g_l^{\text{pseudo}}=\frac{4\pi }{2\mu }\frac{(l+1)}{(2l+1)!}a_l^{2l+1}$$
(44)
for the $`l`$-th partial wave. For this discussion we restrict ourselves to the scattering length approximation of the phase shifts. Due to the fact that the radial derivative acts only to the right hand side the operator of the pseudopotential (43) is not hermitian. This is in contradiction to the basic concept of effective interactions.

A more severe weakness shows up when we evaluate the energy shifts induced by the pseudopotential. As discussed in section II.1 the expectation value of a proper effective interaction with eigenstates $`|nlm>`$ of the noninteracting two-body system should be equal to the energy shift induced by the original potential. For the pseudopotential (43) we obtain an energy shift

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }E_{nl}^{\text{pseudo}}}{E_{nl}}& =\frac{1}{E_{nl}}<nlm\left|𝒗_l^{\text{pseudo}}\right|nlm>\hfill \\ & =\frac{2}{\mathrm{\Lambda }}\frac{(l+1)}{[(2l+1)!!]^2}q_{nl}^{2l}a_l^{2l+1}.\hfill \end{array}$$
(45)
This has to be compared with the full energy shift in scattering length approximation (22) which is by construction reproduced by the ECI. Obviously the energy shift induced by the pseudopotential for states with $`l>0`$ is by a factor $`(l+1)/(2l+1)`$ smaller the the energy shift of the original potential. Thus the pseudopotential underestimates the effect of the two-body interactions beyond s-wave when used in a mean-field framework. For the widely used s-wave part the energy shifts of the pseudopotential agree with the energy shifts of the original potential.

We conclude that the nonhermitian pseudopotential is not a proper effective interaction for a mean-field description of a dilute quantum gases that goes beyond s-wave interactions.

## III Energy Functional of a Trapped Multi-Component Fermi Gas

### III.1 Fundamentals

In the following we investigate the ground state properties of a dilute Fermi gas composed of $`\mathrm{\Xi }`$ distinguishable components that is trapped in an external potential $`U(\stackrel{}{x})`$ at temperature $`T=0`$K. In present experiments DeMarco and Jin (1999) one or two components are used, which belong to the same atomic species but are distinguished by different projections $`M_F`$ of the total angular momentum $`F`$ onto the direction of an external magnetic field. We distinguish the different components by a formal quantum number $`\xi =1,\mathrm{},\mathrm{\Xi }`$. For simplicity we use the same mass $`m`$ of the atoms for all components.

We treat the many-body problem in the framework of density functional theory and construct an energy functional of the inhomogenous multi-component Fermi gas within a proper approximation. The ground state density-distribution of the many-body system is then determined by functional minimization of the energy.

The large particle numbers of the order $`N10^6`$ allow the rather simple Thomas-Fermi approximation for the energy functional. It is assumed that the energy density of the inhomogenous system is described locally by the energy density of the corresponding homogenous system, and higher-order terms which include gradients of the density are small. To check the quality of the Thomas-Fermi approximation we calculated the next order gradient corrections for a trapped noninteracting Fermi gas. For $`N=100`$ particles the relative contribution of the gradient correction to the total energy is of the order of $`10^2`$; for typical particle numbers of $`N=10^6`$ it drops to $`10^5`$ Roth (2000).

As starting point for the Thomas-Fermi approximation we calculate the energy density of the homogenous interacting multi-component Fermi gas in mean-field approximation. The basic restriction of the mean-field picture is that two- and many-body correlations induced by the interaction are not contained in the many-body state. Nevertheless they can be implemented implicitly by using a proper effective interaction that is tailored for the model-space available. In the previous section we constructed the Effective Contact Interaction especially for the mean-field description of dilute not self-bound quantum gases.

A central topic of the following studies is the role of the interaction on the structure and stability of trapped degenerate Fermi gases. Our special interest concerns the p-wave part of the interaction, which contributes even at $`T=0`$K — in contrast to bosonic systems. It will turn out that the p-wave terms can be of substantial importance for the ground state properties of fermionic systems and should not be neglected from the outset.

We write the Hamilton operator of the system as a sum of the external trapping potential $`𝑼`$ and an internal part $`𝑯_{\text{int}}`$

$$𝑯=𝑼+𝑯_{\text{int}}.$$
(46)
The internal part contains the kinetic energy and the Effective Contact Interaction as discussed in section II. We include the s-wave and p-wave terms of the ECI in scattering length formulation as well as the s-wave effective range correction. With equations (24), (28), (34), and (41) the internal Hamiltonian reads

$`𝑯_{\text{int}}`$ $`={\displaystyle \frac{1}{2m}}`$ $`{\displaystyle \underset{i}{}}\stackrel{}{𝒑}_i^2`$
$`+{\displaystyle \frac{4\pi }{m}}a_0`$ $`{\displaystyle \underset{i,j>i}{}}\delta ^{(3)}(\stackrel{}{𝒓}_{ij})`$ (47)
$`{\displaystyle \frac{12\pi }{m}}b_0`$ $`{\displaystyle \underset{i,j>i}{}}{\displaystyle \frac{1}{2}}\left[\delta ^{(3)}(\stackrel{}{𝒓}_{ij})(\widehat{𝒓}_{ij}\stackrel{}{𝒒}_{ij})^2+\text{h.a.}\right]`$
$`+{\displaystyle \frac{12\pi }{m}}a_1^3`$ $`{\displaystyle \underset{i,j>i}{}}(\stackrel{}{𝒒}_{ij}\widehat{𝒓}_{ij})\delta ^{(3)}(\stackrel{}{𝒓}_{ij})(\widehat{𝒓}_{ij}\stackrel{}{𝒒}_{ij}).`$
The summations over the particle indices $`i`$ and $`j`$ range from $`1`$ to the total number of particles. The properties of the two-body interaction are parameterized by the s- and p-wave scattering lengths $`a_0`$ and $`a_1`$, respectively, and by the s-wave effective volume $`b_0`$. In general the interaction parameters depend on the component quantum numbers $`\xi `$ of the interacting particles. In order to discuss the basic phenomena we restrict ourselves to equal interaction parameters for all components. The generalization to scattering length matrices that account for the dependence on the component indices of the two interacting particles is straightforward.

Experimentally each component may experience a different trapping potential $`U_\xi (\stackrel{}{x})`$. For magnetic traps this is due to the different magnetic momenta of the components, which leads to a relative shift of the trapping potentials for the components. Thus the operator of the external potential has the following form

$$𝑼=\underset{i}{}\underset{\xi }{}U_\xi (\stackrel{}{𝒙}_i)𝚷_{\xi ,i}$$
(48)
where $`𝚷_\xi `$ is a projection operator onto states with the component quantum number $`\xi `$.

Many of the results shown in the next sections do not depend on the actual shape of the trapping potential. If the shape enters explicitly we assume a deformed harmonic oscillator potential

$$\begin{array}{cc}\hfill U(\stackrel{}{x})& =\frac{m\omega ^2}{2}\left(\lambda _1^2x_1^2+\lambda _2^2x_2^2+\lambda _3^2x_3^2\right)\hfill \\ & =\frac{1}{2m\mathrm{}^4}\left(\lambda _1^2x_1^2+\lambda _2^2x_2^2+\lambda _3^2x_3^2\right),\hfill \end{array}$$
(49)
where $`\omega =\sqrt[3]{\omega _1\omega _2\omega _3}`$ is the mean oscillator frequency and $`\mathrm{}=(m\omega )^{1/2}`$ the corresponding mean oscillator length, i.e., the mean width of the Gaussian single-particle ground state of the harmonic oscillator potential. The deformation is parameterized by the ratios $`\lambda _i=\omega _i/\omega `$, which fulfill the condition $`\lambda _1\lambda _2\lambda _3=1`$.

### III.2 Energy Density in Thomas-Fermi Approximation

The calculation of the energy density functional of the inhomogenous interacting Fermi gas is performed in two steps: First we calculate the energy density of the corresponding homogenous system in mean-field approximation. In the second step this is translated into an energy density of the inhomogenous system by means of the Thomas-Fermi approximation.

The ground state of a many-fermion system in mean-field approximation is given by an antisymmetrized product of one-body states $`|i>`$. In the case of a homogenous system the one-body states are eigenstate of the momentum operator with eigenvalues $`\stackrel{}{k}_i`$. In addition they are characterized by the component quantum number $`\xi `$

$$|i>=|\stackrel{}{k}_i>|\xi _i>.$$
(50)
Assuming a box of volume $`V`$ with periodic boundary conditions the spatial part of the one-body states is given by

$$<\stackrel{}{x}|\stackrel{}{k}_i>=\frac{1}{\sqrt{V}}\mathrm{exp}\left(\mathrm{i}\stackrel{}{k}_i\stackrel{}{x}\right).$$
(51)
The energy density of the homogenous system is given by the expectation value of the internal part of the Hamilton operator (III.1)

$$_{\text{hom}}=<𝑯_{\text{int}}>/V.$$
(52)
The calculation of the expectation values of the several parts of the Hamiltonian is straightforward Roth (2000). As a function of the Fermi momenta $`\kappa _\xi `$ of the different components the energy density reads

$`_{\text{hom}}(\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }`$ $`)`$
$`={\displaystyle \frac{1}{20\pi ^2m}}`$ $`{\displaystyle \underset{\xi }{}}\kappa _\xi ^5`$
$`+{\displaystyle \frac{a_0}{9\pi ^3m}}`$ $`{\displaystyle \underset{\xi ,\xi ^{}>\xi }{}}\kappa _\xi ^3\kappa _\xi ^{}^3`$ (53)
$`+{\displaystyle \frac{a_1^3}{30\pi ^3m}}`$ $`{\displaystyle \underset{\xi }{}}\kappa _\xi ^8`$
$`+{\displaystyle \frac{a_1^3+b_0}{60\pi ^3m}}`$ $`{\displaystyle \underset{\xi ,\xi ^{}>\xi }{}}[\kappa _\xi ^3\kappa _\xi ^{}^5+\kappa _\xi ^5\kappa _\xi ^{}^3].`$
The summations run over all components $`\xi =1,\mathrm{},\mathrm{\Xi }`$. To avoid fractional exponents we use Fermi momenta $`\kappa _\xi `$ rather than densities $`\rho _\xi =\kappa _\xi ^3/(6\pi ^2)`$.

The basic assumption of the Thomas-Fermi (or local-density) approximation is that the energy density of the inhomogenous Fermi gas is locally given by the energy density of the corresponding homogenous system. Thus the energy density of the inhomogenous system is constructed from (III.2) by replacing $`\kappa _\xi `$ with *local Fermi momenta* $`\kappa _\xi (\stackrel{}{x})`$. In addition the contribution of the external trapping potential has to be included. This results in the following expression for the energy density of the trapped interacting multi-component Fermi gas

$`[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }](\stackrel{}{x})`$
$`={\displaystyle \frac{1}{6\pi ^2}}`$ $`{\displaystyle \underset{\xi }{}}U_\xi (\stackrel{}{x})\kappa _\xi ^3(\stackrel{}{x})`$
$`+{\displaystyle \frac{1}{20\pi ^2m}}`$ $`{\displaystyle \underset{\xi }{}}\kappa _\xi ^5(\stackrel{}{x})`$
$`+{\displaystyle \frac{a_0}{9\pi ^3m}}`$ $`{\displaystyle \underset{\xi ,\xi ^{}>\xi }{}}\kappa _\xi ^3(\stackrel{}{x})\kappa _\xi ^{}^3(\stackrel{}{x})`$ (54)
$`+{\displaystyle \frac{a_1^3}{30\pi ^3m}}`$ $`{\displaystyle \underset{\xi }{}}\kappa _\xi ^8(\stackrel{}{x})`$
$`+{\displaystyle \frac{a_1^3+b_0}{60\pi ^3m}}`$ $`{\displaystyle \underset{\xi ,\xi ^{}>\xi }{}}[\kappa _\xi ^3(\stackrel{}{x})\kappa _\xi ^{}^5(\stackrel{}{x})+\kappa _\xi ^5(\stackrel{}{x})\kappa _\xi ^{}^3(\stackrel{}{x})].`$
The local Fermi momentum is related to the density of particles of the component $`\xi `$ by

$$\rho _\xi (\stackrel{}{x})=\frac{1}{6\pi ^2}\kappa _\xi ^3(\stackrel{}{x}).$$
(55)
Accordingly the number of particles of component $`\xi `$ is given by

$$N_\xi =\mathrm{d}^3x\rho _\xi (\stackrel{}{x})=\frac{1}{6\pi ^2}\mathrm{d}^3x\kappa _\xi ^3(\stackrel{}{x}).$$
(56)
As discussed in section II.3 we can reproduce the two-body energy spectrum with an accuracy of about 5% up to $`a_lq1.5`$. If we take this as a limit for the root mean square of the relative momentum $`𝒒^2^{1/2}=0.53\kappa `$ in the many-body system, we can apply our many-body model up to $`a_l\kappa 3`$.

### III.3 Functional Variation and the Extremum Condition

The ground state density of a system is found by minimizing the energy functional

$$E[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }]=\mathrm{d}^3x[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }](\stackrel{}{x}).$$
(57)
for given particle numbers $`N_\xi `$. This constraint is implemented with help of a set of Lagrange multipliers $`\mu _\xi `$, which are the chemical potentials of the different components. The Legendre transformed functional

$$\begin{array}{cc}& F[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }]\hfill \\ & =E[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }]\underset{\xi }{}\mu _\xi N_\xi \hfill \\ & =\mathrm{d}^3x[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }](\stackrel{}{x})\underset{\xi }{}\frac{\mu _\xi }{6\pi ^2}\kappa _\xi ^3(\stackrel{}{x})\hfill \\ & =\mathrm{d}^3x[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }](\stackrel{}{x}),\hfill \end{array}$$
(58)
has to be minimized by functional variation. A necessary but not sufficient condition for a set of local Fermi momenta $`\{\kappa _1(\stackrel{}{x}),\mathrm{},\kappa _\mathrm{\Xi }(\stackrel{}{x})\}`$ to minimize the transformed energy functional $`F[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }]`$ is stationarity, i.e., that the first variation of $`F[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }]`$ with respect to all $`\kappa _\xi (\stackrel{}{x})`$ vanishes

$$\frac{\delta }{\delta \kappa _\xi }F[\kappa _1,\mathrm{},\kappa _\xi ]=0\text{for all}\xi .$$
(59)
This extremum condition is fulfilled if the derivative of the integrand $`[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }](\stackrel{}{x})`$ with respect to all local Fermi momenta vanishes at each point $`\stackrel{}{x}`$

$$\frac{}{\kappa _\xi (\stackrel{}{x})}[\kappa _1,\mathrm{},\kappa _\mathrm{\Xi }](\stackrel{}{x})=0\text{for all}\stackrel{}{x},\xi .$$
(60)
Inserting expression (III.2) for the energy density and evaluating the derivative results in the extremum condition

$$\begin{array}{cc}& m[\mu _\xi U_\xi (\stackrel{}{x})]\hfill \\ & =\frac{1}{2}\kappa _\xi ^2(\stackrel{}{x})+\frac{2a_0}{3\pi }\underset{\xi ^{}\xi }{}\kappa _\xi ^{}^3(\stackrel{}{x})+\frac{8a_1^3}{15\pi }\kappa _\xi ^5(\stackrel{}{x})\hfill \\ & +\frac{a_1^3+b_0}{30\pi }\underset{\xi ^{}\xi }{}[3\kappa _\xi ^{}^5(\stackrel{}{x})+5\kappa _\xi ^2(\stackrel{}{x})\kappa _\xi ^{}^3(\stackrel{}{x})]\hfill \end{array}$$
(61)
for all $`\stackrel{}{x}`$ and each $`\xi `$. This is a coupled set of $`\mathrm{\Xi }`$ polynomial equations for the local Fermi momenta $`\{\kappa _1(\stackrel{}{x}),\mathrm{},\kappa _\mathrm{\Xi }(\stackrel{}{x})\}`$ at some given point $`\stackrel{}{x}`$. Note that the trivial solutions $`\kappa _\xi (\stackrel{}{x})=0`$ were separated already. Any real solution of the extremum condition (LABEL:eq:ef\_extremumcond) corresponds to a stationary point of the energy functional. In general one has to check explicitly whether they correspond to a minimum of the energy functional or whether they are maxima or saddle points.

All following investigations on the structure and stability of degenerate Fermi gases and on the influence of s- and p-wave interactions are based on the extremum condition (LABEL:eq:ef\_extremumcond). Many physical conclusions can be drawn from its algebraic structure already. We will discuss these questions in detail for the one- and two-component Fermi gas in section IV and V, respectively.

## IV Single-Component Fermi Gas

As a first application of the formalism developed in the preceding section we study the properties of a degenerate single-component Fermi gas.

### IV.1 Effect of the p-Wave Interaction

The energy density of the interacting multi-component Fermi gas (III.2) reduces for the single-component system to the form

$$\begin{array}{cc}\hfill [\kappa ](\stackrel{}{x})& =\frac{1}{6\pi ^2}U(\stackrel{}{x})\kappa ^3(\stackrel{}{x})+\frac{1}{20\pi ^2m}\kappa ^5(\stackrel{}{x})\hfill \\ & +\frac{a_1^3}{30\pi ^3m}\kappa ^8(\stackrel{}{x}),\hfill \end{array}$$
(62)
where $`\kappa (\stackrel{}{x})`$ is the local Fermi momentum. The first term is the contribution of the trapping potential $`U(\stackrel{}{x})`$, the second term is the kinetic energy, and the third term describes the contribution of the p-wave interaction with a p-wave scattering length $`a_1`$. As mentioned earlier the s-wave part of the interaction does not contribute in a system of identical fermions due to the Pauli principle. Therefore the p-wave part is the leading interaction term and there is no reason to neglect it from the outset.

A first hint on the effects of the p-wave interaction are given by the density distributions for different values of the p-wave scattering length. The density distribution is obtained by solution of the extremum condition (LABEL:eq:ef\_extremumcond), which takes the simple form

$$m[\mu U(\stackrel{}{x})]=f_1[\kappa (\stackrel{}{x})]$$
(63)
$$\text{with}f_1(\kappa )=\frac{1}{2}\kappa ^2+\frac{8a_1^3}{15\pi }\kappa ^5.$$

This 5th order polynomial equation for the local Fermi momentum $`\kappa (\stackrel{}{x})`$ is solved numerically for each point $`\stackrel{}{x}`$. The chemical potential $`\mu `$ is adjusted such that particle number (56) assumes the desired value.

Figure 5 shows the resulting radial density profiles $`\rho (x)=\kappa ^3(x)/(6\pi ^2)`$ for a single-component gas of $`N=10^6`$ particles in a spherical trap with oscillator length $`\mathrm{}`$ for different p-wave scattering lengths $`a_1`$. The oscillator length defines the fundamental length scale of the problem and the parameter that determines the strength of the interaction is the ratio of the p-wave scattering length and oscillator length, $`a_1/\mathrm{}`$. To increase the magnitude of this ratio experimentally one can either increase the magnitude of the scattering length or decrease the oscillator length.

For a repulsive p-wave interaction, i.e., $`a_1/\mathrm{}>0`$, of increasing strength (dotted curves) the density distribution flattens and expands radially compared to the noninteracting system (solid line). For a ratio $`a_1/\mathrm{}=0.1`$ the central density has dropped to one half of the density of the noninteracting gas. With a typical experimental oscillator length of $`\mathrm{}=1\mu m`$ this ratio corresponds to a rather large scattering length of $`a_12000a_\text{B}`$, which nevertheless may be within the range of experimental parameters Bohn (2000). For a tightly confining trap with $`\mathrm{}=0.1\mu m`$ a moderate scattering length of $`a_1200a_\text{B}`$ is required to obtain the same ratio.

For an attractive p-wave interaction, $`a_1/\mathrm{}<0`$, the central density increases significantly with increasing interaction strength. If the central density exceeds a certain value, i.e., if $`|a_1/\mathrm{}|`$ exceeds a critical value, then the extremum condition (63) has no real solution any more. Physically this corresponds to a collapse of the dilute gas caused by the attractive mean-field that is generated by the p-wave interaction. We will discuss this question in detail in the following sections.

The dependence of the density distribution on the p-wave scattering length as depicted in Figure 5 already demonstrates that the p-wave interaction may have strong influence on the properties of degenerate Fermi gases.

### IV.2 Mean-Field Instability: A Variational Picture

To illustrate the origin and mechanism of the collapse of the metastable state of the trapped Fermi gas we utilize a simple variational picture. Assume a single-component Fermi gas of $`N`$ particles in a spherical symmetric oscillator potential. The local Fermi momentum of the interacting system is parameterized by the analytic expression for the local Fermi momentum of the noninteracting system

$$\kappa (\stackrel{}{x})=\frac{2(6N)^{1/3}}{X_\text{t}}\sqrt{1\left(\frac{x}{X_\text{t}}\right)^2}\text{for}xX_\text{t},$$
(64)
where the classical turning point $`X_\text{t}`$ is treated as variational parameter. By inserting this parameterization into the energy density (62) and integrating we obtain a closed expression for the energy as function of the parameter $`X_\text{t}`$

$$E(X_\text{t})=C_u\frac{NX_\text{t}^2}{\mathrm{}^4}+C_t\frac{N^{5/3}}{X_\text{t}^2}+C_1\frac{N^{8/3}a_1^3}{X_\text{t}^5},$$
(65)
with constant coefficients

$$C_u=\frac{3}{16m},C_t=\frac{3(9/2)^{1/3}}{2m},C_1=\frac{8^6(4/3)^{1/3}}{1925\pi ^2m}.$$
(66)
Again the first term corresponds to the external potentials, the second to the kinetic energy, and the third term to the p-wave interaction.

Figure 6 shows the dependence of the total energy (65) on the parameter $`X_\text{t}`$ for a system of $`N=10^6`$ particles with different p-wave scattering lengths. For attractive p-wave interactions, i.e., negative scattering length $`a_1`$, the contribution of the interaction in (65) is negative. Due to its $`X_\text{t}^5`$ dependence this interaction contribution leads to a rapid drop of the energy for systems of decreasing spatial extension and thus increasing density. At very high densities (small $`X_\text{t}`$) one formally ends up with states of negative energy, i.e., bound states. One should however keep in mind that the assumptions made for the construction of the Effective Contact Interaction are not valid in this high density regime any more.

The not self-bound metastable state appears as local minimum at positive energies and low densities provided the p-wave attraction is sufficiently weak; the thick solid curve in Figure 6 shows an example. If the strength of the attractive p-wave interaction increases then the local minimum flattens and devolves to a saddle point (dashed curve). From this particular interaction strength on the metastable low-density state does not exist anymore, only the true ground state of the system remains, which is usually a crystal. The system collapses if the barrier caused by the positive kinetic and the attractive mean-field energy vanishes. Since the mean-field attraction grows with increasing density the system is unstable and collapses towards a high-density configuration.

### IV.3 Mean-Field Instability: Stability Conditions

Based on the extremum condition (63) we derive a set of analytic stability conditions that relate the maximum density of a metastable system with the p-wave scattering length. Part of this was already discussed in Roth and Feldmeier (2000).

The mean-field instability of the system occurs at values of the chemical potential $`\mu `$ and the scattering length $`a_1`$, where the extremum condition (63) does not have a real solution any more. This is shown in a pictorial way in Figure 7 where the the right hand side $`f_1(\kappa )`$ of the extremum condition is plotted as function of $`\kappa `$ for different $`a_1/\mathrm{}`$. The solution of the extremum condition at some specific point $`\stackrel{}{x}`$ is given by the value of $`\kappa `$ at which the respective curve reaches the value $`m[\mu U(\stackrel{}{x})]`$. In the minimum of the trapping potential (we assume $`U(\stackrel{}{x})=0`$ in the minimum) the solution is given by the point where $`f_1(\kappa )`$ reaches the value $`m\mu `$. By moving towards the outer regions of the trap $`m[\mu U(\stackrel{}{x})]`$ decreases and one scans $`f_1(\kappa )`$ down along ordinate until one reaches $`m[\mu U(\stackrel{}{x})]=0`$, i.e., the classical turning point.

For repulsive p-wave interactions (dotted curve) the r.h.s. of (63) is a monotonic growing function and solutions exist for arbitrary values of $`m[\mu U(\stackrel{}{x})]`$. If the scattering length $`a_1`$ is negative (dashed curve) then $`f_1(\kappa )`$ exhibits a maximum at a Fermi momentum $`\kappa _{\text{max}}`$ and chemical potential $`\mu _{\text{max}}`$

$$\kappa _{\text{max}}=\frac{\sqrt[3]{3\pi }}{2a_1},\mu _{\text{max}}=\frac{3(3\pi )^{3/2}}{40ma_1^2}.$$
(67)
For values of the chemical potential $`\mu >\mu _{\text{max}}`$ no solution of the extremum condition exists, i.e., there is no metastable low-density state. Equivalently only solutions with local Fermi momenta below $`\kappa _{\text{max}}`$ correspond to minima of the energy functional; those above $`\kappa _{\text{max}}`$ (gray segment of the dashed curve) correspond to maxima of the energy. Thus we get a limiting condition for the local Fermi momentum of the metastable state

$$a_1\kappa (\stackrel{}{x})\frac{\sqrt[3]{3\pi }}{2}$$
(68)
or in terms of the density

$$a_1^3\rho (\stackrel{}{x})\frac{1}{16\pi }.$$
(69)
This is one form of the *stability condition* for the single-component Fermi gas. We note that this condition is completely independent of the trap geometry. As soon as the stability condition is violated somewhere in the trap, in general in the minimum of the trapping potential, the system will become unstable.

For practical purposes we formulate a stability condition in terms of the particle number $`N`$. The maximum particle number $`N_{\text{max}}`$ of the metastable degenerate Fermi gas is directly connected to the maximum chemical potential $`\mu _{\text{max}}`$. This relation is established numerically by solving the extremum condition for the maximum chemical potential $`\mu _{\text{max}}`$ and integrating over the resulting density distribution to obtain the corresponding maximum particle number (56). This is done for several scattering length $`a_1`$ assuming a deformed oscillator potential (49) with mean oscillator length $`\mathrm{}`$. Finally a parametrized form of the stability condition is fitted to this data. The parameterization is motivated by the noninteracting gas, where the maximum local Fermi momentum is proportional to $`\sqrt[6]{N}/\mathrm{}`$. Inserting this into the stability condition (68) leads to the form

$$C(\sqrt[6]{N}\frac{a_1}{\mathrm{}})1,C=2.246.$$
(70)
The parameter $`C`$ is fitted to the numerical results, which are reproduced with a deviation far below 1%. Note that this condition is independent of the deformation of the harmonic oscillator trap Roth (2000).

For an interaction strength of $`a_1/\mathrm{}=0.01`$ which corresponds to a scattering length $`a_1200a_\text{B}`$ for a trap with $`\mathrm{}=1\mu m`$ the maximum particle number is $`N_{\text{max}}=7.8\times 10^9`$. This seems to be out of the range of present experiments. Nevertheless if we increase the strength of the p-wave attraction to $`a_1/\mathrm{}=0.1`$ then the maximum particle number drops to $`N_{\text{max}}=7800`$. Experimentally this could be achieved by utilizing a p-wave Feshbach resonance to increase the p-wave scattering length to $`a_12000a_\text{B}`$ as proposed by J. Bohn Bohn (2000) for the 40K system.

## V Two-Component Fermi Gas

As second application we consider the degenerate two-component Fermi gas.

### V.1 Interplay between s- and p-Wave Interaction

The general energy density of a trapped multi-component Fermi gas in Thomas-Fermi approximation (III.2) takes for the two-component system the following form

$`[\kappa _1,\kappa _2](\stackrel{}{x})`$ $`={\displaystyle \frac{1}{6\pi ^2}}`$ $`\left[U_1(\stackrel{}{x})\kappa _1^3(\stackrel{}{x})+U_2(\stackrel{}{x})\kappa _2^3(\stackrel{}{x})\right]`$
$`+{\displaystyle \frac{1}{20\pi ^2m}}`$ $`\left[\kappa _1^5(\stackrel{}{x})+\kappa _2^5(\stackrel{}{x})\right]`$
$`+{\displaystyle \frac{a_0}{9\pi ^3m}}`$ $`\kappa _1^3(\stackrel{}{x})\kappa _2^3(\stackrel{}{x})`$ (71)
$`+{\displaystyle \frac{a_1^3}{30\pi ^3m}}`$ $`[\kappa _1^8(\stackrel{}{x})+\kappa _2^8(\stackrel{}{x})]`$
$`+{\displaystyle \frac{a_1^3+b_0}{60\pi ^3m}}`$ $`\left[\kappa _1^3(\stackrel{}{x})\kappa _2^5(\stackrel{}{x})+\kappa _1^5(\stackrel{}{x})\kappa _2^3(\stackrel{}{x})\right],`$
where $`\kappa _1(\stackrel{}{x})`$ and $`\kappa _2(\stackrel{}{x})`$ denote the local Fermi momenta of the two components. In contrast to the single-component system, both, s- and p-wave terms of the Effective Contact Interaction contribute. The s-wave interaction acts only between particles of different species and generates a contribution proportional to the product of the densities of both components in the energy density (V.1). The p-wave term acts between particles of different components as well as between particles of the same species. For reasons of simplicity we assume the same p-wave scattering length $`a_1`$ for these different interactions.

Including the constraint of given particle numbers $`N_1`$ and $`N_2`$ of the two components with help of the chemical potentials $`\mu _1`$ and $`\mu _2`$ (see section III.3) leads to the transformed energy density

$$\begin{array}{cc}\hfill [\kappa _1,\kappa _2](\stackrel{}{x})& =[\kappa _1,\kappa _2](\stackrel{}{x})\hfill \\ & \frac{\mu _1}{6\pi ^2}\kappa _1^3(\stackrel{}{x})\frac{\mu _2}{6\pi ^2}\kappa _2^3(\stackrel{}{x}).\hfill \end{array}$$
(72)
Functional variation of the transformed energy functional leads to the extremum condition. For the two-component system the general form (LABEL:eq:ef\_extremumcond) reduces to a coupled set of two polynomial equations

$$\begin{array}{cc}\hfill m[\mu _1U_1(\stackrel{}{x})]& =\frac{1}{2}\kappa _1^2(\stackrel{}{x})+\frac{2a_0}{3\pi }\kappa _2^3(\stackrel{}{x})+\frac{8a_1^3}{15\pi }\kappa _1^5(\stackrel{}{x})\hfill \\ & +\frac{a_1^3+b_0}{30\pi }\left[3\kappa _2^5(\stackrel{}{x})+5\kappa _1^2(\stackrel{}{x})\kappa _2^3(\stackrel{}{x})\right],\hfill \end{array}$$
(73)
where the second equation is generated by the exchange $`\kappa _1(\stackrel{}{x})\kappa _2(\stackrel{}{x})`$ and $`[\mu _1U_1(\stackrel{}{x})][\mu _2U_2(\stackrel{}{x})]`$. Trivial solutions with $`\kappa _1(\stackrel{}{x})=0`$ and $`\kappa _2(\stackrel{}{x})=0`$, resp., are already separated in this expression.

These coupled equations have a great variety of solutions. In order to show the generic phenomena of the two-component system without too many parameters we restrict ourselves to equal numbers of particles in both components $`N=N_1=N_2`$ as well as trapping potentials that differ only by an additive constant, thus $`\mu U(\stackrel{}{x})=\mu _1U_1(\stackrel{}{x})=\mu _2U_2(\stackrel{}{x})`$.

We will concentrate the further studies on the stability of the degenerate two-component Fermi gas against mean-field collapse. For this phenomenon solutions with identical local Fermi momenta for both components, $`\kappa (\stackrel{}{x})=\kappa _1(\stackrel{}{x})=\kappa _2(\stackrel{}{x})`$, are relevant. Under this assumption the extremum condition reduces to a single equation

$$m[\mu U(\stackrel{}{x})]=f_2[\kappa (\stackrel{}{x})]$$
(74)
$$\text{with}f_2(\kappa )=\frac{1}{2}\kappa ^2+\frac{2a_0}{3\pi }\kappa ^3+\frac{4\stackrel{~}{a}_1^3}{5\pi }\kappa ^5.$$

For simplicity we introduce a modified p-wave scattering length 222This can be generalized to include different p-wave scattering lengths for the different combinations of the two species: $`\stackrel{~}{a}_1^3=\frac{1}{3}(a_{1[11]}^3+a_{1[22]}^3+a_{1[12]}^3+b_0)`$.

$$\stackrel{~}{a}_1^3=a_1^3+b_0/3,$$
(75)
which contains the s-wave effective volume parameter. In the following we will discuss the properties of the two-component Fermi gas as function of the s-wave and the modified p-wave scattering length.

For other phenomena, like the separation of the two components due to repulsive interactions, different classes of solutions become important. We will discuss these in a forthcoming publication.

### V.2 Mean-Field Instability: Stability Conditions

The stability of the two-component Fermi gas under the influence of s- and p-wave interactions can be investigated with tools similar to the single-component case. Here the solutions with identical Fermi momenta for both components are of interest.

Similar to the single-component case the right hand side $`f_2(\kappa )`$ of the extremum condition (74) may exhibit a maximum if the s-wave or the p-wave scattering length is negative. Thus the density and particle number of the metastable low-density state may be limited. For a detailed analysis one has to look at all possible combinations of signs of the s- and p-wave scattering lengths separately:

: For a purely repulsive interaction $`f_2(\kappa )`$ is a monotonic growing function and no mean-field induced collapse occurs.

: For purely attractive interactions $`f_2(\kappa )`$ shows a maximum; thus the density of the metastable low-density state is limited.

: The negative contribution of the p-wave interaction dominates $`f_2(\kappa )`$ at high densities and generates a maximum, i.e., the mean-field induced collapse can occur even if the s-wave interaction is repulsive.

: It depends on the relative strength of the s- and p-wave interaction whether the r.h.s. of the extremum condition has a local maximum or grows monotonically.

Especially the stability in the last case depends on a subtle competition between s- and p-wave interaction. Moreover it shows some completely new phenomena, which will be discussed in the following section.

For those cases where $`f_2(\kappa )`$ has a maximum the value of the local Fermi momentum $`\kappa _{\text{max}}`$ at the maximum is given by the equation

$$a_0\kappa _{\text{max}}2[\stackrel{~}{a}_1\kappa _{\text{max}}]^3=\frac{\pi }{2}.$$
(76)
Again $`\kappa _{\text{max}}`$ is an upper limit for the local Fermi momenta, which can occur for a metastable low-density state of the two-component gas. Thus we can formulate the stability condition

$$a_0\kappa (\stackrel{}{x})2[\stackrel{~}{a}_1\kappa (\stackrel{}{x})]^3\frac{\pi }{2}$$
(77)
or equivalently in terms of the density

$$[6\pi ^2a_0^3\rho (\stackrel{}{x})]^{1/3}12\pi ^2\stackrel{~}{a}_1^3\rho (\stackrel{}{x})\frac{\pi }{2}.$$
(78)
If these stability conditions are violated than no metastable low-density state exists for the two-component Fermi gas. For a pure s-wave interaction ($`\stackrel{~}{a}_1=0`$) the stability condition (77) reduces to the form $`a_0\kappa (\stackrel{}{x})\pi /2`$, which was obtained earlier by M. Houbiers *et. al.* Houbiers et al. (1997). Compared to this simple form the inclusion of the p-wave interaction reveals several new effects.

Before we discuss the structure of (77) we formulate an equivalent stability condition in terms of the number of particles $`N=N_1=N_2`$ of each component. For given values of the two scattering lengths the maximum local Fermi momentum and the maximum chemical potential is calculated. From the solution of the extremum condition (74) for these parameters the corresponding maximum particle number is determined. This numerical data is fitted by a suitable parameterization of the stability condition in terms of the particle number and the two scattering lengths $`a_0/\mathrm{}`$ and $`\stackrel{~}{a}_1/\mathrm{}`$

$$\begin{array}{cc}& C_0(\sqrt[6]{N}\frac{a_0}{\mathrm{}})+C_1^3(\sqrt[6]{N}\frac{\stackrel{~}{a}_1}{\mathrm{}})^3\hfill \\ & +C_{01}^{n+1}(\sqrt[6]{N}\frac{a_0}{\mathrm{}})(\sqrt[6]{N}\frac{\stackrel{~}{a}_1}{\mathrm{}})^n1.\hfill \end{array}$$
(79)
This parameterization is constructed in analogy to the single-component case (70); the additional cross-term is necessary to achieve a similar accuracy with typical deviations below 1%. The parameters $`C_0`$, $`C_1`$, and $`C_{01}`$ have to be fitted for each combination of signs of the two scattering length separately. The value of $`n`$ is not included in the fitting procedure but chosen by hand. The resulting values are summarized in Table 1.

Figure 8 illustrates the dependence of the maximum particle number resulting from this stability condition on the s- and p-wave scattering lengths. The contour plot shows the logarithm of the maximum particle number for each component as function of $`a_0/\mathrm{}`$ and $`\stackrel{~}{a}_1/\mathrm{}`$. The first gross observation is that attractive s- and p-wave interactions with similar scattering lengths set similar restrictions to the stability of the two-component Fermi gas. For example, a pure s-wave interaction with $`a_0/\mathrm{}=0.05`$ leads to a maximum particle number of $`N_{\text{max}}1.7\times 10^6`$. In comparison a pure p-wave interaction with same scattering length $`\stackrel{~}{a}_1/\mathrm{}=0.05`$ causes a collapse at even lower particle numbers of $`N_{\text{max}}2.2\times 10^5`$.

If both interaction parts are attractive they cooperate and cause an instability at lower particle numbers or densities.

If the interaction is attractive in one partial wave and repulsive in the other than the repulsive part leads to a stabilization, i.e., it increases the maximum particle number. Here a significant difference between s- and p-wave interactions arises: The stabilization caused by a repulsive s-wave interaction is rather weak. Compared to a pure p-wave interaction with $`\stackrel{~}{a}_1/\mathrm{}=0.05`$ the presence of a s-wave repulsion of the same magnitude $`a_0/\mathrm{}=0.05`$ increases the maximum particle number only from $`2.2\times 10^5`$ to $`8.9\times 10^5`$. In the opposite case of an attractive s-wave interaction a p-wave repulsion of same magnitude will always lead to an absolute stabilization, i.e., there is no collapse for arbitrary large particle number despite of the s-wave attraction. We will study these special effects in detail in the following section.

These results clearly demonstrate that it is necessary to include the p-wave interaction if the scattering length $`\stackrel{~}{a}_1`$ is roughly in the same order of magnitude as the s-wave scattering length. Even if the ratio of the scattering lengths, $`\stackrel{~}{a}_1/a_0`$, is approximately 0.3 dramatic effects like the p-wave stabilization, which is discussed in the next section, can occur. As can be seen from Figure 8 the p-wave interaction may be neglected only if the ratio $`\stackrel{~}{a}_1/a_0`$ is smaller than 0.1.

### V.3 Mean-Field Instability: p-Wave Stabilization

Several new phenomena occur due to the competition between an attractive s-wave ($`a_0<0`$) and a repulsive p-wave interaction ($`\stackrel{~}{a}_1>0`$). To understand the origin of these phenomena, which are a unique property of this type of interactions, we investigate the right hand side $`f_2(\kappa )`$ of the extremum condition (74). Figure 9 depicts the dependence of $`f_2(\kappa )`$ on the local Fermi momentum for a s-wave scattering length $`a_0/\mathrm{}=0.05`$ and three slightly different positive p-wave scattering lengths in the range $`a_0/\mathrm{}=0.014\mathrm{}0.016`$.

Due to the dominant $`\kappa ^5`$-dependence any repulsive p-wave interaction causes $`f_2(\kappa )`$ to grow fast for large Fermi momenta. Thus the maximum is only local and does not determine necessarily the maximum Fermi momentum or chemical potential as in the cases with attractive or vanishing p-wave interaction. If the p-wave repulsion is sufficiently strong the local maximum vanishes completely and $`f_2(\kappa )`$ is a monotonically growing function. Then a solution of the extremum condition exists for any density, chemical potential or particle number. An example is shown by the dashed curve in Figure 9. It can be seen from equation (74) that the local maximum disappears if the ratio of the two scattering lengths fulfills the condition

$$\frac{\stackrel{~}{a}_1}{|a_0|}\frac{2}{3\pi ^{2/3}}0.311.$$
(80)
If this condition is fulfilled the p-wave repulsion causes an *absolute stabilization* of the system against s-wave induced collapse. In this case the total mean-field contribution of the interactions is always repulsive and grows monotonically with density.

Notice that the p-wave scattering length necessary for this stabilization is only approx. $`1/3`$ of the modulus of the s-wave scattering length. Obviously the p-wave interaction may have drastic influence on the stability even if it is significantly weaker than the s-wave interaction in terms of scattering lengths.

For weaker p-wave repulsions $`f_2(\kappa )`$ still shows a local maximum and in addition a (local) minimum at larger Fermi momenta. Examples are shown by the solid and dotted curves in Figure 9. In this case the extremum condition has two separate branches that correspond to local minima of the energy density, which are separated by branch of local maxima (gray segments). The branch at lower Fermi momenta corresponds to the usual family of low-density solutions that were obtained with other types of interactions too. It ends up at the local maximum with $`\kappa _{\text{max}}`$ given by (76) and $`m\mu _{\text{max}}=f_2(\kappa _{\text{max}})`$. The solution branch at higher Fermi momenta gives rise to a new family of high-density solutions, which are unique for this type of interaction. It is bounded from below by the local minimum $`(\kappa _{\text{min}},\mu _{\text{min}})`$ and raises up to arbitrary Fermi momenta and chemical potentials.

For values of $`[\mu U(\stackrel{}{x})]`$ between $`\mu _{\text{min}}`$ and $`\mu _{\text{max}}`$ the extremum condition has two solutions $`\kappa _{\text{low}}`$ and $`\kappa _{\text{high}}`$ with $`f_2(\kappa _{\text{low}})=f_2(\kappa _{\text{high}})`$, see Figure 10. In equilibrium the one with lower free energy density (72) is realized. We define a chemical potential $`\mu _{\text{trans}}`$ at which the free energy densities of both branches are equal; the value of $`\mu _{\text{trans}}`$ can be determined numerically. Since we expect the solution $`\kappa (\stackrel{}{x})`$ to correspond to a minimum of the energy functional at each point $`\stackrel{}{x}`$, for $`[\mu U(\stackrel{}{x})]<\mu _{\text{trans}}`$ the low-density branch gives the equilibrium solution and for $`[\mu U(\stackrel{}{x})]>\mu _{\text{trans}}`$ the high-density branch does. This gives rise to a Maxwell construction for the r.h.s. of the extremum condition (74) as illustrated in the upper plot of Figure 10. The dotted parts of the lower and upper branch in Figure 10 correspond to local minima and may occur as metastable states that eventually undergo a transition to the energetically lower equilibrium solution.

The structure of the density distribution depends crucially on the value of the chemical potential $`\mu `$, i.e., the particle number. The upper plot of Figure 10 shows the r.h.s. of the extremum condition (74) for an interaction with $`a_0/\mathrm{}=0.05`$ and $`\stackrel{~}{a}_1/\mathrm{}=0.015`$. The dashed horizontal lines mark three different chemical potentials. The lower plot shows the radial dependencies of the local Fermi momenta $`\kappa (x)=(6\pi ^2\rho (x))^{1/3}`$ for these three chemical potentials assuming a spherical trap with oscillator length $`\mathrm{}`$. For chemical potentials $`\mu <\mu _{\text{trans}}`$ — case (a) — the equilibrium solution is completely on the low-density branch and we obtain the usual smooth density profile. If $`\mu >\mu _{\text{trans}}`$ — case (b) and (c) — then the equilibrium solution in the center of the trap is given by the high-density branch, while the solution for the outer regions of the trap is given by the low-density branch. Thus the equilibrium density profile shows a jump in density by typically one order of magnitude as one approaches the center of the trap. The location of the discontinuity is always given by the equation $`[\mu U(\stackrel{}{x})]=\mu _{\text{trans}}`$ which reflects the mechanical equilibrium between the low and the high-density phase, i.e., the equality of the pressures $`p=[\kappa ](\stackrel{}{x})`$ at the boundary between the two phases.

If the chemical potential is still below $`\mu _{\text{max}}`$ — case (b) — then a solution with a smooth low-density profile all over the trap may exist as metastable state. Due to density fluctuations this state may undergo a transition to the energetically preferred equilibrium state, which includes the high-density phase.

The physical origin of the high-density phase is quite intuitive: The usual stability is determined by a competition between kinetic energy, which favors low densities, and attractive mean-field, which prefers higher densities. If the attractive mean-field becomes too strong, than the kinetic energy is not able to stabilize the system anymore and the mean-field collapse occurs. In case of the high-density phase the attractive mean-field generated by the s-wave interaction has already overcome the stabilizing effect of the kinetic energy. Nevertheless the collapse is prevented by the repulsive p-wave contribution which grows for higher densities faster than the s-wave attraction and inhibits a further increase of density. We call this new phenomenon *p-wave stabilized high-density phase* — in contrast to the low-density phase stabilized by the kinetic energy, which is still present in the peripheral regions of the trap.

The situation discussed so far assumes a repulsive p-wave interaction that is slightly to weak to cause the absolute stabilization according to (80). If the p-wave strength is decreased further, then the values of $`\mu _{\text{min}}`$, $`\mu _{\text{max}}`$, and $`\mu _{\text{trans}}`$ also decrease. If the ratio of the p-wave and s-wave scattering lengths drops below the limit

$$\frac{\stackrel{~}{a}_1}{|a_0|}<\sqrt[3]{\frac{160}{729\pi ^2}}0.281$$
(81)
then the chemical potential of the minimum is negative, $`\mu _{\text{min}}<0`$. An example is shown by the dotted curve in Figure 9. For even weaker p-wave interactions with

$$\frac{\stackrel{~}{a}_1}{|a_0|}<0.274$$
(82)
$`\mu `$ can be negative. That means that the high-density solution forms a self-bound state independent of the trapping potential. Therefore as soon as the maximum chemical potential $`\mu _{\text{max}}`$ is exceeded the gas collapses into a self-bound high-density state which is independent of the trap.

We summarize the variety of structures that appear for interactions with attractive s-wave ($`a_0<0`$) and repulsive p-wave part ($`\stackrel{~}{a}_1>0`$) in the following list

the p-wave repulsion stabilizes the system for arbitrary densities and particle numbers with a smooth low-density profile. For ratios of the scattering lengths near the limit a smooth but significant increase of the central density occurs.

For chemical potentials below $`\mu _{\text{max}}`$ or particle numbers below the corresponding maximum particle number (LABEL:eq:2comp\_stabilitycond\_N) the usual low density solution exists. Above these values the p-wave stabilized high-density phase appears, i.e., the density in the central region of the trap is increased by typically one order of magnitude compared to the low-density profile in the outer regions.

Below $`\mu _{\text{max}}`$ or $`N_{\text{max}}`$ a stable solution with the regular low-density profile exists. Above the system collapses; the high-density phase is not stable anymore.

This subtle dependence on the ratio of the scattering lengths is illustrated in Figure 8. The white region for “strong” repulsive p-wave interactions shows the domain of absolute stabilization. The solid line corresponds to the condition (80) for absolute stabilization, the dashed line to condition (82) for the stability of the high-density phase. The small area between those lines represents the parameter region where the p-wave stabilized high-density phase exists if the maximum particle number is exceeded.

To conclude this section we perform a *Gedankenexperiment*: Assume a two-component Fermi gas of $`N=N_1=N_2=60000`$ particles in each component trapped in a spherical oscillator potential with $`\mathrm{}=1\mu m`$. May the interaction be composed of a repulsive p-wave part with $`\stackrel{~}{a}_1/\mathrm{}=0.03`$ and an attractive s-wave component which can be tuned within a small range $`a_0/\mathrm{}=0.095,\mathrm{},0.102`$, e.g., by using a Feshbach resonance. Figure 11 shows the evolution of the density profile of the Fermi gas if the strength of the attractive s-wave interaction is increased slowly such that density fluctuations are negligible. For $`a_0/\mathrm{}=0.095`$ (solid curve) and $`0.1`$ (dash-dotted) we observe a smooth low density profile, where the central density increases slightly with increasing s-wave attraction. A dramatic change happens if the attraction is increased to $`a_0/\mathrm{}=0.101`$ (dashed). For this interaction strength the particle number of the system is already above the maximum particle number given by (LABEL:eq:2comp\_stabilitycond\_N) and the high-density phase appears and occupies a rather large volume. Further increase of the s-wave attraction (dotted) causes a growth of the high-density phase. If the limit $`a_0/\mathrm{}=0.11`$ is reached than part of the high-density component is self-bound and the system is expected to collapse.

## VI Summary and Conclusions

We formulated a simple and transparent model to describe the structure and stability of degenerate multi-component Fermi gases trapped in an external potential. In a first step we derived an Effective Contact Interaction (ECI) for all partial waves that reproduces the exact two-body energy spectrum when used in a mean-field model space. Including the s- and p-wave parts of the ECI we constructed the energy density of the inhomogenous Fermi gas in a mean-field calculation using the Thomas-Fermi approximation. By functional minimization of the energy we obtained a set of coupled polynomial equations for the ground state density profile of the system. We showed that the combination of s- and p-wave interactions leads to a rich variety of phenomena in trapped degenerate Fermi gases.

In the single-component system the p-wave part is the leading interaction term since s-wave scatterings are prohibited by the Pauli principle. Attractive p-wave interactions cause a mean-field instability of the one-component gas if a certain maximum density is exceeded. We derived explicit stability conditions in terms of the density or particle number and the p-wave scattering length.

The interplay between s- and p-wave interactions leads to several new effects in the two-component Fermi gas. We discussed the dependence of the mean-field instability on the s- and p-wave scattering lengths and derived also for this case closed stability conditions. It turns out that attractive s-wave as well as attractive p-wave interactions can cause a mean-field collapse. In addition a repulsive interaction part leads to stabilization, i.e., an increase of the maximum possible density of the Fermi gas.

Interactions with attractive s-wave and repulsive p-wave parts show several special properties. If the p-wave scattering length exceeds about $`1/3`$ of the modulus of the s-wave scattering length the s-wave attraction is fully compensated and no mean-field collapse occurs anymore at high densities. In the transition region towards this absolute stabilization a distinct high-density phase may appear in the center of the trap, which is stabilized by the p-wave repulsion alone.

We conclude that the p-wave interaction may have an important influence on the structure and stability of dilute degenerate Fermi gases. Considering the simultaneous s- and p-wave Feshbach resonances predicted for the 40K system Bohn (2000) it can be foreseen that large values of p-wave scattering length will be available experimentally. In a two-component 40K gas this Feshbach resonance would allow to probe nearly the whole stability map shown in Figure 8 by modifying the magnetic field. Alternatively, tightly confining optical traps Viverit et al. (2001) generate large values of the ratio of scattering length and oscillator length, $`a_l/\mathrm{}`$, such that instabilities occur at much lower densities and particle numbers ($`N_{\text{max}}\mathrm{}^6`$; see equations (70) and (LABEL:eq:2comp\_stabilitycond\_N)).

Concerning the envisioned observation of Cooper pairing in trapped dilute Fermi gases two-component systems with strong attractive s-wave interactions are favored DeMarco et al. (2001); O’Hara et al. (2000). Here the mean-field instability limits the density of the normal Fermi gas. With a suitably chosen repulsive p-wave interaction one could use the effect of absolute stabilization, which we discussed, to allow higher densities of the normal Fermi gas and thus increase the transition temperature to a superfluid state Houbiers et al. (1997).