Quench dynamics of an ultracold two-dimensional Bose gas

TL;DR

This paper investigates the non-equilibrium dynamics of a 2D ultracold Bose gas after a rapid quench across the BKT transition, revealing phase ordering and vortex decay consistent with theoretical scaling laws.

Contribution

It provides the first detailed simulation-based analysis of vortex decay and phase ordering kinetics in a quenched 2D Bose gas, aligning with experimental feasibility.

Findings

Vortex number decays as t^{-1}

Dynamical critical exponent z ≈ 2

Phase ordering kinetics confirmed

Abstract

We study the dynamics of a two-dimensional Bose gas after an instantaneous quench of an initially ultracold thermal atomic gas across the Berezinskii-Kosterlitz-Thouless phase transition, confirming via stochastic simulations that the system undergoes phase ordering kinetics and fulfills dynamical scaling hypothesis at late-time dynamics. Specifically, we find in that regime the vortex number decaying in time as $\langle N_v \propto t^{-1}\rangle$, consistent with a dynamical critical exponent $z \approx 2$ for both temperature and interaction quenches. Focusing on finite-size box-like geometries, we demonstrate that such an observation is within current experimental reach.

Equations

$$i\hbar\frac{\partial\psi(x,y,t)}{\partial t}=\hat{\mathcal{P}}\left\{(1-i\gamma)\left[-\frac{\hbar^{2}\nabla^{2}}{2m_{p}}+\right.\right.\\ \left.g_{\textrm{2D}}|\psi(x,y,t)|^{2}-\mu\bigg{]}{\psi(x,y,t)}+\eta(x,y,t)\right\}\,,$$

$${ϵ}_{\mathrm{cut}}(\mu ,T)=k_{\mathrm{B}}T\log (2)+\mu ,$$

$$g^{(1)}(r)=\frac{⟨{\psi }^{\ast }({\vec{r}}_0)\psi ({\vec{r}}_0+\vec{r}){⟩}_{{\vec{r}}_0,\theta ,\mathcal{N}}}{\sqrt{⟨|\psi ({\vec{r}}_0){|}^2{⟩}_{{\vec{r}}_0,\theta ,\mathcal{N}}⟨|\psi ({\vec{r}}_0+\vec{r}){|}^2{⟩}_{{\vec{r}}_0,\theta ,\mathcal{N}}}},$$

$$g^{(1)}(r){|}_c={\left(\frac{r}{{\lambda }_T}\right)}^{-{\alpha }_{\mathrm{c}}},$$

$$\frac{\mu}{k_{\mathrm{B}}T}\bigg{\rvert}_{\mathrm{BKT}}\approx\frac{\tilde{g}}{\pi}\ln\left(\frac{C}{\tilde{g}}\right),$$

$$T_{\mathrm{BKT}}^{\infty }=\frac{{\hslash }^2\pi n}{m_p{k}_{\mathrm{B}}\ln \left(C/\tilde{g}\right)},$$

$$g^{(1)}(r)\equiv \frac{1}{2}[g_x^{(1)}(r)+g_y^{(1)}(r)]$$

$$g_x^{(1)}(r)$$

$$g_y^{(1)}(r)$$

$$g_{\text{exp}}^{(1)}(r)\propto e^{-r/\xi }\text{and}{g}_{\text{alg}}^{(1)}(r)\propto r^{-\alpha }.$$

$${\chi }_{\mathrm{fit}}^2=\sum _i\frac{(g^{(1)}(r_i)-g_{\mathrm{fit}}^{(1)}(r_i){)}^2}{(g^{(1)}(r_i){)}^2},$$

$$n_{\mathrm{qc}}=\frac{\sqrt{2⟨|\psi (\vec{r}){|}^2{⟩}_{\vec{r},\mathcal{N}}^2-⟨|\psi (\vec{r}){|}^4{⟩}_{\vec{r},\mathcal{N}}}}{n},$$

$$n_{\text{sf}}=\frac{1}{{\lambda }_T^2\alpha },$$

$$U=\frac{\langle\big{|}\sum_{i,j}\psi_{i,j}\big{|}^{4}\rangle_{\mathcal{N}}}{\langle\big{|}\sum_{i,j}\psi_{i,j}|^{2}\rangle_{\mathcal{N}}^{2}},$$

$$m=\frac{1}{\sqrt{N_x{N}_y}}\frac{⟨\left|{\sum }_{i,j}{\psi }_{i,j}\right|{⟩}_{\mathcal{N}}}{⟨{\sum }_{i,j}{\left|{\psi }_{i,j}\right|}^2{⟩}_{\mathcal{N}}}.$$

$${⟨N_{\text{v}}⟩}^{-1/2}\propto L,$$

$$g^{(1)}(r,t)\sim g_{\mathrm{eq}}^{(1)}(r)F\left(\frac{r}{L(t)}\right),$$

$$⟨N_{\text{v}}(t)⟩\propto t^{\beta },$$

$$L_{\mathrm{fit}}\propto t^{-\beta /2}\mathrm{and}{⟨N_{\mathrm{v}}⟩}_{\mathrm{fit}}\propto t^{\beta },$$

$${\tilde{g}}_{\mathrm{in}}\approx 0.21\approx 1.19{\tilde{g}}_{\mathrm{crit}}\to {\tilde{g}}_{\mathrm{fin}}\approx 0.12\approx 0.68{\tilde{g}}_{\mathrm{crit}},$$

$${\tilde{g}}_{\mathrm{in}}=0.1\approx 5.46{\tilde{g}}_{\mathrm{crit}}\to {\tilde{g}}_{\mathrm{fin}}=1.83\times 10^{-3}\approx 0.1{\tilde{g}}_{\mathrm{crit}}.$$

Full text

I Introduction

A two-dimensional (2D) Bose gas is known to undergo a Berezinskii-Kosterlitz-Thouless (BKT) phase transition Berezinskii (1972); Kosterlitz and Thouless (1973) between a superfluid and a non-superfluid regime, originating from the binding and unbinding of vortex-antivortex pairs. Such a transition, first observed in thin helium films Bishop and Reppy (1978), has been experimentally investigated in a broad range of systems including ultracold atomic gases in quasi-2D harmonic potentials Stock et al. (2005); Cladé et al. (2009); Hadzibabic et al. (2006); Krüger et al. (2007); Tung et al. (2010); Yefsah et al. (2011) and, more recently, in near-uniform box-like traps Chomaz et al. (2015). Theoretically such settings have been studied with a universal $|\Psi|^{4}$ model on a lattice Prokof’ev, Ruebenacker and Svistunov (2001); Prokof’ev and Svistunov (2002), a semiclassical field approach Giorgetti et al. (2007), Quantum Monte Carlo methods Holzmann and Krauth (2008); Holzmann et al. (2010), classical fields and the stochastic (projected) Gross-Pitaevskii equation applied to both harmonically trapped gases Simula and Blakie (2006); Simula et al. (2008); Bisset et al. (2009); Cockburn and Proukakis (2012); Mathey et al. (2017) and box-like geometries Foster et al. (2010); Karl and Gasenzer (2017); Gawryluk and Brewczyk (2019), and a renormalization group approach Pelissetto and Vicari (2013). In the weakly-interacting regime, such predictions have shown good agreement with experimental findings for the interference fringes in expansion dynamics Hadzibabic et al. (2006); Polkovnikov et al. (2006), the relation between the number of vortices and the emergence of phase coherence Chomaz et al. (2015), the scale invariance and universality Hung et al. (2011); Cockburn and Proukakis (2012), and the propagation of sound near the BKT transition Ville et al. (2018); Ota et al. (2018).

In cold gases, for a given atomic species, the critical temperature of the BKT phase transition is set by the chemical potential and effective interaction strength. These quantities can be experimentally controlled by fixing the atom number and varying either the strength of the transverse confinement or the scattering length by means of a Feshbach resonance, where available not (a). The BKT transition can be crossed by lowering the temperature via evaporating cooling the gas. The dynamics emerging in a temperature quench are sensitive to the quench protocol, and in particular to the quench rate. A quench across a continuous phase transition spontaneously generates defects in the order parameter. A relation between the defect density and the quench rate is provided by the Kibble-Zurek universal scaling law Kibble (1976); Zurek (1985). In the case of an instantaneous (or sufficiently rapid) quench, the emerging dynamics reveal a phase-ordering stage in which self-similar correlation functions collapse onto each other when scaled in terms of a characteristic lengthscale Bray and Rutenberg (1994); Rutenberg and Bray (1995). The growth of this lengthscale, which is directly connected to the defect dynamics, is set by the dynamical critical exponent $z$ which can thus be extracted from the simulations in the appropriate late-time evolution stage.

The purpose of this work is to theoretically study the phase-ordering kinetics of a quenched ultracold atomic Bose gas in a box trap, within existing and envisaged 2D box geometries, mostly inspired by recent experiments realized in the Laboratoire Kastler Brossel (LKB) in Paris with 87Rb Chomaz et al. (2015); Ville et al. (2018); Saint-Jalm et al. (2019) and currently underway at Cambridge (AMOP) with 39K Had . We model the gas by means of the stochastic projected Gross-Pitaevskii equation Stoof and Bijlsma (2001); Stoof (1999); Gardiner and Davis (2003); Bradley and Gardiner and Davis (2008); Blakie et al. (2008); Proukakis and Jackson (2008). Considering the geometry of the LKB box trap, we first perform a detailed analysis of the equilibrium configuration as a function of temperature, with our findings revealing good qualitative agreement with earlier numerical works. Having identified the relevant regimes, we then discuss controlled instantaneous quenches across the BKT phase transition. Firstly, we verify the expected bulk predictions by considering both temperature and interaction quenches in the limit of a large ‘idealized’ box – larger than currently accessible experimentally. Having demonstrated our methodology and confirmed its predictive power both in terms of correlation functions and vortex dynamics, we then specifically address the feasibility of experimental observation. Correlation functions – which are hard to measure directly in 2D box geometries – are likely to be prone to finite size effects in the currently-accessible box sizes, suggesting it may be favorable for observations to rely instead on vortex numbers. Considering the specific finite extent of the recent LKB experimental set-up in Paris Ville et al. (2018), we show that even in such finite-size systems, there is a well-defined temporal window in which one should be able to observe the long-term $t^{-1}$ evolution in the vortex number, associated with vortex-antivortex annihilation processes. We show that the same conclusion is valid also in the different setting inspired by ongoing experiments of AMOP in Cambridge. Similar findings have been previously obtained in the long-term evolution of the closely-related problem of decaying two-dimensional quantum turbulence, connecting ultracold atom experiments Kwon et al. (2014); Kim et al. (2016); Seo et al. (2017); Neely et al. (2013) with numerical studies Baggaley and Barenghi (2018), a problem discussed more generally in the context of dynamical vortex decay via $N$-vortex collision processes Baggaley and Barenghi (2018); Groszek et al. (2016, 2019); Karl and Gasenzer (2017).

II Stochastic Gross-Pitaevskii Equation

We consider a weakly interacting ultracold Bose gas confined in a transverse tight potential sufficiently strong that it creates a uniform 2D system which occupies the lowest energy state in the transverse direction. The stochastic (projected) Gross-Pitaevskii equation (SPGPE) describes its dynamics via the noisy complex field $\psi(x,y,t)$, subject to the equation Stoof and Bijlsma (2001); Stoof (1999); Gardiner and Davis (2003); Bradley and Gardiner and Davis (2008); Blakie et al. (2008); Proukakis and Jackson (2008)

[TABLE]

where $\nabla^{2}$ is the Laplacian in two dimensions and $g_{\textrm{2D}}$ is the 2D coupling constant associated with the s-wave scattering length $a_{s}$, and $m_{p}$ is the mass of the particles. If the transverse confinement is harmonic, $V(z)=m_{p}\omega_{z}^{2}/2$, one can define the harmonic length $\ell_{\perp}=\sqrt{\hbar/m_{p}\omega_{z}}$ and the dimensionless coupling constant $\tilde{g}=m_{p}g_{\textrm{2D}}/\hbar^{2}=\sqrt{8\pi}{a_{s}}/{\ell_{\perp}}$. The projector $\hat{\mathcal{P}}$ constrains the dynamics of the system within a finite number of macroscopically occupied modes, the coherent region, up to an ultraviolet energy cutoff fixed here as

[TABLE]

for which the mean occupation of the last included mode is of order $\sim 1$, assuming a Bose-Einstein distribution in the occupation number spectrum. A Gaussian stochastic noise $\eta(x,y,t)$ with correlation in space and time $\langle\eta^{*}(x,y,t)\eta(x^{\prime},y^{\prime},t^{\prime})\rangle=2\hbar\gamma k_{\mathrm{B}}T\delta(x-x^{\prime})\delta(y-y^{\prime})\delta(t-t^{\prime})$ is added and projected in the coherent region. The cutoff is implemented in Fourier space. The numerical grid for the simulations is determined from $\epsilon_{\mathrm{cut}}$ by the anti-aliasing condition $\Delta x\leq\pi/\sqrt{8m\epsilon_{\mathrm{cut}}}$ Blakie et al. (2008). In comparing numerical results with experiments, one must include the density of atoms in the coherent region $n_{c\textrm{-field}}=\int dxdy\,|\psi(x,y)|^{2}$ as well as in the incoherent region $n_{\mathrm{I}}$, so that the total density is $n=n_{c\textrm{-field}}+n_{\mathrm{I}}$. Here $n_{\mathrm{I}}$ is determined by the density of states $G(\epsilon)$ of the system as $n_{\mathrm{I}}=1/(L_{x}L_{y})\int_{\epsilon_{\mathrm{cut}}}^{\infty}d\epsilon\,G(\epsilon)/(e^{(\epsilon-\mu)/k_{\mathrm{B}}T}-1)$, where the gas in the incoherent region is assumed to be ideal. Finally, the dissipative term $\gamma$, parametrizing the interaction between the high-lying and the low-lying modes of the system, has the practical role of setting the rate at which the system reaches the equilibrium determined by the external parameters temperature $T$ and chemical potential $\mu$. In our simulations we use values of $\gamma$ within a decade centered in the value of $\gamma=0.01$, of the same order of the values used in Liu et al. (2018) where $\gamma$ was fixed to reproduce the time growth of the number of atoms occupying the lowest-momentum mode in a 3D condensate subject to a temperature quench across the BEC transition Donadello et al. (2016).

III Equilibrium: the Berezinskii-Kosterlitz-Thouless phase transition

III.1 Background Theory

The Mermin-Wagner-Hohenberg theorem Mermin and Wagner (1966); Hohenberg (1967) states that, for a system of dimensions 2 or fewer and short range interactions, it is not possible to have the spontaneous breaking of a continuous symmetry at any non-zero temperature. As a notable consequence, there is no Bose-Einstein condensation in these geometries, since the thermal fluctuations at any temperature are strong enough to destroy the long-range coherence in the system. Repulsive interactions, however, may result in the establishment of quasi-long range coherence at sufficiently low temperatures, affecting the behavior of the first-order correlation function defined in an isotropic system as

[TABLE]

where the average $\langle\dots\rangle_{\vec{r}_{0},\theta,\mathcal{N}}$ is performed over the spatial position $\vec{r}_{0}$, the angular part $\theta$ of $\vec{r}$ and over a large number $\mathcal{N}$ of stochastic realizations. In a 2D system, as shown in Kosterlitz and Thouless (1973), the topological properties of the system are determined by the behavior of quantized vortex pairs with respect to a critical temperature $T_{\mathrm{BKT}}$ for the infinite order Berezinskii-Kosterlitz-Thouless phase transition. Specifically:

• •

for $T>T_{\mathrm{BKT}}$ free vortices can exist in the system, there is no superfluid and the first-order correlation function decays as $g^{(1)}(r)\sim e^{-r/\xi}$, where $\xi$ is a correlation length;

• •

for $T<T_{\mathrm{BKT}}$ vortices can only exist in bound pairs, allowing the presence of a superfluid, and the correlation shows an algebraic decay $g^{(1)}(r)\sim r^{-\alpha}$ in terms of an exponent $\alpha$.

Note that the algebraic decay of the correlation function in the degenerate case could lead to very strong requirements in terms of the sample size, for behavior consistent with the Mermin-Wagner-Hohenberg theorem to manifest itself. At the transition ($T=T_{\mathrm{BKT}}$), the correlation function should decay according to Kosterlitz and Thouless (1973); Kosterlitz (1974)

[TABLE]

where $\lambda_{T}=\sqrt{2\pi\hbar^{2}/{m_{p}}k_{\mathrm{B}}T}$ is the thermal de Broglie wavelength, and $\alpha_{\mathrm{c}}=0.25$ Nelson and Kosterlitz (1977). At the thermodynamic limit, when the volume and total number of particles of the system tend to infinity while the density is fixed, the value for the critical temperature is determined by

[TABLE]

where the constant $C$ has been estimated by Monte-Carlo analysis in Prokof’ev and Svistunov (2002) to be $C\sim 13.2$. Thus, by inverting Eq. (5), the temperature at the critical point is estimated as

[TABLE]

where $n$ denotes the system density.

III.2 Numerical Results

We now characterize the equilibrium state of a finite-size uniform 2D Bose gas, focusing on the experimental geometry at LKB. We consider approximately 20,000 87Rb atoms with mass ${m_{p}}=1.4431\times 10^{-25}$\mathrm{kg}$$ in a uniform two-dimensional box of size $L_{x}\times L_{y}=(50\times 50)\mu$m, whose transverse confinement is $\omega_{z}=2\pi(1500)\textrm{Hz}$. The chemical potential is fixed at $\mu/k_{\textrm{B}}=4.8\textrm{nK}$, and $\tilde{g}=9.5\times 10^{-2}$. For our choice of parameters, the system is numerically equilibrated for a time $t_{\mathrm{eq}}=3\times 10^{5}\omega_{z}^{-1}\approx 30$s at a temperature spanning an interval $T/T_{\mathrm{BKT}}^{\infty}=[0.13,1.96]$, where $T_{\mathrm{BKT}}^{\infty}=33.25\mathrm{nK}$. Our simulations were performed by means of the XMDS2 software package described in Dennis et al. (2013), on the High-Performance-Computing cluster at Newcastle University.

To obtain a smooth first-order correlation function, $g^{(1)}(r)$, we average over both $x$ and $y$ directions in our discretized grid (of $N_{x}\times N_{y}$ points, respectively) using the expression

[TABLE]

where

[TABLE]

where $\psi_{i,j}=\psi(x_{i},y_{j})$ and the average $\langle\dots\rangle_{\mathcal{N}}$ is performed over the number $\mathcal{N}$ of stochastic realizations. Once this function is computed, we fit it (as in Dagvadorj et al. (2015); Comaron et al. (2018)) with the functions

[TABLE]

According to the BKT theory, the former should apply above $T_{\mathrm{BKT}}$ and the latter below. As a measure of the quality of the fit we use the quantity

[TABLE]

where the index $i$ accounts again for the spatial discretization.

The results are given in Fig. 1. The equilibrium correlation functions and the corresponding fits are shown in the left panel for the selected temperature range, while the $\chi^{2}$ functions are shown on the right. As expected, the algebraic fit is better than the exponential fit (i.e., has a lower $\chi^{2}$) at low temperatures, while the exponential fit is better at high temperatures. The crossover between the two behaviors occurs in a narrow region close to $T_{\mathrm{BKT}}^{\infty}$, the two values $\chi^{2}_{\rm exp}$ and $\chi^{2}_{\rm alg}$ being equal slightly above $T_{\mathrm{BKT}}^{\infty}$. The shaded area between $0.95<T/T_{\mathrm{BKT}}^{\infty}<1.15$ represents qualitatively the region where the BKT transition occurs.

In the top panel of Fig. 2 we show the parameter $\alpha$ extracted from the algebraic fit. The shaded area is the same as in Fig. 1. At the left border of this area we find $\alpha\approx 0.25$, which corresponds to the prediction of BKT theory at the transition in the thermodynamic limit Kosterlitz (1974); Nelson and Kosterlitz (1977) and in agreement with Nazarenko et al. (2014). At the right border we instead find $\alpha\approx 0.5$. This value agrees with the one obtained in Foster et al. (2010) just above $T_{\mathrm{BKT}}$ in numerical simulations within a classical field approach in the micro-canonical ensemble, reproducing the interference patterns experimentally observed in Hadzibabic et al. (2006). A broad critical region, as identified here, is also consistent with the existence of an “intermediate regime” characterized by $\alpha>0.25$ and a rapidly vanishing zero momentum current-current correlation, as pointed out in Gawryluk and Brewczyk (2019), again within a classical field calculation but in the grand-canonical ensemble. Such a region is predicted to shrink when approaching the thermodynamic limit Foster et al. (2010); Gawryluk and Brewczyk (2019).

Although no true condensation can occur in 2D in the presence of interactions, the quasi-long range coherence leads to the formation of a so called “quasi-condensate”, which can be thought of as a condensate with a fluctuating phase Kagan et al. (1987). Its density can be computed via Prokof’ev and Svistunov (2002)

[TABLE]

where the average of the moduli $\langle\dots\rangle_{\vec{r},\mathcal{N}}$ is performed over the entire spatial grid and all noise realizations. As shown by the green line in Fig. 2(b) – and consistent with earlier works Foster et al. (2010) – the quasi-condensate density has a significant non-zero value even above the critical temperature, with $n_{\mathrm{qc}}(T_{\mathrm{BKT}}^{\infty})\sim 0.7$. In the framework of the BKT theory, the quasi-condensate facilitates the existence of vortices above the critical temperature, where a superfluid is absent. In the same panel of Fig. 2 we also report the normalized superfluid density $n_{\mathrm{sf}}$ (red line), defined as

[TABLE]

as suggested by Nelson and Kosterlitz Nelson and Kosterlitz (1977), which is meaningful only at low $T$ where the algebraic fit to $g^{(1)}$ is reliable. We also plot the condensate fraction $n_{0}$ (blue dashed line), defined as the normalized density of particles which populate the zero momentum ($k=0$) mode. At low $T$, the lowest mode tends to saturate at the same value as $n_{\rm qc}$, while above $T_{\mathrm{BKT}}$ it is just a small fraction, as expected for a strongly fluctuating quasi-condensate. We note in passing that in our simulations the $c$-field fraction lies between $0.84<n_{c\textrm{-field}}/n<0.99$ depending on the temperature. Here, the $c$-field density is calculated as $n_{c\textrm{-field}}=\langle\sum_{i,j}\left|\psi_{i,j}\right|^{2}\rangle_{\mathcal{N}}$. At the critical point we have $n_{c\textrm{-field}}/n(T_{\mathrm{BKT}}^{\infty})=0.92$.

A related quantity characterizing the location of the phase transition is the Binder ratio (or Binder cumulant) Binder (1981); Kobayashi and Cugliandolo (2016), defined as

[TABLE]

and plotted in Fig. 2(c). This quantity is predicted to be a step-function from $1$ (fully coherent system) to $2$ (pure thermal state) in the limit of infinitely large boxes. In finite volumes, it is instead particularly sensitive to finite-size effects, resulting in a smooth function with a slope progressively small with smaller boxes.

Another way to characterize the physical regime of the system is by counting the average number of vortices at equilibrium as a function of temperature; this is done by means of a numerical routine that calculates the phase winding around each grid point, to identify the vortices and their circulation. The results are given in Fig. 2(d). Note that, for the parameters of our simulations, all vortices annihilate for temperatures lower than $\sim 0.75T_{\mathrm{BKT}}^{\infty}$.

Finally, a measure of the degree of degeneracy of the system was introduced in Ref. Kobayashi and Cugliandolo (2016) in the form of an order parameter $m$ which, in our discretized space, is defined as

[TABLE]

This quantity accounts for the off-diagonal terms of the system’s density matrix. Its value at equilibrium, $m_{\mathrm{eq}}$, is plotted in Fig. 2(e) together with the value of $\sqrt{g_{1}(r)}$ calculated at the edges of the system. These quantities should coincide in the limit of large boxes Kobayashi and Cugliandolo (2016). The behaviour of $m_{\mathrm{eq}}$ is qualitatively similar to the one observed in Kobayashi and Cugliandolo (2016) for a planar 2D $XY$ model, showing a rapid decrease above the critical temperature.

All quantities plotted in Fig. 2 demonstrate the occurrence of the BKT transition within a critical range of $T$. Due to finite-size effects the transition is not sharp. Depending on the way one conventionally defines a precise value of the transition temperature $T_{\mathrm{BKT}}$, such a value can be slightly shifted downward, as in Gawryluk and Brewczyk (2019), or upward, as in Foster et al. (2010), with respect to the ideal value $T_{\mathrm{BKT}}^{\infty}$, but the different definitions are expected to converge in the thermodynamic limit. Our results demonstrate clear consistency with such earlier findings, as well as with the results of similar analysis in exciton-polariton 2D gases Comaron et al. (2018); not (a).

IV Quench and phase ordering dynamics

An important part of our analysis is devoted to the study of the relaxation dynamics of coherence and topological defects of the 2D Bose gas across the critical temperature. This study can reveal universal properties through the characterization of critical exponents, which address the universality class of the system considered Bray and Rutenberg (1994). Such a process has been largely studied in conservative Damle et al. (1996); Karl and Gasenzer (2017) and open systems Kulczykowski and Matuszewski (2017); Comaron et al. (2018). The 2D Bose gas is known to belong to the same universality class of the planar 2D $XY$ model, whose dynamics has been studied theoretically Jelić and Cugliandolo (2011) and experimentally Nagaya et al. (1992, 1995).

In Eq. (1), the parameters controlling the physical state of a given atomic species are $T$, $\mu$ and $g_{\mathrm{2D}}$. In theoretical studies assuming a fixed value for the interaction strength $g_{\mathrm{2D}}$, quenching through the critical point of the transition has been implemented by manipulating the temperature parameter $T$ while keeping the chemical potential unaltered Weiler et al. (2008); Jelić and Cugliandolo (2011), by quenching $\mu$ maintaining a fixed temperature Damski and Zurek (2010), and by performing a simultaneous quench of the two quantities Liu et al. (2016, 2018).

IV.1 Background Theory

We start by briefly reviewing the theory of the phase ordering process. The quantities which characterize the nonequilibrium dynamics of a 2D Bose gas are the correlation length $L(t)$, the average length within which coherence is established, and the average number of topological defects detected when crossing the phase transition. In our case, the defects are quantized vortices and their average number at a time $t$ relates to $L$ as

[TABLE]

and the corresponding density is $\left<n_{\mathrm{v}}\right>=\left<N_{\mathrm{v}}\right>/L_{\mathrm{box}}^{2}$, where from now on $\left<\dots\right>$ corresponds to an average over the stochastic realizations. After a sudden quench from a temperature above $T_{\mathrm{BKT}}$, at sufficiently long times, the system is expected to enter a temporal region where the system exhibits universal dynamical scaling Yurke et al. (1993); Bray and Rutenberg (1994), in the sense that the (non-equilibrium) correlation function (3) evolves in time according to the dynamical scaling form Bray and Rutenberg (1994)

[TABLE]

where $g^{(1)}_{\mathrm{eq}}({r})$ is the static (equilibrium) correlation function for the final parameters $\mu,T$ and $\tilde{g}$, decaying algebraically at long distances as $g^{(1)}_{\mathrm{eq}}(r)\propto r^{-\alpha}$. From this scaling law one can extract the correlation length $L(t)$, which is predicted to exhibit a power law dependence on the inverse of time, whose exponent is referred as the dynamical exponent $z$. Accordingly, the average number of vortices is expected to behave as

[TABLE]

where $\beta=-2/z$. The dynamical exponent for the diffusive dynamics of the 2D $XY$ model is predicted to be $z\approx 2$  Jelić and Cugliandolo (2011). Note that due to the presence of steady state vortices, logarithmic corrections to this law are also expected, such that $\left<N_{\textrm{v}}(t)\right>\propto\ \left({t}/{\ln(t/t_{0})}\right)^{\beta}$, where $t_{0}$ is a nonuniversal temperature-dependent timescale (we refer to Refs. Jelić and Cugliandolo (2011); Comaron et al. (2018) for details). As logarithmic corrections only manifest themselves in very large systems and are unlikely to be directly detected in ultracold atomic experiments, we do not focus on such corrections here, simply noting that the numerical results presented here are in principle consistent with the presence of weak logarithmic corrections.

IV.2 Results for a large box

In this section we discuss the rather idealized case of a large 2D box trap, to reveal the key physics expected in the homogeneous thermodynamic limit. Specifically we simulate here a $(200\times 200)\mathrm{\mu m}$ box with periodic boundary conditions, following an instantaneous quench from above to below the BKT phase transition, as a grand canonical evolution.

Firstly, we simulate an infinitely-rapid temperature quench across the critical point of a ${}^{87}\mathrm{Rb}$ gas by evolving Eq. (1) from an equilibrium initial configuration above the critical point to a quasi-ordered state. Similarly to earlier works Liu et al. (2016, 2018), we induce a simultaneous jump in the system chemical potential from $\mu<0$ to symmetrically-located (about zero) $\mu>0$ values. Specifically, we prepare our system in an equilibrated disordered state with temperature $T_{\mathrm{in}}=200\ \mathrm{nK}\gg T_{\mathrm{BKT}}$ and $\mu_{\mathrm{in}}=-2.4\ k_{\mathrm{B}}\ \mathrm{nK}<0$ and induce at $t=0$ a sudden quench by setting values for a chosen final state with $T_{\mathrm{fin}}=5\ \mathrm{nK}\ll T_{\mathrm{BKT}}$ and $\mu_{\mathrm{fin}}=2.4\ k_{\mathrm{B}}\ \mathrm{nK}>0$. The chosen values of $\mu$ and $T$ also set the cutoff for each stage of the system evolution. Interaction strength $\tilde{g}=95\times 10^{-3}$ is fixed by the transverse confinement adopted, namely $\omega_{z}=2\pi(1500)\mathrm{Hz}$.

Typical snapshots of the classical field phase distribution in a single temperature quench are shown in Fig. 3(a),(b) and (c); the three images are taken before (green circle), close to (red circle) and far after (orange circle) the critical point, showing the process of creation and annihilation of vortex-antivortex pairs during the BKT transition. From such distributions, calculated at different times and in many realizations, we extract the correlation function $g^{(1)}$ and the average number of vortices. The results are shown in Fig. 3(d),(e). The green, red and orange circles in Fig. 3(e) correspond to the three snapshots of the top row.

Inspired by Jelić and Cugliandolo (2011); Comaron et al. (2018), we investigate whether the system exhibits universal dynamics in terms of a correlation length $L$. To this aim, we first plot the ratio $g^{(1)}/g^{(1)}_{\mathrm{eq}}$ as a function of $r$ at different times as in the inset of Fig. 3(d). Following Refs. Kobayashi and Cugliandolo (2016); Kulczykowski and Matuszewski (2017), we then extract the length $L(t)$ by imposing $(g^{(1)}/g^{(1)}_{\mathrm{eq}})(L(t),t)=0.2$, and finally we plot $g^{(1)}/g^{(1)}_{\mathrm{eq}}$ again, but as a function of the rescaled distance $r/L(t)$. As a result, all curves nicely collapse onto a single one, except at large distances where boundary effects become relevant. This confirms the universal dynamical scaling. We have also checked that the dependence on the value of $g^{(1)}/g^{(1)}_{\mathrm{eq}}$ chosen for the determination of $L(t)$ is weak and can be neglected.

The curve of $L(t)$ is shown in Fig. 3(e), together with the average number of vortices, as a function of the dimensionless time $\gamma t$. A key result is that both quantities are found to be related by Eq. (15) within the appropriate region where universal scaling is satisfied (grey shaded area in Fig. 3(e)). Such a window is chosen at late enough times to allow the system to enter the universal regime not (b), but before size effects start to play a role in the growth of $L$ Jelić and Cugliandolo (2011); Comaron et al. (2018). By fitting correlation length and average number of vortices with

[TABLE]

respectively, we are able to extract the growth exponents for both quantities and we find $\beta\approx-1$. This is consistent with a value of the dynamical exponent $z=-2/\beta=2$, as for the case of the planar $XY$ model (with non-conserved order parameter) dynamics, which belongs to the same dynamical universality class. It is also in agreement with the results obtained within a microcanonical evolution of the system, solved with classical field methods, from a many-vortices configuration Karl and Gasenzer (2017). The largest uncertainty in our result for $\beta$ comes from the choice of the parameter $\gamma$ of the SPGPE. We have performed simulations for different values of $\gamma$ within a decade centred in the value of $\gamma=0.01$ (which is the range that we think is reasonable for describing our system) and observed $\beta$ to be compatible with $z=2$ within a $10\%$ accuracy. The variations of $\gamma$ also reflect on the time values reported in Figs. 3 and 4. We have also explicitly verified – by changing the size and shape of our vortex counting region within the box – that our numerical extraction of $z=2$ remains unaffected by loss of vortices on the box boundaries.

In order to make more direct contact with controllable experiments, we also implement a quench in the interaction (i.e., in the coupling constant $\tilde{g}$). Experimentally, varying the value of $\tilde{g}$ is possible by changing the transverse confinement $\omega_{z}$ (as in the LKB setup), or by tuning the value of the scattering length $a_{s}$ thanks to the Feshbach resonances Ha et al. (2013); Clark et al. (2017); Fu et al. (2018). By keeping the control parameters of the reservoir $T_{\mathrm{res}}$ and $\mu_{\mathrm{res}}$ fixed, we tune the parameter $\tilde{g}$ such that the system is suddenly quenched from $\tilde{g}_{\mathrm{in}}$ to $\tilde{g}_{\mathrm{fin}}$ with $\tilde{g}_{\mathrm{fin}}<\tilde{g}_{\mathrm{crit}}<\tilde{g}_{\mathrm{in}}$, with $\tilde{g}_{\mathrm{crit}}$ defined by Eq. (5). Note that fixing $T_{\mathrm{res}},\mu_{\mathrm{res}}$ and reducing $\tilde{g}$, causes the critical temperature to increase according to Eq. (6) (distinct to what would happen if one were to instead keep the density constant during the quench). The chemical potential $\mu_{\mathrm{res}}={16.5k_{\mathrm{B}}\mathrm{nK}}$ and the temperature $T_{\mathrm{res}}={70\mathrm{nK}}$ of the reservoir (incoherent region) are chosen in order to have a number of atoms in the box in the order of tens of thousands, consistent with the experimental set-up. Following Eq. (5), the critical value then reads $\tilde{g}_{\mathrm{crit}}=0.18$.

The interaction quench protocol we adopt consists of two stages: first, we prepare our initial equilibrium, disordered state $\psi_{\mathrm{in}}=\psi(\tilde{g}_{\mathrm{in}},T_{\mathrm{res}},\mu_{\mathrm{res}})$ by evolving Eq. (1) starting from random noise. Then we suddenly quench to the quasi-ordered state $\psi_{\mathrm{fin}}=\psi(\tilde{g}_{\mathrm{fin}},T_{\mathrm{res}},\mu_{\mathrm{res}})$.

For the initial and final parameter sets ($\tilde{g},T_{\rm res}$ and $\mu_{\rm res}$), we choose values which lie sufficiently above ($\tilde{g}_{\mathrm{in}}$) and below ($\tilde{g}_{\mathrm{fin}}$) the critical region, in order to avoid any potential weak dependencies associated with the precise location of the critical region on the implemented energy cutoff. Specifically, here we use the values

[TABLE]

which, for the parameters of the LKB experiment, would correspond to relaxing the transverse harmonic confinement $\omega_{z}$ from an initial value $\omega_{z}\leavevmode\nobreak\ \rightarrow\leavevmode\nobreak\ \omega_{z}/3$, where $\omega_{z}\leavevmode\nobreak\ =\leavevmode\nobreak\ 2\pi(7500)\ \mathrm{Hz}$.

Results are shown in Fig. 3(f),(g). As already done for temperature quenches, we can extract the correlation length $L(t)$ such that the curves of $g^{(1)}/g^{(1)}_{\mathrm{eq}}$ collapse onto a universal function of $r/L(t)$. Again we find a time window where the length $L$ and the average number of vortices scale as $t^{-\beta/2}$ and $t^{\beta}$, respectively, with $\beta\approx-1$ as before.

IV.3 Results for experimentally-relevant boxes

In this section we focus on an interaction quench. Contrary to quenching $T$ and $\mu$, the interaction quench is experimentally easier to implement as it requires only a single parameter to be used to drive the system across the phase transition. Interaction quenches are generally performed in the lab by rapidly changing either the scattering length $a_{s}$ on times scales $\approx\leavevmode\nobreak\ 100\ \mathrm{\mu s}$, or the vertical oscillator length $\ell_{\perp}$ (namely changing the transverse confinement potential $\omega_{z}$). This protocol is well-known and already adopted for harmonic trap potentials Fletcher et al. (2015) and uniform systems Ville et al. (2018).

Realistic systems also have finite-size geometries. For this reason in this section we use a box of size $(40\times 40)\mathrm{\mu m}$, consistent with system sizes of the experiment described in Ref. Ville et al. (2018), also displaying results for $(100\times 100)\mathrm{\mu m}$ to highlight the role of finite-size effects in the smaller systems. Given the difficulties in measuring $g^{(1)}$ experimentally, we choose to focus the following analysis on the evolution of vortex number, which is also better tractable in numerics, particularly for small system sizes.

First, we repeat the same sudden interaction quench simulation discussed in Fig. 3 (bottom), but for a smaller system where a more realistic hard bound potential box is implemented. As in the ideal case, we tune the reservoir temperature to $T_{\mathrm{res}}={70\mathrm{nK}}$, and use a modest experimentally-accessible change in transverse confinement, while keeping fixed the value of $\mu_{\mathrm{res}}={16.5k_{\mathrm{B}}\mathrm{nK}}$. Thus, the range of the interaction values is as in Eq. (19). Once the system density has reached saturation, we note that the coherent atom fraction reads $n_{\textrm{c-field}}/n\approx 0.81$, with $N\approx 1.2\times 10^{4}$. The numerical cutoff along the whole evolution is fixed by final values of temperature $T_{\mathrm{{res}}}$, and chemical potential $\mu_{\mathrm{res}}$.

In addition to the set of parameters of the LKB experiments, we also consider here a rather distinct set of parameters inspired by the AMOP experimental set-up Campbell et al. (2010); Fletcher et al. (2015) at Cambridge. In this experiment, the choice of 39K bosons additionally facilitates – at least in principle – the control of the interaction strength through a Feshbach resonance. Similarly to the previous case, the temperature and chemical potential of the reservoir set the critical value of the interaction strength $g_{\mathrm{{crit}}}$. We choose $T=50\,\mathrm{nK}$ and $\mu=1.9\,k_{\mathrm{B}}\mathrm{nK}$, corresponding to values of temperature and chemical potential at the end of the evaporation stage of an AMOP experiment Had . Correspondingly, the critical interaction strength becomes $\tilde{g}_{\mathrm{crit}}\left(\mu,T\right)=1.83\times 10^{-2}$, so that the quench protocol reads

[TABLE]

At the end of the simulation the total number of atoms reads $N\approx 29\times 10^{4}$, where $n_{c\textrm{-field}}/n\approx 0.98$.

In Fig. 4 we show the temporal evolution of the average number of vortices obtained with both sets of parameters. The dashed lines correspond to the universal scaling behavior $\left<N_{\mathrm{v}}\right>\propto t^{-1}$. In order to start observing such a scaling law one has to wait at least a time $t\approx 0.002\gamma^{-1}$ after the sudden quench for the LKB parameters and $t\approx 0.02\gamma^{-1}$ for the AMOP parameters, corresponding to $\approx 0.2$ s and $\approx 2$ s, respectively, based on the value of $\gamma=0.01$ used in our simulations. This difference in the timescale for phase ordering to set in can be explained in terms of the different growth rate of coherence that we predict in the two systems with such a value of $\gamma$. In the inset of Fig. 4 we show the time evolution of the classical field density (red line) and the average number of vortices (blue line) in linear scale. The classical field density is found to grow faster with the LKB parameters; a factor of approximately ten difference in the growth rate is found by fitting the classical field density with an S-shaped growth curve, as done in Ref. Liu et al. (2018). This explains the temporal shift of the window where universal scaling is observed in our simulations. This regime occurs roughly at the point when the classical field density appears to approach saturation to its final value, at which however the system has not yet reached its full coherence for the corresponding reservoir parameters. A precise determination of this window in the experiments, however, should require an independent experimental estimate of the coherence growth rate. Nevertheless, time scales in between fractions of a second to ten seconds, as found here, are well within the typical observation time in current experiments. Also, these results appear to be rather independent of the system size.

V Conclusions

We have performed a detailed analysis of the dynamics following instantaneous temperature and interaction quenches from an incoherent thermal state to a superfluid state below the Berezinskii-Kosterlitz-Thouless phase transition in 2D Bose gases of ultracold atoms. Considering large boxes with periodic boundary conditions we have demonstrated the self-similarity of correlation functions in the phase-ordering regime, and characterized the evolution of the correlation length and vortex number as a function of time. Both were found to be consistent with a dynamical critical exponent $z=2$, as expected for diffusive dynamics of a system in the 2D $XY$ universality class. Using smaller boxes, we have also shown that realistic geometries, experimentally-accessible interactions and small instantaneous quenches across criticality are likely to facilitate the observation of a regime where $\left<N_{\mathrm{v}}\right>\propto t^{-1}$ (i.e. $z\approx 2$), thus providing direct measurement of the dynamics and critical exponent $z$ defining the system universality class. Our approach is also relevant for investigating the 2D/3D crossover in vortex dynamics, related to recent experiments on vortex clustering and turbulence Gauthier et al. (2018); Johnstone et al. (2018); Seo et al. (2017).

Author Contributions

P.C. and F.L. contributed equally to this work, undertaking all numerical simulations and analysis and producing the first draft, in direct consultation with F.D. and N.P.P, coordinating the research and leading interpretations. All authors contributed to discussions, final data analysis and interpretations, and the final form of the manuscript.

Acknowledgements.

We would like to thank C. Barenghi, T. Billam, T. Bland, J. Dalibard, J. Schmitt, Z. Hadzibabic and A. Groszek for fruitful discussions. We acknowledge financial support from the Quantera ERA-NET cofund project NAQUAS through the Consiglio Nazionale delle Ricerche (F.D. and F.L.) and the Engineering and Physical Science Research Council, Grant No. EP/R043434/1 (P.C. and N.P.P.). We acknowledge the Engineering and Physical Science Research Council, Grant No. EP/L504828/1, for DTA support. This work is also supported by Provincia Autonoma di Trento.