Vortex Thermometry for Turbulent Two-Dimensional Fluids

We introduce a new method of statistical analysis to characterise the dynamics of turbulent fluids in two dimensions. We establish that, in equilibrium, the vortex distributions can be uniquely connected to the temperature of the vortex gas, and apply this vortex thermometry to characterise simulations of decaying superfluid turbulence. We confirm the hypothesis of vortex evaporative heating leading to Onsager vortices proposed in Phys. Rev. Lett. 113, 165302 (2014), and find previously unidentified vortex power-law distributions that emerge from the dynamics.

We introduce a new method of statistical analysis to characterise the dynamics of turbulent fluids in two dimensions. We establish that, in equilibrium, the vortex distributions can be uniquely connected to the temperature of the vortex gas, and apply this vortex thermometry to characterise simulations of decaying superfluid turbulence. We confirm the hypothesis of vortex evaporative heating leading to Onsager vortices proposed in Phys. Rev. Lett. 113, 165302 (2014), and find previously unidentified vortex power-law distributions that emerge from the dynamics.
Turbulence arises in chaotic dynamical systems across all scales, from mammalian cardiovascular systems, to climate, and even to the formation of stars and galaxies [1]. The unpredictability inherent to turbulent systems is further confounded by physical properties such as boundaries and spatial dimensionality, and due to its complexity, there is currently no unified theoretical framework to explain turbulence. As such, there is a need to develop new methods to characterise the evolution of turbulent states in order to provide further insights into this important problem.
Onsager developed a model of statistical hydrodynamics to describe turbulence in two-dimensional (2D) flows [2]. In this representation the fluid is modelled by an N-particle gas of interacting point-like vortices which can be characterised by an equilibrium temperature. As the bounded system of vortices has a finite configuration space, the entropy S of the system has a maximum, and hence there is a range of energy E where the absolute Boltzmann temperature T = (∂S /∂E) −1 becomes negative [2][3][4]. These states correspond to large-scale rotational flows known as Onsager vortices [2].
Experimentally, BECs provide unprecedented opportunities to investigate 2D superfluid turbulence due to the high degree of controllability available in these systems. It is now possible to create and image complex vortex configurations such as dipoles [30][31][32][33], few-vortex clusters [34] and quantum von Kármán vortex streets [35]. Many experiments have also been devoted to the study of quantum turbulence in both two- [11,19,20,23] and three-dimensional [22,36,37] geometries. However, the formation of Onsager vortex structures in statistical equilibrium has not yet been reported. Recent theoretical works have suggested that one significant obstacle is the harmonic trapping commonly used in experiments, as vortex clusters appear to be suppressed in this geometry [14,20,38]. In addition, the detection of vortex circulation signs is experimentally difficult, and only recently have techniques been proposed [39] and implemented [23] to achieve this. Analysis of turbulent dynamics is made even more challenging by current condensate imaging methods, which only allow a small number of frames to be captured for a single experimental realisation [40]. As such, it is desirable to be able to characterise the state of a turbulent fluid using a robust method of statistical analysis that links the instantaneous microscopic configuration of the system to its macroscopic behaviour. Onsager's thermodynamical description of turbulence is one such method, and hence we propose to use its central observable-the vortex temperature-for this purpose. In contrast to velocimetry-based observables that require the measurement of the velocities of the atoms, the thermometry presented here only requires the measurement of the positions and circulation signs of the quantised vortices.
We first outline our method for measuring the temperature of the vortex gas, before examining a specific case of decaying superfluid turbulence using mean-field Gross-Pitaevskii simulations. In the dynamics, we observe that the vortex gas undergoes rapid equilibration before settling into a quasiequilibrium state where it continues to heat adiabatically via vortex evaporation [13]. We have discovered that in this evolution, the numbers of clusters, dipoles and free vortices follow robust power-laws with respect to the total vortex number. The existence of this quasi-equilibrium allows quantitative thermometry of the turbulent fluid.
To calibrate the vortex thermometer, we use Monte Carlo (MC) simulations to map out the equilibrium vortex configurations as a function of the inverse temperature β = 1/k B T , where k B is Boltzmann's constant. We do this for a system of lar boundary of radius R [13,43], and set a hard vortex core of radius 0.003 R to prevent energy divergences. As we vary the temperature across both positive and negative regimes, we quantify the effect on the vortex configuration using a vortex classification algorithm [10,44]. The algorithm uniquely divides the vortex gas into three separate components: clusters of ≥ 2 like-sign vortices, closely bound vortex-antivortex dipoles, and relatively isolated free vortices (for further details, see Ref. [44]). We then calculate the number of clusters N c , dipoles N d and free vortices N f as functions of temperature, and the resulting fractional population curves are presented in Fig. 1.
At low positive absolute temperatures (left hand side of Fig. 1), the vortex gas is at its 'coldest', as both the energy and entropy are minimised. In this regime, bound vortexantivortex dipole pairs dominate the configuration, as shown in the schematic inset of Fig. 1. Above the Berezinskii-Kosterlitz-Thouless (BKT) critical temperature β BKT [45][46][47][48][49], the vortex dipoles dissociate, causing an increase in both the energy and entropy. At β = 0, the vortex configuration becomes a disordered arrangement of vortices and antivortices, thereby maximising the entropy. In the negative temperature region, low-entropy clusters of like-sign vortices tend to form (see schematic inset), and because of their high energy, these negative temperature states are 'hotter' than those at positive temperature. Above the critical temperature β EBC , the vortices form an Einstein-Bose condensate (EBC), a state where the Onsager vortex clusters condense, as indicated by the satura-tion of the cluster population in Fig. 1 [13,44,50]. For a neutral vortex gas the two aforementioned critical temperatures are defined as β BKT = 2/E • and β EBC = −4/N v E • [6], respectively, where the energy E • = ρ s κ 2 /4π is defined in terms of the superfluid density ρ s and the quantum of circulation κ = h/m, with m being the mass of the condensed atoms. Figure 1 demonstrates that the dipole and cluster populations are monotonic functions of β-this is the key observation enabling thermometry of the vortex gas. Given an arbitrary vortex configuration in thermal equilibrium, we may determine its temperature by calculating the populations of clusters and/or dipoles and comparing the obtained values to the curves in Fig. 1. Strictly, the p j (β) curves in the negative temperature region of Fig. 1 are dependent on the vortex number. However, we repeated our MC simulations for N v = 100 and 200 vortices and found that, for the vortex numbers relevant to the dynamical simulations presented here, there is no qualitative change to the thermometry curves, and the quantitative change is not significant (see Supplemental Material [41]). The cluster and dipole fractions are not the only observables that vary monotonically with vortex temperature in our MC simulations. For example, both the energy and dipole moment of the vortex gas also fulfill this requirement [13], and could therefore, in principle, be used for thermometry. However, of all variables considered, we have found that the cluster and dipole fractions provide the most robust thermometers.
Also shown in Fig. 1 is the Einstein-Bose condensate fraction, which quantifies the density of vortices in the largest cluster (for details, see Ref. [44]). For β > β EBC , the condensate fraction is zero, but when β < β EBC it rises sharply. In this extreme temperature region, the other thermometers saturate and the condensate fraction becomes the relevant observable for vortex thermometry.
As an application of our vortex thermometer, we use it to characterise decaying turbulence in a disk-shaped BEC as previously studied in Refs. [13,14]. We simulate the two-dimensional time-dependent Gross-Pitaevskii equation where ψ ≡ ψ(r, t) is the classical field of the Bose gas and g 2D is the two-dimensional interaction parameter resulting from the s-wave atomic collisions. To obtain the uniform circular geometry, we use a two-dimensional power-law trapping potential of the form V tr = µ(r/R) 50 , where r = x 2 + y 2 is the radial distance from the axis of the trap, µ is the chemical potential, and R ≈ 171 ξ is the radius of the trap, measured in units of the healing length ξ = 2 /2mµ [14]. The interaction parameter is set to g 2D = 4.6 × 10 4 2 /m. We solve the GPE using a fourth-order split-step Fourier method on a 1024 2 numerical grid with a spacing of approximately ξ/2. Turbulence is generated by imprinting vortices into the phase of ψ and evolving Eq. (1) for a short amount of imaginary time to establish the vortex core structures. We detect vortices and their circulation signs within r < 0.98 R by locating singularities in the phase of the field.
The initial vortex configurations used in our GPE simulations are produced by randomly drawing N v = 800 vortex locations from a uniform distribution, with equal numbers of each circulation sign. The resulting state is well approximated to have β ≈ 0, although the short imaginary time propagation step causes a small amount of cooling towards positive temperatures. As the turbulence decays, the vortices annihilate and the vortex gas evaporatively heats, resulting in the emergence of two large Onsager vortices at late times [13,14]. Three sample frames from a single simulation are shown in Fig. 2, where panels (a)-(c) show the density |ψ| 2 of the fluid, and panels (d)-(f) show the corresponding vortex configuration after the vortex detection and classification steps. A Helmholtz decomposition [8] has been used to extract the divergence-free component of the condensate velocity field, and the resulting streamlines are also shown in the lower panels. The Onsager vortex clusters are clearly visible in panel (f).
The number of clusters, dipoles and free vortices are shown in Fig. 3 as functions of both time t (inset) and the total number of vortices N v (t). The time-dependent populations (inset) do not follow any simple function. However, the populations as functions of the total number of vortices (main frame) show clear power-law scaling behaviour. The corresponding powerlaws are: of vortices per cluster N vc ≡ N c /N cl , where N cl is the total number of clusters of any size at a given time. To study the effects of system size on these power-laws, we have also considered two smaller disk-shaped systems of radii R ≈ 49 ξ and R ≈ 85 ξ respectively, each with N v = 100 vortices initially imprinted. We find that the scaling behaviour is unchanged in these smaller systems, suggesting that the evolution of the vortex gas is underpinned by a universal microscopic process.
In this system, the primary cause of vortex number decay is the annihilation of vortex-antivortex dipoles. Despite this, the populations of dipoles and free vortices follow approximately the same power-law, demonstrating an interconversion between the vortex populations. However, a distinct powerlaw emerges for the vortex clusters. This behaviour points toward a two-fluid model, where the dipoles and free vortices behave as a weakly interacting thermal cloud, while the clusters act as a quasi-condensate whose relative weight grows over time as a result of vortex evaporative heating. Extrapolating the data toward N v → 0 leads to the inevitable decay of all dipoles and free vortices, with only Onsager vortex clusters remaining. At this point, the rate of pair annihilation becomes insignificant in the dynamics due to the very low probability of vortex-antivortex collisions.
In Fig. 3, the N d and N f curves are well described by the N 6/5 v scaling throughout the dynamical evolution. The N c curve, on the other hand, only begins to follow the N 4/5 v powerlaw once the total vortex number has decayed to N v 200, suggesting that the statistical behaviour of the vortex gas changes at this point in the dynamics. In accordance with the existence of power-law scaling, we interpret this change to be the realisation of a state of quasi-equilibrium for the The temperature is measured independently using the populations of both clusters and dipoles. In the inset, the temperature readings from each thermometer is shown as a function of the total vortex number N v (t). As in Fig. 1, the positive and negative temperature regions have been scaled by their respective critical temperatures, and a dashed horizontal line denotes β = 0. The vertical axis of the inset is the same as for the main frame. decaying turbulence. Under this quasi-equilibrium condition, the vortex evaporative heating process becomes adiabatic in the sense that the vortex gas is able to rearrange into a higher entropy configuration between the vortex annihilation events. For N v 200, the vortex number decays too rapidly for this to be possible. This quasi-equilibrium condition is not a true equilibrium of the system, since vortex-antivortex annihilations and vortex-sound interactions are continuously driving energy from the vortices into the sound field. Presumably, the true equilibrium of the condensate will only be realised when all vortices have decayed and the total entropy of the system is maximised. In the Supplemental Material [41], we present vortex number decay data for a range of other initial vortex configurations, observing in all cases evidence for the same power-law and quasi-equilibration behaviour.
We now have an algorithm to assign a vortex temperature to the dynamical GPE simulations. We determine the fractional populations of vortex dipoles and clusters as a function of time, and use each of these to read off a temperature from the curves in Fig. 1. The two resulting measurements of β(t) are presented in the main frame of Fig. 4. Both measurements show that the temperature of the vortex gas is spontaneously increasing as Onsager vortex clusters form, thereby confirming the evaporative heating scenario posited in Ref. [13]. At late times, a small discrepancy between the two temperature readings emerges, which we attribute to the compressibility of the fluid not accounted for in the MC model. The same temperature measurements are plotted as a function of the total vortex number in the inset. Based on our quasi-equilibrium interpretation discussed above, we note that the temperature reading is strictly only valid for N v 200 (t 2000 /µ), since outside of this range the vortices are out of equilibrium and their temperature is not well defined. To obtain these curves, we have applied the thermometer calibrated with N v = 50 vortices, despite the fact that N v varies between ≈ 30 and ≈ 200 throughout the equilibrium dynamical evolution. In the Supplemental Material [41], we show that using a thermometer calibrated with a different number of vortices does not affect the qualitative shape of the β(t) curve.
We have developed a methodology that allows the temperature of point vortices in two-dimensional fluids to be determined using only the information about the vortex positions and their signs of circulation. We have applied the vortex gas thermometers to freely decaying two-dimensional quantum turbulent systems and quantitatively shown the transition to negative temperatures and the emergence of Onsager vortices due to the evaporative heating of the vortex gas [13,14]. Our vortex thermometers may also be useful for characterisation of turbulent classical fluids, as the continuous vorticity distributions can be approximated accurately by a discretised set of point vortices before performing the vortex classification and thermometry. This methodology may therefore open new pathways to quantitative studies of two-dimensional turbulence.
We  To assess the sensitivity of the observed power-laws (Fig. 3  of the main text) to the choice of initial vortex configuration, we have run Gross-Pitaevskii simulations with a diverse range of initial conditions. In addition to the randomly sampled initial condition (case I) described in the main text, we have considered four other types of initial state. The cases II and III are configurations with lower incompressible kinetic energy E i k , created by imprinting the vortices randomly throughout the condensate as dipole pairs with sizes 8 ξ and 12 ξ, respectively (before the imaginary time propagation step). For case IV, motivated by experiments (e.g. [1][2][3]), the vortex creation is simulated dynamically by stirring an initially unperturbed condensate with a repulsive Gaussian potential of waist 30 ξ and amplitude 5 µ. The stirring potential is moved back and forth with centroid position x • (t) = 100ξ cos(2πµt/1050 ) for four periods, and then ramped down to zero over a fifth period. Finally, a large incompressible kinetic energy in case V is initiated by imprinting a periodic square array of vortex clusters with alternating circulation sign, each with a radius of ≈ 43 ξ and containing up to 25 randomly placed vortices. Examples of initial conditions for cases II-V are shown in Fig. 5

FIG. 5. Examples of initial vortex configurations for our turbulent
Gross-Pitaevskii simulations. Panels (a)-(d) correspond to cases II-V, respectively, and show the vortices in positive (negative) clusters as blue (green) squares, dipoles as red triangles, and free vortices as yellow circles. Note that each vortex dipole contains one vortex and one antivortex. The streamlines in each frame are obtained by calculating the incompressible component of the velocity field of the classical field describing the Bose gas. Panel (c) shows the dynamically stirred configuration immediately after the stirrer has been switched off. See also Fig. 2(d) in the main text, which shows the initial condition for case I. Panel (e) shows the incompressible kinetic energy E i k per vortex for each of the five cases, averaged over 80 (10) initial conditions generated for case I (cases II-V). Error bars denote one standard deviation. The shading is indicative of how 'cold' (blue) or 'hot' (red) a given initial state is.
where the vortices have been classified into clusters, dipoles and free vortices as described in the main text. The corresponding mean incompressible kinetic energy for each initial condition is shown in Fig. 5(e). This energy is defined E i k = (m/2) |ψ| 2 |v i | 2 d 2 r, where the v i (r) is the divergencefree component of the total velocity field.
The resulting number decay curves for each vortex type are shown in Fig. 6(a)-(c). The dipole and free vortex decay curves [panels (b) and (c), respectively] remain relatively unchanged across different initial configurations. By contrast, the clusters [panel (a)] show clear variation across the set of initial conditions, suggesting that initially the system is in fact behaving very differently under each constraint. However, despite initial differences (at large N v ), all cluster decay curves eventually exhibit behaviour consistent with the power-laws obtained in Fig. 3 of the main text, demonstrating a loss of memory of the initial vortex configuration. This provides further evidence that these power-laws correspond to a state of quasi-equilibrium in which the vortex gas should have a welldefined temperature, as argued in the main text. In Fig. 6(a), the approximate value of N v at which the N c curve begins to follow the N 4/5 v power-law is also highlighted, which we interpret as being the point at which the vortex gas reaches quasiequilibrium.
Applying our thermometer ( Fig. 1 from main text) to all five initial conditions, we obtain the temperature readings for each, which are presented in Fig. 7. Here, the temperature is calculated from the mean of the cluster and dipole thermometer measurements, which are themselves ensemble averaged over 80 (10) simulations for case I (cases II-V). These curves show a clear dependence on initial condition, with the low energy configurations (cases II and III) being consistently colder than those with high energy (cases IV and V). The random initial configuration (case I) lies between the two extremes. The approximate value of N v at which the vortex gas appears to reach  equilibrium in each case [see Fig. 6(a)] is also shown. Even before this point (i.e. for larger N v ), the vortex thermometer provides a plausible temperature reading, but the measurement is not reliable if the vortex gas is out of equilibrium. In cases IV and V, the equilibration point corresponds to a turning point in the temperature curve, providing further evidence for our interpretation of the vortex gas equilibrium condition.

VORTEX NUMBER DEPENDENCE OF THERMOMETRY
In Fig. 4 of the main text, the temperature measurement of the decaying turbulence for case I was obtained using a single thermometer calibrated with N v = 50 vortices. Strictly, this thermometer is only quantitatively valid for the short time in the dynamical evolution when N v ≈ 50, and additional thermometers should be calibrated to obtain more accurate measurements for other vortex numbers. Here we demonstrate that changing N v in the Monte Carlo simulations has only a small effect on the thermometry curves, and consequentially on the dynamical β(t) measurements.
We have repeated our Monte Carlo simulations with N v = 100 and N v = 200 (as argued in the main text, the vortex gas appears to be out of equilibrium for N v 200 for case I, and hence any vortex numbers beyond 200 are not relevant for thermometry here). The obtained p c (β) and p d (β) curves are shown in Fig. 8 for all three values of N v . The thermometry curves show very little variation as N v is changed, especially in the positive temperature region. Most importantly, the curves are always monotonic, regardless of N v , and hence can always be used for thermometry.
Using the three cluster thermometers in Fig. 8, we have remeasured the dynamical temperature β(t) from our Gross-Pitaevskii simulations, and the resulting curves are presented in Fig. 9. Evidently, our measurement using the N v = 50 FIG. 9. Inverse temperature of the vortex gas as a function of time, averaged over a set of 80 dynamical GPE simulations for case I, and measured using three different thermometers, calibrated using N v = {50, 100, 200}, respectively. The temperature is measured using the populations of clusters only. As in Fig. 8, the positive and negative temperature regions have been scaled by their respective critical temperatures, and a dashed horizontal line denotes β = 0. The mean number of vortices remaining in the system is indicated at particular times by vertical dotted lines. See also Fig. 4 of the main text.
thermometer slightly underestimates the temperature at early times (50 N v 200), and slightly overestimates it at the latest times (N v 50). Despite this minor quantitative correction, the qualitative behaviour of vortex heating-our main conclusion from the data-is unchanged regardless of which thermometer is used. Note that we have measured these temperatures using the spline fits to the data in Fig. 8, as described in the main text. The fluctuations that are visible in the β(t) curves arise from the variations in p c (t), which then appear in each temperature measurement.