PHYSICAL REVIEW A 80, 013607 共2009兲

Repulsive Fermi gas in a harmonic trap: Ferromagnetism and spin textures
1

L. J. LeBlanc,1 J. H. Thywissen,1 A. A. Burkov,2 and A. Paramekanti1

Department of Physics, University of Toronto, Toronto, Ontario, Canada M5S 1A7
Department of Physics, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1
共Received 31 March 2009; published 13 July 2009兲

2

We study ferromagnetism in a repulsively interacting two-component Fermi gas in a harmonic trap. Within
a local density approximation, the two components phase separate beyond a critical interaction strength, with
one species having a higher density at the trap center. We discuss several easily observable experimental
signatures of this transition. The mean-field release energy, its separate kinetic and interaction contributions, aswell as the potential energy all depend on the interaction strength and contain a sharp signature of thistransition. In addition, the conversion rate of atoms to molecules, arising from three-body collisions, peaks at
an interaction strength just beyond the ferromagnetic transition point. We then go beyond the local density
approximation and derive an energy functional that includes a term that depends on the local magnetization
gradient and acts as a “surface tension.” Using this energy functional, we numerically study the energetics ofsome candidate spin textures that may be stabilized in a harmonic trapping potential at zero net magnetization.
We find that a hedgehog state has a lower energy than an “in-out” domain-wall state in an isotropic trap. Upon
inclusion of trap anisotropy we find that the hedgehog magnetization profile gets distorted due to the surface
tension term, this distortion being more apparent for small atom numbers. We estimate that the magnetic dipole
interaction does not play a significant role in this system. We consider possible implications for experiments on
trapped 6Li and 40K gases.
DOI: 10.1103/PhysRevA.80.013607

PACS number共s兲: 03.75.Ss

I. INTRODUCTION

In recent years, a series of beautiful experiments has
shown that a gas of two-component fermions interacting via
a Feshbach resonance exhibits a superfluid state at low temperature 关1–4兴. An exciting new direction for cold atom experiments is the study of ferromagnetism arising from repulsive interactions in a two-component Fermi gas. Such aferromagnetic “Stoner instability” 关5兴 occurs, within meanfield theory of a homogeneous Fermi gas at zero temperature, when the 共repulsive兲 s-wave scattering length, aS, between two spin states is large enough that kFaS ⬎ ␲ / 2, whereបkF is the Fermi momentum of the gas. This condition can besatisfied upon tuning aS to large positive values near a Feshbach resonance provided the system stays stable for a sufficiently long time. Since ferromagnetism arises from twobody interactions, whereas atom loss due to Feshbach
molecule formation is because of three-body collisions 关6兴,
there may be a range of densities where the lifetime is longenough to reach the ferromagnetic state.
The suggestion that ferromagnetism may be achieved as ametastable state in cold Fermi gases is not new. Salasnichand co-workers 关7兴 studied the mean-field theory of a harmonically trapped Fermi gas with repulsive interactions andfound that this should lead to phase separation between thetwo spin species if the net magnetization is zero. A similar
study was carried out by Sogo and Yabu 关8兴 allowing fornonzero spontaneous magnetization. Duine and MacDonald
关9兴 later showed that the ferromagnetic transition in a homogeneous Fermi gas changes from a continuous to a first-ordertransition at low enough temperatures upon going beyond
mean-field theory. They also proposed that an initially magnetized Fermi gas will tend to stay spin-coherent for longtimes even in the presence of magnetic field noise that is1050-2947/2009/80共1兲/013607共10兲

naively expected to cause strong dephasing, provided thesystem is close to the transition into a ferromagnetic state关9兴.
Using an optical lattice and engineering the band structure
to get flat 共dispersionless兲 bands are other interesting routes
to achieving ferromagnetism. Such “flat-band ferromagnetism” 关10兴 has the advantage that the ferromagnetic state
occurs at weak repulsive interactions and can be theoreticallyanalyzed in a reliable fashion; however, an existing theoretical proposal along these lines involves working with fermions in the p band of a honeycomb optical lattice 关11兴. A
recent work 关12兴 has considered the possibility of ferromagnetism for strongly interacting fermions in optical lattices
and studied, within a phenomenological Landau theory, theenergetics of possible spin textures 共such as hedgehog states,
domain walls, and skyrmions兲 that might arise in a trap.
In this paper, motivated by earlier work and by ongoingexperimental efforts, we revisit the problem of ferromagnetism in a harmonically trapped two-component Fermi gas.
We begin by using a “local density approximation” 共LDA兲,
sometimes referred to in the literature as the Thomas-Fermi
approximation, to describe this system. Within the LDA, wefind that the mean-field release energy of the trapped gas 共aswell as the potential energy and the kinetic energy components of the release energy兲 provides a simple, albeit indirect,
diagnostic of the ferromagnetic transition. We find that theformation of nonzero local magnetization in the trap causes a
suppression of the atom loss rate via three-body collisions.
This suppression competes with the growth of the loss rate as
the interactions get stronger, leading to a peak in the atom
loss rate at an interaction strength that is very close to, butslightly beyond, the ferromagnetic transition point. We thenshow how one might incorporate magnetization gradient 共or“surface tension”兲 terms in order to go beyond the LDA. Our
energy functional is akin to an earlier phenomenological013607-1©2009 The American Physical Society

PHYSICAL REVIEW A 80, 013607 共2009兲

LEBLANC et al.

Landau theory for ferromagnetism in an optical lattice 关12兴but explicitly keeps track of the spatial dependence of all the
Landau theory coefficients that arise from density variationsin the trap. Using this extended energy functional, we study
the energetics of various spin textures in the case where thenet magnetization, which does not relax in these quantum
gases, is assumed to be zero. This corresponds to choosing
the initial population to be the same for both hyperfine species of fermions. We show, in this case, that a hedgehog
configuration of the magnetization has a lower energy thanan “in-out” phase-separated configuration with a domain
wall. A similar phenomenon has been predicted for fermionsin optical lattices 关12兴. Finally, we turn to the effect of anisotropic trapping frequencies in a harmonic trap. While suchanisotropies can be incorporated by a trivial rescaling of coordinates in the LDA, this is no longer true in the presence of
surface tension that leads to a breakdown of the LDA. 共A
breakdown of the LDA has been observed 关13兴 and theoretically addressed 关14–16兴 in the context of polarized superfluids in highly anisotropic traps.兲 We use our extended energyfunctional to numerically study how the hedgehog state distorts upon going from a spherically symmetric trap to ananisotropic cigar-shaped trap. We conclude with estimatesthat indicate that the magnetic dipole interaction between
atoms can be neglected for 6Li and 40K.
II. FERROMAGNETISM WITHIN THE LOCAL DENSITYAPPROXIMATIONThe Hamiltonian describing a uniform two-component
Fermi gas interacting through a repulsive s-wave contact interaction is given by
H=兺
␴

冕

d 3K
†
⑀ Kc K␴c K␴ + g
共2␲兲3

冕

共2兲

␴

册

共4兲

A. Rescaling to the isotropic problem for a harmonic trap

Let us assume that the trapping potential is harmonic but
possibly anisotropic, so that V共R兲 = 21 M兺i␻2i R2i . Here ␻x,y,z
are the trapping frequencies along different directions. In order to make the problem appear isotropic we can rescaledistances by setting R̃i = Ri共␻i / ⍀兲, where ⍀ = 共␻x␻y␻z兲1/3 isthe geometric mean of the trap frequencies. With this standard rescaling, we get
E兵关␳␴共R̃兲兴其 =

冕 冋 兺
d3R̃

3
␣ ␳5/3共R̃兲 + g␳↑共R̃兲␳↓共R̃兲
5 ␴ ␴

册

1
+ M⍀2R̃2 兺 ␳␴共R̃兲 − 兺 ␮␴␳␴共R̃兲 . 共5兲
2␴␴

B. Noninteracting unmagnetized gas

For the unmagnetized gas, we have N↑ = N↓ = N / 2 and forthe noninteracting case we can set g = 0. This reduces theenergy functional to
EN0 兵关␳␴共R̃兲兴其 = 兺 EN0 ␴关␳␴共R̃兲兴,

共6兲␴

EN0 ␴关␳␴共R̃兲兴 =共3兲

At this level of treatment the contact interaction strength g isrelated to the two-body scattering length aS in vacuum via
g = 4␲aSប2 / M.
The LDA for a trapped Fermi gas corresponds to simply
assuming that the above results apply locally in the presenceof a trap potential V共R兲. The ground-state energy of thistrapped Fermi gas is then obtained by minimizing the energyfunctional,3␣ ␳5/3共R兲 + g␳↑共R兲␳↓共R兲
5 ␴ ␴

where 兵␳␴共R兲其 denotes the density profile of both spin species 关␳↑共R兲 , ␳↓共R兲兴. Here we have introduced two Lagrange
multipliers ␮␴ that act as chemical potentials for the two spin
species and serve to impose the constraints 兰d3R␳␴共R兲 = N␴.
The separate constraint on each spin component arises fromthe assumption that the two spin components correspond to
the lowest two Zeeman split hyperfine levels of a Fermi gas.
Since the Zeeman splitting near a Feshbach resonance istypically far greater than all other energy scales and the totalenergy must be conserved in these thermally isolated gases,
we arrive at the constraint that the population of the two
Zeeman components cannot change for fermionic atoms
where the only interaction is between different spin components. Thus, unlike in solid state ferromagnets, the magnetization can be conserved on very long time scales.

where V denotes the system volume and EF␴ = ␣␳␴2/3, with
␣ = 共6␲2兲2/3ប2 / 2M, denotes the Fermi energy of particles
with spin ␴. A mean-field theory of the interacting Hamiltonian yields the total interaction energy,
U = gV␳↑␳↓ .

d 3R

␴

2

K = 53 V共EF↑␳↑ + EF↓␳↓兲,冕 冋 兺+ V共R兲 兺 ␳␴共R兲 − 兺 ␮␴␳␴共R兲 ,† †d3RcR↑
cR↓cR↓cR↑ , 共1兲

where ⑀K = ប K / 2M is the kinetic energy of atoms withmass M and momentum បK. For a Fermi gas with N␴ particles of spin ␴, the uniform gas densities of each spin is␳␴ = N␴ / V and the total kinetic energy of the uniform gas isjust
2E兵关␳␴共R兲兴其 =

冕 冋
d3R̃

3 5/3
␣␳ 共R̃兲5 ␴

册

1+ M⍀2R̃2␳␴共R̃兲 − ␮␴␳␴共R̃兲 .
2

共7兲

Setting ␦EN0 / ␦␳␴ = 0 leads to the equation

␣␳␴2/3共R̃兲 = 共␮N0 − 21 M⍀2R̃2兲,

共8兲

where we have used symmetry to set ␮↑ = ␮↓ = ␮N0 . The solu-

tion to this equation is simply

013607-2

PHYSICAL REVIEW A 80, 013607 共2009兲

REPULSIVE FERMI GAS IN A HARMONIC TRAP:…

␳␴共R̃兲 = ␣−3/2关␮N0 − 21 M⍀2R̃2兴3/2 .

共9兲

␧兵关n␴共r兲兴其 =Clearly there is a maximum radius,
RN0 =

冑

+

2␮N0
2,

共10兲

M⍀

beyond which ␳␴共R̃兲 = 0. Integrating up to this maximum radius using spherical symmetry of the density and employing
the constraint, we end up with

␮N0 ⬅ EF0 = ប⍀共3N兲1/3 ,
RN0 =

冑

2EF0
= aHO共24N兲1/6 ,m a⍀ 23
ប⍀共3N兲4/3 ,EN0 = NEF0 =
4
4

冉 冊

2 1/2 1/2
0
−3
␳N,
共0兲
=
=
a
N ,
HO
␴
3␲4
␲2共RN0 兲3
4N␧兵关n␴共r兲兴其 =

n␴ =

␳␴
␳N0 ␴共R̃ = 0兲

,

共16兲

␭ = kF0 共0兲aS ,共17兲␧ = E/EN0 ,共18兲␥␴ = ␮␴/␮N0 .

共19兲

Here ␭ is the dimensionless interaction parameter, ␧ , ␥␴ are
the total energy and chemical potential, respectively, in dimensionless units, and kF0 共0兲 denotes the Fermi wave vector
at the trap center for the unmagnetized noninteracting gas.
In terms of these dimensionless variables, the energyfunctional becomes

冕

drr2

冋

3 5/3
关n 共r兲 + n5/3
↓ 共r兲兴
5 ↑册

4␭n↑共r兲n↓共r兲 − 兺 共␥␴ − r2兲n␴共r兲 . 共21兲3␲␴The variational minimization ␦E / ␦n␴共r兲 leads to the following two equations:n↑共r兲 =n↓共r兲 =冋冉冋冉␥↑ − r2 −4␭n↓共r兲
3␲␥↓ − r2 −4␭n↑共r兲
3␲

冊册
冊册

3/2

,

共22兲

,

共23兲

3/2

subject to the constraints
4␲

共15兲

册

D. Variational minimization

C. Converting the interacting problem to dimensionlessvariables

R̃
,
RN0

3 5/3
关n 共r兲 + n5/3
↓ 共r兲兴
5 ↑4␭n↑共r兲n↓共r兲 − 兺 共␥␴ − r2兲n␴共r兲 .
3␲
␴

64
3␲
+

共13兲

Let us use the noninteracting unmagnetized gas results toconvert to dimensionless variables as follows:d 3r

If we assume that the ground-state solution for the densities
respects the spherical symmetry of this energy functional, we
can further simplify the energy functional to a onedimensional integral,

共12兲

共14兲

冕 冋

共20兲

共11兲

where aHO = 共ប / M⍀兲1/2 is the oscillator length, and we denote the density solution for this noninteracting unmagnetized Fermi gas by ␳N0 ␴共R̃兲. Here ␮N0 is the chemical potentialof the gas, EN0 is the total energy of the gas, and ␳N0 共0兲 is theatom density at the trap center.

r=

16
3␲2

冕

drr2n␴共r兲 =

␲2 N␴
.
4 N

共24兲

These equations can be iteratively solved 共numerically兲 forthe fermion densities given the interaction strength and thetotal fermion numbers for each species. Having solved them
we can use the resulting fermion densities to compute physical observables. We will denote the average magnetizationby m̄ ⬅ 共N↑ − N↓兲 / 共N↑ + N↓兲. For an unmagnetized gas, wefind that increasing the interaction ␭ progressively modifies
the density profile of the gas from that of a noninteracting
Thomas-Fermi profile at ␭ = 0 to that of a fully polarized gas
when ␭ Ⰷ 1.
E. Release energy and “ferromagnetic transition”

The release energy of the trapped atomic gas is measured
by rapidly switching off the trap potential and measuring thetotal kinetic energy of the atoms after some time delay. Assuming that the switch-off process is instantaneous and that
all the interaction energy in the initial state has been converted into the kinetic energy of atoms at the time of measurement, the release energy and its separate kinetic and interaction energy components are given, within the LDA, by

013607-3

␧rel =

64
3␲

冕

drr2

冋兺

册

4
3
n␴5/3共r兲 +␭n↑共r兲n↓共r兲 , 共25兲5 ␴
3␲PHYSICAL REVIEW A 80, 013607 共2009兲

LEBLANC et al.

冕 冋
冕 冋兺
drr2

64
␧kin =
3␲

drr

册

4
␭n↑共r兲n↓共r兲 ,
3␲

2

册

3
n5/3共r兲 .
5 ␴ ␴

0.64

共26兲

0.62
0.60

共27兲

64
3␲

冕

drr2兵r2关n↑共r兲 + n↓共r兲兴其.

0.6

0.56

0.4

0.54

m=0.0

0.52
0.50
0.0

共28兲

0.65

From Fig. 1, we see that the release energy displays asharp transition point, for m̄ = 0, at ␭crit ⬇ 1.84. At this interaction strength, we find that kFaS = ␲ / 2 at the trap center,with kF being the Fermi wave vector at the trap center in theinteracting cloud, which corresponds to the Stoner transition
point in the uniform gas. Further, an examination of the density profile of the two spin species shows that this corresponds to an onset of phase separation in the trap—for ␭⬎ ␭crit, atoms of one spin type tend to have a higher density
at the trap center while atoms of the other spin type arepushed away from the center leading to a nonzero magnetization density near the trap center. Exactly which atom tends
to accumulate at the center is a spontaneously broken symmetry at zero magnetization, and this phase separation is
simply a local manifestation of ferromagnetic ordering. This
result for the ␭crit at m̄ = 0 translates into an estimate for thecritical two-body scattering length,⬇ 0.6␭critaHON−1/6 ⬇ aHON−1/6 ,a共crit兲
S

0.8

0.58

The potential energy of the cloud, due to the confining harmonic trap, can be easily obtained from experimental measurements of the cloud profile and it is given by␧pot =(a) ε rel

0.60

0.5

1.0

共31兲

where a0 ⬇ 0.529 Å is the Bohr radius.
Figure 1 also shows the kinetic energy and the interactionenergy components of the total release energy. Each of theseobservables shows a large and much more dramatic signature
at the transition 共for m̄ = 0兲 than the total release energy and it
is likely to be more promising experimental signature, as
discussed in Ref. 关17兴. In addition, the potential energy of theconfined cloud also shows a maximum at the ferromagnetictransition point.
Strictly speaking, there is no phase transition 共beyond
mean-field theory兲 except in the thermodynamic limit that,for a trapped Fermi gas, is obtained by taking N → ⬁ and⍀ → 0 with N⍀3 held fixed. For nonzero magnetization,
however, there is no phase transition even at mean-fieldlevel; nevertheless, the release energy does display a fairly
sharp crossover at ␭crit for m̄ ⱗ 0.2. The measured releaseenergy can only tell us about the existence of a phasetransition—for ␭ ⬎ ␭crit, in situ measurements of the magnetization profile, which we discuss below, are needed to showthat this transition corresponds to ferromagnetism in the trap.

3.0

3.5

4.0

0.50

2.5

3.0

3.5

4.0

0.6

0.45

0.4
0.2
|

0.40

m=0.0
0.35
0.0
0.20

0.5

1.0

1.5

(c) ε int

2.0

m=0.0
0.2

0.15

0.4
0.10

0.6
0.05

0.8
0.00
0.0
0.66

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

1.5

2.0

2.5

3.0

3.5

4.0

(d) εpot

0.62

0.8
0.58

0.6

0.54

|

共 6Li兲 ⬇ 13800a0,9500a0,6400a0 ,a共crit兲
S

2.5

0.8

0.55

beyond which one expects to see phase separation in the trap.
For ⍀ / 2␲ ⬇ 170 Hz we estimate for N = 104 , 105 , 106 the respective critical scattering lengths,
共30兲

2.0

(b) εkin

共29兲

共 40K兲 ⬇ 5300a0,3600a0,2500a0 ,a共crit兲
S

1.5

|

64
3␲

|

␧int =

m=0.0

0.50
0.0

0.5

1.0

λ
FIG. 1. 共a兲 Dimensionless mean-field release energy, ␧rel, versus
interaction, ␭, within the LDA for indicated values of m̄ = 共N↑− N↓兲 / N. For m̄ = 0, there is a phase separation transition that appears as a sharp kink in the release energy at ␭ ⬇ 1.84. 共b兲–共d兲 Same
as in 共a兲 but for the kinetic energy, interaction energy, and potential
energy of the gas. The energy per particle in physical units may beobtained by multiplying these results by 3E0F / 4, where E0F= ប⍀共3N兲1/3. As shown later, going beyond the LDA leads to negligible quantitative corrections to these results.

013607-4

PHYSICAL REVIEW A 80, 013607 共2009兲|

REPULSIVE FERMI GAS IN A HARMONIC TRAP:…

m=0.0

1.25
1.00

0.4

Γ/Γ

0.75

0

0.6
0.50

0.8

0.25
0.00
0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

λ
FIG. 2. Dimensionless atom loss rate, ⌫ / ⌫0, as a function ofinteraction strength at various average magnetizations.
F. Atom loss rate

Atoms on the repulsive side of the Feshbach resonancetend to be unstable to formation of molecules via three-bodycollisions. Apart from kinematic and statistical constraints onthese processes, there is a simple constraint that two of theseatoms, which eventually form the molecule, must have opposite spins. One consequence of having a nonzero local
magnetization in the trapped gas is a suppression of the probability of finding fermions with opposite spin in the same
region of space, which leads to a strong suppression of such
three-body losses. Up to an unknown prefactor, ⌫0, we canestimate this three-body loss rate as⌫ = ⌫ 0␭ 6冕

d3rn↑共r兲n↓共r兲关n↑共r兲 + n↓共r兲兴,

共32兲

where the ␭6 scaling follows from a study of the threefermion problem 关6兴. Figure 2 depicts a plot of ⌫ / ⌫0 as afunction of the interaction strength ␭. For m̄ = 0, the very
rapid growth of ⌫ / ⌫0 for small interaction strength arises
from the rapid increase in the ␭6 coefficient, while the drop
beyond the ferromagnetic transition point arises from the formation of a nonzero local magnetization in the trap that suppresses the product n↑共r兲n↓共r兲 in the integrand. These twocompeting effects lead to a peak in the rate of atom loss, viaconversion to molecules, at an interaction strength that isslightly beyond the ferromagnetic transition point. A measured peak in the atom loss rate as a function of increasing
interaction strength, as observed in Ref. 关17兴, would thus hint
at local ferromagnetism developing in the trap.
III. BEYOND THE LDA: MAGNETIZATION GRADIENTSThe discussion in Sec. II F has focused on the propertiesof the Fermi gas within the LDA. The energy functional atthis level of approximation does not have any gradient terms.
We will not worry about the shortcomings of this approximation for the density profile—it is well known that the
LDA breaks down near the trap edges—but instead focus on
going beyond the LDA by considering magnetization gradient terms with a view to studying the energetics of spin textures. We begin by noting that although we have been assign-

ing a global spin axis to the magnetization, the LDA energyfunctional would be unchanged if we, in fact, choose thelocal spin quantization axis to vary from point to point; only
the magnitude of the local magnetization plays a role. In
order to go beyond the LDA and to study the energies ofvarious spin textures in such a Fermi gas, we therefore needto extend the energy functional in two respects. First, wehave to promote the local magnetization to a vector quantity
so that the magnetization can point in different directions on
the Bloch sphere at different spatial locations. Second, we
have to include terms in the energy functional that depend on
the local magnetization gradients; this corresponds to adding
a surface tension term to the energy functional. The resultsfrom such an extended energy functional should be compared, in the future, with microscopic Hartree-Fock calculations.
We start with the dimensionless energy functional in Eq.
共20兲 and setn↑共r兲 =n共r兲
关1 + m共r兲兴,2

共33兲

n↓共r兲 =n共r兲关1 − m共r兲兴,2

共34兲

which define the local magnetization density m共r兲. As discussed, the spin quantization axis can be chosen to be different at each space point within the LDA. Let us next expand
the energy functional in powers of m共r兲; we will keep terms
up to m6共r兲 although terminating the expansion at m4共r兲
would not qualitatively affect our results. The energy functional then splits into two parts as␧ = ␧a关n共r兲兴 + ␧b关n共r兲,m共r兲兴,where
␧a关n共r兲兴 =

16
3␲2
+

␧b关n共r兲,m共r兲兴 =冕 冋冉 冊d 3r

6 n共r兲 5/3
5 2

共35兲

册

␭ 2
n 共r兲 − 共␥ − r2兲n共r兲 ,
3␲

16
3␲2

冕

共36兲

d3r关A2共r兲m2共r兲 + A4共r兲m4共r兲

+ A6共r兲m6共r兲 − hn共r兲m共r兲兴.

共37兲

Here ␧a only depends on the density profile that depends on
the interaction ␭ and that we assume is unchanged from thatgiven by the LDA calculation earlier. This is a good approximation since the corrections to the LDA energy are weak fortypical atom numbers used in experiments, as we will see
below. The coefficients of the magnetization-dependent energy functional, ␧b, are013607-5

A2共r兲 =

冉冊

n5/3共r兲 ␭ 2
n 共r兲 ,
−
22/33
3␲

共38兲

n5/3共r兲,
22/381共39兲

A4共r兲 =

PHYSICAL REVIEW A 80, 013607 共2009兲

LEBLANC et al.

A6共r兲 =

7n5/3共r兲
.
22/32187

A2共r兲 depends on ␭ explicitly. In addition, all the coefficients
A2,4,6 depend on the spatial location in the trap through thedensity and thus also depend implicitly on the interaction
strength ␭. This dependence was ignored in earlier phenomenological work on trapped fermions in an optical lattice关12兴. The Lagrange multiplier in the energy functionals are
given by ␥ = 共␥↑ + ␥↓兲 / 2 and h = 共␥↑ − ␥↓兲 / 2. Promoting themagnetization and the Lagrange multiplier h to vectors m
ជ , hជ
and including gradient terms lead to an energy functional of
the form␧b关n共r兲,mជ 共r兲兴 =

16
3␲2

␹−1
0 =

共40兲

冕 再

d3r A2共r兲兩mជ 共r兲兩2 + A4共r兲兩mជ 共r兲兩4

1+ A6共r兲兩m
ជ 共r兲兩6 + ␨s共r兲␣i关ⵜim j共r兲兴22冎− hជ 共r兲 · m
ជ 共r兲 .

Here, ␣i = 共␻i / ⍀兲2, which comes from our rescaling to an
isotropic problem. The stiffness ␨s共r兲 depends on r only
through the density n共r兲, and it can be computed in the uniform Fermi gas assuming that the magnetization variation isslow on the scale of the interparticle spacing but fast on thelength scale over which the total density varies, so that density variations can be ignored in this computation. The
Lagrange multiplier hជ 共r兲 must be chosen to satisfy global
constraints on the magnetization, for instance,兰d3rn共r兲mi共r兲 = 0 for each component i. We next outline thederivation of the stiffness term.

Rescaling distances to get an isotropic harmonic trapping
potential and setting M i = ␳N0 ␴共r = 0兲n共r兲mi共r兲, with r = R / RN0 ,we find␨s共r兲 ==For small magnetization, we can obtain the result for ␨sfrom the result for the magnetic susceptibility of the uniform
Fermi gas. Note that the excess energy in an applied fieldhជ 共q兲 共pointing in any direction兲 is given by ⌬E共q兲
= 21 ␹共q兲hi共q兲hi共−q兲 that defines the wave-vector dependent
magnetic susceptibility. This tells us that the magnetization
M共q兲 in this external field is simply M共q兲 = ␹共q兲h共q兲, so that
we can instead set ⌬E共q兲 = 21 ␹−1共q兲M i共q兲M i共−q兲. Expanding2
␹−1共q兲 = ␹−1
0 共1 + bq 兲 then yields
共42兲

The well-known result for a Fermi gas at T = 0 is that b= 1 / 12kF2 , using which the energy cost becomes, in real
space,

where

1
2␹0

冕 冋

ជ 共R兲兩2 +d3R 兩M1
12kF2

册ជ M i共R兲兴2 ,
关ⵜ

冋

1
1
n−1/3共r兲
2/3
0
2
2/3
2 3 6关3␲ n共r兲兴
␳N␴共RN0 兲3

册

2/3

1
.
72n共r兲共3N兲2/3

共45兲

For general values of the magnetization, higher order gradient terms might also become important. We will focus hereon the effects of this simplest gradient term in the energyfunctional.

Before proceeding to the energetics of various spin textures, let us slightly simplify the energy functional. Noticethat n共r兲 varies over the length scale of 1 in our dimensionless units. For large atom numbers, the stiffness is small as
seen from Eq. 共45兲 and we therefore expect significant variations in the magnetization to occur on length scales ᐉ Ⰶ 1 inour dimensionless units. Making this assumption, we can setជ 关n共r兲mi共r兲兴 ⬇ n共r兲ⵜជ mi共r兲, which results in the slightly simⵜ
plified energy functional
␧b =

共43兲

16
3␲2

冕 冋

d3r A2共r兲兩mជ 共r兲兩2 + A4共r兲兩mជ 共r兲兩4 + A6共r兲兩mជ 共r兲兩6+ hជ 共r兲 · m
ជ 共r兲 +

册n共r兲␣i关ⵜim j共r兲兴关ⵜim j共r兲兴 , 共46兲
144共3N兲2/3

where hជ 共r兲 is chosen to satisfy

冕A. Computation of the stiffness ␨s(r)

⌬E =

共44兲

B. Simplified magnetization energy functional

共41兲

2⌬E共q兲 = 21 ␹−1
0 共1 + bq 兲M i共q兲M i共− q兲.

␲ 2ប 2 ␲ 2ប 2
共3␲2␳兲−1/3 .
=MkF
M

d3rn共r兲mi共r兲 = 0

共47兲

for each component i 共for zero net magnetization兲. Recall
that ␣i = 共␻i / ⍀兲2, where ⍀ = 共␻x␻y␻z兲1/3 is the geometricmean of the trap frequencies.
IV. ENERGETICS OF SPIN TEXTURESThe energy functional we have derived above allows us to
study the energetics of various magnetization patterns in thetrapped Fermi gas. We begin by considering the isotropic
harmonic trap, for which we compare energies of a hedgehog
configuration and a domain-wall configuration of the magnetization. We then turn to an anisotropic cigar-shaped trap andshow how the hedgehog state gets deformed from the isotropic case. In each case, we begin by constructing the appropriate ansatz for the magnetization. We then numericallyminimize the resulting energy functional, by discretizing it
on a fine grid of points and using a simulated annealing
procedure, to obtain the optimal magnetization profile and its
energy. Assuming that the density and magnetization satisfy
the constraints that the total atom number is fixed and the013607-6

PHYSICAL REVIEW A 80, 013607 共2009兲

REPULSIVE FERMI GAS IN A HARMONIC TRAP:…

1.0
m





 

N=104

0.0

DW

1.0


HH



m

total energy



N=106







0.0 0.2 0.4 0.6 0.8 1.0 1.2
rtotal magnetization is zero 共so that the Lagrange multipliers
can be dropped兲, we can express the total energy as a sum¯␧ = ␧1 + ␧2, where

␧2 =

16
3␲2
+

16
3␲2

冕 冋冉 冊
d 3r

册

6 n共r兲 5/3 ␭n2共r兲
++ r2n共r兲 , 共48兲
5 2
3␲

冕 冋

d3r A2共r兲兩mជ 共r兲兩2 + A4共r兲兩mជ 共r兲兩4 + A6共r兲兩mជ 共r兲兩6册n共r兲
␣i关ⵜimជ j共r兲兴2 .
144共3N兲2/3

共49兲

A. Isotropic trap: Hedgehog state

For the isotropic trap, the density profile is sphericallysymmetric, which allows us to set
␧1 =

64
3␲

冕 冋冉 冊
dr r2

+

冕

dr r

2

再

Figure 3 shows the energy ¯␧HH of the hedgehog state obtained by finding the function m共r兲 that minimizes the hedgehog state energy. Figure 4共a兲 shows the magnetization profileof the hedgehog state at two different interaction strengths.
We find that the magnetization is suppressed in a small region around the trap center and vanishes at r = 0. To understand the magnetization profile of the hedgehog near the trap
center, we can focus just on the last two terms in Eq. 共51兲.
Taking a functional derivative with respect to m共r兲 and setting it to zero then suggest that m共r兲 ⬃ r2 at small r, so theenergy density coming from the central region of the hedgehog is finite. Far from the center, we expect the magnetization to be small. These expectations are consistent with themagnetization profiles shown in Fig. 4共a兲.

册

A2共r兲m 共r兲 + A4共r兲m 共r兲 + A6共r兲m 共r兲
2

冋

FIG. 4. 共a兲 Magnetization profiles for the hedgehog state at ␭
= 2.0 共solid兲 and ␭ = 2.4 共dashed兲. 共b兲 Magnetization profiles for thedomain-wall state at ␭ = 2.0 共solid兲 and ␭ = 2.4 共dashed兲. The profiles
have been calculated for 104 atoms in an isotropic trap. The hedgehog state has zero magnetization at the trap center while thedomain-wall state magnetization gets suppressed around the domain
wall but remains nonzero at the trap center. Insets indicate the schematic magnetization plot of the hedgehog state and the domain-wallstate.

B. Isotropic trap: Domain wall state

6 n共r兲 5/3 ␭n2共r兲
+ r2n共r兲 . 共50兲
+
5 2
3␲

For the magnetization-dependent energy functional, we mustset ␣i = 1 in the isotropic trap, and the hedgehog state corresponds to choosing mជ 共r兲 = m共r兲r̂. This leads to
64
␧2 =
3␲

0.5



FIG. 3. 共Color online兲 Dimensionless total energy, ␧1 + ␧2,
shown as a function of interaction strength, ␭, for an isotropic harmonic trap. DW indicates the energy of the domain-wall state for
104 atoms 共solid兲 and 106 atoms 共dashed兲. HH denotes ␧1 + ␧2 forthe hedgehog state that is nearly identical for 104 and 106 atoms.
Also shown 共thin solid line, “unmagnetized”兲 is ␧1, defined in Eq.
共44兲, which depends only on the total density profile.

␧1 =

(b)

0.0




(a)

0.5

4

冉 冊 册冎

n共r兲
dm共r兲 2m2共r兲2+
144共3N兲2/3
r2
dr

6

.

共51兲

We do not have to pay attention to the zero magnetization
constraint since this is guaranteed for any choice of m共r兲 by
the hedgehog ansatz symmetry. For typical particle numbers
in experiments, N ⬃ 104 – 106, the stiffness term has a verysmall coefficient. We will therefore assume that the averagedensity profile n共r兲 obtained from our earlier LDA calculation remains unchanged and only focus on changes in themagnetization profile arising from inclusion of gradientterms.

For the domain-wall state we have, as before,
␧1 =

64
3␲

冕 冋冉 冊dr r2

册

6 n共r兲 5/3 ␭n2共r兲
++ r2n共r兲 . 共52兲
5 2
3␲For the magnetization-dependent energy functional, we set
␣i = 1 in the isotropic trap and choose mជ 共r兲 = m共r兲ẑ for thedomain wall. This is capable of describing a state with spin↑, say, at the trap center with spin ↓ pushed away from thecenter, what we might call an in-out domain wall. We find
␧2 =

64
3␲
+

冕 冋

dr r2 A2共r兲m2共r兲 + A4共r兲m4共r兲 + A6共r兲m6共r兲

冉 冊册

n共r兲
dm共r兲 2,
144共3N兲2/3 dr

共53兲

where, for N↑ = N↓, we must satisfy the constraint
兰dr r2n共r兲m共r兲 = 0. Figure 3 shows the energy ¯␧DW of thedomain-wall state obtained by finding the function m共r兲 thatminimizes its energy subject to the zero magnetization con-013607-7

PHYSICAL REVIEW A 80, 013607 共2009兲

LEBLANC et al.
 

straint. We find that this domain-wall state has a higher energy than the hedgehog state. Figure 4共b兲 shows the magnetization profile of the domain-wall state. As expected, themagnetization is suppressed in a small region around thedomain wall but remains nonzero at the trap center.



C. Cigar-shaped trap: Distorted hedgehog

If we consider a cylindrically symmetric 共cigar-shaped兲
trap, we can look at an ansatz of the formm
ជ 共r兲 = m共␳,z兲

冉

冊

x
y
sin ␺, sin ␺,cos ␺ ,␳␳ 

m2 2
sin ␺
␳2

+ m2关␣⬜共⳵␳␺兲2 + ␣z共⳵z␺兲2兴,共56兲

so the integral 兰d3r → 2␲兰dzd␳␳. We can assume that m isan even function of z and that ␺共␳ , −z兲 = ␲ − ␺共␳ , z兲 关so thatsin2 ␺共␳ , −z兲 = sin2 ␺共␳ , z兲兴 to restrict the energy integration
grid to just z ⬎ 0. These conditions ensure that the total magnetization integrates to zero. The final expression for the energy can thus be recast, with r ⬅ 冑␳2 + z2, as
␧1 =

64
3␲

冕 冋冉 冊
冕 冕冑

␧2 =

64
3␲

册

6 n共r兲 5/3 ␭n2共r兲
++ r2n共r兲 , 共57兲
5 2
3␲

dr r2

Rmax

2
rmax
−z2

dz

0

d␳␳F共␳,z兲,

0F共␳,z兲 = A2m2 + A4m4 + A6m6 +

再

共58兲

m2 2
n
␣
⬜ 2 sin ␺
144共3N兲2/3␳冎+ ␣⬜共⳵␳m兲2 + ␣z共⳵zm兲2 + m2关␣⬜共⳵␳␺兲2 + ␣z共⳵z␺兲2兴 ,共59兲
with ␺共␳ = 0 , z兲 = 0 and ␺共␳ , z = 0兲 = ␲ / 2 by symmetry. For
notational simplicity, we have suppressed the coordinate labels on n , m , ␺ in the above functional.
We find, numerically, that ␺ ⬇ ␪, so, in fact, the ansatz
simplifies to the formm
ជ 共r兲 = m共␳,z兲

冉 冊

x y z
, , .
r r r



共55兲

With the choice of magnetization in Eq. 共54兲, we have兩m
ជ 共r兲兩2 = m2共␳ , z兲, while

␣i共⳵im j兲共⳵im j兲 = ␣⬜共⳵␳m兲2 + ␣z共⳵zm兲2 + ␣⬜





冊

共60兲

The main effect of going from the spherical to the cigarshaped trap is that the magnitude of the magnetization is no





where ␳ = 冑x2 + y 2 and ␺ ⬅ ␺共␳ , z兲. For ␺ = ␪ = tan−1共␳ / z兲 thisreduces to the spherical hedgehog ansatz. Note that the direction of the magnetization on the Bloch sphere is unrelated
to the location in real space. We could equally well havechosen, for instance,

冉





共54兲

x
y
m
ជ 共r兲 = m共␳,z兲 cos ␺, sin ␺, sin ␺ .
␳
␳































FIG. 5. 共Color online兲 Two-dimensional magnetization profile
for the distorted hedgehog showing breakdown of the LDA for themagnetization density for 102 atoms in an anisotropic trap with ␭= 2.4 and ␣⬜ = 2 , ␣z = 0.25共␻⬜ / ␻z ⬇ 2.8兲. 共a兲 Plot of the magnitudeof the magnetization m共␳ , z兲 and equal-magnetization contours displayed in rescaled coordinates in which the trap potential is spherically symmetric. Colorbar to the right indicates the value of m共␳ , z兲.
We see that m共␳ , z兲 is larger in magnitude for larger values of ␳ than
it is for z, indicating that the surface tension makes it easier tochange its value in the weak trapping direction. 共b兲 mជ shown as aquiver plot indicating the magnitude and direction of the magnetization 共plotted in coordinates where the trap anisotropy is explicitly
shown兲. Shaded area indicates the region of the trap where the atom
density is nonzero.

longer just dependent on the radial coordinate r. The magnetization however still points 共in our rescaled coordinates兲
along the radial direction. The plot of the magnetization for␭ = 2.4 in the rescaled and in the original coordinates for a
trap anisotropy corresponding to ␣⬜ = 2 , ␣z = 0.25 共a trap frequency ratio ␻⬜ / ␻z ⬇ 2.8兲 is given in Fig. 5.
V. EFFECT OF DIPOLAR INTERACTIONS

Our results for the spin texture energetics and magnetization profiles have been obtained by neglecting the role of thelong range magnetic dipole interaction between the fermions.
The dipole interaction will add to the magnetic energy ofatoms in the trap. In addition, it will lead to spatial variationsin the magnetic field seen by atoms within the trap and, thus,
tend to cause dephasing as atoms in different regions will
precess at different rates. Such effects are known to be important in 87Rb spinor Bose condensates 关18,19兴. In order toestimate the dipole interaction energy and the time scale ofthis dephasing, we have considered the spatial variations in
the dipole field for the simple case of the spherical trap.
The expression for the precession frequency at distance r
from the center of the spherical trap is, for the hedgehog
state,

013607-8

PHYSICAL REVIEW A 80, 013607 共2009兲

REPULSIVE FERMI GAS IN A HARMONIC TRAP:…

dip
␯HH
共r兲 =

冉 冊冉 冊1 ␮0 2 2冑2␮
h 4␲ B ␲冑3
⫻

冕

冋

ប
M⍀

−3/2

N1/2

dr1d␪r21 sin ␪m共r1兲n共r1兲
共r2 + r21 − 2rr1 cos ␪兲3/2

F共r,r1, ␪兲 = cos ␪ − 3

F共r,r1, ␪兲, 共61兲

共r1 − r cos ␪兲共r1 cos ␪ − r兲
共r2 + r21 − 2rr1 cos ␪兲

册

,
共62兲

where ␮B ⬇ 9.27⫻ 10−24 J / T is the Bohr magneton and ␮0= 4␲ ⫻ 10−7 N / A2 is the permeability of free space. Evaluating this, we find that the typical values of 共and also thevariation in兲 the precession frequency for ␭ = 2.5, N = 106,
dip 6
dip 40
共 Li兲 ⬇ 0.03 Hz and ␯HH
共 K兲
and ⍀ = 2␲共170 Hz兲 are ␯HH⬇ 0.6 Hz. The energy associated with the dipole interactions
is far smaller than our estimate of magnetic exchange energies, ⬃500 Hz, arising from the s-wave contact interaction
between fermions 共in the interaction range where we expect
ferromagnetism兲. At the same time, measurements of thetypical atom lifetime, ␶, on the repulsive side of the Feshbach resonance indicate that ␶ ⬃ 10 ms for 40K 关20兴 and ␶⬃ 100 ms for 6Li 关2,21兴. These are clearly much less thanthe variations in the precession period induced by spatial
variations in the dipolar field as estimated above. Taken together, these estimates show that ignoring the effect of dipole
interactions is a very good approximation in this system.
VI. CONCLUSIONSIn conclusion, we have studied ferromagnetism and spin
textures in ultracold atomic Fermi gases in the regime ofstrongly repulsive interactions using the LDA extended toinclude magnetization gradient corrections. Within the LDA
at zero temperature, we have shown that the release energyof the gas, as well as its separate kinetic energy and interaction energy components, shows a sharp signature of the ferromagnetic transition. We have also shown that the atom loss
rate via three-body collisions has a peak very close to theferromagnetic transition and it provides yet another diagnostic of the transition into the ferromagnetic state. We havegone beyond the LDA by deriving a surface tension correction to the energy functional, which depends on atom number
and the trap geometry and used it to study the energetics of

关1兴 M. Greiner, C. A. Regal, and D. S. Jin, Nature 共London兲 426,
537 共2003兲.
关2兴 T. Bourdel, J. Cubizolles, L. Khaykovich, K. M. F. Magalhäes,
S. J. J. M. F. Kokkelmans, G. V. Shlyapnikov, and C. Salomon,
Phys. Rev. Lett. 91, 020402 共2003兲.
关3兴 S. Gupta, Z. Hadzibabic, M. W. Zwierlein, C. A. Stan, K.
Dieckmann, C. H. Schunk, E. G. M. van Kempen, B. J. Verhaar, and W. Ketterle, Science 300, 1723 共2003兲.
关4兴 M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H.
Schunck, and W. Ketterle, Nature 共London兲 435, 1047 共2005兲.

various spin textures in a two-component trapped Fermi gas.
For a spherically symmetric trap, we find that a hedgehog
magnetization profile has lower energy than a domain-wallstate. For large atom numbers, the small surface tension
leads to a small energy difference between the two spin textures and the results are close to those of the LDA. In thiscase, the surface tension is responsible for selecting thehedgehog state as having the lowest energy but we havechecked that it does not significantly change our results forthe release energy and the atom loss rates. These continue tobe useful, albeit indirect, diagnostics of the transition into theferromagnetic state. For elongated clouds, we have shownthat the surface tension term distorts the hedgehog states, inrescaled coordinates where the trap is isotropic, in such amanner that the magnitude of the magnetization changes
more easily in the weak direction of the trap than would be
expected on the basis of the LDA. Such a breakdown of the
LDA is more apparent for smaller atom numbers. Finally, we
have considered the effect of magnetic dipolar interactions
on our results and find that it is a good approximation toignore dipole interactions in this system.
The typical atom loss rate near the Feshbach resonancesets a constraint that the formation time for the ferromagnetic
state will have to be on the order of tens of milliseconds for40
K 关20兴 and hundreds of milliseconds for 6Li 关2,21兴 in orderfor it to be observed. A direct way to probe the spin textures
discussed here would be through high resolution in situ magnetometry, as has been done for spinor Bose condensates关22兴. An experimental observation of ferromagnetism intrapped Fermi gases would provide impetus for future theoretical work on finite temperature effects and collectivemodes in the strongly interacting regime.

ACKNOWLEDGMENTS

We thank Rembert Duine, Gyu-Boong Jo, Wolfgang
Ketterle, and Allan MacDonald for useful discussions. This
work was supported by the NSERC of Canada 共J.H.T., A.P.,
and A.A.B.兲 and the Canadian Institute for Advanced Research 共J.H.T. and L.J.L.兲. A.P. acknowledges support fromthe Sloan Foundation, the Connaught Foundation, and theOntario ERA. J.H.T. thanks the MIT-Harvard Center for
Ultracold Atoms for hospitality during the preparation of thispaper.

关5兴 E. C. Stoner, Philos. Mag. 15, 1018 共1933兲.
关6兴 D. S. Petrov, Phys. Rev. A 67, 010703共R兲 共2003兲.
关7兴 L. Salasnich, B. Pozzi, A. Parola, and L. Reatto, J. Phys. B 33,
3943 共2000兲.
关8兴 T. Sogo and H. Yabu, Phys. Rev. A 66, 043611 共2002兲.
关9兴 R. A. Duine and A. H. MacDonald, Phys. Rev. Lett. 95,
230403 共2005兲.
关10兴 See, for instance, H. Tasaki, Prog. Theor. Phys. 99, 489
共1998兲.
关11兴 S. Zhang, H.-H. Hung, and C. Wu, e-print arXiv:0805.3031.

013607-9

PHYSICAL REVIEW A 80, 013607 共2009兲

LEBLANC et al.
关12兴 I. Berdnikov, P. Coleman, and S. H. Simon, Phys. Rev. B 79,
224403 共2009兲.
关13兴 G. B. Partridge, W. Li, Y. A. Liao, R. G. Hulet, M. Haque, and
H. T. C. Stoof, Phys. Rev. Lett. 97, 190407 共2006兲.
关14兴 T. N. De Silva and E. J. Mueller, Phys. Rev. Lett. 97, 070402
共2006兲; S. S. Natu and E. J. Mueller, e-print arXiv:0802.2083;
S. K. Baur, S. Basu, T. N. De Silva, and E. J. Mueller, Phys.
Rev. A 79, 063628 共2009兲.
关15兴 A. Imambekov, C. J. Bolech, M. Lukin, and E. Demler, Phys.
Rev. A 74, 053626 共2006兲.
关16兴 R. Sensarma, W. Schneider, R. B. Diener, and M. Randeria,
e-print arXiv:0706.1741.
关17兴 G.-B. Jo, Y. Lee, C. Christensen, J. Choi, J. H. Thywissen, D.

Pritchard, and W. Ketterle 共unpublished兲.
关18兴 M. Vengalattore, S. R. Leslie, J. Guzman, and D. M. StamperKurn, Phys. Rev. Lett. 100, 170403 共2008兲; M. Vengalattore, J.
Guzman, S. Leslie, F. Serwane, and D. M. Stamper-Kurn,e-print arXiv:0901.3800.
关19兴 R. W. Cherng and E. Demler, e-print arXiv:0806.1991.
关20兴 C. A. Regal, M. Greiner, and D. S. Jin, Phys. Rev. Lett. 92,
083201 共2004兲.
关21兴 K. Dieckmann, C. A. Stan, S. Gupta, Z. Hadzibabic, C. H.
Schunck, and W. Ketterle, Phys. Rev. Lett. 89, 203201 共2002兲.
关22兴 L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore, andD. M. Stamper-Kurn, Nature 共London兲 443, 312 共2006兲.

013607-10