First and second sound of a unitary Fermi gas in highly oblate harmonic traps
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] O c t First and second sound of a unitary Fermi gas in highly oblate harmonic traps
Hui Hu , Paul Dyke , Chris J. Vale , and Xia-Ji Liu ∗ Centre for Quantum and Optical Science, Swinburne University of Technology, Melbourne 3122, Australia (Dated: June 19, 2018)We theoretically investigate first and second sound modes of a unitary Fermi gas trapped ina highly oblate harmonic trap at finite temperatures. Following the idea by Stringari and co-workers [Phys. Rev. Lett. , 150402 (2010)], we argue that these modes can be described bythe simplified two-dimensional two-fluid hydrodynamic equations. Two possible schemes - soundwave propagation and breathing mode excitation - are considered. We calculate the sound wavevelocities and discretized sound mode frequencies, as a function of temperature. We find that in bothschemes, the coupling between first and second sound modes is large enough to induce significantdensity fluctuations, suggesting that second sound can be directly observed by measuring in-situ density profiles. The frequency of the second sound breathing mode is found to be highly sensitiveto the superfluid density.
PACS numbers: 67.85.Lm, 03.75.Ss, 05.30.Fk
Contents
I. Introduction II. 2D reduced thermodynamics
III. 2D simplified two-fluid hydrodynamicequations
IV. Free-propagating first and second sound V. Breathing first and second sound modes inharmonic traps
VI. Conclusions Acknowledgments A. Equation of state of a unitary Fermi gas B. Matrix elements of S ( ω ) C. Solving det[ ˜ S (˜ ω )] = 0 References ∗ Electronic address: [email protected]
I. INTRODUCTION
Low-energy excitations of a strongly interacting super-fluid quantum system - in which inter-particle collisionsare sufficiently strong to ensure local thermodynamicequilibrium - can be described by the Landau two-fluidhydrodynamic theory [1–4]. There are two well-knownkinds of excitations, referred to as first and second hydro-dynamic modes, which describe the coupled oscillationsof the superfluid and normal fluid components at finitetemperatures. First sound is an ordinary phenomenon inliquid, describing the propagation of a pressure or den-sity wave that is associated with normal acoustic sound.The motions of the superfluid and of normal fluid com-ponents are locked in phase. In contrast, second sound isa phenomenon characteristic of superfluids. It describesthe ability to propagate undamped entropy oscillations,which are essentially opposite phase oscillations of thetwo components. In the quantum liquid of superfluid he-lium, the velocities of first and second sound were firstmeasured in 1938 by Findlay et al . [5] and in 1946 byPeshkov [6], respectively.Ultracold atomic Fermi gases near a Feshbach reso-nance represent a new type of strongly interacting quan-tum system with unprecedented control over interatomicinteractions, dimensionality and purity [7]. These are an-ticipated to provide an ideal table-top system to deepenour understanding of fermionic superfluidity and Landautwo-fluid hydrodynamics. For a Fermi gas at unitarity,first sound modes have been investigated in great detail,both experimentally and theoretically [8–17], particularlynear zero temperature. In contrast, second sound modeis more difficult to address [18–27] and has only recentlybeen observed in experiments [28]. Key to understand-ing second sound, is knowledge of the superfluid den-sity [20, 29, 30], a quantity which is still a subject of in-tense investigation. Theoretically, first and second soundmodes of the Landau two-fluid hydrodynamic equationswere solved by Taylor et al . for an isotropically trappedunitary Fermi gas [21], using an assumed superfluid den-sity and by developing a variational approach. For gen-eral anisotropic harmonic traps, the Landau hydrody-namic equations are more difficult to solve. However,for a highly elongated configuration, it was shown byBertaina, Pitaevskii and Stringari [24] that, a simplifiedone-dimensional (1D) two-fluid hydrodynamic descrip-tion may be derived, generalizing the dimensional reduc-tion at zero temperature [15]. Following this pioneeringidea, the propagation of a second sound wave was strik-ingly observed in a highly elongated unitary Fermi gasby the the Innsbruck team [28]. Using the measured sec-ond sound velocity data, the simplified 1D hydrodynamicequations were solved and used to extract the superfluiddensity [28]. These results show that studies of secondsound provide a promising route towards an accurate de-termination of the superfluid density.Here, we would consider a unitary Fermi gas in highly oblate harmonic traps, which can be readily realised incurrent experiments. Indeed, such configurations havealready been investigated in laboratory experiments [31–36]. The analysis presented here could therefore be di-rectly tested in future experiments. Our main resultsare briefly summarized as follows. Following the ideasof Stringari and co-workers [24, 26], we derive the sim-plified two-dimensional (2D) Landau two-fluid hydrody-namic equations. Using a variational approach within thelocal density approximation [18, 20, 21], we fully solve thecoupled 2D hydrodynamic equations in the presence of aweak radial harmonic confinement. The discretized modefrequencies of first and second sound are calculated. Wefind that the density fluctuation due to the second soundmode is significant, revealing that second sound is di-rectly observable from in-situ density profiles, after anappropriate modulation of the weak radial trapping po-tential. The mode frequency of second sound is foundto depend very critically on the form of the superfluiddensity. Our quasi-2D analysis is similar to the quasi-1Dstudy of Stringari and co-workers [26], however, there area number of significant differences. Specifically we fullysolve the coupled hydrodynamic equations without as-suming that first and second sound are decoupled. Thisenables us to calculate the density fluctuations associ-ated with the discretized second sound modes. From anexperimental point of view, the possibility of measuringdiscrete breathing mode frequencies may allow a moreaccurate determination of superfluid density.The remainder of the paper is organized as follows.In the next section, we outline the reduced 2D univer-sal thermodynamics, satisfied by the unitary Fermi gasin highly oblate harmonic traps based on the measured3D equation of state [37]. In Sec. III, we derive thesimplified 2D Landau two-fluid hydrodynamic equationsand discuss very briefly their applicability. In Sec. IV, weconsider the propagation of first and second sound, whichmay be excited by creating a short perturbation of thedensity. The nature of first and second sound excitationsin the quasi-2D unitary Fermi gas is described. In Sec. V,we fully solve the simplified 2D hydrodynamic equations and calculate the discretized breathing mode frequenciesof first and second sound. The density fluctuation associ-ated with low-lying second sound modes is discussed andthe dependence of the second sound mode frequency onthe form of superfluid density is examined. Finally, inSec. V we draw our conclusions. The appendix presentsfurther details of our numerical calculations.
II. 2D REDUCED THERMODYNAMICS
Let us consider a unitary Fermi gas trapped in a highlyanisotropic pancake-like harmonic potential, V ext ( r ⊥ , z ) = 12 mω ⊥ r ⊥ + 12 mω z z , (1)with atomic mass m and the radial and axial trappingfrequencies ω ⊥ , ω z . The aspect ratio of the trap, definedby λ = ω z /ω ⊥ ≫ , can be specifically tuned in experi-ments. The number of atoms in the typical sound modeexperiments [9, 10, 13, 16, 28] is about N ∼ . Forsuch a large number of atoms, it is standard to use thelocal density approximation [7], which amounts to treat-ing the atoms in the position ( r ⊥ , z ) locally as a uniformmatter, with a local chemical potential given by, µ ( r ⊥ , z ) = µ − V ext ( r ⊥ , z ) . (2)Here µ is the chemical potential at the trap center. Thevalidity of the local density approximation can be con-veniently estimated by comparing the number of atoms N ∼ with a threshold N D = λ ( λ + 1) , belowwhich there is a dimensional crossover from three- to two-dimensions [33, 34]. At Swinburne, we have previouslydemonstrated pancake traps with λ ∼ . Thus, the ratio N/N D ∼ , indicating that the system could be deeplyin the 3D limit where the local density approximation iswell justified. A. Universal thermodynamics
A remarkable feature of a uniform unitary Fermi gasis that all its thermodynamic functions can be expressedas universal functions a single dimensionless parameter µ/ ( k B T ) [38, 39]. This is simply due to the fact that atunitarity, the s -wave scattering length - the only lengthscale used to characterize the short-range interatomic in-teraction - diverges. Thus, the remaining length scalesare the thermal wavelength λ T ≡ p π ~ / ( mk B T ) (3)and the mean distance between atoms n − / , where n isthe number density. Accordingly, the energy scales aregiven by the temperature, k B T , and by the Fermi energy, k B T F = ~ m (cid:0) π n (cid:1) / , (4) -3 -2 -1 0 1 2 3 41234 -2 0 2 40.00.51.0 MIT expt. smoothed fit g ( x ) = n () / n () x = /k B T dg ( x ) / d x x FIG. 1: (Color online) The universal scaling function ofa strongly interacting unitary Fermi gas, n ( µ, T ) /n ( µ, T ) ,where n ( µ, T ) is the number density of an ideal spin-1/2Fermi gas. The experimental data from the MIT team (redcircles) [37] are compared with a smoothed fit (the black line),as described in the text. The inset shows the comparison forthe derivative of the universal scaling function. The verti-cal grey lines indicate the critical threshold for superfluidity, ( µ/k B T ) c ≃ . [37]. or, alternatively, by the chemical potential µ . Follow-ing the scaling analysis, all the thermodynamic functionscan therefore be expressed in terms of universal func-tions that only depends on the dimensionless parameter x = µ/ ( k B T ) . For the pressure and number density, theseuniversal functions are given, respectively, by f Dp ( x ) = ( k B T ) − λ T P ( µ, T ) , (5) f Dn ( x ) = λ T n ( µ, T ) . (6)Using the thermodynamic relation n = ( ∂P/∂µ ) T , wefind f Dn ( x ) = df Dp ( x ) /dx . It is easy to see that know-ing f Dp and f Dn , we can calculate directly all the ther-modynamic functions of the uniform unitary Fermi gas[26].It is highly non-trivial to theoretically determine theuniversal function f Dp or f Dn , since there are no smallparameters to control the theory of a strongly interact-ing Fermi gas [40–44], except at high temperatures, wherethe virial expansion approach is applicable [45, 46]. For-tunately, accurate experimental data for the equation ofstate are now available from the landmark experimentsperformed by the Massachusetts Institute of Technology(MIT) team [37]. In Fig. 1, we present their data for theuniversal scaling function, g ( x ) ≡ n ( µ, T ) n ( µ, T ) = f Dn ( x ) f Dn, ( x ) , (7)where the subscript “0” indicates the result of an ideal, non-interacting Fermi gas and f Dn, ( x ) = 4 √ π ˆ ∞ dt √ te − t e x e x e − t . (8)For later convenience of numerical calculations, we havefitted the experimental data to an analytic expression, asdetailed in the Appendix A. The fitting curve is smoothlyextrapolated to both low and high temperature regimeswhere the behavior of f Dn ( x ) is well-understood [23, 25,45, 46]. As evident in Fig. 1, the relative error in thefitting is less than . , significantly smaller than thestandard error for the experimental data (i.e., ). B. Local density approximation
Within the local density approximation, we can writethe local pressure and number density using the universalfunctions, P ( r ⊥ , z ) = k B Tλ T f Dp (cid:20) µ ( r ⊥ , z ) k B T (cid:21) , (9) n ( r ⊥ , z ) = 1 λ T f Dn (cid:20) µ ( r ⊥ , z ) k B T (cid:21) , (10)where the local chemical potential is given by Eq. (2).For a trapped Fermi gas in three dimensions, the Fermitemperature T F is defined as [7] k B T F = ~ (cid:0) N ω ⊥ ω z (cid:1) / . (11)Using the number equation N = ´ d r ⊥ dzn ( r ⊥ , z ) , werelate µ /k B T to the reduced temperature T /T F by, TT F = (cid:20) √ π ˆ ∞ dt √ tf Dn ( µ k B T − t ) (cid:21) − / . (12)The onset of superfluidity in a trapped unitary Fermigas occurs at µ /k B T = x c ≃ . . Numerically, wefind that T c ≃ . T F , in agreement with the previousresults [26, 37]. C. 2D reduced thermodynamic functions
In highly oblate harmonic traps, although the cloudis still three-dimensional, its low-energy dynamics aregreatly affected by the tight axial confinement. Thissituation is very similar to a highly elongated unitaryFermi gas considered earlier by Bertaina, Pitaevskii andStringari [24], who showed that with tight radial con-finement, the standard Landau two-fluid hydrodynamicequations could reduce to a simplified 1D form. For thesame reason, which will be explained in greater detailin the next section, first and second sound under thetight axial confinement can be described by a simpli-fied 2D two-fluid hydrodynamic description. In brief, dueto the nonzero viscosity and thermal conductivity, localfluctuations in temperature ( δT ) and chemical potential( δµ ) become independent of the axial coordinate, for anylow-energy excitations at the frequency ω ∼ ω ⊥ ≪ ω z .Therefore, the axial degree of freedom in all the thermo-dynamic variables that enter the Landau two-fluid hy-drodynamic equations becomes irrelevant and can be in-tegrated out. We can then derive 2D reduced thermody-namics and immediately have the reduced Gibbs-Duhemrelation, δP = s δT + n δµ, (13)where the variables P , s and n are the axial integralsof their three-dimensional counterparts, namely the lo-cal pressure, entropy density and number density. Thevariable P is given by, P ( r ⊥ ) ≡ ˆ dzP ( r ⊥ , z ) = m ( k B T ) π / ~ ω z f p ( x ) , (14)where x ≡ (cid:0) µ − mω ⊥ r ⊥ / (cid:1) k B T = µ ( r ⊥ ) k B T (15)and we have introduced the universal scaling function, f p ( x ) ≡ ˆ ∞ dtf Dp (cid:0) x − t (cid:1) . (16)All the 2D thermodynamic variables can be derived fromthe reduced Gibbs-Duhem relation, for example, n = (cid:18) ∂P ∂µ (cid:19) T = m ( k B T ) π / ~ ω z f n ( x ) , (17) s = (cid:18) ∂P ∂T (cid:19) µ = m ( k B T ) π / ~ ω z [3 f p ( x ) − xf n ( x )] , (18)where f n ( x ) ≡ df p ( x ) dx = ˆ ∞ dtf Dn (cid:0) x − t (cid:1) . (19)In addition, it is readily shown that the specific heatper particle at constant column density and pressure aregiven by [25], ¯ c v = T (cid:18) ∂ ¯ s ∂T (cid:19) n = 6 f p ( x ) f n ( x ) − f n ( x ) f ′ n ( x ) , (20) ¯ c p = T (cid:18) ∂ ¯ s ∂T (cid:19) P = ¯ c v " f p ( x ) f ′ n ( x ) f n ( x ) , (21)where ¯ s ≡ s / ( n k B ) is the entropy per particle and f ′ n ( x ) ≡ df n ( x ) /dx . It is also straightforward to checkthe universal relations, (cid:18) ∂P ∂n (cid:19) ¯ s = 3 P n , (22) (cid:18) ∂P ∂s (cid:19) n = T . (23) c p /nk B c v /nk B (b) e n t r opy , s p ec i f i c h ea t s x = /k B Ts/nk B f p (x) f n (x) df n (x)/dx f ( x ) (a) FIG. 2: (Color online) (a) 2D universal scaling functions f p ( x ) , f n ( x ) and df n ( x ) /dx as a function of the dimension-less variable x = µ/ ( k B T ) . (b) 2D entropy ¯ s = s/ ( nk B ) andspecific heats per particle ¯ c v = c v / ( nk B ) and ¯ c p = c p / ( nk B ) ,as a function of x = µ/ ( k B T ) . The vertical grey lines indicatethe critical threshold for superfluidity, x c ≃ . [37]. In Fig. 2, we report the relevant 2D universal scalingfunctions, calculated by using the experimental MIT datafor f Dn ( x ) after smoothing. D. Superfluid density
We can also express the local superfluid density as auniversal (but as yet undetermined) function f Ds : n s ( r ⊥ , z ) = 1 λ T f Ds (cid:20) µ ( r ⊥ , z ) k B T (cid:21) . (24)Integrating out the axial coordinate, we obtain n s ( r ⊥ ) = ˆ dzn s ( r ⊥ , z ) = m ( k B T ) π / ~ ω z f s ( x ) , (25)where the universal scaling function f s ( x ) is given by, f s ( x ) = ˆ ∞ dtf Ds (cid:0) x − t (cid:1) . (26)As mentioned earlier, in contrast to the equation ofstate, the universal function for the superfluid density ofa unitary Fermi gas f Ds ( x ) is not yet known precisely n s / n T/T c x= /k B T F f s (x) FIG. 3: (Color online) Bulk superfluid fraction of superfluidhelium (red circles) [47] and the corresponding 2D reducedsuperfluid fraction (the black line). The inset shows the 2Duniversal scaling function for superfluid density, to be usedfor a unitary Fermi gas. The vertical grey lines indicate thecritical threshold for superfluidity, x c ≃ . [37]. [20, 29, 30]. For illustrative purposes in the present work,we will consider the superfluid fraction of superfluid he-lium [47], (cid:16) n s n (cid:17) D = f He (cid:18) TT c (cid:19) , (27)which is plotted in Fig. 3 by red circles. Our choice is mo-tivated by the similarity for hydrodynamics between su-perfluid helium and unitary Fermi gas, presumably aris-ing from strong correlations in both systems, and justi-fied in part by measurements in Ref. [28]. Using Eq.(6), we find that T /T c = [ f Dn ( x ) /f Dn ( x c ≃ . − / .Thus, from Eq. (27), the universal function f Ds ( x ) canbe calculated using f Ds ( x ) = f Dn ( x ) f He (cid:20) f Dn ( x ) f Dn ( x c ≃ . (cid:21) − / ! (28)and f s ( x ) is then obtained. In the inset of Fig. 3, we show f s ( x ) calculated with the superfluid fraction of superfluidhelium.In a 2D configuration with ω ⊥ = 0 , it is natural to de-fine a Fermi temperature in terms of the column density n : k B T DF = ( ~ ω z ) / (cid:18) π ~ n m (cid:19) / , (29)which coincides with the usual three-dimensional defini-tion of the Fermi temperature, Eq. (4) . The reducedtemperature is then given by, TT DF = (cid:20) √ π f n ( x ) (cid:21) − / . (30) Therefore, the critical temperature of the unitary Fermigas in the pancake geometry is given by, T Dc =[2 f n ( x c ) / √ π ] − / T DF ≃ . T DF . In Fig. 3, we showthe 2D superfluid fraction f s /f n (the black line) as afunction of T /T Dc . It lies systematically below the three-dimensional superfluid fraction (red circles), due to theintegration over the axial degree of freedom.It should be noted that the critical behavior of thesuperfluid density near the phase transition is greatly af-fected by the axial integration. In three dimensions, wemay define the critical exponent α by n s ( T → T c ) ∼ ( T − T c ) α or f Ds ( x → x c ) ∼ ( x − x c ) α . After integratingout the axial degree of freedom, it is easy to show, f s ( x → x c ) ∼ ( x − x c ) α +1 / . (31)As we shall see, the increase in the critical exponent willreduce the second sound velocity and hence the secondsound mode frequency. III. 2D SIMPLIFIED TWO-FLUIDHYDRODYNAMIC EQUATIONS
We now derive the simplified 2D Landau two-fluid hy-drodynamic equations, following the idea by Bertaina,Pitaevskii and Stringari [24]. This 2D picture will gen-erally be valid as we are considering excitations whosewavelength is long compared to the transverse (axial)cloud width.The standard two-fluid equations in the harmonic trap V ext are given by, m∂ t n + ∇ · j = 0 , (32) m∂ t v s + ∇ ( µ + V ext ) = 0 , (33) ∂ t j i + ∂ i P + n∂ i V ext = ∂ t ( η Γ ik ) , (34) ∂ t s + ∇ · ( s v n ) = ∇ · ( κ ∇ T /T ) , (35)where j = m ( n s v s + n n v n ) is the current density, n s and n n are the superfluid and normal density, v s and v n are the corresponding velocity fields, Γ ik ≡ ( ∂ k v ni + ∂ i v nk − δ ik ∂ j v nj / , and finally η and κ are the shearviscosity and thermal conductivity, respectively. In theabove equations, we have kept only linear terms in the ve-locity, as we are interested in small-amplitude dynamicsin the linear response regime. Moreover, we have omit-ted bulk viscosity terms which give smaller contributions.The Landau two-fluid hydrodynamic theory is applicablewhen the mean free path l is much smaller than the wave-length of the sound wave λ . By considering the axial sizeof the Fermi cloud R z ∼ p ~ / ( mω z ) , we further require R z ≪ λ , therefore the hydrodynamic regime in a har-monic trap is achieved when l ≪ R z ≪ λ .As discussed by Stringari and co-workers [24], in thepresence of tight confinement, the nonzero viscosity andthermal conductivity may significantly change the low-energy dynamics at the frequency ω ∼ ω ⊥ ≪ ω z . Thisoccurs when the viscous penetration depth δ - the typicallength scale at which an excitation becomes damped dueto shear viscosity - fulfils the condition, δ ≡ r ηmn n ω ≫ R z . (36)In the case of superfluid helium in a thin capillary, theabove condition makes the normal component of the liq-uid stick to the wall and thus the normal velocity fieldvanishes [48]. With a soft wall caused by the tight axialconfinement, the normal component can move but thenormal velocity field becomes uniform along the axial di-rection [24]. For the same reason, the superfluid velocityalso becomes independent of the axial coordinate. Anal-ogously, a nonzero thermal conductivity in Eq. (35) forthe entropy conservation leads to the independence of thetemperature fluctuation δT on the axial coordinate. Thetight axial confinement also implies that the axial compo-nent of the velocity fields, for both superfluid and normalfluid, must be much smaller than the transverse one andmay be neglected. Using Eq. (33) for the superfluid ve-locity implies that the chemical potential fluctuation δµ isessentially independent of the axial coordinate. In brief,we conclude that, owing to the crucial roles played by theshear viscosity and thermal conductivity, under the tightaxial confinement the velocity fields are independent ofthe axial position and a global thermal equilibrium alongthe axial direction is established.With these observations, it is straightforward to writedown the simplified two-fluid hydrodynamic equations byintegrating out the axial degree of freedom in Eqs. (32),(33), (34) and (35), which take the following forms, m∂ t n + ∇ ⊥ · j ⊥ = 0 , (37) m∂ t v s ⊥ + ∇ ⊥ [ µ ( r ⊥ ) + V ext ( r ⊥ )] = 0 , (38) ∂ t j ⊥ + ∇ ⊥ P + n ∇ ⊥ V ext = 0 , (39) ∂ t s + ∇ ⊥ · ( s v n ⊥ ) = 0 , (40)where the current density now becomes j ⊥ = m ( n s v s ⊥ + n n v n ⊥ ) and we have omitted the residual dissipationterms along the weakly confined direction, as in a uniformfluid the effect of viscosity and thermal conductivity isirrelevant in the long-wavelength limit.We note that, by introducing a characteristic colli-sional time τ related to viscosity, τ ≃ η/ ( mn ¯ v ) ∼ η/ ( n ~ ω z ) , where ¯ v is the average velocity of atoms andis of the order of the Fermi velocity v F ∼ p ~ ω z /m , it iseasy to check that Eq. (36) is equivalent to requiring thelow frequency condition, ω ≪ ω z τ. (41)Recalling that ω ∼ ω ⊥ and ω z = λω ⊥ , the above re-quirement for achieving the reduced 2D hydrodynamicscan be rewritten as λ ωτ ≫ . For a highly oblate trapwith typical aspect ratio λ ∼ ≫ , this conditionis compatible with the condition to enter the hydrody-namic regime, l ≪ λ , which can alternatively be writtenas ωτ ≪ . A. Variational reformulation of the simplifiedLandau hydrodynamic equations
A convenient way to solve the simplified two-fluidhydrodynamic equations is to reformulate them usingHamilton’s variational principle [49]. Following previ-ous work [18, 20], the hydrodynamic modes of theseequations with frequency ω at temperature T can beobtained by minimizing a variational action, which, interms of displacement fields u s ( r ⊥ ) = iωv s ⊥ ( r ⊥ ) and u n ( r ⊥ ) = iωv n ⊥ ( r ⊥ ) [50], takes the following form, S = 12 ˆ d r ⊥ " mω (cid:0) n s u s + n n u n (cid:1) − (cid:18) ∂µ∂n (cid:19) s ( δn ) − (cid:18) ∂T∂n (cid:19) s δnδs − (cid:18) ∂T∂s (cid:19) n ( δs ) . (42)Here, n s ( r ⊥ ) and n n ( r ⊥ ) = n ( r ⊥ ) − n s ( r ⊥ ) are the2D reduced superfluid and normal-fluid densities at equi-librium, as discussed in the previous section Sec. II.The density fluctuation δn and entropy fluctuation δs are given by, δn ( r ⊥ ) ≡ −∇ ⊥ · ( n s u s + n n u n ) (43)and δs ( r ⊥ ) ≡ −∇ ⊥ · ( s u n ) , (44)respectively. The effect of the weak radial trapping po-tential V ext ( r ⊥ ) = mω ⊥ r ⊥ / enters the action Eq. (42)through the coordinate dependence of the equilibriumthermodynamic variables ( ∂µ/∂n ) s , ( ∂T /∂n ) s and ( ∂T /∂s ) n , within the local density approximation. Westress that, all these thermodynamic variables can be ob-tained by using the reduced Gibbs-Duhem relation Eq.(13) and can be expressed by the universal functions f p ( x ) and f n ( x ) .In superfluid helium, the solutions of the Landau two-fluid hydrodynamic equations can be well understood asdensity and entropy (temperature) waves, which are thepure in-phase mode with u s = u n and the pure out-of-phase mode with n s u s + n n u n = 0 , known as first andsecond sound, respectively [2–4]. For a strongly inter-acting unitary Fermi gas, we could use the same classifi-cation [21]. For this purpose, we rewrite the action Eq.(42) in terms of two new displacement fields u a = ( n s u s + n n u n ) /n (45)and u e = u s − u n , (46)considering that the density and temperature fluctua-tions are given by δn = −∇ ⊥ · ( n u a ) (47)and δT = (cid:18) ∂T∂s (cid:19) n ∇ ⊥ · (cid:18) s n s n u e (cid:19) , (48)respectively. Ideally, first sound is characterized by δn =0 but δT = 0 and second sound by δn = 0 but δT = 0 .Using the standard thermodynamic identities derived from the reduced Gibbs-Duhem relation Eq. (13), aftersome straightforward but lengthy algebra, we arrive at S = 12 ˆ d r ⊥ h S ( a ) + 2 S ( ae ) + S ( e ) i , (49)where S ( a ) = mω n u a + ( ∇ ⊥ n · u a ) ( ∇ ⊥ V ext · u a ) + 2 ( n ∇ ⊥ V ext · u a ) ( ∇ ⊥ · u a ) − n (cid:18) ∂P ∂n (cid:19) ¯ s ( ∇ ⊥ · u a ) , (50) S ( ae ) = (cid:18) ∂P ∂s (cid:19) n ( ∇ ⊥ · u a ) (cid:20) ∇ ⊥ · (cid:18) s n s n u e (cid:19)(cid:21) , (51) S ( e ) = mω n s n n n u e − (cid:18) ∂T∂s (cid:19) n (cid:20) ∇ ⊥ · (cid:18) s n s n u e (cid:19)(cid:21) . (52)It is clear that the first and second sound are governedby the actions S ( a ) and S ( e ) , respectively. The couplingbetween first and second sound is controlled by the cou-pling term S ( ae ) , which is in general nonzero. Indeed,in our case, as ( ∂P /∂s ) n = T / , strictly speaking thefirst and second sound are coupled at any finite temper-atures. IV. FREE-PROPAGATING FIRST ANDSECOND SOUND
For a uniform superfluid ( V ext = 0 ), the solutions of S ( a ) and S ( e ) are plane waves of wave vector q with dis-persion ω = c q and ω = c q , where c = s m (cid:18) ∂P ∂n (cid:19) ¯ s (53)and c = s k B Tm ¯ s ¯ c v n s n n . (54)These expressions for first and second sound velocities arethe standard results used to describe superfluid helium,when the corresponding equation of state and superfluiddensity are used [2–4]. In Fig. 4, we show the decoupledfirst and second sound velocities of a highly oblate Fermigas at unitarity, by using dashed lines.Including the coupling term S ( ae ) and using the stan-dard thermodynamic relations, it is straightforward toshow that the solutions for sound velocities u of the sim-plified two-fluid hydrodynamic equations satisfy the fol-lowing equation, u − u (cid:0) c + c (cid:1) + c c γ = 0 , (55) Second sound c / v F D T/T
F2D
First sound
T/T
F2D
LP Ratio
FIG. 4: (Color online) 2D first and second sound velocities(solid lines) as a function of temperature, in units of the Fermivelocity v DF = p k B T DF /m of an ideal Fermi gas at zerotemperature at the trap center. The dashed lines give the de-coupled sound velocities. The inset shows the Landau-Placzekparameter ǫ LP . The vertical grey lines indicate the criticaltemperature for superfluidity, T c ≃ . T DF , in the absenceof harmonic confinement in the transverse direction. where γ ≡ ( ∂P /∂n ) ¯ s ( ∂P /∂n ) T = ¯ c p ¯ c v > . (56)It gives rise to two solutions for the sound velocity, u and u , which in the absence of the coupling term S ( ae ) (i.e., γ = 1 ), coincide with the decoupled first and sec-ond sound velocities, c and c . The numerical resultsfor u and u are shown in Fig. 4 by solid lines. Thetemperature dependence of the reduced sound velocitiesin Fig. 4 is very similar to that of a 3D unitary Fermi gaspredicted in the earlier works [23, 29]. We attribute thisqualitative similarity to the strongly interacting natureof the system.In the case of a small parameter θ ≡ c / ( γc ) ≪ ,which is indeed true for a highly oblate unitary Fermigas, we may solve Eq. (55) perturbatively. We find theexpansions [23], u = c [1 + ( γ − θ + · · · ] , (57) u = c γ [1 − ( γ − θ + · · · ] . (58)To the leading order of θ , the first sound velocity is notaffected by the coupling term, but the second sound ve-locity decreases by a factor of √ γ and is now given by[23, 26], u = s k B Tm ¯ s ¯ c p n s n n . (59)Therefore, quantitatively, the coupling between first andsecond sound can be characterized by the so-calledLandau-Placzek (LP) parameter ǫ LP ≡ γ − [23]. In theinset of Fig. 4, we plot the LP parameter as a function oftemperature. In the superfluid phase, it is always smallerthan . , indicating that the first and second sound cou-ple very weakly in a unitary Fermi gas. A similar situ-ation happens in superfluid helium, where ¯ c p ≃ ¯ c v andhence ǫ LP ≃ [20]. From the above approximate expres-sion for the second sound velocity and Eq. (31) for 2Dsuperfluid density, it is clear the velocity vanishes as u ∼ (cid:18) − TT c (cid:19) α/ / , (60)when approaching to the superfluid phase transition frombelow. Here, α ≃ / is the critical exponent of the su-perfluid density of a unitary Fermi gas in three dimen-sions.Although the coupling between first and second soundis weak in a unitary Fermi gas, it is of significant impor-tance for the purpose of experimental detection of sec-ond sound. Unlike superfluid helium, temperature oscil-lations - the characteristic motion of the second sound -are difficult to observe in ultracold atomic gases. Den-sity measurement is the most efficient way to characterizeany low-energy dynamics. As a result, an observable sec-ond sound mode must have a sizable density fluctuation,which can only be induced by its coupling to first soundmodes. The strength of density fluctuations in secondsound may be conveniently characterized by the follow-ing ratio between the relative density and temperaturefluctuations, (cid:18) δn/n δT /T (cid:19) nd ≃ Tn (cid:18) ∂n ∂T (cid:19) P = 2 − f p f ′ n f n , (61)where the subscript “2nd” indicates the second sound andwe have taken the derivative at const pressure, as the | n | / n T/T
F2D
FIG. 5: (Color online) The density fluctuation of 2D secondsound, calculated in the case of δT /T = 10% . The verticalgrey lines indicate the critical temperature for superfluidity, T c ≃ . T DF , in the absence of harmonic confinement inthe transverse direction. second sound is properly considered as an oscillating waveat constant pressure, rather than at constant density, asindicated by Eq. (59) [26]. In contrast, the ratio betweenthe relative density and temperature fluctuations in firstsound is given by, (cid:18) δn/n δT /T (cid:19) st ≃ Tn (cid:18) ∂n ∂T (cid:19) ¯ s = 2 . (62)In Fig. 5, we report the density fluctuation of sec-ond sound as a function of temperature in the superfluidphase, calculated with an assumed temperature fluctu-ation ratio δT /T = 10% . This is a typical tempera-ture fluctuation, achievable in current experiments. Forexample, in the recent first-sound collective mode mea-surement, it was shown that ( δn/n ) st ∼ over awide temperature window [16, 17]. By using Eq. (62),we therefore assume that a temperature fluctuation at δT /T = 10% can be easily excited. From Fig. 5, it iseasy to see that the density fluctuation of second sound isvery significant when temperature T > . T c ∼ . T DF ,revealing that the propagation of a second sound pulsecould be experimentally detected via density measure-ments. This result suggests that a similar approachto that demonstrated by the Innsbruck experiment [28]could be applied in the oblate geometry considered here,where a focused blue-detuned laser propagating along thetightly confined direction could excite a second soundwave which will propagate radially outwards from thecloud centre. Subsequent absorption images taken atdifferent times after the laser pulse could record the den-sity fluctuation induced by the propagating second soundwave.While the essential physics of our scheme is the same asthat considered by by Stringari and co-workers in highlyelongated trap [26], the quasi-2D geometry may offer ad-vantages for experiments. Specifically, clouds confinedin elongated traps typically have very high peak densi-ties leading to high optical densities that are difficult tomeasure quantitatively using absorption imaging [51, 52].This means quantifying small density fluctuations in suchclouds can prove challenging. In a quasi-2D trap however,such high peak optical densities can be avoided and im-ages taken along the direction of tight confinement canbe integrated over the angular coordinate to obtain anaccurate measure of the radial density. Small densityfluctuations should therefore be more easily detected. V. BREATHING FIRST AND SECOND SOUNDMODES IN HARMONIC TRAPS
We now consider a weak transverse harmonic trap, V ext ( r ⊥ ) = mω ⊥ r ⊥ / and fully solve the simplified Lan-dau two-fluid hydrodynamic equations using a variationalapproach, developed in the previous work [20, 21]. Wefocus on compressional (breathing) modes with the pro-jected angular momentum l z = 0 , since these modes arethe easiest one to excite experimentally. A. Variational approach
For breathing modes, we assume the following polyno-mial ansatz for the displacement fields: u a ( r ⊥ ) = ˆr ⊥ N p − X i =0 A i ˜ r i +1 ⊥ , (63) u e ( r ⊥ ) = ˆr ⊥ N p − X i =0 B i ˜ r i +1 ⊥ , (64)where ˆr ⊥ is the unit vector in the radial direction, { A i , B i } ( i = 0 , · · · , N p − ) are the N p variational pa-rameters, and ˜ r ⊥ ≡ r ⊥ /R F is the dimensionless radialcoordinate with R F being the Thomas-Fermi radius. Byinserting this variational ansatz into the action Eq. (49),we express the action S as functions of the N p varia-tional parameters { A i , B i } . The mode frequencies areobtained by minimizing the action S with respect to these N p parameters. The accuracy of our variational calcu-lations can be improved by increasing the value of N p .It turns out that the calculations converge quickly as N p increases.In more detail, it is straightforward to show that, theexpression of the action can be written in a compact form, S = 12 Λ † S ( ω ) Λ , (65)where Λ ≡ (cid:2) A , B , · · · , A i , B i , · · · , A N p − , B N p − (cid:3) T and S ( ω ) is a N p × N p matrix with block elements ( i, j = 0 , · · · , N p − ), [ S ( ω )] ij ≡ " M ( a ) ij ω − K ( a ) ij − K ( ae ) ij − K ( ae ) ji M ( e ) ij ω − K ( e ) ij . (66)In [ S ( ω )] ij , we have introduced the weighted mass mo-ments, M ( a ) ij = m ˆ d r ⊥ ˜ r i + j +2 ⊥ n ( r ⊥ ) , (67) M ( e ) ij = m ˆ d r ⊥ ˜ r i + j +2 ⊥ (cid:20) n s n n n (cid:21) ( r ⊥ ) , (68)and the spring constants, K ( a ) ij = 32 ( i + 2) ( j + 2) R F ˆ d r ⊥ ˜ r i + j ⊥ P ( r ⊥ ) , (69) K ( ae ) ij = T i ( i + 2) R F ˆ d r ⊥ ˜ r i + j ⊥ (cid:20) s n s n (cid:21) ( r ⊥ ) , (70) K ( e ) ij = ˆ d r ⊥ (cid:18) ∂T∂s (cid:19) n r ⊥ d (cid:2) s n s ˜ r i +2 ⊥ /n (cid:3) d ˜ r ⊥ × d h s n s ˜ r j +2 ⊥ /n i d ˜ r ⊥ . (71)To derive the above equations, we have used the univer-sal relations satisfied by the highly oblate unitary Fermigas: n ( ∂P /∂n ) ¯ s = 3 P / and ( ∂P /∂s ) n = T / .Moreover, we have used integration by parts to simplifythe expressions: for example, the contributions from themiddle two terms in Eq. (50) with V ext ( r ⊥ ) can be shownto cancel with each other with our polynomial ansatz.For a given value of µ /k B T (or T /T F , see Eq.(12)), theweighted mass moments and spring constants can be cal-culated by using local thermodynamic variables in Eqs.(14), (17), (18), (20) and (25). The detailed expressionsfor numerical calculations are listed in Appendix B.It is easy to see that the minimization of the action S is equivalent to solving S ( ω ) Λ = 0 , (72)or det S ( ω ) = 0 . The detailed numerical procedure ispresented in Appendix C. Once a solution (i.e., the modefrequency ω and the coefficient vector Λ ) is found, wecalculate the density fluctuation of the mode, using δn ( r ⊥ ) = − N p − X i =0 A i r ⊥ ddr ⊥ (cid:2) r ⊥ n ( r ⊥ ) ˜ r i +1 ⊥ (cid:3) . (73)We have performed numerical calculations for the num-ber of the variational parameter N p up to , for anygiven chemical potential µ /k B T or temperature T /T F .In the following, we first discuss the decoupled first andsecond sound. Then, we focus on the effect of the cou-pling between first and second sound and the densityfluctuation of second sound modes.0 / T/T F (a) (b) T/T F FIG. 6: (Color online) The mode frequencies of decoupled first(a) and second sound (b). The vertical grey lines show thecritical temperature of a three-dimensional trapped unitaryFermi gas, T c ≃ . T F . B. Decoupled first and second sound
In Figs. 6(a) and 6(b), we present mode frequenciesof decoupled first and second sound modes of a highlyoblate unitary Fermi gas, obtained by solving the indi-vidual action S ( a ) and S ( e ) , respectively.For all the first sound modes, except the lowest-lyingone, the mode frequency decreases monotonically withincreasing temperature, exhibiting the same tempera-ture dependence as the first sound modes in an isotropic[21] or highly elongated harmonic trap [16, 26]. Thelowest-lying breathing mode instead does not depend onthe temperature and takes an invariant mode frequency ω B = √ ω ⊥ . We note that such a temperature indepen-dence is a peculiarity of the unitary Fermi gas due to itsinherent scale invariance. Indeed, our variational ansatz u a ( r ⊥ ) = r ⊥ of the lowest-lying breathing mode is one ofthe few exact scaling solutions exhibited by the Landautwo-fluid hydrodynamic equations at unitarity [25]. Tounderstand this, we simply recall that the spring constant K ( a ) ij satisfies, K ( a ) ij M ( a ) ij = 32 ( i + 2) ( j + 2)( i + j + 2) ω ⊥ . (74)Thus, if i = 0 or j = 0 , we have K ( a ) ij = 3 ω ⊥ M ( a ) ij . Takinginto account the fact that K ( ae ) i =0 ,j = 0 or K ( ae ) j =0 ,i = 0 , itis readily seen from Eq. (66) that ω = 3 ω ⊥ provides anexact solution, regardless the number of the variationalansatz used. We note also that, the first sound solu-tions converge very quickly with the number of the vari-ational parameters, indicating that these solutions areindeed well-approximated by the polynomial function.For second sound modes, we find that the mode fre-quency initially increases with increasing temperatureand then drops to zero very dramatically when the tem- / T/T F FIG. 7: (Color online) Temperature dependence of the fulltwo-fluid hydrodynamic mode frequencies (blue circles). Forcomparison, we show also the mode frequencies of decoupledfirst and second sound, respectively, by black solid and reddashed lines. The vertical grey lines show the critical tem-perature of a three-dimensional trapped unitary Fermi gas, T c ≃ . T F . perature approaches to the critical value. This behav-ior differs from the earlier result in a three-dimensionalisotropic harmonic trap [21] but agrees with a recentprediction made for a high elongated configuration [26].Qualitatively, we may estimate the discretized secondsound mode frequency by using the expression ω ∼ u q with a characteristic wavevector q ∼ /R s , where R s is the size of the superfluid component along the ra-dial direction. Within the local density approxima-tion, we have R s ∼ R F p − T /T c . Using Eq. (60) u ∼ (1 − T /T c ) α/ / , we obtain that ω ∼ (cid:18) − TT c (cid:19) α/ − / . (75)Thus, for the critical exponent α > / , the second soundmode frequency ω vanishes at T c . The superfluid densitydata (i.e., of superfluid helium) that we use have a crit-ical exponent α ≃ / , and consequently, we find thevanishing frequency at the transition. C. Full solutions of 2D two-fluid hydrodynamics
We now include the coupling term S ( ae ) . In Fig. 7,we report the full variational results by blue circles. Forcomparison, the decoupled first and second sound modefrequencies are also shown, by black lines and red dashedlines, respectively. As anticipated, the exact solutionwith frequency ω = √ ω ⊥ remains unchanged with theinclusion of the coupling term.Very similar to the sound velocities in the uniform case(see Fig. 4), the first sound mode frequency is barely1 / T/T F (a) (b) T/T F FIG. 8: (Color online) Enlarged view of the full two-fluidhydrodynamic mode frequencies (blue circles), near the n = 2 first sound mode (a) and the lowest two second sound modes(b). In the left panel, the avoided crossing between first andsecond sound is evident. affected by the coupling term S ( ae ) . This is particularlyevident in Fig. 8(a), in which we show an enlarged viewfor the n = 2 first sound mode (the integer n labels the n -th first sound modes). The correction to the first soundmode frequency due to the coupling is about . andoccurs only at around T ∼ . T F .On the other hand, the frequency of second soundmodes is notably pushed down by the coupling term S ( ae ) , as can be seen from Fig. 8(b). The maximumcorrection is up to when the temperature is about . T F . In the uniform case, i.e., Fig. 4, we find a sim-ilar correction to the second sound velocity in the sametemperature regime. D. Density fluctuations of sound modes
The sizable correction in the mode frequency of a sec-ond sound mode strongly indicates that, its density fluc-tuation, as a result of the coupling to first sound modes,should also be significant. In Fig. 9(b), we show thedensity fluctuations of the lowest two first sound modes( m = 2 and ) and of the lowest two second sound modes( m = 1 and ), at a temperature T = 0 . T F . Here, wehave used the integer m as the index of different modeswhen the coupling term is taken into account.Remarkably, the amplitude of the second sound den-sity fluctuations is just a bit smaller than that of thefirst sound modes, over a wide range of temperatures,revealing that it could be detected experimentally. Inthis respect, we note that, most recently the density fluc-tuations of the low-lying first sound modes have beenmeasured in a highly elongated unitary Fermi gas, to areasonably good precision [16]. Therefore, it is very likelythat a low-lying second sound mode could also be ob-served by looking at its density fluctuation, after proper (b) r n ( r ) r /R F n s n ( r ) / n F n tot T=0.18T F (a) FIG. 9: (Color online) (a) Density distribution (the blackline) and superfluid density distribution (the red dashed line)at T = 0 . T F , in units of the peak linear density of an idealFermi gas at zero temperature and at the trap center ( n F ). (b)Density fluctuations (in arbitrary unit) of the lowest two first-sound (solid lines) and second-sound modes (dashed lines). excitation.Comparing the density and superfluid density profiles,shown in Fig. 9(a), we find that the density fluctuationof second sounds is most significant within the superfluidcore, in accordance with their temperature-wave nature.In contrast, the density fluctuation of first sound modesextends over the whole Fermi cloud. E. Dependence on the superfluid density
We not that, a complete theoretical description of thesuperfluid density for a Fermi gas at unitarity is yet to bedetermined due to the theoretical difficulties in handlingthe strong interaction. In the above studies, we haveconsidered the superfluid fraction of superfluid heliumfor the superfluid density of the unitary Fermi gas, see,for example, Eqs. (27) and (28). As the second soundis arguably the most significant demonstration of super-fluid, it is natural to anticipate that the mode frequencyof second sound modes should depend very sensitively onthe form of superfluid density.We have performed numerical calculations with othermodels of the superfluid density, as shown in Fig. 10.The two additional models we considered are: firstly,2
He II 1 (T/T c ) (1 T/T c ) n s / n T/T c FIG. 10: (Color online) Bulk superfluid fraction: the blackline and blue dashed line correspond to the choices, n s /n =1 − ( T /T c ) and n s /n = (1 − T /T c ) / , respectively. The redcircles are the superfluid fraction of superfluid helium [47].Note that, recently the superfluid fraction of a 3D unitaryFermi gas has been calculated by Salasnich based on the Lan-dau’s expression for superfluid density and physical elemen-tary excitations [29]. The predicted temperature dependenceis similar to what we have shown in the figure, presumablydue to the simialr elementary excitations arising from stronginteractions. ( n s /n ) D = 1 − ( T /T c ) , which is in magnitude simi-lar to the superfluid fraction of superfluid helium butwith a mean-field-like critical exponent α = 1 ; secondly, ( n s /n ) D = (1 − T /T c ) / , which differs significantly fromthe superfluid fraction of superfluid helium but takes intoaccount the correct critical exponent α ≃ / . The pre-dicted mode frequencies for different superfluid densityare shown in Fig. 11.It is readily seen that the mode frequency of secondsound modes displays a very sensitive dependence onboth the magnitude and critical exponent of the super-fluid density. In particular, the mode frequency withthe superfluid-helium-like superfluid fraction (Fig. 11(a))differs significantly from that with the choice ( n s /n ) D =1 − ( T /T c ) , although these two superfluid fractions dif-fer slightly in magnitude. The strong dependence is veryencouraging, indicating that practically, the superfluiddensity of a unitary Fermi gas could be accurately de-termined by measuring the mode frequency of low-lyingsecond sound modes. F. Experimental considerations
Ideally, a compressional second sound mode would beexcited by modulating the trapping frequency close tothe predicted mode frequency. However, in our quasi-2Dconfiguration, this can be problematic due to the factthat the radial confinement is generated by a magnetic field, whose modulation can also change the s -wave scat-tering length and hence the system may be away fromthe unitary regime. We note that this problem does notexist in other setups where it is possible to have opticalradial (and axial) confinement.For our experimental setup at Swinburne, therefore, amore practical scheme would be to locally perturb thecloud and allow it to relax. Taking inspiration from theselective excitation schemes, developed by the Innsbruckexperiments [17], one may direct an appropriately sizedblue-detuned laser beam into the cloud and apply a mod-ulation burst, chosen to provide the best mode matchingwith the calculated second sound mode. The power, du-ration and shape of the laser beam should be optimizedin order to resonantly derive the desired small amplitudeoscillation in the linear response regime. VI. CONCLUSIONS
In conclusion, we have derived the two-dimensionalsimplified Landau two-fluid hydrodynamic equations todescribe the low-energy dynamics of a unitary Fermigas confined in highly oblate harmonic traps. By us-ing a variational approach, these two-fluid hydrodynamicequations have been fully solved. We have discussedin detail the resulting density-wave (firsts sound) andtemperature-wave (second sound) oscillations, in the ab-sence or presence of a weak transverse harmonic trap.First and second sound velocities or discretized mode fre-quencies have been predicted, accordingly.We have found a weak coupling between first and sec-ond sound in highly oblate unitary Fermi gas, very similarto the case of superfluid helium. Though the coupling isweak, it induces significant density fluctuations for sec-ond sound modes, indicating that second sound couldbe potentially observed in the highly oblate configura-tion, by measuring the density fluctuations. Owing tothe strong sensitivity of second sound mode frequenciesto superfluid density, the experimental measurement ofdiscretized second sound modes could provide a promis-ing way of accurately determining the superfluid densityof a unitary Fermi gas.
Acknowledgments
We thank Martin W. Zwierlein and Mark J.-H. Kufor providing us their experimental data. This re-search was supported by the ARC Discovery ProjectsGrant Nos. FT130100815 (HH), DP140103231 (HH),DE140100647 (PD), FT120100034(CJV), DP130101807(CJV) and DP140100637 (XJL), and NFRP-China GrantNo. 2011CB921502 (XJL, HH).3 (c) n s /n = (1 T/T c ) T/T F (b) n s /n = 1 (T/T c ) T/T F / T/T F (a) He II FIG. 11: (Color online) Sensitivity of second sound modes on the bulk superfluid density. The vertical grey lines show thecritical temperature of a three-dimensional trapped unitary Fermi gas, T c ≃ . T F . Appendix A: Equation of state of a unitary Fermigas
In the low temperature superfluid phase, we fit theexperimental data of g ( x ) by using the Padé approximantof order [2/2], g s ( x ) ≃ ξ − / (cid:20) p x − + p x − q x − + q x − (cid:21) , (A1)where ξ ≃ . ± . is the Bertsch parameter [37], andobtain p = − . , p = +5 . , q = − . and q = +5 . .For the normal state, i.e., x ⊂ [ − . , x c ≃ . , we fitthe data with a fourth order polynomial g ( x ) ≃ a + a x + a x + a x + a x (A2)and find that a = +1 . , a = +0 . , a = − . , a = − . and a = 0 . . At hightemperatures where x < − . , we use the virial expan-sion form, g ( x ) ≃ b e x + 3 b e x − − / e x + 3 − / e x , (A3)where b = 3 √ / and b ≃ . are the secondand third virial coefficients of a unitary Fermi gas [45,46], respectively. To make g ( x ) smooth across the wholenormal state, we connect g ( x ) and g ( x ) by using theexpression, g n ( x ) = g ( x )1 + e − x +1 . + g ( x )1 + e +4( x +1 . . (A4)Recalling that g s ( x ) = g n ( x ) at x c ≃ . , we thereforeset x c = 2 . . We note that, we do not needto impose the constraint g ′ s ( x c ) = g ′ n ( x c ) , due to thesuperfluid phase transition. Appendix B: Matrix elements of S ( ω ) In this appendix, we present the weighted mass mo-ments and spring constants in their dimensionless form, ˜ M ij = M ij π / ~ ω z m R F ( k B T ) (B1)and ˜ K ij = K ij π / ~ ω z m R F ( k B T ) ω ⊥ , (B2)respectively, for the purpose of performing numerical cal-culations. Accordingly, we solve the matrix S ( ω ) in itsdimensionless form, h ˜ S (˜ ω ) i ij ≡ " ˜ M ( a ) ij ˜ ω − ˜ K ( a ) ij − ˜ K ( ae ) ij − ˜ K ( ae ) ji ˜ M ( e ) ij ˜ ω − ˜ K ( e ) ij , (B3)where ˜ ω ≡ ω/ω ⊥ is the reduced frequency. For givenchemical potential µ /k B T or temperature T /T F , thesedimensionless matrix elements are given by, after taking y = ˜ r ⊥ , ˜ M ( a ) ij = ˆ ∞ dyy i + j +22 [ f n ] , (B4) ˜ M ( e ) ij = ˆ ∞ dyy i + j +22 (cid:20)(cid:18) − f s f n (cid:19) f s (cid:21) , (B5)for the weighted mass moments and ˜ K ( a ) ij = 32 ( i + 2) ( j + 2)( i + j + 2) ˜ M ( a ) ij , (B6) ˜ K ( ae ) ij = i ( i + 2)4 TT F ˆ ∞ dyy i + j [¯ s f s ] , (B7) ˜ K ( e ) ij = T T F ˆ ∞ dyy i + j (cid:20) ( i + 2) ¯ s f s − T F T (¯ s f s ) ′ y (cid:21) × (cid:20) ( j + 2) ¯ s f s − T F T (¯ s f s ) ′ y (cid:21) f n ¯ c v ] , (B8)4for the spring constants. In the above expressions, theargument for the universal scaling functions is x = µ / ( k B T ) − y/ ( T /T F ) , i.e., [ f n ] is the short-hand nota-tion of f n [ µ / ( k B T ) − y/ ( T /T F )] , and (¯ s f s ) ′ stands forthe derivative d [¯ s f s ] ( x ) /dx. Appendix C: Solving det[ ˜ S (˜ ω )] = 0 To solve the matrix equation ˜ S ( ˜ ω )Λ = 0 , (C1)where the vector of displacement fields, Λ =[ A , B ..., A i , B i , ... ] T , we rewrite the matrix ˜ S (˜ ω ) = M ˜ ω − K . Here, M and K denote collectively the ma-trix of the dimensionless weighted mass moments andthe spring constants, respectively. The matrix M is pos-itively definite, so that we decouple it by a product of alower triangular matrix L and its transpose, M = L · L T . (C2) In terms of this decomposition, the matrix equation, Eq.(C1), becomes h L − · K · (cid:0) L − (cid:1) T i · L T Λ =˜ ω L T Λ . (C3)It is easy to see that the matrix [ L − · K · ( L − ) T ] issymmetric and its eigenvalues give rise to the desired so-lution of mode frequencies. The displacement field foreach eigenvalue can also be calculated accordingly, withknown eigenstate of the matrix [ L − · K · ( L − ) T ] . Moreexplicitly, if we denote an eigenstate by X , the corre-sponding displacement field is given by Λ = (cid:0) L T (cid:1) − X . (C4)Once the displacement fields for a mode are found, wemay calculate its density fluctuation and temperaturefluctuation. [1] L. Tisza, C. R. Phys. , 1035 (1938).[2] L. D. Landau, J. Phys. (USSR) , 71 (1941).[3] I. M. Khalatnikov, An Introduction to the Theory of Su-perfluidity (Westview Press, New York, 2000).[4] A. Griffin, T. Nikuni, and E. Zaremba,
Bose-CondensedGases at Finite Temperatures (Cambridge UniversityPress, Cambridge, 2009).[5] J. C. Findlay, A. Pitt, H. G. Smith, and J. O. Wilhelm,Phys. Rev. , 506 (1938).[6] V. P. Peshkov, J. Phys. (USSR) , 389 (1946).[7] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod.Phys. , 1215 (2008).[8] S. Stringari, Europhys. Lett. , 749 (2004).[9] J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, andJ. E. Thomas, Phys. Rev. Lett. , 150402 (2004).[10] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C.Chin, J. Hecker Denschlag, and R. Grimm, Phys. Rev.Lett. , 203201 (2004).[11] H. Hu, A. Minguzzi, X.-J. Liu, and M. P. Tosi, Phys.Rev. Lett. , 190403 (2004).[12] P. Capuzzi, P. Vignolo, F. Federici, and M. P. Tosi, Phys.Rev. A , 021603(R) (2006).[13] A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R.Geursen, M. Bartenstein, C. Chin, J. Hecker Denschlag,and R. Grimm, Phys. Rev. Lett. , 170401 (2007).[15] S. K. Adhikari and L. Salasnich, New J. Phys. , 023011(2009).[16] M. K. Tey, L. A. Sidorenkov, E. R. Sánchez Guajardo,R. Grimm, M. J. H. Ku, M. W. Zwierlein, Y. H. Hou, L.Pitaevskii, and S. Stringari, Phys. Rev. Lett. , 055303(2013).[17] E. R. Sánchez Guajardo, M. K. Tey, L. A. Sidorenkovand R. Grimm, Phys. Rev. A , 063601 (2013). [18] E. Taylor and A. Griffin, Phys. Rev. A , 053630 (2005).[19] Y. He, Q. Chen, C. C. Chien, and K. Levin, Phys. Rev.A , 051602(R) (2007).[20] E. Taylor, H. Hu, X.-J. Liu, and A. Griffin, Phys. Rev.A , 033608 (2008).[21] E. Taylor, H. Hu, X.-J. Liu, L. P. Pitaevskii, A. Griffin,and S. Stringari, Phys. Rev. A , 053601 (2009).[22] R. Watanabe, S. Tsuchiya and Y. Ohashi, Phys. Rev. A , 043630 (2010).[23] H. Hu, E. Taylor, X.-J. Liu, S. Stringari, and A. Griffin,New J. Phys. , 043040 (2010).[24] G. Bertaina, L. P. Pitaevskii, and S. Stringari, Phys. Rev.Lett. , 150402 (2010).[25] Y.-H. Hou, L. P. Pitaevskii, and S. Stringari, Phys. Rev.A , 033620 (2013).[26] Y.-H. Hou, L. P. Pitaevskii, and S. Stringari, Phys. Rev.A , 043630 (2013).[27] H. Hu and X.-J. Liu, Phys. Rev. A , 053605 (2013).[28] L. A. Sidorenkov, M. K. Tey, R. Grimm, Y.-H. Hou, L.Pitaevskii, and S. Stringari, Nature (London) , 78(2013).[29] L. Salasnich, Phys. Rev. A , 063619 (2010).[30] G. Baym and C. J. Pethick, Phys. Rev. A , 043631(2013).[31] K. Martiyanov, V. Makhalov, and A. Turlapov, Phys.Rev. Lett. , 030404 (2010).[32] B. Fröhlich, M. Feld, E. Vogt, M. Koschorreck, W. Zw-erger, and M. Köhl, Phys. Rev. Lett. , 105301 (2011).[33] P. Dyke, E. D. Kuhnle, S. Whitlock, H. Hu, M. Mark,S. Hoinka, M. Lingham, P. Hannaford, and C. J. Vale,Phys. Rev. Lett. , 105304 (2011).[34] A. A. Orel, P. Dyke, M. Delehaye, C. J. Vale, and H. Hu,New J. Phys. , 113032 (2011).[35] A. T. Sommer, L. W. Cheuk, M. J. H. Ku, W. S. Bakr,and M. W. Zwierlein, Phys. Rev. Lett. , 045302 (2012).[36] E. Vogt, M. Feld, B. Fröhlich, D. Pertot, M. Koschorreck,and M. Köhl, Phys. Rev. Lett. , 070404 (2012).[37] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W.Zwierlein, Science , 563 (2012).[38] T.-L. Ho, Phys. Rev. Lett. , 090402 (2004).[39] H. Hu, P. D. Drummond, and X.-J. Liu, Nature Phys. ,469 (2007).[40] X.-J. Liu and H. Hu, Phys. Rev. A , 063613 (2005).[41] H. Hu, X.-J. Liu, and P. D. Drummond, Europhys. Lett. , 574 (2006).[42] H. Hu, X.-J. Liu, and P. D. Drummond, Phys. Rev. A , 023617 (2006).[43] H. Hu, X.-J. Liu, and P. D. Drummond, Phys. Rev. A , 061605(R) (2008).[44] H. Hu, X.-J. Liu, and P. D. Drummond, New J. Phys. , 063038 (2010).[45] X.-J. Liu, H. Hu, and P. D. Drummond, Phys. Rev. Lett. , 160401 (2009).[46] X.-J. Liu, Phys. Rep. , 37 (2013).[47] J. G. Dash and R. D. Taylor, Phys. Rev. , 7 (1957).[48] K. R. Atkins, Phys. Rev. , 962 (1959).[49] P. R. Zilsel, Phys. Rev. , 309 (1950).[50] In the time domain, the displacement fields u s and u n aredefined by, v s ⊥ ( r ⊥ , t ) = ∂ u s ( r ⊥ , t ) /∂t and v n ⊥ ( r ⊥ , t ) = ∂ u n ( r ⊥ , t ) /∂t .[51] G. Reinaudi, T. Lahaye, Z. Wang, and D. Guéry-Odelin,Opt. Lett. , 3143 (2007).[52] J. Estève, C. Gross, A. Weller, S. Giovanazzi, and M. K.Oberthaler, Nature (London)455