Mechanism for sound dissipation in a two-dimensional degenerate Fermi gas

We numerically study the transport properties of a two-dimensional Fermi gas in a weakly and strongly interacting regimes, in the range of temperatures close to the transition to a superfluid phase. For that we excite sound waves in a fermionic mixture by using the phase imprinting technique, follow their evolution, and finally determine both their speed and attenuation. Our formalism, originated from a density-functional theory, incorporates thermal fluctuations via the grand canonical ensemble description and with the help of Metropolis algoritm. From numerical simulations we extract temperature dependence of the sound velocity and diffusivity as well as the dependence on the interaction strength. We emphasize the role of virtual vortex–antivortex pairs creation in the process of sound dissipation.


Introduction
Exciting sound waves within a substance and studying their propagation allows for the exploration of its equilibrium and dynamical properties.The characteristics like the compressibility and viscosity become available via measuring the speed of sound waves and their attenuation.
The propagation of sound in harmonically trapped and homogeneous systems of ultracold bosonic and fermionic atomic gases has been studied experimentally by many groups [1][2][3][4][5][6][7][8][9][10][11] .In a recent work on a uniform two-dimensional weakly interacting Bose gas a sound wave, moving with velocity close to the Bogoliubov sound speed, was observed below and above the critical temperature for the Berezinskii-Kosterlitz-Thouless (BKT) transition 7 .This phenomenon has been investigated and verified by various theoretical papers [12][13][14][15][16] .In Ref. 10 , the first and the second sound modes were observed for the first time in 2D Bose gas.Based on temperature dependence of sound speeds the superfluid density of bosonic gas was determined.The superfluid density revealed a universal jump at the critical temperature, in accordance with the BKT theory and as already proved numerically for a uniform two-dimensional atomic Bose gas in Ref. 17 .A single density wave was observed also in interacting atomic Fermi gas at temperatures below the superfluid transition temperature 9 (for a theoretical description see Ref. 18 ).The damping of the sound mode was measured and it was found that the damping is minimal in the strongly interacting regime.For the strongly interacting regime the diffusivity approaches the universal quantum limit which is h/m, where m is the mass of an atom.The results of Ref. 9 demonstrate that the strongly interacting regime distinguishes from the BCS/BEC sides of the BEC-BCS crossover with respect to the dissipation properties.The experimental data presented in 8 show that the diffusivity achieves a universal quantum limit also in a three-dimensional homogeneous atomic Fermi gas at the unitarity and below the superfluid transition temperature.When the temperature gets higher the diffusivity increases monotonically.
In this paper we focus on the sound propagation in two-dimensional Fermi-Fermi mixtures at low temperatures (but still above the superfluid transition temperature) in a weakly and strongly interacting regimes.Several approaches to study out-of-equilibrium bosonic and fermionic quantum many-body dynamics has been already developed, including those based on the classical field [19][20][21] and the two-particle irreducible effective action, the renormalisation-group theory, or the Hartree-Fock-Bogoliubov [22][23][24] approximation.Here, we demonstrate yet another attempt to fermionic systems based on density-functional theory in its time-dependent version, which has regard to the thermal fluctuations (for a density-functional theory version valid in a fermionic superfluid phase see Ref. 25 ).It can be viewed as a continuation of our previous studies of thermal fermionic systems 26,27 .In our approach, after excitation, we can follow sound waves propagation and have an access both to their speed and the damping related properties.The damping turns out to be minimal for the values of the interaction parameter close to zero and the diffusivity gets close to the universal quantum limit.Opposite is true for positive and negative values of the interaction parameter, where the diffusivity becomes very large.
The paper is organized as follows.In Sec. 2 we introduce the model of a two-component two-dimensional degenerate Fermi gas at temperatures above the transition to a superfluid phase.Then (Sec.3) we describe the numerical experiment in which we excite the sound waves in the system, allow them to propagate, and at the end determine their speed and attenuation.In Sec. 4  we reveal the mechanism of sound dissipation which turns out to be related to creation of virtual vortex-antivortex pairs in the gas.Finally, we conclude in Sec. 5.

Method
We start with a simple description of a two-component fermionic mixture (a mixture of atoms being in two internal degrees of freedom as in experiment of Ref. 9 ) in terms of semiclassical distribution functions, f ± p (r), already detailed in Ref. 27 .Here, indices " + " and " − " distinguish components and the equilibrium, at a given temperature T and a chemical potential µ, semiclassical distributions are The particle energy ε ± p (r) at position r includes the kinetic energy as well as the potential energy related to external trapping and interactions.According to the meaning of semiclassical distribution functions the integrals f ± p (r)dp and (p 2 /2m) f ± p (r)dp represent the atomic densities and atomic local motion (intrinsic) energies in both components.We further assume that we consider the case of a tightly-confined two-dimensional gas in a box-shaped trap.The atomic and the local motion energy densities then ) are the extended fugacities, the functional V int (n + , n − ) describes the interaction between components, and f 2 (z) = 2 ∞ 0 x 3 /(z −1 e x 2 + 1) dx is one of the standard functions for fermions 28 .The free energy functional (which, according to Ref. 29 , substitutes the energy functional in the case of nonzero temperatures) of a whole system can be written as The energy depends on atomic densities, n ± (r), and macroscopic velocity fields, v ± (r), of both fermionic gases.The density of a local free energy of a two-dimensional gas is and the second integral in Eq. ( 1) represents the energy of a macroscopic flow.
It is convenient to represent the density and velocity fields by a single complex field ψ ± (r) introduced via inverse Madelung transformation [30][31][32] .The density and velocity fields are related to the pseudo-wave functions ψ ± (r) as n ± = |ψ ± | 2 and v ± = (h/m) ∇φ ± , where φ ± (r) are the phases of complex functions ψ ± (r).Then the density of a macroscopic flow energy can be expressed in the following way and the functional Eq. ( 1) is transformed to The equations of motion corresponding to the functional Eq. ( 4) are found in a usual way as ih (∂ /∂t)ψ ± (r,t) = (δ /δ ψ * ± )F tot (ψ ± , ∇ψ ± ) and are given by While evolving pseudo-wave functions according to Eqs. (5), the extended fugacities z ± (r) are found from the self-consistency condition The interaction between components is treated within the beyond mean-field approach, the lowest-order constrained variational (LOCV) method 26,[33][34][35][36][37] .The interaction energy of a uniform balanced system of density n (per component) is calculated as E int /V = (h 2 /m) n 2 A(η), where A(η) is a function of dimensionless parameter η = ln (k F a 2D ) with k F = √ 4πn being the Fermi momentum -see Supplementary materials for a derivation based on LOCV approximation and a pseudopotential technique for the two-dimensional contact interactions 48 .The numerically determined function A(η) for two lowest energy branches (attractive and repulsive ones) is plotted in Fig. 11 (left frame) for the parameter η scanning the BCS to BEC crossover.Since the relative population of two lowest branches is exp and A 1 (η) (A 2 (η)) determines the pair energy in the lower (higher) branch, it is clear that mostly lower branch is occupied for the low temperatures we consider (i.e., T /T F < 0.2).Now we introduce thermal fluctuations into our description.We propose here a different way than the one used in Ref. 26 .The grand canonical ensemble of pseudo-fields ψ ± fulfilling Eqs. ( 5) can be obtained from the free energy functional Eq. ( 4) by using the Metropolis algorithm [40][41][42] .At a given temperature T and a chemical potential µ, the probability of having a system with N particles and the free energy equal to F(N), calculated within the grand canonical ensemble, is e µN/k B T e −F(N)/k B T .For our system the free energy F(N) is given by Eq. ( 4), where N is the total number of atoms.The fields ψ ± are determined by expanding them in a particular set of basis functions.In our case a uniform spatial grid with a given step ∆ is used (i.e. both fields are expanded in a set of Dirac delta functions).The maximal energy available on a grid with spatial step ∆ equals h2 (π/∆) 2 /2m.This, introduced by discretization, cutoff energy should allow to include all energy-relevant many-body states while generating the grand canonical ensemble.For noninteracting mixture of Fermi gases it is meaningful to take the cutoff energy, due to thermal broadening of the Fermi-Dirac distribution, as E cut = E F + α k B T , where E F = 2π(h 2 /m) n is the Fermi energy and α ≳ 1 (we take α = 5).For interacting gases we include the interaction-related chemical potential and modify the cutoff as E cut → E cut + δV int /δ n.Fig. 12, right frame (details in Suplementary materials) shows the cutoff energy as a function of η = ln (k F a 2D ) for the mixture consisted of ⟨N ± ⟩ = 1500 atoms at the temperature T /T F = 0.2.Fig. 12 (right frame) suggests that our description of fermionic mixture breaks down for negative values of interaction parameter η when the system enters the BEC phase of the crossover and the mixture should be described rather in terms of composite bosons (then all treatments of bosonic particles at nonzero temperatures are at hand [19][20][21] ).However, an alternative way to achieve the proper cutoff energy exists.Since the total chemical potential is µ tot (η) ≈ E F + µ int (η) (and µ int = δV int /δ n, shown in the inset of Fig. 12, right frame), we search for the appropriate grid via generating the grand canonical ensemble of pseudo-fields ψ ± at the condition ⟨N ± ⟩ = 1500 for each interaction strength.The cutoff energy found in this way remains close to that obtained as described earlier for any interaction strength provided the absolute value |E F + α k B T + δV int /δ n| is used as the definition for a cutoff energy.

Sound waves propagation
To excite sound waves we disturb the mixture of fermionic atoms with the protocol similar to that applied in experimental work of Ref. 9 .Namely, we imprint a phase profile on the pseudo-wave functions ψ ± (x, y) → ψ ± (x, y) exp iφ (x) along x direction, where φ (x) = φ 0 exp [−(x − 0.3 L) 2 /σ 2 ] (see Fig. 2).Here, φ 0 determines the strength of the phase disturbance, L is the size of a two-dimensional atomic system, and σ = L/13 gives the width of the region with imprinted phase.This disturbance results in increase in the total energy of the system ∆E = E imp − E ini , which depends on the particular value of φ 0 and on the temperature T .For example, ∆E = 16 h2 /mL 2 for the smallest perturbation φ 0 = 0.27 π, ∆E = 28 h2 /mL 2 for φ 0 = 0.8 π, and ∆E = 50 h2 /mL 2 for the strongest perturbation considered φ 0 = 1.3 π, for T = 0.15 T F and η = 1 (here, E ini = 2151 h2 /mL 2 ).Note that we use a profile as in Fig. 2 with a symmetrical, double phase jumps (we work with periodic boundary conditions and the care about phase continuity should be taken).Both phase jumps convert to the density perturbations and result in an appearance of dark and white soliton-like structures traveling in opposite directions, see Refs. 43,44 In our case the density structures are hidden inside a thermal noise, so they are almost not detectable from density profiles.To overcome this problem we repeat the phase imprinting procedure on many grand canonical microstates (typically 200) and average the density over all realizations.As a result we obtain clear picture of a density response of the system to the initial perturbation.To see density time evolution on a single plot we also integrate two-dimensional density along one spatial coordinate (perpendicular to the direction of sound wave motion) obtaining time-dependent density pictures like in Fig. 3. Typical density evolution after phase imprinting is shown in Fig. 3.One can see two wave forms propagating in the opposite directions with respect to the region where the phase profile has the maximum.In fact, each of this structures is made of a dark and a white quasisoliton formed from opposite slopes of the phase profile (they are close to each other and almost form a single entity).We could change this behavior by making the width of a phase profile wider, but then we would see in fact four distinguishable structures propagating in the system, which differs from what is seen in the experiment of Ref. 9 .Also note that the propagating sound mode does not bounce from the walls since we use periodic boundary condition -instead the signal is propagating from the other side after reaching the "edge".As in the experiment 9 , the initial signal is traveling with constant velocity (we work in a weak disturbance regime) and is damped while propagating, and finally disappears.To analyze a sound signal we determine its velocity (v) and the diffusion coefficient (D).For that we calculate the density imbalance defined as ∆n(t) = (n l − n r )/(n l + n r ), where n l (n r ) is density in the left (right) halves of the potential box, as a function of time, just like in Ref. 9 .In the other way we integrate the density along the transverse direction to the propagation and Fourier decompose as n(x,t) = ⟨n⟩ + ∑ j=±1,±2,...A j (t) exp( j 2πx/L), and look at the time-dependence of the Fourier coefficient A 1 (t) (in fact, A 1 (t)/max[A 1 (t)]) as in Ref. 7 .Then we fit both the density imbalance and the A 1 (t) coefficient to an exponentially damped periodic function of the form: a exp(−Γt/2) sin(ωt + ϕ) + bt + c (linear part of the fit is used to neglect influence of constant drift of a signal), see Fig. 4. We calculate the velocity as v = L ω/(2π) and the diffusion coefficient as D = L 2 Γ/(2π) 29 .Generally, both methods give close results, and only in few cases they differ or fail to fit the data (only matching results are accepted).
Densities presented in Fig. 3 are smooth enough to obtain values of Γ and ω, but still these values change depending on the realization, especially for the highest temperature analyzed.This is why we analyzed three sets of real time evolution of the mixture after the phase imprinting, each including 200 realizations (microstates) from the grand canonical ensemble.This allow us to have more data to present, and add the variance as additional characteristics.Results are summarized in Figs. 5  and 6.Two features of dissipation mechanism are clearly visible.The dissipation increases with temperature.It is true for all considered values of the interaction parameter.Similar observation is reported in 8 , where a three-dimensional Fermi gas is studied at unitarity.In two-dimensional system, however, the diffusivity grows with temperature much faster -at T /T F = 0.2 the diffusivity is already as large as in three-dimensional case for T /T F = 1 (see Fig. 4 in 8 ).The reason for that is explained in the next section.The sound dissipation also depends on the interaction strength.It is the weakest in the region of η close to zero, when the diffusivity approaches the quantum limit of h/m at the lowest temperature.Similar behavior of diffusivity coefficient was observed in experiment 9 performed at temperatures below the superfluid transition temperature.As shown in the next section, the results presented in Figs. 5 and 6 can be understood via the dissipation mechanism based on creation of virtual vortex-antivortex pairs.equations corresponding to Eqs. (5), that read (just for one component) can be solved by assuming small deviation of a density from the equilibrium one (which is constant): n + = n + eq + δ n + .At the low-temperature limit one has Leaving only small quantities of the first order, the continuity equation becomes ∂ /∂t δ n + = −n + eq ∇ • v + .In the weak interaction case this equation can be combined with the equation of motion leading to the following wave equation The expression in the bracket gives the square of the sound velocity.Hence, for the parameters as used for Fig. 6 (here, ⟨N + ⟩ = 1500) the estimation of the sound velocity is (2π h2 /m 2 ) n + eq ≈ 97 h/(mL).

Mechanism for sound dissipation
Strong damping of sound waves for positive and negative values of the interaction parameter is strictly related to the dimensionality of the system we study.In two-dimensional case the spectrum of elementary excitations becomes rich, qualitatively new modes appear with respect to what occurs in three-dimensional systems.Since the vortex energy depends logarithmically on the area occupied by the system, its creation is possible only at high enough temperatures, higher than the critical temperature for the Berezinskii-Kosterlitz-Thouless phase transition.However, for lower temperatures it is still energetically allowed to excite pairs of opposite charge vortices.Our hypothesis is that the strong damping of sound waves traveling in a fermionic mixture is caused by the appearance of vortex-antivortex pairs in each fermionic component.Obviously, the vortices in the normal Fermi gas are not quantized (see Ref. 46 ) as they are in the superfluid phase.It is then expected that all fermions in a component, as described generally by individual orbitals, locally share the phase.Then locally each Eq. ( 5) becomes one to one correspondence to the set of the Hartree-Fock equations for a spin-polarized fermions (see Ref. 47 for a derivation).The phase locally shared by individual fermions is just the local phase of the pseudo field while the density of the pseudo field is the sum of the densities of fermionic orbitals.
In Fig. 7 we show a sequence of time snapshots presenting positions of vortex-antivortex pairs in the "+" component while evolving a particular grand canonical microstate for ⟨N + ⟩ = 1500, η = 3, and at the temperature T /T F = 0.1.From Fig. 7, which demonstrates a typical dynamic behavior of fermionic components at nonzero temperatures, it is clear that vortex pairs appear on a time scale of 10 −6 mL 2 /h (which is a fraction of a microsecond for a typical size of the trapping box equal to L ≈ 40 µm 9 ).Careful analysis of lifetimes of vortex-antivortex pairs leads to the histogram, Fig. 8. Since for low temperatures the vortex-antivortex correlations are straightforward and determined just by the distance between the vortex and antivortex, the evolution of pairs can be easily monitored.By looking at all possible separations from vortex to antivortex at each simulation time step, the vortex-antivortex pairs are those which correspond to the shortest distances (there is no ambiguity in this assignment because for low temperatures the number of vortex-antivortex pairs is small, see Fig. 7).Since the energy of vortex pairs is of the order of k B T , the product of their energy and a lifetime is of the order of 10 −3 h, thus breaking the Heisenberg uncertainty principle.The vortex-antivortex pairs must be then the virtual pairs.
Hence, the attenuation of sound modes observed in the simulations is caused by the scattering on virtual vortex-antivortex pairs.The number of virtual vortex pairs is large for positive and negative η (see Fig. 9) and the diffusivity is of the order of ∼ 10 h/m, see Fig. 5 for temperatures T /T F > 0.1.Going to the strongly interacting regime the number of vortex pairs is significantly decreased and the attenuation coefficient gets smaller, for lower temperatures it takes values about 1 h/m, the universal quantum limit.
Finally, to uniquely indicate the origin of sound dissipation we generated an ensemble of microstates by enforcing constant phase for each pseudo-wave function ψ ± .The resulting thermal fields ψ ± do not feature vortices but are still fluctuating density fields.After exciting a sound wave in a fermionic mixture we find that the sound propagates through a medium almost without damping, thus leading to diffusivity D ≈ 1 h/m close to the quantum limit both for positive and negative values of the interaction parameter, see Fig. 10.Hence, the presence of virtual vortex-antivortex pairs must be responsible for strong dissipation of sound waves.Similar conclusion can be found in Ref. 15 , where the authors discuss the propagation of sound waves in a two-dimensional ultracold bosonic gases by using a so called dynamic Kosterlitz-Thouless theory and deduce that the vortex dynamics is crucial to recover experimental data 7 on damping.

Summary
In summary, we have studied the propagation of sound waves in a weakly and strongly interacting two-dimensional Fermi gas, in the range of temperatures close to the transition to a superfluid phase.We find numerically the dependence of the sound velocity and sound diffusivity on both the temperature and interactions.The sound diffusivity monotonically increases while the temperature is getting higher, independently of the interactions.At constant temperature the damping of sound takes the lowest values in the strongly interacting regime, where at temperatures close to the superfluid transition the diffusivity approaches the universal quantum limit.We identify that scattering on virtual vortex-antivortex pairs is responsible for strong dissipation of sound waves.

Suplementary materials: LOCV approximation in two dimensions
The Hamiltonian of a uniform two-component two-dimensional (2D) degenerate Fermi gas is given by  where V FF is a zero-range pseudopotential that lead to the scattering length in two dimensions, a 2D .The many-body ground state of H is assumed to be given in the form of the Jastrow-Slater variational wave function where |Ψ ± S ⟩ is the Slater determinant wave function of a component consisting of N ± fermions and f (r) is the Jastrow function describing the two-body correlations between interacting fermions.The pair correlation function, which is assumed to be spherically symmetric, is determined variationally by minimizing the average value of the energy ⟨Ψ JS |H|Ψ JS ⟩/⟨Ψ JS |Ψ JS ⟩.Within the LOCV approximation the correlation function fulfills the Schrödinger equation in a space region defined by r < d.The healing length d, which is of the order of an average atomic separation, is determined self-consistently from the conditions f (r > d) = 1 and f ′ (r = d) = 0.For r > d, f (r) tends to unity and the correlations disappear from (9).In the LOCV method the interaction energy E int /N is approximated by where n is the atomic density.Hence the interaction energy is E int /N = ξ (on average there are only correlated pairs).For the scattering state (the second lowest branch) ξ = h2 k 2 /m > 0 and in the noninteracting case the solution of Eq. ( 10) is where γ ≈ 0.577 is the Euler's constant.Since f (r = d) = 1, the constant of proportionality equals To determine the healing length (d) and the interaction energy per particle (∼ k 2 ), the constraint (12) and the condition f ′ (r = d) = 0 have to be used.The second condition yields Set of nonlinear equations ( 12) and ( 17) can be solved numerically giving quantities (k, d) for each value of 2D interaction parameter η = ln (k F a 2D ), where k F = √ 4πn is the Fermi momentum of a two-dimensional single-component Fermi gas of density n.Similar considerations can be repeated for the bound state (the lowest energy branch), just by making a replacement k → iκ (then ξ = −h 2 κ 2 /m < 0).Now, the interaction energy density of each uniform component, calculated within the LOCV approximation, is E ± int /V = ξ n ± = 4π(h 2 /m) n ± n ∓ (k/k ∓ F ) 2 .Hence, the total interaction energy density is 9/12

Figure 1 .
Figure 1.Left frame: A(η), where η = ln (k F a 2D ), calculated within 2D LOCV approximation for two lowest energy branches.The solid black lines are just a guide to the eye.Inset shows the interaction energy per particle in units of the energy per particle of the noninteracting gas, V int /N tot /E FG , where E FG = π(h 2 /m)n.The data compare well with the results of the fixed-node diffusion Monte Carlo calculations of Ref. 49 (see Fig. 1).Right frame: Cutoff energy as a function of η = ln (k F a 2D ) for the mixture consisted of ⟨N ± ⟩ = 1500 atoms at the temperature T /T F = 0.2.Inset shows the chemical potential µ int .

Figure 3 .
Figure 3. Evolution of a density (integrated along y direction) following the phase imprinting.Left (right) picture is an average over 10 (200) realizations (see text).Both frames correspond to the strength of the imprinted phase φ 0 = 0.8 π, the temperature T = 0.1 T F , and η = 0.5.

Figure 4 .
Figure 4. Evolution of the density imbalance (top) and the lowest Fourier mode (bottom) after phase imprintings for η = 0.5 and different temperatures: T /T F = 0.1 and φ 0 /π = 0.27, T /T F = 0.15 and φ 0 /π = 0.8, and T /T F = 0.2 and φ 0 /π = 1.3, from left to right.Each picture is an average over 200 realizations.In each frame the extracted values of the velocity field v (in units of h/(mL)) and the diffusivity D (in units of h/m) are given.

Figure 7 .
Figure 7. Time-snapshots showing positions of virtual vortex-antivortex pairs in "+" component for a particular realization, for ⟨N + ⟩ = 1500, η = 3, and at the temperature T /T F = 0.1.Time is scaled in the unit of mL 2 /h.Appearance and disappearance of vortex pairs on a time scale of 10 −6 mL 2 /h is clearly visible.

Figure 8 .
Figure 8. Histogram of vortex-antivortex pairs lifetime.Each bin shows the number of occurrences of vortex pairs with a given lifetime, in units of τ = 2 × 10 −6 mL 2 /h, within a period of duration of 0.05 mL 2 /h.Main plot shows the first 17 bins (that corresponds to lifetimes < 3.4 × 10 −5 mL 2 /h, and the inset presents occupations of the rest of the bins (< 10 −4 mL 2 /h).This histogram was plotted based on the single realization studied in Fig. 7 but represents a typical outcome for a single real-time evolution of the system.

Figure 10 .
Figure 10.Diffusivity as a function of interaction parameter for the temperature T /T F = 0.2 and for thermal fields containing no vortices.Error bars represent the variance calculated based on three different sets of time evolution after imprinting the phase φ 0 /π = 1.3.The diffusivity, which stays close to the quantum limit, should be contrasted with huge values shown in Fig. 5 (red points).

48 f
where a(k) and b(k) are coefficients set by the boundary conditions and J 0 (kr) and Y 0 (kr) are Bessel functions of the first and second kinds, respectively.The two-dimensional contact interactions can be introduced by imposing the Bethe-Peierls boundary conditions at r = 0 (r) ∝ J 0 (kr) − π 2[γ + ln (k a 2D /2)] Y 0 (kr) ,

Figure 11 .
Figure 11.A(η), where η = ln (k F a 2D ), calculated within 2D LOCV approximation for two lowest energy branches.The solid black lines are just a guide to the eye.Inset shows the interaction energy per particle in units of the energy per particle of the noninteracting gas, V int /N tot /E FG , where E FG = π(h 2 /m)n.The data compare well with the results of the fixed-node diffusion Monte Carlo calculations of Ref.49 (see Fig.1).

Figure 12 .
Figure 12.Cutoff energy as a function of η = ln (k F a 2D ) for the mixture consisted of ⟨N ± ⟩ = 1500 atoms at the temperature T /T F = 0.2.Inset shows the chemical potential µ int .