Quantum turbulence simulations using the Gross–Pitaevskii equation: High-performance computing and new numerical benchmarks

Abstract

We present high-performance and high-accuracy numerical simulations of quantum turbulence modelled by the Gross–Pitaevskii equation for the time-evolution of the macroscopic wave function of the system. The hydrodynamic analogue of this model is a flow in which the viscosity is absent and all rotational flow is carried by quantized vortices with identical topological line-structure and circulation. Numerical simulations start from an initial state containing a large number of quantized vortices and follow the chaotic vortex interactions leading to a vortex-tangle turbulent state. The Gross–Pitaevskii equation is solved using a parallel (MPI-OpenMP) code based on a pseudo-spectral spatial discretization and second order splitting for the time integration. We define four quantum-turbulence simulation cases based on different methods used to generate initial states: the first two are based on the hydrodynamic analogy with classical Taylor–Green and Arnold–Beltrami–Childress vortex flows, while the other two methods use a direct manipulation of the wave function by generating a smoothed random phase field, or seeding random vortex-ring pairs. The dynamics of the turbulent field corresponding to each case is analysed in detail by presenting statistical properties (spectra and structure functions) of main quantities of interest (energy, helicity, etc.). Some general features of quantum turbulence are identified, despite the variety of initial states. Numerical and physical parameters of each case are presented in detail by defining corresponding benchmarks that could be used to validate or calibrate new Gross–Pitaevskii codes. The efficiency of the parallel computation for a reference case is also reported.

Introduction

The study of quantum fluids, realized in superfluid helium and atomic Bose–Einstein condensates (BEC), has become a central topic in various fields of physics, such as low temperature physics, fluid dynamics of inviscid flows, quantum physics, statistical physics, cosmology, etc. One of the striking features of quantum fluids is the nucleation of vortices with quantized (fixed) circulation, when an external forcing is applied (rotation, stirring, etc.). The observation of quantized vortices, as a signature of the superfluid (zero-viscosity) nature of these flow systems, was extensively explored in different experimental settings of superfluid helium or BEC. Configurations with a large number of quantized vortices tangled in space can evolve to Quantum Turbulence (QT), generally referred to as vortex tangle turbulence. While QT in superfluid helium has been largely studied in the last two decades (see dedicated volumes [1], [2], [3], [4]), only recent experimental and theoretical studies [5], [6], [7], [8] reported different possible routes to QT in BEC.

A promising path of research in exploring QT is based on the analogy with classical turbulence (CT), observed in conventional viscous fluids and governed by the Navier–Stokes equations. Classical turbulent flows are characterized by the chaotic motion of vortical eddies that populate a continuous hierarchy of intensities and scales, from the large (integral) scale of the flow, down to the Kolmogorov’s viscous length scale. The classical turbulent cascade of energy between scales is characterized by the Kolmogorov’s power-law spectrum [9] in the regime of vanishing viscosity (i. e. large to infinite Reynolds numbers). Quantum turbulence appears then as an equivalent regime, since superfluids are assimilated to flows with zero viscosity. Both finite- and zero-viscosity regimes can be experimentally obtained in liquid helium by changing the temperature: above the lambda transition temperature (2.17 K) the liquid is normal (viscous) and well below it is a pure superfluid. Experimental measurements in superfluid helium-4 at temperatures below 2 K [10] provided indeed evidence of Kolmogorov’s law for the kinetic energy cascade. However, vortex interaction mechanisms are different in the two types of turbulence. Unlike classical vortex eddies, vortices in QT are identical topological line defects in the fluid density field and their circulation is quantized (in units of Planck’s constant over the atomic mass). In QT, bundles of quantum line vortices play the role of classical vortex eddies. Since visualizations in QT experiments are not yet enough precise to provide an accurate image of vortex interactions, numerical simulations are then needed. The vortex filament (VF) and the Gross–Pitaevskii (GP) models are used in the literature to numerically explore vortex interactions mechanisms in QT. The VF model represents quantized vortices as infinitely thin lines and follows their evolution by integrating the Biot–Savart–Laplace law over the vortex filament tangle. This model proved very useful in studying superfluid helium-4 [11]. The GP model is the simplest mathematical model for a superfluid at zero-temperature and it will be the focus of this paper. It can be also regarded as a theoretical and numerical framework used to investigate the inviscid limit of a fully-developed CT. For a comprehensive description of different models of QT, see recent reviews by Halperin and Tsubota [4], Brachet [12], Barenghi et al. [13] and Tsubota et al. [11].

The Gross–Pitaevskii equation (GPE) is a nonlinear Schrödinger equation with cubic nonlinearity. It describes, in the theoretical limit of absolute zero temperature, the time evolution of the macroscopic complex wave function $\psi$ (identical for all particles) defining a weakly interacting Bose system. Consequently, the GP model naturally applies to numerical studies of QT in Bose–Einstein condensates [14]. The relevance of the GPE for describing QT in superfluid helium is discussed in Section 3.3. For general QT flows, the GP model offers the advantage that quantized vortices appear naturally as topological line defects, resulting from the $U\left(1\right)$ symmetry breaking of the phase shift of $\psi$. Subsequent vortex phenomena (reconnection, annihilation) are intrinsically described by the model. Supplementary phenomenological models are needed in the VF approach to take into account the same vortex mechanisms. As a consequence, QT simulations based on the GPE (denoted hereafter as GPE-QT) were largely used in the literature to study vortex interactions and statistical properties of QT in the zero-temperature limit.

There are several challenges when setting a numerical simulation to investigate GPE-QT: (i) generate a physically and mathematically sound initial state with many quantized vortices that finally evolve to a statistically steady state of QT, (ii) use accurate numerical methods that preserve the invariants of the GPE when long time-integration is necessary, (iii) design numerical codes affording large grid resolutions, necessary to accurately capture the dynamics of vortices and (iv) compute appropriate (statistical) diagnostic tools to analyse the superfluid flow evolution.

We use in this contribution a modern parallel (MPI-OpenMP) numerical code satisfying the requirements (ii) and (iii). The code is called GPS (Gross–Pitaevskii Simulator) [15] and is based on a Fourier-spectral space discretization and up-to-date numerical methods: a semi-implicit backward-Euler scheme with Krylov preconditioning for the stationary GP equation [16] and various schemes (Strang splitting, relaxation, Crank–Nicolson) for the real-time GP equation [17]. The GPS code offers a solid framework to address in detail challenges (i) and (iv), for which we review previous models and bring new contributions.

A great deal of attention has been lately devoted to the development of accurate numerical schemes to solve different forms of the GPE, from the classical (stationary or time-dependent) GPE, to systems of coupled GPEs and more recent formulations (e. g. with non-local or high-order interactions). For recent reviews of numerical methods for GPE, see  [17], [18], [19], [20], [21]. Several software packages for solving the GPE were deposited in the CPC Program Library. The spatial discretization is generally based on spectral [16], [22], [23], [24], finite-elements [25], [26] or finite-difference [27], [28], [29], [30], [31], [32] methods. Provided programs are written in Fortran [22], [27], C [28], [29], Matlab [16], [23], [24], [29], [31], FreeFem++ [26] or C and Fortran with OpenMP [32]. None of these programs deal with QT simulations, but can eventually be extended to perform such simulations using the procedures explained in this paper.

As in classical turbulence, the numerical and physical accuracy of the initial condition is crucial in computing properties of numerically generated QT. Using the hydrodynamic analogy for the GP model (through the Madelung transform, as explained below), pioneering numerical simulations of GP-QT [33], [34], [35] suggested initial conditions and statistical analysis tools inspired from CT. A velocity field, derived from the well-known classical flow with Taylor–Green (TG) vortices, was imposed to the superfluid flow. An initial wave function, with nodal lines corresponding to vortex lines of the velocity, was thus generated. This initial wave function was 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 result of the ARGLE procedure was finally used as initial field for the time-dependent GP simulation. A similar approach was more recently used by replacing the TG vortices with the Arnold–Beltrami–Childress (ABC) classical vortex flow [36], [37]. If this approach is well suited to control the hydrodynamic characteristics of the initial superfluid flow (Mach number, helicity), it involves supplemental technicalities and computations through the ARGLE procedure. We suggest in this paper two new approaches to generate the initial condition for the GP-QT simulations, based on the direct manipulation of the wave function. The ARGLE procedure is thus avoided. The first method prescribes a smoothed random-phase (SRP) for the wave function, while the second one generates random vortex rings (RVR). The two new methods, which are simple to implement, are shown to develop QT fields with similar statistical properties as those obtained using the TG or ABC classical initial conditions. Nevertheless, the dynamics of the superfluid flow is different. Compared to TG and ABC cases, in the SRP case the initial field is vortex free and dominated by the compressible kinetic energy; vortices nucleate progressively and do not display long vortex lines. For the RVR flow, evolution is opposite to that observed for the SRP case: in the early stages of the time evolution the incompressible kinetic energy is dominant; the compressible kinetic energy then starts to increase due to sound emissions through vortex reconnections. However, like in well-documented TG and ABC cases, a Kolmogorov-like scaling of the incompressible kinetic energy spectrum is obtained for the new SRP and RVR cases.

Concerning the analysis of the QT field, we present classical diagnostic tools (inspired from CT), as the energy decomposition and associated spectra [33], [34] and also new ones, as the second-order structure function, not reported in the previously cited studies. This supplements the statistical description of the superfluid flow. We also carefully investigate the influence of numerical parameters (as the grid resolution of a vortex, the maximum resolved wave-number and the computed local Mach number) on the characteristic of QT. This topic is generally very briefly addressed in physical 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 also intended to define in detail numerical benchmarks in the framework of parallel computing. We start by revisiting classical GP-QT settings (based on TG and ABC flows). New results obtained with our high-performance/high-accuracy parallel code are compared with available data in the literature. We then present the new numerical benchmarks, based on random phase fields or random vortex rings generation. The new benchmarks offer a new perspective in comparing the classical CT based approaches to more GP-oriented models. For all benchmarks, we offer a comprehensive description of the numerical and physical parameters and give checkpoint values for validating each step of the simulation. This could be useful to assess new numerical methodologies or tune/validate new modern GP numerical codes.

The organization of the paper is as follows. Section 2 introduces the GP mean field equation, its stationary version and the hydrodynamic analogy. Section 3 presents the mathematical and physical QT model used in the present numerical simulations. The main characteristics of the quantum flow (healing length, sound velocity) and the notion of quantized vortex are introduced. The Bogoliubov dispersion relation and the validity of the model are also discussed. In Section 4 we present the main diagnostic quantities generally used to analyse a QT flow: spectra of compressible and incompressible kinetic energies, velocity structure functions and helicity. The numerical methods used to solve the GP equation and the associated GPS code are described in Section 5. An extended subsection is devoted to the derivation of dimensionless equations. The particular numerical methods used in this study to advance the GP wave function in imaginary-time (ARGLE) or real-time (GP) are also described. Section 6 presents in detail four different approaches to generate the initial field for the simulation of decaying GP-QT: Taylor–Green (TG), Arnold–Beltrami–Childress (ABC), smoothed random phase (SRP) and random vortex rings (RVR). Each method is associated to a benchmark. The results obtained for the four benchmarks are discussed in Section 7. We present values, spectra and structure functions of main quantities of interest (energy, helicity, etc.) that could be useful to benchmark numerical codes simulating QT with the GP model. Finally, the main features of the benchmarks and their possible extensions are summarized in Section 8. Appendix presents strong scalability tests of the computation code used to run a specific case (the ABC flow) in both MPI and hybrid MPI-OpenMP configurations.