## Clustering vortices

Many-body systems generally become more disordered as more energy is pumped into them. A curious exception to this rule was predicted in the context of turbulent flow by the physical chemist Lars Onsager. He suggested that the entropy of certain two-dimensional (2D) systems can decrease with increasing energy, corresponding to an effective negative temperature. Using 2D Bose-Einstein condensates of atoms, Gauthier *et al.* and Johnstone *et al.* put Onsager's theory to the test. They provided energy to the system by perturbing the condensate, creating vortices and antivortices. With increasing energy, the system became more ordered as clusters containing either only vortices or only antivortices emerged.

## Abstract

Nonequilibrium interacting systems can evolve to exhibit large-scale structure and order. In two-dimensional turbulent flow, the seemingly random swirling motion of a fluid can evolve toward persistent large-scale vortices. To explain such behavior, Lars Onsager proposed a statistical hydrodynamic model based on quantized vortices. Here, we report on the experimental confirmation of Onsager’s model. We dragged a grid barrier through an oblate superfluid Bose–Einstein condensate to generate nonequilibrium distributions of vortices. We observed signatures of an inverse energy cascade driven by the evaporative heating of vortices, leading to steady-state configurations characterized by negative absolute temperatures. Our results open a pathway for quantitative studies of emergent structures in interacting quantum systems driven far from equilibrium.

Statistical mechanics provides a description of the thermodynamic equilibrium behavior of a system of many interacting particles that individually exhibit seemingly random thermal motion. The system can be taken out of equilibrium by the input of energy, but in isolation is expected to relax back to statistical equilibrium, a process accompanied by an increase in the entropy or randomness of the system. In some situations, however, the system can exhibit a steady-state collective motion of the constituent particles, corresponding to statistical equilibrium in a highly excited state. Examples of such behavior can be found in systems ranging from the flocking of birds (*1*) to labor productivity (*2*) and turbulence in two-dimensional (2D) flow. Turbulence is considered a quintessential example of a system driven far from equilibrium. In hydrodynamic turbulence, the kinetic energy of the fluid is transported without loss across many length scales. However, despite the complexity of the phenomenon, there are some statistical theories that describe the steady-state behavior of turbulent systems, such as Kolmogorov’s power-law scaling of energy flow (*3*). In 3D turbulence, a process known as a Richardson cascade (*4*) results in energy transport to ever smaller length scales, causing vortices to break up over time and the system to appear chaotic. Remarkably, restricting the fluid dynamics to 2D results in an inverse-cascade process: energy flows toward the largest length scales available, resulting in system-scale, persistent vortex flows (*5*). This behavior has been observed in systems ranging in scale from soap films (*6*) to Jupiter’s atmosphere (*7*, *8*).

Onsager proposed an explanation for the appearance of large-scale vortex flow in 2D turbulence in terms of equilibrium statistical mechanics of a model of quantized point vortices (*9*), noting its applicability to superfluids. He assigned an absolute temperature *T* to the point vortex system as *S* is the Boltzmann entropy of the vortex configuration corresponding to the logarithm of the number of possible vortex configurations that result in a flow field with energy *E*, β is the inverse temperature, and *k*_{B} is the Boltzmann constant. The lowest energies are produced by weak flows, which correspond to low-entropy (ordered) states, with β → ∞, in which vortices of opposite sign pair up, whereas at higher energies, β → 0, the vortices arrange into increasingly uncorrelated (high-entropy) configurations. For a confined system, however, the configuration space of vortices is bounded, and *S* must take a maximum value at a finite energy *E*_{0}. Onsager noted that for *E* > *E*_{0}, the absolute temperature will be negative because higher energies can only be obtained at lower entropies—the vortices must become more ordered, with like-sign vortices clustering to produce stronger flows. If energy is continuously injected into a confined 2D flow, then eventually large-scale compound “Onsager vortices” will be the only remaining features, with β → –∞.

Negative absolute temperature states were subsequently used to describe manufactured distributions of nuclear spin systems (*10*) and, more recently, motional degrees of freedom of cold atoms in optical lattices (*11*). In these experiments, changing spin states and the sign of interatom interactions, respectively, resulted in a forced population inversion with the resulting state decaying to lower energy while increasing its entropy. The vortex system described by Onsager is markedly different from the single-particle nuclear spin and optical lattice experiments: energy is injected into the system in a continuous manner and the interactions of the constituent particles (vortices) result in a negative temperature configuration.

In a superfluid, 3D quantum turbulence (3DQT), which manifests as tangles of quantized vortex lines, has been shown to exhibit a direct energy cascade similar to its classical counterpart. The statistical dynamics of 3DQT have been studied over the past three decades both numerically and experimentally in superfluid helium (*12*). More recently, 3DQT has been observed in atomic Bose–Einstein condensates (BECs) (*13*–*15*), where direct imaging is possible owing to the comparatively large vortex cores (micrometer size versus angstrom size in superfluid helium). In 2DQT, an inverse energy cascade is predicted to result from the preferential transport of the energy injected into the superfluid through the creation of vortex–antivortex pairs (*16*, *17*). Because the vortices are quantized, the spatial clustering of like-sign vortices forms the equivalent of classical large-scale flows, as in Onsager’s model. To date, however, the major challenge hindering the understanding of 2DQT has been to devise a method to experimentally measure the velocity field of a superfluid.

Atomic BECs provide an ideal system in which 2DQT can be realized, as they can be readily trapped in highly oblate geometries, where the dynamics of the vortices are restricted to a plane (*18*–*21*). The emergence of persistent currents has been observed in an annular BEC (*19*), and the relaxation of turbulence has been investigated via vortex number statistics (*20*). These early experiments relied heavily on comparisons to numerical simulations, as information about the vortex circulation could not be obtained. This provided motivation for developing new techniques to probe 2DQT, such as multishot vortex tracking (*22*, *23*) and single-shot vortex sign detection (*21*, *24*). In the latter case, a velocity-selective Bragg spectroscopy (*25*) technique was utilized, allowing the sign and position of each vortex (and hence the incompressible velocity field) in a turbulent BEC to be determined (*21*). However, the BEC was harmonically trapped and therefore inhomogeneous, and the vortices preferentially formed vortex–antivortex dipole pairs. Indeed, numerical studies have indicated that the uniformity of the BEC plays a major role in the turbulent dynamics (*26*).

Here, we injected vortices into a uniform, planar BEC by dragging a grid of elliptical obstacles formed by an array of laser beams through the atomic cloud (*18*) and observed the evolution of the resulting states by identifying the sign and location of each vortex in the BEC (*21*). Using velocity-selective Bragg scattering, as implemented in (*21*), we generated a map of the BEC, highlighting regions where the 1D projection of the superfluid velocity is resonant with a stimulated scattering process from a pair of counterpropagating laser beams due to the Doppler shift (Fig. 1, D to F). The resulting momentum transfer to the atoms from the laser beams only scatters a small number of atoms, and density holes visible in the image of the unscattered component of the BEC correspond to the locations of the vortices (Fig. 1, A to C). Because of the circulating flow, the Bragg-scattered signal is antisymmetric across vortices, allowing their signs to be determined from this 1D velocity information (*21*, *27*) and a numerical 2D velocity map to be generated (Fig. 1, G to L). We then used a vortex classification algorithm (*16*, *28*) to identify vortices as clusters, dipole pairs, or free vortices. Any two like-sign vortices are said to belong to the same cluster if they are closer to one another than either is to an opposite-sign vortex, whereas a vortex–antivortex pair is defined as a dipole if they are mutual nearest neighbors. Free vortices are those that are left over after all clusters and dipoles have been assigned. By changing the semimajor axis length *R*_{G} of each grid obstacle (Fig. 2), we were able to control the initial spacing between vortices of opposite sign, which are preferentially shed from opposite sides of each barrier. This in turn allowed us to control the amount of kinetic energy imparted to the BEC, as this energy increases with the increasing separation of vortex-antivortex pairs. For the finest grids used, the resulting vortex distributions are dominated by dipole pairs, whereas for larger grids, clustered vortices form the majority, with the ratio changing monotonically with obstacle size (Fig. 2, A and B).

To analyze the observed vortex distributions in terms of Onsager’s thermodynamic framework, we first assigned to them a vortex temperature. Because the turbulent BEC is in a highly nonequilibrium state, the vortex subsystem itself must be well isolated from the embedding fluid (in this case, the phonon bath), and also be in a state of quasiequilibrium, for the temperature to be a valid observable. It is reasonable to assume that this was the case here because the time scale on which vortices move (and redistribute themselves) is much shorter than the timescales on which vortex–sound interactions become important (*27*). To assign a temperature to our states of 2DQT, we applied a method of vortex thermometry (*29*) that uses a monotonic relationship between the vortex temperature and the (mean) populations of clusters (*p*_{c}) and dipoles (*p*_{d}) at a given vortex number *N*_{v}. By performing Monte Carlo (MC) simulations in the canonical ensemble, we obtained thermometry curves of the mean populations *27*). The vortex temperature of the experimental data was then measured by first taking an ensemble average of *p*_{c} and *p*_{d} for each grid size and hold time and then finding the temperature at which these two fractions simultaneously matched best with the corresponding MC curves [chosen such that *N*_{v}; see (*27*) for details], as demonstrated in Fig. 3A. Figure 3B shows an example of the temperature assignment for all configurations for which the*p*_{c}, *p*_{d}, *p*_{f}} points corresponds to a single grid size and hold time, and the corresponding β axis reading gives the obtained temperature of that configuration (although the free vortex fraction *p*_{f} was not fitted to the thermometer, it was compared as a consistency check). Data for all other configurations are shown in fig. S1.

To corroborate the vortex temperature measurements, we calculated a number of additional independent observables that quantify the vortex clustering (or lack thereof). The first of these is the correlation function *i*th vortex has the same (opposite) sign. For a configuration dominated by dipoles (clusters), the sign of *C*_{1} is negative (positive), and in the limit of large *N*_{v}, *C*_{1} vanishes for a random configuration (near β = 0). We also calculated the mean vortex dipole moment, *s _{i}* = ±1 is the sign of circulation of the

*i*th vortex located at position

*r*relative to the center of the trap. Finally, we numerically generated a spectrum of the incompressible kinetic energy per vortex for each configuration (

_{i}*27*,

*30*), which was integrated to give the total incompressible kinetic energy per vortex. The cluster fraction, correlation function, dipole moment, temperature, and energy per vortex all showed an increasing trend as the grid was coarsened (Fig. 2), indicating an increase in long-wavelength energy in the system (fig. S2).

Time-series data for each grid are shown in Fig. 4 and figs. S3 to S6. For the 6-μm grid (Fig. 4), the distribution began with approximately equal weightings of clusters, dipoles, and free vortices, corresponding to an inverse temperature close to 0 (*29*). During the hold time after the grid sweep, the vortex distribution evolved to a state with a higher clustered fraction (Fig. 4A) and to a negative temperature (Fig. 4D). The correlation function evolved from negative to positive values (Fig. 4B) and the dipole moment grew (Fig. 4C). The incompressible kinetic energy spectrum showed a buildup of energy at scales on the order of the system size (Fig. 4F), leading to an increase in the mean energy per vortex over time (Fig. 4E). These observations point toward the existence of an inverse energy cascade driving the system toward ordered Onsager vortex states. Similar behavior is observed for the finest grid (4.2 μm, fig. S3), though the consistently lower vortex number leads to larger fluctuations in the classified populations and hence vortex temperature.

Although the inverse energy cascade is only predicted to occur in 2D classical systems with continuous driving (*5*, *31*), large-scale flows have been seen to form in numerical simulations of 2DQT with vortex–antivortex annihilation but no continuous driving (*17*, *26*). For a given distribution of vortices and antivortices, the annihilation of vortex–antivortex pairs can be thought of as a quench that causes “evaporative heating” (*17*, *29*)—the removal of the “coldest” vortices, which contribute the least to the total incompressible kinetic energy associated with the vortex flow field. The “hotter” remaining vortices redistribute, with an increased mean energy per vortex. This process is analogous to the mechanism of evaporative cooling, whereby the hottest particles are lost, resulting in a cooler distribution. There is evidence that this heating process occurred for the 6-μm grid (Fig. 4): The clustered vortex fraction, dipole moment, and energy per vortex (in particular at low wavenumbers) all increased over time.

In the configurations that began with an appreciable clustered fraction, we did not observe significant evolution of the classified vortex fractions (Fig. 2B and figs. S4 to S6). We believe that evaporative heating was suppressed in such configurations because of the small number of dipole and free vortices available to annihilate. Indeed, in these cases, the time evolution of the energy per vortex indicates that some cooling of the vortices occurred (see gray lines in Fig. 2F). As there is no appreciable evaporative heating in this regime, the cooling effects of phonons and the thermal background dominate the vortex thermodynamics (the same effect is also evident at the latest times for the 4.2- and 6-μm grid sweeps; see Fig. 4 and fig. S3). We expect that for higher vortex numbers with a similar cluster fraction, evaporative heating would in fact play an important role, as the correspondingly higher number of dipole vortices remaining would allow for further annihilations, resulting in larger Onsager vortex clusters as seen in numerical studies (*17*, *26*).

For all grid sweeps, the vortex number gradually decreased over time during the turbulent decay. Many theoretical models have been put forward to statistically describe this vortex annihilation behavior in terms of *n*-vortex collision events (e.g., *20*, *26*, *32*). For example, it was recently predicted (*32*) that vortex annihilations should occur primarily via two- and three-body loss processes, and that these processes should give rise to number decay equations of the form *N*_{v}(*t*) ∝ *t*^{–1} and *N*_{v}(*t*) ∝ *t*^{–2/5}, respectively. These predictions are displayed alongside the experimental vortex number data in Fig. 4G and figs. S3 to S6, showing reasonable agreement for the 4.2- and 6-μm grids. The coarser grids appear to maintain the slower decay rate to longer hold times. As further support for this theoretical model, we also plotted the mean nearest-neighbor intervortex scale *l*_{v} in Fig. 4F and figs. S3 to S6. All grids appear to be consistent with *l*_{v} ∝ *t*^{1/5} scaling initially, with the finer grids transitioning to *t*^{1/2} scaling, in accordance with the observed number decay behavior (*32*). We emphasize that our data cover only a small range of scales and although they appear consistent with the scaling laws, additional data in larger systems would be required to confirm these processes.

The ability to measure statistical distributions arising from the dynamics of vortices will facilitate studies of nonequilibrium behavior in isolated quantum systems. Such studies may include experimental investigations of quench protocols and various thermalization processes such as the eigenstate thermalization hypothesis (*33*), as well as nonthermal fixed points (*34*) in quantum turbulence and their relation to universal phenomena in far from equilibrium systems. Indeed, the *l*_{v} ∝ *t*^{1/5} scaling consistent with our data are predicted to correspond to an anomalous fixed point, which can be interpreted as a critical slowing of the system’s dynamics (*32*).

Our results corroborate Onsager’s statistical description of 2D turbulence. The group of T. W. Neely has independently and simultaneously observed negative absolute temperature states of quantized vortices in a similar system in the form of large, persistent clusters (*35*).

## Supplementary Materials

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We thank T. W. Neely, M. J. Davis, M. T. Reeves, G. Gauthier, A. S. Bradley, X. Yu, T. A. Bell, H. Rubinsztein-Dunlop, M. Baker, and T. P. Billam for useful discussions.

**Funding:**We acknowledge financial support from the Australian Research Council Discovery Program (project nos. DP130102321, DP170104180, and FT180100020). This research was also partially supported by the Australian Research Council Centre of Excellence in Future Low-Energy Electronics Technologies (project no. CE170100039).

**Author contributions:**S.P.J., T.P.S., and K.H. conceived the experiments. S.P.J. conducted the experiments. S.P.J., P.T.S., C.J.B., and K.H. built the experimental apparatus and developed the experimental control and analysis software. A.J.G. ran the Monte Carlo simulations and assisted with data analysis. S.P.J., A.J.G., T.P.S., and K.H. wrote the manuscript. All authors discussed the experiments and commented on the manuscript.

**Competing interests:**The authors declare no competing interests.

**Data and materials availability:**All experiment data, analysis, and simulation code are available on figshare (

*36*).