First and second sound modes at finite temperature in trapped Fermi gases from BCS to BEC
Yan He, Qijin Chen, Chih-Chun Chien, and K. Levin
James Franck Institute and Department of Physics, University of Chicago, Chicago, Illinois 60637
(Dated: August 11, 2021)
We determine the temperature T dependence of first and second sound mode frequencies for trapped Fermi
gases undergoing BCS to Bose Einstein condensation (BEC) crossover. Our results are based on the two fluid
equations in conjunction with a microscopic calculation of thermodynamical variables. As in experiment and at
unitarity, we show that the lowest radial breathing mode is T independent. At finite T , higher order breathing
modes strongly mix with second sound. Their complex T dependence should provide an alternative way of
measuring the transition temperature, Tc.
PACS numbers: 03.75.Hh, 03.75.Ss, 74.20.-z arXiv:0704.1889
The recent discovery of the superfluid phases of trapped
Fermi gases has led to considerable interest in their collective
mode spectrum [1, 2, 3, 4, 5, 6, 7]. Among the modes of ex-
perimental interest are breathing modes as well as propagating
first sound. While originally theoretical attention [8, 9, 10]
was focused on ground state properties experimental mea-
surements are naturally not restricted to temperature T = 0.
Indeed there is an interesting body of information which is
emerging in these Fermi gases about the finite temperature
behavior [4, 5, 6] of the breathing modes and, more recently
about the propagating sound velocity [7].
The purpose of this paper is to compute sound mode fre-
quencies in spherically trapped Fermi gases undergoing BCS
to Bose-Einstein condensation (BEC) crossover, at general T .
We present a solution of the linearized two fluid equations
and compare with recent experiments. We focus on the radial
breathing modes and present predictions for second sound, as
well. The structure of the two fluid equations for Bose [11]
and Fermi gases [12] has been rather extensively discussed.
In the crossover regime, the normal fluid is novel [13, 14],
containing both fermions and noncondensed pairs, which have
not been systematically included in previous collective mode
literature.
Of great importance to the field as a whole is the future
possibility of second sound observations, in large part because
this may ultimately help assign more precise experimental val-
ues to the transition temperature Tc. While existing experi-
mental techniques such as vortex observation [15] and den-
sity profile features [16] help establish superfluidity they pro-
vide lower bounds on Tc or determine its value for the special
case of a population imbalanced system. Thermodynamical
experiments measure Tc more directly [17, 18] but have been
confined to unitarity. Thus other techniques, such as second
sound observation will be of great value. One of the princi-
pal results of the present paper is an analysis of how second
sound is coupled to the breathing modes. We demonstrate that
higher order breathing modes will reveal Tc through this cou-
pling, and therefore are an alternative to direct second sound
measurements. However, the lowest breathing mode appears
to be remarkably T independent at unitarity. This has been ob-
served experimentally [4, 5] and argued to follow from isen-
tropic considerations [6]. Here we show, that even when we
treat the full coupling between first and second sound, we ob-
tain similar T independent behavior at unitarity.
At the core of the two fluid theory is the assumption that
hydrodynamics is valid and that there are frequent collisions
which produce a state of local thermodynamical equilibrium.
Although there are some exceptions [19], reaching the two
fluid regime has not been easy for atomic Bose gases. Two
fluid dynamics are more readily achieved for Fermi gases,
principally because in the crossover regime the large scatter-
ing lengths produce sufficient collisions. Nevertheless there
has been considerable theoretical interest in setting up [12]
and solving [20] the two fluid equations for Bose conden-
sates. Indeed, hydrodynamical approaches have successfully
addressed both the T = 0 and normal state regimes of the
Bose gases [21] and demonstrated that the breathing mode
frequencies are the same. Here, by contrast we address the
Fermi gas case in a trap. Because they interact more strongly
near unitarity, hydrodynamical descriptions have been argued
quite convincingly [6, 9, 22] to be valid.
Previous theoretical work has been confined to T = 0 treat-
ments of a harmonic trap, or to finite T theories [22] of a
uniform gas. Our work is most similar in spirit to an earlier
Bose gas study [20] although we introduce different numeri-
cal techniques as well as address fermions rather than bosons.
We note that the input thermodynamics of systems undergo-
ing the BCS-BEC crossover which is used in the present pa-
per, has been rather well calibrated against experimental mea-
surements in Ref. [17] and is based on a finite temperature
extension of the simplest (BCS-Leggett) ground state. In the
absence of a trap our results are for the most part similar to
those in Ref. [22].
We begin with the two fluid equations which describe the
dynamical coupling of the superfluid velocity vs and the nor-
mal velocity vn. Just as in the spirit of the original Landau
two fluid equations we ignore viscosity terms. In the presence
of a trap potential Vext =
mωhor
2, the two fluid equations
are given by m∂vs
+∇(µ+Vext+m
) = 0, ∂j
+∇·Π =
−n∇Vext,
+ ∇ · j = 0 and ∂(ns)
+ ∇ · (nsvn) = 0,
with Πij = pδij + nsv
s + nnv
n, n = ns + nn and
http://arxiv.org/abs/0704.1889v2
j = nsvs + nnvn. Here µ is the chemical potential, p the
pressure, and s the entropy per particle. Moreover, we have
nnvn + nsvs = nv. Here ns (vs) and nn (vn) denote the
superfluid and normal densities (velocities), respectively. We
use the subscript “0” to denote equilibrium quantities such that
vs0 = vn0 = 0, ∇(µ + Vext) = 0, ∇p0 = −n0∇Vext, and
n0, s0, µ0, p0 are independent of time t. Combining this with
the thermodynamic relation dµ = −s dT + dp/n, we have
∇T0 = 0, implying that temperature T0 is constant in the trap.
It then follows that in equilibrium µ = µ0 − Vext, consistent
with the Thomas-Fermi approximation.
For small deviations from equilibrium, we may linearize the
two fluid equations. Eliminating the velocities vs and vn, one
finds [12, 20]
= ∇ ·
−∇ · (δT n0∇s0) (1)
ns0n0s
− (∇s0)
2δT +∇s0 · ∇
We will focus on δµ(r) and δT (r) as the principal variables.
This choice, which is different from that in Ref. [20], is made
because both variables are non vanishing at the trap edge so
that in a basis set expansion they will satisfy the same bound-
ary conditions. Moreover, the two fluid equations are simplest
in this form. Expressing δs, δp and δn in terms of δµ, δT , the
two fluid equations can be rewritten as
, (3)
, (4)
A = ∇ · (n∇δµ) +∇ · (ns∇δT ) ,
s2∇δT
+∇s · ∇δµ+ s∇s · ∇δT .
It is understood that all coefficients of δµ and δT are calcu-
lated in equilibrium so that we drop the subscript “0”. The
thermodynamical quantities in equilibrium can be calculated
following Ref. [23], based on the standard local density ap-
proximation (LDA), µ(r) = µ0 − Vext(r). Their derivatives
with respect to T and µ can be calculated analytically, and
their gradients can be obtained via ∇f = −
∇Vext,
where f denotes any of the variables (n, ns, nn, s).
To solve the two coupled differential equations (3) and (4),
we assume a simple harmonic time dependence δµ, δT ∝
e−iωt. We cast the differential two fluid equations into an
eigenfunction problem with ω2 playing the role of eigenvalue
and the eigenfunctions given by the amplitudes of δµ and δT .
Since neither T nor µ depends on the density, they will not
vanish at the trap edge. Our boundary conditions require that
all thermodynamic variables be smooth (but not necessarily
zero) at the trap edge. At finite T , the density in the trap de-
creases exponentially when the local chemical potential be-
comes negative at large radius. We choose, thus, to expand
δµ and δT in terms of Jacobi polynomials. For our numerics
we choose the matrix dimension to be 300; we have similarly
investigated matrices of dimension 200 up to 900, and found
little change in our principal findings.
We now turn to an important aspect of our numerics. Be-
cause we generate some 300 frequencies in our numerical ap-
proach it is essential to establish a mechanism for systemati-
cally identifying first and second sound modes. To help find
such a “fingerprint”, we introduce a “decoupling approxima-
tion” based on reducing Eqs. (3) and (4) to δ̈µ = g1(δµ, δT ),
δ̈T = g2(δµ, δT ), where g1,2 are known functions. We elim-
inate cross terms by setting δT = 0 in g1 and δµ = 0 in
g2. With these two decoupled equations it is then relatively
straightforward to associate a profile plot of the numerically
calculated δp(r) and nδs(r) vs r within a trap with first or sec-
ond sound-like modes. Here we convert δµ and δT to δp and
δs via δp = nδµ + nsδT and δs =
We stress that this procedure differs somewhat from that in
Ref. [20] where it was argued that one could associate second
sound with a mode in which the oscillations in the thermody-
namical variable δT (r) were constrained to be in the conden-
sate region of the trap. As explained below, we do not find
this to be the case. An important check on our procedure is
that we find that there is no sign of second sound above Tc.
Up to this point everything is general applying to both
Fermi and Bose superfluids. All that is needed is a micro-
scopic theory for thermodynamical variables. Here we use a
calculational framework we have developed for treating BCS-
BEC crossover in trapped Fermi gases [23, 24] which em-
phasizes the importance of pseudogap effects or finite mo-
mentum pairs. The local thermodynamical potential (density)
Ω = Ωf + Ωb is associated with a contribution from gapped
fermionic excitationsΩf as well as from non-condensed pairs,
called Ωb. We have
Ωf = −
[(ξk − Ek)− 2T ln(1 + e
−Ek/T )],
Ωb = −Z∆
2µpair +
T ln(1 − e−Ωq/T ) . (5)
Here µpair is the chemical potential of the pairs which is
zero below Tc and the pair dispersion Ωq, along with the (in-
verse) residue Z , can be derived from a microscopic T -matrix
theory, described elsewhere [13, 14]. Using Ω one then ar-
rives at thermodynamical properties such as the entropy den-
sity ns = −∂Ω
as well as self consistent equations for the
total excitation gap ∆, the contribution to ∆ from noncon-
densed pairs (called the pseudogap) and the number equa-
tions. These self-consistent (local) equations are simply given
0 0.1 0.2 0.3 0.4
0 0.2 0.4 0.6 0.8 1
r/Rmax
T = 0.1
n = 1
T = 0.2
n = 1
(a) (b)
(c) (d)
Figure 1: (Color online) Behavior of the first (upper row) and second
sound (lower row) modes for a spherical trap at unitarity within the
“decoupling approximation” (see text). The left column shows the
T -dependence of the frequencies, while the right shows correspond-
ing typical spatial oscillation profiles for δp [(black) solid lines] and
nδs [(red) dashed lines] at the lowest frequencies, which provide fin-
gerprints of first and second sound. For the second sound, nδs dom-
inates and changes sign within the condensate. Here Tc ≈ 0.27TF ,
and the arrows indicate the condensate edge.
by ∂Ω
= 0, ∂Ω
∂µpair
= 0 and n = −∂Ω
, subject to the total
number constraint N =
d3r n(r). When T < Tc, we can
use the gap equation and chain rule to eliminate the variable
∆, which is a function of µ and T , via ∂∆
and ∂∆
. Similarly, when T > Tc we use
the gap and pseudogap equations to eliminate the variables ∆
and µpair to arrive at thermodynamical quantities.
Figure 1 shows the lowest two collective modes at unitar-
ity in a spherical trap, obtained by solving the chemical po-
tential or temperature fluctuation equation in the decoupling
approximation scheme described earlier. It is evident that our
approximated breathing mode frequency [Fig. 1(a)] is inde-
pendent of temperature. We understand this result by noting
that the decoupled equation for the breathing mode is given by
−ω2δµ = Cµµ1∇
2δµ+Cµµ2 · ∇δµ where Cµµ1 = n
and Cµµ2 = ∇µ. At unitarity, n
µ. The only
T -dependence contained in µ0 ≡ µ(r = 0) can be elimi-
nated via a simple re-scaling of r → r
2µ0/mω
ho, yielding
a T -independent breathing mode frequency. These arguments
can be shown to be equivalent to the isentropic assumption of
Ref. [6], where vs = vn is assumed in the two fluid model,
leading to two simpler coupled equations associated with the
Euler and the continuity equations.
By contrast, the second sound mode frequency we obtain
increases rapidly with temperature. Some typical oscillation
profiles of δp(r) and nδs(r) are shown in the right two pan-
els of Figure 1. Although the “entropy density oscillations”
nδs(r) fall off at large r, the entropy per particle δs(r) os-
cillations (not shown) increase very rapidly upon entering the
normal region. Consequently temperature fluctuations δT be-
0 0.2 0.4 0.6 0.8
0 0.2 0.4 0.6 0.8 1
T/TF = 0.01 T/TF = 0.2
1/kFa = 0
T/TF = 0.15
1/kFa = 0
T/TF = 0.08
1/kFa = 1
r/Rmax
1/kFa = 0
(a) (b)
(c) (d)
Figure 2: (Color online) Typical spatial oscillation profiles for δp
[(black) solid lines] and nδs [(red) dashed lines] obtained from the
fully coupled equations for a spherical trap at (a-c) unitarity and (d)
1/kF a = 1 at different T , for the first (top row) and second (bottom
row) sound modes. Also shown in (a) is the T = 0 analytical result
[(green) dot-dashed curve]. The arrows indicate condensate edge.
come large at the trap edge.
Our identification of first sound for the decoupled case leads
us to associate this mode in a coupled situation with a profile
for which at the trap center δp has a large amplitude, while
nδs is almost zero (with a small peak near the trap edge) as
in Fig. 1(b). By contrast, in the trap center the second sound
mode has large entropy fluctuations, while the pressure fluctu-
ations are almost zero [Fig. 1(d)]. These features will serve as
fingerprints for distinguishing (lower order) first and second
sound modes at finite temperatures.
Figure 2 shows some typical eigenfunction profiles of the
lowest modes obtained in the spherical trap upon solution of
the fully coupled two fluid equations. The first row corre-
sponds to the breathing mode in the unitary case. The good
agreement between the very low T result for δp and the T = 0
analytical solution (green dashed line) in Fig. 2(a) helps vali-
date our numerical scheme. Figure 2(b) corresponds to a high
temperature breathing mode. In this regime the pseudogap re-
gion outside the superfluid core is relatively large and the peak
in nδs(r) is accordingly very broad. The lower two panels
correspond to second sound modes for the unitary and BEC
cases. By contrast with the breathing modes, here nδs has a
larger amplitude than δp with an opposite sign. Clearly this
is very similar to what we observed in the decoupled mode
analysis of Fig. 1.
Figures 3(a-c) address the fully coupled equations and show
the behavior of the lowest breathing (upper branch) and sec-
ond sound mode (lower branch) frequencies as a function of
temperature in a spherical trap for 1/kFa = 1, 0, and −0.5,
respectively. For all three values of 1/kFa we find very little
sign of Tc in the lowest breathing mode frequency. In con-
trast, the second sound mode frequencies increase with T and
disappear above Tc. Figure 3(d) presents a more complete
series of modes for the unitary case. Here the lines serve as
guides to the eye for the lowest (blue) and higher order (green)
0 0.1 0.2 0.3
0 0.1 0.2 0.3 0.4
1/kFa = 1
1/kFa = 0
1/kFa = -0.5
1/kFa = 0
(a) (b)
(c) (d)
Figure 3: (Color online) Temperature dependence of breathing mode
and second sound frequencies in a spherical trap. Panels (a-c) are for
the near-BEC, unitary, and near-BCS cases, respectively. The upper
and lower branches in (a-c) represent the lowest breathing mode and
second sound frequencies, respectively. In (d), more complete results
(open circles) are shown at unitarity, where the lines serve as guides
to the eye for the breathing mode (nearly horizontal blue and green
lines) and second sound (increasing red curves) frequencies.
breathing mode, and second sound (red) frequencies.
At unitarity [Fig. 3(d)] one can identify a sequence of
higher order breathing modes which precisely overlap with
analytical calculations for T = 0. Importantly, only the low-
est of these is found to be a constant in temperature; the others
are found to mix with second sound modes, as indicated in the
solid and dashed (red) lines in Fig. 3(d). The behavior of the
lowest mode helps to justify the isentropic assumption made
in Ref. [6]. We understand this by referring back to the de-
coupled profiles in Figs. 1(b) and 1(c), which are seen to be
quite distinct. By contrast the profiles of the decoupled first
and second sound modes at higher order (not shown) appear
more similar to each other than their lowest order counter-
parts. Indeed, the behavior of the profiles at higher order is
associated with an increasing number of nodes in the curves
of δp and nδs (not shown), which leads to a greater similarity
between first and second sound profiles and helps explain why
the higher order modes are more readily mixed [25].
In summary, we have presented predictions for future ex-
periments on higher breathing modes and second sound in a
trap. We find that only the lowest breathing mode frequency
has very weak T dependence. For the unitary case this tem-
perature insensitivity was clearly observed by Thomas and co-
workers [6]. As a result of this experiment it should not be
surprising that we find relatively weak dependencies on ei-
ther side of the Feshbach resonance for this breathing mode
frequency. In the literature there are experimental claims (at
1/kFa = 1.0) which are consistent with a decrease [4] in the
radial breathing mode frequency, as indeed we find here, al-
though ours is probably too weak to reconcile the different
findings in Refs. [4] and [2, 3]. Finally, our more complete
studies at unitarity show that if higher order radial breathing
modes can be accessed, because of their strong hybridization
with second sound, it may be possible to use these breathing
modes, rather than direct second sound to determine the tran-
sition temperatures Tc. This should be of value since there are
currently few experiments which can assign a value to Tc over
the wider crossover regime.
This work was supported by Grant Nos. NSF PHY-
0555325 and NSF-MRSEC DMR-0213745; we thank John
Thomas for suggesting this problem and Cheng Chin for use-
ful conversations.
[1] C. Chin, M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim,
J. Hecker-Denschlag, and R. Grimm, Science 305, 1128 (2004).
[2] J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and J. E.
Thomas, Phys. Rev. Lett. 92, 150402 (2004).
[3] J. Kinast, A. Turlapov, and J. E. Thomas, Phys. Rev. A 70,
051401(R) (2004).
[4] A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R. Geursen,
M. Bartenstein, C. Chin, J. H. Denschlag, and R. Grimm, Phys.
Rev. Lett. 98, 040401 (2007).
[5] J. Kinast, A. Turlapov, and J. E. Thomas, Phys. Rev. Lett. 94,
170404 (2005).
[6] J. E. Thomas, J. Kinast, and A. Turlapov, Phys. Rev. Lett. 95,
120402 (2005).
[7] J. Joseph, B. Clancy, L. Luo, J. Kinast, A. Turlapov, and J. E.
Thomas (2006), cond-mat/0612567.
[8] H. Hu, A. Minguzzi, X.-J. Liu, and M. P. Tosi, Phys. Rev. Lett.
93, 190403 (2004).
[9] H. Heiselberg, Phys. Rev. Lett. 93, 040402 (2004).
[10] G. E. Astrakharchik, R. Combescot, X. Leyronas, and
S. Stringari, Phys. Rev. Lett. 95, 030404 (2005).
[11] A. Griffin and E. Zaremba, Phys. Rev. A 56, 4839 (1997).
[12] E. Taylor and A. Griffin, Phys. Rev. A 72, 053630 (2005).
[13] Q. J. Chen, I. Kosztin, B. Jankó, and K. Levin, Phys. Rev. Lett.
81, 4708 (1998).
[14] Q. J. Chen, J. Stajic, S. N. Tan, and K. Levin, Phys. Rep. 412, 1
(2005).
[15] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, and W. Ket-
terle, Nature 435, 170404 (2005).
[16] M. W. Zwierlein, A. Schirotzek, C. H. Schunck, and W. Ket-
terle, Science 311, 492 (2006).
[17] J. Kinast, A. Turlapov, J. E. Thomas, Q. J. Chen, J. Stajic, and
K. Levin, Science 307, 1296 (2005), published online 27 Jan-
uary 2005; doi:10.1126/science.1109220.
[18] L. Luo, B. Clancy, J. Joseph, J. Kinast, and J. E. Thomas, Phys.
Rev. Lett. 98, 080402 (2007).
[19] D. M. Stamper-Kurn, H.-J. Miesner, S. Inouye, M. R. Andrews,
and W. Ketterle, Phys. Rev. Lett. 81, 500 (1998).
[20] V. B. Shenoy and T.-L. Ho, Phys. Rev. Lett. 80, 3895 (1998).
[21] A. Griffin, W.-C. Wu, and S. Stringari, Phys. Rev. Lett. 78, 1838
(1997).
[22] H. Heiselberg, Phys. Rev. A 73, 013607 (2006).
[23] Q. J. Chen, J. Stajic, and K. Levin, Phys. Rev. Lett. 95, 260405
(2005).
[24] J. Stajic, Q. J. Chen, and K. Levin, Phys. Rev. Lett. 94, 060401
(2005).
[25] It is interesting to note that there exists a lower frequency first-
sound like mode at finite T , which we also find to be present in
the “free” (but hydrodynamic) Fermi gas limit at T 6= 0.

  We determine the temperature $T$ dependence of first and second sound mode
frequencies for trapped Fermi gases undergoing BCS to Bose Einstein
condensation (BEC) crossover. Our results are based on the two fluid equations
in conjunction with a microscopic calculation of thermodynamical variables. As
in experiment and at unitarity, we show that the lowest radial breathing mode
is $T$ independent. At finite $T$, higher order breathing modes strongly mix
with second sound. Their complex $T$ dependence should provide an alternative
way of measuring the transition temperature, $T_c$.

First and second sound modes at finite temperature in trapped Fermi gases from BCS to BEC
Yan He, Qijin Chen, Chih-Chun Chien, and K. Levin
James Franck Institute and Department of Physics, University of Chicago, Chicago, Illinois 60637
(Dated: August 11, 2021)
We determine the temperature T dependence of first and second sound mode frequencies for trapped Fermi
gases undergoing BCS to Bose Einstein condensation (BEC) crossover. Our results are based on the two fluid
equations in conjunction with a microscopic calculation of thermodynamical variables. As in experiment and at
unitarity, we show that the lowest radial breathing mode is T independent. At finite T , higher order breathing
modes strongly mix with second sound. Their complex T dependence should provide an alternative way of
measuring the transition temperature, Tc.
PACS numbers: 03.75.Hh, 03.75.Ss, 74.20.-z arXiv:0704.1889
The recent discovery of the superfluid phases of trapped
Fermi gases has led to considerable interest in their collective
mode spectrum [1, 2, 3, 4, 5, 6, 7]. Among the modes of ex-
perimental interest are breathing modes as well as propagating
first sound. While originally theoretical attention [8, 9, 10]
was focused on ground state properties experimental mea-
surements are naturally not restricted to temperature T = 0.
Indeed there is an interesting body of information which is
emerging in these Fermi gases about the finite temperature
behavior [4, 5, 6] of the breathing modes and, more recently
about the propagating sound velocity [7].
The purpose of this paper is to compute sound mode fre-
quencies in spherically trapped Fermi gases undergoing BCS
to Bose-Einstein condensation (BEC) crossover, at general T .
We present a solution of the linearized two fluid equations
and compare with recent experiments. We focus on the radial
breathing modes and present predictions for second sound, as
well. The structure of the two fluid equations for Bose [11]
and Fermi gases [12] has been rather extensively discussed.
In the crossover regime, the normal fluid is novel [13, 14],
containing both fermions and noncondensed pairs, which have
not been systematically included in previous collective mode
literature.
Of great importance to the field as a whole is the future
possibility of second sound observations, in large part because
this may ultimately help assign more precise experimental val-
ues to the transition temperature Tc. While existing experi-
mental techniques such as vortex observation [15] and den-
sity profile features [16] help establish superfluidity they pro-
vide lower bounds on Tc or determine its value for the special
case of a population imbalanced system. Thermodynamical
experiments measure Tc more directly [17, 18] but have been
confined to unitarity. Thus other techniques, such as second
sound observation will be of great value. One of the princi-
pal results of the present paper is an analysis of how second
sound is coupled to the breathing modes. We demonstrate that
higher order breathing modes will reveal Tc through this cou-
pling, and therefore are an alternative to direct second sound
measurements. However, the lowest breathing mode appears
to be remarkably T independent at unitarity. This has been ob-
served experimentally [4, 5] and argued to follow from isen-
tropic considerations [6]. Here we show, that even when we
treat the full coupling between first and second sound, we ob-
tain similar T independent behavior at unitarity.
At the core of the two fluid theory is the assumption that
hydrodynamics is valid and that there are frequent collisions
which produce a state of local thermodynamical equilibrium.
Although there are some exceptions [19], reaching the two
fluid regime has not been easy for atomic Bose gases. Two
fluid dynamics are more readily achieved for Fermi gases,
principally because in the crossover regime the large scatter-
ing lengths produce sufficient collisions. Nevertheless there
has been considerable theoretical interest in setting up [12]
and solving [20] the two fluid equations for Bose conden-
sates. Indeed, hydrodynamical approaches have successfully
addressed both the T = 0 and normal state regimes of the
Bose gases [21] and demonstrated that the breathing mode
frequencies are the same. Here, by contrast we address the
Fermi gas case in a trap. Because they interact more strongly
near unitarity, hydrodynamical descriptions have been argued
quite convincingly [6, 9, 22] to be valid.
Previous theoretical work has been confined to T = 0 treat-
ments of a harmonic trap, or to finite T theories [22] of a
uniform gas. Our work is most similar in spirit to an earlier
Bose gas study [20] although we introduce different numeri-
cal techniques as well as address fermions rather than bosons.
We note that the input thermodynamics of systems undergo-
ing the BCS-BEC crossover which is used in the present pa-
per, has been rather well calibrated against experimental mea-
surements in Ref. [17] and is based on a finite temperature
extension of the simplest (BCS-Leggett) ground state. In the
absence of a trap our results are for the most part similar to
those in Ref. [22].
We begin with the two fluid equations which describe the
dynamical coupling of the superfluid velocity vs and the nor-
mal velocity vn. Just as in the spirit of the original Landau
two fluid equations we ignore viscosity terms. In the presence
of a trap potential Vext =
mωhor
2, the two fluid equations
are given by m∂vs
+∇(µ+Vext+m
) = 0, ∂j
+∇·Π =
−n∇Vext,
+ ∇ · j = 0 and ∂(ns)
+ ∇ · (nsvn) = 0,
with Πij = pδij + nsv
s + nnv
n, n = ns + nn and
http://arxiv.org/abs/0704.1889v2
j = nsvs + nnvn. Here µ is the chemical potential, p the
pressure, and s the entropy per particle. Moreover, we have
nnvn + nsvs = nv. Here ns (vs) and nn (vn) denote the
superfluid and normal densities (velocities), respectively. We
use the subscript “0” to denote equilibrium quantities such that
vs0 = vn0 = 0, ∇(µ + Vext) = 0, ∇p0 = −n0∇Vext, and
n0, s0, µ0, p0 are independent of time t. Combining this with
the thermodynamic relation dµ = −s dT + dp/n, we have
∇T0 = 0, implying that temperature T0 is constant in the trap.
It then follows that in equilibrium µ = µ0 − Vext, consistent
with the Thomas-Fermi approximation.
For small deviations from equilibrium, we may linearize the
two fluid equations. Eliminating the velocities vs and vn, one
finds [12, 20]
= ∇ ·
−∇ · (δT n0∇s0) (1)
ns0n0s
− (∇s0)
2δT +∇s0 · ∇
We will focus on δµ(r) and δT (r) as the principal variables.
This choice, which is different from that in Ref. [20], is made
because both variables are non vanishing at the trap edge so
that in a basis set expansion they will satisfy the same bound-
ary conditions. Moreover, the two fluid equations are simplest
in this form. Expressing δs, δp and δn in terms of δµ, δT , the
two fluid equations can be rewritten as
, (3)
, (4)
A = ∇ · (n∇δµ) +∇ · (ns∇δT ) ,
s2∇δT
+∇s · ∇δµ+ s∇s · ∇δT .
It is understood that all coefficients of δµ and δT are calcu-
lated in equilibrium so that we drop the subscript “0”. The
thermodynamical quantities in equilibrium can be calculated
following Ref. [23], based on the standard local density ap-
proximation (LDA), µ(r) = µ0 − Vext(r). Their derivatives
with respect to T and µ can be calculated analytically, and
their gradients can be obtained via ∇f = −
∇Vext,
where f denotes any of the variables (n, ns, nn, s).
To solve the two coupled differential equations (3) and (4),
we assume a simple harmonic time dependence δµ, δT ∝
e−iωt. We cast the differential two fluid equations into an
eigenfunction problem with ω2 playing the role of eigenvalue
and the eigenfunctions given by the amplitudes of δµ and δT .
Since neither T nor µ depends on the density, they will not
vanish at the trap edge. Our boundary conditions require that
all thermodynamic variables be smooth (but not necessarily
zero) at the trap edge. At finite T , the density in the trap de-
creases exponentially when the local chemical potential be-
comes negative at large radius. We choose, thus, to expand
δµ and δT in terms of Jacobi polynomials. For our numerics
we choose the matrix dimension to be 300; we have similarly
investigated matrices of dimension 200 up to 900, and found
little change in our principal findings.
We now turn to an important aspect of our numerics. Be-
cause we generate some 300 frequencies in our numerical ap-
proach it is essential to establish a mechanism for systemati-
cally identifying first and second sound modes. To help find
such a “fingerprint”, we introduce a “decoupling approxima-
tion” based on reducing Eqs. (3) and (4) to δ̈µ = g1(δµ, δT ),
δ̈T = g2(δµ, δT ), where g1,2 are known functions. We elim-
inate cross terms by setting δT = 0 in g1 and δµ = 0 in
g2. With these two decoupled equations it is then relatively
straightforward to associate a profile plot of the numerically
calculated δp(r) and nδs(r) vs r within a trap with first or sec-
ond sound-like modes. Here we convert δµ and δT to δp and
δs via δp = nδµ + nsδT and δs =
We stress that this procedure differs somewhat from that in
Ref. [20] where it was argued that one could associate second
sound with a mode in which the oscillations in the thermody-
namical variable δT (r) were constrained to be in the conden-
sate region of the trap. As explained below, we do not find
this to be the case. An important check on our procedure is
that we find that there is no sign of second sound above Tc.
Up to this point everything is general applying to both
Fermi and Bose superfluids. All that is needed is a micro-
scopic theory for thermodynamical variables. Here we use a
calculational framework we have developed for treating BCS-
BEC crossover in trapped Fermi gases [23, 24] which em-
phasizes the importance of pseudogap effects or finite mo-
mentum pairs. The local thermodynamical potential (density)
Ω = Ωf + Ωb is associated with a contribution from gapped
fermionic excitationsΩf as well as from non-condensed pairs,
called Ωb. We have
Ωf = −
[(ξk − Ek)− 2T ln(1 + e
−Ek/T )],
Ωb = −Z∆
2µpair +
T ln(1 − e−Ωq/T ) . (5)
Here µpair is the chemical potential of the pairs which is
zero below Tc and the pair dispersion Ωq, along with the (in-
verse) residue Z , can be derived from a microscopic T -matrix
theory, described elsewhere [13, 14]. Using Ω one then ar-
rives at thermodynamical properties such as the entropy den-
sity ns = −∂Ω
as well as self consistent equations for the
total excitation gap ∆, the contribution to ∆ from noncon-
densed pairs (called the pseudogap) and the number equa-
tions. These self-consistent (local) equations are simply given
0 0.1 0.2 0.3 0.4
0 0.2 0.4 0.6 0.8 1
r/Rmax
T = 0.1
n = 1
T = 0.2
n = 1
(a) (b)
(c) (d)
Figure 1: (Color online) Behavior of the first (upper row) and second
sound (lower row) modes for a spherical trap at unitarity within the
“decoupling approximation” (see text). The left column shows the
T -dependence of the frequencies, while the right shows correspond-
ing typical spatial oscillation profiles for δp [(black) solid lines] and
nδs [(red) dashed lines] at the lowest frequencies, which provide fin-
gerprints of first and second sound. For the second sound, nδs dom-
inates and changes sign within the condensate. Here Tc ≈ 0.27TF ,
and the arrows indicate the condensate edge.
by ∂Ω
= 0, ∂Ω
∂µpair
= 0 and n = −∂Ω
, subject to the total
number constraint N =
d3r n(r). When T < Tc, we can
use the gap equation and chain rule to eliminate the variable
∆, which is a function of µ and T , via ∂∆
and ∂∆
. Similarly, when T > Tc we use
the gap and pseudogap equations to eliminate the variables ∆
and µpair to arrive at thermodynamical quantities.
Figure 1 shows the lowest two collective modes at unitar-
ity in a spherical trap, obtained by solving the chemical po-
tential or temperature fluctuation equation in the decoupling
approximation scheme described earlier. It is evident that our
approximated breathing mode frequency [Fig. 1(a)] is inde-
pendent of temperature. We understand this result by noting
that the decoupled equation for the breathing mode is given by
−ω2δµ = Cµµ1∇
2δµ+Cµµ2 · ∇δµ where Cµµ1 = n
and Cµµ2 = ∇µ. At unitarity, n
µ. The only
T -dependence contained in µ0 ≡ µ(r = 0) can be elimi-
nated via a simple re-scaling of r → r
2µ0/mω
ho, yielding
a T -independent breathing mode frequency. These arguments
can be shown to be equivalent to the isentropic assumption of
Ref. [6], where vs = vn is assumed in the two fluid model,
leading to two simpler coupled equations associated with the
Euler and the continuity equations.
By contrast, the second sound mode frequency we obtain
increases rapidly with temperature. Some typical oscillation
profiles of δp(r) and nδs(r) are shown in the right two pan-
els of Figure 1. Although the “entropy density oscillations”
nδs(r) fall off at large r, the entropy per particle δs(r) os-
cillations (not shown) increase very rapidly upon entering the
normal region. Consequently temperature fluctuations δT be-
0 0.2 0.4 0.6 0.8
0 0.2 0.4 0.6 0.8 1
T/TF = 0.01 T/TF = 0.2
1/kFa = 0
T/TF = 0.15
1/kFa = 0
T/TF = 0.08
1/kFa = 1
r/Rmax
1/kFa = 0
(a) (b)
(c) (d)
Figure 2: (Color online) Typical spatial oscillation profiles for δp
[(black) solid lines] and nδs [(red) dashed lines] obtained from the
fully coupled equations for a spherical trap at (a-c) unitarity and (d)
1/kF a = 1 at different T , for the first (top row) and second (bottom
row) sound modes. Also shown in (a) is the T = 0 analytical result
[(green) dot-dashed curve]. The arrows indicate condensate edge.
come large at the trap edge.
Our identification of first sound for the decoupled case leads
us to associate this mode in a coupled situation with a profile
for which at the trap center δp has a large amplitude, while
nδs is almost zero (with a small peak near the trap edge) as
in Fig. 1(b). By contrast, in the trap center the second sound
mode has large entropy fluctuations, while the pressure fluctu-
ations are almost zero [Fig. 1(d)]. These features will serve as
fingerprints for distinguishing (lower order) first and second
sound modes at finite temperatures.
Figure 2 shows some typical eigenfunction profiles of the
lowest modes obtained in the spherical trap upon solution of
the fully coupled two fluid equations. The first row corre-
sponds to the breathing mode in the unitary case. The good
agreement between the very low T result for δp and the T = 0
analytical solution (green dashed line) in Fig. 2(a) helps vali-
date our numerical scheme. Figure 2(b) corresponds to a high
temperature breathing mode. In this regime the pseudogap re-
gion outside the superfluid core is relatively large and the peak
in nδs(r) is accordingly very broad. The lower two panels
correspond to second sound modes for the unitary and BEC
cases. By contrast with the breathing modes, here nδs has a
larger amplitude than δp with an opposite sign. Clearly this
is very similar to what we observed in the decoupled mode
analysis of Fig. 1.
Figures 3(a-c) address the fully coupled equations and show
the behavior of the lowest breathing (upper branch) and sec-
ond sound mode (lower branch) frequencies as a function of
temperature in a spherical trap for 1/kFa = 1, 0, and −0.5,
respectively. For all three values of 1/kFa we find very little
sign of Tc in the lowest breathing mode frequency. In con-
trast, the second sound mode frequencies increase with T and
disappear above Tc. Figure 3(d) presents a more complete
series of modes for the unitary case. Here the lines serve as
guides to the eye for the lowest (blue) and higher order (green)
0 0.1 0.2 0.3
0 0.1 0.2 0.3 0.4
1/kFa = 1
1/kFa = 0
1/kFa = -0.5
1/kFa = 0
(a) (b)
(c) (d)
Figure 3: (Color online) Temperature dependence of breathing mode
and second sound frequencies in a spherical trap. Panels (a-c) are for
the near-BEC, unitary, and near-BCS cases, respectively. The upper
and lower branches in (a-c) represent the lowest breathing mode and
second sound frequencies, respectively. In (d), more complete results
(open circles) are shown at unitarity, where the lines serve as guides
to the eye for the breathing mode (nearly horizontal blue and green
lines) and second sound (increasing red curves) frequencies.
breathing mode, and second sound (red) frequencies.
At unitarity [Fig. 3(d)] one can identify a sequence of
higher order breathing modes which precisely overlap with
analytical calculations for T = 0. Importantly, only the low-
est of these is found to be a constant in temperature; the others
are found to mix with second sound modes, as indicated in the
solid and dashed (red) lines in Fig. 3(d). The behavior of the
lowest mode helps to justify the isentropic assumption made
in Ref. [6]. We understand this by referring back to the de-
coupled profiles in Figs. 1(b) and 1(c), which are seen to be
quite distinct. By contrast the profiles of the decoupled first
and second sound modes at higher order (not shown) appear
more similar to each other than their lowest order counter-
parts. Indeed, the behavior of the profiles at higher order is
associated with an increasing number of nodes in the curves
of δp and nδs (not shown), which leads to a greater similarity
between first and second sound profiles and helps explain why
the higher order modes are more readily mixed [25].
In summary, we have presented predictions for future ex-
periments on higher breathing modes and second sound in a
trap. We find that only the lowest breathing mode frequency
has very weak T dependence. For the unitary case this tem-
perature insensitivity was clearly observed by Thomas and co-
workers [6]. As a result of this experiment it should not be
surprising that we find relatively weak dependencies on ei-
ther side of the Feshbach resonance for this breathing mode
frequency. In the literature there are experimental claims (at
1/kFa = 1.0) which are consistent with a decrease [4] in the
radial breathing mode frequency, as indeed we find here, al-
though ours is probably too weak to reconcile the different
findings in Refs. [4] and [2, 3]. Finally, our more complete
studies at unitarity show that if higher order radial breathing
modes can be accessed, because of their strong hybridization
with second sound, it may be possible to use these breathing
modes, rather than direct second sound to determine the tran-
sition temperatures Tc. This should be of value since there are
currently few experiments which can assign a value to Tc over
the wider crossover regime.
This work was supported by Grant Nos. NSF PHY-
0555325 and NSF-MRSEC DMR-0213745; we thank John
Thomas for suggesting this problem and Cheng Chin for use-
ful conversations.
[1] C. Chin, M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim,
J. Hecker-Denschlag, and R. Grimm, Science 305, 1128 (2004).
[2] J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and J. E.
Thomas, Phys. Rev. Lett. 92, 150402 (2004).
[3] J. Kinast, A. Turlapov, and J. E. Thomas, Phys. Rev. A 70,
051401(R) (2004).
[4] A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R. Geursen,
M. Bartenstein, C. Chin, J. H. Denschlag, and R. Grimm, Phys.
Rev. Lett. 98, 040401 (2007).
[5] J. Kinast, A. Turlapov, and J. E. Thomas, Phys. Rev. Lett. 94,
170404 (2005).
[6] J. E. Thomas, J. Kinast, and A. Turlapov, Phys. Rev. Lett. 95,
120402 (2005).
[7] J. Joseph, B. Clancy, L. Luo, J. Kinast, A. Turlapov, and J. E.
Thomas (2006), cond-mat/0612567.
[8] H. Hu, A. Minguzzi, X.-J. Liu, and M. P. Tosi, Phys. Rev. Lett.
93, 190403 (2004).
[9] H. Heiselberg, Phys. Rev. Lett. 93, 040402 (2004).
[10] G. E. Astrakharchik, R. Combescot, X. Leyronas, and
S. Stringari, Phys. Rev. Lett. 95, 030404 (2005).
[11] A. Griffin and E. Zaremba, Phys. Rev. A 56, 4839 (1997).
[12] E. Taylor and A. Griffin, Phys. Rev. A 72, 053630 (2005).
[13] Q. J. Chen, I. Kosztin, B. Jankó, and K. Levin, Phys. Rev. Lett.
81, 4708 (1998).
[14] Q. J. Chen, J. Stajic, S. N. Tan, and K. Levin, Phys. Rep. 412, 1
(2005).
[15] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, and W. Ket-
terle, Nature 435, 170404 (2005).
[16] M. W. Zwierlein, A. Schirotzek, C. H. Schunck, and W. Ket-
terle, Science 311, 492 (2006).
[17] J. Kinast, A. Turlapov, J. E. Thomas, Q. J. Chen, J. Stajic, and
K. Levin, Science 307, 1296 (2005), published online 27 Jan-
uary 2005; doi:10.1126/science.1109220.
[18] L. Luo, B. Clancy, J. Joseph, J. Kinast, and J. E. Thomas, Phys.
Rev. Lett. 98, 080402 (2007).
[19] D. M. Stamper-Kurn, H.-J. Miesner, S. Inouye, M. R. Andrews,
and W. Ketterle, Phys. Rev. Lett. 81, 500 (1998).
[20] V. B. Shenoy and T.-L. Ho, Phys. Rev. Lett. 80, 3895 (1998).
[21] A. Griffin, W.-C. Wu, and S. Stringari, Phys. Rev. Lett. 78, 1838
(1997).
[22] H. Heiselberg, Phys. Rev. A 73, 013607 (2006).
[23] Q. J. Chen, J. Stajic, and K. Levin, Phys. Rev. Lett. 95, 260405
(2005).
[24] J. Stajic, Q. J. Chen, and K. Levin, Phys. Rev. Lett. 94, 060401
(2005).
[25] It is interesting to note that there exists a lower frequency first-
sound like mode at finite T , which we also find to be present in
the “free” (but hydrodynamic) Fermi gas limit at T 6= 0.
