Dispersion and Damping of Phononic Excitations in Fermi Superﬂuid Gases in 2D

: We calculate the sound velocity and the damping rate of the collective excitations of a 2D fermionic superﬂuid in a non-perturbative manner. Speciﬁcally, we focus on the Anderson–Bogoliubov excitations in the BEC-BCS crossover regime, as these modes have a sound-like dispersion at low momenta. The calculation is performed within the path-integral formalism and the Gaussian pair ﬂuctuation approximation. From the action functional, we obtain the propagator of the collective excitations and determine their dispersion relation by locating the poles of this propagator. We ﬁnd that there is only one kind of collective excitation, which is stable at T = 0 and has a sound velocity of v F / √ 2 for all binding energies, i.e., throughout the BEC-BCS crossover. As the temperature is raised, the sound velocity decreases and the damping rate shows a non-monotonous behavior: after an initial increase, close to the critical temperature T C the damping rate decreases again. In general, higher binding energies provide higher damping rates. Finally, we calculate the response functions and propose that they can be used as another way to determine the sound velocity.


Introduction
Collective excitations in ultra-cold atomic Fermi gases are a subject of intense experimental [1][2][3][4][5][6][7][8] and theoretical [9][10][11][12][13][14] research, because their spectra provide valuable information on the internal states of the atomic system. In recent works [12][13][14], the eigenfrequencies and damping rates of different collective excitations (phonons, pair-breaking "Higgs"modes, Leggett modes) have been calculated in a non-perturbative way within the Gaussian pair fluctuation (GPF) approach based on the GPF effective bosonic action [15][16][17]. These studies are related to fermionic systems in three dimensions. Because the experimental setup for cold atomic gases made it possible to reach highly anisotropic trapping of atoms, an effectively two-dimensional geometry is at present experimentally achievable. It makes theoretical investigations of collective excitations of ultra-cold quantum gases in two dimensions timely and relevant.
The microscopic GPF method is complementary to the frequently used two-fluid hydrodynamic approximation, which is capable of obtaining the sound velocity in both uniform and trapped quantum gases, see, e.g., References [18][19][20], and is effectively used to study collective properties, such as density distribution, discrete collective modes of a trapped gas, and the sound velocity in the BCS-BEC crossover. These two approaches can be applied in parallel which is significant for a comparison and verification of results. Both methods have their own advantages: the hydrodynamic approach can lead to good results for many observables at a low cost, while the microscopic approach allows us to derive both the eigenfrequencies and damping factors of collective modes mutually consistently.
The purpose of this study is to investigate the dispersion relation of collective excitations of a 2D fermionic superfluid without making the assumption that the damping rate is small enough to treat it perturbatively. In this sense we refer to our method as "non-perturbative", to contrast it with previous treatments. Nevertheless, it must be kept in mind that the Gaussian pair fluctuation approach itself still relies on a fluctuation expansion around a self-consistent mean-field. Here, we particularly focus on phononic (Anderson-Bogoliubov) excitations [9][10][11]13], characterized by their sound-like dispersion relationship at low momenta.

Gaussian Pair Fluctuation Approximation
In the path-integral formalism [21] the partition sum of a quantum field ψ x,τ,σ (that depends on position x, imaginary time τ and spin σ) is obtained by summing all possible field configurations, where and each field configuration is weighted by a factor e −S[{ψ x,τ,σ ,ψ x,τ,σ }] , with S[{ψ x,τ,σ , ψ x,τ,σ }] is the action functional of the system under investigation.
Our calculations start with the action functional for an ultra-cold Fermi gas [15][16][17]. Ultra-cold means that the thermal wavelength of the atomic gas is large compared to the range of the interatomic potential, so that it can be approximated by a contact potential: V(r − r ) = gδ(r − r ). Because of the fermionic anti-symmetry requirement, s-wave interactions can only occur between opposite-spin fermions. Hence, the action potential is given by with m the mass of the atoms and β = 1/(k B T) with T the temperature. The fermionic fields ψ x,τ,↑ and ψ x,τ, ↓ are Grassmann variables, and the number of spin-up and spin-down atoms (σ =↑, ↓) is determined by the chemical potentials µ ↑ and µ ↓ . Here, we restrict ourselves to balanced fermi gases (equal amounts of spin-up and spin-down particles), so that there is only one chemical potential µ = µ ↓ = µ ↑ . The strength of the contact potential is determined by the parameter g. For a two-dimensional atomic gas it can be related to the binding energy E b of the atoms, such that [22]: At low temperature and low binding energy 0 < E b < 2E F , the two-dimensional Fermi gas is in the Bardeen-Cooper-Schrieffer regime of superfluidity, whereas for E b > 2E F it is in the Bose-Einstein quasicondensate regime of tightly bound pairs. Here, we choose units such thath, 2m, k F , k B = 1 (so that E F = 1 as well). These are essentially density-based units: since k F = √ 2πn, the total density n equals 1/(2π) in our units. Please note that the Fermi velocity in these units equals v F = (hk F )/m = 2.
When fermionic pairing is present, it is convenient to introduce a bosonic field representing these pairs, by performing the Hubbard-Stratonovich transformation [23,24]. This transformation converts the partition sum to: with and the Nambu spinors η k,n = ψ k,n,↑ ψ k,n,↓ ,η k,n = ψ k,n,↓ ψ k,n,↑ .
Here, ω n = 2πn/β are the fermionic Matsubara frequencies (n ∈ Z) and k = {k x , k y } = (2π/L){n x , n y } (with n x , n y , n z ∈ Z). The bosonic scalar pair fields are∆ q,m and ∆ q,m . When a uniform pair quasicondensate is present, this implies that a macroscopic number of pairs occupy the same pair state: the quantum field can be approximated by a classical field ∆ that is constant in space and imaginary time: With (6), the partition function becomes [25]: where Ω sp is the saddle-point (or mean-field) grand canonical thermodynamic potential, with ξ k = k 2 − µ and E k = ξ 2 k + ∆ 2 . The gap ∆ must be determined self consistently: as a classical field solution, it satisfies the principle of least action. This leads to the gap equation: Please note that in two dimensions, the gap equation can be solved analytically with respect to ∆: This solution must be supplemented with an equation for the chemical potential. Within the saddle-point approach, the particle density is given by The chemical potential is then obtained from setting n sp = 1/(2π).
To go beyond the saddle-point approximation, we must restore the quantum fluctuations, adding them back to the classical field. The fluctuations around the saddle point can further be split into amplitude fluctuations a q,z and phase fluctuations θ q,z . Assuming that the classical field solution still dominates, the fluctuation fields are considered small so the action can be expanded to quadratic order in the fluctuation fields. This leads to a result for the partition sum that is mathematically similar to the result obtained for a superfluid Fermi gas in 3D [12,13]: As the fluctuation action is quadratic in the fields, the matrix M with elements M ++ , ±iM +− , M −− represents the inverse propagator of the bosonic fluctuations. The matrix elements are given by: with the function The bosonic excitations of the fermionic superfluid determine the poles of the propagator M −1 . For fixed q, the excitation energy is found as the real part of the pole and the damping rate as its imaginary part. Here, our goal is to extract the sound velocity and sound damping, so we focus on the long-wavelength behavior of the Anderson-Bogoliubov mode. The long-wavelength limit is found by setting z = uq for the dispersion relation, and only after this substitution taking the limit q → 0. We arrive at the result: where X (E k ) = dX (E k ) /dE k . In the long-wavelength limit, the matrix elements are a function of the parameter u. The sound mode is then found by identifying for which value of the complex u parameter the propagator has a pole. This comes down to solving det M(u) = 0. Using the long-wavelength expansion (17)- (19), the determinant of the inverse fluctuation propagator can be simplified similarly to Reference [13]: where the coefficients for the 2D case are: with v k = ∂E k /∂k = 2kξ k /E k the group velocity of the fermionic (pair-breaking) excitations. The resulting solution u of Equation 20 is interpreted as the complex sound velocity, as in [13]. Specifically, the real part of is the speed of sound (the sound mode energy divided by q), whereas the negative imaginary part equals the damping rate divided by q.

Sound Mode as Pole of the Propagator
In expressions (21)-(23) for the coefficients U, D and C, there are terms of the form The denominator u 2 − v 2 k results in a branch cuts along the real axis for F (u) in the complex u plane. As we need to obtain U, D and C in the lower complex plane to find solutions of (20), it is necessary to analytically continue these functions through the branch cuts. For this we use the same method as in Refs. [12,13], namely the analytic continuation scheme introduced by Nozières [26]. It is based on computing the spectral function along the real axis Then in the lower complex plane, the continuation of ρ provides the desired analytic continuation F c = F (u) + 2πiρ(u). The spectral function ρ(c) exhibits a kink at The analytic continuation of the spectral function to the left of the kink is different from the continuation of the function to the right of v b . This divides the analytically continued functions in two, as it effectively rotates the branch cut to the line Hence, there are two regions in the lower complex u plane: domain I with Re (u) < v b and domain II with Re (u) > v b . As in the three-dimensional case, the presence of the two domains has a physical origin: the group velocity v k corresponds to a Landau critical velocity for emitting fermionic excitations in the superfluid, and v b is the maximal group velocity for fermionic excitations with k < k F . For more details on the continuation scheme (in 3D), we refer to [13]. An example of the analytic continuation is shown in Figure 1, plotting the numerically computed values of det M in the lower complex plane. A complex root is present at the point where the white line (where Im (det M) = 0) and the black line (where Re (det M) = 0) cross. This root corresponds to the complex sound velocity at the chosen parameters. The vertical red dotted line is the aforementioned branch cut between the two domains. The roots were always sought with an imaginary value of less than 2, because for higher values of Im (u) the collective excitations are overdamped.
Changing parameters such as the temperature, we can track how the root of det M moves across the complex plane. Figure 2 shows the paths of the complex sound velocity u when the temperature passes from zero to T c , the critical temperature at which ∆ → 0. Various paths for u in the complex plane are shown, corresponding to five different values for the binding energy. At T = 0 the results for all binding energies coincide at u = v F / √ 2. This result has been obtained using the mean-field equation of state combined with the mean-field gap equation for background parameters (µ, ∆). The obtained zero-temperature sound velocity reproduces the known zero-temperature limit [27,28] which gives the constant sound velocity in 2D at all coupling strengths. This can also be shown analytically by taking the β → ∞ limit in expressions (21)- (23). For non-zero temperatures, however, the sound velocity differs from the zero-temperature limit. Also taking into account fluctuations beyond the mean-field solutions changes the dependence of the sound velocity on the coupling strength even at zero temperature, so that it is not constant [29]. The present work is focused on the temperature dependence of the sound velocity and damping. Taking into account fluctuations is beyond the scope of this study.  As the temperature is increased towards T c the paths ends up on the real axis. In the BCS regime, the complex sound velocity disappears below the real axis before the critical temperature is reached. The path for E b = 3E F differs from the other paths due to the "bump" at the peak that we could not yet interpret. In general, a higher binding energy gives a higher damping rate, which is possible because the system in the BEC regime is more dependent on the temperature. Figures 3 and 4 show the temperature dependence of the sound velocity c = Re(u) and the damping parameter κ = − Im(u), respectively. Please note that κ equals the damping rate divided by q -just as the energy, the damping rate is linear in q in the long-wavelength limit. The sound velocity decreases monotonically from its zero-temperature value of v F / √ 2, but does not reach zero as the temperature reaches the critical temperature.
The general decreasing behavior of the sound velocity of a superfluid Fermi gas when the temperature rises from zero to T c is common both for 2D and 3D Fermi gases [10]. To understand better the dependence of the sound velocity in 2D on the binding energy and temperature, we can remind that the mean-field result for the sound velocity has been re-derived in Reference [28] using simple thermodynamic relation, which in the chosen units reads: In the mean-field approach, the chemical potential at T = 0 is [28] The chemical potential appears then to be independent on the binding energy and gives the aforesaid sound velocity u = v F / √ 2 in the whole range of E b . This is specific for a 2D gas, as distinct from 3D. Moreover, this result coincides with the zero-temperature sound velocity for an ideal Fermi gas in 2D, which is consistent with the BCS limit of a vanishingly small coupling strength. In the far BCS regime k B T c E F , hence the chemical potential and the sound velocity behave asymptotically close to their zero-temperature expressions. This explains why the decrease of the sound velocity as a function of the temperature becomes less expressed for weaker coupling strengths as can be seen from Figure 3. The smaller the damping parameter κ is, the more stable the collective excitations are, so that at T = 0 we have fully stable excitations. As the temperature rises, as expected, κ also rises, and when the temperature continues to rise, κ goes down again. In the BEC regime, this result makes sense, since one must view κ relative to c and if κ tends to zero, c also goes to zero. The result that one gets in the BCS regime on the other hand is unexpected, namely a stabilization of the mode close to T c . This may mean that our description of the system is not accurate at temperatures close to T c and that the present GPF approximate action for the 2D system is should be restricted to temperatures well below T c . To get a better estimate of the regime of validity, we investigate the response function in the next subsection.

Response Function
An alternative way to determine the sound velocity, without relying on the analytical continuation, consists of using instead the response function. For low momentum, the pair field response function is given by [13]: The dominating term here is proportional to the phase-phase element (M −1 ) 2,2 . We determine the phase-phase response function in the long-wavelength limit as Please note that a presence just below the real axis of a pole (i.e., a node of UC − D 2 ) will lead to a peak of χ(c). The location of the peak approximates the real value of the pole, i.e., to the sound velocity. The height of the peak is determined by the residue, and its width is determined by how far in the lower half plane the pole lies. Hence, the value of c for which the response function achieves a maximum c ≡ c resp can be used as an approximation for the sound velocity.
The numerical solutions for c and c resp lie close to each other, as can be seen from Figure 5. An advantage of this second technique is that it is more reliable to determine the sound speeds close to T c , and hence allows us to estimate the region of validity for our results. From Figure 5, it can be seen that the sound speed starts to rise again close to T c in the BCS regime. This phenomenon is not intuitively expected and deviates from that of Figure 3 but can possibly be explained by second sound. The first sound is the sound wave of the density fluctuations. Normally, temperature differences spread through diffusion, but with second sound we have temperature fluctuations that behave like waves. Second sound only occurs in the Bose condensed state, so it is possible that in the BCS regime the first sound dominates and therefore does not disappear at and above the critical temperature, whereas the BEC regime second sound dominates and therefore disappears at T c . Regardless of this interpretation, we can identify the region where c corresponds to c resp as the region where the method of finding the sound velocity from the poles of the bosonic fluctuation propagator is reliable. This turns out to be a rather large fraction of the temperature range, T 0.95 T c . Figure 5. Sound velocity c resp as a function of temperature for E b = 1.5, 1.8, 2, 2.5, 3 (E F ). Inset: this sound velocity for temperatures close to T C .

Conclusions
The aim of this work is to determine the long-wavelength dispersion of Anderson-Bogoliubov excitations for two-dimensional superfluid Fermi gases by locating the poles of the propagator of pair field fluctuations. The fluctuation propagator is in turn obtained via an expansion of the matrix elements of the GPF action up to the second order in powers of q. As was the case for the 3D system, our method of analytical continuation does not need to assume that the poles lie close to the real axis (i.e., that the damping can be tackled perturbatively). We found that at T = 0, the sound velocity is equal to v F / √ 2 independently on the binding energy, being in line with preceding works, and the damping rate is equal to zero. For non-zero temperatures, the sound velocity differs from this known value. As the temperature rises, the sound velocity decreases and the damping increases. Taking into account fluctuations can of course change the dependence of the sound velocity and damping on the binding energy and temperature. For T = 0 fluctuations can be included in the treatment rigorously, their influence on the equation of state of a Fermi gas in 2D at non-zero temperatures is a nontrivial problem, because the superfluid state in two dimensions is a quasicondensate rather than a true condensate. This is a subject of the future investigation.
The response function has been calculated for various temperatures and binding energies. The maxima of the response functions could then be used to determine the sound rates at different temperatures and binding energies in an alternative way, which gives almost the same results as those obtained through complex poles of the fluctuation propagator, except the close vicinity of the transition temperature, where the second method is more reliable. The estimated sound velocities agree with the previous results, and we get a better view at the sound velocity near the critical temperature T c . In the BCS regime, we still have stable collective excitation at T c , while in the BEC regime this is not the case. The reason for this could be that in the BCS regime first sound dominates and keeps existing past T c , while in the BEC regime second sound dominates and disappears past T c .
Recently, the sound mode was studied experimentally for a uniform Fermi superfluid in a box trap [7,8]. The superfluid was excited at a given frequency by vibrating a wall of the container, and the resulting density response was measured to estimate both the dispersion and the damping. A similar technique could well be employed to study the sound mode in the two-dimensional case, given the versatility of trapping potentials that can be achieved. This could allow experimental verification of the results predicted here. The straightforward extension of the present approach to Fermi gases confined to a frequently used parabolic confinement potential is possible for sufficiently large potential wells. In this case, a change of the sound velocity and damping comes completely from a renormalization of the gap and number equations in a trap [18][19][20] with respect to bulk. In the opposite case of a relatively small-radius confinement potential, the hydrodynamic approach used in Reference [20] can effectively determine spectra of discrete low-lying modes. Another promising development of the present study is related to the collective modes in ultra-cold Fermi gases with dipolar interactions intensely investigated recently [30][31][32]. Different interaction potentials, as well as a deformation of the Fermi surface [32] can be straightforwardly incorporated in the GPF approach. Finally, we note that the path-integral technique used here has as a specific advantage over some other methods that it explicitly allows the studying of damped collective modes at non-zero temperatures, relevant to the experiments.