Non-equilibrium Dynamics of Vortices in Two-Dimensional Quantum Gases: Determining the Dynamical Scaling Region Using the Mahalanobis Distance

When a two-dimensional system undergoes a rapid quench from a disordered to an ordered phase, it does not order instantly but instead relaxes towards equilibrium over time. During this relaxation, the dynamical scaling hypothesis predicts that the length scale of ordered regions increases, with later patterns statistically similar to earlier ones except for a change in global scale. Quantum gases are one of many systems in which such out of equilibrium behaviour is predicted to occur. Here we present a method for systematically testing when dynamical scaling is taking place, by quantifying the similarity of the rescaled two-point correlation function over time using the Mahalanobis distance. Data on the velocity field of a quantum fluid, generated from point vortex simulations, are used to illustrate the application of this method.


Introduction
People search for evidence of dynamical scaling in systems such as the Ising model of magnetic spins, binary alloys [1,2] and Bose-Einstein Condensates (BECs) [3].A key prediction of the dynamical scaling hypothesis is that the correlation function, G(r, t), between two points separated by a distance r, changes in time, t, only by a rescaling in length scale by L c (t), where L c (t) is a correlation length that depends on time only.

G(r, t) = G(r/L c (t)) .
( To test if this prediction is met, the correlation functions at different times are rescaled by L c (t) and plotted on the same axes.If they collapse onto a single curve, this is evidence of dynamical scaling.Typically, this collapse is checked visually; however, a more rigorous, algorithmic test is preferable in order to determine the times at which dynamical scaling applies.In this paper we outline an approach that makes use of the Mahalanobis distance [4] to quantify the deviation of a curve from the mean, eliminate outliers and identify the scaling region.When dynamical scaling holds, the dependence of correlation length on time is predicted to be a power law of the form: where z is the dynamical critical exponent and is expected to have a value of two for a two-dimensional system with a non-conserved order parameter [1,2].A recent numerical study of coarsening dynamics, following the rapid quench of a Bose gas into the ordered phase, found clear evidence of dynamical scaling [3].The values of z found were close to two but dependent upon dissipation, parameterized by a dimensionless constant, γ.This study used the projected Gross-Pitaevskii equation (PGPE) and the stochastic projected Gross-Pitaevskii equation (SPGPE) to model the evolution of the wavefunction ψ(r, t), both in the absence and the presence of dissipation.
The authors observed growing regions of constant phase following the quench, as quantum vortices and antivortices annihilated, as shown in Fig. 1.They found the correlation functions at different times using the wavefunction and observed them to collapse, within a given scaling window, onto a single curve once rescaled by L c (t).Their method for identifying the scaling window was to judge by eye initially and later check that the maximum discrepancy between rescaled correlation functions was sufficiently small.An earlier study [5], simulating a similar system using a truncated Wigner approach, also found evidence of universal scaling behaviour in the time evolution of the single-particle momentum spectrum.They selected a scaling window by inspecting the rescaled spectra and found a value of z = 1.8(3) (avoiding a strongly anomalous non-thermal fixed point associated with vortex clustering).To illustrate our method to test for the collapse of rescaled correlation functions, we also consider a two-dimensional quantum fluid.However, rather than simulate the system using the Gross-Pitaevskii equation (GPE) [6,7], we use the point vortex model (PVM) instead [8,9].The PVM is often used to simulate dynamics in two-dimensional quantum fluids [10][11][12] such as those formed in ultracold atomic gases [13,14].Although vortices in atomic gases are not point-like and the gas itself is compressible, the PVM can give reasonably accurate predictions for vortex behaviour.It also benefits over the more directly applicable GPE because only the coordinates of the vortices need tracking and not the quantum mechanical wavefunction at every grid-point.Consequently, it is computationally less expensive and makes simulations of larger systems possible.The structure of this paper is as follows: Section 2 describes how data for the out of equilibrium dynamics of a two-dimensional quantum fluid is generated using a dissipative PVM.Section 3 details how a two-point correlation function is calculated and rescaled and how the scaling region is determined using a method based on the Mahalanobis distance.Section 4 summarises the main findings and suggests some possibilities for future work using this method.
2 Generating Data Using the Point Vortex Model (PVM)

The PVM on an Infinite Plane
The PVM applies to two-dimensional incompressible fluids with zero viscosity and assumes that vorticity is confined to delta-function-like points [8].Each vortex is advected by its local fluid velocity, which is made up of contributions from each of the other vortices.Consequently, the equations of motion for vortex i, of circulation, κ i , at position (x i , y i ), in a system of N v point vortices are: Where κ j and (x j , y j ) are the circulations and positions of all the vortices other than i and x ij = x i − x j , y ij = y i − y j and r ij = (x ij , y ij ).This result may also be derived from the Hamiltonian [9]: in which the double sum adds contributions from all pairs of vortices.
In the case of quantum fluids, vortex circulations are quantised and take only the values of ±h/m, where h is the Planck constant and m is the mass of an atom.When modelling with the GPE, the natural units to use are the healing length units.Distances are measured in terms of the healing length, ξ = ℏ/ √ mµ, and times in terms of, τ = ℏ/µ, where µ is the chemical potential.In these units, the circulations are ±2π.We refer to a vortex of circulation +2π as a vortex and a vortex of circulation −2π as an antivortex.

A Dissipative PVM in a Doubly-Periodic Square Box
To model the dynamics of the vortices we use a doubly-periodic square box of length L = 2048ξ.To account for interactions between vortices across the periodic boundaries we use the following equations [15]: As above, the circulation of vortex i is κ i but its coordinates, (x i , y i ), are rescaled by 2π/L to give (x i , ỹi ), so that in these units the box has a length of 2π.xi − xj and ỹi − ỹj are written as xij and ỹij .The time units are also scaled to give The first sum is over all vortices except for vortex i and the second is over the integer m.Since the term in this second sum falls with increasing m, an excellent level of precision is achieved by curtailing the sum to be from m = −10 to m = +10.We use a Dormand-Prince 8th Order Runge-Kutta algorithm [16] to solve these equations of motion, with an absolute tolerance of 1 × 10 −8 .This gives us a nondissipative velocity, v i = (dx i /dt, dỹ i /dt), for each vortex.To include the effects of dissipation we add a dissipative component at right angles to this to give the overall velocity, including dissipation, [11].The values of γ we use are {2 −8 , 2 −7 , . . ., 2 −1 , 1, ∞}, where γ = ∞ is the purely dissipative case, in which A final, key aspect of vortex dynamics that we add to the simulation is the annihilation of vortices and antivortices that come into too close proximity.This is observed in experiments and in GPE simulations and is the main mechanism by which phase ordering takes place.To include this in our dissipative PVM, at each time-step we remove any vortex-antivortex pair that are within 1ξ of each other.

Generating Initial Conditions
The initial arrangement of vortices and antivortices affects the dynamics that follow.For example, an isolated vortex-antivortex pair, known as a dipole, will move in a straight line, perpendicular to the line joining the vortices.Alternatively, a cluster made up only of vortices (or only of antivortices) will tend to rotate about its centre, so long as no other vortices are near.The effect of dissipation on dipoles is to draw the constituent vortex and antivortex closer together over time, whereas it tends to scatter clusters apart.Therefore, when generating multiple realizations of the initial conditions for repeat runs of the point vortex simulations it is vital to ensure they have simlar distributions of dipoles, clusters and free vortices.To achieve this, we note that the Hamiltonian, H □ , for a system of point vortices in a doubly-periodic square box [15], gives an unambiguous parameter for measuring vortex configurations: In this equation all symbols have the same meanings as in Eq. 5 & 6.The double sum at the beginning adds contributions from all pairs of vortices and, as before, an excellent level of precision is achieved by curtailing the sum over m to be from m = −10 to m = +10.Calculating this value for the initial configuration and dividing the result by the number of vortices, N v , gives the initial energy per vortex E 0 /N v .A completely random configuration of vortices and antivortices would have E 0 /N v = 0, increasing the amount of clustering of like-signed vortices increases E 0 /N v and pairing vortices and antivortices into dipoles makes E 0 /N v more negative.
To create initial conditions with specified values of E 0 /N v we first scatter 2069 vortices and 2069 antivortices in the box.Then we employ a biased random walk algorithm in which the positions of vortices are updated at random and these updated positions accepted only if they bring E 0 /N v closer to the desired value.Once the desired values is achieved to within a 1% tolerance, further steps are undertaken in order to thermalize the configuration.In each step, two vortices are randomly selected and moved by different, random amounts.If E 0 /N v remains within the tolerance the step is accepted.Throughout all of this process, changes that bring vortices closer than 1ξ are rejected.

Constructing the Velocity Field and Correlation Function
In order to demonstrate dynamical scaling, a two-point correlation function is generally calculated using a relevant scalar order parameter.In GPE simulations this may be the quantum mechanical phase or the fluid density, in the Ising model it is the magnetization, in our case we use the speed of the fluid.To find this we first construct the velocity field from the coordinates of the point vortices.We consider a 100 × 100 grid of points (x g , ỹg ), filling the square box of length 2π and find the velocity of the fluid at each grid point, v = (v x, v ỹ ).This is done by summing the contributions, at that location, of all of the point vortices, whilst accounting for the periodic boundary conditions [15]: Eq. 8 & 9 are almost the same as the non-dissipative point vortex equations of motion, Eq. 5 & 6, except for two differences; the sum is over all of the vortices and the coordinates, (x g , ỹg ), are those of a grid point and not those of a vortex.This is because the point vortex equations arise from considering the contributions of all other vortices to the fluid velocity at that point and asserting that, in the absence of dissipation, the vortex moves at the local fluid velocity.See Fig. 2 (a)-(c) for examples of the velocity field constructed from the point vortex coordinates in this fashion.
Using the velocity, v = (v x, v ỹ ), we find the corrected speed at each grid point, v c , from the difference between the speed, |v|, and the mean speed taken across the whole grid, ⟨|v|⟩: This is the scalar order parameter we use to find the two-point correlation function and ultimately provides the basis for demonstrating dynamical scaling.
The correlation function, G, for the corrected speed, v c , between two points separated by a displacement (x, ỹ), is given by: Where the sum over i is across all grid points in the x direction, the sum over j is for all those in the ỹ direction and the denominator ensures normalization such that G(0, 0) = 1.The Wiener-Khinchin theorem [17] provides a more succint, and computationally more efficient, formula in terms of Fourier transforms: To convert from G(x, ỹ) to a one-dimensional G(r) we consider nested annuli of thickness, δr, where δr is the distance between adjacent grid points and r2 = x2 + ỹ2 .We then find G(r) by calculating the mean value of G(x, ỹ) across all the grid points in the annulus from r to r + δr.Finally, we multiply the distances by L/2π to return to the original length scale used for our system and recover G(r).We repeat this for each time-step, for all 100 stochastic realizations and for the whole ensemble of initial conditions.

Rescaling the Correlation Function and Testing for Dynamical Scaling
The theory of dynamical scaling predicts that the correlation function, G(r, t) should vary only by a change in length scale as time progresses.If this is the case, finding a suitable length at each time-step and rescaling the distances by this length, should cause the plots of G(r) at different times to collapse onto a single curve.
The process we use to rescale the correlation function for a given set of initial conditions is as follows: 1. Consider one time-step and find the mean values of G(r) across the 100 stochastic realizations (from here on G should be understood as the average over these realizations).2. Use cubic spline interpolation to find the radius at which G(r) = 0.1.We define this as the correlation length at this time-step, L c (t). (We use 0.1 for our data as it gives sufficient points, at both smaller and larger radii, to ensure a reliable interpolation).

Rescale the values of radius by dividing them by L c (t).
Fig. 3 Representing curves as points in an Ng-dimensional Euclidean space.(a) Five curves G(r) plotted on a grid of three points r 1 , r 2 , r 3 .(b) The same curves represented by points in a threedimensional Euclidean space.Similar curves, such as those shown with red circles, green squares, blue diamonds and magenta triangles, become points that are clustered close together; whilst a very different curve, like the one shown with orange stars, becomea a point very far from the others.To identify an outlier, we need to quantify the deviation, ∆, of a suspected point from the mean position of all the other points.We use the Mahalanobis distance, in the Ng-dimensional Euclidean space, to do this.

Repeat the steps 1-3 for all time-steps.
This gives arrays of data for G and r/L c at each time-step.Simply plotting G against r/L c allows us to check visually if the curves for different time-steps collapse.However, a more systematic approach is desired to determine the start and end of any period of dynamical scaling in a robust and reliable fashion.The method we implement to test for collapse of the different curves is: 1.For every time-step, use spline interpolation to find the values of G on the same uniform grid of values of r/L c .Use N g points in the range from the smallest value of r/L c (t) for any time-step up to r/L c (t) = 5.By this point the values of G should be very close to zero in all cases, so continuing further is unnecessary.We use N g = 40.2. For each time-step there are now N g values of G. Therefore each time-step can now be represented by a point in an N g -dimensional Euclidean space: where, G i is the vector position representing the correlation function at timestep i and (G i 1 , G i 2 , ..., G i Ng ) are the coordinates in each of the dimensions (i.e. the values of G on the grid of N g radii).The problem is now one of determining which time-steps are outliers and which are not, see Fig. 3. 3. Find the mean, µ, and the covariance, C, of G i across all of the time-steps.µ is a vector of length N g and C is a matrix of size N g × N g .

4.
The squared deviation of a point, G i , from the mean, µ, normalized by the covariance between the data from the orginal curves, is given by the Mahalanobis distance squared, M D 2 i , [4]: Calculate M D 2 i for each time-step and order them from largest to smallest.Then, to identify and eliminate outliers: 5. Choose the time-step with the largest value of M D 2 i and recalculate µ and C without it (i.e.treate it as an outlier).6. Recalculate M D 2 i for the suspected outlier using the new µ and C. If the value is below a cutoff (see discussion below), this timestep is not an outlier and the process is stopped.If it is above the cutoff, this timestep is identified as an outlier and removed from future consideration.

Recalculate M D 2
i for all the other points, with the previously identified outlier(s) excluded and order these values from largest to smallest.8. Identify the time-step with the new largest value of M D 2 i and repeat steps 5-7.9. Continue iteratively until a suspected outlier has M D 2 i less than the cutoff, at which point all outliers have been identified.10.All remaining timesteps are considered to overlap.The scaling region is then the longest stretch of uninterrupted timesteps not containing any outliers.See Fig. 4 for an example.
The choice of cutoff for determining outliers is a slighlty subtle one.If the values of G i at each radius are independent and normally distributed then M D 2 i would qualify as a χ 2 statistic and the critical χ 2 value for a suitable significance level would be an appropriate choice.However, our data are not normally distributed, so instead we use a cutoff of 4×N g , where N g is the number of grid-points.The logic behind this decision is that for the case in which G i at each radius varies independently of every other radius, the covariance matrix is diagonal and M D 2 i becomes a sum of similar terms: So, if a hypothetical curve has every point one standard deviation from the mean, then M D 2 = Ng k=1 1 2 = N g , whereas if each point is two standard deviations from the mean then M D 2 = Ng k=1 2 2 = 4N g .We consider this a reasonable choice of cutoff, whilst recognizing that our values of G i at each radius are not, in general, independent.However, a lower value could be used if a stricter definition of collapse is preferred.The cutoff is the one free parameter in this approach and plays a similar role to the significance level in hypothesis testing.Fig. 5 (a) shows how the choice of cutoff affects the percentage of timesteps identified as within the scaling region and consequently the values of dynamical critical exponent determined from the scaling region data.those time-steps identified as within the scaling region are plotted as solid red lines, getting darker as time progresses, and those outside of the scaling region are plotted as gray dashed lines.In the rescaled plots, (b), all curves in the scaling region overlap.(c) Shows the data for correlation length against time plotted as grey crosses, the bars show the standard error in Lc calculated across the 100 stochastic realizations at each time-step.The red line is a power law Lc ∼ t 1/z , with z = 1.985±0.003determined using a bootstrap method [18] in which we resample Lc(t) 100 times, conduct a power law fit each time and find the mean and standard error in z across the resamplings.(d) Illustrates the time-steps identified as outliers as black blocks in the bottom line and those in the scaling region as red blocks on the top line.

The Dynamical Critical Exponent
To calculate values for the dynamical critical exponent, z, we conduct a least squares fit of the power law L c ∼ t 1/z to the data L c (t) for just the times within the scaling region.By identifying the scaling region in a robust fashion, with a suitable choice of cutoff, prior to conducting the fit, we can be more confident that the predicted power law relationship is valid and the value of z found is an accurate measure of the dynamical critical exponent of the system.Fig. 5 (b) shows how the value of z, determined in this way, varies with the level of dissipation, γ.

Conclusions
We have described a systematic method, based on the concept of Mahalanobis distance, for identifying outliers when comparing graphs of rescaled correlation functions.By eliminating those times for which the curves do not collapse onto one single curve, Fig. 5 Determining the dynamical critical exponent, z, for E 0 /Nv = 4. (a) The black crosses (lefthand axis) show how the value of z varies with the choice of cutoff for γ = 1 (errorbars are not shown as they are either too large, for cutoff < 3Ng, or too small, for cutoff ≥ 3Ng).The red diamonds (right-hand axis) show how the % of timesteps identified as within the scaling region varies.Using a cutoff of less than 3Ng is problematic as only a tiny fraction of timesteps are within the scaling region and the fit is unreliable.(b) The black crosses (left-hand axis) show how the value of z varies with the dissipation, γ, for a cutoff of 4Ng (errorbars would be too small to see on this scale).The red diamonds (right-hand axis) show the corresponding % of timesteps in the scaling region.A large scaling region is identified for all γ and z ≈ 2 but dependent upon γ, a qualitatively similar result to that reported in [3].In all cases we determine z using the same bootstrap method described in Fig. 4.
our approach gives a reliable identification of the scaling region.We have also demonstrated its application to the phase-ordering behaviour a quantum fluid following rapid quench, using data generated using a dissipative PVM.This produced that correspond qualitatively and quantitatively with those from other studies that make use of GPE simulations and identification of scaling region (with subsequent checks) [3,5].In particular, we found values of dynamical critical constant that are close to two but depend on dissipation.By applying this method we can be confident that the data used to determine the dynamical critical exponent came only from times during which the system displays scaling behaviour.Further work could include an application of this approach to the correlation functions for a wider ensemble of initial conditions, to the results of GPE simulations or to studies of non-equilibrium behaviours of systems outside of the field of quantum fluids.

Declarations
• Funding This work was supported by a Lady Bertha Jeffreys PhD Studentship.
• Competing interests The authors have no competing interests to declare that are relevant to the content of this article.• For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

Fig. 1
Fig. 1 Evolution of the phase of the field ψ of a BEC in the case of zero dissipation, reproduced from [3] under CC-BY-4.0international license 1 .Distances are quoted in terms of the healing length, ξ = ℏ/ √ mµ, and times are given in terms of τ = ℏ/µ.These are the characteristic length and time scales used when solving the Gross-Pitaevskii equation, where m is the atomic mass, ℏ is the reduced Planck constant and µ is the chemical potential.The BEC was modelled using the projected Gross-Pitaevskii equation (PGPE) in a doubly-periodic square box of width ≈ 363ξ.The BEC was created in a far from equilibrium configuration with a high density of quantised vortices and antivortices.It then "relaxes", in this case through annihilation of vortex-antivortex pairs.(a), (b) and (c) correspond to times t ≈ {200, 2000, 20000}τ respectively.The open white squares show the positions of vortices (positive circulation), the filled black squares show the positions of antivortices (negative circulation) and the background colour indicates the phase.

Fig. 2
Fig. 2 Fluid velocity calculated from point vortex configurations and two-point correlation functions for the initial conditons γ = 1, E 0 /Nv = 4. (a)-(c) Open (white) circles show the locations of vortices and filled (blue) circles show the locations of antivortices.The pale red streamlines show the fluid velocity at the times t = {4200, 13200, 32700}τ .(d)-(f) Are the two-point correlation function G(r) calculated from the corrected speed at the same times.

Fig. 4
Fig. 4 Identifying the scaling region and determing z for the initial conditions γ = 1, E 0 /Nv = 4. (a) Shows the two-point correlation function G(r) for each time-step.(b) Shows G(r) with the distances rescaled by Lc, where Lc is the value of r when the unscaled G(r) = 0.1.In both (a) and (b)those time-steps identified as within the scaling region are plotted as solid red lines, getting darker as time progresses, and those outside of the scaling region are plotted as gray dashed lines.In the rescaled plots, (b), all curves in the scaling region overlap.(c) Shows the data for correlation length against time plotted as grey crosses, the bars show the standard error in Lc calculated across the 100 stochastic realizations at each time-step.The red line is a power law Lc ∼ t 1/z , with z = 1.985±0.003determined using a bootstrap method[18] in which we resample Lc(t) 100 times, conduct a power law fit each time and find the mean and standard error in z across the resamplings.(d) Illustrates the time-steps identified as outliers as black blocks in the bottom line and those in the scaling region as red blocks on the top line.