Evolution of large-scale flow from turbulence in a two-dimensional superfluid

See allHide authors and affiliations

Science  28 Jun 2019:
Vol. 364, Issue 6447, pp. 1267-1271
DOI: 10.1126/science.aat5793

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.

Science, this issue p. 1264, p. 1267


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 1/T=S/EkBβ, where 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 kB 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 E0. Onsager noted that for E > E0, 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) (1315), 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 (1821). 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 RG 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).

Fig. 1 Vortex configurations in the dipole, random, and clustered regimes.

(A to C) The locations of vortices are visible as dark spots in the optical density images of the BEC. The example distributions shown are dipole dominated, random, and cluster dominated, respectively. (D to F) The corresponding Bragg spectroscopy signals, with vortex locations denoted by crosses. (G to I) The computed point vortex velocity field projected onto the line defined by the directions of the Bragg spectroscopy laser beams. Colors in (D) to (I) indicate projections of the superfluid flow in the direction indicated by the arrow. (J to L) The classification of the vortices based on their signs and positions: vortices (antivortices) are indicated by blue (green) points; clusters by lines of the same color; dipoles are linked by red lines. Streamlines of the computed flow are shown in gray. Data in each row are taken from a single run of the experiment with (top) RG = 4.2 μm and 1-s hold time, (middle) RG = 6 μm and 2.5 s hold time, and (bottom) RG = 9.7 μm and 1 s hold time, where RG is the semimajor axis length of the obstacle grid. Note that vortices on the edge of the BEC cannot be detected reliably; however, these do not contribute appreciably to the superfluidflow, as seen in the Bragg signal.

Fig. 2 Generating grid turbulence at positive and negative absolute temperatures.

Vortex configuration data, time averaged across early (t ≤ 1.5 s, dark points) and late (t ≥ 4 s, light points) hold times after the grids have passed through the BEC, for each grid used. (A) The grid configurations used are shown one-third of the way through the 0.5 s sweep. White areas have high laser intensity, repelling the atoms, and move through the BEC in the directions indicated by yellow arrows on the 7.9-μm grid (see also fig. S7A). (B) Classified populations of clustered (blue), dipole (red), and free (green) vortices. (C) Nearest-neighbor vortex sign correlation function C1 (see main text). (D) Mean vortex dipole moment d (see main text), scaled by the trap radius R. (E) Inverse temperature of the vortices, β. Positive (negative) values are scaled by the critical temperature |βBKT| (|βEBC|) (27). (F) Incompressible kinetic energy (in units of ρκ2/4π, where ρ is the superfluid density, κ = h/m is the unit of circulation, h is the Planck constant, and m is the atomic mass). Data in (B), (C), (D), and (F) are the mean ± SEM calculated from 10 to 25 measurements at each hold time, though most error bars fall within the points plotted.

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 (pc) and dipoles (pd) at a given vortex number Nv. By performing Monte Carlo (MC) simulations in the canonical ensemble, we obtained thermometry curves of the mean populations pcMC(β) and pdMC(β) for NvMC={4,6,8,10,14,20}(27). The vortex temperature of the experimental data was then measured by first taking an ensemble average of pc and pd 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 NvMC is closest to Nv; see (27) for details], as demonstrated in Fig. 3A. Figure 3B shows an example of the temperature assignment for all configurations for which theNvMC=14 thermometer was used. A vertical triad of {pc, pd, pf} 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 pf was not fitted to the thermometer, it was compared as a consistency check). Data for all other configurations are shown in fig. S1.

Fig. 3 Temperature assignment.

Fractional vortex population curves generated by MC simulations (colored solid lines) are used as thermometers to assign a vortex temperature to the experimental data (filled circles). (A) For a given hold time and grid size (0.5 s and 11.5 μm shown here), the inverse temperature reading (βm, vertical gray line) is obtained by finding the minimum of Rc,d(β)=Δc(β)2+Δd(β)2 (black solid line) at which the mean fractions of both clustered (blue) and dipole (red) vortices best match the relevant MC thermometry curves. At each β, the difference between the value of the MC curve and the experimental data are given by Δc(β) and Δd(β) for clustered and dipole vortices, respectively (27). (B) The placement of all data points that fall on the Nv = 14 MC curves (i.e., 12 ≤ Nv < 17) are shown; data for other Nv are shown in fig. S1. Free vortex fractions (green) are not used in the matching process, but typically fall close to the MC curve once the temperature for that configuration has been assigned. The positive (negative) temperature axis is scaled by the critical temperature |βBKT| (|βEBC|). Points show the mean populations ± SEM (vertical bars), with temperature uncertainty (horizontal bars) based on the population uncertainty (see fig. S1H) (27).

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 C1=1Nvi=1Nvci, where ci=1 (1) if the circulation of the nearest neighbor of the ith vortex has the same (opposite) sign. For a configuration dominated by dipoles (clusters), the sign of C1 is negative (positive), and in the limit of large Nv, C1 vanishes for a random configuration (near β = 0). We also calculated the mean vortex dipole moment, d=|1Nvi=1Nvsiri|, where si = ±1 is the sign of circulation of the ith vortex located at position ri relative to the center of the trap. Finally, we numerically generated a spectrum of the incompressible kinetic energy per vortex for each configuration (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.

Fig. 4 Evidence of evaporative heating in the 6-μm grid data.

(A) Evolution of the classified populations of clustered (blue), dipole (red), and free (green) vortices over 5 s, starting from near equal weightings to more clustered states at later times. (B) Nearest-neighbor correlation function. (C) Dipole moment, scaled by the trap radius R. (D) Inverse temperature of the vortices, β (axis scaled as in Fig. 2E). (E) Total incompressible kinetic energy per vortex (in units of ρκ2/4π). (F) Incompressible kinetic energy spectrum (in units of ρκ2/4π). An increase in energy at wavenumbers on the order of the system size (k = π/R, solid gray vertical line) can be associated with the development of large-scale flow. For wavenumbers kξ > 1 (dotted gray vertical line), the spectrum appears to scale as k–3 (dotted black line), corresponding to the shape of vortex cores imprinted in the calculation, whereas k–1 (solid black line) corresponds to the far field of an isolated vortex (30). A k–5/3 scaling (dashed black line), shown for comparison, would be expected for an energy cascade in a driven classical system (3, 5). The grid scale k = π/RG is indicated by the dashed gray vertical line, where RG = 6.0 μm is the semimajor axis of the grid. The inset shows the difference ΔE(k)/Nv from the 0.5 s hold time spectrum. (G) Vortex number decay. Lines show t–2/5 (solid) and t–1 (dotted) power laws for reference. (H) Mean intervortex separation scaled by R. Lines show t1/5 (solid) and t1/2 (dotted) power laws for reference. Points in [(A) to (C)], (E), (G), and (H) are the mean ± SEM. Lines in (F) join the mean k-binned values. The SEM is on the order of the separation between lines.

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 Nv(t) ∝ t–1 and Nv(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 lv in Fig. 4F and figs. S3 to S6. All grids appear to be consistent with lvt1/5 scaling initially, with the finer grids transitioning to t1/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 lvt1/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

Materials and Methods

Figs. S1 to S10

References (3743)

References and Notes

  1. Materials and methods are available in the supplementary materials.
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).

Stay Connected to Science

Navigate This Article