First and second sound in a dilute Bose gas across the BKT transition

We study the propagation of the two sound modes in two-dimensional Bose gases across the Berezinksii-Kosterlitz-Thouless (BKT) transition using classical-field dynamics, which is motivated by recent measurements of Christodoulou et al. Nature 594, 191 (2021). Based on the dynamic structure factor, we identify the two sound modes as the Bogoliubov (B) and the non-Bogoliubov (NB) sound mode below the transition, and as the diffusive and the normal sound mode above the transition. The NB sound mode velocity is higher than the B sound mode velocity, which we refer to as the weak-coupling regime of the sound modes. We excite the sound modes by driving the system as in the experiment and by perturbing the density with a step-pulse perturbation, as a secondary comparison. The driven response depends on the driving strength and results in higher velocities for the B sound mode at high temperatures near the transition, compared to the sound results of the dynamic structure factor and step-pulse excitation. We show that the higher mode velocity has a weak temperature dependence across the transition, which is consistent with the experimental observation.


I. INTRODUCTION
Superfluidity in quantum fluids is in general accompanied by the phenomenon of two sound modes, namely, first and second sound, which is supported by Landau-Tisza's two-fluid hydrodynamic theory [1,2]. The intriguing phenomenon of second sound was first observed in liquid helium and was described as an entropy wave based on the two-fluid hydrodynamic theory [3]. Ultracold atoms expanded the scope of this study by a wide range of trappable quantum liquids including features such as reduced dimensionality and tunable interactions. Two sound modes were measured in a three-dimensional (3D) unitary Fermi gas [4], the BEC-BCS crossover [5], and dilute 2D and 3D Bose gases [6,7].
Contrary to liquid helium, ultracold gases form a wide range of weakly interacting quantum fluids, which undergo an intriguing interplay of sound modes between hydrodynamic and non-hydrodynamic regimes [7][8][9]. Based on the hydrodynamic theory at zero temperature, the second-sound velocity (v 2 = v 1 / √ D) is either below the first-sound velocity v 1 for D = 3 and 2 dimensions or equal to v 1 for D = 1 dimension [10]. Hydrodynamic theory does not support the sound velocity being higher than the Bogoliubov velocity. However, such regimes at low temperatures were predicted in dilute 3D and 2D Bose gases in the weak-coupling regime [11,12].
Sound propagation in 2D quantum fluids is of particular interest because the superfluid density undergoes a universal jump of 4/λ 2 at the Berezinksii-Kosterlitz-Thouless (BKT) transition [13][14][15], where λ is the thermal de Broglie wave length. This has attracted interest to study sound propagation in ultracold 2D quantum gases both experimentally [6,16,17] and theoretically [12,[18][19][20][21][22][23][24][25][26]. In Ref. [12] we discussed and gave numerical evidence for the weak and strong coupling regimes of the sound modes. For the weak-coupling regime, we found that the non-Bogoliubov (NB) sound mode has higher velocity than the Bogoliubov (B) sound mode. We referred to this as a non-hydrodynamic regime. For the strongcoupling regime, we showed that the B sound mode velocity is higher than the NB sound mode velocity. This was consistent with a hydrodynamic scenario. Furthermore, we found that the two sound modes undergo a temperature-dependent hybridization between these two coupling regimes [9]. We note that for a finite-size system the BKT transition manifests itself as a crossover [27][28][29][30], rather than a sudden jump.
Recently, Ref. [6] reported the measurement of the two sound modes in a homogeneous 2D Bose gas of 39 K atoms across the BKT transition. The density response of the driven system is measured as a function of the driving frequency, allowing the detection of both sound modes. The velocity of the lower sound mode decreases with increasing temperature and vanishes above the transition temperature T c , whereas the velocity of the higher sound mode is higher than the Bogoliubov velocity at zero temperature and displays a weak temperature dependence across T c . The measurements are compared to the two-fluid theory in an infinite system [20], which predicts a jump in the velocities at T c and does not describe the temperature dependence of the higher sound mode at T /T c < 0.75.
In this paper, we use classical-field simulations to study the propagation of the sound modes in 2D Bose gases for the experimental parameters of Ref. [6]. The dynamic structure factor (DSF) of the unperturbed cloud shows the Bogoliubov (B) and non-Bogoliubov (NB) sound modes below the transition temperature and the diffusive and normal sound modes above the transition temperature. This allows us to determine the two sound velocities across the BKT transition, independent of an external probe, serving as a benchmark for the results of the density probes. We implement the experimental method to excite the sound modes by driving the system [6]. We find a driving-strength dependent density response and the excitation of two well-resolved sound peaks is observed for strong driving strengths only. As a secondary comparison, we excite the two sound modes using a step-pulse perturbation of the density [5,12]. Finally, we compare the sound velocity results of the driven response, the step-pulse perturbation and the DSF with the measurements [6]. The measured higher-mode velocity within the experimental error agrees with the simulation results at all temperatures, even at temperatures where the hydrodynamic prediction fails. The measured lower velocity shows a shift to higher velocities compared to the results of the DSF and step-pulse perturbation, which is captured by the simulation results of the driven response at high temperatures. This increase due to the nonlinear response of the strong driving partially captures the experimental observations. This paper is organized as follows. In Sec. II we describe the simulation method and the excitation protocols. In Sec. III we calculate the dynamic structure factor to characterize the sound modes. In Sec. IV we present the results of the driven response and the step-pulse perturbation. In Sec. V we compare the simulation and the measurements of the two sound velocities. We conclude in Sec. VI.

II. SYSTEM AND METHODOLOGY
We simulate a bosonic cloud of 39 K atoms confined to 2D motion in a box potential. This geometry was used in Ref. [6]. The system is described by the Hamiltonian (1) ψ (ψ † ) is the bosonic annihilation (creation) operator. g =gh 2 /m is the 2D interaction parameter, withg = √ 8πa s / z being the dimensionless interaction and m the atomic mass. a s is the 3D scattering length and z is the harmonic-oscillator length of the confining potential in the transverse direction. We use the same density n = 3 µm −2 and the sameg = 0.64, as in the experiments [6]. We choose a square box of size L x ×L y = 32×32 µm 2 and various temperatures T /T 0 . We use the temperature T 0 = 2πnh 2 /(mk B D c ), with the critical phasespace density D c = ln(380/g), as the temperature scale, see Refs. [31,32]. This scale gives an estimate of the critical temperature T c of the BKT transition. For the simulations we discretize space on a lattice of size N x ×N y = 64×64 and a discretization length l = 0.5 µm. We note that l is chosen to be smaller than the healing length ξ =h/ √ 2mgn and the thermal de Broglie wave length [33]. We use the classical-field method of Refs. [30,34]. According to this method, we replace the operatorsψ in Eq. 1 and in the equations of motion by complex numbers ψ. We sample the initial states from a grand-canonical ensemble of a chemical potential µ and a temperature T via a classical Metropolis algorithm. We propagate the state using the equations of motion to obtain the many-body dynamics. To excite the sound modes we add the perturbation term where V (r, t) is the perturbation potential that couples to the density n(r, t) at the location r = (x, y) and time t. Within linear response theory, the induced density fluctuation δn(k, ω) is described in terms of the density response function χ nn (k, ω) = δn(k, ω)/V (k, ω), where V (k, ω) is the Fourier transform of V (r, t). This allows us to determine the collective modes of the system. We first implement the experimental method of exciting both sound modes, as in Ref. [6]. We drive the system along the y direction using V (r, t) = V 0 sin(ωt) × (y − L y /2)/L y , where V 0 is the driving strength and ω is the driving frequency. This predominantly excites the longest wavelength sound modes with the wavevector k L = π/L y [35]. For each ω, we calculate the time evolution of the density profile ∆n(y, t) = n(y, t) − n, which is averaged over the ensemble and the x direction. The driving protocol results in center-of-mass oscillations of the cloud, as shown in Fig. 3(a). Assuming that the density fluctuation corresponds to the lowest excitation mode in the box, we fit ∆n(y, t) with the function n(y, t) = b 0 +b(t) sin[π(y −L y /2)/L y ] using b 0 and b(t) as the fitting parameters. From b(t), we obtain the centerof-mass displacement d(t) = (1/L y ) Ly 0 dy yn(y, t) = 2b(t)L y /π 2 . To calculate the steady-state response, we choose the time evolution between 160 − 360 ms, which is fixed for every ω. Fitting d(t) with the function f (t) = [R(ω) sin(ωt) − A(ω) cos(ωt)] exp(−κt) enables us to determine the reactive (R) and absorptive (A) response, where we have included an additional fit parameter κ as global damping rate. Using the Fourier decomposition V (k L , ω) = 4V 0 /π 2 , the A(ω) response yields Imχ nn (k L , ω) = nπ 4 A(ω)/(8V 0 L y ) and thus allows to determine the dynamic structure factor S(k L , ω) [6].
As second method, we employ a step-pulse perturbation of the density, which is motivated by Refs. [5,12].
To perturb the density we use the Gaussian potential is the time-dependent strength and σ is the width. We suddenly turn on and off this potential for a short perturbation time of about 0.5 ms, which excites sound pulses as shown in Fig. 5.

III. DYNAMIC STRUCTURE FACTOR
To characterize the sound modes we calculate the dynamic structure factor (DSF) Excitation spectra of a 2D Bose gas for the density n = 3 µm −2 and the interactiong = 0.64, which are the same as the experimental values in Ref. [6]. Dynamic structure factor S(k, ω) as a function of the wave vector k = ky and the frequency ω is shown for various T /T0 across the BKT transition, where T0 is the estimate of the transition temperature. The white dashed line is the Bogoliubov dispersion determined using the numerically obtained superfluid density ns(T ); see text. where n(k, ω) is the Fourier transform of the density n(r, t) in space and time. We define n(k, ω) as N l = N x N y is the number of lattice sites and T s = 160 ms is the sampling time for the numerical Fourier transform. The DSF gives the overlap of the density degree of freedom with the collective modes. In Fig. 1 we show S(k, ω) as a function of the wave vector k = k y and the frequency ω for various T /T 0 . At low temperatures it primarily shows one excitation branch, while at intermediate temperatures it shows two excitation branches. Above the transition temperature, it displays the diffusive mode at low momenta and the excitation branch of the normal sound mode in a thermal gas. We compare these spectra with the Bogoliubov dispersion hω k = k ( k + 2gn s ), where k = 2J[1 − cos(k y l)] is the free-particle dispersion on the lattice introduced for simulations and J =h 2 /(2ml 2 ) is the tunneling energy.
We numerically determine the superfluid density n s (T ) from the current-current correlations; see Ref. [12]. We show the Bogoliubov prediction in Fig. 1, which agrees with the lower branch at all k for all T /T 0 below the transition temperature. This enables us to identify the lower branch as the Bogoliubov (B) mode and the higher branch as the non-Bogoliubov (NB) mode. At low temperatures the spectral weight is on the B mode, while at intermediate temperatures both B and NB modes are visible. The broadening of the B mode increases with increasing temperature, which occurs due to Landau damping [36]. Above the transition temperature, the B mode transforms into the diffusive mode and the NB mode continuously connects to the normal sound mode.
In experiments, the DSF is measured via the density response χ nn (k, ω) that describes the density fluctuation created by a perturbing potential. Thus, the density response is a useful tool to identify the sound modes using density probes, as we show in Sec. IV. For k B T hω, the density response relates the DSF via S(k, ω) = −k B T Imχ nn (k, ω)/(πnω) [37,38]. The two sound modes are supported by Landau-Tisza's two-fluid hydrodynamic theory, yielding the density response [38] χ nn (k, ω) = nk 2 m which has poles at the velocities v 1 and v 2 corresponding to the two sound modes. v 2 = T s 2 n s /(c v n n ) denotes an additional velocity, where s is the entropy, c v is the heat capacity, n s is the superfluid density, and n n is the normal fluid density. Following Eq. 5 and including linear damping [39], we fit the simulated S(k, ω) with S(ω) = S 1 (ω) + S 2 (ω) for each k = k y , where The amplitudes x 1,2 , mode frequencies ω 1,2 and damping rates Γ 1,2 are the fit parameters. As an example, in 2.8 µm −1 sets the momentum scale, below which the Bogoliubov dispersion has a linear momentum dependence. We fit these spectra with Eq. 6 and show their fits as the continuous lines in Fig. 2. The two-mode feature of the DSF is captured by Eq. 6. At T /T 0 = 0.81, the lower (B) sound peak has higher spectral weight than the higher (NB) sound peak. Above the transition temperature at T /T 0 = 1.22, the B mode vanishes and results in the diffusive mode, while the NB mode becomes the normal sound mode. This numerical observation of the two sound peaks is consistent with the measured spectra [6], which we discuss below. To determine the sound velocities, we perform fits for various k below k ξ and determine the velocities from the linear slope of ω 1,2 /k. When there is mainly one dominant mode at low and high temperatures, we fit with a single function in Eq. 6. The results of the two sound velocities are presented in Fig. 6.

IV. EXCITATION OF DENSITY PULSES
The two sound modes that we find in the dynamic structure factor can be measured by exciting density pulses [4][5][6]. In the following, we present the method of periodic modulation [6] and a potential quench of the local density [5,12].

A. Periodic driving
We first present the method of exciting both sound modes via periodic driving of the center of mass of the cloud, as described in [6] and Sec. II. The driv-ing potential is directed along the y direction and sinusoidally oscillates in time at frequency ω. In Fig.  3(a) we show the resulting time evolution of the density profile ∆n(y, t) for V 0 /µ = 0.8, ω/ω B = 0.74 and T /T 0 = 0.54. µ = gn is the mean-field energy and ω B = v B k L is the Bogoliubov frequency, which results in ω B /(2π) = 34 Hz for v B = gn/m = 2.25 mm/s. The driving protocol excites center-of-mass oscillations of the cloud at ω/ω B , from which we determine the displacement d(t) = 2b(t)L y /π 2 ; see Sec. II. In Fig. 3(b) we show d(t) and the corresponding fit with the function f (t) = [R(ω) sin(ωt) − A(ω) cos(ωt)] exp(−κt), which allows us to determine the reactive (R) and absorptive (A) response. We find that the damping of the oscillations, determined by κ, depends on ω and V 0 of the driving potential. In Fig. 3(c) we show the A(ω) response determined as a function of V 0 /µ and ω/ω B . For low V 0 /µ, A(ω) primarily displays one maximum that corresponds to the B mode. The location of this maximum occurs below the zero-temperature prediction ω/ω B = 1, which is due to the thermal broadening of the phonon modes at nonzero temperatures [12]. For higher V 0 /µ near and above 1, A(ω) shows two maxima corresponding to the B and NB modes. Interestingly, the separation between the peak locations in frequency space increases with increasing V 0 /µ, which is due to the nonlinear response of the system, that sets in for larger perturbations beyond the linear response, as we describe below.
To determine the peak amplitude and the frequency, we fit A(ω) with the function A 1,2 (ω) = ωS 1,2 (ω), based on the relation S(k, ω) ∝ A(ω)/ω, see Sec. III and Eq. 6. In Figs. 3(d) and (e) we plot the frequency and the amplitude of each sound peak, determined via in- dividual fitting, as a function of V 0 /µ. For low driving strengths up to V 0 /µ ∼ 0.6, the B-mode frequency remains qualitatively unchanged and is about ω 1 /ω B ≈ 0.8. In this weak perturbation regime, the B-mode amplitude increases linearly, which is a characteristic of linear response. At higher V 0 /µ, the B-mode amplitude shows nonlinear behavior, where its frequency decreases and drops to ω 1 /ω B = 0.62 for V 0 /µ = 1.5. The higher (NB) sound peak is resolved only above V 0 /µ > ∼ 0.7. Contrary to the B-mode, the NB-mode frequency increases with increasing V 0 /µ, while its amplitude decreases. This reduction of the B-mode frequency and enhancement of the NB-mode frequency derives from the decreasing superfluid and increasing normal fluid density due to stronger probing, respectively. To determine the sound velocities v 1,2 = ω 1,2 /k L , we use the A(ω) response in nonlinear regime at V 0 /µ = 1, where the two sound peaks are well resolved and motivated by the probing regime of Ref. [6]. We note that at this regime the frequencies of the modes are weakly renormalized by the probe, compared to the linear response regime. We show the results of v 1,2 for various T /T 0 in Fig. 6.
We now relate the A(ω) response to the dynamic structure factor (DSF) both below and above the transition temperature. For this, we calculate the A(ω) response using V 0 /µ = 1 at T /T 0 = 0.54 and 1.35. In Fig. 4, we show the dimensionless responseÃ = π 3 mv 2 B A/(8V 0 L y ) and the corresponding DSF S = k B TÃ/(mv 2 B ω) below and above the transition temperature. For T /T 0 = 0.54, A(ω) shows the two sound peaks, which are also visible in the corresponding S(ω). The two-peak structure is captured by the DSF in Eq. 6 and is consistent with the simulated DSF of unperturbed cloud in Fig. 1. Above the transition temperature,Ã(ω) primarily shows the diffusive sound peak and the higher mode is not discernible. This absence of the higher sound peak is in contrast with the measured response [6] as well as the simulated DSF in Fig. 1 and the steppulse excitation in Fig. 5, which show both diffusive and higher sound modes. The reason for this discrepancy can be the decay of the oscillations shown in Fig.  3(b), which is not discernible in the measurements [6]. In Fig. 4(b) we show S(ω) corresponding toÃ(ω) at T /T 0 = 1.35, which results in the diffusive sound mode centered at ω = 0. We fit S(ω) with the Lorentzian f (ω) = x T Γ T /(ω 2 + Γ 2 T ) using x T and Γ T as the fitting parameters. The width of the diffusive mode yields the thermal diffusivity D T = Γ T /k 2 L = (7.2±0.2)h/m, which agrees with the measured D T = (5 ± 2)h/m [6]. Furthermore, our result of D T is above the sound diffusivities measured in strongly interacting 2D Fermi gases [17].

B. Step-pulse perturbation
In this section we demonstrate the method of steppulse perturbation, which excites sound modes by locally perturbing the density. The perturbation sequence is described in Sec. II and the perturbation potential is a Gaussian, for which we use the strength |V 0 |/µ in the range 0.2 − 1.0 and the width σ/ξ = 2.9, where ξ = 0.51 µm is the healing length. In Fig. 5 we show the time evolution of the density profile ∆n(y, t). The perturbation potential was turned on for about 0.5 ms, which excites two (fast and slow) sound pulses below the transition temperature and one sound pulse above the transition temperature. The increased density at the location of the perturbation results in an additional pulse after the perturbation is turned off. Below T /T 0 ∼ 1, the fast and slow pulse initially overlap and then eventually separate into two pulses propagating at different velocities, as shown in Fig. 5. The fast pulse corresponds to the NB mode, whereas the slow pulse represents the B mode. In the long-time evolution the fast pulse bounces off the wall of the box and continues propagating towards the center. At low temperatures, the amplitude of the B mode is higher than the NB mode, whereas at higher temperature T /T 0 = 0.81, the NB mode has higher amplitude than the B mode. This is consistent with the spectral weights of the modes in the dynamic structure factor in Fig. 1. Near the transition temperature T /T 0 ∼ 1, the propagation of the B mode vanishes and results in the diffusive sound mode. Above the transition temperature, the time evolution primarily shows the propagation of the normal sound mode and also diffusive dynamics at the location of the perturbation. To obtain the sound velocities, we fit the density profile with one or two Gaussians to determine the locations of one or two density pulses. The time dependence of the locations gives the sound velocities, which we show in Fig. 6.

V. COMPARISON TO EXPERIMENTS
In Fig. 6, we combine our simulation results of the sound velocities and compare them with the measurements [6]. The temperature dependence of the two mode velocities determined from the dynamic structure factor (DSF) of unperturbed cloud serves as a benchmark for the results based on the density response involving external perturbations. The higher-mode velocity displays a weak temperature dependence at all temperatures and no signature of a jump at the transition. On the other hand, the lower-mode velocity decreases with increasing temperature and vanishes above the transition temperature. We compare this result with the nonzero temperature Bogoliubov estimate v B,T = gn s (T )/m, which is obtained using the numerically determined superfluid density n s (T ); see, for details, Ref. [12]. v B,T agrees with the lower-mode velocity and shows the crossover behavior at the transition, rather than a jump, which is expected for a finite-size system [27][28][29][30].
In Fig. 6 we now compare the DSF results with the sound velocities determined using density probes in the experiments and in our simulations. Overall, the measured higher-mode velocity agrees with the DSF highermode velocity. However, the measured lower velocity is higher than the DSF lower velocity at T /T 0 < 1, except for the measurements at T /T 0 ≥ 1, which show agreement. We note that the measured T c /T 0 = 1 is determined based on the disappearance of lower sound peak [6]. We also present the simulation velocities obtained by imitating the experimental protocol described in Sec. IV A. We used the driving strength V 0 /µ = 1 in line with the used values in the experiments between V 0 /µ = 0.47 − 1 [6]. The simulated higher-mode velocity agrees with the measured higher-mode velocity and shows deviations from the DSF higher-mode velocity at low temperatures. On the other hand, the simulated lower velocity is smaller than the measured lower velocity at T /T 0 < 0.75 and agrees with the measurements at T /T 0 > 0.75. In Fig. 6 we also show the two sound velocities determined via a step-pulse perturbation in Sec. IV B, which agree with the results of the DSF throughout the transition for both the higher and lower velocity.
The lower-velocity results of the DSF and step-pulse perturbation are in agreement with the Bogoliubov estimate v B,T , confirming that the sound velocity is smaller than the measurements and there is smooth crossover occurring near the transition point. The deviation near the transition temperature is reproduced by the simulations of the driven response, which yield similar velocities as in the measurements. These simulations then systematically deviate from the measurements at intermediate and low temperatures, which seems to occur due to a variation in the used value of the driving strength in the measurements and the corresponding change in nonlinear response. Furthermore, the measurement uncertainty of the box length and the density can also affect the magnitude of the measured sound velocities [6].

VI. CONCLUSIONS
We have determined and discussed the propagation of first and second sound in homogeneous 2D Bose gases across the BKT transition using classical-field simulations for the experimental parameters of Ref. [6]. We have identified the two sound modes based on the dynamic structure factor, which are the Bogoliubov (B) and non-Bogoliubov (NB) sound modes below the transition and the diffusive and normal sound modes above the transition. We have excited the sound modes using the experimental method of periodic driving [6] and the method of step-pulse perturbation [5,12]. We have determined the sound velocities from the dynamic structure factor (DSF), the driven response and the step-pulse excitation and compared them with the measurements of Ref. [6]. While the sound results of the DSF and step-pulse excitation show excellent agreement, the results of the driven response show a systematic deviation, compared to the DSF results, due to nonlinear response. If the probing strength is small, the driven response recovers the Bogoliubov mode velocity. However, for small probing strengths, the signal of the non-Bogoliubov mode vanishes. It only becomes measurable for intermediate probing strengths. At these probing strengths, the nonlinear character of the probe influences the frequencies of the measured sound modes. Therefore, this approach only gives an approximate value of the sound velocities, for these probing strengths. Overall, the simulated higher-mode velocity is above the B mode velocity and displays a weak-temperature dependence across the transition, which is in agreement with the corresponding measurements of the higher mode velocity. On the other hand, the measured lower mode velocity is below the simulated velocities of the DSF and step-pulse excitation but agrees with the results of the simulated density response at high temperatures across the transition.
Our results give insight into the temperature dependence of the two sound modes of dilute 2D Bose gases, and the signature of these modes in the measurement techniques of Refs. [6] and [5]. Our results reproduce largely the measurement results of [6], while demonstrating that the strong-driving results are subject to drivinginduced frequency shifts. These shifts might obscure the measurements of these frequencies. Furthermore, generating a signal for the non-Bogoliubov mode requires strong driving, making this probe technique more suitable for a qualitative investigation of the modes, in contrast to the step-pulse technique of [5]. For increasing interactions or densities, these modes undergo hybridization and crossover to the strong-coupling regime occurs [9,12], which warrant further experimental investigations. The two sound modes and their coupling can affect the dynamics, such as the propagation of deterministic vortex colliders [40]. Our results enable the further study of these phenomena, as they provide an in-depth insight into the key probing techniques of the field.