Quantum turbulence simulations using the Gross-Pitaevskii equation: high-performance computing and new numerical benchmarks
Michikazu Kobayashi, Philippe Parnaudeau, Francky Luddens, Corentin Lothode, Luminita Danaila, Marc Brachet, Ionut Danaila
QQuantum turbulence simulations using theGross-Pitaevskii equation: high-performancecomputing and new numerical benchmarks
Michikazu Kobayashi , Philippe Parnaudeau , Francky Luddens ,Corentin Lothod´e , Luminita Danaila , Marc Brachet and Ionut Danaila ∗ Department of Physics, Kyoto University, Japan Institut Pprime, CNRS,Universit´e de Poitiers - ISAE-ENSMA - UPR 3346,86962 Futuroscope Chasseneuil Cedex, France Universit´e de Rouen Normandie,Laboratoire de Math´ematiques Rapha¨el Salem,F-76801 Saint- ´Etienne-du-Rouvray, France Universit´e de Rouen Normandie,CORIA, CNRS UMR 6614,F-76801 Saint- ´Etienne-du-Rouvray, France Laboratoire de Physique de l’ ´Ecole Normale Sup´erieure,ENS, Universit´e PSL, CNRS, Sorbonne Universit´e,Universit´e de Paris, F-75005 Paris, France ∗ Corresponding author: [email protected] 10, 2020 a r X i v : . [ phy s i c s . f l u - dyn ] M a r bstract This paper is concerned with the numerical investigation of Quan-tum Turbulence (QT) described by the Gross-Pitaevskii (GP) equation.Numerical simulations are performed using a parallel (MPI-OpenMP)code based on a pseudo-spectral spatial discretization and second ordersplitting for the time integration. We start by revisiting (in the frame-work of high-performance/high-accuracy computations) well-knownGP-QT settings, based on the analogy with classical vortical flows:Taylor-Green (TG) vortices and Arnold-Beltrami-Childress (ABC) flow.Two new settings are suggested to build the initial condition for theQT simulation. They are based on the direct manipulation of thewave function by generating a smoothed random phase (SRP) field, orseeding random vortex rings (RVR) pairs. The new initial conditionshave the advantage to be simpler to implement than the TG andABC approaches, while generating statistically equivalent QT fields.Each of these four GP-QT settings is described in detail by definingcorresponding benchmarks that could be used to validate/calibratenew GP codes. We offer a comprehensive description of the numericaland physical parameters of each benchmark. We analyze the resultsin detail and present values, spectra and structure functions of mainquantities of interest (energy, helicity, etc.) that are useful to describethe turbulent flow. Some general features of QT are identified, despitethe variety of initial states.
Keywords: Quantum Turbulence, Gross-Pitaevskii equation, Taylor-Green,ABC, parallel computing, spectral method.
The study of quantum fluids, realized in superfluid helium and atomic Bose-Einstein condensates (BEC), has become a central topic in various fields ofphysics, such as low temperature physics, fluid dynamics of inviscid flows,quantum physics, statistical physics, cosmology, etc. One of the strikingfeatures of quantum fluids is the nucleation of vortices with quantized (fixed)circulation, when an external forcing is applied (rotation, stirring, etc). Theobservation of quantized vortices, as a signature of the superfluid (zero-viscosity) nature of these flow systems, was extensively explored in differentexperimental settings of superfluid helium or BEC. Configurations with a2arge number of quantized vortices tangled in space can evolve to QuantumTurbulence (QT), generally referred to as vortex tangle turbulence . WhileQT in superfluid helium has been largely studied in the last two decades(see dedicated volumes Vinen and Niemela (2002); Barenghi et al. (2001);Barenghi and Sergeev (2008); Halperin and Tsubota (2009)), only recentexperimental and theoretical studies (Henn et al., 2010; Seman et al., 2011;Kwon et al., 2014; Navon et al., 2016) reported different possible routes toQT in BEC.The focus of this paper is on the numerical study of QT in superfluidhelium. A large body of literature in this field is based on the macroscopicdescription of the superfluid using the Gross-Pitaevskii (GP) mean-fieldequation. The GP equation is a nonlinear Schr¨odinger (NLS) equation withcubic nonlinearity. It describes, in the theoretical limit of absolute zerotemperature, the time evolution of the single-particle complex wave function ψ (identical for all particles). Quantized vortices appear then as topologicalline defects, resulting from the U (1) symmetry breaking of the phase shiftof ψ . QT in superfluid helium based on the GP model (denoted hereafteras GP-QT) implies complicated space and time interactions between a largenumber of quantized vortices. For a comprehensive description of differentmodels of QT, see recent reviews by Halperin and Tsubota (2009); Brachet(2012); Barenghi et al. (2014); Tsubota et al. (2017).There are several challenges when setting a numerical simulation to inves-tigate GP-QT: (i) generate a physically and mathematically sound initial state with manyquantized vortices that finally evolve to a statistically steady state of QT, (ii) use accurate numerical methods that preserve the invariants of the GPequation when long time-integration is necessary, (iii) design numerical codesaffording large grid resolutions, necessary to accurately capture the dynam-ics of vortices and (iv) compute appropriate (statistical) diagnostic tools toanalyze the superfluid flow evolution.We use in this contribution a modern parallel (MPI-OpenMP) numericalcode satisfying the requirements (ii) and (iii) . The code is called GPS(Gross-Pitaevskii Simulator) (Parnaudeau et al., 2015) and is based on aFourier-spectral space discretization and up-to-date numerical methods: asemi-implicit backward-Euler scheme with Krylov preconditioning for thestationary GP equation (Antoine and Duboscq, 2014) and various schemes(Strang splitting, relaxation, Crank-Nicolson) for the real-time GP equation(Antoine et al., 2013). The GPS code offers a solid framework to address in3etail challenges (i) and (iv) , for which we review previous models and bringnew contributions.As in classical turbulence (CT), the numerical and physical accuracyof the initial condition is crucial in computing properties of numericallygenerated QT. Using the hydrodynamic analogy for the GP model (through theMadelung transform, as explained below), pioneering numerical simulations ofGP-QT (Nore et al., 1993, 1997a; Abid et al., 2003) suggested initial conditionsand statistical analysis tools inspired from CT. A velocity field, derived fromthe well-known classical flow with Taylor-Green (TG) vortices, was imposed tothe superfluid flow. An initial wave function, with nodal lines correspondingto vortex lines of the velocity, was thus generated. This initial wave functionwas then used in the Advective Real Ginzburg-Landau equation (ARGLE),equivalent to the imaginary-time GP equation with Galilean transformation(see also below), to reduce the acoustic emission of the initial field. The resultof the ARGLE procedure was finally used as initial field for the time-dependentGP simulation. A similar approach was more recently used by replacing theTG vortices with the Arnold-Beltrami-Childress (ABC) classical vortex flow(Clark di Leoni et al., 2016, 2017). If this approach is well suited to controlthe hydrodynamic characteristics of the initial superfluid flow (Mach number,helicity), it involves supplemental technicalities and computations throughthe ARGLE procedure. We suggest in this paper two new approaches togenerate the initial condition for the GP-QT simulations, based on the directmanipulation of the wave function. The ARGLE procedure is thus avoided.The first method prescribes a smoothed random-phase (SRP) for the wavefunction, while the second one generates random vortex rings (RVR). Thetwo new methods, which are simply to implement, are shown to develop QTfields with similar statistical properties as those obtained using the TG orABC classical initial conditions. Nevertheless, the dynamics of the superfluidflow is different. Compared to TG and ABC cases, in the SRP case the initialfield is vortex free and dominated by the compressible kinetic energy; vorticesnucleate progressively and do not display long vortex lines. For the RVRflow, evolution is opposite to that observed for the SRP case: in the earlystages of the time evolution the incompressible kinetic energy is dominant;the compressible kinetic energy then starts to increase due to sound emissionsthrough vortex reconnections. However, like in well-documented TG andABC cases, a Kolmogorov-like scaling of the incompressible kinetic energyspectrum is obtained for the new SRP and RVR cases.Concerning the analysis of the QT field, we present classical diagnostic4ools (inspired from CT), as the energy decomposition and associated spectra(Nore et al., 1993, 1997a) and also new ones, as the second-order structurefunction, not reported in the previously cited studies. This supplements thestatistical description of the superfluid flow. We also carefully investigatethe influence of numerical parameters (as the grid resolution of a vortex, themaximum resolved wave-number and the computed local Mach number) onthe characteristic of the QT. This topic is generally very briefly addressed inphysical papers on GP-QT.Starting from the observation that in previously published studies of GP-QT, the focus was mainly given to the physics of turbulence, this paper is alsointended to define in detail numerical benchmarks in the framework of parallelcomputing. We start by revisiting classical GP-QT settings (based on TG andABC flows). New results obtained with our high-performance/high-accuracyparallel code are compared with available data in the literature. We thenpresent the new numerical benchmarks, based on random phase fields orrandom vortex rings generation. The new benchmarks offer a new perspectivein comparing the classical CT based approaches to more GP-oriented models.For all benchmarks, we offer a comprehensive description of the numericaland physical parameters and give checkpoint values for validating each stepof the simulation. This could be useful to verify new numerical methodologiesor tune/validate new modern GP numerical codes.The organization of the paper is as follows. Section 2 introduces the GPmean field equation, its stationary version and the hydrodynamic analogy,derived using the Madelung transform. The main conservation laws areemphasized and the Bogoliubov-de Gennes (BdG) linearized system for theresponse of the system to small oscillations is derived. Section 3 introduces themain characteristics of the QT flow that are used to set numerical simulations:the healing length, the sound velocity, the Bogoliubov dispersion relation,the compressible and incompressible kinetic energy and the helicity. Thenotion of quantum vortex and its representation on a computational gridare also discussed. The numerical method used to solve the GP equationand the associated GPS code are described in Section 4. An extendedsubsection is devoted to the derivation of dimensionless equations, by settinga unified framework covering different existing approaches in the literature.The particular numerical methods used in this study to advance the GP wavefunction in imaginary-time (ARGLE) or real-time (GP) are also described.Section 5 presents in detail four different approaches to generate the initial fieldfor the simulation of decaying GP-QT: Taylor-Green (TG), Arnold-Beltrami-5hildress (ABC), smoothed random phase (SRP) and random vortex rings(RVR). Each method is associated to a benchmark. The results obtained forthe four benchmarks are discussed in Section 6. We present values, spectraand structure functions of main quantities of interest (energy, helicity, etc.)that could be useful to benchmark numerical codes simulating QT with theGP model. Finally, the main features of the benchmarks and their possibleextensions are summarized in Section 7. We present in this section the GP model used for the numerical simulationof QT. The equations in this part apply to both BEC and superfluid heliumQT. We first introduce the GP mean field equation and its stationary version.The hydrodynamic analogy, important for relating QT to CT, is then derivedusing the Madelung transform. The main conservation laws are emphasizedas important checkpoints for numerical simulations. Finally, we derive theBogoliubov-de Gennes (BdG) linearized system for the response of the systemto small oscillations. The BdG formalism will be used in the next section toderive the dispersion relation for the QT flow.
In the zero-temperature limit, the superfluid system of weakly interactingbosons of mass m , is described by the Gross-Pitaevskii mean field equation(Pitaevskii and Stringari, 2003): i (cid:126) ∂∂t ψ ( x , t ) = − (cid:126) m ∇ + V trap ( x ) + g | ψ ( x , t ) | ! ψ ( x , t ), (1)where V trap is the external trapping potential and g the non-linear interactioncoefficient g = 4 π (cid:126) a s m , (2)with a s the s-wave scattering length for the binary collisions within the system.The complex wave function ψ is generally represented as (Madelungtransformation): ψ = q n ( x , t ) e iθ ( x , t ) , n ( x , t ) = | ψ ( x , t ) | = ψ ( x , t ) ψ ∗ ( x , t ), (3)6ith n the atomic density and θ the phase of the order parameter. We denoteby ψ ∗ the complex conjugate of the wave function. Inserting (3) in (1) weobtain by separating the imaginary and real parts, the time-evolution equationfor the atomic density: ∂n∂t + (cid:126) m (cid:16) ∇ n · ∇ θ + n ∇ θ (cid:17) = 0, or ∂n∂t + (cid:126) m ∇ · ( n ∇ θ ) = 0, (4)and the equation for the phase: (cid:126) ∂θ∂t + (cid:126) m ( ∇ θ ) + V trap + gn − (cid:126) m √ n ∇ ( √ n ) = 0. (5)The last term in the left-hand side of (5) is called quantum pressure andis a direct consequence of the Heisenberg uncertainty principle (Pitaevskiiand Stringari, 2003). It depends on the gradient of density, suggesting thatquantum effects are important in non-uniform gases and, for uniform systems,close to vortex cores. The system of equations (4)-(5) is equivalent to theoriginal GP equation (1) and will be used in the next section to derive ahydrodynamic analogy.It is important to define the integral quantities conserved by the GPequation (1). First, after multiplying (1) by ψ ∗ , integrating in space andtaking the imaginary part, we obtain that the number of atoms N in thesystem is conserved ( ∂N/∂t = 0), with N = Z | ψ | d x = Z n d x . (6)Second, after multiplying (1) by ∂ψ ∗ /∂t , integrating in space and taking thereal part, we obtain that the energy of the system is also conserved: E ( ψ ) = Z (cid:126) m |∇ ψ | + V trap | ψ | + 12 g | ψ | ! d x . (7)One can also obtain that ∂E/∂t = 0 by directly differentiating (7).Finally, it is useful to introduce stationary solutions of the GP equation.These are obtained by considering that the wave function evolves in time as: ψ ( x , t ) = Ψ( x ) exp( − iµt/ (cid:126) ), (8)7ith µ the chemical potential. Note that | ψ | = | Ψ | and the number of atomsin (6) is conserved by the stationary field. The time-evolution GP equation(1) then reduces to the stationary (time-independent) GP equation: − (cid:126) m ∇ Ψ + V trap Ψ + g | Ψ | Ψ = µ Ψ. (9)The chemical potential µ is fixed by the normalization condition (6) andexpressed from (8) and (9) as: µ = 1 N Z (cid:126) m |∇ Ψ | + V trap | Ψ | + g | Ψ | ! d x = 1 N (cid:18) E (Ψ) + Z g | Ψ | d x (cid:19) .(10) The hydrodynamic analogy of the GP equation (1) is obtained by relatingthe wave function ψ to a superfluid flow of mass density ρ ( x , t ) = m n ( x , t ) = m | ψ ( x , t ) | , (11)and velocity v ( x , t ) = (cid:126) m ∇ θ ( x , t ) = (cid:126) ρ ψ ∗ ∇ ψ − ψ ∇ ψ ∗ i . (12)For a flow of non-vanishing density ( ρ = 0), we infer from (12) that thesuperfluid is irrotational: ∇ × v = 0. (13)The velocity v is usually related to the current density (Pitaevskii andStringari, 2003): j ( x , t ) = n ( x , t ) v ( x , t ) = 1 m ρ ( x , t ) v ( x , t ). (14)Taking the divergence of (14) and using (12) we obtain the expression: ∇ · j = (cid:126) m ψ ∗ ∇ ψ − ψ ∇ ψ ∗ i = (cid:126) m (cid:16) ∇ n · ∇ θ + n ∇ θ (cid:17) = (cid:126) m ∇ · ( n ∇ θ ). (15)We infer from (15) and (4) that the equation for the time evolution of thedensity could be now written as: ∂n∂t + ∇ · j = 0 ⇐⇒ ∂ρ∂t + ∇ · ( ρ v ) = 0. (16)8quation (16) is the continuity equation of the superflow, expressing theconservation of the number of particles N given by (6).To complete the hydrodynamic equations, we take the gradient of theequation (5) and use (12) to obtain the equation for the time evolution of thevelocity: ∂ v ∂t + 12 ∇ ( v ) = − m ∇ ( gn + V trap ) + (cid:126) m ∇ √ ρ ∇ ( √ ρ ) ! . (17)Using the identity 12 ∇ ( v ) = ( v · ∇ ) v + v × ( ∇ × v ),and (13), we infer from (17) that, after neglecting the last (quantum pressure)term and the trapping potential ( V trap = 0): ∂ v ∂t + ( v · ∇ ) v = − ρ ∇ gρ m ! . (18)Equations (16) and (18) are similar to the Euler equations describingthe evolution of a compressible, barotropic and inviscid classical flow withpressure: P = gρ m . (19)Notice that (17) is exactly the non-conservative form for the time evolutionof the momentum ρ v = m j . As a consequence, equations (16) and (18) showthat the total momentum of the superfluid, p = Z ( ρ v ) d x = m Z j d x = (cid:126) Z I mag ( ψ ∗ ∇ ψ ) d x , (20)is also conserved. Note that the momentum conservation can be directlyderived from (1). We recall that the energy (7) is conserved; it can berewritten as: E = Z ρ v + (cid:126) m |∇√ ρ | + ρm V trap + gρ m ! d x . (21)To summarize, the three main integral invariants of the GP equation (1) are:the number of particles N (6), the energy E , in the form (7) or (21), and,when V trap = 0, the momentum p (20).9 .3 The Thomas-Fermi limit The Thomas-Fermi regime is characterized by strong interactions (the kineticenergy is negligible compared to the interaction energy). If the kinetic energy(first term) is neglected in (9), the stationary solution Ψ takes a simple form,that can be expressed from the equation: V trap ( x ) + g | Ψ | = µ . (22)The Thomas-Fermi limit (22) can also be derived from (17) by considering v = 0 and neglecting the quantum pressure term. For a harmonic trappingpotential with V trap = mω x /
2, with a reference trapping frequency ω , theThomas-Fermi regime is attained when N a s /a h (cid:29)
1, with a h = q (cid:126) / ( mω )the harmonic oscillator length. In the absence of a trapping potential, weobtain the Bogoliubov relation µ = gn (Pitaevskii and Stringari, 2003). Since the GP equation (1) sustains wave solutions, it is interesting to estimatethe response of the system to small perturbations. The Bogoliubov-de Gennesmodel is based on the linearisation of (1) assuming that: ψ ( x , t ) = [Ψ( x ) + δψ ( x , t )] exp( − iµt/ (cid:126) ), (23)where Ψ( x ) is a stationary solution satisfying (9), and δψ a perturbationof small amplitude. After discarding the terms of order δψ , we obtain thefollowing linearized equation: µδψ + i (cid:126) ∂δψ∂t = − (cid:126) m ∇ δψ + V trap δψ + g Ψ δψ ∗ + 2 g | Ψ | δψ . (24)We consider the following form of the perturbation: δψ = (cid:16) a ( x ) e − iωt + b ∗ ( x ) e iω ∗ t (cid:17) , (25)where a , b are complex functions (small amplitudes) and ω is a complexfrequency. From (24) we separate the terms in e − iωt and e iω ∗ t and obtain thefollowing Bogoliubov-de Gennes equations: − (cid:126) m ∇ + V trap − µ + 2 g | Ψ | ! a + g Ψ b = (cid:126) ω a , (26) g (Ψ ∗ ) a + − (cid:126) m ∇ + V trap − µ + 2 g | Ψ | ! b = − (cid:126) ω b . (27)10olving the system of differential equations (26)-(27) provides the eigenfre-quencies ω and the modes amplitudes a and b . Note from (25) that if ( a , b , ω )is a solution of this system, then ( b ∗ , a ∗ , − ω ∗ ) is also solution (representingthe same physical oscillation). Moreover, due to the Hamiltonian nature ofthe system (26)-(27), ω ∗ is also an eigenfrequency. Consequently, modes cor-responding to eigenfrequencies with a non-zero imaginary part are amplifiedin time. Such modes are associated to the dynamic instability of the system. Since the focus of this study is the QT in superfluid helium, we set to zero thetrapping potential ( V trap = 0) in the remainder of the paper. In this section,we introduce the main characteristics of the QT flow that have to be resolvedin numerical simulations. We first analyze the constant density uniformsuperflow, which is the background of the QT field. Important characteristicscales for properly designing numerical simulations are introduced: the healinglength and the sound velocity. The response of this uniform flow to smallperturbations is also presented by deriving the Bogoliubov dispersion relation.This relation is an important guide in setting the numerical resolution inspectral methods. It also defines the domain of the validity of the GP modelin describing superfluid helium QT. We then introduce the notion of quantumvortex and its representation on a computational grid. Finally, we presentthe main integral quantities that will be used to analyze the QT field: thecompressible and incompressible kinetic energy and the helicity. An elementary solution to the stationary GP equation (9), in absence of thetrapping potential, is a flow with constant density ρ = ρ . From (9) and (8),we infer that the corresponding wave functions and chemical potential arerelated by µ = g | Ψ | = g | ψ | = gn = g ρ m . (28)The solution Ψ could be taken as real and represents a first approximationof a quantum flow developing in a container of volume V , far from the walls.The number of atoms in the container and the energy of the system follow11rom (6) and (7), respectively: N = | Ψ | V = n V = ρ m V . (29) E = 12 g | Ψ | V = gN V . (30)It is now possible to determine the constitutive relation of the ideal barotropicfluid: P = − ∂E ∂V = gN V = 12 gn = 12 gρ m , (31)which was put into evidence in the momentum equation (17) and also definedin (19). Note that the pressure of the superflow does not vanish at zerotemperature, as in an ideal gas (Pitaevskii and Stringari, 2003). We also inferfrom (31) that the flow is compressible, with the sound velocity c defined asin classical hydrodynamics: c = s ∂P ∂ρ = √ gρ m = r gn m . (32)Combining (28) and (32) we obtain the equality: µ = mc . (33)The sound velocity gives a characteristic velocity of the system. To havea complete space-time description of the system we need to introduce acharacteristic length scale. The healing length indicates the distance overwhich density variations take place in the system. It can de derived in severalways. If in the stationary GP equation (9) we assume that the first (kinetic-energy) term balances the non-linear interaction term over the length ξ , wecan use the following approximation in the stationary limit: (cid:126) m Ψ ξ ≈ g | Ψ | Ψ = ⇒ ξ = (cid:126) q mg | Ψ | . (34)For a constant density flow, using (28) and (33), the expression of the healinglength becomes: ξ = (cid:126) √ mgn = (cid:126) √ mµ = 1 √ (cid:126) mc . (35)12he same expression for the healing length could be obtained from themomentum equation (17) by imposing that the pressure term balances thequantum pressure term over the healing length. It results that the quantumpressure is negligible for distances R (cid:29) ξ , which is exactly the domain ofvalidity of the hydrodynamic analogy with the Euler equations. The Bogoliubov-de Gennes linearized system derived in §2.4 is greatly sim-plified in the case of the flow of uniform density described by (28), with realwave function Ψ = √ n . The perturbations are then taken as plane waves: a ( x ) = ue i k · x , b ( x ) = ve i k · x , (36)with k the wave number vector. The BdG system (26)-(27) becomes for thiscase: (cid:126) m k + gn − (cid:126) ω ! u + ( gn ) v = 0, (37)( gn ) u + (cid:126) m k + gn + (cid:126) ω ! v = 0, (38)with a non-trivial solution if( (cid:126) ω ) = (cid:126) m k ! + ( gn ) (cid:126) m k . (39)Using (32) to express the sound velocity and (35) for the healing length, theBogoliubov dispersion relation (39) becomes: ω = ( ck ) s ξ k kξ (cid:28)
1) and the excitations in thisregime are called phonons (sound waves). Going back to the momentumequation (17), we can easily see that in the phonons regime, the quantumpressure is negligible in front of the hydrodynamic pressure. Consequently,the validity of the hydrodynamic analogy is limited to the phonons excitations,a regime where the quantum pressure could be neglected. The transition13etween the phonon regime and the particle regime takes place at ( kξ ∼ kξ (cid:29) ω → cξk / √ ω → (cid:126) k / m ). The information on the separation ofdifferent regimes is quite useful for numerical simulations, since the quantity( kξ ) must be carefully assessed to represent correctly the density waves onthe the computational grid. Note that in helium II, due to strong interactionsbetween particles, the dispersion relation has a different shape, with a linearregime followed by a quadratic regime with a maxon (local maximum) anda roton (local minimum) (Barenghi and Parker, 2016). The excitations inthe quadratic region near the minimum of the dispersion curve are calledrotons. Consequently, using the GP equation allows us to capture only thephonons regime of excitations in a quantum flow. Note that GP equationcan be modified by including non-local terms to model any dispersion curveBerloff et al. (2014). The definition (12) of the superflow velocity becomes singular along lines withvanishing density ( ρ = 0). The lines along which both real and imaginary partthe order parameter are zero define topological defects, known as quantumvortices . The hydrodynamic analogy through the Madelung transformationbecomes singular when vortices are present in the superflow. A detailed reviewof mathematical problems related to the Madelung transformation in presenceof quantum vortices is offered in Carles et al. (2012). Another intriguingconsequence of the hydrodynamic analogy appears when calculating thecirculation Γ along a regular path C surrounding a simply connected domain S : Γ = I C v · d l . (41)From (13) we infer by applying the Stokes theorem that Γ = 0 if S is a simplyconnected domain without vortices. When a vortex line crosses the domain S , the curl of the velocity expressed by (13) becomes a Dirac function andthe circulation Γ takes non zero values. It is important to note, however, thatvortex solutions are not singular solutions of the GP equation (1). Indeed, astraight-line vortex solution Ψ v is represented using cylindrical coordinates( r , ϕ , z ) by: Ψ v = | Ψ v ( r ) | e iκϕ , (42)14here κ is necessarily an integer to ensure that the wave function is singlevalued. Using (12) we infer that the velocity around the vortex line istangential and singular for r = 0: v v = (cid:126) m r ∂ ( κϕ ) ∂ϕ e ϕ = (cid:126) m κr e ϕ . (43)The circulation around this vortex line follows from (41):Γ v = (2 π ) κ (cid:126) m = κ hm . (44)This quantification of the vortex line circulation is the outstanding differencebetween quantum and classical hydrodynamics. The integer κ is usuallyreferred to as winding number or charge of the vortex.The vortex solution (42) can be further developed by assuming that: | Ψ v ( r ) | = √ n f ( η ), η = rξ , (45)with n the far-field background solution with constant density and ξ thehealing length. Injecting (42) with (45) in the stationary GP equation (9) andusing definitions (28) and (35), we obtain the following ordinary differentialequation for f : − η ddη η dfdη ! + f κ η − ! + f = 0, (46)with limit conditions: f → η → f → η → ∞ . Theasymptotic behaviour of this solution near the origin ( r = 0) is well-known(Neu, 1990): f ( η ) ∼ η | κ | + O ( η | κ | +2 ), η →
0, (47)suggesting that the vortex core, i. e. the region near the vortex line where thedensity is varying in a significant way, is of the order of the healing length ξ .Using (43), the kinetic energy ( R ρ v v d x ) of the vortex solution, which isthe main contribution to the total energy (21), results to be proportional to κ ( e. g. Barenghi and Parker, 2016). This implies that a multiply quantizedvortex with κ > κ -singly quantizedvortices in GP-QT. As shown in Takeuchi et al. (2018), the multiply quantizedvortex is also dynamically unstable, since the complex frequency ω obtainedin the Bogoliubov-de Gennes equations (26)-(27) has a non-zero imaginarypart. 15 .4 Energy decomposition The energy is the first integral quantity that will be used to characterizethe QT flow field. The accuracy of numerical simulations in conserving thisquantity is an important checkpoint to validate the numerical scheme and thegrid resolution. As in CT, energy spectra will be used to identify different(Kolmogorov) regimes/ranges in the structure of the turbulent field. Startingfrom the observation that the QT field can be viewed as a background uniformflow to which a large number of quantum vortices are superimposed, the totalenergy of the system in QT studies ( e. g.
Nore et al., 1997a,b) is generallycomputed using the form: E T ( ψ ) = Z (cid:126) m |∇ ψ | + 12 g (cid:16) | ψ | − | Ψ | (cid:17) ! d x , (48)where | Ψ | = n is the atomic density of the uniform flow. This expression isstrictly equivalent to the form (7) of the energy because of the conservationof the number of atoms (6). The corresponding GP equation, equivalent to(1) is then: i (cid:126) ∂∂t ψ ( x , t ) = δE T δψ ∗ = − (cid:126) m ∇ + g (cid:16) | ψ ( x , t ) | − | Ψ | (cid:17)! ψ ( x , t ), (49)Using the hydrodynamic analogy presented in §2.2, the total energy (48) canbe also presented as: E T ( ρ , v ) = Z ρ v + (cid:126) m |∇√ ρ | + 12 m g ( ρ − ρ ) ! d x . (50)Again (50) is strictly equivalent to (21). The three terms is (50) correspondto (Nore et al., 1997a,b):– the kinetic energy E kin = Z |√ ρ v | d x , (51)– the so-called quantum energy (expressed using (35)) E q = Z (cid:126) m |∇√ ρ | d x = Z c ξ |∇√ ρ | d x . (52)– and the internal energy (expressed using (32)): E int = Z m g ( ρ − ρ ) d x = Z c ( ρ − ρ ) ρ d x . (53)16he kinetic energy E kin can be further decomposed (Nore et al., 1997a,b) asthe sum of a compressible part E ckin and an incompressible part E ikin : E ckin = Z | ( √ ρ v ) c | d x , E ikin = Z | ( √ ρ v ) i | d x , (54)owing to the Helmholtz decomposition:( √ ρ v ) = ( √ ρ v ) c + ( √ ρ v ) i , with ∇ × ( √ ρ v ) c = 0, and ∇ · ( √ ρ v ) i = 0.(55) Spectra of the different components of the energy and structure functions ofvelocity will be used to analyze the QT field, as in CT. The energy spectra arecomputed using the following expressions resulting after applying Parseval’stheorem for the Fourier transform: E i kin ( k ) = 12(2 π ) Z | k | = k |F k ( √ ρ v ) i | d Ω k , E c kin ( k ) = 12(2 π ) Z | k | = k |F k ( √ ρ v ) c | d Ω k , E int ( k ) = c ρ (2 π ) Z | k | = k |F k ( ρ − ρ ) | d Ω k , E q ( k ) = c ξ (2 π ) Z | k | = k |F k ( ∇√ ρ ) | d Ω k , (56)where F k is the Fourier transform F k ( f ( x )) = Z f ( x ) e − i k · x d x , F − x ( g ( k )) = 1(2 π ) Z g ( k ) e i k · x d k , (57)and Ω k is the solid angle in the spectral space. Note that the components ofthe energy (51)-(53) can also be represented and computed in the spectralspace following: E kin = 12(2 π ) Z |F k ( √ ρ v ) | d k = Z E kin ( k ) dk , E int = c ρ (2 π ) Z |F k ( ρ − ρ ) | d k = Z E int ( k ) dk , E q = c ξ (2 π ) Z |F k ( ∇√ ρ ) | d k = Z E q ( k ) dk . (58)17he structure function for the velocity following the x -direction (withunitary vector e x ) is defined as: S p// ( r ) = Z (( v ( x + r e x ) − v ( x )) · e x ) p , (59)where p is the order of the structure function and r the length scale. Similarexpressions are used for the structure functions following the y and z directions.Assuming a homogeneous and isotropic distribution of the QT velocity fieldstatistics, averaging over different directions should give the same results. Asa verification, for p = 2 and large length scale r , the structure function couldbe reasonably approximated by:lim r →∞ S // ( r ) ’ Z | v ( x ) · e x | = 2 Z v x . (60) The helicity is another important integral quantity characterizing the QTflow field. The definition of helicity in a classical flow is: H = Z v · ω d x , (61)where ω = ∇× v is the vorticity. In a quantum fluid, the vorticity concentratesin vortex cores as ω ( r ) = hm Z d r ds δ ( r − r ( s )) ds , (62)where r ( s ) denotes the position of the vortex line, δ is the Dirac delta functionand s the arclength. Therefore, only quantized vortices bring a non-zerocontribution to the helicity. We consider in the following the case in which allvortices form closed loops. Then, the formalism described by (61) and (62) istopologically equivalent to (Moffat, 1969): H = hm ! X i = j Lk ij + X i SL i , (63)where Lk ij ∈ Z denotes the Gauss linking number Lk ij = 14 π Z C i Z C j ( X i − X j ) | X i − X j | · ( d X i × d X j ), (64)18etween two vortex loops i and j . X i and X j denote the position vectorsof points on the centerlines C i and C j ( i = j ) for vortex loops i and j ,respectively. SL i ∈ Z denotes the C˘alug˘areanu-White self-linking number ofeach single vortex loop. SL i can be also written as the sum SL i = W r i + T w i of the writhe W r i ∈ R and the total twist T w i ∈ R . The writhe is defined by W r i = 14 π Z C i Z C i ( X i − Y i ) | X i − Y i | · ( d X i × d Y i ), (65)where X i and Y i denote two distinct points on the same curve C i . The totaltwist T w i can be further decomposed in terms of normalized total torsion T i ∈ R given by T i = 12 π Z C i τ ( s ) ds , (66)and intrinsic twist N i ∈ Z of the ribbon R ( C i , C ∗ i ), where τ ( s ) is the torsionof C i . The second curve C ∗ i of the ribbon R ( C i , C ∗ i ) is fixed in order forthe ribbon to be identified by the points of constant phase that lie on the (cid:15) -portion of the iso-phase surface.Within the definition (63), the helicity always vanishes due to the identity(Laing et al., 2015): X i = j Lk ij + X i W r i = − X i T w i . (67)As a consequence, the centerline helicity H = hm ! X i = j Lk ij + X i W r i = h πm X i , j Z C i Z C j ( X i − X j ) | X i − X j | · ( d X i × d X j ), (68)has also been introduced as a nonzero helicity in quantum fluid. Note thatthe double integral in Eq. (68) means Eq. (64) for i = j and Eq. (65)for i = j . Although the centerline helicity is just the quantification of howknotted vortex lines are, its calculation requires detailed extraction of allcenterlines of the quantized vortices. This can be done by using, for example,the interpolation of the computed wave function at sub-grid scales (Scheeleret al., 2014). 19 new method which yields the same results as the centerline helicity wasintroduced in Clark di Leoni et al. (2016). Because the superfluid velocity v diverges at the vortex cores as shown in Eq. (43), the direct calculation of thehelicity looks ill-defined. However, only the superfluid velocity perpendicularto the quantized vortex has a singularity, while the component parallel tothe quantized vortex remains regular. Using this observation, a regularizedvelocity is defined as: v reg = v k w / √ w j w j , (69)where w = (cid:126) im ∇ ψ ∗ × ∇ ψ , (70)and v k = (cid:126) w i [( ∂ i ∂ j ψ ) ∂ j ψ ∗ − ( ∂ i ∂ j ψ ∗ ) ∂ j ψ ]2 im √ w k w k ( ∂ l ψ )( ∂ l ψ ∗ ) . (71)The resultant regularized helicity H reg = Z v reg · ω d x , (72)can be well defined as another formalism for the centerline helicity. Thisexpression was proven useful and efficient in computing the helicity of flowswith hundreds of thousands of knots (Clark di Leoni et al., 2016). Numerical simulations were performed using the parallel code called GPS(Gross-Pitaevskii Simulator) (Parnaudeau et al., 2015). The code is based on aFourier-spectral space discretization and recent up-to-date numerical methods:a semi-implicit backward-Euler scheme with Krylov preconditioning for thestationary GP equation (Antoine and Duboscq, 2014) and various schemes(Strang splitting, relaxation, Crank-Nicolson) for the real-time GP equation(Antoine et al., 2013). GPS is written in Fortran 90 and uses a two-levelcommunication scheme based on MPI across nodes and OpenMP within nodes.Only one external library, FFTW (Frigo and Johnson, 2005), is required forthe computation. Initially designed to simulate BEC configurations (with orwithout rotation), the GPS code was adapted in this study for the simulationof QT flows. We present in this section the main features of the numericalsystem: the particular scaling used to obtain the GP dimensionless equations,20nd the particular numerical methods used to prepare the initial state andthen to advance in real-time the GP wave function.
For the numerical resolution of the GP equation (1), it is convenient to usea dimensionless form obtained after scaling all physical quantities with thecharacteristic scales of the QT field introduced in §3. We present below ageneral formalism that covers the various forms of scaling used in the physicaland mathematical GP literature. We start by considering general referencescales ( L ref , v ref ) for length and velocity, respectively. A natural scale for thewave function ψ is ψ ref = √ n . With the scaling:˜ x = x L ref , ˜ t = v ref L ref t , ˜ ψ = ψ √ n , (73)the dimensionless GP equation (49) (with V trap = 0) becomes: i ∂∂ ˜ t ˜ ψ ( x , t ) = (cid:16) − α ˜ ∇ + β (cid:16) | ˜ ψ (˜ x , t ) | − (cid:17)(cid:17) ˜ ψ ( x , t ), (74)with non-dimensional coefficients: α = (cid:126) m L ref v ref = √ ξc L ref v ref = 1 √ ξL ref ! (cid:18) cv ref (cid:19) , (75) β = gn (cid:126) L ref v ref = mc (cid:126) L ref v ref = c √ ξc L ref v ref = 1 √ L ref ξ ! (cid:18) cv ref (cid:19) . (76)From (75)-(76) we infer that non-dimensional coefficients α and β are relatedto physically relevant scales through:˜ ξ = ξL ref = s αβ , ˜ c = cv ref = q αβ = 1 M ref , (77)where ( ξ/L ref ) represents the non-dimensional healing length and ( c/v ref ) thenon-dimensional sound velocity. M ref is the reference Mach number, definedas the ratio between the reference velocity and the sound velocity.The last important parameters to define when working with non-dimensionalequations are the size of the computational box and the grid resolution. Ifthe physical GP equation (1) is defined in a cubic computational domain of21hysical size L , the non-dimensional size L of the computational box used todiscretize the non-dimensional equation (74) is then: L = LL ref = Lξ ! ξL ref ! = Lξ ! s αβ . (78)We recall that ξ is a good approximation of the radius of a quantum vortex(see §3.3). It follows that the ratio ( L/ξ ) in (78) is physically important sinceit indicates how many vortices the computational domain can accommodatein one direction: N dv = L ξ = L s βα . (79)Thus, increasing the value of L (for fixed α and β ) will result in a highernumber of vortices present in the computational box.When defining the grid resolution, it is important to control the numberof grid points inside the vortex core. If the numerical simulation uses N x grid points in each direction, the physical grid spacing is δx = L/N x , or innon-dimensional units δ ˜ x = ( δx/L ref ) = L /N x . It is important to quantifythe grid spacing with respect to the healing length by defining: χ = δxξ = (cid:18) L N x (cid:19) L ref ξ ! = (cid:18) L N x (cid:19) s βα . (80)The parameter χ defined in (80) is also important when analyzing the dis-persion relation (40) presented in §3.2 to assess on the validity of the GPmodel. Indeed, the maximum wave-number represented on a grid of size δx is k max = (2 π ) / (2 δx ) and, consequently, the non-dimensional quantity ( k max ξ )is expressed as: k max ξ = π ξδx = πχ . (81)The numerical resolution will be then fixed in order to keep ( k max ξ ) ≈ β = 1. The value of α is fixed to 1 or 1 /
2, as commonlyused in the classical GP equation. The reference scales result then from2277) and the size of the computational domain from (78). For the firstchoice we obtain: β = α = 1 = ⇒ ˜ ξ = ξL ref = 1, v ref = c √ M ref = 1 √ Lξ = L .(82)In this case, a quantum vortex is of size 1, the velocity of sound is √ L is a large integer (to resolve alarge number of vortices). A second common choice is: β = 1, α = 12 = ⇒ ˜ ξ = ξL ref = 1 √ v ref = c , M ref = 1, Lξ = √ L .(83)In this case, the velocity of sound is 1 and the size of the vortex oforder of (1 / √ ξ , the size of the vortex in non-dimensionalunits. In exchange, the size of the domain L has to be adapted infunction of the resolution N x , in order to keep constant the parameter χ in (80) and, implicitly, ( k max ξ ) in (81).• A second approach (Nore et al., 1997a,b) is used in the present paper.The size of the non-dimensional computational box is first set to L = 2 π ,which is convenient for spectral methods. Moreover, instead of settingindependently the constants α and β , only the value of the referenceMach number M ref is fixed to a relatively low value. This is equivalentto impose the value of the product αβ . From previous relations we inferthat: L = 2 π , αβ = 12 M = ⇒ ξL ref = √ αM ref , v ref = M ref c , Lξ = 2 π √ αM ref , k max ξ = (cid:18) N x (cid:19) √ αM ref . (84)We note from (84) that the parameter α can be used to control the non-dimensional size of the vortex, while the grid resolution N x can be setto control the parameter ( k max ξ ). We generally set for QT simulations M ref = 0.5, equivalent to αβ = 2.23articular care has to be devoted when computing non-dimensional valuesof different quantities appearing in integral invariants or in the hydrodynamicanalogy. If the non-dimensional wave function is computed from (74) as˜ ψ (˜ x , ˜ t ) = | ˜ ψ | exp( iθ (˜ x , ˜ t )), the scaled number of atoms results from (6) and(29): ˜ N = NN = 1 L Z D | ˜ ψ | d ˜ x , (85)where D is the non-dimensional computation domain. The scaled total energy(per volume unit) results from (48):˜ E ( ˜ ψ ) = E T ( ψ ) E ref = 2 α L Z D α | ˜ ∇ ˜ ψ | + β (cid:16) | ˜ ψ | − (cid:17) ! d ˜ x , (86)with energy units E ref = ρ v L = (cid:126) n L L / (4 mα ).For the hydrodynamic analogy developed in §2.2, taking as reference thedensity of the background uniform flow, i. e. ρ ref = ρ = mn , results in:˜ ρ = ρρ ref = mn | ˜ ψ | ρ ref = | ˜ ψ | . (87)The momentum is derived from (12) and thus computed in the non-dimensionalcode as: ˜ ρ ˜ v (˜ x , ˜ t ) = 2 α ˜ ψ ∗ ˜ ∇ ˜ ψ − ˜ ψ ˜ ∇ ˜ ψ ∗ i = (2 α ) I mag ( ˜ ψ ∗ ˜ ∇ ˜ ψ ). (88)The non-dimensional superflow velocity also results from (12):˜ v (˜ x , ˜ t ) = v ( x , t ) v ref = (cid:126) mv ref L ref ˜ ∇ θ (˜ x , ˜ t ) = 2 α ˜ ∇ θ (˜ x , ˜ t ), (89)and the non-dimensional circulation of a vortex of winding number ( κ = 1)from (44): ˜Γ = Γ v v ref L ref = 2 π (cid:126) mv ref L ref = 4 πα . (90)Finally, the hydrodynamic expression (50) of the total energy becomes˜ E ( ˜ ψ ) = E T ( ψ ) E ref = 1 L Z D (cid:18)
12 ˜ ρ ˜ v + (2 α ) |∇ ( √ ˜ ρ ) | + ( αβ ) ( ˜ ρ − (cid:19) d ˜ x , (91)with the same reference energy as in (86) E ref = ρ v L . In (91) the firstterm represents the kinetic energy ˜ E kin ( ˜ ψ ), the second the quantum energy24 E q ( ˜ ψ ) and the third the interaction energy ˜ E int ( ˜ ψ ). Note that 2 α = ˜ ξ ˜ c and αβ = 1 / (2 M ).To simplify the presentation, we drop in the following the tilde notation.All the developments and results in the remaining of the paper concernnon-dimensional quantities. To find stationary solutions to (74), a very popular numerical method is the normalized gradient flow (Bao and Du, 2004). The idea is to propagate thewave function following the gradient flow corresponding to the minimizationof the energy (86). In the original method, the solution is subsequentlynormalized to satisfy the constraint of the conservation of the number ofatoms (equivalent to imposing the L -norm of the solution). In our case, wewant to find a stationary state that mimics a classical flow with prescribedvelocity v ext . Assuming that ∇ · v ext = 0, after applying a local Galileantransformation, the non-dimensional energy of the driven field becomes (seeNore et al., 1997a, for details): E v = 2 α L Z D α (cid:12)(cid:12)(cid:12)(cid:12) ∇ ψ − i v ext α ψ (cid:12)(cid:12)(cid:12)(cid:12) + β | ψ | − ! d x , (92)or, using the hydrodynamic analogy: E v ( ψ ) = 1 L Z D (cid:18) ρ | v − v ext | + (2 α ) |∇ ( √ ρ ) | + ( αβ ) ( ρ − (cid:19) d x . (93)In this setting, we are searching a unconstrained minimizer of E v . Owing to theprevious decomposition, there is a competition between the background uni-form distribution | ψ | = 1 and a phase accommodating to v ext . Numerically,we solve the gradient descent equation (or Advective Real Ginzburg-LandauEquation, ARGLE): ∂∂τ φ ( x , t ) = α ∇ − i v ext · ∇ − | v ext | α + β − β | φ ( x , t ) | ! φ ( x , t ), x ∈ D ,(94)with initial condition φ ( x , 0 + ) = φ ( x ). Note that τ is here a pseudo-timeused to propagate the solution until a stationary state is reached. Hence, thismethod belongs to the class of so-called imaginary time propagation methods.25e use a semi-implicit Backward Euler scheme to advance the solution in thepseudo-time interval ( τ n , τ n +1 ):˜ φ n +1 ( x ) − φ n ( x ) δτ = α ∇ ˜ φ n +1 + α ∇ − i v ext · ∇ − | v ext | α + β − β | φ n ( x ) | ! ˜ φ n ( x ).(95)where δτ := τ n +1 − τ n . The resulting system is solved with spectral accuracyusing FFTs.The ARGLE procedure stops either after a pre-definite number of pseudo-time steps or when the convergence criterion is reached: || φ n +1 − φ n || ∞ δτ ≤ (cid:15) , (96)where (cid:15) is a user defined parameter. If this criterion is not satisfied at theend of the computation, it is still possible to check the energy convergencecondition: | E v ( φ n +1 ) − E v ( φ n ) | δτ E v ( φ n ) ≤ (cid:15) , (97)which is generally less constraining than (96). Note that, even when theconvergence is achieved, we can only guarantee that the Backward Eulermethod provides a local minimum of E v . The simulation of QT consists of solving the GP equation (74) using apseudo-spectral scheme in space and a second order splitting for the timediscretization (ADI, Alternating Direction Implicit or Strang splitting). Letus rewrite (74) as: ∂∂t ψ = iα ∆ ψ − iβ | ψ | ψ = L x ψ + L y ψ + L z ψ + N ( ψ ), (98)with the following definitions: L x ψ = iα∂ xx ψ , L y ψ = iα∂ yy ψ , L z ψ = iα∂ zz ψ , N ( ψ ) = − iβ | ψ | ψ . (99)26f H denotes one of the previous operators ( L x , L y , L z or N ), and φ is agiven field, we denote by S ( s , H ) φ := ψ ( s ) the solution at time t = s of thefollowing Cauchy problem: ∂∂t ψ = H ( ψ ), ψ ( t = 0) = φ . (100)Then the second order ADI time scheme could be presented as: ψ n +1 = S δt L x ! S δt L y ! S δt L z ! S ( δt , N ) S δt L z ! S δt L y ! S δt L x ! ψ n .(101)This scheme is second order accurate, provided that we solve the partialproblems exactly, i. e. each term S ( s , H ) is computed exactly. This isachieved using the spectral representation. For j ∈ { x , y , z } , using F j theFourier transform in the j direction, we obtain S ( s , L j ) φ = F − j (cid:16) e − iαk j s F j φ (cid:17) . (102)For the non linear operator S ( δt , N ), we notice that, if ψ N is such that ∂ t ψ N = − iβ | ψ N | ψ N , then ∂ ( | ψ N | ) /∂t = 0. Consequently, this step can besolved analytically and: S ( s , N ) φ = e − iβ | φ | s φ . (103)In conclusion, using the spectral discretization we obtain a second orderaccurate scheme for the time integration. As in numerical studies of classical turbulence, the preparation of the initialstate is crucial in investigating statistical properties of QT. We describe in thissection four different approaches to generate the initial field for the simulationof decaying GP-QT. Each method is associated to a benchmark for the GP-QTsimulation. The first two methods are classical (Nore et al., 1997a,b) andinspired from CT. They start from defining a velocity field containing vortices.The Taylor-Green or the Arnold-Beltrami-Childress (ABC) model flows areused for this step. A wave function field is then constructed such that its27odal lines correspond to vortex lines of the velocity field. This initial wavefunction is then used in the ARGLE procedure described in §4.2 to generatean initial field for the real-time GP simulations. The role of the ARGLE stepis to reduce the acoustic emission of the initial field. The last two methods arenew and based on the direct manipulation of the wave function. We prescribeeither a random phase field or we manufacture an initial field containing manyquantum vortex rings. The four methods are described in detail below.
The velocity v TG of the Taylor-Green (TG) three-dimensional vortex flow isdefined as: v TG, x ( x , y , z ) = sin( x ) cos( y ) cos( z ), v TG, y ( x , y , z ) = − cos( x ) sin( y ) cos( z ), v TG, z ( x , y , z ) = 0. (104)To create a wave function field ψ TG with zeros along vortex lines of v TG ,we make use of the Clebsch representation of the velocity field (Nore et al.,1997a,b): ∇ × v TG = ∇ λ × ∇ µ , (105)with Clebsch potentials λ ( x , y , z ) = cos( x ) q | cos( z ) | , µ ( x , y , z ) = cos( y ) q | cos( z ) | sgn(cos( z )), (106)where sgn is the sign function. Note that a zero in the ( λ , µ ) plane correspondsto a vortex line of v TG (see Nore et al. (1997a,b) for details).In practice, we start by defining in the ( λ , µ ) plane a complex field ψ e with a simple zero at the origin: ψ e ( λ , µ ) = ( λ + iµ ) tanh( √ λ + µ / √ ξ ) √ λ + µ . (107)When replacing (106) into (107), a three-dimensional complex field is obtained,with one nodal line. We can further define on [0, π ] : ψ ( x , y , z ) = ψ ( λ ( x , y , z ), µ ( x , y , z )) = ψ e ( λ − √ µ ) ψ e ( λ , µ − √ × ψ e ( λ + 1 √ µ ) ψ e ( λ , µ + 1 √ ψ e is extendedby mirror reflection to the entire domain [0, 2 π ] , the obtained wave functionfield contains closed rings inside the domain (see Fig. 1 a, right).The last manipulation of the wave function is intended to match thecirculation of the velocity field v TG . From (105) and (106) we compute thecirculation on the face z = 0, ( x , y ) ∈ [0, π ] × [0, π ] using the Stokes’ theorem:Γ z =0 = Z π Z π ( ∇ × v TG ) · e z dx dy = Z π Z π x ) sin( y ) dx dy = 8. (109)Defining the ratio of the total circulation to the circulation (90) of a single vor-tex as γ d = Γ / Γ v = 2 / ( πα ), the wave-function field matching the circulationof the TG velocity field is (Nore et al., 1997a) is ψ ARGLE ( x , y , z ) = ψ ( λ ( x , y , z ), µ ( x , y , z )) [ γ d / , (110)where [.] denotes the integer part. In this setting, each vortex line correspondsto a multiple zero line (see Fig. 1 a). The next step in the preparation ofthe initial field is to use the ARGLE imaginary time procedure (94) with v ext = v TG and initial condition φ ( t = 0 + ) = ψ ARGLE . During the ARGLEdynamics the multiple zero lines in ψ ARGLE will spontaneously split into[ γ d /
4] = [1 / (2 πα )] single zero lines (see Fig. 1 b). The system will finallyconverge to initial conditions for the GPE, compatible with the TG flow, andwith minimal sound emission. We denote the resulting converged state as φ TG (see Fig. 1 c). 29)b)c)Figure 1: Illustration of the initial field preparation using the Taylor-Greenvortex flow. Imaginary time evolution of quantized vortices (iso-surfaces oflow φ TG ) during the ARGLE calculation. Case N x = 128 with [ γ d /
4] = 3(see Table 2). Two different views, on the subdomain [0, π ] (left) and theentire domain [0, 2 π ] (right), illustrating the symmetry of the flow. Panelsfrom top to bottom: (a) τ = 0, the initial condition φ ( t = 0 + ) = ψ ARGLE , Eq.(110) with multiply quantized (thick) vortices, (b) τ = 1 when each initialvortex line splits in 3 singly quantized vortices and (c) τ = 60 for the finalconverged ARGLE field, with closed loops inside the domain.30 .2 Arnold-Beltrami-Childress (ABC) flow The TG vortex flow (104) has zero helicity. To obtain a a helical flow atlarge scales, we use the method suggested by Clark di Leoni et al. (2017) intheir study of helical quantum turbulence at zero temperature. The externalvelocity field is defined as the superposition of two Arnold-Beltrami-Childress(ABC) flows: v ABC = v (1)ABC + v (2)ABC , (111)with v ( k )ABC, x ( x , y , z ) = B cos( ky ) + C sin( kz ), v ( k )ABC, y ( x , y , z ) = C cos( kz ) + A sin( kx ), v ( k )ABC, z ( x , y , z ) = A cos( kx ) + B sin( ky ). (112)Unless stated otherwise, we use ( A , B , C ) = (0.9, 1, 1.1) / √
3. As for theTaylor-Green flow, we use the ARGLE procedure (94) with v ext = v ABC andinitial condition: φ ( t = 0 + ) = ψ (1)ABC × ψ (2)ABC . (113)The wave functions ψ ABC are defined as: ψ ( k )ABC = ψ x , y , zA , k × ψ y , z , xB , k × ψ z , x , yC , k , ψ x , y , zA , k = exp i " A sin( kx )2 α y + i " A cos( kx )2 α z ! ,(114)where [ a ] stands for the nearest integer to a . The ARGLE procedure hasthe role to minimize the amount of energy of acoustic modes in the initialcondition. Details of the quantum ABC flow are discussed by Clark di Leoniet al. (2016, 2017). We denote the resulting converged state as φ ABC . Previous initial fields for the simulation of the QT were built based on theanalogy with classical flows (TG and ABC) with vortices. We present inthis section the first method to set an initial field by direct manipulationof the wave function. A smoothed random phase (SRP) is assigned to theinitial wave function ψ SRP . Initially, there are no vortices present in the field.Vortices nucleate during the time evolution and their interaction generate a31T field. In practice, to obtain the nucleation of enough vortices for QT, weinitialize the field as follows: ψ SRP = e iθ ( x ) , (115)where θ is a smooth random periodic function in the computational box.To create this initial phase, we first generate the random phase θ i , j , k ∈ [ − K , K ] at N r points x i , j , k = N s × ( i , j , k ) where N s = N/N r and i , j , k ∈{
0, 1, 2, · · · , N r − } . Then, θ is obtained by cubic spline interpolation (withperiodicity) using the points ( x i , j , k , θ i , j , k ). The one-dimensional cubic (anduniform) spline interpolation is expressed as: θ N s i r + i s , N s j r , N s k r = A i s θ N s i r , N s j r , N s k r + B i s θ N s ( i r +1), N s j r , N s k r + C i s θ N s i r , N s j r , N s k r + D i s θ N s ( i r +1), N s j r , N s k r , A i s = N s − i s N s , B i s = i s N s , C i s = ( A i s − A i s ) N s D i s = ( B i s − B i s ) N s i r , j r , k r ∈ {
0, 1, · · · , N r − } and i s ∈ {
1, 2, · · · , N s − } . The secondderivative θ is obtained by solving the following linear system with tridiagonalmatrix: N s ( θ N s i r − N s j r , N s k r + 4 θ N s i r , N s j r , N s k r + θ N s i r +1, N s j r , N s k r )= 6 ( θ N s i r − N s j r , N s k r − θ N s i r , N s j r , N s k r + θ N s i r +1, N s j r , N s k r ) . (117)After the interpolation along the i -direction, we compute the spline interpola-tion along the j -direction θ i , N s j r + j s , N s k r = A j s θ i , N s j r , N s k r + B j s θ i , N s ( j r +1), N s k r + C j s θ i , N s j r , N s k r + D j s θ i , N s ( j r +1), N s k r , (118)for i ∈ {
0, 1, · · · , N − } , j r , k r ∈ {
0, 1, · · · , N r − } , and j s ∈ {
1, 2, · · · , N s − } ,and, finally, that along the k -direction θ i , j , N s k r + k s = A k s θ i , j , N s k r + B k s θ i , j , N s ( k r +1) + C k s θ i , j , N s k r + D k s θ i , j , N s ( k r +1) ,(119)for i , j ∈ {
0, 1, · · · , N − } , k r ∈ {
0, 1, · · · , N r − } , and k s ∈ {
1, 2, · · · , N s − } .With this method, the characteristic variation of the phase θ is KN r /π .The characteristic velocity results from (89): v = 2 α ( KN r /π ). The Machnumber of the system is computed using (77) as M = v/c = √ αKN r /π √ β .We denote the resulting converged state as ψ SRP . An example of theresulting flow is shown in Fig. 2. 32 KK K − K Figure 2: Illustration of the initial field preparation using the SRP (smoothedrandom phase) method. Spline interpolation in one dimension using randomvalues for the phase (left) and density contours (right) of the final 3D wavefunction ψ SRP . The main idea for this last initial condition is to prepare an initial statecontaining enough vortices to lead to QT. We derive in the following amethod to fill the computational bow with vortex rings. The challenge is toobtain a physically acceptable ansatz. We start from the single vortex ringsolution to the GP equation (Pitaevskii and Stringari, 2003). A vortex ring33f radius R and constant translational speed can be approximated as: ψ VR ( x , y , z , R ) = f (cid:18)q ( r − R ) + ˜ z (cid:19) e ± i tan − ( ˜ zr − R ), (120) f ( r ) = vuut a ( r/ξ ) + a ( r/ξ ) b ( r/ξ ) + a ( r/ξ ) , (121) a = 73 + 3 √ a = 6 + √ b = 21 + √ x = x − π , (123) r = q ˜ x + ˜ y , (124) ξ = q α/β (125)where f ( r ) is the solution to the GP equation (74) written in cylindricalcoordinates for ψ = f ( r ) e iκ tan − ( y/x ) with κ = 1: − αr ddr r dfdr ! + ακ fr + β ( f − f = 0. (126)The form in (121) is obtained as the Pad´e approximation of this solution.Coefficients a , a , and b in (122) are fixed by satisfying (126) to the orderof ( r/ξ ) for both r/ξ (cid:28) r/ξ (cid:29)
1. This expression stands for a vortexring centered in the origin. Note that this definition is consistent with avortex core size of the order of ξ .The vortex ring ansatz ψ R has the finite net momentum j (see Eq. 14).To eliminate this momentum, we add an opposite-symmetrical ring by settingthe wave function for a vortex-ring pair (VRP) as: ψ VRP ( x , y , z , R , d ) = ψ VR ( x , y , z − d/ R ) ψ ∗ VR ( x , y , z + d/ R ), (127)where d is the inter-vortex distance. Because the ansatz ψ VRP for a vortex-ringpair does not satisfy the periodic boundary condition, we rewrite it as ψ VRP ( x , y , z , R , d ) → ψ VRP ( x , y , z , R , d ) × ψ ∗ VRP (2 L − x , y , z , R , d ) ψ ∗ VRP ( − x , y , z , R , d ) × ψ ∗ VRP ( x , 2 L − y , z , R , d ) ψ ∗ VRP ( x , − y , z , R , d ) × ψ ∗ VRP ( x , y , 2 L − z , R , d ) ψ ∗ VRP ( x , y , − z , R , d ). (128)34he last step to prepare the initial state ψ RVR (random vortex rings) isobtained by randomly putting vortex-ring pairs in the domain. First, werandomly translate the ansatz ψ VRP (128) as ψ RVR ( x , y , z , R , d ) ≡ F − x (cid:16) e i k · X F k ( ψ VRP ( x , y , z , R , d )) (cid:17) , (129)where X = ( X , Y , Z ) ∈ [0, 2 π ] are uniform random numbers. After that, werandomly rotate the ansatz by: ψ RVR ( x , y , z , R , d ) → ψ RVR ( x , y , z , R , d ) ψ RVR ( x , z , y , R , d ) ψ RVR ( y , x , z , R , d ) ψ RVR ( y , z , x , R , d ) ψ RVR ( z , x , y , R , d ) ψ RVR ( z , y , x , R , d ) . (130)Finally, the initial state ψ RPR is obtained by preparing N V different ansatze ψ RVR and multiplying them. Changing the radius of the ring R , the inter-vortex distance d or the number of vortex rings pairs N V will impact thebehaviour of QT. An example of the resulting flow is shown in Fig. 3.Figure 3: Illustration of the initial field preparation using random vortex ringpairs. Vortex lines (iso-surfaces of low ρ ) for the wave function ψ RVR (Eq.130) with, from left to right, N V =1, 20 and 50 vortex ring pairs.35 Numerical results
We describe in this section quantum turbulence flows simulated using thefour initial conditions described in the previous section: Taylor-Green (TG),Arnold-Beltrami-Childress (ABC), smoothed random phase (SRP) and ran-dom vortex rings (RVR). We present values, spectra and structure functionsof main quantities of interest (energy, helicity, etc.) that could be useful tobenchmark numerical codes simulating QT.The main physical and numerical parameters of the runs were fixedfollowing the scaling analysis provided in §4.1 and are summarized in Table1. Runs are identified using the abbreviation of the initial condition, followedby the identifier of the space resolution, e. g.
TG a is the run using theTaylor-Green initial condition and a 128 grid. Resolutions up to 1024 gridpoints (runs ” d”) were considered for some cases. For all simulations, thegrid is equidistant in each space direction, covering a domain of the same size[0, L ] , with L = 2 π . We recall that a Fourier spectral spatial discretizationwith periodic boundary conditions is used in the GP solvers. Run L N x M ref c α β k max ξ δx/ξ ξ N dv a 2 π
128 0.5 2.0 0.05000 40 2.26 1.388 0.035355 88b 2 π
256 0.5 2.0 0.02500 80 2.26 1.388 0.017678 177c 2 π
512 0.5 2.0 0.01250 160 2.26 1.388 0.008839 355d 2 π Table 1: Numerical and physical parameters used in the QT simulations.Using the contribution by Nore et al. (1997a) as guideline, the referenceMach number was fixed to M ref = 0.5, equivalent to a non-dimensional speedof sound c = 2. Consequently, αβ = 2 for all cases. Following (84), whenincreasing the grid resolution N x by a factor of 2, the value of the parameter α is diminished by the same factor in order to keep constant the value of k max ξ = 8 √ / ξ = √ αM ref diminisheswhen N x is increased, while the grid resolution of a vortex is kept constant δx/ξ = π/ ( k max ξ ) = 1.388. We check for the TG case that this grid resolutionis enough to accurately resolve vortices in our QT simulations. Since the size L of the computational box is kept constant, the higher the grid resolution N x , the larger is the number of vortices present in the domain (see values of N dv in Table 1). 36 .1 Benchmark The Taylor-Green initial field was prepared as described in §5.1. We displayin Table 2 the values of the time step δt used in the GP solver (see §4.3) andthe final time T f of each simulation. The parameters of the correspondingimaginary-time (IT) run cases preparing the initial condition using the ARGLEsolver are also presented, with δτ and τ f the imaginary-time step and finalvalue at convergence, respectively, and [ γ d /
4] the winding number of initialTG vortices seed at τ = 0 (see Eq. (110) and Figure 1). Run N x δt T f TG a 128 1.250e-3 12TG b 256 6.250e-4 12TG c 512 3.125e-4 12TG d 1024 3.125e-4 10 Run N x [ γ d / δτ τ f TG aIT 128 3 1.2500e-2 60TG bIT 256 6 6.2500e-3 60TG cIT 512 12 3.1250e-3 60TG dIT 1024 25 1.5625e-3 60
Table 2: Runs for the TG-QT case. Parameters used in the GP solver (casesTG a to TG d) and the imaginary-time (IT) ARGLE solver (cases TG aITto TG dIT). For each space resolution N x , the corresponding physical andnumerical parameters are displayed in Table 1. In defining this benchmark, it is important to describe in detail the initialfield obtained after the imaginary time (ARGLE) procedure. We recall thatthis procedure starts from the ansatz ψ ARGLE (110) containing multiple zeroTG vortices that split into [ γ d /
4] = [1 / (2 πα )] singly quantized vortices duringthe imaginary time propagation (see Fig. 1 illustrating the case TG aIT).Note from Table 2 that when increasing the grid resolution N x , the ansatzTG vortices split in a larger number of individual quantized vortices (up to25 for N x = 1024). This is illustrated in Fig. 4 showing vortex configurationsobtained at the end of the ARGLE procedure for runs TG aIT, TG bIT andTG cIT.To validate the ARGLE runs, we report in Table 3 the values of differentenergies (see §3.4) computed for the final field (at τ f ). The results are in goodagreement with those reported by Nore et al. (1997a,b). The values of thehelicity (see §3.6) are also reported in Table 3. Note that in this particular37igure 4: TG-QT. Initial condition computed with the imaginary-time AR-GLE solver. Vortex lines (iso-surfaces of low ρ ) of the converged wave function φ TG at final imaginary-time τ f . From left to right: grid resolutions N x = 128,256, 512 (corresponding to runs TG aIT, TG bIT and TG cIT in Table 2).case, we expect the helicity to be zero, which is satisfied for regularized helicity.As already stated by Clark di Leoni et al. (2017), the regularized helicity issmoother and less noisy, which explains the discrepancies for helicity in runsTG cIT and TG dIT. Run E ikin E ckin E q E int H H reg
TG aIT 0.12901707 4.8667051e-04 7.9239425e-03 1.2995235e-02 1.37e-13 -1.87e-11TG bIT 0.11334487 2.2334712e-04 4.0373670e-03 6.8223665e-03 2.96e-07 -5.52e-07TG cIT 0.12884207 1.5065059e-04 2.4895687e-03 4.2864757e-03 3.83e-03 -5.63e-07TG dIT 0.12968555 9.5590716e-05 1.3476259e-03 2.3466209e-03 -3.94e-04 -9.71e-08
Table 3: TG-QT. Values of different energies and helicity at τ f for theruns preparing the Taylor-Green initial condition, using the imaginary-timeARGLE solver. Starting from the initial condition presented in Fig. 4, we used the Strang–splitting GP solver (see §4.3) to advance the wave function in real time. Thefinal (at t = T f ) QT field is displayed in Fig. 5 for runs TG a, TG b andTG c. As explained before, when the grid resolution N x is increased, the sizeof a vortex core ξ diminishes and, consequently, the density of the tangled38ortex lines is increased in the computational box. Meanwhile, we recall thatthe grid resolution of the vortex core ( δx/ξ ) is the same for all simulations.Figure 5: TG-QT. Instantaneous fields computed with the real-time GPsolver, starting from the initial condition presented in Fig. 4. Vortex lines(iso-surfaces of low ρ ) of the wave function at final time T f . From left to right:grid resolutions N x = 128, 256, 512 (corresponding to runs TG a, TG b andTG c in Table 2).To compare our results with those reported by Nore et al. (1997a,b), theTG-QT fields are analyzed by providing in Fig. 6 the time evolution of theincompressible ( E ikin ) and compressible ( E ckin ) parts of the kinetic energy (91)for cases TG a to TG d. For each case, the incompressible kinetic energy isdominant at the beginning of the simulation, and slowly decreases in time,while the compressible part increases. We report in Fig. 7 (a) the spectrumof E ikin for the case TG c at different time instants of the computation. Forsmall k , the spectrum follows a (Kolmogorov-like) power law E ikin ( k ) ∼ k − / (dashed line in Fig. 7 a), especially for early times of the simulation. Theseresults concerning the incompressible energy evolution and its spectrum arein very good agreement with the numerical results reported by Nore et al.(1997a,b) for the grid resolution N x = 512. As a novel diagnostic tool of theturbulent field (not presented in Nore et al. (1997a,b)), we show in Fig. 7 (b)the time-evolution of the second-order structure function S // ( r ) (see Eq. (59)).For a developed QT field at t = 12, the slope of the structure function curveat the origin is close to 2, while for length scales larger than 0.1, the slopeevolves to 1/3. Using Eq. (60) to check the structure function calculation,we also plot in Fig. 7 (b) as a dotted line the value 2 R v x which is reachedfor large length scales (see Eq. (60)).39) t E i k i n TG_aTG_bTG_cTG_d b) t E c k i n TG_aTG_bTG_cTG_d
Figure 6: TG-QT. Time evolution of incompressible kinetic energy E ikin (a)and compressible kinetic energy E ckin (b) for runs TG a to TG d (see Table2).a) k k E ikin t = 0 E ikin t = 4 E ikin t = 6 E ikin t = 12 b) r S // t=0t=4t=6t=12 Figure 7: TG QT. Spectrum of E ikin (a) and second–order structure function(b) for the case TG c (see Table 2). There are two important check–points in validating a QT-GP simulation: theaccuracy in verifying conservation laws (see §2.1) and the grid convergence.In our numerical simulations, we monitor the time variation of the number ofparticles N (see Eq. (85)) and total energy per volume unit E (see Eq. (86)).These two quantities should be conserved by the GP solver. For the TG-QTsimulations, we report in Table 4 initial and final values for the norm and40ormalized energy, as well as their relative maximum variation during thetime evolution. Note from Table 4 that N is perfectly conserved, and energyrelative fluctuations δ ( E ) are less than 0.02%, which is sufficiently small valueto guarantee the validity of the computation. Run N | t =0 N | t = T f δ ( N ) E | t =0 E | t = T f δ ( E )TG a 0.9789997 0.9789997 0.0 0.1504230 0.1504111 7.97e-05TG b 0.9873160 0.9873160 0.0 0.1244280 0.1244235 3.60e-05TG c 0.9920083 0.9920083 0.0 0.1357688 0.1357597 6.67e-05TG d 0.9955732 0.9955732 0.0 0.1275738 0.1275550 1.47e-04 Table 4: TG-QT. Conservation of the number of particles N and energy pervolume unit. Initial (at t = 0) and final values (and t = T f ) and relative maxi-mum variation, defined following e. g. δ ( E ) = max t ∈ [0; T f ] | E ( t ) − E t =0 | /E t =0 .The second important check-point is the grid convergence. To correctlycapture vortices of radius ξ , we need enough discretization points in eachvortex core. We recall that the grid step size was fixed to δx/ξ = 1.338for all runs, corresponding to ξk max = √ ’ δx/ξ . Run N x α β ξk max δx/ξ TG a 128 0.05 40 2.26 1.388TG g 64 0.05 40 1.13 2.776TG h 256 0.05 40 4.52 0.694
Table 5: TG-QT. Supplementary runs used to check the grid convergence ofthe results (to be compared to runs in Table 1).In Table 6, we report the values of different energies obtained at the endof the imaginary-time ARGLE procedure for these new cases. Relative errorswere computed with respect to reference values of the case TG a. We concludethat a value of ξk max ’ i. e. δx/ξ = π/
2, is sufficient to ensure the gridconvergence and good accuracy of numerical results. To further check thisassessment, we also simulated the QT evolution starting from these runs. Wereport in Fig. 8 the time evolution of incompressible and compressible kineticenergies. The similarities between cases TG a and TG h suggest that the41esolution used for the case TG a is fine enough to capture the vortices inthe QT field. This validates the choice of parameters in Table 1.
TG a TG g (rel. err.) TG h (rel. err.) E ikin E ckin E q E int E v Table 6: TG-QT. Energies computed from φ T G , the wave function obtainedat the end of the imaginary-time ARGLE procedure for cases TG a, TG gand TG h. Relative errors (rel. err.) were computed with respect to referencevalues of the case TG a.a) t E i k i n TG_aTG_gTG_h b) t E c k i n TG_aTG_gTG_h
Figure 8: TG-QT. Time evolution of incompressible E ikin (a) and compressible E ckin (b) energies for cases TG a, TG g, TG h, used to check grid convergence.42 .2 Benchmark The ABC initial field was prepared as described in §5.2. We display in Table7 the values of the time step δt used in the GP solver (see §4.3) and the finaltime T f of each simulation. The parameters of the corresponding imaginary-time (IT) run cases preparing the initial condition using the ARGLE solverare also presented, with δτ and τ f the imaginary-time step and final value atconvergence, respectively. Run N x δt T f ABC a 128 8.0e-4 10ABC b 256 4.0e-4 10ABC c 512 2.0e-4 10 Run N x δτ τ f ABC aIT 128 4.0e-3 30ABC bIT 256 2.0e-3 30ABC cIT 512 1.0e-3 30
Table 7: Runs for the ABC-QT case. Parameters used in the GP solver(cases ABC a to ABC c) and the imaginary-time (IT) ARGLE solver (casesABC aIT to ABC cIT). For each space resolution N x , the correspondingphysical and numerical parameters are displayed in Table 1. Following (113) and (114), the initial condition for the imaginary–time AR-GLE procedure is obtained only by phase manipulations of the wave function.Therefore, vortices are not present at τ = 0, but they nucleate during theimaginary-time evolution, which a dissipative process. The obtained fieldswith vortices at the end of the ARGLE procedure are illustrated in Fig. 9.Note that, compared to the TG fields in Fig. 4, the distribution of vortices inthe computational box displays no symmetries with respect to central planes.This is the first feature that makes the ABC case different from the TG case.The second differentiating feature is the presence of helicity in the ABCflow obtained after the ARGLE procedure. We recall that the helicity of theTG flow is strictly zero (see Table 3). We report in Table 8 the values ofdifferent energies (see §3.4) and helicity computed for the final ABC field (at τ f ). As expected, the value of the incompressible kinetic energy is close to 1,which corresponds to the energy of the classical ABC flow. For the helicity,the theoretical value for the classical ABC flow is 3. A close value to 3 is43igure 9: ABC-QT. Initial condition computed with the imaginary-timeARGLE solver. Vortex lines (iso-surfaces of low ρ ) of the converged wavefunction φ ABC at final imaginary-time τ f . From left to right: grid resolutions N x = 128, 256, 512 (corresponding to runs ABC aIT, ABC bIT and ABC cITin Table 7).obtained only for large grid resolutions ( N x ≥ i. e. for sufficiently smallvalues of the vortex size ξ . Run E ikin E ckin E q E int H H reg
ABC aIT 0.9485114 0.0014218 0.0277287 0.0430379 2.4106982 2.4726091ABC bIT 0.9792042 0.0008142 0.0144931 0.0237201 2.7439625 2.6532860ABC cIT 0.9884992 0.0006486 0.0073802 0.0124975 2.7217161 2.7365301
Table 8: ABC-QT. Values of different energies and helicity at τ f for the runspreparing the ABC initial condition, using the imaginary-time ARGLE solver.For these computations, the ARGLE procedure required a significantcomputational time and was therefore stopped before the convergence criterion(96) was satisfied. To ensure the validity of the ARGLE solution, we estimatedthe criterion (97) by monitoring in Fig. 10a the imaginary-time evolutionof energy fluctuations defined as | E v ( φ n +1 ) − E v ( φ n ) | / ( δτ E v ( φ n )), with E v expressed by (92). The convergence criterion (97) is satisfied to a fairlygood degree of precision (10 − ). Figure 10b shows the spectra of E ikin , theincompressible kinetic energy of ARGLE solutions. This is an importantbenchmark verification, since E ikin represents the most important part in thetotal energy of the ABC super-flow (see Table 8). The similar slopes of thespectra for large wave numbers k indicate that the energy distribution of the44hree ABC flows are similar at small scales. The low- k part of the spectrum( k (cid:28) k ξ ) reproduces the classical spectrum which has only 2 nonzero modes, k = 1, 2. The slope for k (cid:29) k ξ is − i. e. E ikin ( k ) ∼ k − ). This feature ofthe high- k spectrum is detailed in Krstulovic and Brachet (2010).a) E n e r g y f l u c t u a t i o n ABC_aITABC_bITABC_cIT b) k E i k i n ABC_aITABC_bITABC_cIT
Figure 10: ABC-QT. (a) Relative fluctuation of the total energy( | E v ( φ n +1 ) − E v ( φ n ) | / ( δτ E v ( φ n ))) during the ARGLE computation. (b)Spectrum of E ikin , the incompressible kinetic energy of ARGLE solutions.Results for cases ABC aIT, ABC bIT and ABC cIT described in Table 7. Starting from the initial condition presented in Fig. 9, we used the Strang–splitting GP solver (see §4.3) to advance the wave function in real time. Thefinal (at t = T f ) QT field is displayed in Fig. 11 for runs ABC a, ABC b andABC c. As for the TG case, when the grid resolution N x is increased, the sizeof a vortex core ξ diminishes and, consequently, the density of the tangledvortex lines is increased.The time evolution of the incompressible kinetic energy E ikin and theregularized helicity H reg (see Eq. (72)) are shown in Fig. 12. These resultsare in good agreement with those reported by Clark di Leoni et al. (2017).To analyze the turbulent super-flow, we plot in Fig. 13 spectra for theincompressible kinetic energy E ikin (panel a) and the regularized helicity H reg (panel b) at different time instants and the second-order structure function S // ( r ) (panel c). We plot with dashed lines in Figs. 13a and 13b the reference45igure 11: ABC-QT. Instantaneous fields computed with the real-time GPsolver, starting from the initial condition presented in Fig. 9. Vortex lines(iso-surfaces of low ρ ) of the wave function at final time T f . From left to right:grid resolutions N x = 128, 256, 512 (corresponding to runs ABC a, ABC band ABC c in Table 7).a) t E i k i n ABC_aABC_bABC_c b) t H r e g ABC_aABC_bABC_c
Figure 12: ABC-QT. Time evolution of incompressible kinetic energy E ikin (a) and regularized helicity H reg (b) (see Eq. (72)).(Kolmogorov-like) power laws ε / k − / for E ikin and ηε − / k − / for helicity,respectively. The constants ε and η were computed as: ε = − dE ikin dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t =10 , η = − dHdt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t =10 . (131)We note from Figs. 13a and 13b that for t >
5, both energy and helicity H reg spectra exhibit at large scales a power law variation with exponent − / k k E ikin t= 0.0 E ikin t= 5.0 E ikin t= 7.0 E ikin t= 10.0 b) k k H reg t=0.0 H reg t=5.0 H reg t=7.0 H reg t=10.0 c) r S // t=0.0t=5.0t=7.0t=10.0 Figure 13: ABC-QT. Analysis of the turbulent super-flow. Spectrum ofincompressible kinetic energy E ikin (a) and regularized helicity H reg (b) atdifferent time instants for case ABC c. Dashed lines represent reference powerlaws with slope − /
3. Panel (c) displays the second order structure functionfor the same case ABC c and same time instants.compatible with a dual energy and helicity cascade. Again, this result isin good agreement with the results of Clark di Leoni et al. (2017). Thenovel diagnostic tool introduced in the previous section for the TG flow isalso performed with the ABC flow by computing the second–order structurefunction S // ( r ) (see Eq. (59)). Figure 13c displays the structure function forthe same case ABC c and same time instants considered for plotting spectra.A similar evolution as noted for the TG case (see Fig. 7) is observed: theslope of the structure function curve at the origin is close to 2, and, for largelength scales, the asymptotic value 2 R v x (dotted line).47 .2.3 Accuracy of numerical results and influence of the Machnumber As for the TG case, we monitor the time variation of the number of particles N (see Eq. (85)) and total energy per volume unit E (see Eq. (86)). Theaccuracy to which these two quantities are conserved by the GP solver isreported in Table 9 initial and final values for the norm and normalized energy,as well as their relative maximum variation during the time evolution. Notefrom Table 4 that N is perfectly conserved, and energy relative fluctuations δ ( E ) are less than 0.01%, which is sufficiently small value to guarantee thevalidity of the computation.Run N | t =0 N | t = T f δ ( N ) E | t =0 E | t = T f δE ABC a 0.9410138 0.9410138 0.0 1.0206996 1.0206400 5.98e-05ABC b 0.9647786 0.9647786 0.0 1.0182316 1.0181689 6.18e-05ABC c 0.9796053 0.9796053 0.0 1.0090255 1.0089631 6.19e-05Table 9: ABC-QT. Conservation of the number of particles N and energy pervolume unit. Initial (at t = 0) and final values (and t = T f ) and relative maxi-mum variation, defined following e. g. δ ( E ) = max t ∈ [0; T f ] | E ( t ) − E t =0 | /E t =0 .Another interesting question that can be addressed using the ABC flow isthe influence of the Mach number on the QT dynamics. Since the velocity v is singular at the vortex center r = 0, we considered in defining the localMach number the quantity √ ρ v which is not singular ( v ∼ /r and √ ρ ∼ r ,see §3.3). We thus computed two representative values: a maximum Machnumber M max based on the maximum superfluid velocity, and a Mach number M rms based on averaged values: M max := k√ ρ v k L ∞ ( D ) c , M rms := k√ ρ v k L ( D ) c √L = √ E kin c . (132)Keeping c and ξ constant, one can change the Mach number in the ABC flowby tuning the values of the parameters A , B , C in (114). Using as reference thecase ABC c ( N x = 512) we performed two new runs for which the parametersare displayed in Table 10. The values of constants A , B , C were divided(ABC c1) or multiplied (ABC c2) by a factor of 2. As a result, comparedto case ABC c, the velocities are divided (resp. multiplied) by 2 for caseABC c1 (resp. ABC c2). The values for the Mach number reported in Table480 were computed at the end of the ARGLE procedure preparing the initialcondition. Figure 14 shows the time evolution for the two values of the Machnumber, M max and M rms computed by the (real-time) GP solver. The ratioof 2 is well conserved in time, though the values are varying significantly. Thisproves that tuning the values of constants A , B , C is a simple and practicalapproach in modifying the Mach number of the QT super-flow.name ( A , B , C ) N x M max M rms ABC c (0.9, 1, 1.1) / √ / (2 √
3) 512 0.836509 0.357021ABC c2 2(0.9, 1, 1.1) / √ A , B , C in definingthe ABC flow were modified (see Eq. (114)).a) t m a x ABC_cABC_c1ABC_c2 b) t r m s ABC_cABC_c1ABC_c2
Figure 14: ABC-QT. Time evolution of the Mach number M max (a) and M rms (a) for cases ABC c, ABC c1 and ABC c2 (see Table 10).We present in Fig. 15 the time evolution of incompressible kinetic energy E ikin and regularized helicity H reg for new cases with different Mach numbers.As expected from the analysis above, the energy and helicity associatedwith the classical flow v ABC are divided (resp. multiplied) by 4 for caseABC c1 (resp. ABC c2). We note that the time evolution of these mainquantities depends on the Mach number. To assess on the distribution of theincompressible kinetic energy among scales, we plotted in Fig. 16 spectra of49 ikin at significant time instants, t = 5 and final time t = T f = 10. The spectrafor the three cases are quite similar showing that the obtained dynamics ofthe QT is equivalent when varying the Mach number of the flow.a) t E i k i n / E i k i n | t = ABC_cABC_c1ABC_c2 b) t H r e g / H r e g | t = ABC_cABC_c1ABC_c2
Figure 15: ABC-QT. Influence of the Mach number. Time evolution of theincompressible kinetic energy E ikin (a) and regularized helicity H reg (b). Tobe compared with curves in Fig. 12.a) k E i k i n ABC_cABC_c1ABC_c2 b) k E i k i n ABC_cABC_c1ABC_c2
Figure 16: ABC-QT. Influence of the Mach number. Spectrum of the in-compressible kinetic energy E ikin for case ABC c, ABC c1, ABC c2, at timeinstants t = 5 (a) and t = T f = 10 (b).50 .3 Benchmark The SRP initial field was prepared as described in §5.3. The advantage ofthis new initial condition is that the time-imaginary ARGLE simulation is nolonger necessary in the preparation of the initial field. We display in Table 11the values of the time step δt used in the GP solver (see §4.3), the final time T f of each simulation, and the parameters K (maximum amplitude of thephase) and N r (number of random values) of the method generating the phasefield (see Fig. 2). We recall that the characteristic velocity of the generatedflow field results is v = 2 α ( KN r /π ), and the corresponding theoretical Machnumber M = √ αKN r /π √ β . Run N x δt T f K N r SRP a 128 1 / π / π / π Table 11: Runs for the SRP-QT case. For each space resolution N x , thecorresponding physical and numerical parameters are displayed in Table 1.Figure 17 illustrates the vortex structures in the QT super-flow generatedwith this method. Compared to TG and ABC cases, in the SRP case vorticesnucleate progressively and do not display long vortex lines. A very fine grainstructure of vortices is observed in all SRP runs.To analyze the SRP-QT flow we plotted in Fig. 18a the time evolution ofthe compressible E ckin and incompressible E ikin kinetic energies. An ensembleaverage for 10 different (random) initial conditions was taken to displaythe results. Since the initial filed (at t = 0) does not contain vortices, theincompressible kinetic energy E ckin is initially zero and subsequently increasesdue to vortex nucleations. After reaching the maximum value at t ∼ E ckin gradually decreases to the end of the simulation ( t = T f ). During the entiretime evolution, the dynamics of the flow is dominated by the compressiblekinetic energy E ckin , which is always larger than E ikin . Figure 18b shows thespectrum of E ikin . As for TG and ABC cases, a Kolmogorov-like scaling isobtained, with a − / k . Hereof, the SRP-QTflow is statistically similar to the TG and ABC QT flows and can be usedin a detailed parametric study of the decay of quantum turbulence (which isbeyond the scope of this contribution).51igure 17: SRP-QT. Instantaneous fields computed with the real-time GPsolver, starting from the initial condition presented in Fig. 2. Vortex structures(iso-surfaces of low ρ ) of the wave function at final time T f . Grid resolution N x = 512, corresponding to run SRP c in Table 11).a) t E ckin SRP_a E ckin SRP_b E ckin SRP_c E ikin SRP_a E ikin SRP_b E ikin SRP_c b) k E i k i n t = 0.5 t = 1 t = 1.5 t = 2 k Figure 18: SRP-QT. (a) Time evolution of the compressible E ckin and incom-pressible E ikin kinetic energies. (b) Spectrum of E ikin at different time instants.Case SRP c ( N x = 512). In both panels, the results represent an ensembleaverage for 10 different (random) initial conditions.52 .4 Benchmark The RVR initial field was prepared as described in §5.4. Like in the SRPcase, building this new initial condition avoids the use of the time-imaginaryARGLE computation. We display in Table 12 the values of the time step δt used in the GP solver (see §4.3), the final time T f of each simulation, theparameter N V representing the number of pairs of vortex rings seeded in theinitial field, the radius R of a vortex ring, and the distance d between thevortex rings forming a pair (see Eq. (130)). Run N x δt T f N V R d
RVR a 128 1 / π/ π RVR b 256 1 / π/ π RVR c 512 1 / π/ π Table 12: Runs for the RVR-QT case. For each space resolution N x , thecorresponding physical and numerical parameters are displayed in Table 1.Note that in Fig. 3 we represented, to illustrate the method, a few numberof vortex pairs ( N V =1, 20 and 50). In the GP calculations we used a muchlarger value for N V , up to 800 for the case RVR c. The initial field for thethree considered cases is displayed in Fig. 19. Like in the TG and ABC cases,when the grid resolution N x is increased, ξ diminishes and, consequently,thinner vortex rings are seeded in the initial field.The obtained RVR-QT flow is illustrated in Fig. 20. Multiple vortex ringreconnections lead to a dense vortex distribution in the QT field, similar tothat obtained for the ABC flow (see Fig. 11).For the analysis of the RVR-QT flow we provide in Fig. 21(a) the timeevolution of the compressible E ckin and incompressible E ikin kinetic energies forthe case RVR c. Since the initial distribution of vortex rings pairs is randomin the computational box, we present the ensemble average results for 10 runswith random positions of the same number of vortex ring pairs ( N V = 800).In the early stages of the time evolution ( t < E ikin is dominant. Thecompressible kinetic energy E ckin starts to increase at t ∼
1, due to soundemissions through vortex reconnections. This evolution is opposite to thatobserved for the SRP-QT cases. Figure 21b shows the spectrum of E ikin . Likein the SRP cases (see Fig. 18), we note a Kolmogorov-like scaling of thespectrum, with a − / k .53igure 19: RVR-QT. Initial field containing N V randomly distributed vortexring pairs. Vortex lines (iso-surfaces of low ρ ) of the wave function. Fromleft to right: grid resolutions N x = 128, 256, 512 and N V = 200, 400, 800(corresponding to runs RVR a, RVR b and RVR c in Table 12).Figure 20: RVR-QT. Instantaneous fields computed with the real-time GPsolver, starting from the initial condition presented in Fig. 19. Vortex lines(iso-surfaces of low ρ ) of the wave function at final time T f . From left to right:grid resolutions N x = 128, 256, 512 (corresponding to runs RVR a, RVR band RVR c in Table 12). We simulated in this paper quantum turbulence superfluid flows describedby the Gross-Pitaevskii equation. Numerical simulations were performedusing a parallel (MPI-OpenMP) code based on a pseudo-spectral spatialdiscretization and second order splitting for the time integration. As expected54) t E ckin RVR_a E ckin RVR_b E ckin RVR_c E ikin RVR_a E ikin RVR_b E ikin RVR_c b) k E i k i n t = 0 t = 1 t = 2 t = 3 k Figure 21: RVR-QT. (a) Time evolution of the compressible E ckin and incom-pressible E ikin kinetic energies. (b) Spectrum of E ikin at different time instants.Case RVR c ( N x = 512). In both panels, the results represent an ensembleaverage for 10 different runs, with random initial distribution of N V = 800vortex ring pairs in the computational domain.from the theoretical numerical analysis, this approach ensured an accuratecapture of the dynamics of the flow, with a perfect conservation of thenumber of particles and a negligible drift in time of the total energy. Severalconfigurations of QT were simulated using four different initial conditions:Taylor-Green (TG) vortices, Arnold-Beltrami-Childress (ABC) flow, smoothedrandom phase (SRP) fields and random vortex rings (RVR) pairs. Each ofthese case was described in detail by setting corresponding benchmarksthat could be used to validate/calibrate new GP codes. Particular carewas devoted in describing dimensionless equations, characteristic scales andoptimal numerical parameters. We presented values, spectra and structurefunctions of main quantities of interest (energy, helicity, etc.) that are usefulto describe the turbulent flow. Some general features of QT were identified,despite the variety of initial states: the spectrum of the incompressible kineticenergy exhibits a Kolmogorov-type − / Acknowledgements
The authors acknowledge financial support from the French ANR grantANR-18-CE46-0013 QUTE-HPC. Part of this work used computational re-sources provided by IDRIS (Institut du d´eveloppement et des ressources eninformatique scientifique) and CRIANN (Centre R´egional Informatique etd’Applications Num´eriques de Normandie).
A Parallel performance of the code
A.1 Execution time
Execution times for runs ABC a to ABC c and ABC aIT to ABC cIT arereported in Table 13. 56ase N x Iterations MPI proc. Execution time (s) ratioABC aIT 128 7500 56 539.126 0.000257075ABC bIT 256 15000 112 5940.076 0.000354056ABC cIT 512 30000 224 53066.445 0.000395375ABC a 128 12500 56 700.524 0.000334036ABC b 256 25000 112 7364.754 0.000438973ABC c 512 50000 224 63397.341 0.000472346Table 13: Execution time for ABC runs. The last column reports the executiontime divided by the number of degrees of freedom ( N x ).When switching from one case to the next one, we doubled the totalnumber of iterations and also the number of processes. We expected a smallvariation of the value of the execution time divided by the grid resolution.For the ARGLE procedure, we monitored an efficiency of 65% from caseABC aIT to case ABC cIT. For the time-dependent GP simulation, weobtained an efficiency of 70% from case ABC a to case ABC c. Note that themeasured time is the total time for the execution of the program, not solelythe computational part of the code. A.2 Strong scalability of GPS
Strong scalability results of the GPS code are presented in figure 22. A3D test case (with grid resolutions up to 2048 ) was performed using adifferent number of processes (up to 64536), and the execution time wasmonitored. Strong scalability tests using only MPI (Fig. 22a) or the hybridcode MPI-OpenMP (Fig. 22b) show scalability and speed-up close to idealperformances. BibliographyReferences
Abid, M., Huepe, C., Metens, S., Nore, C., Pham, C., Tuckerman, L., Brachet,M., 2003. Gross-Pitaevskii dynamics of Bose-Einstein condensates andsuperfluid turbulence. Fluid Dynamics Research 33 (5-6), 509–544.57 (~ 1/np) (~ np) ( ~ np) Figure 22: Parallel performance of the GPS code when computing 3D caseswith 512 , 1024 and 20483