Berezinskii-Kosterlitz-Thouless phase induced by dissipating quasisolitons

We theoretically study the sound propagation in a two-dimensional weakly interacting uniform Bose gas. Using the classical fields approximation we analyze in detail the properties of density waves generated both in a weak and strong perturbation regimes. While in the former case density excitations can be described in terms of hydrodynamic or collisionless sound, the strong disturbance of the system results in a qualitatively different response. We identify observed structures as quasisolitons and uncover their internal complexity for strong perturbation case. For this regime quasisolitons break into vortex pairs as time progresses, eventually reaching an equilibrium state. We find this state, characterized by only fluctuating in time averaged number of pairs of opposite charge vortices and by appearance of a quasi-long-range order, as the Berezinskii-Kosterlitz-Thouless (BKT) phase.


Introduction
Sound waves carry information on both thermodynamic and transport properties of a medium they propagate through. In classical hydrodynamics, measuring the speed of sound waves and their attenuation gives an access to characteristics of the medium such as the compressibility and viscosity. In quantum hydrodynamics, with superfluids present, the picture is more complex 1 . For example, liquid helium and weakly interacting Bose gas respond to local perturbation in a qualitatively different ways.
Many experiments exploring the phenomenon of sound propagation in ultracold atomic systems have been already performed. Sound waves were studied in harmonically trapped three-dimensional bosonic gases at very low [2][3][4] as well as at higher 5 temperatures. Sound velocity was measured at resonance in a mixture of fermionic lithium atoms 6 . Two sound modes propagating at different speeds, according to predictions of two-fluid hydrodynamics, were observed in a resonant Fermi gas 7 . The sound diffusion in a unitary three-dimensional Fermi gas was investigated and a universal quantum limit of diffusivity was observed in Ref. 8 . Recently, sound propagation and damping were studied in two-dimensional systems -a weakly interacting Bose 9, 10 and strongly interacting Fermi 11 gases.
In 9 , a gas of rubidium-87 atoms is confined inside a quasi-2D rectangular potential with hard walls. A density perturbation is introduced by applying a repulsive potential to the cloud of atoms. This additional potential creates a density dip along one direction which propagates at constant speed when the laser is switched off. During this evolution the density perturbation bounces several times off the walls of the box and its velocity is found to be close to the Bogoliubov sound speed. Several attempts have been already undertaken to reproduce results of this experiment [12][13][14] . They all support experimental observation that below the BKT transition the generated sound waves propagate with velocities close to the speed of second sound, predicted by two-fluid hydrodynamic model. They also predict that sound waves can propagate above the BKT transition, in agreement with experiment and in opposition to what is predicted by two-fluid hydrodynamic model with respect to second sound 15,16 .

Sound waves propagation
In our numerical simulations we first start with analyzing the response of two-dimensional Bose gas to weak perturbation, as in the experiment of Ref. 9 . With the help of Metropolis algorithm 14,[17][18][19][20][21][22] we create an ensemble of initial states, ψ(r), for a two-dimensional Bose gas confined in a box potential with periodic boundary conditions. We work within the grand canonical ensemble, so the temperature T and the chemical potential µ are control parameters. Members of an ensemble of classical fields are drawn with the probability ∼ exp((µN − E)/k B T ), where N = ∑ k |α k | 2 and E = (h 2 /2m) ∑ k k 2 |α k | 2 + (g/2L 2 ) ∑ k,j,l α * k α * j α l α k+j−l are the number of particles and the energy of the state, respectively. Here, g/L 2 is the two-body interaction energy with L being the length of a square box potential and α k are the amplitudes of expansion of the classical field in the basis of the box eigenfunctions.
Next, we disturb the cloud of atoms with the protocol similar to that applied in 9 . We create a dip in density in one direction, just by multiplying the initial classical field by an appropriate steplike function. The dip is characterized by two parameters: the width w and the depth d. Such disturbance, a kind of density imprinting, is then evolved according to the classical field approximation (CFA) prescription. An equation for a disturbed classical field ψ(r,t) is 23 ih ∂ /∂t ψ(r,t) = (−h 2 /2m ∇ 2 + g|ψ(r,t)| 2 ) ψ(r,t) . (1) In our case the perturbation is initially located at the center of the box. As a consequence we have two density patterns propagating symmetrically outward as in the experiment of 2 . Integrating density along the direction transverse to the sound propagation we are able to identify traveling waves and find their speeds. To get the velocities of density waves we decompose the density integrated along the direction perpendicular to the wave propagation, which is y axis in our case, via the Fourier transform at each instant of time: n(x,t) = n +∑ j=±1,±2,... A j (t) exp( j 2πx/L). Here, n is the average density along x axis. Then we do time analysis of A j (t)/A j (0) coefficients -we fit real (or imaginary) part of them to an exponentially damped sinusoidal function e −Γ j t/2 [Γ j /(2ω j ) sin(ω j t) + cos(ω j t)]. Focusing on lower energy modes we obtain the velocity of density waves as v = ω j /k j , where k j is the wave vector (k j = ± j2π/L, since periodic boundary conditions are applied).
In Fig. 1 we summarize our results, showing velocities of density waves in units of the Bogoliubov speed c B = gN/L 2 m and assuming the interaction strength is g = 0.16h 2 /m. The temperature is given in units of T c = 2πnh 2 /[mk B ln (380h 2 /mg)], which is the calculated critical temperature for the BKT phase transition 24 . Two sets of equilibrium states, corresponding to different values of chemical potential (as indicated in the legend), are used for the analysis. Numerical values of the speed of sound waves (black and brown bullets, black circles) remain in a good agreement with experimental data (red squares 9 ) in the whole range of studied temperatures. Most of the data in Fig. 1 were obtained assuming a density imprinting parameters: w/L = 0.25 and d = 5/16. Some of the points (brown bullets) were calculated in the regime of very weak perturbation, d = 1/16. Close to the transition and for the normal phase we double results showing also the response of the system to stronger initial disturbance (black circles, for which d > 5/16). Numerical data clearly demonstrate that the response of the system depends on the strength of initial density imprinting. Similar observation was reported in Ref. 4 , where the wave fronts propagation through the condensate was studied after it was split by a strong laser beam.

2/11
Additionally, we put in Fig. 1 values of the speed of second sound, coming from the two-fluid model 1,15,16 . Within this framework one can calculate a value of the speed of second sound, c 2 , from the equation where n s , n n , and n (= n s + n n ) are superfluid, normal, and total densities, respectively. κ T (κ S ) is the isothermal (adiabatic) compressibility,s is the entropy, andc V -the specific heat at constant volume, both per unit mass. For temperatures such that k B T > g n, a good approximation to the velocity of second sound is given by the expression c 2 = n s /n/ √ mnκ T 15 . Within CFA we can calculate both the superfluid density (from current-current correlations 22 ) and the isothermal compressibility. The latter is obtained from the fluctuation-dissipation relation ( The speed of second sound, as a function of T /T c , is plotted in Fig. 1 by black crosses and the values are close to experimental (and CFA) results.
In CFA approach the classical field is expanded in the basis of modes which are macroscopically occupied -in our case these are simply the plane waves 23 . Dynamical equations for modes amplitudes (expansion coefficients) can be written in the form of a set of nonlinear differential equations. At each instant of time, these amplitudes determine fully the classical field. Hence, we can watch the density (square modulus of the classical field) with arbitrary high spatial resolution. In particular, this resolution can be high enough to resolve the healing length scale.

Figure 2.
Time evolution of initially perturbed equilibrium state. Two-dimensional density of the system is first integrated along one direction, next averaged over approximately 100 realizations, and then shown as a function of time. Perturbation parameters, the width w and depth d are given in the frames. Temperature is T = 0.37T c . Total evolution time t = 0.025 mL 2 /h corresponds to 50 ms for L = 38µm. Here, the spatial resolution is 1.64 of the healing length.

Quasisolitons
New observations are possible in a strong perturbation regime with high resolution employed. Interesting things happen when the density imprinting parameters, the width w and the depth d of the initial depletion region, are being increased. First of all, we find that the visible density waves, in fact, consist of several thinner structures, which are characterized by slightly different velocities, see Fig. 2. For weaker initial perturbations they rather form a single beam (Fig. 2, upper row), whereas for stronger ones we clearly see signs of nonlinear behavior -lines in time-position plots (Fig. 2, lower row) are no longer straight. In both cases, internal structures, traveling with slightly different velocities, separate each other. For stronger initial perturbation the multiple density dips accelerate during the motion. At the same time, larger the width of initial perturbation w bigger the number of created waves. The number of internal thin waves visible inside broad dips in Fig. 2 is well described by the formula w/(4ξ ), where w is the width of initial perturbation and ξ is the healing length. Factor 4 arises from a typical width of these structures, see Fig. 2, which turns out to be about 4ξ . Structures presented in Fig. 2 survive several crossings and exhibit a moderate change of shapes during the evolution, so one can try to assign some solitonic properties to them. Below we exam numerically obtained densities by comparing their profiles with analytical formulas for a single soliton solution of a nonlinear Schrödinger equation, known as Zakharov-Shabat solutions with the nonlinear term corresponding to repulsive interactions 26 Here, n 0 is the density of the medium far away from the soliton dip, q andq are the position and velocity of the soliton, respectively. The healing length ξ =h/ √ mn 0 g and c B = n 0 g/m. The width and the depth of solitonic solution (3) is solely determined by soliton velocity,q. For weak perturbations (top left frame in Fig. 2) density dips propagate with constant velocity but they get shallower in time (due to snake instabilities [27][28][29][30][31] ). Hence, additional damping factor has to be included in Eq. (3) and then a dissipative Zakharov-Shabat solution takes the form On the other hand, for strong perturbations (lower row in Fig. 2) the density waves accelerate in time. In this case, the parameteṙ q in Eq. (3) should depend on time.
To verify this concept we analyze first the simplest case -a small initial perturbation in both w and d which leads to single wave moving with constant velocity. In Fig. 3, frames (a) and (b), we compare numerically obtained density cuts for various times with Eq. (4), whereq is fixed by the value of v we get from the time-dependent analysis of A 1 (t) coefficient (actually, the numerical results are fitted to the sum of two one-soliton solutions of Eq. (4)). It turns out that both densities match perfectly up to first collision, after which the agreement slowly degrades. Surprisingly, Eq. (3) works well also for strong perturbation, although in this case one has to consider rather an accelerating two-soliton solutions of (3), i.e. a solution with time-dependentq. Frames (c) and (d) prove that an agreement is good, here w/L = 2/50 and d = 14/16. Hence, we will be naming observed structures as quasisolitons as opposed to true solitonic objects present in 1D but also in 2D, where they are called the Jones-Roberts (JR) solitons 32,33 . Similarly to our case (as discussed in the next section), the JR solitons can transform into vortex-antivortes pairs. However, this happens in the limit of vanishing velocity which critically distinguish JR structures from ours. Surprisingly, Eq. (3) also works for weak perturbations at longer times, i.e. after first collision of density waves. Density dips become, as time progresses, broader and shallower and their velocity increases towards speed of Bogoliubov sound (similarly to dissipative dynamics of solitons already observed in early experiments on dark solitons created by phase imprinting technique 34 and theoretically discussed in 35 ).
The excitation protocol we use in numerical simulations is rather a kind of ideal one. It differs from that applied in Ref. 9 in two ways. First, we generate the system at equilibrium in a box potential, i.e. without an additional potential repelling atoms out of desired space in the box, and only after that part of the atoms is removed from the sample. Second, the change of the density at the border of perturbed-unperturbed regions is extremely sharp. The second difference, as we checked, does not influence the results in weak perturbation regime. However, when the system is strongly disturbed, multi-quasisoliton structures visible in Fig. 2 get less pronounced. The same happens when the system is disturbed and evolved after it is equilibrated in the presence of repelling barrier.

Emergence of the BKT phase
As observed in simulations, quasisolitons dissipate and eventually breake into elementary vortex pairs (i.e. pairs of vortices with winding numbers equal to ±1). It is illustrated in Fig. 4, second and fifth rows, where densities for a single realization are shown at different times. Except most left frames almost no trace of density waves is visible, although, while averaged over a hundred realizations, density dips are still detectable for short times (first left three columns in first and fourth rows). Additionally, positions of pairs of opposite charge vortices are marked by color (white and red) circles for a single realization cases. We identify vortices by calculating the phase winding around each lattice plaquette as the system evolves. The number of vortex pairs changes in time but finally, for times longer than 0.1 mL 2 /h, an equilibrium is established -see plateau in Fig. 5 (middle frame), exhibiting only fluctuating, after averaging over a hundred realizations, number of vortex pairs. This feature sustains over hundreds of milliseconds as shown in inset in Fig. 5. The final number of vortex pairs depends on the strength of initial density perturbation, the stronger disturbance the larger number of pairs. For longer times no trace of quasisolitons' positions is visible (see most right column in Fig. 4) and motion of pairs of vortices resembles rather that of an ideal gas particles.
Production of vortex-antivortex pairs in two-dimensional condensates, due to snake instability of initially imprinted dark soliton, was thoroughly discussed in 36 . As demonstrated theoretically in 37 , creation of vortex-antivortex pairs is expected to occur in a laser-stirred two-dimensional Bose gas, as in the experiment of 38 . There, an additional energy is pumped into the system due to stirring the gas at high enough velocity. This energy dissipates via mechanism of creation of vortex-antivortex pairs. Vortex dynamics itself features fast (related to annihilation of vortices of opposite charges) and slow (related to drifting vortices out to the thermal region of the cloud) decay in number of vortices. Remarkably, the first channel is quickly closed and the system exhibits very slow relaxation, on time scales of seconds. The somehow opposite route, which is the crossover from the BKT phase 39,40 containing pairs of vortices to a vortex free Bose-Einstein condensate (BEC) in weakly interacting quasi-two-dimensional Bose gas was experimentally realized in 41 .
In numerical simulations, initially prepared sample of the Bose gas at equilibrium is modified by immediate removing some part of the atomic cloud. Then, in our case the total energy of the system is decreased and preserved during further evolution as the system is isolated. How the energy per atom is changed by such a perturbation is shown in Fig. 6. What is most important is the behavior of kinetic energy per atom, since it determines the temperature of the system. Fig. 6 tells us that this quantity does not change much. Since the atomic density gets lower, the transition temperature is reduced (see 24 ). Then effectively, a two-dimensional Bose gas is shifted towards the transition point because T /T c gets larger. This supposition is confirmed also by analyzing the properties of the system in plateau region in Fig. 5, middle frame. It turns out that there exist such values of the chemical potential and temperature that the grand canonical ensemble formalism predicts the average energy and particles number just as those exhibited by the system in plateau region in Fig. 5.
As visible in Fig. 5, the number of vortex pairs increases with the strength of perturbation. Left frame in Fig. 5 shows effective ratio T /T c after disturbance. In the background a superfluid fraction is plotted as a function of relative temperature T /T c (black dots) for the Bose gas of the chemical potential µ = 2500h 2 /mL 2 (see 22 ). The most right line in Fig. 5 (left frame) indicates the transition temperature as argued in 22 , based on consideration of current-current correlations. The superfluid density at the transition temperature represents the superfluid density jump characteristic for the BKT theory as shown by Nelson and Kosterlitz 42 . What we observe in our simulations is then an opposite route to that realized in 41 , it is rather the crossover from a vortex free BEC to the BKT phase containing pairs of vortices.
Presence of vortices in plateau region, Fig. 5 (middle frame), could be a signature of appearance of the BKT phase 43,44 .  Fig. 5). The spatial resolution is 1.12µm which is 1.39 of the healing length.

6/11
Here, we actually observe an appearance and persistence of pairs of vortices which is a strong signature for the BKT phase 39,40 .
As we checked, the average number of vortex pairs does not change over a time of hundreds of milliseconds (see inset in Fig.  5), suggesting that the system remains frozen in the BKT state. To clarify the situation we calculate the first-order correlation function g (1) (r, r ,t), defined as a normalized average g (1) (r, r ,t) = ψ (r,t) ψ(r ,t) / ψ (r,t) ψ(r ,t) over the grand canonical ensemble, while the system evolves. According to Kosterlitz and Thouless 39,40 it should decay algebraically with a distance in a uniformn two-dimensional Bose gas. We plot the results in Fig. 7 for a strong perturbation with w/L = 12/34 and d = 14/16. We find that already when the number of vortex pairs approaches plateau (see Fig. 5), the first-order correlation function decays algebraically (Fig. 7, left and middle frames in top row). However, the correlations are not yet spatially uniform. They behave as ∼ 1/r α with a distance but with different exponents in x and y directions. However, for t 0.10 mL 2 /h the system starts exhibiting the quasi-long-range order, i.e. g (1) (r, 0) drops algebraically with a distance, here with exponent equal to ≈ 0.13 for t > 0.20 mL 2 /h. Hence, the system enters the BKT phase. To check the consistency of our simulations we calculate the superfluid density of the system which, according to the BKT theory 42     Details of the relaxation dynamics of our system, i.e. those exhibited in a period before the plateau region is reached, are depicted in Fig. 5. We find the number of vortices is decreasing in time. Several experiments, studying two-dimensional turbulent flow in bosonic superfluids, have already reported a decay of vortices after an energy was pumped to the system through transient stirring [45][46][47][48] . In 48 two distinct regimes have been identified -the one showing t −2/5 and the other exhibiting t −1 decay in time of number of vortices. Such vortex annihilation behavior can be understood in terms of collisions involving few vortices. In Ref. 49 it was predicted that vortices annihilate mainly through two-and three-body loss events. These processes result in specific time decay of the number of vortices, governed by ∝ t −1 and ∝ t −2/5 dependence, respectively. Experimental results of 48 show reasonable agreement with this prediction for some kinds of disturbance protocols. Fig. 5 (right frame) clearly shows the existence of two different regimes while the system approaches an equilibrium. We uncover faster (∝ t −1 ) and slower (∝ t −2/5 ) decrease of number of vortices as well. In our case, a slow decay is succeeding the fast one. This must be related to the perturbation protocol we use. It leads to creation of pairs of vortices (vortex dipoles) due to snake instability of quasisolitons rather than to spatially random distributed vortices. Evidently, in the first stage two-body head-to-tail collisions of vortex dipoles are dominant, whereas later on all other relative orientations of vortex dipoles get possible which, effectively, corresponds to the three-body collision events.
Yet another transition can be realized with the perturbation protocol we apply. Now, we start with the two-dimensional Bose gas at higher temperature, T /T c = 0.8. The system remains in the BKT phase 22 , see Fig. 7, most left frame in bottom row. The first-order correlation function decays algebraically. Then we strongly disturb the Bose gas with w/L = 10/36 (middle frame) and w/L = 20/36 (right frame). After stronger perturbation the correlation function starts to decrease exponentially with a distance, the system looses the quasi-long-range order. Hence, we observe the BKT to normal phase transition in this case.
Moreover, we checked that also in the case of less ideal excitation protocols the long-term evolution exhibits plateau in the number of vortex pairs as in Fig. 5. Such behavior is expected because the transition temperature T c is again reduced since the average density is lower than the unperturbed one, due to the presence of additional atoms-repelling barrier.

Summary
In summary, we have studied the propagation of density waves in a two-dimensional weakly interacting uniform Bose gas. In a weak perturbation regime, the velocities of sound waves agree with the experimentally measured. When the system is strongly disturbed, its response becomes complex. We then identify patterns of multiple density dips, propagating with slightly different velocities. Size of these structures is of submicron and their shape coincides with dissipative Zakharov-Shabat profile. These

8/11
structures dissipate, changing into pairs of opposite charge vortices. On a scale of hundreds of milliseconds an equilibrium is reached, characterized by only fluctuating in time averaged number of pairs of opposite signs vortices and by appearance of a quasi-long-range order -the system enters the BKT phase. The BKT phase appears as a result of an increase of system's relative temperature T /T c .