Vinen turbulence via the decay of multicharged vortices in trapped atomic Bose-Einstein condensates

We investigate a procedure to generate turbulence in a trapped Bose-Einstein condensate which takes advantage of the decay of multicharged vortices. We show that the resulting singly-charged vortices twist around each other, intertwined in the shape of helical Kelvin waves, which collide and undergo vortex reconnections, creating a disordered vortex state. By examining the velocity statistics, the energy spectrum, the correlation functions and the temporal decay, and comparing these properties with the properties of ordinary turbulence and observations in superfluid helium, we conclude that this disordered vortex state can be identified with the `Vinen' regime of turbulence which has been discovered in the context of superfluid helium.

The singular nature of quantized vorticity (concentrated along vortex lines) and the absence of viscosity make superfluids remarkably different from ordinary fluids. Nevertheless, recent studies [1] have revealed that superfluid helium, when suitably stirred, shares an important property with ordinary turbulence: the same Kolmogorov energy spectrum [2], describing a distribution of kinetic energy over the length scales which signifies an energy cascade from large length scales to small length scales. This finding suggests that turbulence of quantized vortices (quantum turbulence) may represent the 'skeleton' of ordinary (classical) turbulence [3]. A puzzle arises however: experiments [4,5] show that there are other regimes in which turbulent superfluid helium lacks the Kolmogorov spectrum.
Trapped atomic Bose-Einstein condensates (BECs) are ideal systems to tackle this puzzle. This is because, unlike helium, the physical properties of BECs (e.g. the strength of atomic interactions, the density, the vortex core radius) can be controlled. Moreover, individual quantized vortices are more easily nucleated, manipulated [6,7] and observed [8][9][10] in BECs than in helium.
A second motivation behind our work is that the study of three-dimensional turbulence in BECs [11] is held back by the lack of a standard method to excite turbulence in a reproducible way, so that experiments and numerical simulations can be compared with each other and any generality of results can be more easily recognized. In classical turbulence, standard benchmarks are flows driven along channels or stirred by propellers, flows around welldefined obstacles (e.g. cylinders, spheres, steps) and wind tunnel flows past grids. Similar techniques are used for superfluid helium [12,13], which is mechanically or thermally driven along channels, or stirred by oscillating wires, grids, forks, propellers and spheres. Vortices and turbulence in BECs have been generated by moving a laser beam across the BEC [9,[14][15][16][17][18][19][20], by shaking [21] or stirring the trap rotating it around two perpendicular axes [22], by phase imprinting staggered vortices [23], or by thermally quenching the system (Kibble-Zurek mechanism) [24][25][26][27]. This variety of techniques and the arbitrarily chosen values of physical parameters means that comparisons are difficult. Moreover, the disadvantage of some of these techniques is that they tend to induce large surface oscillations or even fragmentation [28] of the condensate which complicate the interpretation of results and the comparison with classical turbulence.
The aim of this report is to propose a technique to induce turbulence in BECs based on the decay of multicharged vortices (a more controlled and less forceful technique than the above mentioned methods), and to characterize the turbulence which is produced. For this purpose, we shall compare two disordered vortex states, which we shall call 'anisotropic' and 'quasi-isotropic' for simplicity, resulting respectively from the decay of a single j = 4 charged vortex and from the decay of two antiparallel j = 2 vortices. We shall reveal a new connection between BEC turbulence and the so-called 'Vinen' turbulent regime discovered in superfluid helium.

II. MULTICHARGED VORTICES
In a superfluid, the circulation around a vortex line is an integer multiple (j = 1, 2, · · · ) of the quantum of circulation [29] where C is a closed path around the vortex line, h is Planck's constant and m the atomic mass. The velocity field around an isolated vortex line is therefore constrained to the form v = jκ/(2πr) where r is the distance to the vortex axis. This property is in marked contrast to ordinary fluids where the velocity of rotation about an axis (e.g. swirls, tornadoes, galaxies) has arbitrary strength and radial dependence.
The angular momentum and the energy of an isolated vortex in a homogeneous superfluid grow respectively with j and j 2 [29]. Therefore, for the same angular momentum, multicharged (j > 1) vortices carry more en-ergy, and, in the presence of thermal dissipative mechanisms, tend to decay into singly-charged vortices [30][31][32][33]. Besides the energy instability, there is also a dynamical instability [34][35][36] which would destabilize a multicharged vortex. The time-scale of these effects has been investigated [30,32,33,35,[37][38][39][40]. The technique of topological phase imprinting [41] has allowed the controlled generation of multicharged vortices [30,32,39,42] in atomic condensates. The decay of a doubly-quantized vortex into two singly-quantized vortices has been studied [30,37,38] in a Na BEC. Quadruply-charged quantized vortices are also theoretically predicted [43] to decay. Recent work has determined that the stability of such vortices is affected by the condensate's density [30] and size [38], and by the nature of the perturbations [43].

III. MODEL
We model the condensate's dynamics using the 3D Gross-Pitaevskii equation (GPE) [29] for a zerotemperature condensate where ψ(r, t) is the condensate's wavefunction, r the position, t the time, and the harmonic trapping potential. The parameter g = 4πh 2 a s /m characterizes the strength of the inter-atomic interactions, a s is the s-wave scattering length. The normalization is V |ψ| 2 dV = N where V is the BEC's volume and N the number of atoms. We cast the GPE in dimensionless form using τ HO = ω −1 r , HO = h/mω r and hω r as units of time, distance and energy respectively. The inter-atomic interaction parameter g is chosen to describe a typical BEC with N ≈ 1 × 10 5 atoms of 87 Rb trapped harmonically in a cigar-shaped BEC with radial and axial frequencies such that ω z /ω r = λ = 0.129. It is of our particular interest to study properties of the condensate's velocity field components, which are computed from the definition v(r) = (ψ∇ψ * − ψ * ∇ψ) /2i |ψ| 2 . The dimensionless GPE is solved numerically in the 3D domain −10 HO ≤ x, y ≤ 10 HO and −40 HO ≤ z ≤ 40 HO on a 128 × 128 × 512 grid (keeping the same spatial discretization in the three directions) with time-step ∆t = 10 −3 using the 4th order Runge-Kutta method with XMDS [44]. We have performed tests with different grid sizes and verified that our results are independent of the discretization.

IV. DECAY OF SINGLE QUADRUPLY-CHARGED VORTEX
The shape of the singly-charged vortex lines emerging from the decay of a multicharged vortex depends on where and when the decay starts. The singly-charged lines may be straight or intertwined (as reported here) depending on the perturbation's symmetry and the local density homogeneity. If the perturbation is uniform and the density does not vary much in the z-direction, every point on the vortex unwinds at the same rate, and straight singly-charged vortex lines will emerge. However, if the density changes significantly along z, the unwinding takes place at different times at different positions, inducing intertwining, as discussed in [37]. The twisted vortex decay is shown in Fig. 1 and in the first movie in the Supplementary Material [45]. The main feature which is visible during the evolution are Kelvin waves. Kelvin waves consist of helical displacements of the vortex core axis, and play an important role in quantum turbulence [12]; they have been recently experimentally identified in superfluid helium [46] and their presence has been recognized in atomic condensates [47,48]. Considering Fig. 1, it is worth distinguishing the Kelvin waves which, in our case, emerge on parallel vortices from the decay of a multicharged vortex [35,38] in a confined geometry, from the Kelvin waves generated by the Crow instability [49] on anti-parallel vortices in a homogeneous condensate.
Our numerical experiments suggest that the decay of the multicharged vortex can be sped up. Imposing random fluctuations (≤ 10% of |ψ|) to the initial j = 4 wave function does not significantly change the decay time scale, probably because the symmetry of the initial condition is not completely broken. A small displacement of the vortex core axis (≈ 0.5 HO ) is more efficient, triggering the onset of the twisted unwinding in about 2.4 τ HO ; a larger displacement (≈ 0.2 HO ) reduces this time to 2.0 τ HO . Among the other methods which we have investigated, the most efficient is to gently squeeze the harmonic potential in the xy plane by an amount ω x /ω y = 0.9 when preparing the initial state in imaginary time, then resetting ω x /ω y = 1 when propagating the GPE in real time; the squeeze triggers the onset of decay in only 1.0 τ HO . In the experiments, it is usually difficult to control perturbations well enough to reproducibly determine the time scale of decay.

V. DECAY OF TWO ANTIPARALLEL DOUBLY-CHARGED VORTICES
Now we exploit the twisted unwinding of a multicharged vortex as a convenient technique to generate turbulent vortex tangles which are relatively free of large density perturbations. We start by numerically imprinting anti-parallel, doubly-charged vortices as the initial state. One vortex is centered at position (x, y) = (1.8 HO , 1.5 HO ) and the other at (x, y) = (1.0 HO , −1.3 HO ), as shown in Fig. 4(a). During the evolution, the vortices unwind, twist, move slightly forward due to the self-induced velocity field, and then reconnect, generating a turbulent state with only moderate density oscillations, as shown in Fig. 4(c).
We analyze the turbulent state in terms of the spatial averages of the velocity components. At the beginning of the decay, we find |v y | / |v x | ≈ 1 and (as expected, as vortices are initially aligned in the z-direction) |v z | / |v x | ≈ 0. Fig. 5 (top) shows that, as the time proceeds, the vortex configuration becomes almost isotropic, indeed at t = 10 τ HO we have |v y | / |v x | ≈ 1.00 and |v z | / |v x | ≈ 0.77. For simplicity of reference, we call a 'quasi-isotropic state' the resulting disordered vortex configuration.
The oscillations of the transverse velocity components v x and v y appear practically in phase for the anisotropic case (see Fig. 2 (bottom)) and almost completely out of phase for the quasi-isotropic case (see Fig. 5 (bottom)), and these suggest the existence of collective modes. In order to properly identify these modes we have evaluated the time evolution of the condensate's transverse widths w x and w y (found by adjusting Gaussian fits on x-and ydirections over the z-integrated density) for both cases. As can be seen in Fig. 3, the anisotropic case exhibits a quadrupolar mode, in which w x and w y oscillate out of phase in time. After performing a Fourier analysis of these quadrupolar oscillations we verified the mode's frequency to be ω ∼ √ 2ω r , in agreement with theoretical predictions for a vortex-free cloud [50]. This mode is reminiscent of the unstable (quadrupolar) Bogoliubov mode that drives the decay of the initial multicharged vortex, as identified for an analogous 2D trapped system in [40]. However, in the quasi-isotropic case, Fig. 6 shows that after t ≈ 10 τ HO there is a small (∼ 0.2 HO ) in-phase oscillation of the widths (as opposed to the larger oscillations for the anisotropic case in Fig. 3). The latter can be identified as a breathing mode and (again, through a Fourier analysis) it was found to exhibit a characteristic frequency of ω ∼ 2ω r , also in agreement with theoretical predictions for a vortex-free cloud [50]; this value is expected to hold even for rapidly rotating trapped systems [51]. These results, alongside a visual comparison of Figs. 1 and 4, show that the outer surface of the condensate is actually slightly less disturbed in the quasiisotropic case. Our turbulent condensate (Fig. 3(c)) is clearly less 'wobbly' than condensates made turbulent via other stirring methods, notably Refs. [23] and [28]).
We proceed and analyze the distribution of values of the turbulent velocity components (Fig. 7). We find that the probability density functions (PDFs, or nor- malized histograms) display the typical power-law scal- where α x ≈ −2.97, α y ≈ −2.95, α z ≈ −3.20, in agreement with findings in larger condensates [23]. Such power-law scaling, characteristic of quantum turbulence and observed in helium experiments [52], contrasts the Gaussian PDFs which are typical of ordinary turbulence. The difference between power-law and Gaussian statistics is important at

VI. IDENTIFICATION OF THE TURBULENCE
The two disordered vortex states, anisotropic and quasi-isotropic, produced respectively by the decay of a single quadruply-charged vortex (Section IV) and by decay of two antiparallel doubly-charged vortices (Section V) are clearly different. In the anisotropic state, all vortex lines are aligned in the same direction, and the net nonzero angular momentum constrains the flow. In the quasi-isotropic state, the oscillations of the transverse velocity components are three times larger, suggesting the presence of high-velocity events (vortex reconnections between Kelvin waves growing on opposite-oriented vortices) which are the hallmarks of turbulence. Moreover, in the quasi-isotropic state, the zero angular momentum of the initial configuration allows a redistribution of the velocity field, making the amplitudes of the three velocity components almost equal; indeed, after the initial vortex has split (t ≈ 10 τ HO ), the axial velocity |v z | ≈ 0.9 HO /τ HO , is not much smaller than the transverse velocity |v x | ≈ 1.2 HO /τ HO and |v y | ≈ 1.1 HO /τ HO . On the contrary, in the anisotropic state, at the same stage (t ≈ 10 τ HO ), the axial velocity component is much smaller than the transverse components. In other words, the velocity field which results from the decay of the antiparallel doubly-charged vortex state is indeed almost isotropic.
Since there is not yet a precise definition of turbulence, in principle both disordered states investigated here could be considered somewhat 'turbulent'. However at this early stage of investigation, we want to make conceptual connections to the simple isotropic cases known in the literature (in particular, the Kolmogorov and Vinen types of turbulence). Therefore hereafter we concentrate on the quasi-isotropic state.
The next question is whether the quasi-isotropic state is 'turbulent' in the sense of ordinary turbulence or in comparison with turbulent superfluid helium. To find the answer, we turn to the workhorse of statistical physics: the correlation function.
The quantity which measures the degree of randomness of a turbulent flow is the (normalized) longitudinal correlation functions [54,55], defined by where the symbol · · · denotes average over x, andê i is the unit vector in the direction i = x, y, z.
In fluid dynamics, c represents the size of the large eddies. In our case, c is the length scale over which the velocity field is highly correlated. Fig. 8 shows that the correlation functions drop to only ≈ 10% at distances of the order of the average separation between the vortex lines, ≈ 6.1 HO ; the last quantity is estimated as ≈ L −1/2 from the measurement of the vortex line density L (the vortex length per unit volume). Physically, this lack of correlation means that the vortex lines are randomly oriented with respect to each other. The analysis of the correlation function therefore suggests that the turbulent velocity field arising from the decay of antiparallel doubly-charged vortices is essentially a random flow.
This result implies that the distribution of the kinetic energy over the length scales, or energy spectrum E(k) (where the wavenumber k represents the inverse length scale), should be very different from the celebrated Kolmogorov scaling, E(k) ∼ k −5/3 , which is observed in classical turbulence and implies a particular structure of the flow. The importance of the Kolmogorov scaling is that it is the signature of a nonlinear cascade mechanism which transfers energy from large length scales to small length scales. Besides ordinary turbulence, the classical Kolmogorov scaling has been observed in helium experiments [56,57] when the turbulence is generated by grids in wind tunnels or by counter-rotating propellers. Numerical simulations [58] revealed that the Kolmogorov scaling is associated with the presence of large-scale polarization (bundles) of vortex lines which become locally parallel to each other; the bundles locally create a net average rotation over length scales larger than , thus building up energy at wavenumbers k < k = 2π/ (the hydrodynamical range).
To compute the energy spectrum E(k) of our turbulent condensate, we use a standard procedure [59] to extract the incompressible kinetic energy from the total energy, obtaining Fig. 9. To interpret the figure, we mark with vertical lines the wavenumbers k ξ = 2π/ξ, k a = 2π/a, k = 2π/ , k d = 2π/d TF and k D = 2π/D TF corresponding to the healing length ξ, the vortex core radius a, the average vortex separation , and the radial and axial Thomas-Fermi radii, d TF and D TF . Fig. 9 clearly shows that, at the hydrodynamical length scales k < k the energy spectrum is not of Kolmogorov type: under the classical k −5/3 scaling, substantially more energy would be contained in the small k region. Instead, over a wide range of wavenumbers up to k a , we observe the E(k) ∼ k −1 spectrum which is characteristic of an isolated straight vortex line. This means that, at distances less than , the velocity field is dominated by the nearest vortex in the vicinity of the point of observation -the effects of all the other vortices in a random tangle canceling each other out. Fig. 9 is therefore consistent with the random flow interpretation which results from the analysis of the correlation functions. It is also worth noticing that the E(k) ∼ k −3 scaling which appears in the region k a < k < k ξ is characteristic of the vortex core [60,61].
Besides the correlation function and the energy spectrum, further insight into the nature of our turbulence is acquired by measuring the temporal decay of the vortex line density L(t). This decay is caused by sound radiated away by vortices as they accelerate about each other [62] or reconnect [63] with each other. Fig. 8 shows that, at large t, the decay is consistent with the form L(t) ∼ t −1 , as reported for a larger spherical condensate [23]; the decaying turbulence is shown in the second movie in the Supplementary Material [45]. The L(t) ∼ t −1 decay is characteristic of the 'Vinen' turbulence regime (also called 'ultra-quantum regime') which was generated in superfluid helium by Walmsley & Golov [4] by short injections of vortex rings; long injections generated a 'Kolmogorov' (or 'quasi-classical') regime which decayed as L(t) ∼ t −3/2 . Numerical simulations [64] of Walmsley & Golov's experiment revealed that (at large k within the hydrodynamical range k < k ) 'Vinen turbulence' (decaying as L(t) ∼ t −1 ) has spectrum E(k) ∼ k −1 and 'Kolmogorov turbulence' (decaying as L(t) ∼ t −3/2 ) has The distinction between 'Vinen turbulence' and 'Kolmogorov turbulence', first proposed by Volovik [65], seems to be characteristic of superfluids. It is now understood that Vinen's pioneering experiments [66] on counterflow heat transfer in superfluid helium generated 'Vinen turbulence': this was demonstrated by numerical models of counterflow turbulence [67] which produced the expected E(k) ∼ k −1 energy spectrum and the L(t) ∼ t −1 vortex line density decay, which is the larget decaying solution of the equation dL/dt ∼ −L 2 proposed by Vinen on simple physical arguments to model a random-like flow.
A recent work [54] has examined the properties of the turbulence following a thermal quench of a Bose gas. It found that topological defects created by the Kibble-Zurek mechanism evolve into a turbulent vortex tangle [68] which eventually decays into a vortex free state. During the decay, which has the form L(t) ∼ t −1 , the energy spectrum is E(k) ∼ k −1 and the correlation functions drop to few percent over the distance . The thermal quench is therefore another clear example of 'Vinen turbulence'.

VII. CONCLUSION
In this work we have explored the decay of initially imprinted multicharged vortices as a method to generate turbulence in a trapped Bose-Einstein condensate which is relatively free from large surface oscillations and fragmentation. We have examined the decay of two multicharged vortex systems, in a typical cigar-shaped harmonically confined atomic Bose-Einstein condensate.
The first (a quadruply charged vortex) led to the disordered, anisotropic vortex state. The second (two antiparallel doubly-charged vortices) generated helical Kelvin waves on opposite oriented vortex lines which reconnected, creating a second disordered, quasi-isotropic vortex state. Looking for similarities with ordinary turbulence in its simplest possible forms -in particular with the property of isotropy -we have concentrated the attention on the quasi-isotropic state, and carefully considered in which sense it is turbulent.
This question is subtle. In classical physics, turbulence implies a large range of length scales which are all excited and interact nonlinearly. In ordinary fluids, the presence and the intensity of turbulence is inferred from the Reynolds number Re, which must be sufficiently large (typically few thousands, depending on the problem) for turbulence to exist. But the Reynolds number has two definitions. The first definition is where D is the system's large length scale, i.e. the system's size or the size of the energy containing eddies, u is the flow's velocity at that large scale, and ν is the kinematic viscosity; this definition follows directly from the Navier-Stokes equation, and measures the ratio of the magnitudes of inertial and viscous forces acting on a fluid parcel. The definition makes apparent why large scale (e.g. geophysical) flows are always turbulent, and why microfluids flows are not (indeed, with microfluids one has to rely on chaos, not turbulence, for achieving any desired mixing). The second definition assumes Kolmogorov theory, and is (where η is the length scale of viscous dissipation). This definition measures the degree of separation between the large length scale (at which energy is typically injected) and the small length scale (at which energy is dissipated).
In the context of condensates, the two definitions clash with each other: The first definition implies that Re is infinite (because the viscosity is zero), and the second that Re is not much larger than unity (because the size of a typical condensate is larger, but not orders of magnitude larger, than the healing length, which can be considered the length scale at which acoustic dissipation of kinetic energy occurs). Although interesting work is in progress to identify a definition of Reynolds number suitable for a superfluid system (for example exploring dynamical similarities [69]), to answer the question which we asked, at this stage, we have to leave the Reynolds number and proceed in other ways. We have therefore carefully examined the properties of the disordered, quasi-isotropic state in terms of velocity statistics, energy spectrum, correlation function and temporal decay, and compared them to the properties of ordinary turbulence and of turbulent superfluid helium. Clearly, the quasi-isotropic state does not compare well with the properties of ordinary turbulence. Despite the limited range of length scales available in a small BEC, we conclude that, in the decay of the two antiparallel doubly-charged vortices, the quasi-isotropy of our disordered state, the short correlation length, the E(k) ∼ k −1 scaling of the energy spectrum at the hydrodynamical length scales, and the L(t) ∼ t −1 temporal behavior of the decay, identify our disordered, quasi-isotropic vortex state as an example of the 'Vinen turbulent regime' discovered in superfluid helium at low temperatures.
The nature of the disorder, or turbulence, in the anisotropic case generated by the decay of the single quadruply-charged vortex will be the subject of future investigations: one should vary the amount of polarization and study fluctuations of the velocity field over the mean rotating flow. There are no numerical studies yet of such turbulence in trapped Bose systems, and (because of the role played by boundaries in the spin-down of viscous flows) no immediate classical analogies, so this case is less straightforward to analyze; spin-down dynamics experiments in superfluid helium [70][71][72] should be the main reference systems.
Future work will also address the next natural question: can the classical Kolmogorov regime be achieved in much larger condensates under suitable forcing at the largest length scale, as suggested by some numerical simulations [22], thus identifying the cross over between Vinen and Kolmogorov turbulence?