Universal dynamics in the expansion of vortex clusters in a dissipative two-dimensional superﬂuid

A large ensemble of quantum vortices in a superﬂuid may itself be treated as a novel kind of ﬂuid that exhibits anomalous hydrodynamics. Here we consider the dynamics of vortex clusters under thermal friction and present an analytic solution that uncovers a new universality class in the out-of-equilibrium dynamics of dissipative superﬂuids. We ﬁnd that the long-time dynamics of the vorticity distribution is universal in the form of an expanding Rankine vortex (i.e., top-hat distribution) independent of initial conditions. This highlights a fundamentally different decay process to classical ﬂuids, where the Rankine vortex is forbidden by viscous diffusion. Numerical simulations of large ensembles of point vortices conﬁrm the universal expansion dynamics and further reveal the emergence of a frustrated lattice structure marked by strong correlations. We present experimental results of expanding vortex clusters in a quasi-two-dimensional Bose-Einstein condensate that are in excellent agreement with the vortex ﬂuid theory predictions, demonstrating that the signatures of vortex ﬂuid theory can be observed with as few as N ∼ 11 vortices. Our theoretical, numerical, and experimental results establish the validity of the vortex ﬂuid theory for superﬂuid systems.


I. INTRODUCTION
A defining feature of quantum fluids is that they exhibit quantized vortices.These stable topological defects have a circulation that is quantized in units of = h/m, where h is Planck's constant and m is the mass of a fluid particle.Despite this key difference from classical viscous fluids, many features of turbulence, i.e., the complex, collective behavior of many vortices, are common to classical and quantum fluids.In threedimensional quantum turbulence, which has been extensively studied in bulk superfluid helium [1,2], examples include the Kolmogorov energy cascade [3,4], the dissipation anomaly [5], and boundary layers [6].More recently, experimental advances in quasi-two-dimensional (2D) ultracold atomic gases [7][8][9][10][11] and superfluid optomechanical systems with thin-film helium [12] have renewed interest in turbulence and vortex dynamics in two dimensions, where markedly different behavior to three dimensions is often observed.Here, 2D quantum fluid Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
One might expect the emergence of classical phenomena from quantum vortex dynamics to follow from Bohr's correspondence principle; provided many quantum vortices of the same sign are bundled together, collectively they should mimic classical vortex tubes.In two dimensions, however, recent theoretical work has shown that a dense system of chiral (i.e., same sign) quantum vortices at large scales may be treated as a kind of fluid in its own right [22].In such a vortex fluid, the dynamics are governed by a hydrodynamic equation that contains anomalous stress terms absent in the standard Euler equation, allowing for phenomena such as analog edge states of the fractional quantum Hall effect [23].This theory was recently extended to describe dissipative effects [24], accounting for mutual friction due to the interaction between the superfluid and a stationary thermal component present in experiments.However, exact solutions to the vortex fluid theory equations, and by extension the collective dynamics of many-vortex systems, are still lacking.Furthermore, an understanding of how the anomalous stresses affect large ensembles of quantum vortices remains unexplored.Despite the relevance of the theory to many experimental systems, neither the conservative nor dissipative vortex fluid theory have been demonstrated experimentally.
Here we consider the dynamics of a 2D chiral vortex cluster within the dissipative vortex fluid theory [24].Generally, solving for the out-of-equilibrium dynamics of many-body systems poses significant challenges.We provide an analytical solution to this theory, demonstrating a universality class in the out-of-equilibrium dynamics of dissipative superfluids.We show that dense vortex clusters evolve into a Rankine vortex (i.e., top-hat distribution) at long times independent of the initial vorticity distribution.This behavior is markedly different to the case of a classical viscous fluid, where not only is the Rankine vortex forbidden, but the expansion of a vortex is governed by ordinary viscous diffusion [25].
To corroborate the theory, we simulate large collections of point vortices and demonstrate that any distribution of vorticity evolves into a Rankine vortex.Beyond the vortex fluid theory, at long times we observe frustrated ordering of the vortices that becomes highly correlated at both short and long distances.Finally, we experimentally observe the emergence of a Rankine vortex in a quasi-2D 87 Rb Bose-Einstein condensate, demonstrating the utility of vortex fluid theory for a system with as few as N ∼ 11 vortices.Our findings establish a connection between the abstract concepts of the vortex fluid theory developed in Refs.[22,24] and their physical realizations.Through our numerical and experimental results, we demonstrate a platform for further experiments investigating the vortex fluid theory.

II. POINT VORTEX MODEL
We consider the motion of N vortices in a homogeneous quasi-2D superfluid characterized by healing length ξ , speed of sound c, and vertical thickness d.We assume that d ∼ ξ so that the longitudinal excitations of the vortex line are suppressed and hence the vortex motion is effectively 2D, with all vortex cores parallel to ẑ [26,27].Each vortex carries singly quantized circulation i = κ i h/m with κ i = ±1.In the limit where the velocity of the vortices is far below the speed of sound (v c) and the intervortex separation is greater than the healing length ( ξ ), we are able to approximate vortices at positions r i = (x i , y i ) as pointlike particles generating vorticity ω(r) = i i δ(r − r i ) and a fluid velocity field u(r) Hence we can describe the motion of the vortices with a point vortex model [28,29].
The incompressible kinetic energy of a 2D fluid can be expressed in terms of the relative vortex positions.In free space (i.e., in the absence of boundaries), the Hamiltonian is where ρ s is the 2D superfluid density, r i j = |r i − r j |, and L is an arbitrary length scale [30].Using Hamilton's equations for the point vortex system, and setting ρ s = L = 1 for simplicity, the velocity of vortex i is expressed in terms of the other vortex positions as where x i j = x i − x j , y i j = y i − y j , and For systems at finite temperature, interactions between the superfluid and thermal component result in the dissipation of energy proportional to the relative velocity of the two components [31].In two dimensions, the thermal component is typically stationary due to viscous clamping in superfluid helium [12] or strong trap anisotropy in atomic gas superfluids [32].The effect of dissipation is described with the equation of motion where the dimensionless mutual friction coefficient γ (typically 1) characterizes the strength of the dissipation.In the context of superfluid helium, Eq. ( 4) and the mutual friction coefficient can be rigorously derived from the interactions between a vortex and a thermal phonon bath, where the mutual friction coefficient is temperature dependent, i.e., γ ∝ T [33][34][35].In the context of ultracold atomic gases and assuming superfluid density gradients are negligible, Eq. ( 4) can be rigorously derived from the dissipative Gross-Pitaevskii equation (dGPE) [36], where is the condensate wave function, m is the atomic mass, V (r) is the trapping potential, g is the two-body contact interaction parameter, and μ is the chemical potential of the system.The dGPE describes a weakly interacting Bose-Einstein condensate (BEC) coupled to a uniform stationary thermal reservoir.The parameter γ in Eq. ( 4) is the same parameter that describes the dissipation in Eq. (5).

III. DISSIPATIVE VORTEX FLUID THEORY
A system containing a large number of well-separated 2D quantum vortices can be viewed as a fluid in its own right and its dynamics described by a set of anomalous hydrodynamic equations [22,24].For a chiral system where all vortices have the same-sign circulation ( i κ i = ±N), the collective dynamical variables are vortex density ρ ≡ i δ(r − r i ) and vortex fluid velocity field ρv ≡ i v i δ(r − r i ).Note that here ω = ρ.
A complete description of the dynamics of dissipative chiral vortex fluids is presented in Appendix A. The full equations of motion are complex and ultimately not needed for the purposes of this paper.Here, we consider only the evolution of the vortex fluid density.For a vortex fluid governed microscopically by Eq. ( 4), the anomalous hydrodynamical equation for the vortex fluid density is where is the material derivative, and γ v × ∇ρ describes transverse convection [24].
In Eq. ( 6) the term ∝ ∇ 2 ρ describes uphill diffusion of ρ, which, in contrast to ordinary viscous diffusion, serves to steepen local vorticity gradients.The damping term −γ ρ 2 strives to suppress regions of high density and, together with the nonlinear term, balances the uphill diffusion to prevent a singular solution.
While Eq. ( 6) is intractable, by considering the competition of terms we assume a uniform vortex fluid density, i.e., ∇ρ = 0.This simplifies Eq. ( 6) to ∂ t ρ = −γ ρ 2 , which has the solution where ρ 0 = ρ(0) is the initial density.This solution describes an expanding Rankine vortex, where the density distribution is uniform within the cluster and zero outside [37].
The Rankine vortex expansion is characterized by the mean radius which shows that the cluster exhibits diffusive-type growth.
The canonical angular momentum of the fluid reads L f = − dr r 2 ω/2 = −9π N r(t ) 2 /8, and hence its dynamics is fully determined by Eq. (8).
The energy of the cluster, calculated from H = ρ s dr |u| 2 /2 with the velocity field of the Rankine vortex u φ (r < r c (t )) = r ρ/2; u φ (r > r c (t )) = N/(2π r), evolves as To demonstrate the stability of the Rankine vortex, we show its solution is robust against density fluctuations within perturbation analysis.Consider a local density perturbation ρ = ρ a + δρ to the universal expansion solution, where ρ a is given by Eq. ( 7), |δρ| |ρ a |, and δρ(|r| → ∞) = 0.For simplicity, we assume that δρ has cylindrical symmetry, namely, δρ(r) = δρ(r).Substituting the perturbed Rankine vortex solution into Eq.( 6) and keeping the leading-order terms in δρ, we obtain where v a = r ρ a /2.The last two terms in Eq. ( 10) become less relevant at long times as the coefficients in front of δρ and ∂ r δρ are proportional to ρ a , which tend to zero in the longtime limit.Hence the evolution of the vorticity distribution will eventually be dominated by a diffusion equation with a negative diffusion coefficient.The only physical solution (i.e., no singularities during the time evolution) is therefore constant density.We have solved Eq. ( 10) numerically for several cylindrically symmetric perturbations and find each decays to zero as the system is evolved (see Appendix A for details).
We ultimately find that the Rankine vortex is an asymptotic solution of the full anomalous hydrodynamical equation [Eq.( 6)] independent of initial conditions.The combination of the damping and negative-viscosity terms in Eq. ( 6) suppress density fluctuations in the cluster and yield the formation of a Rankine vortex that is an attractor of the dissipative dynamics.These results suggest that an initially nonuniform density will evolve towards the universal Rankine scaling solution described by Eq. (7).In classical viscous fluids, an axisymmetric vortex expands via diffusion, and the latetime profile of an isolated line vortex instead tends to the Lamb-Oseen vortex [25].This distinct behavior highlights the fundamentally different dissipation mechanisms between finite-temperature superfluids and classical viscous fluids.

A. Universal expansion regime
We further demonstrate the universality of the Rankine vortex by numerically simulating large ensembles of vortex clusters.In Fig. 1 we show the results of simulating the expansion of five ensembles of N = 1000 vortices according to the point vortex model [Eq.( 4)] with γ = 0.01.The initial conditions are drawn from three axisymmetric states: a top-hat, Gaussian, and ring distributions, as well as two nonaxisymmetric states: a one-dimensional line and a random distribution of multiple clusters.Without dissipation, the top-hat and Gaussian distributions (in the limit of N → ∞) are thermal equilibrium solutions for the vortex distribution [38].The ring distribution and the nonaxisymmetric states are highly nonequilibrium initial conditions [38,39].In each row of Fig. 1, we plot the cluster at a different time through the expansion indicated on the left (we have chosen γ τ = 4.8 × 10 −2 2π R 2 / as a convenient unit of time), with the relative size of the cluster labeled on the right.After sufficient evolution (t ∼ 1000γ τ ), each of the vortex clusters are well within the universal expansion regime.
In the final row of Fig. 1 [t = 1.3 × 10 5 γ τ ], we have plotted the position of the vortices and highlighted dislocations in the cluster to show the emergence of crystallization.It is evident that beyond this local structure, the final states are identical and the universal expansion regime holds regardless of initial condition.In Fig. 2(a), we find excellent agreement in comparing the mean radius of the clusters from the simulations and the analytic solutions of the vortex fluid theory [Eq.( 8 (d) and 1(e)] behave very similarly to the ring distribution to their highly nonequilibrium nature.Thus for the remainder of the paper, we focus only on the top-hat, Gaussian, and ring cases.In the bottom panel of Fig. 2(a), we plot the error between the theory and simulation, defined as ε = ( r VFT − r PVM )/ r PVM , and observe it to approach zero as the cluster expands.We similarly find strong agreement for the energy of the cluster, as seen in Fig. 2(b).
To further emphasize the uniformity of the expanded clusters and indeed the emergence of the Rankine vortex as predicted by the vortex fluid theory, we plot the average radial density (averaged over n = 100 realizations of each distribution) of each ensemble throughout the expansion.Figures 3(a)-3(d) show that the cluster evolves towards uniformity, with the density of each cluster becoming more uniform at different rates.The ring initial condition exhibits a slower timescale in evolving towards a top hat due to its highly nonequilibrium initial state.The specific value of γ we have chosen here is unimportant; γ simply determines how quickly the cluster reaches the regime of universal dynamics, after which it is simply a scaling factor in time.To understand and demonstrate the effect of the dissipation, i.e., the cause of the expansion, we also provide comparisons with the conservative dynamics (γ = 0) for all scenarios in Fig. 9 (see Appendix B).The universal dynamics we observe in the expanding regime is in contrast to the nonergodicity observed in conservative point vortex dynamics, where large ensembles of point vortices can get trapped in asymmetric, nonequilibrium stationary states [40].
At sufficiently late times, however, the cluster begins to crystallize; additional features emerge at length scales comparable to the typical intervortex spacing, going beyond the vortex fluid theory predictions.Vortices at the cluster edge gradually organize into concentric circles, leading to a density  3(d), as demonstrated by the black curve, yields quantitative agreement with the vortex fluid theory.We note these oscillations may be related to the predictions of Ref. [23], where it was found that the superfluid Rankine vortex supports an edge layer with a number of interesting properties, including a density overshoot and soliton solutions with quantized charge.The emergence of the short-length-scale features in the density can better be seen in Fig. 3(e), where we plot the density as a function of radius and time for the top-hat initial distribution.At early times (i.e., t 100γ τ ), the density is uniform.As the cluster evolves, discrete peaks form at the outer edges of the distribution and eventually in the bulk of the vortex cluster.

B. Beyond vortex fluid theory: Crystallization
To investigate this emerging density structure and quantify the ordering and crystallization of the vortex lattice, we calculate the vortex-vortex correlation function, where s is the spatial separation of a pair of vortices, N is the number of vortices contributing to separation s, and θ i is the accessible angle at distance s for a given vortex assuming a circular distribution [41].The correlation functions of the top-hat, Gaussian, and ring initial distributions are shown in Fig. 4(a).No correlation is indicated by g(s) = 1, as is the case for the initial top-hat distribution [blue curve in Fig. 4(a)].As the cluster expands, periodic correlations emerge for all ensembles during expansion, as seen in Fig. 4(c) at t = 1000γ τ .In the long-time limit, the separation of vortices on opposite sides of the cluster also become correlated due to the formation of concentric rings as discussed in Sec.IV A. The correlation functions for the three ensembles exhibit the same period (i.e., there is a characteristic nearest-neighbor distance, next-nearest-neighbor distance, etc.) but differ in the strength of the correlation for the same evolution time.
To understand the magnitude of the correlations in the vortex cluster, we quantify the dynamics of the crystallization in Fig. 4(e).We plot the evolution of the geometric disorder parameter σ g = σ nn /μ nn , where σ nn and μ nn are the standard deviation and mean of the nearest-neighbor distances of vortices in the cluster, respectively [42].For a perfectly ordered Abrikosov lattice, σ g = 0.It is clear the strength of the correlation function depends upon how ordered the cluster is, e.g., the top-hat expansion (the most ordered state given by the lowest σ g ) is marked by the strongest correlations.
We find there is a period where all three ensembles evolve approximately as σ g ∼ t −α with fitting parameter α 4/9 [gray region in Fig. 4(e)].Although we do not have a theoretical explanation for the observed scaling, we note that it was found to be independent of γ for γ 0.01, suggesting it is a robust feature of the expansion.Additionally, there is a late time plateau where the disorder parameter ceases to decay.This persists for a wide range of N (albeit fluctuating slightly), suggesting lim t→∞ σ g = 0 is not an artifact due to a particular choice of vortex number.At this stage of the expansion, where all the vortex clusters are equally ordered, the vortex-vortex correlation function for each of the three cases exhibits approximately the same magnitude in the peaks [Fig.4(d)].The tendency for the cluster to reach a steady state where the disorder is nonzero suggests the outer concentric rings of vortices prevent the emerging Abrikosov lattice from spreading throughout the entire cluster.Despite clear signatures of an Abrikosov lattice emerging in patches of the clusters [Fig. 1, t = 1.3 × 10 5 γ τ ], they are broken by dislocations that arise in the dynamics and persist.The frustrated structure observed is distinct from the familiar vortex lattices seen in rotating superfluid systems [43][44][45] as this system is out of equilibrium.

A. Experimental procedure
In this section, we compare the results of the vortex fluid theory to data from experiments observing the expansion of vortex clusters in a quasi-2D BEC.A planar 87 Rb BEC of N c ∼ 2.2 × 10 6 atoms is trapped in a gravity-compensated optical potential.The horizontal confinement of the atoms in the x-y plane is provided by a repulsive blue-detuned optical dipole potential that results from the direct imaging of a digital micromirror device; here it is configured to produce a disk-shaped trap with radius R = 50 μm.This results in the approximately hard-walled confinement of the BEC with a near-uniform density.In the vertical direction, a red-detuned optical dipole potential provides harmonic trapping with frequency ω z = 2π × 108 Hz, leading to a vertical Thomas-Fermi radius of 6 μm.The healing length of the BEC is ξ ∼ 500 nm [46].
A schematic of the preparation of the initial vortex cluster is shown in Fig. 5(a).We form the BEC in the disk trap before transferring it to an annular trap by ramping-on an additional repulsive central barrier with a radius of R 0 = 15 μm over 200 ms.Simultaneously, an elliptical stirring barrier is introduced with a major and minor axis of 50 and 2 μm, respectively.This barrier crosses the annulus, which results in a split ring.Following previous techniques for stirring persistent currents in ring-trapped BECs [47,48], the stirring barrier is linearly accelerated at 980 μm s −2 for a time of 400 ms around the annulus.While still accelerating, the barrier height is then linearly decreased to zero by reducing both the barrier width and length over 100 ms, effectively removing the stirring barrier through to the central barrier.After a 400-ms period of equilibration in the annular trap, the central barrier is removed over 200 ms by linearly reducing its radius to zero.
This procedure results in a high-energy cluster of N ∼ 11 vortices within a radius of r c ∼ 8 μm near the trap center.During the removal of the stirrer, occasionally one or two stray vortices (of the same sign) are produced away from the main cluster [see left image in Fig. 5(b)].The system has a condensate fraction of N c /N ∼ 80%, and the small thermal cloud leads to weakly dissipative vortex dynamics.
We observe the expansion of the vortex cluster for 7 s, destructively sampling ∼40 independent images at intervals separated by 250 ms.High-resolution images of the BEC and vortex cores are obtained as in Ref. [19], where a short 3-ms time of flight expansion allows the vortex cores to expand and become visible using darkground Faraday imaging [49].The radial distribution of the condensate is essentially unchanged during this expansion.Two example images are shown in Fig. 5(b) at t = 0 s and t = 2 s.The vortex positions were identified using a Gaussian fitting algorithm [19,42], which can fail to detect vortices at early times when vortices are not well separated [e.g., t = 0 s in Fig. 5(b)].In Fig. 5(c), we show the two-dimensional histograms of observed vortex positions at 1-s intervals to demonstrate the expansion of the vortex cluster in the BEC.

B. Experimental results
To better visualize the vortex cluster expansion, we plot the average radial density of the cluster for three different time intervals in Figs.6(a)-6(c).In each case, we have fitted a top-hat distribution where the radius is found from measuring the cluster size in the experiment.We find that as the cluster expands, the density evolves towards a top-hat In Fig. 7(a) we plot the mean radius of the central cluster as a function of time and find excellent agreement with the r ∼ √ t prediction of the vortex fluid theory.To emphasize the strength of the agreement, we fit a linear curve (dashed line) to the late-time data.We fit a dissipation constant of γ = 3.2 × 10 −3 ; this is consistent with previous experiments of atomic BECs [19,32,50] and sets the experimental time in units of γ τ to be ∼4.1 × 10 4 γ τ (for comparison with simulation).
Figure 7(b) shows the energy decay of the vortex cluster.Given the experiment is in a circular bounded domain of radius R, the Hamiltonian for the vortex cluster is now where r i = |r i | is the magnitude of the radial position of the vortex [51].The first term of Eq. ( 12) describes the energy due to vortex-vortex interactions and is simply Eq. ( 1).The second and third terms describe the interactions with fictitious image vortices that enforce the boundary condition u • n| r=R = 0, which ensures the superfluid flow normal to the trap boundary is zero.As the vortex cluster is located towards the center of the trap (even at late times), the additional energy due to the confinement is negligible (i.e., the second and third term have little contribution to the energy).
As mentioned in Sec.V A, we encounter difficulty in imaging vortices when the cluster radius is small, which manifests as an artificial growth in N, as seen in Fig. 7(c).Due to the strong dependence of energy upon N (H ∼ N 2 ), we plot H/N 2 to factor out the vortex number dependence upon the energy.We observe a discrepancy in energy between experiment and theory due to only having N ∼ 11 vortices in the experiment as the vortex fluid theory assumes N 1.However, we find that applying a simple scaling factor N 2 /(N 2 − 2N ) fits the theory well [52].Performing point vortex simulations of varying N in a bounded disk, we find this scaling to yield good agreement with the vortex fluid theory for all N 10 (see Appendix C for further details).
Finally, in Fig. 7(d) we measure the disorder of the vortex cluster σ g experimentally and compare it to simulations of the bounded point vortex model using the ensemble of initial vortex positions taken from the experiment at t = 1.0 s.We simulate the entire ensemble of vortices observed (i.e., stray vortices included) but only plot the ordering of the central cluster.We find strong agreement between the point vortex model and the experiment, which both show the cluster does not become more ordered over time.If we only simulate the dynamics of the vortices in the central cluster (i.e., do not include the stray vortices in the initial conditions), we find the simulations predict a decrease in σ g with time, as was found in Sec.IV B. This suggests that the stray vortices significantly suppress the ordering of the vortex cluster.Figure 7(c) shows that at the end of the experiment, there tends to be only one stray vortex in the system.Remarkably, this suggests a single stray vortex is sufficient to destroy the ordering of the cluster.Without stray vortices, our simulations show that dissipation can provide sufficient ordering to study the crystallization dynamics in future experiments with larger vortex numbers.

VI. CONCLUSIONS
We have analytically shown within dissipative vortex fluid theory that any dense vortex cluster in a finite-temperature quantum fluid evolves to form a Rankine vortex, confirming a new universality class in dissipative superfluids.Numerical simulations of the microscopic dynamics of the vortices confirms this universal expanding regime, as well as revealing phenomena beyond the predictions of vortex fluid theory.We presented the emergence of a frustrated lattice structure within the vortex cluster, which is approached through a power-law decay in the geometric disorder and is marked by strong vortex-vortex correlations.Finally, we have presented experimental observations of vortex cluster expansions in a quasi-2D BEC and found they are in good agreement with the vortex fluid theory and our numerical simulations, despite being for a cluster of only 11 vortices.Our experimental results validate the vortex fluid theory, paving the way to better understand dissipative mechanisms in quantum fluids.
Whereas the Rankine vortex is forbidden in viscous classical fluids [25], our results suggest it may be highly relevant in finite-temperature dissipative superfluid systems.Our findings also suggest that recently predicted phenomena associated with the superfluid Rankine vortex, such as quantized edge solitons [23], may be within reach experimentally provided larger vortex numbers can be achieved.Superfluids with a larger ratio between the system size and healing length, such as strongly interacting Fermi gases [11] or thin-film superfluid helium [21], may be promising alternative platforms to test these predictions.Beyond the vortex fluid theory, the emerging fractured lattice structure in the vortex clusters could have Kibble-Zurek-type behavior, as well as demonstrating qualities reminiscent of the hexatic phase observed in systems such as liquid-crystal films [53].Hence the complete set of equations for dissipative chiral vortex fluids reads

Dynamics of density perturbations
We demonstrate that the Rankine vortex solution [Eq.(7)] is an asymptotic solution of Eq. ( 6) of the main text (anomalous hydrodynamical equation for the vortex density) in the long-time limit, independent of initial vortex distributions.Consider a local density perturbation ρ = ρ a + δρ to the universal expansion solution, where ρ a is given by Eq. ( 7), |δρ| |ρ a |, and δρ(|r| → ∞) = 0.For simplicity, we assume that δρ has cylindrical symmetry, namely, δρ(r) = δρ(r).Keeping the leading-order terms in δρ, we obtain where v a = r ρ a /2.Solving Eq. (A14) numerically for a variety of initial conditions, we find that for all cases the perturbation decays over time and reduces to zero across the fluid δρ(t → ∞) → 0. Figure 8 shows four examples where the initial perturbation takes different forms, as specified in the figure caption.Here the boundary conditions are δρ regular at r = 0 and δρ(r → ∞) → 0.
This result suggests that the top-hat distribution is an attractor of the dissipative vortex fluid dynamics and hence further supports our conclusion that any initial distribution of vortex density will eventually tend towards a Rankine vortex.

APPENDIX B: COHERENT EVOLUTION OF VORTEX CLUSTERS
To emphasize the dramatic difference between the dynamics of large clusters of vortices in the presence or absence of dissipation, we plot the evolution of the five different initial distributions when the dissipation is zero [i.e., solve Eq. ( 4) for γ = 0].In Fig. 9, it is clear that the Rankine vortex does not emerge in these situations and instead the final vortex distributions are strongly dependent upon the initial conditions.

APPENDIX C: POINT VORTEX SIMULATIONS FOR FINITE-SIZED SYSTEMS
The vortex cluster expansion experiments we describe in the main text are performed with ∼ 11 vortices in an atomic gas BEC.The superfluid is confined in a flat-bottomed, disk trap of radius R using the experimental techniques described in the main text and Ref. [54].Similar systems have been realized for superfluid helium [21,55].
The vortex fluid hydrodynamic theory, however, is for an infinite system with no boundaries in which the vortex density is coarse grained and treated as a continuous quantity.Thus an important question is, How applicable is the vortex fluid theory to the experimental observations?In this section we compare simulations of point vortex dynamics in the disk trap with the predictions of the time dependence of the cluster radius and energy from the vortex fluid theory.
The Hamiltonian for N vortices at positions r i in a circular domain of radius R is given in Eq. ( 12).Hamilton's equations of motion for this system are given by Eq. ( 2), which yields the equation of motion for vortices in a circular disk, where x i j = x i − x j , y i j = y i − y j , and r 2 i j = x 2 i j + y 2 i j .The first term of Eq. (C1) arises from the flow field of the other vortices, while the second term arises from the image vortices.The barred terms xij = x i − x j , r2 i j = x2 i j + ȳ2 i j , etc. correspond to the positions of the image vortices, with circulation ¯ j = − j , located outside the disk at the inverse point ri = R 2 r i /|r i | 2 .Using Eq. (C1) for the bounded vortex velocities in Eq. ( 4) of the main text, we solve for the dynamics of bounded vortices in the dissipative regime.
When a vortex approaches the boundary and pairs with its image, the energy of the system reduces while the velocity of the vortex increases dramatically.This can be seen in Fig. 10(a), where we plot the energy (left axis) and radius (right axis) of two positive vortices in the bounded point vortex model.Unsurprisingly, the hydrodynamic approximation is a poor description of the radius of the vortex pair (plotted as the solid curve), as the core assumption is that vortices are densely packed so the vorticity can be coarse grained.Comparing this result with Fig. 10(b), where we plot the energy and radius of an expanding cluster of N = 5 vortices (the initial condition is drawn from a uniform distribution and sits within r = 0.1R), we see that the boundary has much less of an effect upon the dynamics.As there are more vortices within the center of the disk, the energy is dominated by vortex-vortex interactions, unlike the N = 2 case where the images strongly influence the dynamics.As image vortices are positioned at the inverse point, their effect upon the system reduces significantly for vortices far from the boundary.As a result, we begin to see signatures of the anomalous hydrodynamics [i.e., r ∝ √ t and H ∝ − ln(t )] despite there being only N = 5 vortices.
Upon further increasing the number of vortices to N = 10 in Fig. 10(c), the difference between the bound simulations and the free-space anomalous hydrodynamics is further reduced.This particular result is close to the experimental system, and we see good agreement between the vortex fluid theory and bound point vortex simulation.For N = 20 [Fig.10(c)], the expansion of the cluster in the point vortex model is close to the vortex fluid prediction.For even larger numbers the difference continues to decrease, and so in conclusion, it seems that for a sufficiently large vortex number (N > 20), an expanding cluster in a bounded circular domain can be closely approximated as a free-space expansion until the cluster radius nears the boundary.
There is a distinct difference, however, between the energy in the point vortex simulation and that of the anomalous hydrodynamic solution.Given the anomalous hydrodynamics assumes N 1, we find there are finite N effects for a small number of vortices.For N > 10, we find that scaling the point vortex simulation energy (and indeed the exper-iment, seen in the main text) by N 2 − 2N gives excellent agreement with the scaled hydrodynamic energy H/N 2 .It can be seen that the difference vanishes in the N → ∞ limit.

FIG. 1 .
FIG. 1. (a)Exemplar vortex distributions throughout the expansion in free space with γ = 0.01.We show the (a) top-hat, (b) Gaussian, and (c) ring initial distributions as well as two nonaxisymmetric initial conditions in (d) and (e).The time in the simulation is indicated on the left.The radii of the dotted lines in the first three rows, emphasizing the expansion, have a constant radius of r = 0.2.The r = 0.2 dashed line is omitted in the final row as it cannot be seen.In the bottom row we have highlighted lattice dislocations (defined when vortices have either five or seven nearest neighbors) with black points.

FIG. 2 .
FIG. 2. Comparison between vortex fluid theory (VFT) and numerical results for (a) average radius of vortex cluster and (b) energy of the vortex cluster.In the bottom panel we plot the error between the hydrodynamic solution and the numerical simulation.

FIG. 3 .
FIG. 3. Normalized ensemble radial vortex density for the tophat, Gaussian, and ring initial conditions at times (a) t = 0, (b) t = 10γ τ , (c) t = 100γ τ , and (d) t = 1000γ τ averaged over n = 100 different realizations for each distribution.The black curve in (d) is the coarse-grained density for the top-hat initial distribution.(e) Normalized average density as a function of radius and time for the top-hat initial distribution.

FIG. 5 .
FIG. 5. (a) Schematic of the procedure used to create the initial vortex cluster.An elliptical stirrer is swept around the annulus as its length is decreased over time, creating same-signed vortices (denoted by a +) which become pinned to the central barrier.However, some vortices can become stray and remain in the bulk of the fluid.The radius of the central barrier is then shrunk to zero, and the pinned circulation (denoted by large + at the center) splits into singly charged vortices which are free to expand.(b) Experimental image after 0-and 2-s hold time and 3-ms time-of-flight expansion showing resolved vortex cores as dark holes in the atom density.White circles indicate a vortex detected by the Gaussian fitting algorithm.(c) Two-dimensional histograms of the vortex positions over time, averaged over ∼40 runs.The red dashed circle indicates the cutoff radius for cluster analysis (see text).

FIG. 6 .
FIG. 6. (a-c) Radial density of vortex clusters averaged over 1 s (period given in top right) along with top-hat distribution fit.Insets: corresponding 2D histograms of vortex positions.The red dashed line represents the cluster cutoff radius.(d) Error between histograms and top-hat fit.

FIG. 7 .
FIG. 7. (a) Comparison of vortex fluid theory and experimental measurements of the average vortex radius r .The dashed line is a linear fit to the late-time data.The gray region is prior to the removal of the pinning potential.(b) Energy of experiment and vortex fluid theory.Experimental energy is also shown after being scaled (hollow points) by a constant factor (see text).(c) Average measured number of vortices in the system (i.e., cluster and stray vortices) and the number of vortices measured in the cluster.Each exhibits an artificial growth due to imaging techniques.(d) Observed geometric disorder σ g of the vortex cluster compared with point vortex simulations with and without stray vortices used to calculate the dynamics.Legend in (a) applies to all four subfigures.Uncertainty in the experiment is given by the standard error.

FIG. 9 .
FIG. 9. Evolution of vortex clusters without dissipation (i.e., γ = 0) for initial distributions corresponding to (a) top-hat, (b) Gaussian, (c) ring, (d) one-dimensional line, and (e) multiple random clusters.The upper and middle row are exemplar snapshots of the cluster at the beginning and end of the simulation (as labeled).The bottom row is the average radial density (averaged over n = 20 realizations) of distributions before and after evolution.

FIG. 10 .
FIG. 10.Comparison of radius and energy of point vortex cluster expansion with the predictions of the vortex fluid hydrodynamics.(a) A system of N = 2 vortices initially equidistant from origin separated by θ = π .The black curves correspond to the left y axis and are the energy of the vortex cluster, while the orange curves correspond to the right y axis and indicate the average radius of the vortex cluster.Dashed lines are from simulations of the point vortex model within a disk, and the solid lines are the solution to the vortex fluid hydrodynamics.(b) N = 5 vortices.(c) N = 10 vortices.(d) N = 20 vortices.