Collective modes and superfluidity of a two-dimensional ultracold Bose gas

The collective modes of a quantum liquid shape and impact its properties profoundly, including its emergent phenomena such as superfluidity. Here we present how a two-dimensional Bose gas responds to a moving lattice potential. In particular we discuss how the induced heating rate depends on the interaction strength and the temperature. This study is motivated by the recent measurements of Sobirey {\it et al.} arXiv:2005.07607 (2020), for which we provide a quantitative understanding. Going beyond the existing measurements, we demonstrate that this probing method allows to identify first and second sound in quantum liquids. We show that the two sound modes undergo hybridization as a function of interaction strength, which we propose to detect experimentally. This gives a novel insight into the two regimes of Bose gases, defined via the hierarchy of sounds modes.


INTRODUCTION
The emergent phenomena of quantum liquids such as superfluidity and sound modes depend on a multitude of system features, such as interaction strength, dimensionality and temperature. These two classes of phenomena are linked in an intricate manner, as exemplified by the Landau criterion, which predicts dissipationless flow for a perturbation moving with a velocity below a critical velocity. This critical velocity in turn is associated with the creation of elementary excitations [1][2][3][4][5]. While study of superfluids was first motivated by the properties of helium-II, ultracold atoms have expanded the scope of the study of superfluidity by a wide range of trappable quantum liquids including superfluids having reduced dimensionality and tunable interactions. Measurements of the critical velocity have been performed by perturbing ultracold atom clouds with a moving laser beam [6][7][8] or lattice potential [9,10].
Superfluidity in two-dimensional (2D) systems is a particularly intriguing case, due to their critical properties, which differ from higher dimensional systems. Although 2D systems have no long-range order due to increased thermal fluctuations, they can become a superfluid via the Berezinksii-Kosterlitz-Thouless (BKT) mechanism [11][12][13]. The superfluid phase is a quasicondensate characterized by an algebraically decaying phase coherence. Ultracold atom systems provide unprecedented control and tunablity, allowing the study of superfluidity in 2D. This led to the observation of pair condensation [14], phase coherence [15][16][17], critical velocity [7,18], and sound propagation [19,20]. Two sound modes were recently detected in Ref. [21]. However, their hybridization and sublinear dissipation for small perturbation velocities, which we identify in this paper, have not been detected yet.
In this article, we use both classical-field dynamics and analytical estimates to investigate the induced heating rate as a function of the velocity of a moving lattice po- For the weak-coupling regime, for which the velocity of the non-Bogoliubov (NB) mode is larger than the velocity vB of the Bogoliubov (B) mode, the Landau criterion predicts the critical velocity vc = vB, as indicated. (C) Simulated heating rate R(v) for k0/kc = 0.6, where kc = 2.2 µm −1 is the wave vector above which the Bogoliubov dispersion approaches a quadratic momentum dependence. The fit (red curve) yields a sharp increase at vc = 4.0 mm/s, see text. The two maxima of the heating rate correspond to the two modes, where the first maximum is close to vB = 5.8 mm/s. tential. The simulated heating rate shows two maxima corresponding to two sound modes, and sublinear dissipation at low velocities. By changing the periodicity of the lattice potential we examine the heating rate at varying wave vector. We find phononic excitations for low wave vectors and free-particle excitations for high wave vectors. This is in excellent agreement with Bogoliubov theory and the measurements of Ref. [18]. We determine a critical velocity from the sharp onset of dissi-pation and compare it with the measurements for various interactions [18]. Below the critical velocity, we find that the heating rate has a power-law dependence in the velocity and its magnitude increases with the temperature, which is supported by the analytic estimate of the quasicondensed phase. Finally, we determine the two mode velocities from the heating rate and identify their hybridization dependent on interaction and temperature.

System and dynamical response
We simulate bosonic clouds of 6 Li 2 molecules confined to 2D motion in a box potential. This geometry offers tunability of the effective interaction strength, and was used in Ref. [18], and depicted in Fig. 1A. The system is perturbed by a lattice potential moving at a constant velocity v. The unperturbed system is described by the Hamiltonian (1) ψ (ψ † ) is the bosonic annihilation (creation) operator. The 2D interaction parameter is g =gh 2 /m, wherẽ g = √ 8πa s / z is the dimensionless interaction, m the molecular mass, a s the 3D molecular scattering length, and z the harmonic oscillator length in the transverse direction [22]. We consider a rectangular system with dimensions L x × L y = 100 × 100 µm 2 and a density n = 1.2 µm −2 , comparable to the experimental parameters [18]. We solve the dynamics using the classical-field method described in Refs. [23,24], see also Methods. According to this methodology, we replace the operatorŝ ψ in Eq. 1 and in the equations of motion by complex numbers ψ. We sample the initial states from a grandcanonical ensemble with a chemical potential µ and a temperature T via a classical Metropolis algorithm. To obtain the many-body dynamics of the system, we propagate the state using the classical equations of motion. We model the lattice perturbation via the additional term where n(r) is the density at the location r = (x, y). The lattice potential V (r, t) is directed along the x direction: V 0 is the strength, k 0 = 2π/λ l the wavevector, and v the velocity, where we define λ l as the distance between two maxima of the potential. We move this potential for a fixed perturbation time t p of 100 ms. We calculate an ensemble average of the energy change ∆E = H 0 (t p ) − H 0 (0) using Eq. 1 and ψ(r, t). From this change of energy, we determine the heating rate ∆E/t p for various sets of parameters v, k 0 , V 0 ,g, and T /T 0 . Throughout this paper, we will use the temperature T 0 = 2πnh 2 /(mk B D c ), with the critical phase-space density D c = ln(380/g), as the temperature scale, see Refs. [25,26]. This scale gives an estimate of the critical temperature T c , which is the temperature of the BKT transition. We define a dimensionless heating rate R =h∆E/(t p V 2 0 N ), where N is the number of molecules. This heating rate and its velocity, temperature and interaction dependence are central quantities of this paper. As an example, in Fig. 1C we show R(v) for g = 1 and T /T 0 = 0.3. We used V 0 /µ = 0.05 and k 0 /k c = 0.6. µ = gn is the mean-field energy and k c ≡ √ 2/ξ is the wave vector above which the Bogoliubov dispersion approaches a quadratic momentum dependence, with ξ =h/ √ 2mgn being the healing length.
As we depict in Fig. 1C, the heating rate is small at low velocities. Below, we comment on the velocity dependence in this regime in more detail. The heating rate increases rapidly around a velocity which we refer to as the critical velocity of the condensate. As a simple estimate of this velocity, we fit the heating rate below the Bogoliubov velocity v B = gn/m = 5.8 mm/s to the , with A 0 and v c as fitting parameters. The fit gives v c = 4.0 mm/s, as depicted in Fig. 1C.
In addition to the sudden increase of the heating rate, captured by the critical velocity, the heating rate displays two maxima. These two maxima derive from the two excitation branches of the system. In the example shown in Fig. 1C, the maximum at lower velocities is close to the Bogoliubov estimate v B , as shown. We give further evidence for this interpretation below. We refer to this branch as the Bogoliubov (B) mode. Additionally, there is a second maximum at a higher velocity, which we refer to as the non-Bogoliubov (NB) mode. This is consistent with an excitation spectrum sketched in Fig. 1B. Characterizing these two modes further is the second objective of this paper, in addition to the critical velocity.

Analytical estimates
We present two analytical estimates for properties of the heating rate. The first estimate uses Bogoliubov theory, and the second estimate uses the quasicondensate properties of 2D quantum gases. We derive the heating rate perturbatively at second order in the probing term. The Bogoliubov estimate of the heating rate is [27] hω k = k ( k + 2mv 2 B ) is the Bogoliubov dispersion, u k and v k are the Bogoliubov parameters, with (u k + v k ) 2 = Critical velocity. Measurements of the critical velocity (diamonds) of Ref. [18] are compared to the simulations of T /T0 = 0.1 (blue crosses) and 0.3 (red crosses) for various interactions ln(kFa2D). kF is the Fermi wave vector, a2D is the 2D scattering length, and vF is the Fermi velocity. We also show the propagation velocity of a density wave (circles) for T /T0 = 0.1 and the Bogoliubov velocity vB (continuous line).
is the time independent part of the lattice potential in momentum space, and N 0 is the number of condensed atoms. In dimensionless form R ≡h(dE/dt)/(N 0 V 2 0 ) and after simplifying Eq. 4 we obtain where v 0 is v 0 =hk 0 /m. The normalized heating rate R B has a maximum at v max = (v 2 0 /4+v 2 B ) 1/2 . The location of the maximum moves to higher velocities with increasing k 0 or v 0 . We show v max as a red line in Fig. 2, which cap-tures the trend of the measurement and the simulation. We include the thermal damping of the phonon modes by considering a Landau-type damping Γ k = v Γ k, where v Γ is the damping velocity. This results in [27] We show this estimate in Fig. 2C.
To describe the heating rate at low velocities we relate the heating rate to the momentum distribution n k , and arrive at the expression [27] where E k =hv B |k| is the phononic dispersion at long wavelengths. The momentum density n k follows a powerlaw dependence in the wave vector k as where r 0 is the short-range cutoff of the order of ξ and τ is the algebraic scaling exponent. τ increases monotonously from 0 to 1 as the temperature is increased from 0 to T c .
In dimensionless form and for v < v B , we obtain the heating rate of a quasicondensate at low temperatures [27] R qc = π 32 which scales as R qc ∝ τ v 2 for low v and yields a rapid increase for v close to v B . It vanishes, as τ approaches 0.
Low-velocity dependence of the heating rate.

Comparison
In this section, we compare the analytical estimates with the simulation results and the experimental results. In Fig. 2A we show the measurements of Ref. [18]. The measurement was performed at ln(k F a 2D ) = −2.8, where k F = √ 4πn is the Fermi wave vector and a 2D is the 2D scattering length. The maximum of this measured response is close to the phonon velocity for small k 0 /k c and shifts to larger velocities with increasing k 0 /k c . We perform a simulation for the same system parameters, which results in the heating rate shown in Fig. 2B. The heating rate displays the same overall dependence on the lattice wave vector k 0 . In particular, for vanishing k 0 the maximum of the heating rate approaches the sound velocity. Furthermore, the simulation recovers the measurements for intermediate and high wave vectors. We note that for k 0 /k c > 1 the lattice wave vector k 0 approaches the maximum momentum set by the system discretization length, so that the simulation becomes quantitatively unreliable. In Fig. 2C we show the Bogoliubov estimate R B of Eq. 6. We set the value of v Γ /v B = 0.03, that we have obtained numerically by determining the damping of a density wave [27]. The magnitude of the maximum of R B increases with increasing k 0 as in the simulation in Fig. 2B, and provides a good estimate for both the measurement and the simulation. We note that the simulation also displays a faint second branch above the Bogoliubov branch. As we elaborate below, this is due to the non-Bogoliubov mode. This feature becomes more pronounced and realistically detectable for higher temperatures and for an interaction strength near the hybridization interaction, as we discuss below.
Next, we determine the critical velocity for the same value k 0 /k F = 0.3 and for varying interactions as in the experiment. In Fig. 3 we compare the simulation results of the critical velocity v c for T /T 0 = 0.1 and 0.3 with the measurements of Ref. [18]. The measurement agrees with the simulation of T /T 0 = 0.3 for the interaction strengths ln(k F a 2D ) < ∼ −1.5, except for the data point at ln(k F a 2D ) = −3.2. For strong interactions ln(k F a 2D ) > ∼ −1.5, the measurements are closer to the simulations of T /T 0 = 0.1. This might suggest that the experimental results were obtained at temperatures that varied with varying interaction strengths. We also show the simulation results of the sound velocity v s for T /T 0 = 0.1, which is determined from the propagation of a density wave [27]. For all interactions, this sound velocity is slightly below the Bogoliubov velocity v B and above the results of v c . The difference between v s and v c is higher for strong interaction and high temperature, which is due to the broadening of the heating rate.
Finally, we examine the low-velocity behavior of the heating rate at v/v B 1. In Fig. 4 we show the simulated heating rate R(v) for k 0 /k c = 0.4 and various T /T 0 . R(v) shows a power-law dependence at low velocities, visible as a linear dependence on v on a log-log scale. More specifically, we observe a quadratic dependence on v, as supported by fitting with the function f (ṽ) = αṽ 2 , where α is the fitting parameter.ṽ = v/v B is the scaled velocity. We refer to this power-law dependence as sublinear dissipation, because the power-law dependence is sublinear. We note that above the critical temperature, where the bosons form a thermal cloud, the dissipation is linear in v [28]. The quantity α(T ) increases with the temperature, as we show in the inset of Fig. 4. Both the power-law dependence and α(T ) are consistent with the estimate R qc in Eq. 9, where α(T ) is related to the BKT exponent τ . Since the magnitude of the heating rate is small, it would be demanding to extract τ experimentally. However, in principle this result relates the low-velocity heating rate to the quasicondensate scaling exponent.

Two sound modes
In this section we characterize the two sound modes that we observe in the heating rate. In Fig. 5 we show the simulation results of R(v) forg = 1.4, k 0 /k c = 0.4, and various T /T 0 . R(v) displays two maxima. The first maximum corresponds to the B mode and the second maximum to the NB mode. At low temperatures the B mode has a higher magnitude than the NB mode. This changes as the temperature is increased. For high temperatures the NB mode has a higher magnitude than the B mode. This shift of the magnitude is consistent with the shift of spectral weights of the two modes in the dynamic structure factor [29]. Thus, the heating rate provides direct access to the relative amplitude of two modes [30]. The first maximum disappears and turns into a broad background of the diffusive sound mode above the critical temperature T c /T 0 = 0.86, where the value of T c is identified by the critical velocity going to zero [27].
Next, we analyze the interaction dependence by varyingg in the range 0.6 − 3.4 and using the lattice wave vector fixed at k 0 = 1.1 µm −1 . In Figs. 6(A-C) we show R(v) as a function ofg for T /T 0 = 0.1, 0.3, and 0.6. R(v) displays a hybridization of the two sound modes as a function ofg. This hybridization is more visible at intermediate and higher temperatures, due to the larger weight of the non-Bogoliubov mode, which we observed in Fig. 5. For low temperatures, as in Fig. 6A, this hybridization occurs at a high value ofg, which shifts to low values ofg for high temperatures in Figs. 6B and 6C. Similarly, the magnitude of the velocity difference at the avoided crossing increases with increasing temperature. To study this hybridization we determine the two mode velocities by fitting R(v) with a single and a bi-Lorentzian function. In Figs. 6(D-F) we present the interaction dependence of the two velocities v 1/2 obtained from the fit. We indicate the hybridization point at the interaction at which the magnitude of the heating rate of the high-velocity (v 1 ) mode exceeds the magnitude of the heating rate of the low-velocity (v 2 ) mode by a vertical dashed line in Figs. 6E and 6F. We compare v 1/2 with the velocity v s obtained from the propagation of a density wave [27]. We also show the Bogoliubov velocities v B and v B,T = gn s (T )/m determined using either the total density or the numerically obtained superfluid density n s (T ) [27]. The lower velocity v 2 agrees with v s and v B,T for all interactions and all temperatures. The upper velocity v 1 agrees with v B above the hybridization point only. This suggests that for weak interaction and low temperature the system is in the non-hydrodynamic regime, where the lower velocity is described the Bogoliubov (B) approach and the upper mode is the non-Bogoliubov (NB) mode, as we pointed out in Refs. [29,31]. This scenario changes for strong interaction and high temperatures, where the system enters the hydrodynamic regime and a hydrodynamic two-fluid model describes the two sound modes. Here, the upper mode propagates with the total density and the lower mode with the superfluid density [32,33]. DISCUSSION We have determined and discussed the heating rate of a superfluid 2D Bose cloud by perturbing it with a moving lattice potential, using classical-field simulations and analytical estimates. This study is primarily motivated by the experiment reported in Ref. [18]. Indeed, we find, firstly, that the signature of the Bogoliubov mode in the heating rate is well-reproduced in our study, for which we present a numerical result as well as an analytical estimate. Secondly, we show that the critical velocity that emerges in this type of stirring experiment is consistent with the experimental findings of Ref. [18].
However, our study also suggests to broaden the scope of this type of stirring experiment. The results that we report here, elucidate the general behavior of sound excitations in 2D Bose gases, in particular they give access to the two-mode structure of the excitation spectrum, and their interaction and temperature dependence. We show that the heating rate has two maxima, as a function of its velocity and for fixed lattice wave length, corresponding to the two sound modes of the fluid. The relative weight of these modes changes significantly as a function of temperature, and reverts its hierarchy. Furthermore, we show that the two modes undergo hybridization as a function of the interaction strengths. At interaction strengths below the hybridization strength, the non-Bogoliubov mode has a higher velocity than the Bogoliubov mode, whereas for interaction strengths above the hybridization strengths the hierarchy is reversed, which is consistent with the scenario described in [29]. This provides a new insight into the collective modes of Bose gases, which we put forth here, and which we propose to be tested experimentally.

METHODS
To perform numerical simulations we discretize space using a lattice of size N x × N y = 200 × 200 and a discretization length l = 0.5 µm. We note that l is chosen to be smaller than or comparable to the healing length ξ =h/ √ 2mgn and the de Broglie wavelength λ = 2πh 2 /(mk B T ), see Ref. [34]. This results in the cloud size L x × L y = 100 × 100 µm 2 . For the system parameters we use the density n = 1.2 µm −2 , and various interactions and temperatures. The initial ensemble consists of 100 − 1000 states according to the temperature and the interaction. We propagate the state using the equations of motion. For perturbing the cloud, we lin-early turn on the lattice potential over 100 ms and then move the potential at a velocity v and a lattice wave vector k 0 for a fixed perturbation time t p of 100 ms.
Supplemental Materials: Collective modes and superfluidity of a two-dimensional ultracold Bose gas

ANALYTIC HEATING RATE
We determine the heating rate peturbatively by considering a weak perturbation term. To second order, the heating rate is given by [28] H 0 is the unperturbed Hamiltonian.Ĥ s,I is the perturbation term in the interaction picture. In momentum space the perturbation term is described asĤ where V k (t) is the Fourier transform of the potential V (r, t) andn k = q a † q a k+q is the Fourier transform of the density operatorn(r). a k (a † k ) is the bosonic annihilation (creation) operator. V (r, t) is the lattice potential directed along the x direction: V (x, t) = V 0 cos 2 [k 0 (x + vt)/2], where V 0 is the strength, v the velocity, and k 0 the wave vector. We Fourier transform V (x, t) and obtain We transform Eq. S2 to the interaction picture viaĤ s,I (t) = exp(iĤ 0 t)Ĥ s (t) exp(−iĤ 0 t).

Bogoliubov heating rate
Here we derive the heating rate for a condensate at zero temperature. We use the Bogoliubov approximation and obtain the diagonalized Hamiltonian whereb k andb † k are the Bogoliubov operators, andhω k = k ( k + 2mv 2 B ) is the Bogoliubov dispersion. k = h 2 k 2 /(2m) is the free-particle dispersion and v B is the Bogoliubov velocity. We expand the momentum occupation around the condensate mode asn k = √ N 0 (u k +v k )(b † −k +b k ), where N 0 is the number of condensed atoms, and u k and v k are the Bogoliubov parameters, with (u k +v Using Eqs. S4 and S5 we solve the commutator in Eq. S1 and arrive at the heating rate [28] V k is the time-independent part of the lattice potential in Eq. S3. Using the expression of |V k | 2 and ω k (u k + v k ) 2 = hk 2 /(2m), we obtain The delta term in the heating rate gives an onset of dissipation at the velocity v = ω k0,0 /k 0 , where ω k0,0 = k 0 v 2 0 /4 + v 2 B and v 0 =hk 0 /m. In dimensionless form R ≡h(dE/dt)/(N 0 V 2 0 ) we have We extend this result to nonzero temperatures by including the thermal damping of the phonon modes. We consider a Landau-type damping Γ k = v Γ k, where v Γ is the damping velocity. We replace the delta distribution by a Lorentzian distribution, i.e πδ(x) = lim →0 /(x 2 + 2 ). This results in Quasicondensate heating rate To derive the heating rate for a quasicondensate, we consider a Hamiltonian of the form where E k ≡hω k is the excitation spectrum. For the perturbation term, we transform Eq. S2 to the interaction picture by using a k → exp(−iω k t)a k . This results in We now use Eqs. S10 and S11 to calculate the commutator [H s,I (t), H 0 ], which gives Using Eqs. S11 and S12 we calculate the commutator Integrating Eq. S13 over time t 1 yields the heating rate which at long times approaches This result relates the heating rate to the momentum distribution n k . We consider the phononic dispersion at long wave lengths, i.e., E k =hv B |k|. The momentum distribution is given by where n is the real-space density, r 0 is the short-range cutoff of the order of ξ, and τ is the algebraic scaling exponent. We simplify Eq. S15 and obtain with can be solved, giving where x = (1 + τ /4) and 2 F 1 (a, b, c, d) are the hypergeometric functions. We expand I(φ) as I(φ) = 4πṽ + O(τ ). The dimensionless heating rate R =h(dE/dt)/(N V 2 0 ) is This low-temperature estimate scales linearly in τ , while it vanishes for τ approaching zero. At this order the dependence of k 0 and r 0 drops out. For low velocity the heating rate scales quadratically in v as R qc = πτṽ 2 /32. In Fig. S1(a) we show the result of Eq. S17 as a function ofṽ and τ for k 0 r 0 = 2. The dissipation is nonzero below the Bogoliubov velocity and increases with increasing bothṽ and τ . This nonzero dissipation at low velocities is peculiar to 2D superfluids. In Fig. S1(b) we compare the results of Eqs. S17 and S20, which show good agreement for allṽ and all τ . The low-velocity sublinear behavior is also captured well by the estimate R qc = πτṽ 2 /32.

INFLUENCE OF TEMPERATURE AND LATTICE STRENGTH ON THE CRITICAL VELOCITY
As we show in the main text, the critical velocity is smaller for high temperatures. In this section we analyze this thermal reduction systematically by varying the temperature T /T 0 in the range 0.1 − 0.9, where T 0 is the estimate of the critical temperature. We simulate the heating rate R(v) for n = 1.2 µm −2 ,g = 1.6, and k 0 /k c = 0.4. k c = √ 2/ξ is determined by the healing length ξ. In Fig. S2A we show R(v) as a function of T /T 0 for the lattice strength V 0 /µ = 0.01. With increasing temperature, the broadening of the heating rate increases and the onset of rapid increase occurs at a lower velocity. At low temperatures the heating rate primarily shows one maximum close to the Bogoliubov velocity v B = 7.3 mm/s. At high temperatures the heating rate shows two maxima corresponding to the two sound modes. We fit the heating rate below v B to the function fit, we also extrapolate the temperature T /T 0 = 0.86, for which v c is zero. We denote this temperature as the critical temperature. In Fig. S2B we show R(v) as a function of the lattice strength V 0 /µ for T /T 0 = 0.1. µ is the mean-field energy. The lattice strength introduces an additional broadening of the heating rate. This broadening is higher for higher V 0 /µ. We determine the value of the critical velocity v c , which is shown in Fig. S2B. The critical velocity decreases in a non-linear fashion as a function of V 0 /µ.

DETERMINING THE SOUND VELOCITY FROM THE PROPAGATION OF A DENSITY WAVE
To determine the sound velocity we excite a density wave following the method of Refs. [17,24]. We imprint a phase difference on one-half of the system along x direction. This sudden imprint of the phase results in an oscillation of the phase difference between the two subsystems, ∆φ(t), as shown in Fig. S3A. We fit ∆φ(t) with a damped sinusoidal function f (t) = A 0 exp(−Γt) sin(2πf + φ 0 ) to determine the amplitude A 0 , the damping rate Γ, the frequency f , and the phase shift φ 0 . From the fit to the results in Fig. S3A we obtain f = 35.8 Hz and Γ/(2π) = 1.29 Hz. The sound velocity is determined by v s = 2f L x , where L x is the system length in the x direction. This results in v s = 7.16 mm/s and the damping velocity v Γ = 0.26 mm/s, which are v s /v B = 0.98 and v Γ /v B = 0.03, respectively. The reduction of the sound velocity is small for T /T 0 = 0.1, as expected. In Fig. S3B we show the spectrum of ∆φ(t), which yields the peak at the sound frequency and is the same as the one determined from the time evolution in Fig. S3A.
We choose the gradient along x/y directions and calculate the Fourier transform of the current density (j k ) x/y in the x and y directions. This allows us to determine (j * k ) x (j k ) y , which in the limit of k → 0 are approximated by (j * k ) l (j k ) m = k B T m A n s k l k m k 2 + n n δ lm .
n s and n n are the superfluid and the normal fluid density, respectively. A is the system area. By taking a cut along the line k x = k y = k/ √ 2 and by estimating the k = 0 value using a linear fit, we determine the superfluid density based on Eq. S22. In Fig. S4 we show the interaction dependence of the numerically determined superfluid density. While the superfluid density doesn't show a qualitative dependence on the interaction strength for low and intermediate temperatures, it decreases with increasing interaction strength for high temperature.