Higher-Harmonic Collective Modes in a Trapped Gas from Second-Order Hydrodynamics

Utilizing a second-order hydrodynamics formalism, the dispersion relations for the frequencies and damping rates of collective oscillations as well as spatial structure of these modes up to the decapole oscillation in both two- and three- dimensional gas geometries are calculated. In addition to higher-order modes, the formalism also gives rise to purely damped"non-hydrodynamic"modes. We calculate the amplitude of the various modes for both symmetric and asymmetric trap quenches, finding excellent agreement with an exact quantum mechanical calculation. We find that higher-order hydrodynamic modes are more sensitive to the value of shear viscosity, which may be of interest for the precision extraction of transport coefficients in Fermi gas systems.


I. INTRODUCTION
Strongly interacting quantum fluids (SIQFs) such as high T c superconductors [1], clean graphene [2], the quark-gluon plasma [3], and Fermi gases tuned to a Feschbach resonance [4] seem to lack a description in terms of quasi-particle degrees of freedom. This has fueled interest in developing new tools to understand the transport properties of these fluids, as well as trying to experimentally determine those properties more precisely.
As of yet, one of the cleanest experimental realization of SIQFs is a Fermi gas tuned to a Feschbach resonance. Fermi gases offer unprecedented control of a multitude of properties such as interaction strength, system geometry, spin imbalance [5,6], and mass imbalance [7]. In the case of a spin and mass balanced gas, there have been a number experiments aimed at the extraction of shear viscosity [8][9][10][11], and -to a lesser extent -bulk viscosity [12].
The Navier-Stokes equations provide a relatively straightforward model for the dependence of cloud expansion and collective oscillation phenomena on transport coefficients, making them a seemingly ideal candidate for extraction of such coefficients. Yet, in the low density corona of trapped atom gases the local mean free path becomes large, and hence one cannot expect the Navier-Stokes equations to apply. This, as well as uncertainties arising e.g. from trap averaging, gives rise to a large systematic error in transport coefficients thus extracted from experimental data. Hence a theory which can address both the hydrodynamic behavior of the high density region as well as the low density corona of the cloud is desirable. It has been shown that by including extra "non-hydrodynamic" degrees of freedom in a fluid dynamical description, termed anisotropic fluid dynamics, one can obtain a smooth crossover between Navier-Stokes dynamics in the high density core of the gas cloud and kinetic theory in the low density corona [13]. This theory has been recently used to determine the shear viscosity in the high temperature regime with an error of five percent by comparing experimental data for an expanding cloud to an anisotropic hydrodynamic description [14]. Similar precision determinations for transport properties at lower temperatures, e.g. close to the superfluid transition, are still outstanding.
This present work is related to studies using anisotropic hydrodynamics in the sense that we will also employ a hydrodynamic description beyond Navier-Stokes ("second-order hydrodynamics") in order to study collective oscillations of harmonically trapped Fermi gases above the super-fluid transition (T > T c ). In the linear response regime we are considering in this work, it turns out that the second-order and anisotropic hydrodynamic equations of motion are identical. We will refer to our approach as second-order hydrodynamics to simplify the discussion, but the only difference to an anisotropic hydrodynamics framework will be in name.
Our work is closely related to Ref. [15], and in many aspects is complementary to the results therein. In this work, the focus is on the effects arising from a non-vanishing shear viscosity, and we limit our consideration to ideal equation of state, whereas in Ref. [15] collective modes for polytropic equations of state and zero shear viscosity were studied.
The outline of the paper is as follows: We begin by describing our theoretical framework of second-order hydrodynamics in Sec. II. We then proceed to calculate the frequencies, damping rates, and spatial structures for the collective modes of harmonically trapped gases in both two and three dimensions in Secs. III,IV. We calculate mode excitation amplitudes for experimentally relevant conditions in Sec. V and offer our conclusions in Sec. VI. Detailed results on the spatial mode structure and mode amplitude calculations can be found in three appendices.

II. SECOND-ORDER HYDRODYNAMICS
The Navier-Stokes equations are conservation equations for mass, momentum, and energy. To close the system of equations, constitutive relations between the viscous stresses and fluid variables need to supplied. For Navier-Stokes, the viscous stress tensor π ij is set to first-order gradients of the fluid dynamic variables, mass density ρ, flow velocity u, and temperature T . While widely successful in many fluid dynamics applications, such a first order gradient expansion suffers from certain problems, in particular, in systems where the fluid speed approaches the speed of light [16]. Thus, more recently a second-order hydrodynamic framework has been developed which -true to its nameincludes second-order gradients in the expansion of the stress tensor, with appropriate new second-order transport coefficients. Unlike the similar framework of the Burnett equations, second-order hydrodynamics in addition contains a resummation procedure which ensures that it is a consistent, causal and instability-free generalization of Navier-Stokes (cf. the reviews in Refs. [4,17]). For the case of a unitary Fermi gas, scale-invariance seems to be a good symmetry, and the resulting non-relativistic form of the second-order hydrodynamic equations has been derived in Ref. [18].

II.1. Basic Equations
In the following we consider what is maybe the simplest possible second-order hydrodynamics formalism to describe a fluid in d spatial dimensions with a trapping force F. Namely, we utilize a relaxation equation for the stress tensor.
In this case our fluid equations are given by where is the energy density, η is the shear viscosity, and τ π is the relaxation time for the stress tensor. In the above equations, and σ ij are specified in terms of the fluid velocity, mass density and pressure P as Note that Eq. (5) corresponds to the equation of state for a scale-invariant system. It is easy to show that the familiar Navier-Stokes equations are recovered upon taking the limit τ π → 0 in Eq. (4).

II.2. Assumptions
For simplicity, we have assumed the bulk viscosity and heat conductivity coefficients to vanish. The assumption of vanishing bulk viscosity is consistent with measurements in two dimensions. [19,20]. Furthermore, calculations of bulk viscosity in d = 3 imply that the value of bulk viscosity near unitarity in the high temperature limit should be small [21]. Since we will consider a Fermi gas in the normal phase, i.e. above the superfluid transition temperature T c , taking the bulk viscosity to vanish should be a good approximation in the case d = 3 as well. The assumption of vanishing thermal conductivity is justified as it is already a second-order gradient effect as discussed in Ref. [22]. Hence we assume the gas is isothermal, but it is straightforward to see how the procedure below can be extended to the non-isothermal case. As a consequence, the temperature is a function of time only and not of spatial coordinates. In order to obtain analytically tractable results, we additionally make the approximation that the gas may be described with an ideal equation of state: where n is the number density of particles (we let = k B = 1 throughout). The effects of a realistic non-ideal equation of state on collective mode behavior in a viscous fluid typically require numerical treatments such as those presented in Ref. [23]. Moreover, we assume η/P to be constant. While this assumption is not expected to hold in the low density corona, it will allow analytic access to the spatial structure, frequency and damping rates of collective modes using a secondorder hydrodynamics framework. More accurate numerical studies including temperature and density effects on the shear viscosity are left for future work.
Finally, in order to access collective mode behavior of the gas, we will assume small perturbations around a time independent equilibrium state characterized by ρ 0 (x), u 0 = 0, and T 0 which are solutions to Eqs. (1)- (7). Thus, we set ρ = ρ 0 (1 + δρ), u = δu, and T = T 0 + δT with δρ, δu, δT assumed to be small. Working in the frequency domain we have δρ(t, x) = e −iωt δρ(x), with similar expressions holding for δu and δT . To simplify notation, from now on perturbations such as δρ denote quantities where the time dependence has been factored out, unless otherwise stated.

II.4. Configuration Space Expansion
For a harmonic trapping potential with trapping frequency ω ⊥ , the solution for the equilibrium density configuration is given by ]. In the following, we will be using dimensionless units such that all distances are measured in units of (T 0 /(mω 2 ⊥ )) 1/2 , times are measured in units of ω −1 ⊥ , temperatures in units of T 0 , and densities in units of 0 . In these units the equilibrium solution is given by where A 0 is a dimensionless positive number setting the number of particles (cf. the discussion in App. C).
In the absence of a trapping potential, it is usually convenient to perform a spatial Fourier transform of Eqs. (8)-(11) in order to obtain the collective modes of the system. However, here we are interested in a harmonic trapping potential (linear trapping force) which breaks translation symmetry. Thus, it is more convenient to use a different expansion basis for the perturbations. Here we choose to expand perturbations in tensor Hermite polynomials, though any complete basis of linearly independent polynomials will do. The N th order tensor Hermite polynomials in d spatial dimensions are given by the Rodrigues formula [24] where i k ∈ {1, 2, ...d} for k = {1, 2, ..., N }. The tensor Hermite polynomials are orthogonal with respect to a Gaussian weight which makes them particularly useful for the case of a harmonic trapping potential. In particular, they satisfy the orthonormality condition Assuming translational invariance along the z-axis in d = 3 spatial dimensions, the expansion in both d = 2 and d = 3 will involve only the tensor Hermite polynomials for d = 2. Recalling the assumption that the gas is isothermal, the polynomial expansion of perturbations is then given by δT where, in the sum over j, m j is understood to run over all combinations of indices unique up to permutations. For example, if M = 2 the second sum runs over m = {(1, 1), (1, 2), (2, 2)}, while (2, 1) is excluded. The reason for this restriction is that H (M ) m (x) is fully symmetric in the indices as can be seen from Eq. (12). One should also note that in Eqs. (14), b (M ) m is used as shorthand for the polynomial coefficients of all components of δu, and for a given M and m is a column vector with d components.
Let us now discuss the details of accessing the collective modes whose spatial structure is associated with polynomials of low degree ("low-lying modes"). Substituting Eqs. (14) truncated at polynomial order N into the linearized secondorder hydrodynamics equations and taking projections onto different tensor Hermite polynomials of order K ≤ N we obtain a matrix equation for the polynomial coefficients in Eqs. (14). The (complex) collective mode frequenciesω are then obtained from requiring a non-trivial null-space of this matrix, and subsequently the spatial structures are obtained from the corresponding null-vectors.

III. COLLECTIVE MODE SOLUTIONS IN d = 2
Results for the density and velocity of low-lying collective modes in d = 2 are shown in Fig. 1. In particular, we find a breathing (monopole) mode which corresponds to a cylindrically symmetric oscillatory change in cloud volume, a sloshing (or dipole) mode where the center of mass of the cloud oscillates about the trap center, a quadrupole mode which is elliptical in shape, and higher-order modes corresponding to higher-order geometric shapes. Note that the spatial structure of these collective modes are similar to those reported in Ref. [15]. More detailed information about the d = 2 collective modes can be found in App. A.
FIG. 1. Time snap shots of density profiles and subsequent momentum density (ρu) for the oscillatory modes in d = 2. Note that the center of the monopole mode is at a lower density than the centers of the other modes since it is volume changing and has a larger radius than the equilibrium configuration. The damping rate of higher-order modes is more sensitive to η/P as discussed in the text. Also note that non-hydrodynamic modes share the same spatial structure as their hydrodynamic counterpart.
The collective mode frequencies ω and damping rates Γ are given as the real and imaginary parts of roots of polynomials, which generally do not admit simple closed form expressions. Hence, in Tab. I we choose to report expressions for the complex frequencies and spatial mode structure from second-order hydrodynamics for the lowlying modes in the hydrodynamic limit η/P 1 and τ π 1 (assuming that τ π and η/P are of the same order of magnitude), in which case simple analytic expressions can be obtained. In addition to the modes shown in Fig. 1 there are three modes in Tab. I which have zero complex frequency. The first corresponds to a change in total particle number, the second corresponds to a change in temperature and width of the cloud, and the third "zero mode" is simply a rotation of the fluid about the central axis. While they are required for the mode amplitude analysis (see Sec. V), the role of the first two of these zero frequency modes is relatively uninteresting. Hence, we relegate detailed discussion of these modes to App. C.
The rows of Tab. I starting with the number mode and ending with the decapole mode are all hydrodynamic modes. We note that at order O(η/P ) the results for these modes match those from an analysis of the mode frequencies of the Navier-Stokes equations at the same order. However, for values of η/P where corrections to the hydrodynamic limit become significant, the frequencies found from the Navier-Stokes equations and second-order hydrodynamics disagree. In Fig. 2 we show the full dependence of the hydrodynamic mode frequencies and damping rates on η/P (assuming τ π = η/P based on kinetic theory [22,25,26]). Note that the result of second-order hydrodynamics for the quadrupole mode exactly matches the result from kinetic theory when setting τ π = τ R = η/P [23,27]. I. Frequencies and damping rates in d = 2 from linearized second-order hydrodynamics assuming η P , τπ 1. The hydrodynamic mode damping rates depend on η/P times a prefactor which increases with mode order. Note that for d = 2 there is no non-hydrodynamic sloshing or breathing mode.
Furthermore, results shown in in Tab. I demonstrate that the hydrodynamic mode damping rates depend on η/P times a prefactor which increases with mode order. This is completely analogous to what has been observed in experiments on relativistic ion collisions, where simultaneous measurements of multiple modes have been used to obtain strong constraints on the value of η/s, cf. Ref. [28]. While higher-order modes have not yet been studied in experiment, it is conceivable that measuring their damping rates could lead to a similarly strong experimental constraint on shear viscosity in the unitary Fermi gas. We are not aware of this approach having been suggested elsewhere in the literature. When aiming for using higher-order modes to analyze shear viscosity in Fermi gases we recall that the present analysis is based on a linear response treatment. Quantitative analysis of higher-order flows will, however, require the inclusion of nonlinear effects, especially for analysis of flows beyond hexapolar order due to mode mixing. For this reason, we suggest the hexapolar mode as a prime candidate for the use of higher-order modes to extract shear viscosity.
Finally, Tab. I also indicates the presence of non-hydrodynamic modes (e.g. modes not present in a Navier-Stokes description). The physics of non-hydrodynamic modes is largely unexplored (cf. Refs. [29,30] for a brief discussion of the topic in the context of cold quantum gases). Results shown in Tab. I imply that several such non-hydrodynamic modes exist, all of which are purely damped in second-order hydrodynamics. The non-hydrodynamic mode damping rates are sensitive to τ π and η/P . Thus the value of τ π could be extracted experimentally by measuring any of the non-hydrodynamic mode damping rates in combination with a hydrodynamic mode damping rate required to determine η P . In Fig. 3, non-hydrodynamic damping rates are shown as a function of η/P when setting τ π = η/P .  3. Two-dimensional non-hydrodynamic collective mode damping rates Γ as a function of η/P (using τπ = η/P ). Subscripts denote mode name (quadrupole "Q"; hexapole "H"; octupole "O"; decapole "D"). Dotted line is Γ nh = 1/τπ (cf. Ref. [29]). Note that in d = 2 there is no non-hydrodynamic sloshing or breathing mode.

IV. COLLECTIVE MODE SOLUTIONS IN d = 3
In the case of a three-dimensional gas in a harmonic trap with trapping frequencies ω z ω x = ω y , the resulting gas cloud takes on an elongated cigar shaped geometry. For ω z = 0, the configuration space expansion Eqs. (14) can be applied because there is no dependence on the coordinate z if we assume a translationally invariant system along the z-axis. In this case, the collective mode structures in d = 3 are qualitatively similar to those obtained in the two-dimensional case, cf. the discussion in Refs. [31,32].
We report results for the low-lying modes in the limit η/P, τ π 1 in Tab. II whereas the full dependence of frequencies and damping rates on η P is shown in Figs. 4,5 for the case τ π = η/P . The only qualitative difference with respect to the d = 2 case is that the breathing mode in d = 3 has a different frequency, a non-zero damping rate, and there is now a non-hydrodynamic breathing mode. See App. B for more details about the spatial structure of the d=3 collective modes.
It should be pointed out that, while second-order hydrodynamics predicts purely damped non-hydrodynamic modes for both d = 2, 3, more general (string-theory-based) calculations suggest that there should be a non-vanishing frequency component in the case of d = 3 [30]. It would be interesting to measure non-hydrodynamic mode frequencies and damping rates in order to describe transport beyond Navier-Stokes on a quantitative level.

V. MODE AMPLITUDES CALCULATIONS
In this section, experimentally relevant scenarios to excite the collective modes of the previous sections are discussed, and the corresponding mode amplitudes are calculated. For simplicity, we assume τ π = η/P in the following. In particular, we focus on studying the excitation of the non-hydrodynamic quadrupole (in d = 2, 3) and non-hydrodynamic breathing (in d = 3) modes, leaving a study of higher-order modes for future work. For simplicity, only simple trap quenches (rapid changes in trap configuration) are considered. We will assume the gas cloud to start in an equilibrium configuration of a (possibly biaxial, i.e. ω x, init = ω y, init ) harmonic trap. At some initial time, a rapid quench will bring the trap configuration into a final harmonic form, which is assumed to be isotropic in the x-y plane with trapping frequency ω x, f inal = ω y, f inal = 1 in our units.
In the case of Navier-Stokes equations, initial conditions are fully specified through the initial density ρ init , velocity u init , and temperature T init or appropriate time derivatives of such quantities. However, second-order hydrodynamics treats the stress tensor π ij as a hydrodynamic variable, so, in addition, an initial condition π ij, init or its time derivative needs to be specified.
For equilibrium initial conditions of a general biaxial harmonic trap with trapping force given by F = −γ x x − γ y y  TABLE II. Frequencies and damping rates in d = 3 from linearized second-order hydrodynamics assuming η P , τπ 1. The hydrodynamic mode damping rates depend on η/P times a prefactor with increases with mode order. Note that there is no non-hydrodynamic sloshing mode, but, unlike for d = 2, there is a non-hydrodynamic breathing mode for d = 3.
we have where T init also needs to be specified. Initial equilibrium implies the condition γ i = ω 2 i, init = σ i /T init so that the cloud width is fully specified once γ i for i = x, y and T init are fixed. In addition, equilibrium of the initial trap allows us to take π ij, init = 0. The mode amplitudes can then be obtained by projecting initial conditions onto the collective modes found in the preceding sections (see App. C for details of the calculation).
Isotropic Trap Quench in d = 2 We first consider the case of an isotropic trap quench γ x = γ y ≡ γ in d = 2 and assume A i /A 0 = 1 and T init = 1 for simplicity. Although this case does not exhibit non-hydrodynamic or higher-order collective mode excitation, it does allow us to make direct comparison to results from the literature for the breathing mode excitation amplitude. This type of initial condition corresponds to a rotationally symmetric trap quench with no initial fluid angular momentum. Symmetry then implies that only the number, temperature, and breathing modes Three-dimensional non-hydrodynamic collective mode damping rates Γ as a function of η/P (using τπ = η/P ). Subscripts denote mode name (monopole "B"; quadrupole "Q"; hexapole "H"; octupole "O"; decapole "D"). Dotted line is Γ nh = 1/τπ (cf. Ref. [29]). Note that in d = 2 there is no non-hydrodynamic sloshing mode, but there is a non-hydrodynamic breathing mode.
can be excited (cf. Tab. I), and the initial amplitude for these modes are readily calculated. Fig. 6 displays the (dimensionless) breathing mode amplitude as a function of the quench strength γ. (Note that the amplitude of the temperature mode is identical to the breathing mode amplitude in this case.) The number mode is not excited since the number of atoms taken in the initial condition match the number of atoms we assumed in our final trap equilibrium (A i /A 0 = 1). The amplitude of the breathing mode for the isotropic trap quench is compared to the results from an exact quantum mechanical scaling solution by Moroz [33] in Fig. 6. As can be seen from this figure, there is exact agreement between the calculations for all strength values γ. Note that the amplitudes in this case are independent of η/P since for d = 2, the breathing mode does not couple to the shear stress tensor π ij .  [33]. Right: Absolute value of the (dimensionless) breathing ("B"), hydrodynamic ("Q h ") and non-hydrodynamic ("Q nh ") quadrupole mode amplitudes as a function the quench strength parameter γy for an anisotropic trap quench in d=2. Results shown are for η P = 0.5. Note that the temperature mode amplitude (not shown) matches the breathing mode amplitude for both the isotropic and anisotropic trap quench in d=2.
Anisotropic Trap Quench in d = 2 We perform a similar analysis to that above, considering the case A i /A 0 = 1, and T init = 1, but now taking γ x γ y = 1, which corresponds to an anisotropic trap quench. The mode amplitudes in this case depend on the value of η/P . In this case, in particular the temperature, breathing and quadrupole modes are excited. Fig. 6 shows the absolute value of the mode amplitudes for the hydrodynamic breathing and quadrupole modes, as well as the non-hydrodynamic quadrupole mode as a function of the quench strength γ y . Not surprisingly, Fig. 6 shows that the anisotropic trap quench gives rise to a considerably larger quadrupole mode amplitude (both hydrodynamic and non-hydrodynamic) than the amplitude of the breathing mode.
For a potential experimental observation of the non-hydrodynamic quadrupole mode, it is interesting to consider the relative amplitude of this non-hydrodynamic mode to the (readily observable) hydrodynamic quadrupole mode. The (absolute) amplitude ratio calculated using the above anisotropic trap quench initial condition is plotted in Fig. 7 as a function of η/P . One finds that the non-hydrodynamic mode amplitude is monotonically increasing as a function of η/P . This is plausible given that for small viscosities one expects the hydrodynamic mode to be dominant, whereas one expects the non-hydrodynamic mode to dominate in the ballistic η/P → ∞ limit.
The present calculation is compared to mode amplitude ratios extracted from experimental data [20] in Ref. [29]. To compare non-hydrodynamic damping rate data and theory, we follow the procedure used in Ref. [29] by employing the approximate kinetic theory relation where K ≈ 0.12 in order to relate the experimentally determined k F a to η/P (see discussion in Refs. [23,29] and references therein for more details on this relation). Using this procedure, one observes qualitative agreement of the amplitude ratios between calculation and experimental data in Fig. 7 (left panel). In addition, one can compare the non-hydrodynamic quadrupole mode damping rate, finding reasonable agreement (cf. right panel of Fig. 7). Right: non-hydrodynamic quadrupole mode damping rate Γ. Experimental data is from the reanalysis done in Ref. [29].
There are several possible reasons why a quantitative agreement between second-order hydro and experiment in Fig. 7 should not be expected. For instance, the present theory calculations neglect the presence of a pseudogap phase and pairing correlations (see e.g. Refs. [27,[34][35][36][37] on this topic). Furthermore, it is likely that the quantitative disagreement in Fig. 7 is at least in part due to the assumptions discussed in Sec. II, such as small perturbations, constant η/P and ideal equation of state. In particular, for the strongly interacting quasi-two-dimensional Fermi gas near the pseudogap temperature T * , significant modification of the equation of state have been predicted and observed, cf. Ref. [34]. Studies aiming for achieving a quantitative agreement most likely will have to rely on full numerical solutions, such as e.g. those discussed in Refs. [14,23], which we leave for future work. In addition, our framework only admits a single non-hydrodynamic mode for each of the collective modes. While this may be appropriate in the kinetic theory regime, other approaches such as that of Ref. [30] indicate that our model may be too simple to capture quantitative features of early time dynamics. Finally, we note that the data shown in Fig. 7 was extracted from experiments that were not designed with the purpose of considering early time dynamics. This may contribute to the large uncertainty of the existing data, as well as possibly introducing significant systematic error.
Isotropic Trap Quench in d = 3 For the case of an isotropic trap quench (A i /A 0 = 1, T init = 1) in d=3, results for the (hydrodynamic and non-hydrodynamic) breathing mode and temperature mode amplitudes are shown in Fig. 8. Unlike the case of d=2, the three-dimensional geometry is capable of supporting a non-hydrodynamic breathing mode; furthermore, we find the temperature mode amplitude to differ from the (hydrodynamic) breathing mode amplitude.
The right panel of Fig. 8 shows the ratio of non-hydrodynamic to hydrodynamic breathing mode amplitudes, which reaches up to about 20% for large enough η/P . It is interesting to note that the apparent saturation of this ratio at about 20% is consistent with the amplitude ratio from Ref. [29], extracted from experimental data in Ref. [31]. (Note that in the experiment of Ref. [31], the gas was released from a symmetric trap, allowed to expand for a short period, and then recaptured in a symmetric trapping potential, which is a different protocol than the trap quench considered here. For this reason we do not attempt a direct comparison to mode amplitudes of Ref. [29] in this case.) , non-hydrodynamic breathing ("B nh ") and temperature ("T") modes as a function the quench strength parameter γ for an isotropic trap quench in d=3. Results shown are for η P = 0.5. Right: Ratio of the non-hydrodynamic to hydrodynamic breathing mode amplitudes as a function of η/P for an isotropic trap quench in d=3 (result independent of γ). The maximal ratio of about 20% is roughly consistent with the results of Ref. [29], and suggests the possibility of experimental observation.

VI. CONCLUSION
We have utilized a second-order hydrodynamics framework in order to gain analytic insight into the collective response of Fermi gases in both two-dimensional ("pancake") and three-dimensional ("cigar") trap geometries. Our results demonstrate a number of interesting features which may be reasonably expected to qualitatively describe experimental results. In some cases we even expect quantitative agreement, such as was the case for a Fermi gas undergoing an isotropic trap quench where we found our results to match those obtained from an exact quantum mechanical scaling solution [33].
For instance, our analysis demonstrates that the damping rate of the volume conserving higher-harmonic modes is proportional to the shear viscosity times the harmonic mode number (i.e. the mode winding number w, see App. A). Similar features have been predicted and experimentally observed in the context of relativistic heavy-ion collisions, cf. Refs. [28,38]. While higher-order mode excitations likely result in weaker signal-to-noise ratio, our work suggests a potential experimental avenue towards a precision extraction of shear viscosity in cold Fermi gases by analogy with the relativistic heavy-ion results.
Our study also discusses the presence, damping rates, and expected mode amplitudes of non-hydrodynamic modes in trapped Fermi gases in detail. Recent studies on these non-hydrodynamic modes suggest they could provide information about the presence of quasi-particles in strongly coupled Fermi gases as well as exhibit deep connections between atomic physics and string theory [29,30]. The present results for the non-hydrodynamic mode amplitudes suggest that non-hydrodynamic modes should be accessible by state-of-the-art cold Fermi gas experiments.
There are a number of ways that the results presented here can be generalized and improved. For instance, we plan to study the case of collective modes for uniform density and temperature equilibrium configurations ("Fermi gas in a box"), as well as anisotropic trap geometries in the future. Furthermore, the availability of fully numerical second-order hydrodynamic algorithms [13,23] will allow relaxing the present assumptions of constant η/P , ideal equation of state, and small perturbations around equilibrium in order to aim for a fully quantitative exploration of the collective modes of trapped Fermi gases.

ACKNOWLEDGMENTS
This work was supported in part by the Department of Energy, DOE award No. de-sc0008132. Publication of this article was funded by the University of Colorado Boulder Libraries Open Access Fund. PR would like to thank John Thomas for fruitful discussions as well as the organizers of the workshop on "Non-Equilibrium Physics and Holography" at Oxford University in July 2016 for providing a stimulating environment for many interesting discussions on this topic. WL would like to thank Jasmine Brewer and Roman Chapurin for useful comments regarding clarity of the manuscript.
Appendix A: Spatial Mode Structure for d = 2 In this appendix we collect the spatial mode structure for the low lying modes in d = 2 (see Tab. III). It is interesting to note that, of the modes found, only the temperature mode and breathing mode are associated with a non-zero value of δT . This is because they are the only two modes which change the volume of the cloud, and hence can lead to heating and cooling of the gas. Also note that each irrotational volume conserving mode (here: dipole, quadrupole, hexapole, octupole, and decapole modes) has two independent realizations related by an appropriate coordinate transformation. For example, the quadrupole mode and tilted quadrupole mode are related by a rotation of 45 o . In general, we can assign a mode the winding number w of the associated velocity field u on a circle centered on the origin. This quantity merely counts the number of full rotations made when following a vector around the prescribed circle. For example, w Dip = 0, w Quad = 1, w Hex = 2,... so that the angle of coordinate rotations to get the second independent mode for a given mode is conveniently given by ∆φ = π 2(w + 1) .
Of course, any rotation through an angle in the range ∆θ ∈ (0, π/(w + 1)) will produce an equally valid independent mode, but the angle in Eq. (A1) provides a uniform approach to finding an independent mode from one already found. This provides a connection for example to the approach in Ref. [29] where inequivalent polynomials under rotation up to quadrupole mode were considered in the second-order hydrodynamics framework used here. one axis and isotropic along the other two in d = 3, the breathing mode now couples to shear stresses. Hence there is now an associated non-hydrodynamic breathing mode as well as a difference in the corresponding temperature perturbation associated with the volume change of the cloud. Since there is no associated velocity field for the time independent temperature zero mode, this mode is associated with a vanishing stress tensor, and hence the mode structure is the same as in the d = 2 case. All other modes also exhibit the same spatial and frequency structure.
Non-hydro Quad. Note that the tilted modes denoted Tilted-or T-for short in the table can be found by an appropriate rotation of coordinates.
It should be noted that as can be seen from Tabs.III and IV, not all of the individual perturbations are orthogonal  (the full mode structures are, however, independent). For example, while δρ T emp = δρ B , it is not possible to construct the full mode structure δ T emp = {δρ T emp , δu T emp , δT T emp } as a linear combination of the full mode structure of the other modes. We note that as a result, care should be used when obtaining the system of equations for the amplitudes give by evaluating Eqs. (C1)-(C4) not to miss contributions from all the important modes. We also point out that in general our process for finding mode structure, it is found that mode frequencies come in pairs, one with positive and the other with negative real part. However, in allowing for a complex amplitude and taking the real part in Eqs. (C1)-(C4) we need only consider modes to have positive real part of the frequency. Let us consider a generic isotropic trap quench in d = 2 for multiple initial conditions in order to demonstrate the role of the temperature and number modes. In this case, the amplitudes take on a fairly simple form Note that the above expressions contain both phase and magnitude information, as they may be negative. We will plot phase and amplitude separately below. Additionally, we see from Eqs. (C5)-(C8) several features which should be expected. To explore this, we will break up our analysis into several cases.
This is the case discussed in the main text (see Sec. V).
Case 2: A i /A 0 = 1, T init = 1 For this case we see from Eqs. (C5)-(C8) that the ratio A i /A 0 gives rise to a non-zero amplitude of the number mode, but leaves the location of the zero of the other two modes at γ = 1. This should be expected since this merely means that at γ = 1 we have more (A i /A 0 > 1) or less (A i /A 0 < 1) atoms in the trap than what was assumed in the equilibrium we expanded about. This should make the role of the number mode more clear and is demonstrated for the case A i /A 0 > and T init = 1 in Fig. 9. Particularly, it is only important if for some reason we chose to expand our dynamics about an equilibrium with a different number of particles than given by our initial conditions. Case 3: A i /A 0 = 1, T init = 1 For this case we see from Eqs. (C5)-(C8) that the number mode is not excited, while the value of T init alters the location of the zero for the temperature and breathing modes. This should be expected since at γ = 1 the breathing mode is excited through the temperature difference. Fig. 10 shows the case where T init < 1, demonstrating that the phase of the breathing mode vanishes at γ = 1 while the magnitude is non-zero. A positive amplitude at γ = 1 is expected since the temperature is below its equilibrium value for the given cloud radius so the cloud will reduce its size to try to reach equilibrium.