Sound propagation in a two-dimensional Bose gas across the superfluid transition

Motivated by recent experiments in Phys. Rev. Lett. 121, 145301 (2018), we study sound propagation in a two-dimensional (2D) Bose gas across the superfluid-thermal transition using classical field dynamics. Below the transition temperature we find a Bogoliubov and a non-Bogoliubov mode, above it we find the normal sound mode and the diffusive mode, as we determine from the dynamical structure factor. Our simulations of the experimental procedure agree with the measured velocities, and show that below the transition temperature the measurements detect the Bogoliubov mode. Above the transition, they either detect the normal sound mode for low densities or weak interactions, or the diffusive mode for high densities or strong interactions. As a key observation, we discuss the weak coupling regime in which the non-Bogoliubov mode has a higher velocity than the Bogoliubov mode, in contrast to a hydrodynamic scenario. We propose to detect this regime via step-pulse density perturbation, which simultaneously detects both sound modes

Sound modes in 2D Bose gases are of special interest, as an interacting 2D system undergoes a superfluid transition via the Berezinskii-Kosterlitz-Thouless mechanism [28]. At the transition, and in the thermodynamic limit, the superfluid density vanishes with a universal jump of 4/λ 2 , where λ is the de Broglie wavelength. Furthermore, 2D systems exhibit a universal scale invariance: the dimensionless thermodynamic quantities, such as the phase-space density and the entropy, depend only on a single dimensionless parameter µ/k B T or equivalently T /T c , where µ is the mean-field energy, T the temperature, and T c the critical temperature. This is confirmed by Refs. [29,30], where no jump in the thermodynamic quantities is observed. Refs. [31,32] studied the sound modes of 2D quasi-condensates using the two-fluid model, which show a jump at the transition.
Recently, Ref. [33] reported on the measurements of the sound propagation in a uniform 2D Bose gas of 87 Rb atoms across the superfluid-thermal transition. The temperature dependence of the measured sound velocity shows no discernible jump in the crossover regime and a nonzero velocity above the transition. Theoretical studies of this measurement were reported in Refs. [34,35].
In this paper, we investigate sound mode dynamics of a uniform 2D Bose gas of 87 Rb atoms across the superfluidthermal transition using c-field simulations. We determine a sound velocity c by exciting running and standing waves with a weak Gaussian potential. These results show good agreement with the measurements of Ref. [33]. Below T c , the temperature dependence of c is captured by a Bogoliubov estimate that includes the superfluid density at nonzero temperature. Near and above T c , c displays a temperature dependence that depends on the density in a qualitative manner: c increases and decreases for low and high density, respectively. This is also reflected in the dynamic structure factor, showing the density-dependent interplay between two sound modes that we refer to as the Bogoliubov and the non-Bogoliubov mode below T c , and the diffusive and the normal sound mode above T c . The results of c show a breaking of the universal scale invariance at nonzero temperature due to Landau damping. Going beyond the experimental work of Ref. [33], we propose to excite the two modes using a step-pulse density perturbation. The results of the step-pulse excitation across the transition show excellent agreement with the results of the dynamic structure factor, which provides a simultaneous measurement of both sound velocities. This paper is organized as follows. In Sec. II we illustrate the terminology of first and second sound. In Sec. III we describe our simulation method. In Sec. IV we determine sound velocities by exciting running waves. In Sec. V we analyze the scale invariance of the sound velocity. In Sec. VI we compare the running-wave velocity with the standing wave velocity and the Bogoliubov estimate. In Sec. VII we show the dynamic structure factor. In Sec. VIII we excite two sound modes with a step-pulse perturbation, and in Sec. IX we conclude.
FIG. 1. Sketch of the qualitative temperature dependence of first sound (upper curve) and second sound (lower curve), for (a) weak and (b) strong interactions. B labels the Bogoliubov sound mode, NB the non-Bogoliubov mode, D the diffusive mode, and N the normal sound mode. The line color represents the spectral weight of the modes in the dynamic structure factor. For strong interactions the two modes undergo an avoided crossing at a hybridization temperature below the critical temperature, where the hybridization point is indicated by the crossing of the decoupled modes.

II. WEAK AND STRONG COUPLING REGIME
In this paper, we refer to the faster mode as first sound, and the slower mode as second sound. We note that this terminology is inherited from the study of superfluid helium, and that the application of its terminology to cold atom systems could be done in several ways. As we describe in this paper, for cold atom systems, we find two regimes. For weak interactions, or small densities, the temperature dependence of the sound velocities is sketched qualitatively in Fig. 1(a). At low temperatures, one sound mode is well described by the Bogoliubov approximation (B), and one mode that we refer to as a non-Bogoliubov (NB) mode. Here, the non-Bogoliubov mode is the first sound mode, in the sense that it is the faster mode. In Ref. [27] we have given a weak coupling description of this mode as a squeezing mode. For the interaction g → +0, the ratio of the sound mode veloci-ties approaches two. As the temperature is raised above the critical temperature, the NB mode continuously connects to the normal sound mode of a thermal gas. The velocity of the Bogoliubov mode undergoes a universal jump to zero at the critical temperature in the thermodynamic limit, and becomes the diffusive mode. We note that this sudden jump is replaced by a crossover regime for finite systems. The spectral weight of the modes in the dynamic structure factor has been indicated by the line color in the sketch.
For strong interactions the temperature dependence of the mode velocities is sketched in Fig. 1(b). At low temperatures the Bogoliubov mode is the faster mode, which we refer to as first sound in this regime. The non-Bogoliubov mode is the slower mode, and connects to the diffusive mode. The two modes display an avoided crossing at a hybridization temperature below the critical temperature. We give numerical evidence in support of these scenarios below.

III. SIMULATION METHOD
We simulate the dynamics of a 2D quasi-condensate using the c-field method of Ref. [36]. We describe the system with the Hamiltonian (1) ψ andψ † are the bosonic annihilation and creation operator, respectively. The 2D interaction parameter is given by g =gh 2 /m, whereg = √ 8πa s / z is the dimensionless interaction, a s the s-wave scattering length, and z the harmonic oscillator length in the transverse direction. We useg = 0.167, as in Ref. [33].
Inspired by the experimental setup of Ref. [33], we consider a 2D Bose cloud of 87 Rb atoms confined in a rectangular box geometry of dimensions L x × L y = 34 × 39 µm. For the numerical simulations, we discretize space with a lattice of size N x × N y = 68 × 78 and the discretization length l = 0.5 µm, where l is chosen to be smaller than or comparable to the healing length ξ and λ (see Ref. [37]). In our c-field approach, we replace the operatorsψ in Eq. 1 and in the equations of motion by complex numbers ψ. We sample the initial states in a grand-canonical ensemble of chemical potential µ and temperature T via a classical Metropolis algorithm. We propagate this initial state according to the equations of motion. For each trajectory, we calculate the desired observables, and average over the initial thermal ensemble. The density of the atoms is in the range n 2D = 3.0 − 53 µm −2 . For each n 2D , we choose several temperatures across the transition. The critical temperature T c is estimated by the critical phase-space density D c = ln(380/g) [38], which results in T c = 2πn 2Dh 2 /(mk B D c ). To excite sound modes we add the perturbation H ex = drV (r, t)n(r), where n(r) is the density at the location r = (x, y). The excitation potential V (r, t) is given by where V 0 is the time-dependent strength and σ the width. This potential is used along the upper edge of the box at the location y 0 = 36 µm to excite running and standing waves. For all simulations, σ is 5 µm and V 0 is typically in the range V 0 /µ = 0.1 − 0.4, where µ = gn 2D is the mean-field energy. A running wave is excited using the following scheme. We slowly turn on the potential over t on = 200 ms, i.e. V 0 (t on ) = V 0 , wait for t w = 100 ms, and then suddenly turn it off. This excites a sound wave propagating in space along y direction as a function of t, see Fig. 2(a). To excite a standing wave, the following scheme is used. We slowly turn on the potential in the manner described above, and then sinusoidally mod- , where t = t on + t w and f is the modulation frequency. We perform this modulation at various frequencies f . After 1 − 2s excitation time, we analyze the density modulation following Ref. [33]. For each f , we explore one oscillation by recording the density profiles n i (y, f ) at four different times t i , i ∈ {1, 2, 3, 4}, where t i are chosen according to ωt i = ωt 1 + (i − 1)π/2, with ω = 2πf . The amplitude of the standing waves is calculated by the quantity q 2 (y, f ) = q 2 1 (y, f ) + q 2 2 (y, f ), where q 1 (y, f ) = n 3 (y, f )−n 1 (y, f ) and q 2 (y, f ) = n 4 (y, f )−n 2 (y, f ). The squared amplitudes determined at various f are shown in Fig. 5(a). We determine the sound velocities by exciting running and standing waves in Secs. IV and VI, respec-tively.

IV. SOUND PROPAGATION
In this section, we present the results of running-wave excitation for various combinations of n 2D and T /T c . As an example, we first choose n 2D = 29.2 µm −2 and T /T c = 0.37, which is one parameter set used in the experiment. We excite a running wave following the sequence described in Sec. III. In Fig. 2(a) we show the time evolution of the density profile ∆n(y, t) = n(y, t) − n 2D , which is averaged over the x direction and the ensemble. The excited sound wave is indicated by the density dip propagating in space as a function of time. The sound wave travels back and forth between the edges at a constant velocity and forms a triangular pattern. We fit the locations of the sound wave with a triangular-wave function to determine its velocity c. From the fit, we obtain c = 1.47 mm/s, which is in excellent agreement with the measured c = 1.49 mm/s. The simulated c is slightly below the Bogoliubov estimate of the sound velocity at zero temperature c 0 = gn 2D /m = 1.61 mm/s. Furthermore, we examine the damping of the sound mode in Fig. 2(a). We fit the density profile with the function n(y, t) =n + A 1 (t) cos(πy/L y ) to determine the amplitudes A 1 . This function represents the lowestenergy density mode, and the functional form is motivated by the experiments.n is a fitting parameter and represents the average density in y direction. In Fig. 2 we show the extracted amplitudesÃ 1 (t) = A 1 (t)/A 1 (0) as a function of the propagation time t. We fitÃ 1 with an exponentially damped sinusoidal function [33] f (t) = e −Γt/2 Γ/(2ω) sin(ωt) + cos(ωt) to determine the frequency ω and the damping rate Γ. From these, we determine the sound velocity c = ωL y /π and the quality factor Q = 2ω/Γ. ForÃ 1 in Fig. 2(b), the fit yields ω = 119.7 s −1 and Γ = 11.2 s −1 , and we obtain c = 1.47 mm/s and Q = 21.4. The value of c is the same as for the triangular pattern fit, and the high value of Q implies weak damping of the sound mode. As the main origin of the damping of the sound modes, we identify Landau damping as we explain below.
We now consider the three other sets of n 2D and T that are used in the experiment, which are (53 µm −2 , 0.21 T c ), (52 µm −2 , 0.95 T c ), and (11 µm −2 , 1.38 T c ). For each set, we repeat the running-wave excitation and determine ω and Γ, as above. We show the extracted amplitudesÃ 1 in Figs. 2(c), 2(d), and 2(e), respectively. The values of c, Γ, and Q are given in Table I, where we compare them with their corresponding measured values. They are in agreement below T c , while they deviate for the parameter set above T c . We link this deviation to the measurement uncertainty and possibly different values of V 0 between experiment and simulation.
We now analyze the temperature dependence of c across the transition systematically. We choose the three densities n 2D ≈ 3, 12, and 27 µm −2 . We refer to them as low, moderate, and high density, respectively. For each n 2D , we determine c, Γ, and Q at various T /T c , with the running-wave excitation described above. We use the same V 0 ≈ 0.2 µ for all simulations. We show the normalized results of c/c 0 as a function of T /T c in Fig. 3(a). The temperature range includes the superfluid, crossover, and thermal regime. In the superfluid regime, c overall decreases with increasing T . The reduction in c/c 0 is higher for low n 2D as compared to high n 2D . In the crossover and the thermal regime the temperature dependence of c/c 0 depends on the density in a qualitative manner. With increasing T /T c , c/c 0 increases for a small density n 2D , but decreases for large density. Note that c eventually vanishes in the thermal regime for high n 2D . This result indicates that at all densities, the running-wave measurement primarily excites the Bogoli- ubov mode at temperatures below the transition temperature. However, above the transition temperature, the potential quench primarily excites the normal sound mode at low density, or weak-coupling, and the diffusive mode at high densities. As we describe below, the same trend is visible in the standing-wave experiment. Furthermore, the dynamic structure factor that we discuss in Sec. VII supports this scenario as well. We emphasize that, in general, both modes are excited in these exper-iments. However, the amplitudes of the excited states are in general very different so that only one mode is detectable. In Sec. VIII we present a proposal for exciting both modes simultaneously with detectable amplitudes. In Figs. 3(b) and 3(c) we show the damping rate Γ and the quality factor Q, respectively. Γ shows a densitydependent behavior as a function of T /T c , which translates into a density dependence of the temperature dependence of Q. As a comparison we depict the prediction for the Q factor, which assumes that Landau damping is the primary mechanism for the line broadening, see Refs. [33,39]. The comparison shows good agreement.

V. SCALE INVARIANCE
Here we examine the scale invariance of c across the transition. We first demonstrate the scale invariance of the phase-space density D = n 2D λ 2 . We calculate D at various T /T c for the same three densities as before. In Fig. 4(a) we show the results of the inverse phasespace density D −1 determined at various T /T c . The different n 2D results collapse on a single line all across the transition. We compare them with the scaling prediction where D c is the critical phase-space density [38]. The simulations are in excellent agreement with the prediction. This confirms the universal scale invariance of the phase-space density.
We now test the scale invariance of the dimensionless sound velocities c 0 /c T and c/c T , where we refer to c T = k B T /m as the thermal velocity. The scaling prediction for c 0 /c T is (c 0 /c T ) scale = D cg T c /(2πT ), which depends only on T /T c , whileg is a fixed parameter. We show the results of c 0 /c T and c/c T in Figs. 4(b) and 4(c), respectively. The results of c 0 /c T collapse on a single line and agree very well with (c 0 /c T ) scale . This is a direct consequence of the data collapse shown in Fig. 4(a). However, for c, obtained from the simulation, the different n 2D results do not collapse on a single line, which shows a breaking of scale invariance regarding the sound velocity. The results of c/c T and the prediction (c 0 /c T ) scale agree only at low T , whereas they deviate at intermediate and high T . At low T , the damping of the sound mode is small compared to the mode frequency, i.e. Γ ω. However, at high T , Γ is comparable to ω and the deviation from the scaling prediction increases. The magnitude of the damping is expressed as a velocity, and shown as errorbars in Fig. 4(c). The deviation from the scaling prediction is comparable to the errorbars, suggesting that this breaking of scale invariance is due to the damping of the sound mode. Near and above T c , c undergoes the density-dependent changes that we have pointed out in the previous section.  Fig. 3(b).

VI. STANDING WAVES
As a second measurement, we analyze standing waves for the same system parameters as in Sec. IV. As an illustration, we choose n 2D ≈ 27 µm −2 and T /T c = 0.23, and create standing waves by periodically modulating the excitation potential, following the scheme described in Sec. III. After 1 s excitation time, we calculate the squared amplitude q 2 (y, f ) of the density modulation at varying modulation frequency f , see Sec. III for details. We show the results of q 2 (y, f ) determined as a function of f in Fig. 5(a). This response demonstrates the excitation of the first three standing waves at their mode frequencies. We fit the spatial dependence of q 2 (y, f ) with the function q 2 (y, ω j ) = j B 2 j cos 2 (k j,ω y/2) to determine the amplitudes B 2 j , where j is the mode index. We show the extracted amplitudes B 2 j of the standing waves in Fig.  5(b). We fit B 2 j with a Lorentzian function to determine the mode frequency f j and the damping rate Γ j . We show the determined f j and Γ j in the insets of Fig. 5(b). f j increases linearly with j, which demonstrates that the simulated standing waves correspond to the first three lowest-energy spatial modes. Γ j also increases linearly with j, which is a feature that is consistent with Landau damping.
We use the lowest-energy standing wave to determine the sound velocity c = ω 1 L y /π and compare it to the running wave measurement in Sec. IV. For the example given in Fig. 5, we obtain c = 1.48 mm/s, which agrees very well with c = 1.49 mm/s of the running wave measurement. We extend the comparison between the two measurements to the low and high n 2D systems across the transition. We use V 0 in the range V 0 /µ = 0.1 − 0. for all simulations. We present the results of standing and running wave simulations in Fig. 6.
In addition, we compare the simulation results to the Bogoliubov estimate of the sound velocity at nonzero temperature. We expressψ in the density-phase representation asψ(r) = n + δn(r) exp iφ(r) , where δn andφ are the density and phase fluctuations, respectively. From Eq. 1 we obtain the linearized Hamiltonian where n s is the superfluid density. The long-wavelength excitations are sound waves with velocity Following our description in Sec. II, this is the second (first) sound estimate for weak (strong) interactions. We calculate c B (T ) by numerically determining n s (T ) using the current-current correlations, see Appendix A. In Fig.  6 we present the results of c B (T ) determined for lowand high-n 2D systems. c B (T ) shows a density-dependent behavior and is nonzero above the transition. As mentioned above, the sudden jump of the superfluid density is replaced by a crossover regime due to the finite size of the system. Both the running-wave and the standingwave measurement are consistent with the Bogoliubov estimate below the transition. For low densities, both measurements show an upward trend as the temperatures approach the crossover regime. As it was demonstrated for the running-wave measurement earlier, this upward trend continues at temperatures above the critical temperature. This again suggests the interpretation that for the low density regime the normal sound mode is excited at higher temperatures. For high densities, the measured velocities both show a downward trend above the transition temperature. The standing wave measurement follows the Bogoliubov estimate closely, while the running wave measurement stays at a slightly higher value before it approaches zero as well. These measurements indicate that the primary excitation is the diffusive mode, while being slightly sensitive to the specific excitation method.

VII. DYNAMIC STRUCTURE FACTOR
We calculate the dynamic structure factor where n(k, ω) is the Fourier transform of the density n(r, t) in space and time. We determine n(k, ω) via N l is the number of lattice sites and T s = 328 ms is the sampling time for the numerical Fourier transform. The dynamic structure factor displays the overlap of the density degree of freedom with the collective excitations. We calculate S(k, ω) at various T /T c for low and high n 2D . In Fig. 7(a) we show S(k, ω) as a function of the wavevector k = k y and frequency ω for low n 2D across the transition. At low T , S(k, ω) has most of its weight at the Bogoliubov branch. At intermediate T , an additional branch with higher velocity appears. For comparison, we plot the Bogoliubov spectrum  8. S(k, ω) plots at k = 0.8 µm −1 and 1.7 µm −1 , for low n2D (upper row) and high n2D (lower row). The color scheme is the same as Fig. 7. The vertical dashed lines mark the frequencies of the Bogoliubov dispersion shown in Fig. 7.
, where c B (T ) is determined numerically, as above. k = 2J 1 − cos(kl) is the freeparticle spectrum on the lattice that is introduced to perform the numerical work, and J =h 2 /(2ml 2 ) is the tunneling energy. This dispersion recovers the continuum dispersion for l → 0. The Bogoliubov spectrum agrees well with the lower excitation branch at all k, for all T below T c . With this, we identify the lower branch as the Bogoliubov (B) mode and the upper branch as the non-Bogoliubov (NB) mode. This additional peak is also visible in Fig. 8(a), where the dynamic structure factor is depicted at two fixed values of the momentum. As illustrated in Sec. II, the faster mode is the NB mode and the slower mode is the B mode, for this density regime. Near T c , the B mode vanishes and becomes the diffusive mode, while the NB mode continuously connects to the normal sound mode of a thermal gas. Furthermore, the broadening of the B mode is visible, which corresponds to Landau damping, discussed before.
In Fig. 7(b) we show S(k, ω) for a high density n 2D . At low T , the weight is again mainly on the Bogoliubov branch, similar to the case of low density. At intermediate T , an additional branch with a lower velocity appears, in contrast to the case of low density where the velocity was higher. This corresponds to the second scenario described in Sec. II. These two branches are also visible in Fig. 8(b). We note that the dispersion of the Bogoliubov mode is renormalized to slightly higher values due to level repulsion between the two branches. Furthermore, both branches are broadened more strongly than for low densities, due to the higher interaction. This results in overlapping branches. At the transition, the B mode crosses over into the normal sound mode while the second sound mode transforms into the diffusive mode. The diffusive mode is broader than for low densities, and has higher weight. This leads to the previous observation that for this regime it is the diffusive mode that is primarily excited with a perturbation of the density.

VIII. EXCITATION OF BOTH SOUND MODES
We propose to excite both sound modes simultaneously by using a step-pulse density perturbation which is created by suddenly turning on and off the Gaussian potential at the location L y /2, see [40]. We choose the excitation time to be about 1 ms. For low-n 2D and various T /T c , we excite sound modes using both attractive and repulsive potentials. For all simulations, we use σ = 2 µm and V 0 in the range V 0 /µ = 0.25−2. We show the results in Fig. 9. At T /T c = 0.5, the time evolution of the density profile ∆n(y, t) shows primarily the excitation of the B mode. We do not observe significant NB mode excitation at and below T /T c = 0.5 compared to the numerical noise. At higher T /T c , the time evolution shows both B and NB mode excitations which are characterized by two density pulses traveling at different velocities. The NB mode travels faster than the B mode. At T /T c = 1, the B mode transforms into the diffusive mode and the NB mode into the normal sound mode of a thermal gas. Above T c , the time evolution shows primarily the normal sound propagation, as well as diffusive dynamics at the location of the perturbation. We fit the density profile with one or two Gaussians to determine the locations of one or two density pulses. From these locations we determine the sound velocities. To cancel out nonlinear effects due to the perturbation potential, we estimate the average squared velocity c 2 = (c 2 att + c 2 rep )/2, where c att (c rep ) corresponds to the attractive (repulsive) potential.
In Fig. 10 we show the temperature dependence of the two mode velocities of the step-pulse excitation. For comparison, we determine the mode velocities from the dynamic structure factor of low-n 2D shown in Fig. 7(a). We fit the excitation spectrum in the low-energy regime with a Lorentzian function to determine the mode frequencies. The NB mode frequency is determined after subtracting the background of the B mode. From the frequencies of the NB and B mode, we determine the first and second sound velocity, respectively. We show these results for various T /T c in Fig. 10. The results of the dynamic structure factor show excellent agreement to those of the step-pulse excitation. Overall, the first sound velocity shows a weak temperature dependence across the transition and is in the range c/c 0 = 2.6 − 2.8. The second sound velocity decreases with increasing temperature and vanishes above T c . The second sound results are in good agreement with the Bogoliubov estimates and the running-and standing-wave velocities. This set of results correspond to the first scenario of Sec. II.

IX. CONCLUSIONS
We have studied the propagation of sound in a 2D quasi-condensate of 87 Rb atoms across the superfluidthermal transition using c-field dynamics. We have identified two sound modes. To determine one or both velocities of these modes, we employ several methods. The first two methods are inspired by Ref. [33]: we excite running and standing waves with a weak Gaussian potential, from which we obtain a single velocity. Our simulations are in good agreement with the measurements of Ref. [33]. Furthermore, we have determined the dynamic structure factor. It displays two sound modes, and provides information about the overlap of these modes with the density degree of freedom. Below the critical temperature, one of the modes is the Bogoliubov mode. We refer to the other mode as the non-Bogoliubov mode. Above the critical temperature, we find the normal sound mode and the diffusive mode. The modes that are detected in Ref. [33] are the Bogoliubov mode below the critical temperature and the normal and the diffusive mode above the critical temperature.
As a key observation, we find that the non-Bogoliubov mode can have a higher or a lower velocity than the Bogoliubov mode. For weak interactions or low densities, the non-Bogoliubov mode has a higher velocity than the Bogoliubov mode, while for stronger interactions or higher density, the Bogoliubov mode has the higher velocity. While the strongly interacting regime is consistent with a hydrodynamic two-fluid approach, the weakly interacting regime provides a non-hydrodynamic scenario for the collective modes of Bose-Einstein condensates. We propose to measure the two sound modes simultaneously via a step-pulse density perturbation. By choosing the weak and the strong coupling regime of a condensate, these two regimes can be identified, which provides insight into this dynamical regime of condensates.
In the limit k → 0, these correlations are approximated by ( [41,42]) (j * k ) l (j k ) m = k B T m A n s k l k m k 2 + n n δ lm . (A2) n s and n n are the superfluid and the normal fluid density, respectively. A is the system area. We analyze the correlations along the line k x = k y = k/ √ 2 and determine the k = 0 value using a linear fit in the low-k regime. This allows us to determine n s at temperature T following Eq. A2. In Fig. 11 we show the determined n s as a function of T /T c for low and high n 2D . n s /n 2D shows a densitydependent behavior and no jump at the transition, due to the finite size of the system. For comparison, we numerically determine the condensate density n 0 and show this result in Fig. 11. n s and n 0 show good agreement for low density, while they deviate for high density. We note that for finite systems the condensate density scales algebraically with the system size, where the scaling ex-ponent is associated with the superfluid density [43].