Unitary Qubit Lattice Gas Representation of 2 D and 3 D Quantum Turbulence

Turbulence is of vital interest and importance to the study of fluid dynamics Pope (1990). In classical physics, turbulence was first studied carefully for incompressible flows whose evolution was given by the Navier-Stokes equations. One of the most celebrated results of incompressible classical turbulence (CT) is the existence of an inertial range with the cascade of kinetic energy from large to small spatial scales until one reaches scale lengths on the order of the dissipation wave length and the eddies/vortices are destroyed. The Kolmogorov kinetic energy spectrum in this inertial range follows the power law in wave number space.


Introduction
Turbulence is of vital interest and importance to the study of fluid dynamics Pope (1990).In classical physics, turbulence was first studied carefully for incompressible flows whose evolution was given by the Navier-Stokes equations.One of the most celebrated results of incompressible classical turbulence (CT) is the existence of an inertial range with the cascade of kinetic energy from large to small spatial scales until one reaches scale lengths on the order of the dissipation wave length and the eddies/vortices are destroyed.The Kolmogorov kinetic energy spectrum in this inertial range follows the power law in wave number space.
(1) Independently, quantum turbulence (QT) was being studied in the low-temperature physics community on superfluid 4 He Pethick & Smith (2009); Tsubota (2008).However, as this QT dealt with a two-component fluid (an inviscid superfluid and a viscous normal fluid interacting with each other) it considered phenomena not present in CT and so no direct correspondence could be made.With the onset of experiments in the Bose Einstein condensation (BEC) of dilute gases, we come to a many body wave function that at zero temperature reduces to a product of one-body distributions.The evolution of this one-body distribution function ϕ(x, t) is given by the Gross-Petaevskii (GP) equation Gross (1963); Pitaevskii (1961): with the nonlinear term |ϕ| 2 ϕ arising from the weak boson-boson interactions of the dilute BEC gas at temperature T = 0. Eq. ( 2) is ubiquitous in nonlinear physics: in plasma physics and astrophysics it appears as the envelope equation of the modulational instability while in nonlinear optics it is known as the Nonlinear Schrodinger (NLS) equation Kivshar & Agrawal (2003).
The quantum vortex is a topological singularity with the wave function |ϕ|→0 at the vortex core Pethick & Smith (2009).A 2π circumnavigation about the vortex core leads to an integer multiple of the fundamental circulation about the core: i.e., the circulation is quantized.Thus the quantum vortex is a more fundamental quantity than a classical vortex.In CT one has at any time instant a sea of classical vortices with continuously varying circulation and sizes, making it difficult to define what is a vortex.Thus there was the hope that QT (which consists of entangled quantum vortices Feynman (1955) , while of fundamental importance in its own right, might shed some light on CT Barenghi (2008) .Under the Madelung transformation the GP Eq. ( 2) can be rewritten in conservation fluid form (ρ the mean density and v the mean velocity) ). (5) Thus QT in a BEC gas will occur in a compressible inviscid quantum fluid whose pressure now has two major contributions: a barotropic pressure ∼ gρ 2 and a s-called quantum pressure ∼− √ ρ ∇ 2 √ ρ.
In the seminal paper on QT in a BEC , Nore et.al. Nore et al. (1997) stress the Hamiltonian nature of the GP system and performed 3D QT simulations on a 512 3 grid using Taylor-Green vortices as initial conditions.Nore et.al. Nore et al. (1997) find a transient glimmer of a k −5/3 incompressible energy spectrum in the low-k wave number range, but later in time this disappeared.Other groups Kobayashi & Tsubota (2007); Machida et al. (2010); Numasato et al. (2010); Tsubota (2008)] attacked QT in GP and just recently achieved a simulation on a 2048 3 grid using pseudo-spectral methods Machida et al. (2010) .However, most the algorithms of the Tsubota group introduce an ad hoc dissipative term to the GP equation.Ostensibly, this added dissipative term damps out the very large k-modes and a hyperviscosity effect is commonly used in (dissipative) Navier-Stokes turbulence to suppress numerical instabilities.However the introduction of dissipation into the GP Eq. ( 2) may be more severe since it destroys its Hamiltonian nature -while in CT this dissipation just augments the actual viscosity in the Navier-Stokes equation.There are a considerable number of studies of the effects of hyperviscosity on backscatter and bottlenecks in the kinetic energy spectrum.The presence of ad hoc dissipation in QT can have important effects on the energy backscatter from the large k to smaller k that is even seen in CT simulations [ref..].They Machida et al. (2010) find a pronounced small-k region that exhibited a Kolmogorov-like incompressible energy spectrum of k −5/3 , which after some time decayed away.We now briefly outline the sections in this chapter.In Sec. 2 the difference between 2D and 3D CT are briefly summarized.In Sec. 3 we introduce our novel unitary quantum lattice gas (QLG) algorithm to solve the GP equation.The untiary nature of our algorithm has several important features: (a) it completely respects the Hamiltonian structure of the T = 0 BEC dilute gas, and (b) it permits a determination of a class of initial conditions that exhibit very short Poincare recurrence times.This has not been seen before in other simulations of the GP equation because either this class of initial conditions has been missed or because the introduction of any ad hoc dissipative terms destroy the existence of any Poincare recurrence.
In Sec. 4 we discuss our QLG simulations on 2D QT using a variety of initial conditions and compare our findings to the recent papers of Horng et. al. Horng et al. (2009) and Numasato et. al. Numasato et al. (2010) In Sec. 5 QLG simulations of 3D QT are presented with particular emphasis on the energy spectra.Finally in Sec.6 we discuss future directions of QLG.

Coherence length, ξ
Before we summarize the differences between 2D and 3D CT we shall consider the concept of coherence length that is very commonly introduced in low-temperature physics.In particular we consider some of the ways the coherence length ξ is introduced for BECs and for QT.The first is a back-of-the-envelope definition of the coherence length Pethick & Smith (2009) : it is the length scale at which the kinetic energy is of the same order as the nonlinear interaction term in the GP Eq. ( 2).Approximating ∇∼ξ −1 and the nonlinear term by the asymptotic density |ϕ|∼ A more quantitatively derived expression for ξ for1D vortices is given by Pethick & Smith Pethick & Smith (2009) who readily show that the exact steady state solution to the GP Eq. ( 2) under the boundary conditions: ϕ(0)=0 and ϕ(x → ∞)= √ ρ 0 .The coherence length is thus the minimal distance from the singular core to the asymptotic value of the vortex wave function.A similar result, under the same boundary conditions, hold for 2D and 3D line vortices when one uses the Pade approximant expression of the radial part of the wave function following Berloff Berloff (2004) .Altenatively, Nore et.la.Nore et al. (1997) where < ρ 0 > is the spatially averaged BEC density.The concept of coherence length becomes much more subtle for strongly nonlinear flows and when the (initial) wave function in the simulation is not a quasi-solution of the GP equation.
For example, we shall consider random initial conditions or rescaled wave functions which thus are no longer solutions to the GP equation because the GP equation is nonlinear.For such problems the coherence length is at best a qualitative concept with possibly ξ = ξ(t).

2D CT
In 2D incompressible CT it is convenient to introduce the scalar vorticity and write the 2D Navier-Stokes equation in the form In the inviscid limit, ν → 0, there are two constants of the motion -the energy E and the enstrophy Z Under the assumption of incompressibility, isotropy and self-similarity, Batchelor (1969); Kraichman (1967) have argued for the existence of a dual cascade in the inertial range: an inverse cascade of (incompressible) kinetic energy to small k while a direct cascade of enstrophy to large k.From dimensional arguments they then predicted the spectral exponents of the energy spectrum E(k), where E = dkE(k) : from the energy inverse cascade to small k: E(k) ∝ k −5/3 ; (12) from the enstrophy direct cascade to large k: E(k) ∝ k −3 .( 13)

3D CT
In 3D incompressible CT there is only one inviscid constant of the motion -the total energy.This leads to a direct cascade of energy from large to small k with the well-known Kolmogorov spectrum When one considers compressible classical turbulence in real flows one is typically forced into some subgrid large eddy simulation (LES) modeling.It is interesting to note that in these models one typically requires closure schemes that will produce a Kolmogorov k −5/3 spectrum for the subgrid total kinetic energy Menon & Genin (2010) 3. Unitary quantum lattice gas algorithm for the Gross-Pitaevskii equation

The unitary algorithm
We briefly outline the unitary quantum lattice gas (QLG) algorithm for the solution of the GP equation.Instead of employing a direct scheme to solve the GP equation we move to a more fundamental mesoscopic level and introduce qubits on a lattice.In particular, to recover the scalar GP equation we need to just introduce 2 qubits at each lattice site.A classical bits can take on only one of 2 possible values, which shall be denoted by |0 or |1 .However a qubit can take on an arbitrary superposition of these two states: |q = α q |0 + β q |1 .To take advantage of quantum entanglement we need to introduce at least 2 qubits per lattice site if the unitary collision operator is restricted to entangle only on-site qubits.To recover the GP Eq. ( 2), it will turn out that we will need to consider just 2 qubits per lattice site and of the basis set of the 4 posisble states |00 , |01 , |10 and |11 , only the subset |01 , |10 are required.We introduce the complex probability amplitudes α and β for these states |01 , |10 .Thus at each cubic lattice position x we introduce the (complex) two-spinor field and construct the evolution operator U[Ω] -consisting of an appropriate sequence of non-commuting unitary collision and streaming operators -so that in the continuum limit the two spinor equation will reduce to the GP equation for the 1-particle boson wave function ϕ under the projection Consider the unitary collision operator C that locally entangles the complex amplitudes α and β

242
Advanced Fluid Dynamics www.intechopen.com where the σ are the Pauli spin matrices Now C 2 turns out to be just the swap operator because of its action on the amplitudes α and β : Hence C is typically known as the square-root-of-swap gate.
The streaming operators shift just one of these amplitudes at x to a neighboring lattice point at x + Δx: The subscript γ = 0 on the streaming operator S Δx,γ refers to shifting just the amplitude α while the subscript γ = 1 refers to just shifting the amplitude β.These streaming operators can be expressed in an explicit unitary form by using the Paul spin matrices where n =(1 − σ z )/2, n =(1 + σ z )/2.Note that the collision and streaming operators do not commute: [C, S] = 0. Now consider the following interleaved sequence of unitary collision and streaming operators Since |Δx|≪1 and C 4 = I, Eq.( 23) yields J 2 xγ = I + O(Δx), where I is the identity operator.This is because the streaming operators are O(Δx) deviations from the identity operator I.W e first consider the effect of the evolution operator acting on the γ component of the 2-spinor ψ.Here ε ≪ 1 is a perturbative parameter and Ω is a function to be specified later.
Using perturbation theory, it can be shown that the time advancement of ψ Yepez et al. (2009;2010) yields the spinor equation with γ = 0or1andΔx = O(ε).Since the order ε 3 term in Eq. ( 26) changes sign with γ, one can eliminate this term by introducing the symmetrized evolution operator rather than just a single U γ operator.Under diffusion ordering, Δt = O(ε 2 )=Δx 2 , the evolution equation leads to the spinor equation where the function Ω is still arbitrary.To recover the scalar GP Eq.2), one simply rescales the spatial grid ∇→a −1 ∇, contracts the 2-component spinor field ψ to the (scalar) BEC wave function ϕ ϕ =(1, 1) and chooses Ω = g|ϕ| 2 − 1: Thus the QLG algorithm that recovers the scalar GP Eq,(31) in the mesoscopic 2-spinor representation Eq.( 28) with the evolution operator Eq.( 27) and the component unitary operators defined by Eq.( 24).It is important to realize that there is nothing per se in the QLG that enforces diffusion ordering -which is critical for QLG to recover the GP Eq.( 2).This diffusion ordering must be recovered numerically by appropriate choices of available parameters and initial amplitudes.Since the QLG is unitary, the norm of the spinor ψ is automatically conserved.However, one finds (small) fluctuations in the mean density due to the overlap of the components of the 2-spinor If α is kept purely imaginary and β is kept purely real (or vice versa), the overlap between the two components vanishes and the mean density is conserved exactly.This is achieved by introducing two modifications to the QLA algorithm described above : • initialize the 2-spinor so that α = Re[ϕ], β = iIm[ϕ]; • replace the scalar potential function Ω by the unitary non-diagonal matrix: Now QLG, like its somewhat distant mesoscopic cousin the lattice Boltzmann algorithm Succi (2001), are explicit second order accurate schemes in both space and time.We have found that QLG and lattice Boltzmann have similar numerical accuracy in 1D soliton simulations.Moreover, it has been shown that lattice Boltzmann approaches the accuracy of spectral methods Succi (2001).This has been attributed to the extremely small coefficient multiplying the second order error term.All mesoscopic algorithms have to be carefully benchmarked -particularly as there can be some scaling requirements in order that the more general mesoscopic algorithm actually simulates the physics of interest.In particular, we have benchmarked our 1D QLG code against exact (theoretical) soliton solutions of the nonlinear Schrodinger equation -both scalar Vahala et al. (2003) and vector Vahala et al. (2004) soliton collisions -of nonlinear optics.

Parallelization performance on various architectures
One important property of QLG is that at this mesoscopic level the collision operator acts only on local data stored at the lattice site, while the streaming operator moves the post collision data to the nearest neighbor cell.This leads to ideal parallelization on supercomputers.Its unitary properties will permit direct encoding onto quantum computers.The scalings we report here are for recent codes that solve the coupled set of GP equations (BEC2) or for the T > 0 system which involves a coupling between the BEC ground state and the Bogoliubov modes (BdG).While the BEC2 code is an immediate vector generalization of the QLG of Sec. 2, the BdG code is a significant variation.Yepez has developed a 2-qubit QLG that now involves all 4-qubit states with 4 × 4 coupling matrices rather than the 2 × 2 coupling matrices in BEC2.The physics from these codes will be discussed elsewhere.In strong scaling, one considers the code's performance on a fixed grid.One increases the number of processors and checks the wallclock time needed to do the same fixed number of time steps as before.Ideal scaling occurs if the wallclock time decreases by 0.5 if the number of processors are doubled.We see that the strong scaling of the BEC2-code is excellent.The  dashed curve is the ideal scaling while the solid curve is our QLA code scaling.On the IBM BlueGene/P Intrepid this has been tested to 131072 cores and on larger grids (6400 3 ) it shows superlinear scaling (Fig. 1).Similar results are achieved on the CRAY XT5Jaguarpf, tested to 110592 cores and grids 6000 3 .The BdG-code, with its 4 × 4 entangling matrices on the 2 qubits, shows quite strong superlinear scaling on both the CRAY XT5 and XE6, testing to 216000 cores on Jaguarpf and to 150000 cores on HopperII (Fig. 2).
In weak scaling, one keeps the work done by each processor fixed.Thus if one doubles the grid in each direction, then one would need to increase the number of processors by 2 3 in a 3D simulations.For ideal scaling the wallclock time to completion should remain invariant.The weak scaling of these codes are given in TablesI-IV ,andagain are excellent, with fluctuations in timings typically below 5% as one moves from 216 cores to 110592 cores.

2D QT
We first consider 2D quantum turbulence.From the GP Eq.2, the total energy in 2D QT is conserved Since the GP equation also conserves the number of particles, the last term in Eq.( 32) yields an ignorable constant.The other terms in Eq.( 32) can be categorized as Nore et al. (1997) kinetic energy: To examine the effect of compressibility on the GP Eq.( 2), E K can be further decomposed into its incompressible and compressible components via Helmholtz decomposition Nore et al. (1997).Because a quantum vortex is a topological singularity, one needs to regularize the using definitions of velocity and vorticity.In particular, one defines a density weighted velocity q = √ ρv and its Fourier transform to be q.This differs from the standard Favre averaging used in standard compressible computational fluid dynamics (CFD) in that in QT one uses √ ρ Unitary Qubit Lattice Gas Representation of 2D and 3D Quantum Turbulence www.intechopen.comrather than ρ.Unlike CFD, in QT at the vortex core the density ρ → 0 with √ ρv → const..The longitudinal component and the transverse component of q are Correspondingly the compressible energy and incompressible energy can be defined as compressible : Consequently the energy densities of compressible and incompressible energy become using polar coordinates k, θ.For 2D QT it is useful to introduce a renormalized vorticity Horng et al. ( 2009) The two components of ω q are: Consequently, the time evolution of the density weighted enstrophy, Z q = dx 2 |ω q | 2 ,i sa n excellent measure of the time variation in the total number of vortices present in the 2D QT.
Incompressibility and the conservation of enstrophy are crucial elements for the existence of the dual cascade in 2D CT.However, in 2D QT, one has neither incompressibility nor enstrophy conservation.Bogoliubov elementary excitations, permitted by the compressibility of the quantum fluid, can create quantum vortices.On the other hand, counter-rotating vortex pairs can annihilate each other (in CT one has the merging of like-rotating vortices).Lacking the conservation of enstrophy due to compressibility, 2D QT does not necessitate dual cascade.
Another important feature of the GP Eq.( 2) is that it is Hamiltonian with total energy conserved.Thus we are guaranteed that after a sufficiently long time, the system will return to a state that is very close to its initial state.For nearly all continuous Hamiltonian systems, this recurrence time T P is effectively infinite.However Tracy et al. (1984) has demonstrated that in the 1D NLS (i..e the GP) system the Poincaré recurrence time can be unexpectedly short for certain initial manifolds.Arnold's cat map is an invertible chaotic map of a torus onto itself which is ergodic, mixing and structurally stable and it exhibits a short Poincaré recurrence time.To visualize such a recurrence under the Arnold cat map, we discretize a sample 2D square picture into m × m segments.Each segment is indexed with its 'coordinate' [q t , p t ], with t = 1....m.The Arnold cat map is applied to these segments: These segments are then patched together to form a new picture according to their updated coordinate [q t+1 , p t+1 ].We now consider the effect of pixel resolution on the Arnold cat map by considering m = 74 and m = 300.For m = 74 we find Poincaré recurrence after T P = 114 iterations, Fig. 3, while for the higher resolution m = 300, T P = 300 iterations, Fig. 4. What is quite interesting is the image after t = T P /2 iterations: for the low resolution run, m = 74 there is a point inversion, Fig. 3 (b) , while there is no point inversion for the high pixel resolution m = 300, Fig. 4 (b).However, the Arnold cat map exhibits quite complex characteristics with pixel resolution.For example, at resolution m = 307, one does see a point inversion of the image at t = T P /2 = 22.In our simulations of 2D QT we encounter unexpected short Poincaré recurrence time provided that the ratio between the time-averaged internal energy and averaged kinetic energy is sufficiently small: The plots of amplitude, vorticity and phase in the following sections will adopt the "thermal" scheme: blue stands for low values while red stands for high values.
is the approximated wave function of a single 2D vortex with winding number n. Vortices with opposite charge, e.g.n = 1 and n = −1, are interleaved to form a vortex array that satisfies the periodic boundary conditions.Fig. 5 illustrates the amplitude, density weighted enstrophy and phase of the initial wave function.The vortices can be identified by two methods: a) peaks in the vorticity plots or b) the end points of the branch cuts in the phase plots.For example, in Fig. 5(c), vortices with positive charge (phase varying from −π to π counter-clockwise around the singularity) are encircled in black while those with negative charge (phase varying from −π to π clockwise around the singularity) are encircled in white.In this simulation, grid size is 512 and the total number of iteration Vortices with winding number n = 2 are energetically unstable and will rapidly split into two n = 1 vortices.The Poincaré recurrence for winding number n = 2 case should be viewed as reproducting the state immediately following the vortex splitting.In this simulation, the energy ratio γ = 0.0036, which is much smaller than for the winding number n = 1 case due to the rotation induced by more vortices.The recurrence structure is most clearly visible via the evolution of the internal energy, Fig. 8.With more vortices present in the system, there is more energy exchange between vortices and sound waves.Therefore when the Poincaré recurrence occurs, more background noise is to be expected, and is seen in Fig. 9.One interesting observation in this simulation is that Kelvin-Helmholtz instability is important for vortices generation as suggested in Blaauwgeers et al. (2002); Henn et al. (2009).
When counter-rotating vortices approaches each other, in the regions between the vortex cores the velocity gradient can reach a certain critical value and trigger the Kelvin-Helmholtz instability.This instability will create pairs of counter-rotating vortices.As the number (f) ω q ,t=0 (g) t=10500 (h) t=21000 (i) =31300 (j) t=41900 Fig. 9. Poincaré recurrence for winding numbern =2.At the semi-Poincaré time, the wave function approximates the initial condition with shift in the origins.At t = 41900, the wave function bears a similar structure as at t = 0, but with more noise.
of vortices increases in the system, the probability for Kelvin-Helmholtz triggered vortex generation increases.Fig. 10 demonstrates such a process.

Random phase initial condition
In this simulation, the initial wave function features a constant amplitude and random phase ψ(r)=he iθ(r) .Bicubic interpolation is adopted to produce random phase which satisfies periodic boundary condition, c.f. Keys (1981).Under this interpolation, a 2D function x i y j .The coefficients a i,j are determined by the enforced continuity at the corners.Since there are 16 unknown coefficients, 16 equations are needed to determine these a i,j .Usually one can enforce continuity for p(x, y), ∂ x p(x, y), ∂ y p(x, y), ∂ x,y p(x, y).To generate random phase in the domain L 2 , the following procedure is followed: • Discretize the domain into m × m squares.m here is the level of randomness; • Generate 4 pseudo-random numbers at the 4 corners of each sub-square.These random numbers correspond to p, ∂ x p, ∂ y p, ∂ x,y p at the corners of this sub-square; • Enforce periodicity at the domain boundaries; • Solve a i,j belonging to each single sub-square.
For a grid of 512 2 we consider a randomness level m = 8.The random phase wave function is dynamically unstable.Sound waves are immediately emitted and create quantum vortices.Typically these vortices will decay away and the GP system tends to a thermal equilibrium, as demonstrated by Numasato et al. (2010).However, if the initial condition is chosen such that energy ratio γ ≪ 1, a Poincaré recurrence emerges.In our simulation, γ = 0.00287 and number of iteration is 100000.From the energy evolution plot, Fig. 11, the Poincaré recurrence time can be clearly identified by the abrupt energy exchanges (i.e., the spikes).
The first spike in Fig. 11 appears around t = 21000, with the second at t = 41900.These are just the semi-Poincaré and Poincaré times.One thus expects the phase distribution of the wave function at t = 0, t = 10500, t = 21000 and t = 41900 to illustrate the recurrence, c.f. Fig. 12.What is remarkable is that the randomly distributed vortices suddenly disappear from the system at T P /2 and T P .A tt = 41900, the phase distribution is very close to the initial state despite a constant shift in the central region.In Fig. 13 we have plotted the density-weighted vorticity at various times: there are no vortices at t = 0oratT P /2, (e), although a considerable amount of sound waves.
As energy ratio γ increases, the strength of the Poincaré recurrence is weakened by noise.Fig. 14 demonstrate how the Poincaré recurrence is lost as γ increases.When γ = 0.0567, one still can observe the depletion of vortices from the system at T P /2 = 21000, however, at T P = 41900, the initial condition, which is vortex free, can not be reproduced.For γ = 0.133, Fig. 12. Poincaré recurrence with random phase initial condition.At t = 10500, many randomly distributed vortices can be identified via the branch cuts.At t = 21000 and t = 41900, no vortices exist in the system since there are no branch cuts.There is an induced phase shift seen in the color scheme of the phase plots at t = 0 and t = 41900, but the geometric patterns are the same.
no trace of a short Poincaré recurrence can now be found.It needs to be pointed out that in our simulations, the Poincaré recurrence is characterized by abrupt energy exchange E K and E Q as well as among the compressible and incompressible components of the kinetic energy, E C and E IC .Therefore such phenomena can not be analyzed via standard turbulence theories invoking things like inverse cascades....

Energy cascade in 2D QT
For the simulations with vortices initially embedded in an Guassian BEC background, a k −3 power law is found ubiquitously in the compressible, incompressible and quantum energy spectra whenever vortices are present in the system.Fig. 15 describes the time evolution of the incompressible energy spectrum k s IC .The linear regression fit is over k ∈ [50, 100].Simulation grid 512 2 .In the time interval around t ∼ 24500 for winding number n = 1 embedded vortices, we examine the sudden drop in the spectral exponent s IC .In Fig. 16, the linear regression fit for the incompressible kinetic energy spectrum is made over the wave number interval k ∈ [50, 100].At t = 24400 and t = 24600, vortices are present in the system and the spectrum exponent s IC ∼−3.At t = 24500, when all the vortices are depleted, the incompressible kinetic energy spectral exponent decreases to s IC = −5.828.This could well indicate that the existence in     On this 32768 2 grid we performed a relatively short run to t max = 15000.The evolution of the quantum, internal and total energies are shown in Fig. 17.Based on the time evolution of the energies, the dynamics can be broadly categorized into two stages: (I) generation of vortices and (II) decay of vortices.In stage (I), the compressible energy decreases rapidly while the incompressible energy increases rapidly.Thus a significant amount of energy in the sound waves is transformed into incompressible energy induced by the rotational motion of vortices.In stage (II), the randomly distributed vortices disappear.The energy of the vortices is transferred into sound waves through vortex-vortex annihilation.Note that in this stage, the only major energy exchange occurs between the incompressible and compressible energies while the quantum and internal energies remain almost constants.The spectra for incompressible and compressible energy at t = 8000 is given in Fig. 18.At large-k region (k > 3000), a k −3 power law is present which can be interpreted as result of FFT of quantum vortices.It is interesting to notice that at k ∼ k ξ , e.g.k ∈ [700, 1200],ak −4 law is observed.We sampled the incompressible energy spectra every 50 iteration steps between 6000 < t < 10000 within a wave number window k ∈ [800, 1200].The time averaged slope s IC = −4.145± 0.066.This is in good agreement with the results obtained in Horng et al. (2009).This k −4 power law can be interpreted as the result of dissipation of randomly distributed vortices, as suggested in Horng et al. (2009).However, in low-k region (region I and II in Fig.18(b) and Fig. 18(c)) where semi-classical Kolmogorov cascade is expected, we did not observe the k −5/3 power law.This could be attributed to the compressibility of quantum fluid.Finally, we consider the case of 12 vortices of winding number n = 6 on an 8192 2 grid.The wave function is rescaled so that it is not a quasi-eigenstate of the GP Eq. ( 2).This will lead more quickly to turbulence.The compressible kinetic energy spectrum is shown in Fig. 19(a) and clearly exhibits a triple cascade: a small-k region with spectral exponent s C = −1.83,an intermediate-k region with exponent s C = −8.1, and a large-k region with exponent s C = −2.85.We will comment on this triple cascade spectrum more when we discuss QT in 3D.In Fig. 19(b) we plot the incompressible kinetic energy spectrum and notice the dual cascade region: for high-k we see the ubiquitous exponent s IC = −2.93 while for lower-k, the exponent becomes similar to the Saffman exponent, s IC = −4.0.The k ξ at which we have this dual spectrum meet is around the join of the steep intermediate range spectrum with the −3 spectral tail.We shall see this also in 3D QT.

3D QT
The major differences between 2D and 3D CT is in the behavior of the vorticity vector.In 2D, the vorticity is always perpendicular to the plane of motion while in 3D the vorticity vector can have arbitrary orientation.Here, in 3D we will employ variants of a set of linear vortices following the Pade approximant methods of Berloff Berloff (2004).For winding number n = 1 , using cylindrical polar coordinates (r, φ, z), a linear vortex that lies along the z-axis (and centered at the origin) is given by ϕ(r)=g −1/2 e iφ 11ar 2 (12 + ar 2 ) 384 + ar 2 (128 with |ϕ|→1/ √ g, and |ϕ 0 |→1a sr→ ∞, and |ϕ|∼r a/g as r → 0. Eq.( 45) is an asymptotic solution of the GP Eq.( 2).For this isolated linear vortex, the coherence length, from Eq. ( 6), ξ ∼ 4/ √ a.This is one of the reasons for introducing the factor into the GP Eq.( 2): a small a permits resolution of the vortex core in the simulations.If one starts with a periodic set of well-spaced non-overlapping Pade asymptotic vortices (clearly, of course, this will be dependent on the choice of the parameter a and the grid size L of the lattice) an asymptotic solution of the GP Eq.( 2) is simply a product of the shifted ϕ 0 's , weighted by g −1/2 .The system will evolve slowly into turbulence because this initial state is very weakly unstable.For these wave functions the coherence length is initially fairly well defined.On the other hand, in most of our runs, we just simply rescaled the asymptotic basis vortex function ϕ → g σ ϕ, for some σ.Because the GP Eq.((2) is nonlinear, g σ ϕ is no longer an asymptotic solution and the definition of coherence length becomes fuzzy.In Fig. 20 we show a somewhat complex initial vortex core situation.The initial wave function has winding number n = 5 and the positions of the initial line vortices are chosen so that there is considerable overlap of the wave functions around the center of the domain.In this plot we show not only the phase information on the vortex core isosurfaces but also the phase information on the boundary walls.On the vortex isosuface at t = 0 one can distinguish the 5 periods around the core.On the boundary walls, the intersection of the cores with the walls gives the location of the 4 branch-point like topological singularities.Emanating from each of these singularities are 5 branch cuts because of the chosen winding number n = 5.These branch cuts then join the branch points.Because the n = 5 state is energetically unfavorable, the initial state rapidly decays into 5× winding number n = 1 vortices, Fig. 20(b).It is tempting to identify the wave structures on the vortices as quantum Kelvin waves.Sound waves can also be identified on the boundaries.Near the center of the lattice, where there was initially considerable overlap of the vortices, many vortex loops have now formed.

Poincare recurrence for certain classes of initial conditions
As in 2D, a class of initial conditions will also be found for which the Poincare recurrence time is very short.The definitions of the incompressible and compressible kinetic energy, the quantum energy and the internal energy are immediate generalizations of those given in .Phase information is given on both the vortices and the boundary walls.The winding number n = 5 is evident from both the 5-fold periodicity around the each vortex as well as the 5 branch cuts emanating from each branch point on the boundary.By t = 3000, (b), the 5-fold degeneracy is removed with what seem like quantum Kelvin waves on the n = 1 cores.Basic phase coding : φ = 0 in blue, φ = 2π in red.Grid 2048 3 2D.As in the 2D GP case, short Poincare recurrence will be found for initial conditions such that E int (0), E qu (0) << E kin (0) with E comp (0) << E incomp (0).These conditions are readily satisfied in 3D by considering localized quantum line vortices so that ρ ∼ const.with ∇ √ ρ ∼ 0 much of the lattice.The evolution of a set of 3 × 16 linear vortices in the 3 planes and are examined for winding numbers n = 1 and n = 2, Fig. 21.Phase information is shown on the boundary walls: φ = 0-blue, φ = 2π -red.First consider the case of winding number n = 1.At t = 0, the phase information on the straight line vortex cores are clearly identified on their intersection with the boundaries.The corresponding branch cuts join the 48 branch points.A snapshot of the vortex isosurfaces is shown at t = 84000 and shows strong vortex entanglement with many vortex loops -basically a snapshot of a quantum turbulence state.However by t = 115000 we see a point inversion of the Poincare recurrence of the initial line vortices at t = 0, as was also seen in 2D GP flows.The full Poincare recurrence occurs around t = 230000, (d).The kinks along the the vortex cores may be quantum Kelvin waves: since one outputs at discrete time intervals, one is not at the exact T P .The robustness of the Poincare recurrence time is further exhibited by considering the evolution of 48 quantum vortices with winding number n = 2.The two-fold degeneracy translates into a more complex phase on the boundary walls, as can be seen on comparing (e) with (a) in Fig. 21.At t = T P = 230000, one sees considerable very small scale vortex loops that have arisen from the splitting of the confluent degeneracy although the overall structure of the line vortices can be seen globally seen.The phase information on the boundaries are only a slight perturbation from those initially.The grid for these runs was 1200 3 .For winding number n = 2 we show the details of the isosurfaces at the semi-Poincare recurrence time t = 115000 with phase information on the boundaries, and from a slightly different perspective.Also we show a detailed zoomed-in isoruface plot of the vortex cores at the Poincare recurrence time t = 230000 with phase information on the vortices themselves, Fig. 22.As in 2D, the signature of the occurrence of the Poincare recurrence can be seen in the evolution of the kinetic and quantum energies as a function of time, Fig. 23.The total energy is very well conserved throughout this run, t max = 250000, by our unitary algorithm on a grid 1200 3 : E TOT = const.For the parameters chosen here the internal energy is negligible.Note that the peaks in the kinetic energy are well preserved for vortices with winding number n = 1:E kin (0) ≈ E kin (t = 115000) ≈ E kin (t = 230000) ≈•••.However, for line vortices with winding number n = 2, there is a gradual decay in the peak in the kinetic energy E kin (0) ≥ E kin (t = 115000) ≥ E kin (t = 230000) ≥••• .This also explains why the Poincare recurrence in the isosurfaces for winding number n = 2 is not as clean as for winding number n = 1.Also it can be seen that in the evolution of E kin the vortex motion is much more turbulent for winding number n = 2. Since the internal energy for these runs is so low, the quantum energy evolution is the complement of kinetic energy (E TOT = const.= E kin (t)+E qu (t)+E int (t)).It should be noted that the time evolution of E kin (t) and E qu (t) (and, of course, E int (t)) are determined directly from their definitions, Eq.( 33) and their sum then gives us the (conserved) total energy.We note the loss of the semi-Poincaré time as the pixel resolution of the Arnold cat is increased from 74 to 300 × 300 yet find the persistence of the semi-Poincaré time for grid resolution from 512 3 to 1200 3 .Presumably this is because the QLG algorithm strictly obeys diffusion ordering so that T P on a 512 3 grid occurs at T P = 41775 and at T P = 230000 on a 1200 3 grid.(Diffusion ordering would give t = 229477).There is no such physics scaling laws in the Arnold cat map.(typically α > 7) joins to the large k spectrum, we find the switch over in the incompressible spectral exponents.We have noticed this in basically every simulation we have performed (and grids up to 4096 3 ).For vortices with winding number n = 2, the kinetic energy spectra much cleaner, presumably because the vortex entanglements are stronger and hence the QT is stronger.For the incompressible spectrum E inc (k) one again sees two spectral cascade regions: for k > k ξ the spectral exponent is α ∼ 3.07 while for k < k ξ , the exponent α ∼ 3.90.For the compressible spectrum, we find three spectral energy cacades: for very low k (5 < k < 30) a Kolmogorov-like cascade with exponent α ∼ 1.67, with a steep spectra decay followed for k > k ξ a compressible kinetic energy with exponent α ∼ 3.28.(The total kinetic energy spectrum has the exponents α ∼ 1.64 for low k, and α ∼ 3.17 for large k).The crossover k ξ ∼ 70 − 90.Somewhat surprising, we find a time interval during which we loose the incompressible kinetic energy spectrum E inc (k) ∼ k −3 for winding number n = 1 vortices.In particular, in the time interval 81400 < t < 84300 -except for very brief transient reestablishment of the k −3 spectrum -we find spectra as shown in  there is a sharp drop in the incompressible energy spectrum for wave numbers k > 100, except for a very brief transient recovery around t ∼ 83000.There is also a sharp cutoff in the compressible spectrum for k > 500.Around these intermittencies, the incompressible kinetic energy spectrum also To investigate the cause of this intermittent loss of the incompressible k −3 spectrum, we then examined the vortex isosurfaces around this time interval, Fig. 27.One notices that the loss of the k −3 corresponds to the apparent loss of vortex loops, i.e., of vortices.This would be consistent with the assumption that the incompressible kinetic energy spectrum of k −3 in the very large-k regime is due to the Fourier transform of an isolated Nore et al. (1997).As the vortex loops are reestablished, so is the incompressible k −3 kinetic energy spectrum.An alternative but somewhat more speculative explanation rests on the assumption that the incompressible k −3 spectrum is due to the quantum Kelvin wave cascade on the quantum vortices.As the quantum vortex loop shrink topologically, the Kelvin waves are inhibited and hence the loss of the k −3 spectrum.Moreover, if one looks at the time evolution of the mean kinetic E kin (t) and quantum E qu (t) energies one notices that this loss of the vortex loops occurs around the t ∼ 82000 around which the E kin (t) has a secondary  peak.As there is another secondary peak in E kin (t) around t ∼ (82000 + T P /2) one expects another transient loss in the k −3 spectrum and in the vortex loops.This is indeed found to occur around 196400 < t < 199300.For winding number n = 2 vortices, we do not find such intermittent loss of vortex loops or any intermittent loss of the k −3 spectrum.These results are in agreement with those found earlier in 2D QT.Moreover, there is not a similar intermittent loss of vortex loops for 48 linear vortices with winding number n = 1 (c.f., Fig. 21(b))

Total kinetic energy spectrum for large grid simulations on 5760 3
We have performed simulations on 5760 3 grid using winding number n = 6 straight line vortices as initial conditions.By t = 40000 one obtains the following total kinetic and quantum energy spectra, Fig. 28.Due to an oversight, we did not retain data for the incompressible In Yepez et al. (2009), we tried to identify these 3 regions as the small-k classical Kolmogorov casacde, followed by a semi-classical intermediate-k band (with non-universal exponent α) which is then followed a quantum Kelvin wave cascade for the very large-k band.Objections were raised against this interpretation Krstulovic & Brachet (2010);L'vov & Nazarenko (2010) based on (a) the kinetic energy spectrum of a single isolated vortex is k −3 for all k (for a straight line vortex all the kinetic energy is incompressible); and (b) using the standard definition of the coherence length ξ for the parameters chosen in our simulations, ξ > 2000 -i.e., it is claimed that we are investigating the physics of very strong vortex-vortex core overlapping wave functions.We counter that the definition of ξ is based on a boundary value solution of the GP equation for an isolated single vortex under the condition that the wave function asymptotes to the background value as one moves away from the vortex core.Our simulations are with periodic boundary conditions and we have no pointwise convergence of our wave function to some nice smooth 'background' value.While it can be argued that one should simply replace the usual background density ρ 0 by its spatial average < ρ >, the definition now of ξ becomes qualitative and does not handle large fluctuations about < ρ >.It is clear that we cannot categorically claim that the ubiquitous k −3 spectrum for the large-k band is due to quantum Kelvin wave cascade on single vortices -especially as this k −3 spectrum is also seen in our 2D QT simulations and in 2D there are no quantum Kelvin waves since the vortex core is just a point singularity.It is possible, however, that with the co-existence of this triple cascade region in the kinetic energy spectrum, and with the small-k region exhibiting a quasi-classical Kolmogorov k −5/3 spectrum, that this large-k band k −3 spectrum could be dominated by the spectrum of a single vortex as we are now at very small scales in the problem.
We have investigated in some detail the 3D QT for winding number n = 2 on a 3076 3 grid.One again finds the triple cascade region in the compressible kinetic energy spectrum, Fig. 29, (red squares), as well as in the quantum energy spectrum (gold diamonds).However the

Conclusion
Here we have discussed a novel unitary qubit algorithm for a mesoscopic study of quantum turbulence in a BEC dilute gas.Since it requires just 2 qubits/lattice site with unitary collision 270 Advanced Fluid Dynamics www.intechopen.comoperators that entangle the qubit probabilities and unitary streaming operators that propagate this entanglement throughout the lattice, QLG has a small memory imprint.This permits production runs on relatively few processors: e.g., production runs on 5760 3 grids using just over 10000 cores.Using standard pseudo-spectral codes, the largest grids achieved so far have been 2048 3 by Machida et al. (2010) -and these codes required the ad hoc addition of dissipative terms to damp out high-k modes and for numerical stability.For BEC turbulence this destroys the Hamiltonian structure of the GP Eq. ( 2) and destroys any Poincare recurrence phenomena.In its current formulation, it is critical that parameters are so chosen that the mesoscopic QLG algorithm yields diffusion ordering at the macroscopic GP level.This does restrict the choice in the values of kinetic, quantum and internal energies.If parameter choices are made that violate the diffusion ordering then the QLG algorithm will not be simulating the GP Eq. ( 2).There are various tests for the validity of our QLG solution of the GP equation: the conservation of total energy, and the fact the the Poincare recurrence time scales as grid 2whether in 2D or 3D GP.This replaces the naive thought that the time for QLG phenomena would necessarily scale as grid D , where D is the dimensionality of the macroscopic problem.On the other hand, standard CFD algorithms have the spatial and time step independent which, of course, is quite beneficial.We have presented significant spectral results for both 2D and 3D QT -although their interpretation is not straightforward and much still needs to be done in this area.Much of the controversy surrounds the significance of the coherence length ξ and the k −3 spectrum in the high-k region.We believe there is much new physics occurring in our QT simulations even if the k −3 spectrum is attributed to the dominance of the spectrum of an isolated vortex: (a) there is clear evidence in both 2D and 3D QLG of a dual cascade in the incompressible kinetic energy spectrum, with a spectrum of k −4 followed by the k −3 .This k −4 spectrum has also been seen by Horng et al. (2009) and connection implied with the Saffman spectrum arising from vorticity discontinuities; (b) the triple cascade in the total kinetic energy spectrum with the small-k regime yielding a quasi-Kolmogorov k −5/3 spectrum.It is a bit strange that the QT community is that concerned with achieving the k −5/3 energy spectrum in the incompressible energy E inc (k) in the small k-region where the quantization of the vortex core circulation becomes unimportant.The dynamics of the GP Eq.( 2) is fundamentally compressible.In compressible CFD simulations the emphasis changes to the total kinetic energy spectrum and that its power law is k −5/3 .However much work remains to be done on clarifying the role of quantum Kelvin waves on the energy spectrum.Finally, we comment on future directions of QLG.In this chapter we have restricted ourselves to the scalar GP equation.This is appropriate for a BEC gas with spin f confined in a magnetic well.The spin of the atom is aligned to the magnetic field and so the BEC dynamics is given by just one scalar GP equation.However if this BEC gas is confined in an optical lattice, the spin is no longer constrained and one now must work with 2 f + 1 GP equations.These so-called spinor BECs yield an enormous field of future research Ueda & Kawaguchi (2010).Moreover, the quantum vortices of a scalar BEC are necessarily Abelian vortices, but the vortices of spinor BECs can be non-Abelian in structure.Since QT is driven by vortex-vortex interactions, research needs to be performed to ascertain the role played by the non-Abelian in the energy cascades.(These non-Abelian vortices have non-integer multiples of the base circulation -not dissimilar to fractional quantum electron charge in the fractional quantum Hall effect).Other interesting vortices are skyrmions -used by high energy physicists to model baryons.

Fig. 1 .
Fig. 1.Strong scaling of BEC2 code on IBM BlueGene/Intrepid and CRAY XT5/Jaguarpf for various grids.The solid (blue) line is wallclock time for the BEC2 code while the dashed (red) line is ideal timing.100 iterations.

Fig. 2 .
Fig. 2. Strong scaling of the T > 0 BdG code on CRAY XT5/Jaguarpf and CRAY XE6/Hopper II.The solid (blue) line is wallclock time for the BdG code while the dashed (red) line is ideal timing.Some superlinear scaling is apparent.Grid 8400 3 , 100 iterations.
with major contributions coming from density variations near vortices and from sound and shock waves;•√ ρ∇×v • e z ∝ δ(r − r 0 ), which pinpoints the locations of the vortices.
Fig. 3. Poincaré recurrence under the Arnold cat mapping (I).Picture resolution: 74 × 74 pixels.Poincaré recurrence time T P = 114.Notice at the semi-Poincaré time t=57, the picture is almost a point inversion of the initial condition.

Fig. 4 .
Fig. 4. Poincaré recurrence under the Arnold cat mapping (II).Picture resolution: 300 × 300 pixels.Poincaré recurrence time T P = 300.At the semi Poincaré time, the picture shows no strong symmetry to the initial condition.

4. 1
Poincaré recurrence in 2D QT 4.1.1Four vortices with winding number 1 embedded in a Gaussian background BEC For the first set of simulations, we study the evolution of a set of vortices embedded in a Gaussian BEC background.Periodic boundary conditions are assumed.The initial wave function takes the form of products of individual quantum vortices, modulated by the inhomogeneous Gaussian background: Fig. 5. Four vortices with winding number 1 embedded in a Gaussian background.h = 0.05, a = 0.01, w g = 0.01, g = 5.0.Distance between vortices is L/4 with L being grid size 512.Since winding number is 1, the wave function's phase undergoes ±2π change around the singularities.steps is 50000.The ratio parameter γ = 0.018.A recurrence structure is clearly visible via the evolution of internal energy, E int (t), c.f. Fig.6.To further investigate such structure, we

Fig. 8 .
Fig.8.Evolution of internal energy for winding number 2. Around t ∼ 21000 and t ∼ 41900, internal energy reaches peaks.Comparing with winding number 1 case, more fluctuation is observed and the peaks are much narrower.This is caused by the high density of number of vortices and the frequent annihilation and creation of counter-rotating vortex pairs.

Fig. 11 .
Fig. 11.Energy evolution of the GP system with random phase initial condition.(a): evolution of energies (E K , E Q , E T ).Internal energy is negligible compared to the kinetic energy E K .(b) evolution of incompressible (red), E IC , and compressible kinetic energy (blue), E C .
Fig. 13.The evolution of the vorticity |ω(x, t))|.At a quantum vortex ω(x) ∼ ffi(x − x i ).At t = 0, the initial conditions for the wave function ϕ = √ ρ exp(iθ) are ρ = const and θ random.Thus there are no quantum vortices at t = 0. Very rapidly vortices are born.The vortices are annihilated at t = T P /2.Grid 512 2 .T P scales as L 2 , diffusion ordering, independent of whether 2D or 3D.

Fig. 14 .
Fig. 14.Loss of Poincaré recurrence.(a) evolution of E T , E K , E Q and E I , γ = 0.0567; (b) evolution of E T , E K , E Q and E I , γ = 0.133; (c) evolution of enstrophy Z q , γ = 0.0567; (d) evolution of enstrophy Z q , γ = 0.133.Depletion of vortices can be identified from the sharp decrease of enstrophy.
Fig. 15.Time evolution of the incompressible kinetic energy spectrum s IC .The red horizontal line indicates the k −3 power law.For winding number n = 1, there are spikes in the slope with s IC >> 3 in many instances.While for winding number n = 2 the variation in s IC is greatly reduced.
Fig.16.Phase plot and the incompressible energy spectrum round ∼ 24500.At t = 24400 and t = 24600 there are branch cuts and vortices in the BEC with spectral exponent s IC ∼−3.But at t = 24500, no branch cuts exist in the phase plot, indicative of no vortices in the system.The incompressible spectrum exhibits a discontinuity in the high-k region with a strong decrease in the exponent, s IC ∼−5.8.the incompressible kinetic energy spectrum of k −3 power law in the high-k region could be the by-product of the spectrum of a topological singularity -at least in 2D QT.It should be remembered that in 2D QT there can be no quantum Kelvin wave cascade since the quantum vortex core is just a point singularity, unlike 3D QT where the vortex core is a line or loop.To examine the spectral exponents of the compressible and incompressible kinetic energies in more detail we now discuss some high grid resolution runs: (a) grid 32768 2 with random phase initial conditions, and (b) grid 8192 2 with winding number n = 6 linear vortices in a uniform BEC background.(a) For the 32768 2 run, we choose initial conditions similar toNumasato et al. (2010), with parameters yielding a 'coherence length' ξ =(ag|ψ 0 | 2 ) −1/2 = 33.33 -even though, of course, initially there are no vortices in the system because of the random phase initial condition for the wave function.
Fig. 17.(a) The time evolution of the kinetic, quantum, internal and total energies on a 32768 2 grid with random phase initial condition.(b) Evolution of theincompressible energy (red) and compressible energy (blue).
Fig. 18. (a): The incompressible spectrum E inc (k) (red) and compressible spectrum E comp (k) (blue) on a 32768 2 grid t = 8000.The dashed vertical line indicates the location of k ξ , based on the qualitative notation of the coherence length ξ.The encircled dip in the compressible energy propagates towards the lower-k region, resembling a backward propagating pulse.(b) Incompressible kinetic energy spectrum ain Region I (k 0.01k ξ ) and Region II (0.01k ξ k 0.1k ξ ).The spectral exponents are s IC =+2.34 (red line) and s IC =+0.65 (green line).(c) Incompressible energy spectrum in Region II and Region III (0.1k ξ k k ξ ).Spectral exponents are: s IC = 0.65 (green line); s IC = −4.17(purple line).(d) Incompressible energy spectra in Region III and Region IV (k ξ k).Spectral exponents are: s IC = −4.17(purpleline) and s IC = −3.03(black line).
Fig. 18. (a): The incompressible spectrum E inc (k) (red) and compressible spectrum E comp (k) (blue) on a 32768 2 grid t = 8000.The dashed vertical line indicates the location of k ξ , based on the qualitative notation of the coherence length ξ.The encircled dip in the compressible energy propagates towards the lower-k region, resembling a backward propagating pulse.(b) Incompressible kinetic energy spectrum ain Region I (k 0.01k ξ ) and Region II (0.01k ξ k 0.1k ξ ).The spectral exponents are s IC =+2.34 (red line) and s IC =+0.65 (green line).(c) Incompressible energy spectrum in Region II and Region III (0.1k ξ k k ξ ).Spectral exponents are: s IC = 0.65 (green line); s IC = −4.17(purple line).(d) Incompressible energy spectra in Region III and Region IV (k ξ k).Spectral exponents are: s IC = −4.17(purpleline) and s IC = −3.03(black line).
Fig. 19.The (a) compressible kinetic energy spectrum and (b) incompressible kinetic energy spectrum for winding number n = 6 vortices in a uniform BEC gas.Grid 8192 2 Notice the triple cascade region for the compressible energy spectrum and the dual cascade spectrum for the incompressible energy.
Fig. 20.The |ϕ| isosurfaces very near the vortex core singularity for winding number n = 5:.Phase information is given on both the vortices and the boundary walls.The winding number n = 5 is evident from both the 5-fold periodicity around the each vortex as well as the 5 branch cuts emanating from each branch point on the boundary.By t = 3000, (b), the 5-fold degeneracy is removed with what seem like quantum Kelvin waves on the n = 1 cores.Basic phase coding : φ = 0 in blue, φ = 2π in red.Grid 2048 3 Fig. 21.The evolution of quantum core singularities from an initial set of 48 straight line vortices.(a) winding number n = 1 and the corresponding wall phase information at t = 0, (b) winding number n = 1 isosurface cores at t = 84000.(c) winding number n = 1 isosurface cores at t = 115000 = 0.5T Poin .The 2π phase changes at the core singularity intersections at the walls is very evident.(d) winding number n = 1 isosurface cores at t = 230000 = T Poin showing only small perturbative changes from the initial state given in (a).(e) winding number n = 2 and the corresponding wall phase information at t = 0 showing the confluence degeneracy.(f) The corresponding isosurface cores at t = 230000 = T Poin for winding number n = 2.The wall phase information is a simple perturbative change to that at t = 0, (e) -but there is much small scale vortex loops that has evolved due to the initial confluent degeneracy .Phase coding : φ = 0 blue, φ = 2π red.Grid 1200 3 Fig. 23.The time evolution of the E kin (t) (in blue) and E qu (t) (in red) for 0 ≤ t ≤ 250000 for (a) Winding Number n = 1, and (b) Winding Number n = 2. Grid 1200 3 .exhibits a triple cascade k −α with α ∼ 3.7 for small k,anα ∼ 6 for the intermediate cascade, and α ∼ 3.0 for the large-k cascade.At the intermittency, the large-k exponent increases to a noisy α ∼ 5.2 as well as a steeped intermediate wave number exponent.To investigate the cause of this intermittent loss of the incompressible k −3 spectrum, we then examined the vortex isosurfaces around this time interval, Fig.27.One notices that the loss of the k −3 corresponds to the apparent loss of vortex loops, i.e., of vortices.This would be consistent with the assumption that the incompressible kinetic energy spectrum of k −3 in the very large-k regime is due to the Fourier transform of an isolatedNore et al. (1997).As the vortex loops are reestablished, so is the incompressible k −3 kinetic energy spectrum.An alternative but somewhat more speculative explanation rests on the assumption that the incompressible k −3 spectrum is due to the quantum Kelvin wave cascade on the quantum vortices.As the quantum vortex loop shrink topologically, the Kelvin waves are inhibited and hence the loss of the k −3 spectrum.Moreover, if one looks at the time evolution of the mean kinetic E kin (t) and quantum E qu (t) energies one notices that this loss of the vortex loops occurs around the t ∼ 82000 around which the E kin (t) has a secondary Fig. 24.The initial incompressible (blue) E inc (k,0) and compressible (red) E comp (k,0) kinetic energy spectra for (a) Winding Number n = 1, and (b) Winding Number n = 2.The linear regression (blue dashed line) fit to the incompressible kinetic energy, E inc (k,0) ∼ k −α is: (a) α = 3.15, (b) α = 3.30.Grid 1200 3 .
Fig. 27. 4 snapshots of the winding number n = 1 singular vortex core isosurfaces at times (a) t = 78000, (b) t = 81000, (c) t = 82000, and (d) t = 88000.The phase information (blue: φ = 0, red: φ = 2π on the vortex core singularities clearly shows the 2π phase change in circumnavigating the vortex core loops.At t = 82000 there is a different morphology in the |ϕ| -isosurfaces.Grid 1200 3 Fig. 28.(a) Total kinetic energy spectrum, E TOT (k) − (blue dots), and quantum energy spectrum, E qu (k) − (red dots),att = 40000 for winding number n = 6 on a large 5760 3 grid.3 energy cascade regions can be readily distinguished, with both the total and quantum energy spectra being very similar.Dashed curves -linear regression fits.(b) Linear regression fits for different k-bands in the small-k region: for 30 < k < 200 (green dashed line), the spectra exponent α ∼ 1.30, while for 100 < k < 250 (red dashed line), the spectral exponent α ∼ 1.68.The standard deviation error is 0.06 in both cases.andcompressible components of the kinetic energy spectrum.Both E TOT (k) and E qu (k) have basically the same spectral properties.One sees 3 distinct energy cascade regions k −α :a small-k band with α ∼ 1.30, an intermediate-k band with steep slope α ∼ 7.76 and a large-k band with α ∼ 3.00.InYepez et al. (2009), we tried to identify these 3 regions as the small-k classical Kolmogorov casacde, followed by a semi-classical intermediate-k band (with non-universal exponent α) which is then followed a quantum Kelvin wave cascade for the very large-k band.Objections were raised against this interpretationKrstulovic & Brachet (2010);L'vov & Nazarenko (2010) based on (a) the kinetic energy spectrum of a single isolated vortex is k −3 for all k (for a straight line vortex all the kinetic energy is incompressible); and (b) using the standard definition of the coherence length ξ for the parameters chosen in our simulations, ξ > 2000 -i.e., it is claimed that we are investigating the physics of very strong vortex-vortex core overlapping wave functions.We counter that the definition of ξ is based on a boundary value solution of the GP equation for an isolated single vortex under the condition that the wave function asymptotes to the background value as one moves away from the vortex core.Our simulations are with periodic boundary conditions and we have no pointwise convergence of our wave function to some nice smooth 'background' value.While it can be argued that one should simply replace the usual background density ρ 0 by its spatial average < ρ >, the definition now of ξ becomes qualitative and does not handle large fluctuations about < ρ >.It is clear that we cannot categorically claim that the ubiquitous k −3 spectrum for the large-k band is due to quantum Kelvin wave cascade on single vortices -especially as this k −3 spectrum is also seen in our 2D QT simulations and in 2D there are no quantum Kelvin waves since the Fig. 29.(a) Energy spectra at t = 48000 for winding number n = 2. Blue circleskinetic energy, red squares -compressible kinetic energy, gold diamondsquantum energy.A triple cascade is quite evident in both the quantum and compressible kinetic energy spectra.These two spectra only deviate around the transition from the medium k to large k cascade, i.e., around k ∼ 300.(b) A blow-up of the transition in the incompressible kinetic energy pectrum from k −3.6 to k −3.0 around k ∼ 300.Grid 3072 3 incompressible kinetic energy (blue circles) has a slight bend in its spectral exponent around the wave number k ξ ∼ 300 from the large-k exponent of α ∼ 3.0 for k > k ξ to α ∼ 3.6 for k < k ξ .This bend in the incompressible kinetic energy spectrum occurs at the k ξ where the compressible and quantum energy spectra make their transition from the intermediate-k band large spectral exponent to the large-k band spectral exponent of α ∼ 3. It is interesting to note similar behavior in 2D QT, Fig.18(d) and Fig.19(b).The spectral exponents for E TOT (k) are: a Kolmogorov α ∼ 1.66 for the small-k band 15 < k < 90, a steep α ∼ 8.53 for the intermediate-k band 180 < k < 280 and exponent α ∼ 3.04 in the large-k band.

Table 4 .
Weak scaling of the BdG code on CRAYXE6/Hopper II (100 iterations)