Onsager vortex formation in Bose–Einstein condensates in two-dimensional power-law traps

Onsager vortex formation in Bose–Einstein condensates in two-dimensional
power-law traps

Andrew J. Groszek, Tapio P. Simula, David M. Paganin and Kristian Helmerson
School of Physics and Astronomy, Monash University, Victoria 3800, Australia

We study computationally dynamics of quantised vortices in two-dimensional superfluid Bose–Einstein condensates confined in highly oblate power-law traps. We have found that the formation of large scale Onsager vortex clusters prevalent in steep-walled traps is suppressed in condensates confined by harmonic potentials. However, the shape of the trapping potential does not appear to adversely affect the evaporative heating efficiency of the vortex gas. Instead, the suppression of Onsager vortex formation in harmonic traps can be understood in terms of the energy of the vortex configurations. Furthermore, we find that the vortex–antivortex pair annihilation that underpins the vortex evaporative heating mechanism requires the interaction of at least three vortices. We conclude that experimental observation of Onsager vortices should be the most apparent in flat or inverted-bottom traps.

03.75.Lm, 03.75.Kk, 67.85.-d, 67.25.dk
preprint: DOI: 10.1103/PhysRevA.93.043614

I Introduction

Non-equilibrium physics of quantum gases has attracted significant activity recently Giamarchi et al. (2012). Quantum turbulence (QT) is an archetype of non-equilibrium dynamics, which manifests as a chaotic motion of large numbers of quantised vortices and features an intriguing interplay between chaos and order. In three-dimensional QT, vortex filaments form tangles—a phenomenon which has been widely studied but only imaged directly in recent years in superfluid helium Bewley et al. (2006); Paoletti et al. (2008); Guo et al. (2014). Remarkably, despite the fact that the microscopic behaviour of three-dimensional QT is driven by Kelvin waves Thomson (1880); Bretin et al. (2003); Simula et al. (2008a); Koens et al. (2013); Fonda et al. (2014), Crow instabilities Crow (1970); Berloff and Roberts (2001); Simula (2011), vortex reconnections Feynman (1955); Schwarz (1985); Fonda et al. (2014), phonon radiation Samuels and Barenghi (1998); Barenghi et al. (2005) and mutual friction between the normal and superfluid components Hall and Vinen (1956), statistically the dynamics is thought to yield the same Kolmogorov scaling of incompressible kinetic energy as in classical fluid turbulence.

Restricting the motion of the quantised vortices in one of the spatial dimensions results in two-dimensional (2D) quantum turbulence, in which the vortex tangle reduces to a chaotic configuration of point-like vortices. Recent studies have focused on observing the decay of QT in Bose–Einstein condensates, both experimentally Neely et al. (2013); Kwon et al. (2014); Seman et al. (2011) and using computer simulations Parker and Adams (2005); Nazarenko and Onorato (2006); Numasato and Tsubota (2009); Numasato et al. (2010); Schole et al. (2012); Reeves et al. (2013); Simula et al. (2014); Billam et al. (2014, 2015); Stagg et al. (2015). Two-dimensional systems have attracted particular interest due to a prediction of an inverse energy cascade Kraichnan (1967, 1971) from small to large spatial scales, which originates from the theory of classical fluid turbulence. Using a statistical model of point-vortices, Onsager Onsager (1949) predicted for such systems the emergence of large-scale vortex structures, such as those seen in geophysical systems Eyink and Sreenivasan (2006). Similarly, in 2D QT, the inverse cascade is anticipated to lead to the clustering of like-sign vortices into large-scale Onsager vortex structures.

Recent experimental advances in producing Seman et al. (2011); Wilson et al. (2013); Kwon et al. (2015); Kang et al. (2015), imaging Wilson et al. (2015) and controlling Samson et al. (2016) quantised vortices in ultracold atomic dilute gas Bose–Einstein condensates (BECs) have resulted in detailed measurements of vortex dynamics in these superfluid systems Freilich et al. (2010); Kwon et al. (2014); Navarro et al. (2013); Neely et al. (2010). However, the debate continues regarding whether or not an inverse cascade and associated Onsager vortices should emerge in compressible 2D QT Neely et al. (2013); Billam et al. (2014, 2015); Reeves et al. (2013); Simula et al. (2014); Reeves et al. (2013); Chesler and Lucas (2014); Kwon et al. (2014); Stagg et al. (2015); Numasato et al. (2010); White et al. (2012). Numasato et al. Numasato et al. (2010) simulated quantum turbulence in a uniform 2D superfluid and found evidence of a direct cascade pushing incompressible kinetic energy towards small length scales. In accordance with this finding, a recent experiment Kwon et al. (2014) and simulation Stagg et al. (2015) of a turbulent harmonically trapped highly oblate BEC did not find evidence for the formation of Onsager vortices. By contrast, Simula et al. Simula et al. (2014) observed strong evidence of vortex clustering in their quasi-2D simulations in a flat trap with steep walls.

One key difference between these studies which could explain the disparity between their findings is the trapping potential used for confining the condensate. The aim of this paper is therefore to investigate the role of the trap geometry with regard to the emergence of Onsager vortices. We focus on numerical studies of decaying two-dimensional quantum turbulence in power-law traps, with a particular emphasis on comparing harmonically trapped condensates to those in uniform disk potentials with steep walls. A variety of techniques exist for producing such steep-walled trapping potentials experimentally Gaunt et al. (2013); Chomaz et al. (2015); Lee et al. (2015).

We simulate BEC dynamics using a Gross–Pitaevskii model and also study their thermodynamic properties using a Markov Chain Monte Carlo technique, interpreting the vortex dynamics in each trap in terms of vortex evaporative heating Simula et al. (2014). In addition, we examine in detail the microscopic process of vortex–antivortex annihilation, an essential aspect of the decaying turbulence in these systems. In Sec. II, we introduce the numerical model and computational techniques. In Sec. III, we present the key findings from our simulations of decaying superfluid turbulence in different trapping potentials and interpret our observations using a statistical mechanics framework. We then examine the vortex dynamics on a microscopic scale, focusing in particular on vortex–antivortex annihilation in 2D QT, showing it to be a many-vortex process. Finally, Sec. IV is devoted to discussion.

Ii Model

ii.1 System parameters

The vortex dynamics of two-dimensional dilute gas Bose–Einstein condensates are inherent in the time-dependent condensate wavefunction , whose evolution is modelled here using the Gross–Pitaevskii equation (GPE):


where is the mass of an atom, and is the trapping potential which radially confines the condensate. The effective interaction parameter accounts for reduction of dimensionality of the Gross–Pitaevskii equation from three to two, where is the normalised axial wavefunction, is the total number of atoms in the condensate, and is the interaction parameter for the three-dimensional system being modeled, defined in terms of the -wave scattering length . For a uniform cylindrical trap, , where is the axial length of the three-dimensional condensate. In a sufficiently tight axial harmonic trap, the axial wave function is well approximated by a Gaussian, and hence , where is the axial harmonic oscillator length with an effective harmonic trapping frequency . The wavefunction is normalized such that .

We consider general power-law trapping potentials


where is the radial harmonic trapping frequency, is a parameter which defines the steepness of the trap walls, and is the effective system radius. For this potential is a standard harmonic trap with Thomas–Fermi radius . In the limit of infinite steepness () it approaches a cylindrically symmetric well of radius .

Our system parameters correspond to a two-dimensional BEC with a radial trapping frequency of , and a Thomas–Fermi radius of , where the radial harmonic oscillator length scale is defined as . To this end we choose . Hence, the radial extent of our system is similar to those used in the recent experiment by Kwon et al. Kwon et al. (2014) and simulations by Stagg et al. Stagg et al. (2015).

ii.2 Numerical techniques

In the beginning of our simulations we solve for the approximate ground state of the system using imaginary time evolution of the GPE. We then imprint vortices in the condensate by multiplying the ground state wavefunction by a phase factor , where . Here, the co-ordinate defines the position of the th vortex, whose circulation sign is . We imprint equal numbers of vortices () and antivortices ().

We choose initial conditions which approximate high entropy, highly randomised states which could be produced by stirring the condensate. To this end, we first construct a density of states distribution for our chosen vortex number by iteratively generating random vortex configurations and calculating their energy using a point-vortex Hamiltonian Simula et al. (2014). The maximum entropy state corresponds to the peak of this distribution; hence, we ensure that all initial conditions generated have an energy lying within 10% of this maximum entropy value 111The evaporative heating mechanism does not rely on starting with a specific vortex configuration—the initial condition simply determines how much heating is required to reach the clustered Onsager vortex states..

After the vortex imprinting step, the wavefunction is evolved further in imaginary time for to establish the structure of the vortex cores. This can lead to the annihilation of vortices near the boundary, as well as vortex–antivortex pairs if they were imprinted very close together. The number of vortices at the start of the real-time evolution is therefore on average seven fewer than the that were originally imprinted.

We solve the GPE using a fourth-order split-step Fourier method on a spatial grid with spacing (approximately 0.65 condensate healing lengths) unless otherwise stated. The locations of the vortices in the system are detected by measuring the positions of the phase singularities in the wavefunction at predetermined time intervals. The direction of the phase winding about each singularity determines the circulation sign of the vortex. The vortex locations are only measured in a region in order to avoid detection of ghost vortices Tsubota et al. (2002) in the low density region of the traps with lower values.

Iii Results

iii.1 Macroscopic dynamical behaviour

We first compare the results of decaying turbulence in the two traps discussed in the literature: a harmonic trap () and a uniform trap with steep walls (), which has constant density to within healing lengths of the boundary. For each simulation, we monitor the number of vortices , which decreases over time due to vortex annihilation events. We also measure the dipole moment of the vortex distribution, defined as , where is the position of the ith vortex, and is its charge. For the confined systems being studied here, it is convenient to scale with the system size and the number of vortices . If the vortices are randomly distributed, will approach zero for large systems. A large , on the other hand, signals the presence of two Onsager vortex clusters in our system.

Figure 1 shows the characteristic time evolution of the vortex distribution in the two traps, along with the respective dipole moments. In agreement with previous simulations and experiments Stagg et al. (2015); Kwon et al. (2014), we observe no significant vortex clustering in the harmonic trap. However, and also in agreement with previous 3D simulations Simula et al. (2014), the uniform trap exhibits a strong tendency to form Onsager vortices, as indicated by the increasing dipole moment. Thus, we conclude that the shape of the trapping potential has a strong influence on the vortex clustering behaviour, partially resolving the apparent contradiction in the existing literature.

Figure 1: Comparison of the time evolution of the vortex configuration between the harmonic trap (a)–(c) and uniform trap (d)–(f). The grayscale value represents the superfluid density, and the colorbars are normalised to the maximum density: and for the top and bottom rows, respectively. Vortices and antivortices are denoted by blue (dark) and green (light) circles, respectively. The red line represents the effective dipole moment of the vortex distribution. Movies S1 and S2 in the Supplemental Materials sup show the dynamics of each simulation.

iii.2 Statistical mechanics interpretation

The spontaneous formation of Onsager vortices found in Ref. Simula et al. (2014) was attributed to the evaporative heating mechanism of vortices. When a vortex–antivortex pair annihilates, the incompressible kinetic energy of the system is redistributed among the vortices remaining in the system. This process can lead to evaporative heating of the vortex gas, whereby the mean energy per vortex increases each time an annihilation occurs. When the mean energy per vortex crosses a critical value, a transition into the Onsager vortex state is possible Simula et al. (2014).

Figure 2: Comparison of the vortex number decay and dipole moment evolution (inset) for the harmonic (red/light) and uniform (blue/dark) traps. The black circles in the inset correspond to the timeframes displayed in Fig. 1. The fluctuations in the vortex number are due to vortices crossing the counting radius of , in addition to occasional vortex–antivortex pair creation.

The absence of strong clustering in the harmonic trap could be due to (i) the rate of evaporative heating per annihilation event being too low, leading to inefficient evaporative heating of the vortex gas, (ii) the critical energy per vortex for the Onsager vortex transition in a harmonic trap being out of reach despite the vortices being evaporatively heated, or (iii) the critical value of the dipole moment for harmonic traps being too small to allow a clear distinction to be made between the disordered and clustered vortex configurations. In the following we argue that the combined effect of (ii) and (iii) may explain the observed behaviour.

iii.2.1 Dynamical statistical behaviour

Figure 3: Comparison of statistical behaviour between (a) the harmonic trap and (b) the uniform trap. For each dynamical simulation, the dipole moment of the vortex configuration is shown as a function of the incompressible kinetic energy per vortex number squared. The initial state in each plot is the bottom-left corner, and the evaporative heating increases the energy per vortex number squared over time. The data appears as columns because each vortex annihilation increases the energy per vortex number squared by a discrete amount.
Figure 4: Statistical data obtained from 100,000-step Markov Chain Monte Carlo simulations for a harmonic (red/light) and a uniform (blue/dark) trap for a total of 12 vortices with equal numbers of vortices and antivortices. The subfigures show (a) the specific heat, (b) the incompressible kinetic energy per vortex number squared, and (c) the dipole moment of the configuration, each plotted as a function of the statistical temperature. The shaded regions in (b) and (c) correspond to the standard deviation of each observable at a given temperature. The maximum in the specific heat indicates the transition to the Onsager vortex state in each trap, and is accompanied by an increase in both the energy per vortex number squared and the dipole moment. Frames (d) & (e) and (f) & (g) show typical vortex configurations at the temperature extremes in the harmonic and uniform traps, respectively, with labelling as in Fig. 1. The temperatures shown in these frames are indicated in (c) with vertical dashed lines.

Figure 2 shows little difference between the vortex number decay in the two traps. This suggests that the evaporation of vortices is only weakly affected by the details of the trapping potential. However, the dipole moment shows quantitatively different behaviour between the two traps, and indicates strongly enhanced clustering in the uniform trap. To better understand this difference, we construct a probability distribution of different vortex configurations generated by the dynamics in the space spanned by the dipole moment and energy per vortex number squared by taking the vortex configuration at each time-step to correspond to an independently sampled microstate. We choose to normalise the energy to the square of the vortex number to cancel out the scaling which occurs in the Onsager limit, as the system tends towards a multi-quantum vortex dipole configuration. Figure 3 shows the resulting histograms for each trap. In the harmonic trap (a), the dipole moment shows no significant variation over the measured range of energy per vortex number squared, and hence there is no evidence that the system crosses the Onsager vortex transition. Conversely, the trend in the uniform trap (b) is a clear indication that the evaporative heating is on average increasing the dipole moment, causing the system to evolve towards the Onsager vortex state.

iii.2.2 Monte Carlo thermodynamics

In order to determine the statistical behaviour of the vortex gas beyond the range accessible via the dynamics, we implement a Markov Chain Monte Carlo (MCMC) algorithm for the two traps on a grid. The algorithm is initialised by imprinting a random configuration of vortices into the condensate ground state using the imaginary time propagation method described in Sec. II.2. We set (six vortices of each sign) to approximate the late time configurations of the dynamical simulations presented in Figs. 1 and 2. Keeping fixed, each step in the algorithm shifts a single randomly chosen vortex in the configuration and calculates the value of a predetermined weighting function . This new configuration is then either accepted or rejected based on the change in the weighting function. Here, we use a Boltzmann factor as our weighting function, defining as the statistical temperature of the vortex gas (which in this case is negative—see Refs. Kraichnan and Montgomery (1980); Simula et al. (2014); Viecelli (1995) for details). The choice of determines the most probable vortex configuration. Hence, we can vary to alter the statistical behaviour of the system—this is the basis of the evaporative heating interpretation of the dynamics. To characterise the temperature dependence, we measure three observables: the energy per vortex number squared, the dipole moment and the specific heat, defined as . The system is evolved for 110,000 Monte Carlo steps, the first 10,000 of which are disregarded as the initial condition is, in general, unrepresentative of the chosen temperature. The results for both traps are shown in Fig. 4. This MCMC data shows the transition from the disordered state to the Onsager vortex state in each trap, characterised most obviously by a maximum in the respective specific heat curves in Fig. 4(a). In addition, both the energy per vortex number squared and dipole moment begin to rapidly climb around this critical temperature, signalling the formation of vortex clusters. For a uniform system with superfluid density , the critical temperature is predicted to be Kraichnan and Montgomery (1980); Simula et al. (2014). For vortices, this yields a critical temperature of , which agrees well with our data. In a harmonically trapped system, Fig. 4 shows that will be shifted towards lower temperatures compared to the uniform system.

The key differences between the two traps are evident in Fig. 4. Figure 4(c) shows that the dipole moment climbs to a significantly higher value at the highest temperatures in the uniform trap compared to the harmonic trap—the respective vortex configurations are displayed in frames (e) and (g). In fact, the dipole moment shows only a weak temperature dependence in the harmonic trap, the most marked effect being a decrease in its variance at high temperatures. This suggests that, even if the harmonically trapped system transition to the Onsager state, the resulting dipole moment would remain relatively small when compared to the steeper traps. Figure 4(b) also shows that the energy per vortex number squared required to cross the transition is significantly higher in the harmonic trap. This provides further support for the absence of clustering in the GPE dynamics in the harmonic trap, as the evaporative heating does not supply enough energy to drive the system to these temperatures.

Figure 5: Incompressible kinetic energies of a vortex–antivortex pair for a range of power-law traps computed using the GPE. The pair is placed symmetrically in the trap and both vortices are an equal radial distance from the center. In order of peak location from left to right, the power-law exponents are (red), (thin black lines) and (blue). In addition, the dipole energy for an inverted trap (as described in the main text) is shown in light green (far right peak). The maximal separation is indicated on each curve with a circle, and is emphasised further on the two extreme power-law traps, as well as the inverted trap, with a vertical dashed line.

iii.2.3 Maximum achievable dipole moment

We can predict numerically the maximal separation of the two Onsager vortex clusters in a given system by calculating the energy of a vortex dipole as a function of the separation of the vortex and the antivortex. This yields further insight as to why the two traps show different clustering behaviour. In an infinite system, increasing the dipole separation will logarithmically increase the energy of the Onsager dipole without bound. However, for a bounded system, there exists a separation which maximises the energy. For a harmonic trap, this maximum energy configuration also corresponds to a stationary state Freilich et al. (2010); Middelkamp et al. (2010); Möttönen et al. (2005); Kuopanportti et al. (2011). The dipole energy landscapes obtained for various trap steepnesses are presented in Fig. 5, showing that the energy maximising separation increases as a function of the steepness. This result explains why the MCMC dipole moments in Fig. 4(b) asymptote to different values in the high temperature limit, as the two systems reach their highest energy at differing cluster separations. In addition to various power law traps, Fig. 5 shows the energy in an ‘inverted’ trap. This trap consists of a steep wall () and an additional repulsive Gaussian potential with a width of which causes the condensate density to dip in the center. By pushing the fluid radially outwards, the energy maximising separation of a vortex dipole increases significantly, suggesting that an Onsager state in this trap should have an even greater dipole moment than that in the trap. We have confirmed this prediction with a dynamical GPE simulation presented as Movie S3 in the Supplemental Materials sup .

iii.3 Vortex annihilation is a many-vortex process

The microscopic underpinning of the evaporative heating mechanism of vortices is vortex–antivortex annihilation Simula et al. (2014). Scalar Bose–Einstein condensates with quantised vortices have two types of low-lying excitations—Bogoliubov phonons and vortex waves Simula et al. (2008b); Simula and Machida (2010); Simula (2013). Such modes can resonate, mediating vortex–sound interactions Leadbeater et al. (2001); Zuccher et al. (2012). In principle such vortex–phonon interactions could cause vortex–antivortex pairs to annihilate via soundwave emission, which would account for the conservation of energy and momentum. However, for a single vortex–antivortex pair this does not occur as has been supported experimentally Neely et al. (2010) and shown theoretically Lucas and Surówka (2014). If such vortex–antivortex pair annhilations are forbidden, this raises the question of how the vortex number can decay over time as observed both in the simulations and experiments Kwon et al. (2014).

Figure 6: A vortexonium state (formed from a vortex–antivortex pair) highlighted in (a) with a dashed oval colliding with an antivortex and dissipating into fluid soundwaves which disperse radially, indicated in (b) and (c) with dashed circles. The vortices and colorbar are labelled as per Fig. 1. Supplemental Movie S4 shows the dynamics of this event sup , including the wavefunction phase.

The answer must be that vortex–antivortex annihilation is a many-vortex process. Figure 6 shows a three-vortex process whereby a vortex–antivortex pair has formed a neutral vortexonium state (a rarefaction pulse), in which the individual vortex phase singularities are no longer discernible yet the excitation retains its identity as a spatially localised bound state. This excitation is reminiscent of positronium—a neutral bound state of an electron and a positron. The vortexonium, which is identifiable by a phase step, travels close to the speed of sound until it eventually scatters off an additional vortex or antivortex, as shown in Fig. 6(b) and (c). This decay process irreversibly disperses the energy and momentum of the vortexonium into sound waves Nazarenko and Onorato (2006); Smirnov and Smirnov (2015). Until this secondary process occurs, the vortexonium can also re-form as a vortex–antivortex pair, an event which frequently occurs when a vortexonium state travels into the low density region near the boundary of the trap. The formation of vortexonium as a precursor to the vortex–antivortex annihilation process in 2D BECs has been discussed previously Prabhakar et al. (2013); Kwon et al. (2014); Stagg et al. (2015); Du et al. (2015). Here, we identify the three-vortex collision to be an essential part of the annihilation process in 2D superfluid turbulence.

Figure 7: Feynman diagram depicting the entire vortex–antivortex annihilation process observed, with time flowing from left to right. The straight lines represent vortices () and antivortices (), the double line represents vortexonium (), and the wavy line denotes the sound waves emitted at each vertex (the magnitude of the second burst of sound is far greater than the first). The light blue lines indicate participating catalyst vortices, which are not annihilated during the process.

The question remains how these vortexonium states form to begin with. In a uniform system free from dissipation, an isolated vortex–antivortex pair will travel with constant velocity and inter-vortex separation. Therefore, some mechanism other than sound induced interaction must be responsible for reducing the pair’s separation and forming vortexonium. In our simulations, we observe two ways this bound state can form. Firstly, a vortex–antivortex pair travelling towards a higher density region will reduce its separation in order to satisfy energy conservation, often forming a vortexonium state. However, this process only occurs in traps with shallow walls, where the density variation is significant. The second process we observe is the shrinking of a vortex–antivortex pair via a long-range interaction with a third vortex. By giving up some of its energy to this catalyst vortex, the pair can reduce their separation, ultimately resulting in a vortex–antivortex fusion event and the formation of a vortexonium state. We note that this latter process is ubiquitous in all traps studied. However, in the presence of dissipation, both the formation and annihilation of vortexonium would be possible without additional interactions, as the loss of energy would gradually drive vortex dipoles closer together regardless (see Sec. III.4).

Combining these observations, we obtain a picture of the vortex–antivortex annihilation process, depicted as a Feynman diagram in Fig. 7 which shows how four vortices are involved in the annihilation. Movie S4 in the Supplemental Materials sup shows one such four-vortex process. In the first stage a vortex–antivortex pair interacts with a catalyst vortex to produce a vortexonium state and in the second stage the vortexonium scatters off a catalyst vortex leading to the ultimate destruction of the vortex–antivortex pair and the emission of sound. The catalysts can be any vortex or antivortex in the system.

iii.4 Rate equation for evaporative heating of vortices

Attempts have previously been made to fit a universal law to the vortex number decay Kwon et al. (2014); Stagg et al. (2015); Schole et al. (2012); Cidrim et al. (2016). Kwon et al. Kwon et al. (2014) proposed a phenomenological model of the form , comprised of a linear term to model vortex drift out of the condensate and a nonlinear term to account for vortex–antivortex annihilation, where the and are the one-body and two-body decay constants, respectively.

We find that, due to the zero temperature of the GPE simulations, this equation does not provide an adequate fit to our vortex number decay curves. Instead, for , the vortex number decay is well described by a power law of the form in all traps. This is evident in Fig. 8, which shows the number decay in a harmonic trap averaged over five simulations at resolution. This power law was also observed by Schole et al. Schole et al. (2012), who further suggested that the vortex number rate equation should have the form . This would reflect the importance of a four-body loss process at zero temperature, in contrast to the one- and two-body loss observed in Kwon et al.’s experiments Kwon et al. (2014). The four-vortex annihilation events discussed in Sec. III.3 are consistent with this four-body loss mechanism.

To study the effect of the thermal cloud, we model non-zero temperatures using a damped Gross–Pitaevskii equation Gardiner et al. (2001):


where is the temperature dependent dimensionless damping parameter, and is the chemical potential. We propose a general rate equation for vortex loss at all temperatures:


where is the decay constant corresponding to a particular -body loss mechanism. This model combines the one- and two-body loss processes observed in experiments Kwon et al. (2014) with the higher order three- and four-vortex loss processes observed in our zero temperature simulations. Strictly, a three-vortex decay process is not possible since it would violate the vortex charge conservation law. We instead interpret the three-body term as the loss of two vortices arising from the collision of three (i.e. a vortexonium colliding with a catalyst vortex, as discussed in Sec. III.3).

We have chosen the damping parameter to study the vortex number decay behaviour at non-zero temperature. Figure 8 shows the decay curves for zero temperature () and non-zero temperature (), each averaged over five simulations in a harmonic trap using a numerical grid. We model both cases using Eq. (4). For the case, we find that the decay is best described by a one- and two-body model, with , and . These values are in good agreement with those found by Kwon et al. Kwon et al. (2014). By contrast, the case is best described by a three- and four-body decay model with decay constants , and . We conclude that the three- and four-body vortex loss processes are characteristic of zero temperature systems, and that one- and two-body events become dominant at sufficiently high temperature. Quantifying the transition between these two behaviours at intermediate temperatures is left for future study.

Figure 8: Ensemble averaged vortex number decay curves for harmonically trapped systems at zero temperature (, solid red/dark line) and non-zero temperature (, solid green/light line). The fits for each curve to Eq. (4) are shown as black dashed lines, with , , for the non-zero temperature case, and , , for the zero temperature case.

iii.5 Interaction between vortices and boundaries

In our harmonic trap simulations, the multi-vortex collision process described in Sec. III.3 is the only mechanism of vortex annihilation, excluding a small proportion of vortices which drift out of the condensate. By contrast, the presence of a hard boundary in the steeper traps allows for a number of additional phenomena relating to the dynamics and decay of vortex–antivortex pairs. In particular, we observe three distinct vortex–boundary collision processes, two of which give rise to additional vortex decay branches.

When a single vortex is near the boundary, it will pair up with its image vortex of opposite sign beyond the wall and travel around the circumference of the trap at high velocity. If the separation reduces sufficiently, this vortex–image pair can form a vortexonium with a phase step along the tangent of the wall. As this bound state travels around the boundary, it can either unbind and reform the vortex–image pair, or it can annihilate in much the same way as a vortexonium in the fluid bulk—by colliding with another vortex.

We observe a similar process involving the collision of a vortex–antivortex pair with the boundary. When the pair collides with the wall, it unbinds into two separate vortex–image pairs, which then travel around the boundary in opposite directions, as shown in Fig. 9(a)–(c). If travelling at high enough velocity, one or both of these new vortex–image pairs can form vortexonium excitations, which can then decay as described above. Often, the collision will be violent enough to cause one of the vortices in the initial pair to annihilate immediately, while the other one is left to travel around the boundary.

Figure 9: (a)–(c) Unbinding of a vortex pair at the boundary in the uniform trap; (d)–(f) reflection of a vortexonium state at the boundary in the uniform trap. The green arrows show the direction each excitation is travelling. The insets in (a) and (d) show the phase of the wavefunction in the corresponding frame, showing the two singularities in (a) and the phase step in (d). The soundwave produced by each collision event propagates outwards in (c) and (f). The colorbar is normalised to the maximum condensate density, as in Fig. 1. Movies S6 and S7 in the Supplemental Materials sup show each event in full.

If the initial conditions are such that the vortex–antivortex pair which is incident on the boundary has already fused and formed a vortexonium excitation, the collision dynamics become markedly different. Figure 9(d)–(f) shows that the vortexonium will not separate at the boundary, but rather reflect from it, reversing its propagation direction. This effectively changes the sign of the vortices in the bound state, and can be understood as an exchange of locations with the image vortices beyond the boundary. Effectively, the image vortexonium travels into the condensate, while the real vortexonium leaves.

Remarkably, for the steepest potentials, the proportion of vortices annihilated at the boundary (i.e. via one of the first two processes described above) accounts for approximately half of the total vortex loss. Despite this clear spatial dependence of annihilation behaviour which is absent in the harmonic trap, the vortex number decays at the same rate (see Fig. 2) and the evaporative heating does not appear to be any more efficient. It seems plausible that boundary annihilations would increase evaporative heating efficiency, as less energy should be lost per annihilation (as the energy of a vortex in the low density close to the system’s boundary is less than in the fluid bulk), leaving more for the remaining vortices. However, we have found no strong evidence of this effect.

iii.6 Onsager vortex formation as a function of trap steepness

We repeated our Gross–Pitaevskii simulations of decaying turbulence for a number of trap steepnesses ranging between the two extremes examined in Sections III.1 and III.2 by varying the value of in Eq. 2. Five GPE simulations were performed in each of the chosen trap geometries using a grid, and the dipole moment curves obtained for each steepness were combined by taking averages at each point in time. These averaged dipole moment curves are shown in Fig. 10. On average, a steeper trap produces a larger dipole moment and thus a greater separation of vortex charge. As predicted from energy considerations in Sec. III.2, an inverted trap produces even stronger clustering than any of the power-law traps. For the power-law traps, it appears that the clustering behaviour saturates beyond a steepness of . The dipole moments in Fig. 10 should be compared with their predicted maximum values shown in Fig. 5.

Figure 10: Comparison of dipole moment evolution in traps of varying steepness. Each curve is averaged over five simulations. The trap steepnesses are (from bottom to top) (solid red line), (dot-dashed black line), (dotted black line), (dashed black line), (solid blue line) and an inverted trap (solid green line).

Iv Discussion

We have studied decaying two-dimensional quantum turbulence using the Gross–Pitaevskii model. We have considered Bose–Einstein condensates confined in generic power-law traps which, in particular, enables a comparison to be made between vortex dynamics in harmonically trapped condensates and in condensates confined in (nearly) uniform density disk traps. When an initial random vortex configuration is left to decay, we find that in uniform traps the vortices and antivortices arrange into Onsager vortex clusters due to the evaporative heating mechanism posited in Ref. Simula et al. (2014). However, when a harmonic trapping potential is used, the emergence of Onsager vortices is not obvious—a finding which agrees with experimental observations Kwon et al. (2014). To verify that these results are not specific to our randomly sampled initial vortex configurations, we repeated our simulations in both traps using a repulsive Gaussian laser potential to stir the fluid and produce the initial state vortex configuration, as in Ref. Stagg et al. (2015). Considering both lateral and circular stirring motions, the qualitative vortex evaporative heating behaviour in the harmonic and uniform traps was unaffected. This result was expected since a turbulent system should rapidly forget its history, washing out any initial state dependence.

We also performed Monte Carlo calculations to study equilibrium vortex configurations in harmonic and uniform traps. These calculations showed that the transition from disordered vortex configurations to the clustered Onsager vortex states exist also in harmonic traps but the resulting vortex dipole moment is significantly smaller than for uniform traps, which partly explains why the Onsager vortex clusters have not been observed to emerge in harmonically trapped Bose–Einstein condensates.

To obtain an improved understanding of the vortex evaporative heating mechanism Simula et al. (2014), we carefully tracked the vortex–antivortex annihilation events in the simulations. At zero temperature, we found that vortex–antivortex pair annihilation in these quantum turbulent systems occurs via a combination of three- and four-body processes, involving one or two catalyst vortices, respectively, in addition to the annihilating pair. Firstly, a vortex-antivortex pair forms a vortexonium bound state by either travelling through varying density (three-body) or interacting with a catalyst vortex (four-body). In both cases, this bound state then has to interact with a catalyst vortex for it to irreversibly decay into phonons. Indeed, it has been shown both experimentally Neely et al. (2010) and theoretically Lucas and Surówka (2014) that an isolated vortex–antivortex pair is resistant to sound induced decay. By adding dissipation to the Gross-Pitaevskii model, we simulated a non-zero temperature system and found that the three- and four-body annihilation mechanisms become less important, and instead one- and two-body annihilation events begin to dominate, in agreement with experimental observations Kwon et al. (2014).

By considering power-law traps of varying steepnesses, we found that the vortex clustering tendency becomes stronger as the trap steepness is increased. Finally, we found that a locally and weakly anti-trapping potential Lewandowski et al. (2003); Coddington et al. (2004); Simula et al. (2005) should provide the most promising route to experimental observation of the emergence of the Onsager vortices, which could possibly be detected using the vortex gyroscope imaging technique proposed in Ref. Powis et al. (2014).

We are grateful for useful discussions with Matthew Davis, Thomas Gasenzer and Andy Martin. We acknowledge support from an Australian Postgraduate Award (A. G.), the Australian Research Council via Discovery Project No. DP130102321 (T. S., K. H.) and the nVidia Hardware Grant Program (T. S.).


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description