**1. Introduction**

238 Advanced Fluid Dynamics

Tran-Thu, T. & Temperville, A., 1994. Numerical model of sediment transport in the wave-

Wilcox, D.C., 1993. *Turbulent modelling of CFD*. DCW Industries, La Canada, Calif.

Santos & A. Temperville, 271-281.

current interaction, in *Modelling of Coastal and Estuarine Processes*, Eds. F. Seabra-

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.

$$E(k) \sim k^{-5/3} \tag{1}$$

Independently, quantum turbulence (QT) was being studied in the low-temperature physics community on superfluid <sup>4</sup>*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):

$$i\partial\_t \varphi = -\nabla^2 \varphi + a(g|\varphi|^2 - 1)\varphi,\tag{2}$$

with the nonlinear term |*ϕ*| <sup>2</sup>*ϕ* 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

**1.1 Coherence length,** *ξ*

density <sup>|</sup>*ϕ*| ∼ <sup>√</sup>*ρ*<sup>0</sup>

is

**2.1 2D CT**

enstrophy *Z*

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 ∇ ∼ *<sup>ξ</sup>*−<sup>1</sup> and the nonlinear term by the asymptotic

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

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: *<sup>ϕ</sup>*(0) = 0 and *<sup>ϕ</sup>*(*<sup>x</sup>* <sup>→</sup> <sup>∞</sup>) = <sup>√</sup>*ρ*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) base their definition of coherence length on a *linear* perturbuation dispersion relation about uniform density. As Proment et. al.Proment et al. (2010) mention, the coherence length is strictly

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

In the *inviscid* limit, *ν* → 0, there are two constants of the motion - the energy *E* and the

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:

2, 2*Z* =

*d***r**|*ω*|

*d***r**|**v**|

√

*<sup>ϕ</sup>*(*x*) = <sup>√</sup>*ρ*<sup>0</sup> tanh(*x*/

defined only for a single isolated vortex and to extend this to QT one assumes that

problems the coherence length is at best a qualitative concept with possibly *ξ* = *ξ*(*t*).

where < *ρ*<sup>0</sup> > is the spatially averaged BEC density.

and write the 2D Navier-Stokes equation in the form

**2. Incompressible Classical Turbulence: 2D and 3D CT**

2*E* = 

In 2D incompressible CT it is convenient to introduce the scalar vorticity

*<sup>ξ</sup>* = (√*agρ*0)−<sup>1</sup> (6)

*<sup>ξ</sup>* = (√*ag* <sup>&</sup>lt; *<sup>ρ</sup>*<sup>0</sup> <sup>&</sup>gt;)−<sup>1</sup> (8)

*ω* = (∇ × **v**) · **ez** (9)

*<sup>∂</sup>t<sup>ω</sup>* <sup>+</sup> **<sup>v</sup>** · ∇*<sup>ω</sup>* <sup>=</sup> <sup>−</sup>*ν*∇2*<sup>ω</sup>* (10)

<sup>2</sup> (11)

2*ξ*) (7)

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

$$
\varphi = \sqrt{\rho} \, e^{i\theta/2}, \qquad \text{with} \qquad \mathbf{v} = \nabla \theta \tag{3}
$$

the GP Eq. (2) can be rewritten in conservation fluid form (*ρ* the mean density and **v** the mean velocity)

$$
\partial\_t \rho + \nabla \cdot (\rho \mathbf{v}) = 0 \tag{4}
$$

$$
\rho(\partial\_t \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{v}) = -2\rho \nabla(\mathbf{g} \,\rho - \frac{\nabla^2 \sqrt{\rho}}{\sqrt{\rho}}).\tag{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 <sup>∼</sup> *<sup>g</sup>ρ*<sup>2</sup> and a s-called quantum pressure ∼ −√*<sup>ρ</sup>* <sup>∇</sup>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<sup>3</sup> 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<sup>3</sup> 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.

#### **1.1 Coherence length,** *ξ*

2 Will-be-set-by-IN-TECH

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

Thus QT in a BEC gas will occur in a *compressible inviscid* quantum fluid whose pressure now has two major contributions: a barotropic pressure <sup>∼</sup> *<sup>g</sup>ρ*<sup>2</sup> and a s-called quantum pressure

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<sup>3</sup> 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<sup>3</sup> 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

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.

*∂tρ* + ∇ · (*ρ***v**) = 0 (4)

*<sup>i</sup>θ*/2, with **<sup>v</sup>** <sup>=</sup> <sup>∇</sup>*<sup>θ</sup>* (3)

<sup>√</sup>*<sup>ρ</sup>* ). (5)

*<sup>ϕ</sup>* = √*<sup>ρ</sup> <sup>e</sup>*

*<sup>ρ</sup>*(*∂t***<sup>v</sup>** <sup>+</sup> **<sup>v</sup>** · ∇**v**) = <sup>−</sup>2*ρ*∇(*<sup>g</sup> <sup>ρ</sup>* <sup>−</sup> <sup>∇</sup>2√*<sup>ρ</sup>*

spectrum of *k*−5/3, which after some time decayed away.

velocity)

∼ −√*<sup>ρ</sup>* <sup>∇</sup>2√*ρ*.

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 ∇ ∼ *<sup>ξ</sup>*−<sup>1</sup> and the nonlinear term by the asymptotic density <sup>|</sup>*ϕ*| ∼ <sup>√</sup>*ρ*<sup>0</sup>

$$\mathfrak{F} = (\sqrt{\mathfrak{a}\mathfrak{g}\rho\_0})^{-1} \tag{6}$$

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) is

$$\varphi(\mathbf{x}) = \sqrt{\rho\_0} \tanh(\mathbf{x}/\sqrt{2}\xi) \tag{7}$$

under the boundary conditions: *<sup>ϕ</sup>*(0) = 0 and *<sup>ϕ</sup>*(*<sup>x</sup>* <sup>→</sup> <sup>∞</sup>) = <sup>√</sup>*ρ*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) base their definition of coherence length on a *linear* perturbuation dispersion relation about uniform density. As Proment et. al.Proment et al. (2010) mention, the coherence length is strictly defined only for a single isolated vortex and to extend this to QT one assumes that

$$\xi = \left(\sqrt{\arg < \rho\_0 >}\right)^{-1} \tag{8}$$

where < *ρ*<sup>0</sup> > 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*).

#### **2. Incompressible Classical Turbulence: 2D and 3D CT**

#### **2.1 2D CT**

In 2D incompressible CT it is convenient to introduce the scalar vorticity

$$
\omega = (\nabla \times \mathbf{v}) \cdot \mathbf{e}\_\mathbf{z} \tag{9}
$$

and write the 2D Navier-Stokes equation in the form

$$
\partial\_t \omega + \mathbf{v} \cdot \nabla \omega = -\nu \nabla^2 \omega \tag{10}
$$

In the *inviscid* limit, *ν* → 0, there are two constants of the motion - the energy *E* and the enstrophy *Z*

$$2\mathbf{E} = \int d\mathbf{r} \left| \mathbf{v} \right|^2, \qquad 2\mathbf{Z} = \int d\mathbf{r} \left| \omega \right|^2 \tag{11}$$

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:

where the *σ* are the Pauli spin matrices

*β* :

at *x* + Δ*x*:

commute: [*C*, *S*] �= 0.

*σ<sup>x</sup>* = 0 1 1 0

Hence *C* is typically known as the square-root-of-swap gate.

*<sup>S</sup>*Δ*x*,0

*<sup>S</sup>*Δ*x*,1

*S*Δ*x*,0 = *n* + *e*

first consider the effect of the evolution operator *Uγ*[Ω(*x*)]

Since <sup>|</sup>Δ*x*| � 1 and *<sup>C</sup>*<sup>4</sup> <sup>=</sup> *<sup>I</sup>*, Eq.(23) yields *<sup>J</sup>*<sup>2</sup>

a function to be specified later.

yields the spinor equation

2010)

*σ<sup>y</sup>* = 0 −*i i* 0

*C*2 *α β* = *β α* 

*α*(*x*, *t*) *β*(*x*, *t*)

*α*(*x*, *t*) *β*(*x*, *t*)

can be expressed in an explicit unitary form by using the Paul spin matrices

*Uγ*[Ω(*x*)] = *J*

*ψ*(*x*, *t* + Δ*t*) = *ψ*(*x*, *t*) − *iε*

(−1)*γε*<sup>3</sup>

 ≡ 

 ≡ Now *C*<sup>2</sup> turns out to be just the swap operator because of its action on the amplitudes *α* and

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

The streaming operators shift just one of these amplitudes at *x* to a neighboring lattice point

The subscript *γ* = 0 on the streaming operator *S*Δ*x*,*<sup>γ</sup>* refers to shifting just the amplitude *α* while the subscript *γ* = 1 refers to just shifting the amplitude *β*. These streaming operators

where *n* = (1 − *σz*)/2, *n*¯ = (1 + *σz*)/2. Note that the collision and streaming operators do not

Now consider the following interleaved sequence of unitary collision and streaming operators

This is because the streaming operators are *O*(Δ*x*) deviations from the identity operator *I*. We

2 *<sup>x</sup><sup>γ</sup> J* 2 *<sup>y</sup><sup>γ</sup> J* 2 *<sup>z</sup>γe*

acting on the *γ* component of the 2-spinor *ψ*. Here *ε* � 1 is a perturbative parameter and Ω is

Using perturbation theory, it can be shown that the time advancement of *ψ* Yepez et al. (2009;

2 −1 2

<sup>4</sup> (*σ<sup>y</sup>* <sup>+</sup> *<sup>σ</sup>z*)∇3*ψ*(*x*, *<sup>t</sup>*) + <sup>O</sup>(*<sup>ε</sup>*

<sup>Δ</sup>*x∂<sup>x</sup> n*¯, *S*Δ*x*,1 = *n*¯ + *e*

*σ<sup>z</sup>* = 1 0 0 −1

*α*(*x* + Δ*x*, *t*) *β*(*x*, *t*)

 *α*(*x*, *t*) *β*(*x* + Δ*x*, *t*)

*Jx<sup>γ</sup>* = *<sup>S</sup>*−Δ*x*,*γCS*Δ*x*,*γ<sup>C</sup>* (23)

<sup>−</sup>*iε*2Ω(*x*)

*<sup>σ</sup>x*∇<sup>2</sup> <sup>+</sup> <sup>Ω</sup>

*ψ*(*x*, *t* + Δ*t*) = *Uγ*[Ω] *ψ*(*x*, *t*). (25)

*ψ*(*x*, *t*)+

4),

*<sup>x</sup><sup>γ</sup>* = *I* + *O*(Δ*x*), where *I* is the identity operator.

. (20)

. (21b)

<sup>Δ</sup>*x∂<sup>x</sup> n*, (22)

, (24)

. (19)

(21a)

(26)

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*) :


#### **2.2 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

$$E(k) \propto k^{-5/3} \tag{14}$$

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**

#### **3.1 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

$$\psi(x,t) = \begin{pmatrix} \alpha(x,t) \\ \beta(x,t) \end{pmatrix} \tag{15}$$

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

$$
\psi(x, t + \Delta t) = \mathcal{U}[\Omega] \,\psi(x, t). \tag{16}
$$

will reduce to the GP equation for the 1-particle boson wave function *ϕ* under the projection

$$(1,1)\cdot\psi = \varrho.\tag{17}$$

Consider the unitary collision operator *C* that locally entangles the complex amplitudes *α* and *β*

$$\mathbb{C} \equiv e^{i\frac{\pi}{4}\sigma\_{\mathbb{x}}(1-\sigma\_{\mathbb{x}})} = \begin{pmatrix} \frac{1-i}{2} & \frac{1+i}{2} \\ \frac{1+i}{2} & \frac{1-i}{2} \end{pmatrix} \, \tag{18}$$

where the *σ* are the Pauli spin matrices

4 Will-be-set-by-IN-TECH

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

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

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

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

**3. Unitary quantum lattice gas algorithm for the Gross-Pitaevskii equation**

spectrum for the subgrid *total* kinetic energy Menon & Genin (2010)

each cubic lattice position *x* we introduce the (complex) two-spinor field

*C* ≡ *e i π*

*ψ*(*x*, *t*) =

 *α*(*x*, *t*) *β*(*x*, *t*)

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

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

1−*i* 2 1+*i* 2 1+*i* 2 1−*i* 2

<sup>4</sup> *<sup>σ</sup><sup>x</sup>* (1−*σ<sup>x</sup>* ) =

*ψ*(*x*, *t* + Δ*t*) = *U*[Ω] *ψ*(*x*, *t*). (16)

(1, 1) · *ψ* = *ϕ*. (17)

, (18)

(15)

from the energy inverse cascade to small *k*: *E*(*k*) ∝ *k*<sup>−</sup>5/3; (12) from the enstrophy direct cascade to large *k*: *E*(*k*) ∝ *k*<sup>−</sup>3. (13)

*E*(*k*) ∝ *k*−5/3 (14)

of the energy spectrum *E*(*k*), where *E* = *dkE*(*k*) :

**2.2 3D CT**

spectrum

**3.1 The unitary algorithm**

the two spinor equation

*β*

$$
\sigma\_{\mathcal{X}} = \begin{pmatrix} 0 \ 1 \\ 1 \ 0 \end{pmatrix} \qquad \sigma\_{\mathcal{Y}} = \begin{pmatrix} 0 \ -i \\ i \ 0 \end{pmatrix} \qquad \sigma\_{\mathcal{z}} = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} . \tag{19}
$$

Now *C*<sup>2</sup> turns out to be just the swap operator because of its action on the amplitudes *α* and *β* :

$$\mathbb{C}^2 \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \begin{pmatrix} \beta \\ \alpha \end{pmatrix} . \tag{20}$$

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*:

$$S\_{\Delta x,0} \begin{pmatrix} \mathfrak{a}(x,t) \\ \beta(x,t) \end{pmatrix} \equiv \begin{pmatrix} \mathfrak{a}(x+\Delta x,t) \\ \beta(x,t) \end{pmatrix} \tag{21a}$$

$$S\_{\Delta x,1} \begin{pmatrix} \mathfrak{a}(x,t) \\ \beta(x,t) \end{pmatrix} \equiv \begin{pmatrix} \mathfrak{a}(x,t) \\ \beta(x+\Delta x,t) \end{pmatrix} . \tag{21b}$$

The subscript *γ* = 0 on the streaming operator *S*Δ*x*,*<sup>γ</sup>* 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

$$S\_{\Delta \mathbf{x},0} = \mathfrak{n} + e^{\Delta \mathbf{x} \partial\_{\mathbf{z}}} \,\, \mathfrak{n}, \qquad S\_{\Delta \mathbf{z},1} = \mathfrak{n} + e^{\Delta \mathbf{z} \partial\_{\mathbf{z}}} \,\, \mathfrak{n}, \tag{22}$$

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

$$J\_{\mathbf{x}\gamma} = \mathcal{S}\_{-\Delta \mathbf{z}, \gamma} \mathcal{C} \mathcal{S}\_{\Delta \mathbf{z}, \gamma} \mathcal{C} \tag{23}$$

Since <sup>|</sup>Δ*x*| � 1 and *<sup>C</sup>*<sup>4</sup> <sup>=</sup> *<sup>I</sup>*, Eq.(23) yields *<sup>J</sup>*<sup>2</sup> *<sup>x</sup><sup>γ</sup>* = *I* + *O*(Δ*x*), where *I* is the identity operator. This is because the streaming operators are *O*(Δ*x*) deviations from the identity operator *I*. We first consider the effect of the evolution operator *Uγ*[Ω(*x*)]

$$\mathcal{U}\_{\gamma}[\Omega(x)] = \mathcal{J}\_{x\gamma}^{2}\mathcal{J}\_{y\gamma}^{2}\mathcal{J}\_{z\gamma}^{2}e^{-i\varepsilon^{2}\Omega(x)},\tag{24}$$

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)

$$
\psi(x, t + \Delta t) = \mathcal{U}\_{\uparrow}[\Omega] \,\psi(x, t). \tag{25}
$$

yields the spinor equation

$$\begin{split} \psi(\mathbf{z}, t + \Delta t) &= \psi(\mathbf{z}, t) - i\varepsilon^2 \left[ -\frac{1}{2} \sigma\_{\mathbf{x}} \nabla^2 + \Omega \right] \psi(\mathbf{z}, t) + \\ &\frac{(-1)^\gamma \varepsilon^3}{4} (\sigma\_y + \sigma\_z) \nabla^3 \psi(\mathbf{z}, t) + \mathcal{O}(\varepsilon^4), \end{split} \tag{26}$$

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

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

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.

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 (*BEC*2) or for the *T* > 0 system which involves a coupling between the BEC ground state and the Bogoliubov modes (*BdG*). While the *BEC*2 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

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)

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 *BEC*2-code is excellent. The

the second order error term.

line is ideal timing. 100 iterations.

**3.2 Parallelization performance on various architectures**

in *BEC*2. The physics from these codes will be discussed elsewhere.

with *γ* = 0 or 1 and Δ*x* = *O*(*ε*). Since the order *ε*<sup>3</sup> term in Eq. (26) changes sign with *γ*, one can eliminate this term by introducing the symmetrized evolution operator

$$\mathcal{U}[\Omega] = \mathcal{U}\_1 \left[\frac{\Omega}{2}\right] \mathcal{U}\_0 \left[\frac{\Omega}{2}\right]. \tag{27}$$

rather than just a single *U<sup>γ</sup>* operator. Under diffusion ordering, Δ*t* = *O*(*ε*2) = Δ*x*2, the evolution equation

$$
\psi(x, t + \Delta t) = \mathcal{U}[\Omega(x)]\,\psi(x, t) \tag{28}
$$

leads to the spinor equation

$$i\partial\_t \psi(x,t) = \left[ -\frac{1}{2}\sigma\_\mathcal{X}\nabla^2 + \Omega \right] \psi(x,t) + \mathcal{O}(\varepsilon^2),\tag{29}$$

where the function Ω is still arbitrary. To recover the scalar GP Eq.2), one simply rescales the spatial grid ∇ → *<sup>a</sup>*−1∇, contracts the 2-component spinor field *<sup>ψ</sup>* to the (scalar) BEC wave function *ϕ*

$$
\varphi = (1,1) \cdot \psi = \mathfrak{a} + \mathfrak{z} \tag{30}
$$

and chooses Ω = *g*|*ϕ*| <sup>2</sup> <sup>−</sup> 1 :

$$i\partial\_l \varphi = -\nabla^2 \varphi + a(g|\varphi|^2 - 1)\varphi + \mathcal{O}(\varepsilon^2). \tag{31}$$

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

$$\delta \bar{\rho} = \int d\mathbf{x} (|\varphi|^2 - |\psi|^2) = \int d\mathbf{x} (\alpha^\dagger \beta + \alpha \beta^\dagger).$$

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 :


$$
\Omega\_N = \begin{pmatrix}
\cos[\Omega \Delta t] & -i\sin[\Omega \Delta t] \\
\end{pmatrix},
$$

$$
\text{such that}
\sum\_{\gamma} (\Omega\_N \cdot \psi)\_{\gamma} = e^{-i\Omega \Delta t} \varphi\_\gamma \text{ with } \Delta t = \varepsilon^2.
$$

6 Will-be-set-by-IN-TECH

with *γ* = 0 or 1 and Δ*x* = *O*(*ε*). Since the order *ε*<sup>3</sup> term in Eq. (26) changes sign with *γ*, one

 Ω 2 *U*<sup>0</sup> Ω 2 

*<sup>σ</sup>x*∇<sup>2</sup> <sup>+</sup> <sup>Ω</sup>

where the function Ω is still arbitrary. To recover the scalar GP Eq.2), one simply rescales the spatial grid ∇ → *<sup>a</sup>*−1∇, contracts the 2-component spinor field *<sup>ψ</sup>* to the (scalar) BEC wave

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

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

> <sup>2</sup>) =

 cos[ΩΔ*t*] <sup>−</sup>*<sup>i</sup>* sin[ΩΔ*t*] −*i* sin[ΩΔ*t*] cos[ΩΔ*t*]

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

<sup>2</sup> − |*ψ*<sup>|</sup>

. (27)

<sup>2</sup>), (29)

<sup>2</sup>). (31)

*ψ*(*x*, *t* + Δ*t*) = *U*[Ω(*x*)] *ψ*(*x*, *t*) (28)

*ϕ* = (1, 1) · *ψ* = *α* + *β* (30)

*ψ*(*x*, *t*) + O(*ε*

<sup>2</sup> <sup>−</sup> <sup>1</sup>)*<sup>ϕ</sup>* <sup>+</sup> <sup>O</sup>(*<sup>ε</sup>*

*d***x**(*α*†*β* + *αβ*†).

 ,

can eliminate this term by introducing the symmetrized evolution operator

Under diffusion ordering, Δ*t* = *O*(*ε*2) = Δ*x*2, the evolution equation

 −1 2

*<sup>i</sup>∂t<sup>ϕ</sup>* <sup>=</sup> −∇2*<sup>ϕ</sup>* <sup>+</sup> *<sup>a</sup>*(*g*|*ϕ*<sup>|</sup>

*d***x**(|*ϕ*|

introducing two modifications to the QLA algorithm described above :

Ω*<sup>N</sup>* =

−*i*ΩΔ*t*

• replace the scalar potential function Ω by the unitary non-diagonal matrix:

*ϕ*, with Δ*t* = *ε*2.

*i∂tψ*(*x*, *t*) =

<sup>2</sup> <sup>−</sup> 1 :

*δρ*¯ = 

• initialize the 2-spinor so that *α* = *Re*[*ϕ*], *β* = *i Im*[*ϕ*];

rather than just a single *U<sup>γ</sup>* operator.

leads to the spinor equation

function *ϕ*

2-spinor

such that ∑

*γ*

(Ω*<sup>N</sup>* · *ψ*)*<sup>γ</sup>* = *e*

and chooses Ω = *g*|*ϕ*|

parameters and initial amplitudes.

*U*[Ω] = *U*<sup>1</sup>

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.

#### **3.2 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 (*BEC*2) or for the *T* > 0 system which involves a coupling between the BEC ground state and the Bogoliubov modes (*BdG*). While the *BEC*2 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 *BEC*2. The physics from these codes will be discussed elsewhere.

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.

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 *BEC*2-code is excellent. The

GRID CORES WALLCLOCK (IDEAL: 209.2S) % PEAK <sup>3</sup> 216 209.2 9.5% <sup>3</sup> 1728 212.7 9.4% <sup>3</sup> 13824 212.1 9.3% <sup>3</sup> 46656 212.7 9.0% <sup>3</sup> 64000 213.3 9.1%

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

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 (64003) it shows superlinear scaling (Fig. 1). Similar results are achieved on the CRAY *XT*5*Jaguarp f* , tested to 110592 cores and grids 60003. The *BdG*-code, with its 4 <sup>×</sup> 4 entangling matrices on the 2 qubits, shows quite strong superlinear scaling on both the CRAY *XT*5 and *XE*6, testing to

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<sup>3</sup> in a 3D simulations. For ideal scaling the wallclock time to completion should remain invariant. The weak scaling of these codes are given in Tables I - IV, and again are excellent, with fluctuations

We first consider 2D quantum turbulence. From the GP Eq.2, the total energy in 2D QT is

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)

To examine the effect of compressibility on the GP Eq.(2), *EK* 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 √*ρ*

kinetic energy: *Ekin* <sup>=</sup> <sup>1</sup>

internal energy: *Eint* <sup>=</sup> *ag*

quantum energy: *Eqnt* = 2

<sup>2</sup> <sup>+</sup> *agρ*<sup>2</sup> <sup>+</sup> <sup>2</sup>|∇√*ρ*<sup>|</sup>

2 

<sup>2</sup> <sup>−</sup> *<sup>a</sup><sup>ρ</sup>* 

*dx*2*ρ*|**v**<sup>|</sup>

*dx*2|∇√*ρ*<sup>|</sup>

= *const*. (32)

<sup>2</sup> (33)

2. (35)

*dx*2*ρ*<sup>2</sup> (34)

4800<sup>3</sup> 110592 214.0

216000 cores on *Jaguarp f* and to 150000 cores on *HopperI I* (Fig. 2).

**4. 2D QT**

conserved

*ETOT* =

 *dx*<sup>2</sup> 1 2 *ρ*|**v**|

in timings typically below 5% as one moves from 216 cores to 110592 cores.

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

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 84003, 100 iterations.


Table 1. Weak scaling of the BEC2 code on IBM/BlueGene Intrepid (100 iterations)


Table 2. Weak scaling of the BEC2 code on CRAY XT5/Jaguarpf (100 iterations)


Table 3. Weak scaling of the BdG code on CRAY XT5/Jaguarpf (100 iterations)


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

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 (64003) it shows superlinear scaling (Fig. 1). Similar results are achieved on the CRAY *XT*5*Jaguarp f* , tested to 110592 cores and grids 60003. The *BdG*-code, with its 4 <sup>×</sup> 4 entangling matrices on the 2 qubits, shows quite strong superlinear scaling on both the CRAY *XT*5 and *XE*6, testing to 216000 cores on *Jaguarp f* and to 150000 cores on *HopperI I* (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<sup>3</sup> in a 3D simulations. For ideal scaling the wallclock time to completion should remain invariant. The weak scaling of these codes are given in Tables I - IV, and again are excellent, with fluctuations in timings typically below 5% as one moves from 216 cores to 110592 cores.

## **4. 2D QT**

8 Will-be-set-by-IN-TECH

Fig. 2. Strong scaling of the *T* > 0 BdG code on CRAY XT5/Jaguarpf and CRAY

Table 1. Weak scaling of the BEC2 code on IBM/BlueGene Intrepid (100 iterations)

Table 2. Weak scaling of the BEC2 code on CRAY XT5/Jaguarpf (100 iterations)

Table 3. Weak scaling of the BdG code on CRAY XT5/Jaguarpf (100 iterations)

<sup>3</sup> 216 246.8 <sup>3</sup> 1728 245.6 <sup>3</sup> 13824 249.6 <sup>3</sup> 110592 266.1

<sup>3</sup> 216 261.0 <sup>3</sup> 1728 261.8 <sup>3</sup> 13824 263.9 <sup>3</sup> 46656 268.4

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 84003, 100 iterations.

> GRID CORES WALLCLOCK (IDEAL: 143.6S) % PEAK <sup>3</sup> 64 143.6 19.5% <sup>3</sup> 512 144.5 19.4% <sup>3</sup> 4096 144.8 19.4% <sup>3</sup> 32768 150.8 18.6%

GRID CORES WALLCLOCK (IDEAL: 246.8S)

GRID CORES WALLCLOCK (IDEAL: 261.0S)

We first consider 2D quantum turbulence. From the GP Eq.2, the total energy in 2D QT is conserved

$$E\_{TOT} = \int d\mathbf{x}^2 \left(\frac{1}{2}\rho|\mathbf{v}|^2 + ag\rho^2 + 2|\nabla\sqrt{\rho}|^2 - a\rho\right) = const.\tag{32}$$

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)

$$\text{kinetic energy: } E\_{kin} = \frac{1}{2} \int d\mathbf{x}^2 \rho |\mathbf{v}|^2 \tag{33}$$

$$\text{internal energy: } E\_{\text{int}} = ag \int d\mathbf{x}^2 \rho^2 \tag{34}$$

$$\text{quantum energy: } E\_{\text{qut}} = 2 \int d\mathbf{x}^2 |\nabla \sqrt{\rho}|^2. \tag{35}$$

To examine the effect of compressibility on the GP Eq.(2), *EK* 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 √*ρ*

These segments are then patched together to form a new picture according to their updated coordinate [*qt*+1, *pt*+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 *TP* = 114 iterations, Fig.3, while for the higher resolution *m* = 300, *TP* = 300 iterations, Fig.4. What is quite interesting is the image after *t* = *TP*/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

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

(a) t=0 (b) t=57 (c) t=114

(a) t=0 (b) t=150 (c) t=300

Fig. 4. Poincaré recurrence under the Arnold cat mapping (II). Picture resolution: 300 × 300 pixels. Poincaré recurrence time *TP* = 300. At the semi Poincaré time, the picture shows no

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

The plots of amplitude, vorticity and phase in the following sections will adopt the "thermal"

**4.1.1 Four 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

scheme: *blue* stands for low values while *red* stands for high values.

*<sup>γ</sup>* <sup>=</sup> �*Eint*(*t*)� / �*Ekin*(*t*)� <sup>O</sup>(10−1) (44)

Fig. 3. Poincaré recurrence under the Arnold cat mapping (I). Picture resolution: 74 × 74 pixels. Poincaré recurrence time *TP* = 114. Notice at the semi-Poincaré time t=57, the picture

image at *t* = *TP*/2 = 22.

is almost a point inversion of the initial condition.

strong symmetry to the initial condition.

energy is sufficiently small:

**4.1 Poincaré recurrence in 2D QT**

rather than *<sup>ρ</sup>*. Unlike CFD, in QT at the vortex core the density *<sup>ρ</sup>* <sup>→</sup> 0 with <sup>√</sup>*ρ***<sup>v</sup>** <sup>→</sup> *const*.. The longitudinal component and the transverse component of **q**˜ are

$$\tilde{\mathbf{q}}\_l = \frac{\mathbf{k} \cdot \tilde{\mathbf{q}}}{|\mathbf{k}|^2} \tag{36}$$

$$
\tilde{\mathbf{q}}\_l = \tilde{\mathbf{q}} - \frac{\mathbf{k} \cdot \tilde{\mathbf{q}}}{|\mathbf{k}|^2} \tag{37}
$$

so that **k** × **q**˜*<sup>l</sup>* = 0, **k** · **q**˜*<sup>t</sup>* = 0. Correspondingly the compressible energy and incompressible energy can be defined as

$$\text{compressible}: E\_{\mathbb{C}} = (2\pi)^2 \int dk^2 |\vec{\mathbf{q}}\_l|^2 \tag{38}$$

$$\text{incompressible}: E\_{IC} = \left(2\pi\right)^{2} \int dk^{2} \left|\tilde{\mathbf{q}}\_{t}\right|^{2}. \tag{39}$$

Consequently the energy densities of compressible and incompressible energy become

$$
\varepsilon\_{c,ic}(k) = k \int\_0^{2\pi} d\theta |\mathbf{q}\_{\mathbf{l},\mathbf{t}}(k,\theta)|^2 \,\prime \,\tag{40}
$$

using polar coordinates *k*, *θ*. For 2D QT it is useful to introduce a renormalized vorticity Horng et al. (2009)

$$
\omega\_{\mathfrak{q}} = (\nabla \times [\sqrt{\rho} \mathbf{v}]) \cdot \mathbf{e}\_z. \tag{41}
$$

The two components of *ω<sup>q</sup>* are:


Consequently, the time evolution of the density weighted enstrophy, *Zq* <sup>=</sup> *dx*2|*ωq*<sup>|</sup> 2, is an 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 *TP* 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' [*qt*, *pt*], with *t* = 1....*m*. The Arnold cat map is applied to these segments:

$$q\_{t+1} = \text{Mod}(2q\_t + p\_{t\prime}m) \tag{42}$$

$$p\_{t+1} = Mod(q\_t + p\_{t\prime}m).\tag{43}$$

10 Will-be-set-by-IN-TECH

rather than *<sup>ρ</sup>*. Unlike CFD, in QT at the vortex core the density *<sup>ρ</sup>* <sup>→</sup> 0 with <sup>√</sup>*ρ***<sup>v</sup>** <sup>→</sup> *const*.. The

**<sup>q</sup>**˜*<sup>t</sup>* <sup>=</sup> **<sup>q</sup>**˜ <sup>−</sup> **<sup>k</sup>** · **<sup>q</sup>**˜

so that **k** × **q**˜*<sup>l</sup>* = 0, **k** · **q**˜*<sup>t</sup>* = 0. Correspondingly the compressible energy and incompressible

*dθ*|**q**˜**l**,**t**(*k*, *θ*)|

*dk*2|**q**˜*l*<sup>|</sup>

*dk*2|**q**˜*t*<sup>|</sup>

*qt*+<sup>1</sup> = *Mod*(2*qt* + *pt*, *m*) (42) *pt*+<sup>1</sup> = *Mod*(*qt* + *pt*, *m*). (43)

compressible : *EC* = (2*π*)<sup>2</sup>

incompressible : *EIC* = (2*π*)<sup>2</sup>

*ω<sup>q</sup>* = (∇ × [

Consequently, the time evolution of the density weighted enstrophy, *Zq* <sup>=</sup> *dx*2|*ωq*<sup>|</sup>

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 *TP* 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' [*qt*, *pt*],

*εc*,*ic*(*k*) = *k*

• <sup>√</sup>*ρ*∇ × **<sup>v</sup>** · **<sup>e</sup>***<sup>z</sup>* <sup>∝</sup> *<sup>δ</sup>*(**<sup>r</sup>** <sup>−</sup> **<sup>r</sup>**0), which pinpoints the locations of the vortices.

with *t* = 1....*m*. The Arnold cat map is applied to these segments:

Consequently the energy densities of compressible and incompressible energy become

 2*π* 0

using polar coordinates *k*, *θ*. For 2D QT it is useful to introduce a renormalized vorticity Horng

• (∇√*ρ*) <sup>×</sup> **<sup>v</sup>** · **<sup>e</sup>***z*, with major contributions coming from density variations near vortices and

<sup>|</sup>**k**|<sup>2</sup> (36)

<sup>|</sup>**k**|<sup>2</sup> (37)

<sup>2</sup> (38)

2. (39)

2, is an

2, (40)

<sup>√</sup>*ρ***v**]) · **<sup>e</sup>***z*. (41)

**<sup>q</sup>**˜*<sup>l</sup>* <sup>=</sup> **<sup>k</sup>** · **<sup>q</sup>**˜

longitudinal component and the transverse component of **q**˜ are

energy can be defined as

The two components of *ω<sup>q</sup>* are:

from sound and shock waves;

et al. (2009)

These segments are then patched together to form a new picture according to their updated coordinate [*qt*+1, *pt*+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 *TP* = 114 iterations, Fig.3, while for the higher resolution *m* = 300, *TP* = 300 iterations, Fig.4. What is quite interesting is the image after *t* = *TP*/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* = *TP*/2 = 22.

Fig. 3. Poincaré recurrence under the Arnold cat mapping (I). Picture resolution: 74 × 74 pixels. Poincaré recurrence time *TP* = 114. Notice at the semi-Poincaré time t=57, the picture is almost a point inversion of the initial condition.

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

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:

$$\gamma = \langle \mathcal{E}\_{\text{int}}(t) \rangle \, / \, \langle \mathcal{E}\_{\text{kin}}(t) \rangle \lesssim \mathcal{O}(10^{-1}) \tag{44}$$

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.

#### **4.1 Poincaré recurrence in 2D QT**

#### **4.1.1 Four 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

the appearance of up-side-down point inversion pattern of the Arnold's cat mapping at a

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

(a) |*ψ*|,t=0 (b) t=10500 (c) t=21000 (d) t=31300 (e) t=41900=*TP*

(f) *ωq*,t=0 (g) t=10500 (h) t=21000 (i) t=31300 (j) t=41900=*TP*

(k) *θ*,t=0 (l) t=10500 (m) t=21000 (n) t=31300 (o) t=41900=*TP*

Fig. 7. Poincaré recurrence for vortices with winding number *n* = 1. At Poincaré time *TP* = 41900, all three distributions replicate the initial state except for some background noise, such as the small vortices at the boundaries appearing on the phase plot (o). At the semi-Poincaré time, due to periodic boundary conditions, the system is almost identical to its initial state provided the origins shift by (−*L*/2, −*L*/2), i.e., a point inversion. Grid length

**4.1.2 Four vortices with winding number** *n* = 2 **embedded in a Gaussian BEC background** 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é

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

recurrence occurs, more background noise is to be expected, and is seen in Fig.9.

resolution of 74-pixels.

*L* = 512.

function takes the form of products of individual quantum vortices, modulated by the inhomogeneous Gaussian background:

$$
\psi(\mathbf{r}) = \hbar \, e^{-a \, w\_{\S} \, r^2} \prod\_{i} \phi\_i(\mathbf{r} - \mathbf{r}\_i).
$$

*φ<sup>i</sup>* = *tanh*( <sup>√</sup>*a*|**<sup>r</sup>** <sup>−</sup> **<sup>r</sup>***i*|)*ei n Arg*(**r**−**r***i*) 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

Fig. 5. Four vortices with winding number 1 embedded in a Gaussian background. *h* = 0.05, *a* = 0.01, *wg* = 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, *Eint*(*t*), c.f. Fig.6. To further investigate such structure, we

Fig. 6. Evolution of internal energy for winding number 1. The similarities exhibited by the two peaks around t=21000 and t=41900 clearly indicates a periodical structure.

plot the wave function's amplitude, vorticity and phase distribution. At Poincaré recurrence time *TP* = 41900, the wave function recovers the initial condition except for some background noise. At semi-Poincaré recurrence time, *t* = *TP*/2 = 21000, the wave function would be the same as initial condition if the origin of domain was shifted by (*L*/2, *L*/2). This resembles 12 Will-be-set-by-IN-TECH

function takes the form of products of individual quantum vortices, modulated by the

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

(a) *amplitude*|*ψ*| (b) *vorticity ω<sup>q</sup>* (c) *phase θ*

*h* = 0.05, *a* = 0.01, *wg* = 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

steps is 50000. The ratio parameter *γ* = 0.018. A recurrence structure is clearly visible via the evolution of internal energy, *Eint*(*t*), c.f. Fig.6. To further investigate such structure, we

Fig. 6. Evolution of internal energy for winding number 1. The similarities exhibited by the

plot the wave function's amplitude, vorticity and phase distribution. At Poincaré recurrence time *TP* = 41900, the wave function recovers the initial condition except for some background noise. At semi-Poincaré recurrence time, *t* = *TP*/2 = 21000, the wave function would be the same as initial condition if the origin of domain was shifted by (*L*/2, *L*/2). This resembles

two peaks around t=21000 and t=41900 clearly indicates a periodical structure.

Fig. 5. Four vortices with winding number 1 embedded in a Gaussian background.

∏ *i*

<sup>√</sup>*a*|**<sup>r</sup>** <sup>−</sup> **<sup>r</sup>***i*|)*ei n Arg*(**r**−**r***i*) is the approximated wave function of a single 2D vortex with

*φi*(**r** − **r***i*).

*ψ*(**r**) = *h e*−*a wg <sup>r</sup>*<sup>2</sup>

inhomogeneous Gaussian background:

*φ<sup>i</sup>* = *tanh*(

the singularities.

the appearance of up-side-down point inversion pattern of the Arnold's cat mapping at a resolution of 74-pixels.

Fig. 7. Poincaré recurrence for vortices with winding number *n* = 1. At Poincaré time *TP* = 41900, all three distributions replicate the initial state except for some background noise, such as the small vortices at the boundaries appearing on the phase plot (o). At the semi-Poincaré time, due to periodic boundary conditions, the system is almost identical to its initial state provided the origins shift by (−*L*/2, −*L*/2), i.e., a point inversion. Grid length *L* = 512.

### **4.1.2 Four vortices with winding number** *n* = 2 **embedded in a Gaussian BEC background**

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

(a) *θ*, t=10500 (b) *θ*, t=10600

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

coefficients *ai*,*<sup>j</sup>* are determined by the enforced continuity at the corners. Since there are 16 unknown coefficients, 16 equations are needed to determine these *ai*,*j*. Usually one can enforce continuity for *p*(*x*, *y*), *∂<sup>x</sup> p*(*x*, *y*), *∂<sup>y</sup> p*(*x*, *y*), *∂x*,*<sup>y</sup> p*(*x*, *y*). To generate random phase in

• Generate 4 pseudo-random numbers at the 4 corners of each sub-square. These random

For a grid of 5122 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

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 *TP*/2 and *TP*. At *t* = 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* = 0 or at *TP*/2, (e), although a considerable

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 *TP*/2 = 21000, however, at *TP* = 41900, the initial condition, which is vortex free, can not be reproduced. For *γ* = 0.133,

3 ∑ *i*=0

3 ∑ *j*=0

*ai*,*jx<sup>i</sup> yj* . The

Fig. 10. Vortices generation via the Kelvin-Helmholtz instability. The region depicted in these plots is a blow-up of [−256, −128] × [−128, 128]. (a) phase *θ* at t=10500; (b) *θ* at t=10600. At *t* = 10600, a new pair of counter-rotating vortices (pointed out by the black arrows) are

generated between neighboring counter-rotating vortices.

the domain *L*2, the following procedure is followed:

• Enforce periodicity at the domain boundaries; • Solve *ai*,*<sup>j</sup>* belonging to each single sub-square.

amount of sound waves.

defined on a unit square is approximated by polynomials: *p*(*x*, *y*) =

• Discretize the domain into *m* × *m* squares. *m* here is the level of randomness;

numbers correspond to *p*, *∂<sup>x</sup> p*, *∂<sup>y</sup> p*, *∂x*,*<sup>y</sup> p* at the corners of this sub-square;

time can be clearly identified by the abrupt energy exchanges (i.e., the spikes).

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. 9. Poincaré recurrence for winding number*n* =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.
