First and second sound of a unitary Fermi gas in highly oblate harmonic traps

We theoretically investigate first and second sound modes of a unitary Fermi gas trapped in a highly oblate harmonic trap at finite temperatures. Following the idea by Stringari and co-workers {[}Phys. Rev. Lett. \textbf{105}, 150402 (2010){]}, we argue that these modes can be described by the simplified two-dimensional two-fluid hydrodynamic equations. Two possible schemes - sound wave propagation and breathing mode excitation - are considered. We calculate the sound wave velocities and discretized sound mode frequencies, as a function of temperature. We find that in both schemes, the coupling between first and second sound modes is large enough to induce significant density fluctuations, suggesting that second sound can be directly observed by measuring \textit{in-situ} density profiles. The frequency of the second sound breathing mode is found to be highly sensitive to the superfluid density.


I. INTRODUCTION
Low-energy excitations of a strongly interacting superfluid quantum system -in which inter-particle collisions are sufficiently strong to ensure local thermodynamic equilibrium -can be described by the Landau two-fluid hydrodynamic theory [1][2][3][4]. There are two well-known kinds of excitations, referred to as first and second hydrodynamic modes, which describe the coupled oscillations of the superfluid and normal fluid components at finite temperatures. First sound is an ordinary phenomenon in liquid, describing the propagation of a pressure or density wave that is associated with normal acoustic sound. The motions of the superfluid and of normal fluid components are locked in phase. In contrast, second sound is a phenomenon characteristic of superfluids. It describes the ability to propagate undamped entropy oscillations, which are essentially opposite phase oscillations of the two components. In the quantum liquid of superfluid helium, the velocities of first and second sound were first measured in 1938 by Findlay et al. [5] and in 1946 by Peshkov [6], respectively.
Ultracold atomic Fermi gases near a Feshbach resonance represent a new type of strongly interacting quantum system with unprecedented control over interatomic interactions, dimensionality and purity [7]. These are anticipated to provide an ideal table-top system to deepen our understanding of fermionic superfluidity and Landau two-fluid hydrodynamics. For a Fermi gas at unitarity, first sound modes have been investigated in great detail, both experimentally and theoretically [8][9][10][11][12][13][14][15][16][17], particularly near zero temperature. In contrast, second sound mode is more difficult to address [18][19][20][21][22][23][24][25][26][27] and has only recently been observed in experiments [28]. Key to understanding second sound, is knowledge of the superfluid density [20,29,30], a quantity which is still a subject of intense investigation. Theoretically, first and second sound modes of the Landau two-fluid hydrodynamic equations were solved by Taylor et al. for an isotropically trapped unitary Fermi gas [21], using an assumed superfluid den-sity and by developing a variational approach. For general anisotropic harmonic traps, the Landau hydrodynamic equations are more difficult to solve. However, for a highly elongated configuration, it was shown by Bertaina, Pitaevskii and Stringari [24] that, a simplified one-dimensional (1D) two-fluid hydrodynamic description may be derived, generalizing the dimensional reduction at zero temperature [15]. Following this pioneering idea, the propagation of a second sound wave was strikingly observed in a highly elongated unitary Fermi gas by the the Innsbruck team [28]. Using the measured second sound velocity data, the simplified 1D hydrodynamic equations were solved and used to extract the superfluid density [28]. These results show that studies of second sound provide a promising route towards an accurate determination of the superfluid density.
Here, we would consider a unitary Fermi gas in highly oblate harmonic traps, which can be readily realised in current experiments. Indeed, such configurations have already been investigated in laboratory experiments [31][32][33][34][35][36]. The analysis presented here could therefore be directly tested in future experiments. Our main results are briefly summarized as follows. Following the ideas of Stringari and co-workers [24,26], we derive the simplified two-dimensional (2D) Landau two-fluid hydrodynamic equations. Using a variational approach within the local density approximation [18,20,21], we fully solve the coupled 2D hydrodynamic equations in the presence of a weak radial harmonic confinement. The discretized mode frequencies of first and second sound are calculated. We find that the density fluctuation due to the second sound mode is significant, revealing that second sound is directly observable from in-situ density profiles, after an appropriate modulation of the weak radial trapping potential. The mode frequency of second sound is found to depend very critically on the form of the superfluid density. Our quasi-2D analysis is similar to the quasi-1D study of Stringari and co-workers [26], however, there are a number of significant differences. Specifically we fully solve the coupled hydrodynamic equations without assuming that first and second sound are decoupled. This enables us to calculate the density fluctuations associated with the discretized second sound modes. From an experimental point of view, the possibility of measuring discrete breathing mode frequencies may allow a more accurate determination of superfluid density.
The remainder of the paper is organized as follows. In the next section, we outline the reduced 2D universal thermodynamics, satisfied by the unitary Fermi gas in highly oblate harmonic traps based on the measured 3D equation of state [37]. In Sec. III, we derive the simplified 2D Landau two-fluid hydrodynamic equations and discuss very briefly their applicability. In Sec. IV, we consider the propagation of first and second sound, which may be excited by creating a short perturbation of the density. The nature of first and second sound excitations in 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 frequencies of first and second sound. The density fluctuation associated with low-lying second sound modes is discussed and the dependence of the second sound mode frequency on the form of superfluid density is examined. Finally, in Sec. V we draw our conclusions. The appendix presents further details of our numerical calculations.

II. 2D REDUCED THERMODYNAMICS
Let us consider a unitary Fermi gas trapped in a highly anisotropic pancake-like harmonic potential, with atomic mass m and the radial and axial trapping frequencies ω ⊥ , ω z . The aspect ratio of the trap, defined by λ = ω z /ω ⊥ ≫ 1, can be specifically tuned in experiments. The number of atoms in the typical sound mode experiments [9,10,13,16,28] is about N ∼ 10 5 . For such a large number of atoms, it is standard to use the local density approximation [7], which amounts to treating the atoms in the position (r ⊥ , z) locally as a uniform matter, with a local chemical potential given by, Here µ 0 is the chemical potential at the trap center. The validity of the local density approximation can be conveniently estimated by comparing the number of atoms N ∼ 10 5 with a threshold N 2D = λ (λ + 1), below which there is a dimensional crossover from three-to twodimensions [33,34]. At Swinburne, we have previously demonstrated pancake traps with λ ∼ 60. Thus, the ratio N/N 2D ∼ 30 , indicating that the system could be deeply in the 3D limit where the local density approximation is well justified.

A. Universal thermodynamics
A remarkable feature of a uniform unitary Fermi gas is that all its thermodynamic functions can be expressed as universal functions a single dimensionless parameter µ/(k B T ) [38,39]. This is simply due to the fact that at unitarity, the s-wave scattering length -the only length scale used to characterize the short-range interatomic interaction -diverges. Thus, the remaining length scales are the thermal wavelength and the mean distance between atoms n −1/3 , where n is the number density. Accordingly, the energy scales are given by the temperature, k B T , and by the Fermi energy, or, alternatively, by the chemical potential µ. Following the scaling analysis, all the thermodynamic functions can therefore be expressed in terms of universal functions that only depends on the dimensionless parameter x = µ/(k B T ). For the pressure and number density, these universal functions are given, respectively, by Using the thermodynamic relation n = (∂P/∂µ) T , we find f 3D n (x) = df 3D p (x)/dx. It is easy to see that knowing f 3D p and f 3D n , we can calculate directly all the thermodynamic functions of the uniform unitary Fermi gas [26].
It is highly non-trivial to theoretically determine the universal function f 3D p or f 3D n , since there are no small parameters to control the theory of a strongly interacting Fermi gas [40][41][42][43][44], except at high temperatures, where the virial expansion approach is applicable [45,46]. Fortunately, accurate experimental data for the equation of state are now available from the landmark experiments performed by the Massachusetts Institute of Technology (MIT) team [37]. In Fig. 1, we present their data for the universal scaling function, where the subscript "0" indicates the result of an ideal, non-interacting Fermi gas and For later convenience of numerical calculations, we have fitted the experimental data to an analytic expression, as detailed in the Appendix A. The fitting curve is smoothly extrapolated to both low and high temperature regimes where the behavior of f 3D n (x) is well-understood [23,25,45,46]. As evident in Fig. 1, the relative error in the fitting is less than 0.5%, significantly smaller than the standard error for the experimental data (i.e., 2%).

B. Local density approximation
Within the local density approximation, we can write the local pressure and number density using the universal functions, where the local chemical potential is given by Eq. (2). For a trapped Fermi gas in three dimensions, the Fermi temperature T F is defined as [7] Using the number equation N =´dr ⊥ dzn(r ⊥ , z), we relate µ 0 /k B T to the reduced temperature T /T F by, The onset of superfluidity in a trapped unitary Fermi gas occurs at µ 0 /k B T = x c ≃ 2.49. Numerically, we find that T c ≃ 0.223T F , in agreement with the previous results [26,37].

C. 2D reduced thermodynamic functions
In highly oblate harmonic traps, although the cloud is still three-dimensional, its low-energy dynamics are greatly affected by the tight axial confinement. This situation is very similar to a highly elongated unitary Fermi gas considered earlier by Bertaina, Pitaevskii and Stringari [24], who showed that with tight radial confinement, the standard Landau two-fluid hydrodynamic equations could reduce to a simplified 1D form. For the same reason, which will be explained in greater detail in the next section, first and second sound under the tight axial confinement can be described by a simplified 2D two-fluid hydrodynamic description. In brief, due to the nonzero viscosity and thermal conductivity, local fluctuations in temperature (δT ) and chemical potential (δµ) become independent of the axial coordinate, for any low-energy excitations at the frequency ω ∼ ω ⊥ ≪ ω z . Therefore, the axial degree of freedom in all the thermodynamic variables that enter the Landau two-fluid hydrodynamic equations becomes irrelevant and can be integrated out. We can then derive 2D reduced thermodynamics and immediately have the reduced Gibbs-Duhem relation, where the variables P 2 , s 2 and n 2 are the axial integrals of their three-dimensional counterparts, namely the local pressure, entropy density and number density. The variable P 2 is given by, where and we have introduced the universal scaling function, All the 2D thermodynamic variables can be derived from the reduced Gibbs-Duhem relation, for example, where In addition, it is readily shown that the specific heat per particle at constant column density and pressure are given by [25], wheres 2 ≡ s 2 /(n 2 k B ) is the entropy per particle and f ′ n (x) ≡ df n (x)/dx. It is also straightforward to check the universal relations, In Fig. 2, we report the relevant 2D universal scaling functions, calculated by using the experimental MIT data for f 3D n (x) after smoothing.

D. Superfluid density
We can also express the local superfluid density as a universal (but as yet undetermined) function f 3D s : Integrating out the axial coordinate, we obtain where the universal scaling function f s (x) is given by, As mentioned earlier, in contrast to the equation of state, the universal function for the superfluid density of a unitary Fermi gas f 3D s (x) is not yet known precisely  [47] and the corresponding 2D reduced superfluid fraction (the black line). The inset shows the 2D universal scaling function for superfluid density, to be used for a unitary Fermi gas. The vertical grey lines indicate the critical threshold for superfluidity, xc ≃ 2.49 [37]. [20,29,30]. For illustrative purposes in the present work, we will consider the superfluid fraction of superfluid helium [47], which is plotted in Fig. 3 by red circles. Our choice is motivated by the similarity for hydrodynamics between superfluid helium and unitary Fermi gas, presumably arising from strong correlations in both systems, and justified in part by measurements in Ref. [28]. Using Eq.
and f s (x) is then obtained. In the inset of Fig. 3, we show f s (x) calculated with the superfluid fraction of superfluid helium.
In a 2D configuration with ω ⊥ = 0, it is natural to define a Fermi temperature in terms of the column density n 2 : which coincides with the usual three-dimensional definition of the Fermi temperature, Eq. (4) . The reduced temperature is then given by, Therefore, the critical temperature of the unitary Fermi gas in the pancake geometry is given by, Fig. 3, we show the 2D superfluid fraction f s /f n (the black line) as a function of T /T 2D c . It lies systematically below the threedimensional superfluid fraction (red circles), due to the integration over the axial degree of freedom.
It should be noted that the critical behavior of the superfluid density near the phase transition is greatly affected by the axial integration. In three dimensions, we may define the critical exponent α by After integrating out the axial degree of freedom, it is easy to show, As we shall see, the increase in the critical exponent will reduce the second sound velocity and hence the second sound mode frequency.

III. 2D SIMPLIFIED TWO-FLUID HYDRODYNAMIC EQUATIONS
We now derive the simplified 2D Landau two-fluid hydrodynamic equations, following the idea by Bertaina, Pitaevskii and Stringari [24]. This 2D picture will generally be valid as we are considering excitations whose wavelength is long compared to the transverse (axial) cloud width.
The standard two-fluid equations in the harmonic trap V ext are given by, 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, , and finally η and κ are the shear viscosity and thermal conductivity, respectively. In the above equations, we have kept only linear terms in the velocity, as we are interested in small-amplitude dynamics in the linear response regime. Moreover, we have omitted bulk viscosity terms which give smaller contributions. The Landau two-fluid hydrodynamic theory is applicable when the mean free path l is much smaller than the wavelength of the sound wave λ. By considering the axial size of the Fermi cloud R z ∼ / (mω z ), we further require R z ≪ λ, therefore the hydrodynamic regime in a harmonic trap is achieved when l ≪ R z ≪ λ.
As discussed by Stringari and co-workers [24], in the presence of tight confinement, the nonzero viscosity and thermal conductivity may significantly change the lowenergy dynamics at the frequency ω ∼ ω ⊥ ≪ ω z . This occurs when the viscous penetration depth δ -the typical length scale at which an excitation becomes damped due to shear viscosity -fulfils the condition, In the case of superfluid helium in a thin capillary, the above condition makes the normal component of the liquid stick to the wall and thus the normal velocity field vanishes [48]. With a soft wall caused by the tight axial confinement, the normal component can move but the normal velocity field becomes uniform along the axial direction [24]. For the same reason, the superfluid velocity also becomes independent of the axial coordinate. Analogously, a nonzero thermal conductivity in Eq. (35) for the entropy conservation leads to the independence of the temperature fluctuation δT on the axial coordinate. The tight axial confinement also implies that the axial component of the velocity fields, for both superfluid and normal fluid, must be much smaller than the transverse one and may be neglected. Using Eq. (33) for the superfluid velocity implies that the chemical potential fluctuation δµ is essentially independent of the axial coordinate. In brief, we conclude that, owing to the crucial roles played by the shear viscosity and thermal conductivity, under the tight axial confinement the velocity fields are independent of the axial position and a global thermal equilibrium along the axial direction is established. With these observations, it is straightforward to write down the simplified two-fluid hydrodynamic equations by integrating out the axial degree of freedom in Eqs. (32), (33), (34) and (35), which take the following forms, where the current density now becomes j ⊥ = m(n s2 v s⊥ + n n2 v n⊥ ) and we have omitted the residual dissipation terms along the weakly confined direction, as in a uniform fluid the effect of viscosity and thermal conductivity is irrelevant in the long-wavelength limit. We note that, by introducing a characteristic collisional time τ related to viscosity, τ ≃ η/(mnv 2 ) ∼ η/(n ω z ), wherev is the average velocity of atoms and is of the order of the Fermi velocity v F ∼ ω z /m, it is easy to check that Eq. (36) is equivalent to requiring the low frequency condition, Recalling that ω ∼ ω ⊥ and ω z = λω ⊥ , the above requirement for achieving the reduced 2D hydrodynamics can be rewritten as λ 2 ωτ ≫ 1. For a highly oblate trap with typical aspect ratio λ ∼ 60 ≫ 1, this condition is compatible with the condition to enter the hydrodynamic regime, l ≪ λ, which can alternatively be written as ωτ ≪ 1.

A. Variational reformulation of the simplified Landau hydrodynamic equations
A convenient way to solve the simplified two-fluid hydrodynamic equations is to reformulate them using Hamilton's variational principle [49]. Following previous work [18,20], the hydrodynamic modes of these equations with frequency ω at temperature T can be obtained by minimizing a variational action, which, in terms of displacement fields u s (r ⊥ ) = iωv s⊥ (r ⊥ ) and u n (r ⊥ ) = iωv n⊥ (r ⊥ ) [50], takes the following form, Here, n s2 (r ⊥ ) and n n2 (r ⊥ ) = n 2 (r ⊥ ) − n s2 (r ⊥ ) are the 2D reduced superfluid and normal-fluid densities at equilibrium, as discussed in the previous section Sec. II. The density fluctuation δn and entropy fluctuation δs are given by, and respectively. The effect of the weak radial trapping potential V ext (r ⊥ ) = mω 2 ⊥ r 2 ⊥ /2 enters the action Eq. (42) through the coordinate dependence of the equilibrium thermodynamic variables (∂µ/∂n 2 ) s2 , (∂T /∂n 2 ) s2 and (∂T /∂s 2 ) n2 , within the local density approximation. We stress that, all these thermodynamic variables can be obtained 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 twofluid hydrodynamic equations can be well understood as density and entropy (temperature) waves, which are the pure in-phase mode with u s = u n and the pure out-ofphase mode with n s u s + n n u n = 0, known as first and second sound, respectively [2][3][4]. For a strongly interacting unitary Fermi gas, we could use the same classification [21]. For this purpose, we rewrite the action Eq. (42) in terms of two new displacement fields u a = (n s2 u s + n n2 u n ) /n 2 (45) and considering that the density and temperature fluctuations are given by and 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), after some straightforward but lengthy algebra, we arrive at where S (e) = mω 2 n s2 n n2 It is clear that the first and second sound are governed by the actions S (a) and S (e) , respectively. The coupling between first and second sound is controlled by the coupling term S (ae) , which is in general nonzero. Indeed, in our case, as (∂P 2 /∂s 2 ) n2 = T /2, strictly speaking the first and second sound are coupled at any finite temperatures.

IV. FREE-PROPAGATING FIRST AND SECOND SOUND
For a uniform superfluid (V ext = 0), the solutions of S (a) and S (e) are plane waves of wave vector q with dispersion ω 1 = c 1 q and ω 2 = c 2 q, where and These expressions for first and second sound velocities are the standard results used to describe superfluid helium, when the corresponding equation of state and superfluid density are used [2][3][4]. In Fig. 4, we show the decoupled first and second sound velocities of a highly oblate Fermi gas at unitarity, by using dashed lines. Including the coupling term S (ae) and using the standard thermodynamic relations, it is straightforward to show that the solutions for sound velocities u of the simplified two-fluid hydrodynamic equations satisfy the following equation, It gives rise to two solutions for the sound velocity, u 1 and u 2 , which in the absence of the coupling term S (ae) (i.e., γ = 1), coincide with the decoupled first and second sound velocities, c 1 and c 2 . The numerical results for u 1 and u 2 are shown in Fig. 4 by solid lines. The temperature dependence of the reduced sound velocities in Fig. 4 is very similar to that of a 3D unitary Fermi gas predicted in the earlier works [23,29]. We attribute this qualitative similarity to the strongly interacting nature of the system. In the case of a small parameter θ ≡ c 2 2 /(γc 2 1 ) ≪ 1, which is indeed true for a highly oblate unitary Fermi gas, we may solve Eq. (55) perturbatively. We find the expansions [23], To the leading order of θ, the first sound velocity is not affected by the coupling term, but the second sound velocity decreases by a factor of √ γ and is now given by [23,26], Therefore, quantitatively, the coupling between first and second sound can be characterized by the so-called Landau-Placzek (LP) parameter ǫ LP ≡ γ − 1 [23]. In the inset of Fig. 4, we plot the LP parameter as a function of temperature. In the superfluid phase, it is always smaller than 0.4, indicating that the first and second sound couple very weakly in a unitary Fermi gas. A similar situation happens in superfluid helium, wherec p ≃c v and hence ǫ LP ≃ 0 [20]. From the above approximate expression for the second sound velocity and Eq. (31) for 2D superfluid density, it is clear the velocity vanishes as when approaching to the superfluid phase transition from below. Here, α ≃ 2/3 is the critical exponent of the superfluid density of a unitary Fermi gas in three dimensions.
Although the coupling between first and second sound is weak in a unitary Fermi gas, it is of significant importance for the purpose of experimental detection of second sound. Unlike superfluid helium, temperature oscillations -the characteristic motion of the second soundare difficult to observe in ultracold atomic gases. Density measurement is the most efficient way to characterize any low-energy dynamics. As a result, an observable second sound mode must have a sizable density fluctuation, which can only be induced by its coupling to first sound modes. The strength of density fluctuations in second sound may be conveniently characterized by the following ratio between the relative density and temperature fluctuations, where the subscript "2nd" indicates the second sound and we have taken the derivative at const pressure, as the second sound is properly considered as an oscillating wave at constant pressure, rather than at constant density, as indicated by Eq. (59) [26]. In contrast, the ratio between the relative density and temperature fluctuations in first sound is given by, In Fig. 5, we report the density fluctuation of second sound as a function of temperature in the superfluid phase, calculated with an assumed temperature fluctuation ratio δT /T = 10%. This is a typical temperature fluctuation, achievable in current experiments. For example, in the recent first-sound collective mode measurement, it was shown that (δn/n 2 ) 1st ∼ 20% over a wide 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 is easy to see that the density fluctuation of second sound is very significant when temperature T > 0.7T c ∼ 0.14T 2D F , revealing that the propagation of a second sound pulse could be experimentally detected via density measurements. This result suggests that a similar approach to that demonstrated by the Innsbruck experiment [28] could be applied in the oblate geometry considered here, where a focused blue-detuned laser propagating along the tightly confined direction could excite a second sound wave which will propagate radially outwards from the cloud centre. Subsequent absorption images taken at different times after the laser pulse could record the density fluctuation induced by the propagating second sound wave.
While the essential physics of our scheme is the same as that considered by by Stringari and co-workers in highly elongated trap [26], the quasi-2D geometry may offer ad-vantages for experiments. Specifically, clouds confined in elongated traps typically have very high peak densities leading to high optical densities that are difficult to measure quantitatively using absorption imaging [51,52]. This means quantifying small density fluctuations in such clouds can prove challenging. In a quasi-2D trap however, such high peak optical densities can be avoided and images taken along the direction of tight confinement can be integrated over the angular coordinate to obtain an accurate measure of the radial density. Small density fluctuations should therefore be more easily detected.

V. BREATHING FIRST AND SECOND SOUND MODES IN HARMONIC TRAPS
We now consider a weak transverse harmonic trap, V ext (r ⊥ ) = mω 2 ⊥ r 2 ⊥ /2 and fully solve the simplified Landau two-fluid hydrodynamic equations using a variational approach, developed in the previous work [20,21]. We focus on compressional (breathing) modes with the projected angular momentum l z = 0, since these modes are the easiest one to excite experimentally.

A. Variational approach
For breathing modes, we assume the following polynomial ansatz for the displacement fields: wherer ⊥ is the unit vector in the radial direction, {A i , B i } (i = 0, · · · , N p − 1) are the 2N p variational parameters, andr ⊥ ≡ r ⊥ /R F is the dimensionless radial coordinate with R F being the Thomas-Fermi radius. By inserting this variational ansatz into the action Eq. (49), we express the action S as functions of the 2N p variational parameters {A i , B i } . The mode frequencies are obtained by minimizing the action S with respect to these 2N p parameters. The accuracy of our variational calculations 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, the expression of the action can be written in a compact form, where Λ ≡ A 0 , B 0 , · · · , A i , B i , · · · , A Np−1 , B Np−1 T and S(ω) is a 2N p × 2N p matrix with block elements (i, j = 0, · · · , N p − 1), In [S (ω)] ij , we have introduced the weighted mass moments, and the spring constants, To derive the above equations, we have used the universal relations satisfied by the highly oblate unitary Fermi gas: n 2 (∂P 2 /∂n 2 )s 2 = 3P 2 /2 and (∂P 2 /∂s 2 ) n2 = T /2. Moreover, we have used integration by parts to simplify the expressions: for example, the contributions from the middle two terms in Eq. (50) with V ext (r ⊥ ) can be shown to cancel with each other with our polynomial ansatz. For a given value of µ 0 /k B T (or T /T F , see Eq.(12)), the weighted mass moments and spring constants can be calculated by using local thermodynamic variables in Eqs. (14), (17), (18), (20) and (25). The detailed expressions for numerical calculations are listed in Appendix B.
It is easy to see that the minimization of the action S is equivalent to solving or det S(ω) = 0. The detailed numerical procedure is presented in Appendix C. Once a solution (i.e., the mode frequency ω and the coefficient vector Λ) is found, we calculate the density fluctuation of the mode, using We have performed numerical calculations for the number of the variational parameter N p up to 10, for any given chemical potential µ 0 /k B T or temperature T /T F . In the following, we first discuss the decoupled first and second sound. Then, we focus on the effect of the coupling between first and second sound and the density fluctuation of second sound modes.

B. Decoupled first and second sound
In Figs. 6(a) and 6(b), we present mode frequencies of decoupled first and second sound modes of a highly oblate unitary Fermi gas, obtained by solving the individual action S (a) and S (e) , respectively.
For all the first sound modes, except the lowest-lying one, the mode frequency decreases monotonically with increasing temperature, exhibiting the same temperature dependence as the first sound modes in an isotropic [21] or highly elongated harmonic trap [16,26]. The lowest-lying breathing mode instead does not depend on the temperature and takes an invariant mode frequency ω B = √ 3ω ⊥ . We note that such a temperature independence is a peculiarity of the unitary Fermi gas due to its inherent scale invariance. Indeed, our variational ansatz u a (r ⊥ ) = r ⊥ of the lowest-lying breathing mode is one of the few exact scaling solutions exhibited by the Landau two-fluid hydrodynamic equations at unitarity [25]. To understand this, we simply recall that the spring constant K (a) ij satisfies, ij . Taking into account the fact that K (ae) i=0,j = 0 or K (ae) j=0,i = 0, it is readily seen from Eq. (66) that ω 2 = 3ω 2 ⊥ provides an exact solution, regardless the number of the variational ansatz used. We note also that, the first sound solutions converge very quickly with the number of the variational parameters, indicating that these solutions are indeed well-approximated by the polynomial function.
For second sound modes, we find that the mode frequency initially increases with increasing temperature and then drops to zero very dramatically when the tem- perature approaches to the critical value. This behavior differs from the earlier result in a three-dimensional isotropic harmonic trap [21] but agrees with a recent prediction made for a high elongated configuration [26]. Qualitatively, we may estimate the discretized second sound mode frequency by using the expression ω ∼ u 2 q with a characteristic wavevector q ∼ 1/R s , where R s is the size of the superfluid component along the radial direction. Within the local density approximation, we have R s ∼ R F 1 − T /T c . Using Eq. (60) u 2 ∼ (1 − T /T c ) α/2+1/4 , we obtain that Thus, for the critical exponent α > 1/2, the second sound mode frequency ω vanishes at T c . The superfluid density data (i.e., of superfluid helium) that we use have a critical exponent α ≃ 2/3, and consequently, we find the vanishing 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. For comparison, the decoupled first and second sound mode frequencies are also shown, by black lines and red dashed lines, respectively. As anticipated, the exact solution with frequency ω = √ 3ω ⊥ remains unchanged with the inclusion of the coupling term.
Very similar to the sound velocities in the uniform case (see Fig. 4), the first sound mode frequency is barely affected by the coupling term S (ae) . This is particularly evident in Fig. 8(a), in which we show an enlarged view for the n = 2 first sound mode (the integer n labels the n-th first sound modes). The correction to the first sound mode frequency due to the coupling is about 0.5% and occurs only at around T ∼ 0.14T F . On the other hand, the frequency of second sound modes is notably pushed down by the coupling term S (ae) , as can be seen from Fig. 8(b). The maximum correction is up to 20% when the temperature is about 0.18T F . In the uniform case, i.e., Fig. 4, we find a similar correction to the second sound velocity in the same temperature regime.

D. Density fluctuations of sound modes
The sizable correction in the mode frequency of a second sound mode strongly indicates that, its density fluctuation, as a result of the coupling to first sound modes, should also be significant. In Fig. 9(b), we show the density fluctuations of the lowest two first sound modes (m = 2 and 4) and of the lowest two second sound modes (m = 1 and 3), at a temperature T = 0.18T F . Here, we have used the integer m as the index of different modes when the coupling term is taken into account.
Remarkably, the amplitude of the second sound density fluctuations is just a bit smaller than that of the first sound modes, over a wide range of temperatures, revealing that it could be detected experimentally. In this respect, we note that, most recently the density fluctuations of the low-lying first sound modes have been measured in a highly elongated unitary Fermi gas, to a reasonably good precision [16]. Therefore, it is very likely that a low-lying second sound mode could also be observed by looking at its density fluctuation, after proper

excitation.
Comparing the density and superfluid density profiles, shown in Fig. 9(a), we find that the density fluctuation of second sounds is most significant within the superfluid core, in accordance with their temperature-wave nature. In contrast, the density fluctuation of first sound modes extends over the whole Fermi cloud.

E. Dependence on the superfluid density
We not that, a complete theoretical description of the superfluid density for a Fermi gas at unitarity is yet to be determined due to the theoretical difficulties in handling the strong interaction. In the above studies, we have considered the superfluid fraction of superfluid helium for the superfluid density of the unitary Fermi gas, see, for example, Eqs. (27) and (28). As the second sound is arguably the most significant demonstration of superfluid, it is natural to anticipate that the mode frequency of second sound modes should depend very sensitively on the form of superfluid density.
We have performed numerical calculations with other models of the superfluid density, as shown in Fig. 10. The two additional models we considered are: firstly,  [47]. Note that, recently the superfluid fraction of a 3D unitary Fermi gas has been calculated by Salasnich based on the Landau's expression for superfluid density and physical elementary excitations [29]. The predicted temperature dependence is similar to what we have shown in the figure, presumably due to the simialr elementary excitations arising from strong interactions.
(n s /n) 3D = 1 − (T /T c ) 4 , which is in magnitude similar to the superfluid fraction of superfluid helium but with a mean-field-like critical exponent α = 1; secondly, (n s /n) 3D = (1−T /T c ) 2/3 , which differs significantly from the superfluid fraction of superfluid helium but takes into account the correct critical exponent α ≃ 2/3. The predicted mode frequencies for different superfluid density are shown in Fig. 11. It is readily seen that the mode frequency of second sound modes displays a very sensitive dependence on both the magnitude and critical exponent of the superfluid density. In particular, the mode frequency with the superfluid-helium-like superfluid fraction ( Fig. 11(a)) differs significantly from that with the choice (n s /n) 3D = 1 − (T /T c ) 4 , although these two superfluid fractions differ slightly in magnitude. The strong dependence is very encouraging, indicating that practically, the superfluid density of a unitary Fermi gas could be accurately determined by measuring the mode frequency of low-lying second sound modes.

F. Experimental considerations
Ideally, a compressional second sound mode would be excited by modulating the trapping frequency close to the predicted mode frequency. However, in our quasi-2D configuration, this can be problematic due to the fact that the radial confinement is generated by a magnetic field, whose modulation can also change the s-wave scattering length and hence the system may be away from the unitary regime. We note that this problem does not exist in other setups where it is possible to have optical radial (and axial) confinement.
For our experimental setup at Swinburne, therefore, a more practical scheme would be to locally perturb the cloud and allow it to relax. Taking inspiration from the selective excitation schemes, developed by the Innsbruck experiments [17], one may direct an appropriately sized blue-detuned laser beam into the cloud and apply a modulation burst, chosen to provide the best mode matching with the calculated second sound mode. The power, duration and shape of the laser beam should be optimized in order to resonantly derive the desired small amplitude oscillation in the linear response regime.

VI. CONCLUSIONS
In conclusion, we have derived the two-dimensional simplified Landau two-fluid hydrodynamic equations to describe the low-energy dynamics of a unitary Fermi gas confined in highly oblate harmonic traps. By using a variational approach, these two-fluid hydrodynamic equations have been fully solved. We have discussed in detail the resulting density-wave (firsts sound) and temperature-wave (second sound) oscillations, in the absence or presence of a weak transverse harmonic trap. First and second sound velocities or discretized mode frequencies have been predicted, accordingly.
We have found a weak coupling between first and second sound in highly oblate unitary Fermi gas, very similar to the case of superfluid helium. Though the coupling is weak, it induces significant density fluctuations for second sound modes, indicating that second sound could be potentially observed in the highly oblate configuration, by measuring the density fluctuations. Owing to the strong sensitivity of second sound mode frequencies to superfluid density, the experimental measurement of discretized second sound modes could provide a promising way of accurately determining the superfluid density of a unitary Fermi gas. Here, M and K denote collectively the matrix of the dimensionless weighted mass moments and the spring constants, respectively. The matrix M is positively definite, so that we decouple it by a product of a lower triangular matrix L and its transpose, In terms of this decomposition, the matrix equation, Eq. (C1), becomes It is easy to see that the matrix [L −1 · K · (L −1 ) T ] is symmetric and its eigenvalues give rise to the desired solution of mode frequencies. The displacement field for each eigenvalue can also be calculated accordingly, with known eigenstate of the matrix [L −1 · K · (L −1 ) T ]. More explicitly, if we denote an eigenstate by X, the corresponding displacement field is given by Once the displacement fields for a mode are found, we may calculate its density fluctuation and temperature fluctuation.