Periodic quenches across the Berezinskii-Kosterlitz-Thouless phase transition

The quenched dynamics of an ultracold homogeneous atomic two-dimensional Bose gas subjected to periodic quenches across the Berezinskii-Kosterlitz-Thouless (BKT) phase transition are discussed. Specifically, we address the effect of periodic cycling of the effective atomic interaction strength between a thermal disordered state above, and a highly ordered state below the critical BKT interaction strength, by means of numerical simulations of the stochastic projected Gross-Pitaevskii equation. Probing the emerging dynamics as a function of the frequency of sinusoidal driving from low to high frequencies reveals diverse dynamical features, including phase-lagged quasi adiabatic reversible condensate formation, resonant excitation consistent with an intrinsic system relaxation timescale, and gradual establishment of dynamically-recurring or time-averaged non-equilibrium states with enhanced coherence which are neither condensed, nor thermal. Our study paves the way for experimental observation of such driven non-equilibrium ultracold superfluid states.

The related question of how the presence of a sinusoidally-modulated driving across a phase transition may affect the system dynamics arises naturally. Such a periodic cycling was experimentally investigated in the context of three-dimensional (3D) harmonically trapped ultracold atoms, where a time-dependent dimple microtrap was used to controllably induce a periodic phase space modulation, leading to a reversible cycling across the thermal and the Bose-Einstein condensation phase [43]. Condensate formation was found to lag behind the applied constant-frequency modulation of the laser power. Such findings were subsequently reproduced qualitatively by means of numerical simulations based on the stochastic Gross-Pitaevskii equation [44]. Motivated by the above pioneering works [43,44], and by our recent studies of instantaneously quenched two-dimensional quantum gases [26,45], here we address the corresponding periodically driven phase transition crossing in the context of (quasi) two-dimensional homogeneous ultracold atomic Bose gases. Two-dimensional systems are interesting in their own right, due to the different nature of the underlying Berezinskii-Kosterlitz-Thouless (BKT) phase transition [46,47], associated with binding-unbinding of vortex-antivortex pairs, previously observed in diverse physical contexts [48][49][50][51]. Such properties have also been studied in ultracold atomic systems both experimentally [38,[52][53][54][55][56][57][58][59][60], and theoretically [26,[61][62][63][64][65][66][67][68][69][70][71][72][73][74][75] with the transversal (typically harmonic) confinement in the other direction offering a way to control the effective two-dimensional interaction experienced between the atoms.
In this work, we perform a detailed quantitative analysis of the role of driving frequency on the cyclic phase transition crossing, focussing on the particular case of a periodically-driven 2D homogeneous ultracold atomic Bose gas. This is facilitated by a periodically-modulatedin-time interaction strength, between an initial incoherent state close to, but above, the BKT phase transition, and a state with lower interaction strength below the BKT phase transition exhibiting a high degree of coherence. Our numerical study, performed at fixed system temperature, focusses in parallel on the effect of external driving on density, vortex number, momentum spectrum and coherence. Our parameter choice is based on accessible experimental regimes, building on our earlier work based on an instantaneous interaction quench from above to below the BKT threshold, which focussed on equilibrium properties and late-time phase-ordering dynamics [26].
Deep in the superfluid regime, where the system is highly condensed, modulations of the interaction strength have been studied both experimentally and theoretically, in different contexts. Periodic modulations of the interaction strength between two values in the superfluid regime -conducted within the context of the arXiv:2006.16683v1 [cond-mat.quant-gas] 30 Jun 2020 Gross-Pitaevskii equation, for which the system is assumed to be coherent, and often termed 'Feshbach Resonance Management' [76] -have demonstrated the emergence of Faraday patterns [77] and have been used to study aspects of condensate stability and soliton dynamics [78,79]. Condensate experiments based on a sudden interaction modulation found an interesting analogy with Sakharov oscillations in the early Universe [80], while periodic interaction strength modulations led to matterwave jet emission (Bose fireworks) [81,82] and other interesting patterns [83].
The novel feature of our present study is that our periodic interaction quenches are neither conducted between two states deep in the superfluid regime, nor do they start from a well-formed condensate -but instead the interaction strength is modulated across the phase transition multiple times, in direct analogy to the work of [43]. For this reason, our analysis is based on the stochastic (projected) Gross-Pitaevskii equation [84][85][86][87].
Beyond the expected regimes of extremely slow pumping (which allows the system to proceed adiabatically through instantaneous equilibrium states) and extremely rapid pumping which only mildly perturbs the initial state, we find the periodic driving to be intrinsically resonant with a characteristic relaxation time, at which the periodic modulation of the scattering length causes strongly-non-equilibrium features to emerge. This resonance separates two interesting distinct and experimentally relevant physical regimes: Driving frequencies lower than the resonant value lead to reproducible dynamics which are independent of the quench cycle, and -while not fully equilibrated -resemble certain features of the corresponding-parameter equilibrium states. This regime bears close analogy to the experimental 3D findings of Ref. [43] (upon excluding observed atom losses). In the opposite regime of frequencies exceeding the resonant value, the system grows gradually over multiple quench cycles, and can -for frequencies few times the resonant frequency -accommodate significant coherence. This paper is structured as follows: After presenting a brief overview of the 3D experiment of Ref. [43], we outline the numerical scheme and quench protocol to be used in our work (Sec. II). We start by discussing the dependence of the emerging density and vortex dynamics on driving frequency, and present a detailed quantitative analysis of their respective dynamical time-delay with respect to the external periodic driving (Sec. III): this allows us to identify an intrinsic system frequency, and thus to focus our subsequent analysis on the physically interesting and experimentally relevant regime of driving frequencies which are up to 20 times slower, or faster, than the identified resonant frequency. Within this experimentally interesting range where most novel dynamical features emerge, we analyse the dependence of the momentum spectrum on the driving frequency and identify parameter regimes of high coherence in spite of fast driving (Sec. IV). We then study the maximum attained coherence as a function of quench frequency, and the emerging non-equilibrium steady-state under rapid pumping (Sec. V), before concluding.

II. QUENCH PROTOCOL AND MODELLING SCHEME
A. Experiment of Stamper-Kurn et al. [43] In a ground-breaking paper, Stamper-Kurn et al. [43] achieved reversible condensate formation in a threedimensional gas of sodium atoms. Initially, this gas was confined to a broad harmonic trap and cooled to just above the critical temperature, T c , at which condensation onsets. Subsequently, a thin well (dimple trap) was introduced at the centre of the trap, using an infrared laser. This granted the sodium atoms access to a new, lower-energy state. The well was made so steep that it could only contain a single new energy level. As well depth was increased, by increasing the strength of the laser, the energy of this single state decreased to the extent that a condensate formed in this well. The condensate fraction grew with well depth, and hence with laser power. The power of the laser beam was then sinusoidally modulated between 0 and 7mW at a frequency of 1Hz, which resulted in the periodic cyclic growth and decay of the condensate fraction, N 0 /N , from 0 to 6%. A non-zero condensate fraction was observed to recur for 15 oscillations, even though the peak condensate fraction decayed slightly with each oscillation. This decay was shown to arise from atom loss, rather than as a consequence of this periodic crossing of the phase transition. The latter was confirmed by the fact that the same decay in the peak values was observed even when the power of the laser beam was held constant. Importantly, a ∼ 70ms time delay was observed between the times of maximum laser power and corresponding time of peak condensate fraction, giving some insight into the condensate formation time under this protocol.
To interpret the observed findings, Stoof and Bijlsma [44] pioneered the use of a stochastic Gross-Pitaevskii equation as a model to study reversible condensate formation. Introducing a steep well in the potential, as done experimentally, is equivalent to altering the effective chemical potential, µ eff , of the system. Hence, in this numerical investigation, the reversible formation of a one-dimensional condensate was simulated by modulating the value of µ eff as µ eff (t) = µsin(ω D t), where µ is the chemical potential of the system in the absence of a well and ω D is the frequency of the imposed periodic modulation in µ eff (t). Stoof and Bijlsma analysed the evolution of the condensate density at the trap centre and qualitatively reproduced the delay in condensate growth observed by Stamper-Kurn et. al [43] (but observing a slightly longer time of ∼100ms, which could be attributed to the different probed dimensionality and potential uncertainties of exact experimental parameters). That study provided the first numerical evidence of a repeated phase transition crossing in ultracold atomic gases. Since that work, the stochastic Gross-Pitaevskii model introduced by Stoof [44], and its closely-related stochastic projected Gross-Pitaevskii equation [84,85], have been used to study a broad range of dynamical ultracold phenomena [85][86][87][88], including the study of quenched phase transitions across different geometries, dimensionalities and mixtures [18, 20, 22-24, 26, 30, 32-34, 89-98], with direct successful modelling of phase transition experiments [18,24,34].

B. Numerical Model and Parameter Regime
In this work we analyse the quenched system dynamics in terms of the Stochastic Projected Gross-Pitaevskii Equation (SPGPE), which describes the 'classical' field, Φ(r, t), of all highly-populated modes up to a fixed energy cut-off. Its dynamics are governed by [85]: Here −iγ corresponds to a dissipative/growth term, arising from the coupling of the 'classical field' modes Φ to the high-lying modes, which are treated as a heat bath. Consistent with other treatments [26,98], we set γ = 0.01 in our current calculations, noting that typical physical evolution timescales are set by the scaled time γt, controlling system dynamics and growth. The long-term evolution/steady-state is determined by the balancing of the kinetic energy contribution and the nonlinear interaction term g|Φ(r, t)| 2 with the bath chemical potential µ. The presence of the dynamical noise term η, associated with the collisional randomness of growth/decay processes differentiates each numerical run in a manner analogous to experimental shot-to-shot fluctuations. The noise correlations are given by η * (r, t)η(r , t ) = 2 γk B T δ(t − t )δ(r − r ). Finally,P is a projector which constrains the dynamics of the system within a finite number of macroscopically occupied modes, up to an ultraviolet energy cutoff cut (µ, T ) = k B T log(2) + µ, with the mean occupation of the last included mode set to be ∼ 1. This relation sets our numerical grid spacing to ∆x ≤ π/ √ 8m cut [85]. In this work, we consider a two-dimensional, weakly interacting, homogeneous Bose gas, confined to a square box of sides L x = L y = 100µm. Considering experimentally accessible regimes, our presented analysis is based on a gas of ultracold 39 K bosons with mass m = 6.47 × 10 −26 kg, a chemical potential µ = 1.92k B nK and a temperature, T = 50nK [99]. The corresponding 2D interaction strength is given in terms of the tight transversal harmonic confinement of frequency ω ⊥ by g 2D = g 3D / √ 2π ⊥ , where g 3D = 4π 2 a s /m, a s = 7.36nm is the background s-wave scattering length [100], and ⊥ = /mω ⊥ the transverse harmonic length. This yields a dimensionless coupling constant g = mg 2D / 2 = √ 8πa s / ⊥ . The location of the BKT phase transition is fixed by the relation [61] where g c denotes the critical value of the dimensionless strength at the BKT transition, and Monte-Carlo analysis gives the constant C ∼ 13.2 [62]. This relation sets the value of g c = 1.83 × 10 −2 for our parameters. This coupling constant can be varied by changing the transverse confinement, and the parameters considered in this work relate to the range 1.83×10 −3 < g < 3.48×10 −2 . Implicit in our description is the modulation of the scattering length through Feshbach resonance [101], or the transverse confining harmonic oscillator frequency ω ⊥ . All presented results correspond to an average over N = 50 stochastic realizations, which is large enough such that the error bars on all figures are of the same size order as the marker itself.

C. Quench Protocol
Stimulated by the work of Refs. [43,44], we impose here a symmetric interaction quench, as follows: Firstly, we allow the system to dynamically equilibrate close to, and on the disordered side of, the BKT critical region. Specifically we choose to equilibrate initially to a value g = 1.9g c (which corresponds to an initial number of N = 3 × 10 4 atoms in the classical field). We verify that this state is in equilibrium via calculation of the first order correlation function, which exhibits the correct exponential decay law, and ensuring the vortex number reaches a steady state value. We then initiate a periodic quench in the interaction strength which is symmetric about the critical point, via where, ω D is the driving frequency of the quench and α is the constant periodic amplitude. Based on our chosen initial state, we fix the modulation amplitude to α = 0.9, implying that the quench proceeds between an initial state with g = 1.9g c (corresponding to T > T BKT ), and a final highly-ordered state with g = 0.1g c (T < T BKT ) [102]. Throughout this work, we examine how the system behaviour varies for different driving frequencies in the range ω D = 2π × [1/40, 500]Hz, which cover the entire behaviour from quasi-adiabatic (ω D = 2π × (1/40)Hz) to extremely rapid (ω D = 2π × 500Hz), which fully encompass the entire range of potential parameters for experimental investigation.
In presenting our results, we introduce the time-  Evolution of the c-field density averaged over space and stochastic realisations, compared to its initial value n(t = 0) ≡ n bg for various driving frequencies ωD. Also shown are the equilibrium densities for the initial thermal state (eq th , lower purple horizontal line) and the final-parameter part-superfluid system eq c (upper blue horizontal line) respectively, found by fixing g(t) = gmax (∆g = 0.9) and g(t) = gmin (∆g = −0.9), and the equilibrium density n0(t) = µ/g(t). The slowest quench and condensate equilibrium data are reduced by a factor of 4 for visual aid. For a given driving frequency, the phase lag between the i th trough of ∆g and the i th peak of the average density is denoted φ n i lag . The vertical grey dashed line throughout all three panels highlights the time of the first trough in ∆g. (c) Evolution of the average number of vortices, Nv , for the same quench frequencies and equilibria as above. dependent deviation from criticality which is constrained to oscillate between +α (disordered) and −α (ordered). This quench protocol can be easily visualized in Fig. 1(a). Positive values of ∆g indicate that g > g c , with the system in the disordered (thermal) phase, while ∆g < 0 indicate that g < g c , which initiates the (quasi)condensate formation process. For better understanding, in what follows we also briefly compare our dynamical periodic quench results to the corresponding equilibrium regimes in the two limiting cases of our quench cycle.

III. DYNAMICAL RESULTS
The system dynamics under periodic quenching can be characterized in terms of its density, vortex number, occupation spectrum and coherence, to which we next turn our attention.

A. Density and Vortex Delayed Oscillations
We start our analysis by considering the impact of the periodic interaction strength modulation shown in Fig. 1 on density and vortex number.
Specifically, Fig. 1(b) shows the evolution of the average modulation n − n bg of the c-field density n , compared to its initial value n bg , averaged over space and stochastic realisations, for a range of driving frequencies, while Fig. 1(c) shows the corresponding average evolution of the vortex number in the system. It is wellknown that above the BKT phase transition (∆g > 0) the phase of the system is random, corresponding to an exponentially-decaying phase correlation function: this incoherent state can be thought of as a system with a large number of 'free' vortices, whose exact number is fixed by the system size, grid resolution, atomic mass, temperature and interaction strength (the number itself is not important). As we scan the system across different equilibrium configurations, we find that the number of vortices at equilibrium decreases abruptly across the BKT critical region, reaching a very low, or even zero value as the system approaches its pure superfluid limit (g g c ) [26,69,75]. Let us first consider the limiting (idealized) cases of extremely slow and rapid driving: For extremely slow driving, we expect the system to evolve adiabatically through corresponding equilibrium states, oscillating between the initial purely incoherent state with n = n bg which is filled by vortices (here N v ∼ 160), and a higher value ( n − n bg ≈ 80µm −2 ), which corresponds to the equilibrium density value at g = g min = 0.1g c , and T = 50nK, with such state having N v = 0. The timescale for such periodic change would be set by the slow drive frequency 2π/ω D , in the limit ω D → 0. In fact, the equilibrium density is simply given by the time-dependent homogeneous solution n 0 (t) = µ/g(t), which the idealized case ω D → 0 follows in Fig.1(b). The corresponding averaged densities and vortices in the two limiting cases are shown Fig.1(b)-(c) by the horizontal purple (initial incoherent state with g = g max ) and blue (part-superfluid state with g = g min ) lines.
In the opposite extreme limit of very rapid driving (ω D → ∞), the system has no time to adjust to the driving parameters, and remains in an effective incoherent steady-state close to the initial state, i.e. its density and vortex number closely mimic the flat behaviour of the purple horizontal lines, as shown here for ω D /2π = 500Hz.
Neither of those regimes is intrinsically interesting to study, so our numerical analysis focuses on a broad range of intermediate driving frequencies, from the quasiadiabatic to the dynamically-driven regimes identified below.
The slowest quench analyzed in detail in this work is shown by the red curves in Fig. 1(b)-(c), and corresponds to a drive frequency ω D = 2π × (1/40) Hz. This has been selected such that both the density and vortex number closely resemble the adiabatic case (with N v ≈ 0), but strictly speaking the frequency of the drive is such that the state reached at g min is not yet fully equilibrated, consistent with the state likely to arise under experimental driving conditions: we label such state here as quasi-adiabatic. We have explicitly verified that driving the system with such ω D for half a period (i.e. up to t = π/ω D ), and subsequently removing the drive and waiting sufficient time for the system to relax indeed recovers the corresponding equilibrium state at g min .
For such slow quenches, the maximum system density is achieved almost exactly at each local minimum ∆g(t), while density minima temporally coincide with maxima in ∆g(t). In other words, the density modulation oscillations are in phase with the underlying drive. Moreover, the vortex number rapidly decreases acquiring its minimum value at the local ∆g minimum when the density modulation is also maximised.
As the oscillation frequency increases beyond the adiabatic limit, the overall dynamical evolution becomes faster. To best characterize the underlying dynamics, we consider the dynamical density and vortex number modulation as a function of time scaled to the drive timescale 2π/ω D : this enables evolution graphs to be compared against the same ∆g(t) graph of Fig. 1(a). In such scaled time units, the evolution is comparatively slowed down with increased driving frequency. As a result of the faster driving, the system no longer has sufficient time to adjust to the same density maximum within each cycle, with the oscillating density modulation going out of phase (lagging behind) the external periodic driving, by an amount φ ni lag , corresponding to the phase lag of the i th trough in ∆g(t). The local density maximum is thus reached at a time after ∆g(t) has reached its minimum, when the interaction strength is g min < g(t) < g c .
This effect becomes pronounced with faster driving, which leads both to a lower local density maximum being reached, and a longer time delay between driving and density growth. For intermediate external driving frequencies, the first density maximum reached is lower than the subsequent ones, with the system evidently requiring multiple driving cycles for the temporally-local density maximum to plateau to its overall maximum attainable value under periodic driving. Moreover, the rapid evolution implies that the overall observed density modulation amplitude per cycle is decreasing with increasing ω D , so that -as the density maximum is gradually increasing -the corresponding density modulation minimum is no longer constrained to coincide with the initial density value: such features are most evident when comparing ω D /2π = 1/2Hz (blue) and ω D /2π = 5/3Hz (brown) curves in Fig. 1(b). This effect becomes more pronounced with higher values of ω D , until the driving becomes so rapid that the system has practically no time to adjust.
Such dynamical behaviour can be better understood by looking at the corresponding evolution of mean vortex number N v in the system. As the driving frequency increases, the rate of decrease of vortices in a given driving cycle becomes significantly lower, with the local vortex minimum being significantly delayed from the corresponding ∆g minimum. The delay of the first local vortex minimum φ V1 lag as measured from the ∆g minimum at t = π/ω D is in fact significantly pronounced, and occurs clearly after the corresponding density modulation maximum (which was itself found to phase-lag behind the external driving). This can be understood in terms of the additional phase-ordering time required for the annihilation of a vortex-antivortex pair, which is a competing effect to the external drive. For faster quenches the gradual vortex number decrease occurs over multiple quench cycles and does not necessarily reach N v ∼ 0 as evident from the brown quench in Fig. 1(c). (In the idealized limiting case of extremely rapid quenches (pink), the system has insufficient time to react, exhibiting only a marginal periodic decrease in vortex number).
The above discussion has identified an evident difference in the evolutions of density and vortex number during the driven quench, identifying two rather different behaviours with respect to the driving frequency. This points to the existence of a critical (resonant) frequency which separates these two regimes. We can identify a characteristic relaxation, or 'resonant' frequency through the identification of a characteristic timescale in the system. A relevant timescale in the system is [32] This timescale corresponds to a relaxation time towards the symmetric equilibrium above the critical point, while below the critical point it marks a timescale for dynamical instability, through exponential growth of small fluctuations. This timescale in turn implies a characteristic system frequency ω 0 , defined as This frequency is marked across Fig. 2 (and relevant subsequent figures) as a vertical red dashed line, and clearly marks the boundary between quasi-adiabatic driving (ω D < ω 0 ) and dynamically-driven systems (ω D > ω 0 ). The role of this intrinsic 'resonant' frequency on momentum spectrum and coherence is further highlighted below.
To characterize the above behaviour quantitatively, we next consider the mean phase lag of both density and vortex evolution as a function of the driving frequency ω D , as shown in Fig. 2(a). The average vortex number phase lag is greater than the density phase lag for all driving frequencies considered, in keeping with observations from Fig. 1. A phase lag of φ lag /2π = 0.25 corresponds to the external driving parameter crossing of the BKT threshold, which the density phase lag never exceeds, implying the maximum density per cycle is achieved during the coherent phase. However, the minimum in N v occurs later whilst in the thermal phase, suggesting that the vortex number is not in equilibrium before the next cycle occurs.
Given that the vortex number decay exhibits oscillatory behaviour, and that for a large range of probed driving frequencies it actually requires multiple such oscillations before the vortex number decreases to its steadystate (potentially non-zero) values, it is also interesting to characterize the time taken for the vortex to decay to a near-zero final vortex number. Similarly, the average density also requires several cycles to saturate to its maximum value. To quantify these effects, we have here chosen arbitrarily to define a timescale, n cycle , as the number of quench cycles at which the measured quantity reaches 90% saturation of its steady-state values. This corresponds to the time the vortex number first reaches 10% of its initial value, or the density reaches 90% of its final maximum value.
Apart from the very slow quenches, Fig. 2(b) shows a near-linear increase in the number of cycles required for the system to reach the steady-state as a function of driving frequency. Moreover, even though it may take multiple cycles to reach that value, the linear-like curve implies a near constant duration of time to saturation (2πn cycle /ω D ) for ω D ω 0 . This figure thus confirms the two distinct regimes either side of ω D ∼ ω 0 probed here in the interesting range 0.05 (ω D /ω 0 ) 10. For relatively slow driving ω D < ω 0 , the vortex number and average density saturate within a single quench cycle, suggesting maximum coherence has been reached within one cycle and the vortex number and density evolution are only dependent on the time during a cycle when they are measured, and not the total number of cycles completed. However, when driving faster than ω 0 , the system saturates after multiple quench cycles, with extreme limits of (ω D /ω 0 ) 10 (not shown here) leading to only very small variations from the initial (incoherent) state.
We now turn our attention to the evolution of the momentum spectrum as a function of different driving frequencies. The mean phase shift for the density (φ n ) and vortex number (φ V ) between the i th peaks of ∆g and the peak average density, as a function of ωD. Error bars are smaller than the marker size. Points corresponding to the quench frequencies considered previously are coloured according to the legend of Fig. 1. (b) Average number of oscillation cycles, n cycle , between the first quench and density (vortex number) reaching 90% of its maximum value (10% of its initial value). The relaxation frequency ω0 (Eq. (6)) is marked by a vertical dashed red line, with the grey-shaded area corresponding to the 'quasiadiabatic' regime ωD < ω0. Horizontal black lines in (a) and (b) respectively denote the values of φ lag /2π = 1/4, marking the transition from the observed phase lag occurring within a quarter of the driving cycle (i.e. while the system is within the 'coherent' phase), and n cycle = 1, marking the transition between vortex number counts that are independent and dependent on cycle number. For ωD/ω0 100 (not shown here) the extremely modest response of the system to the drive leads to practically no change in density, or phase.

B. Momentum Evolution and Condensate Mode Dynamics
The density analysis considered above focussed on the entire classical field region, which encompasses both the condensate-which requires macroscopic occupation of the lowest lying momentum mode-and those higher-lying macroscopically occupied modes affected by its presence. To quantify the emerging (quasi)condensate dynamics, we thus consider the evolution of the momentum distribution.
Firstly, we define the fractional momentum occupation showing the profile ofn k at different times, given as vertical lines in (a) (i) and (iv). Also shown is the initial t = 0 distribution prior to the initiation of our periodic quenching (hollow purple circles), with error bars indicative of the magnitude of the error in all of the data shown, and the distribution at equilibrium after a half-quench cycle, followed by a constant g = gmin post-quench relaxation period (dot-dashed blue line). The expected behaviour in the two limiting cases is shown by solid grey lines: exponential behaviour for a thermal cloud ∼ e −ξk [with a fitted coherence length ξ = 0.91µm Lx, Ly] (see also Sec. III C [103]) and a k −2 scaling for a superfluid system. (c)(i) Occupation of then k=0 ≡ N0/N mode for the first density peak (circles) and a late density peak corresponding to the time γt * = 0.12s (squares) [see later]. (c)(ii) Total atom number evaluated at the same times as (c)(i). A vertical red line depicting the relaxation frequency ω0 (Eq. (6)) appears in both panels of (c).

of mode k vian
where · N denotes averaging over stochastic realisations. The evolution of the normalized momentum densityn k is plotted in Fig. 3(a) for different driving frequencies from slow (left, ω D /ω 0 ∼ 0.05) through to fast (right, ω D /ω 0 ∼ 10) quenches. For a truly adiabatic evolution, proceeding slowly enough that the system passes through fully equilibrated states during its entire evolution, one would expect to see the periodic emergence of a 'condensate' mode, in the sense of a highly populated k = 0 mode Indeed, this emerges for the slowest quasi-adiabatic periodic quench [panel (a)(i)], whose modes up to k ∼ 0.1µm are maximally populated at times t = (2π/ω D )(2m+1)/2 (where m is an integer), when ∆g has a local minimum.
This regime clearly shows a periodic reversible condensate formation process, qualitatively similar to the behaviour observed in harmonic dimple microtraps cycled across the phase transition [43,44]. As the quench frequency increases beyond ω 0 ∼ 2π × 0.4 Hz we expect coherence to form over multiple quench cycles. Indeed for ω D ∼ 1.25ω 0 [panel (a)(ii)], we find the initial momentum peak after a single quench cycle to be significantly broadened around k ∼ 0, consistent with the existence of mul-tiple highly-occupied modes (and the gradual emergence of a quasi-condensate). Further increasing the driving frequency leads to coherence only appearing after multiple cycles, consistent with earlier findings, due to the more gradual decrease of vortex numbers which-for faster quenches-becomes significant only after multiple driving cycles through the critical point.
Example slices ofn k vs. k (for k > 0), at extreme limits of ∆g are shown in Fig. 3(b) for both a relatively slow (red; corresponding to a(i)) and a relatively fast (green; corresponding to a(iv)) drive. As a guide for the eye, solid grey lines show the following limiting behaviour: a thermal field is expected to follow an e −ξk scaling, where ξ is the coherence length from the first order correlation function g (1) (r) for a thermal cloud [103], and a fully equilibrated 2D Bose gas is expected to exhibit a k −2 scaling [109]. Our findings show that the driven states oscillate between these limiting cases. As a comparison, we show the extracted momentum distribution for the thermal field at t = 0 (purple circles), and of the Bose gas in equilibrium showing the expected k −2 behaviour (blue dash dotted), and find excellent agreement to the expected scaling.
Considering first the slow driving case [case (a)(i)], we note that the system exhibits a monotonically-decayinĝ n k [ Fig. 3(b)] with more than 95% condensate fraction [ Fig. 3(c)] both at the first (solid red line) and subsequent ∆g minima (dotted red line). This is in contrast to the corresponding behaviour at ∆g maxima (dashed red line) which shows decreasing occupation as k → 0.
Considering now a much faster quench [case a(iv)], the first ∆g minimum [solid green line] is clearly not condensed, as evident from the low k = 0 occupation in Fig. 3(c)(i) [green open circle], and still resembles the initial thermal field distribution at t = 0. After many driving cycles, such a rapidly-driven system develops clear evidence of non-equilibrium condensation, with increasing occupation into the lower k modes [dotted and dashed green lines] and ∼ 50% occupation into the k = 0 mode (Fig. 3(c)(i)), distinct from the expected thermal distribution and with little variation between the maximum (dotted) and minimum (dashed) cases of ∆g. Remarkably, for fast quenches the system maintains a coherent fraction regardless of the drive.
In all cases, it's worth noting that the spectrum at large k acquires a power-law behaviour with an exponent resembling ∼ −1 [104], clearly distinct from both expected thermal and Bose gas scalings. We attribute such numerically observed 'anomalous' scaling to the highly non-equilibrium nature of the driven system, a feature that could be tested in experiments.
The maximum of the condensate fraction as a function of driving frequency is shown in Fig. 3(c)(i), and the total atom number in (c)(ii). The open circles show the condensate fraction measured at the first peak in the total density, whereas the open squares show the condensate fraction measured at an arbitrary peak at sufficiently late times. For faster driving frequencies the condensate fraction increases over time, as a result of the repeated crossing of the phase transition, thus leading to a nonequilibrium steady-state with irreversible growth. This can be attributed to the total time spent above the phase transition, π/ω D , being smaller than the relaxation time τ 0 , and hence a small fraction of the atoms remain Bose condensed after each repeated crossing of the phase transition. This effect disappears both for small driving frequencies, where the condensate fraction follows the equilibrium value expected from an instantaneous evaluation of the dissipative GPE, and for very rapid frequencies (ω D /ω 0 100 -beyond the scale of this graph) which do not allow for any noticeable coherence to form even after multiple quench cycles. In this limit the condensate feels a time averaged drive, g(t) t = g c , and the c-field remains incoherent.
Interestingly, the smallest condensate fraction occurs around ω D ≈ ω 0 , when the competing driving and relaxation are in resonance with one another, leading to a highly non-equilibrium partly-superfluid state, an effect we explore in more detail in Sec. III D.

C. Evolution of Phase Coherence
We now discuss the effect of driving on the spatial coherence of the system.
For a homogeneous 2D system, the spatial correlation function, or normalised first order correlation function, defined by [109]: is known (at equilibrium) to exhibit a transition from exponential (above the critical region) to algebraic-order decay (below the critical region), a characteristic of the BKT phase transition [26,53,69,75,105,106]. We have indeed verified such behaviour in the steady-state profiles for the two limiting equilibrium cases g = g max and g = g min . For a thermal cloud g (1) (r) ∼ e −r/ξ , where ξ is the coherence length [see also Fig. 3(b)]. Here we are instead interested in investigating the effect of the periodic driving on the emerging maximum spatial phase coherence. With this in mind, we perform a detailed analysis of the system coherence when ∆g = −0.9. To account for the fact that the system does not necessarily relax after a single driving cycle (for faster driving), we consider the spatial phase coherence at two times. These correspond to: (a) the end of the first driving half-period (t p1 ) when the system reaches its minimum ∆g value after a single phase transition crossing -this is shown in Fig. 4(a). This corresponds to the same point in the driving cycle irrespective of ω D , thus amounting to different physical times. (b) We also measure g (1) at a fixed absolute late time for all ω D , corresponding to a sufficiently long time such that all considered cases have already acquired their maximum steady-state values. Specifically: Error bars are located where g (1) is measured at a peak in the average density. Data for ωD/2π = 1/40 (red) and 5/33 (grey) Hz extend far beyond γt = 0.12s, with corresponding red/grey bands denoting error bars over the total range covered. Plotted error bars indicate one standard deviation of g (1) (r).
we investigate g (1) (r) at the ∆g trough at, or just after the time γt * ≈ 0.12s: such a time corresponds to the first trough in the limiting case of quasi-adiabatic driving (ω D = 2π(1/40) Hz = 0.05ω 0 ) and the 60 th trough, at which a non-equilibrium steady-state has been reached for the ultrafast quench ω D = 2π × 5Hz = 10ω 0 ]. This is shown in Fig. 4(b). Fig. 4(a) shows that for the quasi-adiabatic driving [ω D = 2π × (1/40)Hz], the system exhibits very strong spatial coherence at the first ∆g minimum, which is perfectly fit with an algebraic-decay function g (1) (r) ∼ r −α , where α = 6.8 × 10 −3 . This correlation is actually identical in all subsequent peaks with the same driving. Although this is highly coherent, consistent with the complete absence of vortices in most numerical realizations, we nonetheless note here once again that this system is not yet fully equilibrated. This is because the emerging state is found to exhibit a momentum spectrum distinct from the expected k −2 equilibrium spectrum, whereas a k −2 spectrum does emerge upon allowing this same system to relax (without further driving) for a significantly longer timescale at ∆g = −0.9.
Increasing the driving frequency leads to the establishment of less coherence in the system after a single cycle, with coherence nonetheless building up gradually after multiple cycles. This can be seen by comparing the computed g (1) (r) between Figs. 4(a) and (b). Although very fast quenches lead to very small coherence after a single drive half-cycle, probing the correlation function instantaneously at g = g min g c nonetheless still facilitates a very good algebraic-decay fit (but with a higher value of α, implying a highly non-equilibrium state). Thus, in order to compare with the maximal coherence imparted at steady-state, we analyze all plots here in terms of an algebraic-decay fit.
The increase of the (algebraic) decay rate of g (1) (r) with increasing ω D seen in Fig. 4(a) should come as no surprise, since the build-up of significant coherence requires the decay of all vortices in the system, whereas all probed quenches except the quasi-adiabatic (red) one still have at least some vortices present at the first halfcycle [see Fig. 1(c)]. Remarkably, even the presence of a few vortices is enough to significantly destroy the system's overall phase coherence, as evident from the blue curves [ω D = 2π × (1/2) Hz] in Fig. 4(a) and Fig. 1(c).
An interesting and rather distinctive feature arises when looking at corresponding late-time peaks shown in Fig. 4(b), after allowing the periodically-driven system to reach its non-equilibrium steady-state, following a sufficiently large number of drive cycles, and repeated crossings of the phase transition in both directions. We have already shown earlier that for sufficiently fast driving, the vortex number curve decreases periodically in an oscillating manner, but with a rapidly decaying envelope, requiring multiple cycles for a significant decrease of its mean vortex number N v . As a result, coherence builds up gradually in such systems, and this was to be expected. The interesting emerging feature here is that at steady-state, there is an optimum driving frequency, above which the final acquired system coherence (at one of the ∆g peaks) at steady-state grows over longer distances again. This is easily seen by comparing Fig. 4(a)-(b). Focussing at the late time γt * ∼ 0.12s, we see that increasing the frequency (in Hz) from ω D /2π = 1/40 (red) through 5/33 (brown) to 1/2 ∼ O(ω 0 /2π) (blue), g (1) (r) decays faster spatially with increasing drive frequency. However, driving the system even faster than ω 0 periodically across a sufficiently large number of phasetransition cycles leads to an increase in the system phase coherence, with the spatial coherence at the maximum distance (half the box length) of the ω D /2π = 5 Hz (green) being comparable to that of ω D /2π = 5/33 Hz (grey), which is driven 33 times slower (and whose coherence does not change between consecutive peaks, having already saturated at its maximal coherence at the end of the first half-cycle).
Such behaviour becomes more evident in Fig. 4(c), when examining the temporal evolution of the value of the phase correlation function near the edge of the system, at g (1) (r = 40µm). Slow quenches (ω D /2π < 0.5 Hz) allow maximal coherence to be effectively established already after a single quench half-cycle: this is indicated by the red/grey horizontal lines/bands in Fig. 4(c). For faster driving than that, the system initially reacts rather slowly, with the coherence growth rate gradually picking up around γt ∼ 0.01s, eventually saturating at a higher value (well before the γt ∼ 0.12s used for the comparison in Fig. 4(b)). Interestingly the extra growth of coherence after multiple cycles is very small around the frequency ω D /2π ∼ 1/2 Hz (blue curve), but increases noticeably for quenches few times faster than that, as evident from the enhanced values shown by the yellow, brown and green curves.

D. Resonances with intrinsic timescales
The behaviour of the spatial correlation function is further characterized in comparison to the relaxation frequency ω 0 = 2π × 0.4 Hz in Fig. 5. Here we plot (a) the steady-state value (t = t * ) of g (1) near the box edge (r = r ), (b) the fractional increase of its value at t * , compared to that at its first peak, (c) the value of the power-law decay exponent α, and (d) the density phase lag γτ n delay = 2πφ n lag /ω D (a) (corresponding to Fig. 2(a), but now plotted in terms of average delay time which offers a different perspective). This figure clearly highlights the importance of the critical driving frequency, corresponding to 'resonant' driving distinguishing the slow and fast driving regimes. This resonant driving frequency is well explained in terms of the frequency ω 0 . Looking at the various subplots, we easily infer that maximal coherence is achieved in the quasi-adiabatic regime (leftmost, red points in all subplots), with g (1) (r = 40µm, γt * = 0.12s) ∼ 0.96 and α = 6.8 × 10 −3 , which is, as expected for such a lowtemperature state with T /T BKT ∼ 0.15, much less than the value α = 0.25 occurring at equilibrium for g = g c [107].
As the system is driven faster, the late-time coherence at the system edge rapidly decreases (down to ∼ 0.11) [ Fig. 5(a)], with such value being already reached after a single quench half-cycle: the latter can be inferred by looking at the fractional change in the value of the correlation function at r from the first peak up to a late converged peak at t * . Correspondingly, faster driving leads to an increasing delay (phase lag) of density growth, as shown in Fig. 5(d). Interestingly, an algebraic power-law fit g (1) (r, t) ∼ r −α of the correlation function at γt * ∼ 0.12s -when the system clearly possesses a non-zero occupation of the k = 0 mode -leads to a significant increase in the value of α, much exceeding the equilibrium value of (1/4). This latter observation pro- vides strong evidence of the non-equilibrium nature of the achieved steady-state, and points to the strong interplay between the system's external driving across the phase transition (quantified in ω D ), and its ability to adjust (quantified in ω 0 ) to an averaged behaviour across two vastly different equilibrium states located well above and well below the critical region.
When driving the system even faster than that, the system loses its ability to adjust rapidly, as it starts experiencing an averaged quench effect. The overall steadystate coherence at the trap edge starts increasing again with increasing frequency [ Fig. 5(a)], with an increase of as much as 50% over its initial (but rather low) acquired value at the first quench half-cycle [ Fig. 5(b)]. As the system has less and less time to react to the external driving, it effectively only exhibits a time-averaged effect, and so becomes less non-equilibrium, consistent with a decreas-ing value of α with faster driving ω D > ω 0 [Fig. 5(c)]. As a result of the rapid driving, the overall density evolution time-delay now also decreases significantly (in 'absolute' time γt [Fig. 5(d)]), which is testament to the very tiny increase of the system density facilitated.
Nonetheless, it is important to note that this delay time is actually a significant fraction of the drive time, as evident from Fig. 2(a), revealing that the density delay time for the fastest quench considered here is 20% of the quench period, while the vortex delay time [not shown here, see Fig. 2(a)] exceeds 60% of the period.
Much faster driving than that (ω D /ω 0 100) [not shown here] leads to a practically monotonic timeaveraged evolution and a minor decrease of about 30% in vortex number, with an associated increase of less than 10%in the population of the k = 0 mode, indicating a vortex-filled, low-coherence, non-equilibrium state. Moreover, driving at an idealized 'infinite' frequency ω D /ω 0 ∼ 1000 was found to lead to practically no change to the system's initial configuration.

IV. CONCLUSIONS
We have considered the dynamical response of a homogeneous two-dimensional ultracold Bose gas under periodic quenches of its interaction strength through the Berezinskii-Kosterlitz-Thouless phase transition at a driving frequency ω D . We have identified an intrinsic system response frequency ω 0 and demonstrated that resonant driving leads to a highly non-equilibrium state exhibiting only limited coherence growth, compared to when driving the system on either side of this resonance. Focusing on the most interesting regime of driving at a frequency within one order of magnitude from the resonant value, we characterized the system response in terms of the driving frequency, by analysing the evolution of densities, vortices, condensate fractions, spectra and coherence.
Specifically, we identified two distinct driving regimes of experimental relevance: Driving at a frequency ω D much smaller than ω 0 gives rise to quasi-adiabatic dynamics, with maximum coherence achieved after a single quench half-cycle, even if the system has not had sufficient time to fully equilibrate in momentum space (with the limiting case of extremely slow drive, ω D → 0, corresponding to the adiabatic regime). Under (quasi)adiabatic conditions it is possible to analytically describe the evolution of the average density.
In the opposite regime, driving the system faster than the system response frequency leads to a dynamicallydriven non-equilibrium system, whose coherence grows gradually through multiple quench cycles. Interestingly, such driving can lead to a significant condensate fraction and enhancement in coherence: for example, even for large driving frequencies ω D /ω 0 ∼ 10, the coherence can grow up to 50% over the range of the box, despite only achieving < 1% coherence after the first quench. As expected, quenching much faster than that (e.g. ω D /ω 0 100) leads to practically no generation of coherence, with the system entering an incoherent, vortex-filled steady-state.
Our study extends earlier works on cyclic phase transition crossing, and may help guide future non-equilibrium driven quench experiments and lead to the generation of interesting non-equilibrium steady-states of ultracold atoms with partial superfluidity.
Data supporting this publication is openly available under a Creative Commons CC-BY-4.0 License on the data.ncl.ac.uk site [108].