Crossover in the dynamical critical exponent of a quenched two-dimensional Bose gas

We study the phase ordering dynamics of a uniform Bose gas in two dimensions following a quench into the ordered phase. Within the classical field methodology, we perform numerical simulations exploring the crossover between dissipative and conservative evolution. Regardless of the dissipation strength, we find clear evidence for universal scaling, with dynamical critical exponent $z$ characterising the growth of the correlation length. In the dissipative limit we find growth consistent with the logarithmically corrected law $[t/\log(t/t_0)]^{1/z}$, and exponent $z=2$, in agreement with previous studies. For decreasing dissipation we observe a smooth crossover to the conservative limit, where we find strong numerical evidence for the expected growth law $t^{1/z}$, but with anomalously low exponent $z \approx 1.7$. We show that this lower exponent may be attributable to a power-law vortex mobility arising from vortex-sound interactions.

We study the phase ordering dynamics of a uniform Bose gas in two dimensions following a quench into the ordered phase. Within the classical field methodology, we perform numerical simulations exploring the crossover between dissipative and conservative evolution. Regardless of the dissipation strength, we find clear evidence for universal scaling, with dynamical critical exponent z characterising the growth of the correlation length. In the dissipative limit we find growth consistent with the logarithmically corrected law [t/ log(t/t 0 )] 1/z , and exponent z = 2, in agreement with previous studies. For decreasing dissipation we observe a smooth crossover to the conservative limit, where we find strong numerical evidence for the expected growth law t 1/z , but with anomalously low exponent z ≈ 1.7. We show that this lower exponent may be attributable to a power-law vortex mobility arising from vortex-sound interactions.
Introduction-A many-body system quenched from a disordered to an ordered phase has long been a topic of interest in nonequilibrium physics. Following the quench the system relaxes toward a new equilibrium configuration via a process of domain coarsening, with an associated growth of the correlation length L c (t). The dynamical scaling hypothesis posits that at sufficiently late times the system should approach a statistically invariant state in which L c becomes the only relevant length scale. In this state, the correlation length is predicted to grow ∼ t 1/z , where z is the dynamical critical exponent [1]. In the classical theory of phase ordering kinetics, coarsening is described in terms of the dynamics and annealing of topological defects, with conservation laws and dimensionality playing a key role in determining z [1]. Extensive numerical studies in two-dimensional (2D) systems such as Ising [2,3] and XY [4][5][6] models have provided broad support for this simple physical picture.
In recent years, ultracold Bose gases have become an established platform for the exploration of nonequilibrium dynamics in a quantum setting. Owing to their exquisite tunability, experiments have been able to use these gases to probe physics such as the Kibble-Zurek mechanism [7][8][9][10] and quantum turbulence [11][12][13][14], as well as scale invariant dynamics following a quench, similar to the scenario described above [15][16][17][18]. In this last context, the concept of nonthermal fixed points [19][20][21][22][23][24][25] has emerged as a powerful theoretical description of scaling behaviour in which topological defects appear to play a less crucial role than in phase ordering kinetics [25].
Theoretical studies have addressed coarsening following an instantaneous quench in 2D bosonic systems such as binary [26] and spinor [27,28] condensates and drivendissipative systems [29][30][31][32]. However, in the apparently simple case of a quenched scalar 2D Bose gas there remain open questions regarding the link between coarsening behaviour and conservation laws in the dynamics. These stretch back to the well-known classification of dynamical universality classes established in Ref. [33]. For Bose gases, there exist both conservative and nonconservative classical field descriptions of the dynamics [34]. The former conserve energy and particle number [35,36]. The latter include dissipation [37][38][39][40][41]; they have no conserved quantities, and in the dissipative limit they reduce to a purely relaxational time-dependent Ginzburg-Landau equation. Hence, for nonconservative dynamics the relevant dynamical universality class would appear to be Model A [33]. There is general theoretical and numerical agreement that for Model A, z = 2 with logarithmic corrections [5,6,[42][43][44]. This scaling can be explained in terms of defect dynamics, under the assumption that the vortex density and correlation length are related, n v ∼ L −2 c [1]. However, as noted in Ref. [33], the coarsening behaviour of a Bose gas with conservative dynamics is theoretically less tractable. Previous numerical studies in this scenario measured contradictory exponents z ∼ 1 [35] and z = 1.8(3) [24,45], while the related conservative XY model yielded z = 1.80(5) [46]. On the other hand, simulations of coarsening in a nonconservative Bose gas have yielded exponents z = 2.0(2) [47] and z = 1.9(2) [24,45], although the weak dissipation included in such works places these results between the purely dissipative Model A and the conservative limit. Additionally, while none of these studies included logarithmic corrections, they also did not rule out their relevance. As such, an overall picture remains elusive. What the precise value of z is for conservative dynamics, and how z varies with the dissipation strength, remain important open questions.
In this Letter, we revisit the problem of coarsening in a 2D scalar Bose gas after an instantaneous quench. By varying the dissipation rate γ in the stochastic projected Gross-Pitaevskii equation (SPGPE) we explore the crossover between the dissipative limit and conservative dynamics. In the dissipative limit, we observe an exponent consistent with z = 2 with logarithmic corrections, in good agreement with previous results for Model A. For decreasing dissipation we observe a crossover in the value of z, which reduces to z ≈ 1.7 in the conservative limit. We analyse the vortex motion and find that this decrease may be attributable to a power-law vortex mobility resulting from vortex-sound interactions.
Simulations-To describe a Bose gas at finite temperature, we adopt a classical field model [34], where L GP = −( 2 /2m)∇ 2 + g|ψ| 2 . In this model, the gas is represented with a complex scalar field ψ(r, t), which includes contributions from all highly occupied single-particle modes of the system, up to some chosen cutoff in the single-particle energy spectrum. The gas is considered to be in contact with a thermal reservoir at temperature T and with chemical potential µ (corresponding to the above-cutoff atoms), with which it can exchange both energy and particles. The projection operator P ensures that no population is transferred outside the chosen subset of single-particle modes during the evolution, while the constants m and g correspond to the particle mass and the 2D interaction strength, respectively. The dimensionless dissipation rate γ controls the strength of the coupling between the system and the bath, and dW(r, t) is a complex Gaussian noise term satisfying dW * (r, t)dW(r , t) = (2γk B T/ )δ(r − r )dt. With α = 1, this model is known generally as the SPGPE [48], and good quantitative agreement with ultracold Bose gas experiments typically requires γ 1 [7,34,[49][50][51]. In the limit γ → 0, the coupling is removed, and Eq. (1) reduces to the projected Gross-Pitaevskii equation (PGPE), for which both the energy E = ( 2 |∇ψ| 2 /2m + g|ψ| 4 /2)dr and norm N = |ψ| 2 dr are conserved under time evolution. For γ 1, on the other hand, the first term on the right hand side of Eq. (1) becomes negligible, and the dissipative Model A is recovered. In practice, we set α = 0 to access Model A, and use α = 1 otherwise.
In this work, we consider a system in a doubly periodic square domain of size L × L. The properties of the thermal bath (µ and T ) are held fixed, and only the dissipation rate γ is varied. Two quench protocols are used, depending on the choice of γ. For γ > 0 we begin with ψ = 0 and instantaneously switch on the reservoir coupling at time t = 0, forcing the classical field density to grow nonadiabatically. For γ = 0, a more careful choice of initial condition must be made, because of the conservation of energy and particle number. To facilitate direct comparison between these cases, the mean energy-and particle-densities for the γ = 0 initial states are chosen to be equal to their ensemble-averaged values, =Ē/L 2 andn =N/L 2 , in the γ > 0 system after equilibration (i.e., at t → ∞). This is achieved by initiating the wavefunction as a populated disk of radius k d in wavenumber The populations n k are chosen to be uniform and to ensure the correct mean densitȳ n. The phase φ k of each mode is initially randomised, and a Powell minimisation algorithm [52] is subsequently used to adjust the phases to achieve mean energy-density¯ . Both types of quench initialise the system far from equilibrium, with a high density of quantised vortices and antivortices. The single-particle modes for our chosen geometry are plane waves satisfying |k| < k cut for some wavenumber cutoff k cut . To prevent aliasing, the wavenumber cutoff is set to k cut = π/(2∆x) (half the Nyquist wavenumber of the grid), with a numerical grid spacing of ∆x ≈ 0.7 ξ, where ξ = /(mµ) 1/2 is the healing length. This choice of k cut results in an occupation of ∼ 1 particle per mode at the cutoff [34]. For the γ = 0 initial condition, we use k d = 0.2 k cut . We set T ≈ 2.1 µ/k B and g ≈ (1). The γ > 0 SPGPE is solved numerically using XMDS2 [53], while the γ = 0 PGPE is parallelised on nVidia Tesla V100 GPUs using CUDA [54], allowing for the significantly longer evolution time required in this limit. The system and ensemble sizes are respectively chosen to be L ≈ 262 ξ and N = 400 for γ > 0.1, and L ≈ 363 ξ and N = 256 for γ ≤ 0.1.
Analysis-Following the quench, the Bose gas begins to relax and the vortex number decreases via the annihilation of vortex-antivortex pairs. This process is illustrated in Fig. 1, where the phase, arg{ψ}, is shown at three times during the evolution of the PGPE. The vortex density is seen to decrease over time, allowing regions of phase coherence to develop.
As a measure of the spatial coherence of the field at a given time, we calculate the first-order correlation function, The angular brackets in this expression correspond to an average over both stochastic realisations and the co-ordinate r . The scaling hypothesis asserts that a universal form for this correlator, G(r, t) = G eq (r)F(r, t), should emerge at late times [1].
Here, G eq (r) = G(r, t → ∞) is the equilibrium correlation function, and F is a scaling function that should have the form F(r, t) = F(r/L c (t)), with F(0) = 1. The correlation length L c (t) in this expression is defined as the average distance over which equilibrium correlations have been established at time t, corresponding to the average size of the phase domains. We calculate F(r, t) from our simulations by measuring both G(r, t) and G eq (r), where the latter is obtained from a temporal and ensemble average of G(r, t) once the system has equilibrated (the onset of equilibrium can be inferred, for example, from the stabilisation of the k = 0 mode population at late times). Below the Berezinskii-Kosterlitz-Thouless transition [56][57][58], we expect that G eq (r) ∼ r −η for r ξ, where 0 ≤ η(T ) ≤ 0.25 is a temperature-dependent exponent [59]. We find that η ≈ 0.06 for our parameters, using a method described in Ref. [60]. We measure the correlation length L c (t) as the radial distance satisfying F(L c , t) = F 0 , where we set F 0 = 0.5 [55]. This choice assists in excluding discretisation effects at small scales and finite-size effects at large scales. (a,i) Scaling behaviour with α = 1, and dissipation rates γ = 0 (a) and γ = 1 (b). Ensemble averaged scaling function F(r, t) plotted against radial distance r (i), both before (inset) and after (main frame) rescaling by the correlation length L c (t). The time at which each curve has been sampled is denoted by the colorbar in the insets, and a grey dot signifies the threshold F 0 = 0.5 used to define L c . Evolution of the mean correlation length L c (t) (ii) and vortex density n v (t) (iii). In columns (ii,iii), the chosen scaling region is highlighted, and the best power-law fit to the data within that region is shown as a black dashed line (offset for visibility). In these four panels, the left insets are histograms showing the distribution of exponents z measured from fits to subsets of the data within the highlighted region, while the right inset shows the compensated correlation length (ii) and vortex density (iii) as a function of time (horizontal axis same as for main frames).
Similarly, we should extract L c (t) over a scaling window in time that both suppresses finite-size effects at long times and excludes initial transients. In practice we end the window as soon as L c (t) > L/4 in any one of the simulations in the ensemble; this stringent condition generally corresponds to an average L c (t) of ∼ L/10 at the end of the window. We start the window as early as possible while ensuring that there are minimal deviations from the unique scaling function F(r/L c (t)). To check for finite-size effects, we have also repeated the γ = 0 simulations with L ≈ 726 ξ using a smaller ensemble of N = 64 trajectories.
Results-The evolution of F(r, t) is displayed in Fig. 2(a,i) and (b,i) for γ = 0 and γ = 1, respectively, and a collapse of the data onto a single unique curve is evident over the time windows shown. In Fig. 2(a,ii) and (b,ii), the evolution of L c (t) is shown for the same two γ values, and in both cases the data are well described by L c (t) ∼ t 1/z within the highlighted scaling window. To measure the exponent z in each case, we fit a power-law to the data across all possible subintervals of ≥ 8 consecutive points within the scaling window. This yields a distribution of z values characterising the statistical uncertainty associated with temporal variations in the scaling of L c (t) (left inset of each frame). We measure z and its uncertainty as the mean and standard deviation of this distribution, respectively.
For these two cases, we obtain z = 1.68(4) (γ = 0) [z = 1.69 (3) with L ≈ 726 ξ] and z = 2.19(3) (γ = 1). As an illustration of the goodness-of-fit, we also plot the compensated correlation lengthL c (t) = L c (t)/L fit c (t) [right insets of column (ii)], which remains close to unity within the highlighted interval. In fact, both curves are well described by the same power-law growth even outside of the chosen windows, suggesting that the scaling may be continuing beyond our estimated regions. It is interesting to note the difference in shape between the L c (t) curves; while there is a long delay before scaling emerges in the γ = 0 case, the strong dissipation in the γ = 1 system causes the power-law growth to begin almost immediately.
As the correlation length grows and phase coherence develops across the system, there is a concurrent decrease in the mean vortex density n v (t). It therefore seems reasonable to expect that dynamical scaling should also manifest in the decay of these defects [24,61,62]. Under the assumption that the vortex distribution is random and approximately spatially uniform at all times, we can predict power-law scaling of the form n v (t) ≈ L −2 c (t) ∼ t −2/z . The mean vortex density is shown in Fig. 2(a,iii) and (b,iii), and a fit to the data within the highlighted window is found in the same way as for the L c (t) curves. The left and right insets, as in column (ii), correspond respectively to the histogram of z values from repeated fits,  (5) with L ≈ 726 ξ] and z = 2.32(7) for γ = 1; both of these values are slightly larger than those obtained from the corresponding fits to L c (t). We note, however, that measuring z from n v (t) is a less rigorous approach, because the vortices are not guaranteed to remain uniformly distributed. Indeed, within the scaling windows we measure negative nearest-neighbour vortex correlations, indicating a tendency toward dipole pairing of vortices and antivortices [55]. This effect is stronger for the larger γ.
We have repeated the analysis shown in Fig. 2 for a range of γ values, and find equally clear evidence for dynamical scaling in all cases. As above, the exponent z is measured from fits to both L c (t) and n v (t) in each case. We find a smooth crossover from z ≈ 2.3 in Model A (γ → ∞) to z ≈ 1.7 in the PGPE (γ = 0), as shown in Fig. 3 [63]. This clearly shows that the coarsening behaviour of the system changes as one crosses between Model A and conservative dynamics. We have performed several additional simulations and analyses to verify the robustness of this result [55].
In the context of Model A dynamics, it has long been argued that z = 2 with logarithmic corrections [43,44]. These corrections are predicted to modify the scaling such that L c ∼ [t/ log(t/t 0 )] 1/z [5,6,30,47], and hence L c ∼ t 1/z only in the limit t t 0 for some microscopic timescale t 0 . We find that a fit to the above expression using our Model A L c (t) data gives an exponent of z = 2 if we choose t 0 = 0.5 /µ, although we note that the fit quality is no better than an uncorrected power-law. Nonetheless, this establishes consistency between our results and the predicted behaviour for Model A. Log-corrected exponents fitted using the same t 0 at other values of γ are shown in Fig. 3, although existing arguments for log-corrections only apply to Model A [64].
To elucidate the origin of the measured crossover in z, we assemble movies of the evolution with vortex tracking [55]. We observe a distinct qualitative change in the dynamics of vortices as γ is varied, suggesting that the crossover in z may be explained in terms of vortex motion. We start by assuming that: (i) the correlation length grows at a rate determined by the mean vortex velocityū v ∼ dL c /dt; (ii) the characteristic vortex velocity isū v ∼ µ v (L c )/L c , where µ v (L c ) is the vortex mobility. We therefore must have dL c /dt ∼ µ v (L c )/L c , which can be integrated to yield z. A simple description of vortex motion in a Bose gas is a weakly damped point-vortex model [55,65,66]. In this idealised case, µ v = const. and hence L c ∼ t 1/2 , i.e. z = 2. The same argument with mobility µ v ∼ 1/ log(L c /ξ) results in the log-corrected expression quoted above, also with z = 2 [43,44]. A possible origin for z < 2 in the conservative limit is a power-law mobility µ v ∼ L ε c , yielding L c ∼ t 1/(2−ε) , i.e. z = 2 − ε [67]. Our γ = 0 data allows us to test this possibility. We directly measureū v from the tracked vortices and extract ε = 0.27(9) [55], in good agreement with the measured z ≈ 1.7. We have also checked that assumption (i) is reasonably satisfied. This supports a direct relation between z and the vortex dynamics, but raises the question of the origin of the power-law mobility. To investigate further, we perform weaklydamped point-vortex simulations, using as initial conditions the vortex configurations extracted from the γ = 0 simulations at the beginning of the scaling window. We find that these simulations exhibit scaling, but with ε ≈ 0 (i.e. µ v ≈ const.) as expected [55]. This analysis strongly suggests that the powerlaw mobility results from vortex-sound interactions, which are absent in the point-vortex model. Turning to our Model A data, we are unable to reliably measure the mobility directly, because the slow drift velocity resulting from vortex interactions is overwhelmed by a fluctuating fast velocity arising from the noise. This issue also arises in the SPGPE for γ 0.1.
Finally, we have also compared our Model A results to the vortex-dipole model of a temperature quench developed in Refs. [68,69], and find them to be broadly consistent. We observe the sampled dipole distribution function Γ(r, t) to follow the predicted form L −ζ c (t)F d (r/L c (t)), with scaling function F d , for ζ ≈ 4 [55].
Conclusions-We investigated the coarsening of a scalar 2D Bose gas following an instantaneous quench into the ordered phase. By varying the dimensionless dissipation rate γ, we explored the crossover between purely relaxational (Model A, γ → ∞) and conservative (γ = 0) dynamics. Our results for Model A are consistent with dynamical critical exponent z = 2 with logarithmic corrections, the generally agreed result in the literature. Our central result is that for decreasing dissipation rate γ we continue to observe universal scaling in time, but with a smooth reduction in exponent towards the value z ≈ 1.7 in the conservative limit. We found evidence that this deviation from z = 2 may be attributed to an anomalous power-law vortex mobility that arises from interactions between vortices and sound waves. For 2D Bose gas experiments, our results imply that measured exponents will depend on what dissipation rate is appropriate to a particular experimental set-up; however, γ 1 is typical [7,34,[49][50][51], suggesting that deviations from z = 2 may be experimentally relevant. In future, it will be interesting to further investigate the source of the power-law vortex mobility identified here.
Data supporting this publication is openly available under a  For the γ > 0 movies, the ψ = 0 initial condition is used, and the mean density is seen to grow rapidly, plateauing by t ∼ 100 /µ. Note that there is no input temperature parameter T for the γ = 0 simulation, but we have chosen the effective equilibrium temperature to be the same as for γ > 0 by restricting the number-and energy-density, as described in the main text. In all movies, tightly bound thermal dipoles [S1] are seen to spontaneously appear in the superfluid from time to time, surviving only briefly before annihilating again.
We note that changing γ gives rise to a change in the direction of preferred vortex motion. For γ 1, nearest-neighbour vortex dipoles travel approximately perpendicular to their separation vector, and can traverse many times the average interdefect distance before annihilating. By contrast, dipoles experience mutual attraction when γ 1, and hence vortices rarely travel beyond their closest neighbours before annihilating. We additionally observe that for γ 1 the vortices rapidly evolve to a state where mutual attraction becomes overwhelmed by fluctuations.

MEASUREMENT OF THE VORTEX MOBILITY
We compare our γ = 0 (projected Gross-Pitaevskii equation; PGPE) results with those of a weakly dissipative point-vortex (PV) model. In this model, the velocity u k of vortex k is: where is the conservative equation of motion for a configuration of point-vortices at locations {x k , y k } in a square domain of sidelength L with periodic boundary conditions [S2]. Here, x jk = x j − x k (likewise for y jk ), and γ PV 1 is a phenomenological damping parameter that models the loss of energy to sound waves. As initial conditions, we use the vortex positions extracted from our 256 PGPE simulations at t = 3400 /µ, which is the beginning of the scaling window identified in the main text (≈ 130 vortices remain at that time). We then solve the above equations using a semi-implicit integration scheme, with the inner sum truncated to −3 ≤ n ≤ 3 [S3]. To best model the PGPE dynamics, we remove vortex-antivortex pairs if they come within ξ of one another.
Under time evolution, the point-vortex density is seen to decay. We find that a choice of γ PV = 0.01 results in almost immediate power-law scaling n v (t) ∼ t −2/z , in agreement with the PGPE. The vortex configuration also remains similar to the PGPE throughout the evolution, as evidenced by C nn (t) (see next section and Fig. S2). However, the exponent is significantly different from z = 1.74(3), as measured from the PGPE data. A fit within the window 6 × 10 3 ≤ µt/ ≤ 10 4 gives z = 1.99(1), consistent with z = 2 (other values of γ PV delay the onset of power-law scaling, but eventually also result in z ≈ 2). The vortex density from the PV simulations is shown in Fig. S1(a), alongside the PGPE data within the scaling window [the data is the same as in Fig. 2(a,iii)  To extract the exponent ε described in the main text, we measure the mean vortex velocityū v as a function of time, and compare it to n v (t) [since it is difficult to extract L c (t) from the PV simulations, we instead assume that n v ∼ L −2 c ]. With a power-law mobility µ v (L c ) ∼ L ε c , we expectū v ∼ n (1−ε)/2 v . For the PV data,ū v is calculated by averaging Eq. (S1) over all vortices. For the PGPE data, we instead measure the velocity of each vortex by tracking defects between adjacent time samples (separated by ∆t), and calculate a finite difference velocity u k (t) = [r k (t) − r k (t − ∆t)]/∆t. Although simplistic, this measurement appears to capture the overall behaviour ofū v . The results are shown in Fig. S1(b), in terms of the sound speed c s = µξ/ . Power-law fits give ε = 0.27 (9) for the PGPE data, consistent with z ≈ 1.7, and ε = 0.04(4) for the PV data (in the same temporal window stated above), in agreement with z ≈ 2. We therefore conclude that the exponent z < 2 measured in the PGPE is consistent with a power-law vortex mobility, with ε > 0. Furthermore, since ε ≈ 0 (and z ≈ 2) in the PV model, the exponent ε > 0 must be a beyond-point-vortex effect.

ANALYSIS OF THE VORTEX CONFIGURATION
As an indicator of the vortex configuration at a given time, we calculate the nearest-neighbour correlator, where N v is the total number of vortices, s j = ±1 is the circulation sign of vortex j, and s nn j is the circulation sign of its nearest neighbour. A value of −1 ≤ C nn 0 indicates a vortex configuration predominantly paired into dipoles, 0 C nn ≤ 1 indicates clustering of same-sign vortices, and C nn ≈ 0 corresponds to an approximately random vortex distribution [S3].   . S2 shows the evolution of the mean C nn (t) as measured from the γ = 0 and γ = 1 simulations (for comparison with Fig. 2 of the main text). Throughout the evolution, C nn < 0 in both cases, indicating that the vortices are in the dipolar regime. Interestingly, the vortices are distributed quite differently within the respective scaling windows, as evidenced by the significantly different values of C nn .

DIPOLE PAIR DISTRIBUTION
In Refs. [S4, S5], an analytic model was introduced for describing the relaxational dynamics of a two-dimensional superfluid following a temperature quench. In this description, the system is represented as a vortex-dipole pair distribution function Γ(r, t), whose evolution is governed by a Fokker-Planck equation. The distribution Γ(r, t) represents the probability of finding a vortex dipole of separation r in the vortex configuration, and can be integrated to give the mean vortex density, n v (t) = 2 Γ(r, t)d 2 r (the prefactor here accounts for the two vortices per dipole). In Ref. [S5], it was predicted that the dipole distribution should obey a scaling form Γ(r, t) ∼ L −ζ c (t)F d (r/L c (t)) with scaling function F d , assuming that L c ∼ t 1/z . The exponent ζ is predicted to depend on the initial and final temperatures of the quench; but in particular, ζ = 4 for quenches from the Berezinskii-Kosterlitz-Thouless critical temperature.
Given that relaxational dynamics are assumed in Refs. [S4, S5], our γ → ∞ (Model A) simulations are the closest point of comparison to the predictions stated above. Although there is an inherent ambiguity in assigning a configuration of vortices into a set of dipoles, we nonetheless attempt a measurement of Γ(r, t) using the same prescription as in Ref. [S6]. We first rank in increasing order the distance between every possible pair of opposite circulation vortices in the system (taking into account the periodic boundary conditions). Beginning with the smallest pair, we proceed sequentially through this list, assigning a pair as a dipole only if neither of its constituent vortices has already been assigned to another dipole. In this way, we assemble a unique list of N v /2 pairs, where every vortex has been paired exactly once. We then construct a histogram n d (r, t) of all dipole sizes r in the system, and average the distribution over all simulations in the ensemble. Finally, we define Γ(r, t) = n d (r, t)/(2πrL 2 ); here, the factor of L 2 converts to a spatial density, while the factor of 2πr accounts for the number of ways of configuring a dipole of size r. Figure S3 shows the resulting distribution Γ(r, t) for times within our identified scaling window (determined in the same way as discussed in the main text). Upon rescaling the data according to the above prediction with ζ = 4, a convincing collapse is obtained. We also observe a power-law tail in the dipole distribution, Γ(r, t) ∼ r −3 , for r L c . The peak visible in the smallest radial bin can be attributed to the aforementioned thermal dipoles that appear throughout the system during the dynamics.

POSSIBLE SYSTEMATIC EFFECTS
Choice of threshold-We find that the choice of scaling function threshold F 0 (used for defining the correlation length) has a slight systematic effect on the power-law scaling of L c (t). In Fig. S4, the critical exponent z measured from power-law fits to L c (t) is shown as a function of dissipation rate γ, with the correlation length extracted using three different choices of F 0 . The fitting is carried out in the same way as described in the main text. Evidently, a larger F 0 gives rise to a slightly larger measured exponent z (on average). As noted in the main text, F 0 should be chosen to minimise both discretisation effects at small scales and finite size effects at large scales. As such, we do not explore thresholds outside of 0.3 ≤ F 0 ≤ 0.7 here.
Lack of finite size effects-We check for finite size effects by repeating our γ = 0 simulation in a system of size L ≈ FIG. S5. Measured critical exponent z as a function of dimensionless dissipation rate γ. The blue circles and red diamonds are the same as in Fig. 3 of the main text, while the additional data (labelled) corresponds to a quench at γ = 0 with a doubled system size (black circle/green diamond), and multiple quenches at γ > 0 using our nonzero density initial condition (purple circles/orange diamonds). Where multiple sets of data fall onto a single γ value, the points have been symmetrically offset along the horizontal axis for clarity. 726 ξ, for an ensemble of N = 64 trajectories. The measured exponents z from fits to both L c (t) and n v (t) are shown in Fig. S5 (as a black circle and green diamond, respectively) alongside the data from Fig. 3 of the main text. To within our estimated uncertainties, we obtain the same exponents in the larger system, suggesting that finite size effects are negligible in our simulations.
Lack of initial condition dependence-In the theory of phase ordering kinetics [S7], it is generally expected that the precise form of initial conditions used to initiate the coarsening behaviour should play no role in determining the observed scaling. To confirm that this holds for our system, we have repeated our γ = {0.1, 0.5, ∞} simulations using the nonzero density initial condition described in the main text, where it was used for γ = 0 simulations. We do this using a system size of L ≈ 363 ξ, and ensembles of N = 64 (N = 256) trajectories for γ = {0.5, ∞} (γ = 0.1). In Fig. S6, we compare the evolution of the scaling function [column (i)], correlation length [column (ii)] and vortex density [column (iii)] for the two initial conditions, with γ = {0.1, 0.5, ∞} in rows (a,b,c), respectively. The procedures used to identify the scaling windows and perform fits to the data have been carried out as described in the main text. For all γ values, the rescaled F(r, t) data [main frames of column (i)] are almost indistinguishable between the two initial conditions. Likewise, the observed power-law scaling in columns (ii) and (iii) appears to be almost unaffected by the initial condition, despite substantial differences in the curves at early times. We note that for the zero (nonzero) initial density configurations, the curves approach the scaling law from steeper (shallower) evolution. This provides strong evidence that the system has reached a universal scaling regime. Figure S5 displays the critical exponents z measured from  Comparison of scaling behaviour between the zero (orange/triangles) and nonzero (blue/circles) density initial conditions, for dissipation rates γ = 0.1 (a), γ = 0.5 (b), and γ → ∞ (c). Ensemble averaged scaling function F(r, t) plotted against radial distance r (i), both before (inset) and after (main frame) rescaling by the correlation length L c (t). The time at which each curve has been sampled is denoted by the color bar in the insets, and a grey dot signifies the threshold F 0 = 0.5 used to define L c . Evolution of the mean correlation length L c (t) (ii) and vortex density n v (t) (iii). In columns (ii,iii), the scaling window is highlighted in the appropriate colour, and the best power-law fit to the data is shown as an orange dotted (blue dashed) line for the zero (nonzero) density initial condition (offset for clarity). The initial and final times of the scaling window, t i and t f , can be read off the horizontal axis in columns (ii,iii).
fits to both L c (t) and n v (t) in these simulations (purple circles and orange diamonds, respectively). In all cases, these exponents are consistent (within our estimated uncertainties) with the exponents measured for the same γ value using the zero density initial condition.