# Turbulence in the two-dimensional Fourier-truncated Gross-Pitaevskii equation

Vishwanath Shukla, Marc Brachet, Rahul Pandit

► To cite this version:

Vishwanath Shukla, Marc Brachet, Rahul Pandit. Turbulence in the two-dimensional Fourier-truncated Gross-Pitaevskii equation. New Journal of Physics, 2013, 15, pp.113025. 10.1088/1367-2630/15/11/113025. [hal-01587337](https://hal-fr.sorbonne-universite.fr/hal-01587337)

HAL Id: hal-01587337

https://hal.sorbonne-universite.fr/hal-01587337

Submitted on 14 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Turbulence in the two-dimensional Fourier-truncated Gross-Pitaevskii equation

To cite this article: Vishwanath Shukla et al 2013 *New J. Phys.* **15** 113025

View the [article online](#) for updates and enhancements.

## Related content

- Quench dynamics of the 2d XY model
  Asja Jeli and Leticia F Cugliandolo

- Systematics of the magnetic-Prandtl-number dependence of homogeneous, isotropic magnetohydrodynamic turbulence
  Ganapati Sahoo, Prasad Perlekar and
  Rahul Pandit

- Berezinskii-Kosterlitz-Thouless transition with a constraint lattice action
  Wolfgang Bietenholz, Urs Gerber and
  Fernando G Rejón-Barrera

## Recent citations

- Delayed collapses of Bose-Einstein condensates in relation to anti-de Sitter gravity
  Anxo F. Biasi et al

- The onset of thermalization in finite-dimensional equations of hydrodynamics: insights from the Burgers equation
  Divya Venkataraman and Samriddhi
  Sankar Ray

- Multiscaling in superfluid turbulence: A shell-model study
  Vishwanath Shukla and Rahul Pandit
# Turbulence in the two-dimensional Fourier-truncated Gross-Pitaevskii equation

Vishwanath Shukla<fnref n="1" />, Marc Brachet<fnref n="2" /> and Rahul Pandit<fnref n="1" />,<fnref n="3" />

E-mail: vishwanath@physics.iisc.ernet.in, brachet@physique.ens.fr and rahul@physics.iisc.ernet.in

New Journal of Physics 15 (2013) 113025 (32pp)

Received 4 February 2013

Published 11 November 2013

Online at http://www.njp.org/

doi:10.1088/1367-2630/15/11/113025

**Abstract.** We undertake a systematic, direct numerical simulation of the two-dimensional, Fourier-truncated, Gross-Pitaevskii equation to study the turbulent evolutions of its solutions for a variety of initial conditions and a wide range of parameters. We find that the time evolution of this system can be classified into four regimes with qualitatively different statistical properties. Firstly, there are transients that depend on the initial conditions. In the second regime, power-law scaling regions, in the energy and the occupation-number spectra, appear and start to develop; the exponents of these power laws and the extents of the scaling regions change with time and depend on the initial condition. In the third regime, the spectra drop rapidly for modes with wave numbers $k > k_c$ and partial thermalization takes place for modes with $k < k_c$; the self-truncation wave number $k_c(t)$ depends on the initial conditions and it grows either as a power of $t$ or as $\log t$. Finally, in the fourth regime, complete thermalization is achieved and, if we account for finite-size effects carefully, correlation functions and spectra are consistent with their nontrivial Berezinskii-Kosterlitz-Thouless

<footnote n="1">¹Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India</footnote>

<footnote n="2">²Laboratoire de Physique Statistique de l'Ecole Normale Supérieure, associé au CNRS et aux Universités Paris VI et VII, 24 Rue Lhomond, F-75231 Paris, France</footnote>

<footnote n="3">³ Also at Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore, India.</footnote>
forms. Our work is a natural generalization of recent studies of thermalization in the Euler and other hydrodynamical equations; it combines ideas from fluid dynamics and turbulence, on the one hand, and equilibrium and nonequilibrium statistical mechanics on the other.

S Online supplementary data available from stacks.iop.org/NJP/15/113025/mmedia

## Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Model, initial conditions and numerical methods</td><td>4</td></tr><tr><td>2.1. The Gross-Pitaevskii equation</td><td>4</td></tr><tr><td>2.2. Numerical methods and initial conditions</td><td>6</td></tr><tr><td>3. Results</td><td>8</td></tr><tr><td>3.1. Evolution of energies, velocity probability distribution functions and the zero-wave-number population</td><td>8</td></tr><tr><td>3.2. Initial transients and the onset of thermalization</td><td>13</td></tr><tr><td>3.3. Partial thermalization and self-truncation</td><td>15</td></tr><tr><td>3.4. Complete thermalization</td><td>20</td></tr><tr><td>4. Conclusions</td><td>23</td></tr><tr><td>Acknowledgments</td><td>24</td></tr><tr><td>Appendix A</td><td>24</td></tr><tr><td>Appendix B</td><td>25</td></tr><tr><td>References</td><td>29</td></tr></table>

## 1. Introduction

The elucidation of the nature of superfluid turbulence, which began with the pioneering studies of Feynman [1] and Vinen and Hall [2–6], has continued to engage the attention of experimentalists, theoreticians and numerical simulators [7–13] and has experienced a renaissance over the last few years. Experimental systems, in which such turbulence is studied, include the bosonic superfluid $^4$He, its fermionic counterpart $^3$He and Bose-Einstein condensates (BECs) of cold atoms in traps and their optical analogues; for representative studies, we refer the reader to [14–23]. Theoretical and numerical studies have used a variety of models to study superfluid turbulence; these include the two-fluid model [24, 25], Biot-Savart-type models with [26, 27] or without [28, 29] the local-induction approximation and the Gross-Pitaevskii (GP) or nonlinear Schrödinger equations [30, 31]. These models have been studied by a combination of theoretical methods, such as wave-turbulence theory [30–33] and numerical simulations [34–40]. Most of these studies have been carried out in three dimensions (3D); numerical simulations of two-dimensional (2D) models for superfluid turbulence have been increasing steadily over the last few years [41–44]. Here we undertake a systematic direct numerical simulation (DNS) of the dissipationless, unforced, Fourier-truncated, 2D GP equation with a view to identifying what, if any, features of the turbulent evolution of the solutions of this equation are universal, i.e. they do not depend on initial conditions. Some, although not
all, of our results are contained in earlier simulations [41–48]. The perspective of our study is different from earlier studies of the 2D GP equation; in particular, we elucidate in detail the dynamical evolution of this system and examine the various stages of its thermalization; in this sense our work is akin to recent studies of thermalization in Euler and other hydrodynamical equations [49–51], which combine ideas from fluid dynamics and turbulence, on the one hand, and equilibrium and nonequilibrium statistical mechanics on the other. Recent studies [49–51] of the dynamics of spectrally truncated, 3D, incompressible Euler flows and related systems have shown that the inviscid and conservative Euler equation, with a high-wave-number spectral truncation, has long-lasting transients that behave just as those of the 3D dissipative Navier–Stokes equation, with generalized dissipation. This is so because the thermalized modes, between some transition wave number and the maximum wave number, act as an effective microworld that provides viscosity to the modes, with wave numbers below the transition wave number; a similar study for the 3D GP equation has been carried out by Krstulovic and Brachet [38, 52].

In a recent review on quantum turbulence, Paoletti and Lathrop [12] write, ‘Despite the abundant examples of turbulence, there is no consensus definition of the term. Here, we define turbulence as a dynamic field that is spatially complex, aperiodic in time, and involves processes spanning several orders of magnitude in spatial extent and temporal frequency’. It is in this sense that we use the term turbulence in our study of the dynamical evolution of solutions of the 2D, Fourier-truncated GP equation.

It is useful to begin with a qualitative overview of our principal results. We find that the dynamical evolution of the dissipationless, unforced, 2D, Fourier-truncated GP equation can be classified, roughly, into the following four regimes, which have qualitatively different statistical properties. (i) The first is the region of initial transients; this depends on the initial conditions. (ii) This is followed by the second regime, in which we see the onset of thermalization; here the energy and occupation-number spectra begin to show power-law-scaling behaviour, but the power-law exponent and the extents of the scaling regions change with time and depend on the initial conditions. (iii) In the third regime, which we call the region of partial thermalization, these spectra show clear, power-law, scaling behaviour, with a power that is independent of the initial conditions, and, at large wave vectors, an initial-condition-dependent, self-truncation regime, where spectra drop rapidly. (iv) Finally, in the fourth regime, the system thermalizes completely and exhibits correlation functions that are consistent with the predictions of the Berezinskii–Kosterlitz–Thouless (BKT) theory [47, 53–55], if the simulation domain and simulation time are large enough. Although some of these regimes have been seen in some earlier numerical studies of the 2D GP equation, we are not aware of any study that has systematized the study of these four dynamical regimes. In particular, regime 3, which shows partial thermalization and self-truncation in spectra, has not been identified in the 2D, Fourier-truncated, GP equation, even though its analogue has been investigated in the 3D case [32, 38, 52].

The remaining part of this paper is organized as follows. In section 2, we describe the 2D GP equation and the different statistical measures we use to characterize turbulence in the Fourier-truncated, 2D GP equation (section 2.1); the details of our numerical methods and initial conditions are given in section 2.2. In section 3, we present our results; these are described in the four sections 3.1–3.4 that are devoted, respectively, to the following: (a) the temporal evolution of the energy components, velocity-component probability distribution functions (PDFs) and the population $N_0$ in the zero-wave-number mode; (b) the statistical characterization of the first
two regimes of the dynamical evolution (by using various energies and the occupation-number spectra for different initial conditions); (c) a similar statistical characterization, as in section 3.2, but for the regime with partial thermalization, and the study of the nature of the growth of the self-truncation region; (d) the final, completely thermalized state of the Fourier-truncated, 2D GP equation. Section 4 contains our conclusions. A note on the units used for the GP equation and details of some analytical calculations are presented in appendices A and B, respectively.

## 2. Model, initial conditions and numerical methods

In this section, we describe the 2D GP equation. We define all the statistical measures that we use to characterize the time evolution of this equation, given the three types of initial conditions that we describe below. We also describe the numerical methods and computational procedures that we use to solve this equation.

### 2.1. The Gross-Pitaevskii equation

The GP equation, which describes the dynamical evolution of the wave function $\psi$ of a weakly interacting 2D Bose gas at low temperatures, is

$$i \frac{\partial \psi(\mathbf{x}, t)}{\partial t} = -\nabla^2 \psi(\mathbf{x}, t) + g |\psi|^2 \psi(\mathbf{x}, t), \quad (1)$$

$\psi(\mathbf{x}, t)$ is a complex, classical field and $g$ is the effective interaction strength [56, 57]. This equation conserves the total energy

$$E = \int_{\mathcal{A}} \left[ |\nabla \psi|^2 + \frac{1}{2} g |\psi|^4 \right] d^2 x \quad (2)$$

and the total number of particles

$$N = \int_{\mathcal{A}} |\psi|^2 d^2 x, \quad (3)$$

where $\mathcal{A} = L^2$ is the area of our 2D, periodic, computational domain of side $L$. From (1) we obtain the continuity equation

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \quad (4)$$

where $\rho = |\psi|^2$ is interpreted as the particle density and the velocity is

$$\mathbf{v}(\mathbf{x}, t) = \frac{\psi^* \nabla \psi - \psi \nabla \psi^*}{i |\psi|^2}. \quad (5)$$

We can use the Madelung transformation $\psi(\mathbf{x}, t) = \sqrt{\rho} e^{i\theta(\mathbf{x}, t)}$, where $\theta(\mathbf{x}, t)$ is the phase of $\psi(\mathbf{x}, t)$, to write $\mathbf{v}(\mathbf{x}, t) = 2\nabla\theta(\mathbf{x}, t)$, whence we obtain [35]

$$E = \int_{\mathcal{A}} \left[ \frac{1}{4} \rho v^2 + \frac{1}{2} g |\psi|^4 + [\nabla \rho^{1/2}]^2 \right] d^2 x = E_{\text{kin}} + E_{\text{int}} + E_{\text{q}}, \quad (6)$$

where the kinetic, interaction and quantum-pressure energies are defined, respectively, as

$$E_{\text{kin}} = \frac{1}{4} \int_{\mathcal{A}} |\sqrt{\rho} v|^2 d^2 x, \quad (7a)$$
$$E_{\text{int}} = \frac{1}{2} \int_A g |\psi|^4 d^2x, \qquad (7b)$$

$$E_q = \int_A |\nabla \rho^{1/2}|^2 d^2x. \qquad (7c)$$

We separate the compressible (superscript c) and the incompressible (superscript i) parts of the kinetic energy by making use of the decomposition

$$\rho^{1/2} \mathbf{v} = (\rho^{1/2} \mathbf{v})^i + (\rho^{1/2} \mathbf{v})^c, \qquad (8)$$

where $\nabla \cdot (\rho^{1/2} \mathbf{v})^i = 0$ and $\nabla \times (\rho^{1/2} \mathbf{v})^c = 0$, whence we obtain the following:

$$E_{\text{kin}}^i = \frac{1}{4} \int_A |(\sqrt{\rho}v)^i|^2 d^2x, \qquad (9a)$$

$$E_{\text{kin}}^c = \frac{1}{4} \int_A |(\sqrt{\rho}v)^c|^2 d^2x. \qquad (9b)$$

The spectra for these energies are defined as follows:

$$E_{\text{kin}}^i = \frac{1}{4} \int |(\widehat{\rho^{1/2} \mathbf{v}})^i|^2 d^2k \equiv \int E_{\text{kin}}^i(k) dk, \qquad (10)$$

$$E_{\text{kin}}^c = \frac{1}{4} \int |(\widehat{\rho^{1/2} \mathbf{v}})^c|^2 d^2k \equiv \int E_{\text{kin}}^c(k) dk, \qquad (11)$$

$$E_{\text{int}} = \int |\sqrt{g/2}| |\psi|^2 d^2k \equiv \int E_{\text{int}}(k) dk \qquad (12)$$

and

$$E_q = \int |\nabla \widehat{\rho^{1/2}}|^2 d^2k \equiv \int E_q(k) dk. \qquad (13)$$

Furthermore, we define an occupation-number spectrum $n(k)$ via

$$N = \int |\hat{\psi}|^2 d^2k \equiv \int n(k) dk, \qquad (14)$$

where we denote the Fourier transform of $A(x)$ by $\hat{A}$; and, for notational convenience, we do not show explicitly the dependence of these spectra on time $t$. In any computational study, we must limit the number of Fourier modes that we use in our study of the GP equation; we refer to such a GP equation as a Fourier-truncated GP equation (cf [49, 50] for studies of the Fourier- or Galerkin-truncated Euler equation).

The Bogoliubov dispersion relation $\omega_B(k)$ is obtained by linearizing (1) around a constant $\psi$. For total number of particles (3) $N=1$, it is

$$\omega_B(k) = kc\sqrt{1 + \frac{\xi^2 k^2}{2}}, \qquad (15)$$

where the sound velocity is $c = \frac{\sqrt{2g}}{L}$ and the coherence length is

$$\xi = \frac{L}{\sqrt{g}}. \qquad (16)$$
We investigate thermalization in the 2D GP equation, so it is useful to recall that a uniform, interacting, 2D Bose gas has a high-temperature disordered phase and a low-temperature BKT phase [58–61], which shows quasi-long-range order with an algebraic decay of the spatial correlation function [53]

$$c(r) = \langle [e^{-i\theta(\mathbf{x})} - \langle e^{-i\theta(\mathbf{x})} \rangle][e^{i\theta(\mathbf{x}+\mathbf{r})} - \langle e^{i\theta(\mathbf{x}+\mathbf{r})} \rangle] \rangle \quad (17)$$

for temperatures T below the transition temperature $T_{BKT}$ (or energy $E_{BKT}$ in the microcanonical ensemble)

$$c(r) \sim r^{-\eta}, \quad (18)$$

where $r \equiv |\mathbf{r}|$ and the critical exponent $\eta < 0.25$ for $T < T_{BKT}$; and $\eta = 0.25$ at $T = T_{BKT}$ [54]. The BKT phase shows bound vortex–antivortex pairs; these unbind above $T_{BKT}$, so

$$c(r) \sim e^{-r/\ell} \quad (19)$$

in the disordered phase, with $\ell$ the correlation length.

## 2.2. Numerical methods and initial conditions

To perform a systematic, pseudospectral DNS of the spatiotemporal evolution of the 2D, Fourier-truncated, GP equation, we have developed a parallel message passing interface (MPI) code in which we discretize $\psi(\mathbf{x}, t)$ on a square simulation domain of side $L = 32$ with $N_c^2$ collocation points. We use periodic boundary conditions in both spatial directions, because we study homogeneous, isotropic turbulence in this 2D system, and a fourth-order, Runge–Kutta scheme, with time step $\Delta t$, for time marching. We evaluate the linear term in (1) in Fourier space and the nonlinear term in physical space; for the Fourier-transform operations we use the FFTW library (www.fttw.org). Thus, the maximum wave-number $k_{\max} = (N_c/2)\Delta k$, where $\Delta k = 2\pi/L$, and

$$\xi k_{\max} = \frac{\pi N_c}{\sqrt{g}}. \quad (20)$$

We have checked that, for the quantities we calculate, dealiasing of our pseudospectral code does not change our results substantially; here we present the results from our pseudospectral simulations that do not use dealiasing. For a general reference on numerical methods for quantum fluids, see [62].

To initiate turbulence in the 2D GP equation we use three types of initial conditions IC1 [41], IC2 and IC3 [52], always normalized to correspond to total number of particles (3) $N = 1$. The first of these is best represented in Fourier space as follows:

$$\hat{\psi}(\mathbf{k}, t=0) = \frac{1}{\sqrt{\pi^{1/2}\sigma}} \exp\left(-\frac{(k-k_0)^2}{2\sigma^2}\right) \exp(i\Theta(k_x, k_y)), \quad (21)$$

where $k = \sqrt{k_x^2 + k_y^2}$, $\Theta(k_x, k_y)$ are random numbers distributed uniformly on the interval $[0, 2\pi]$; $k_0 = N_0\Delta k$ and $\sigma = B\Delta k$, where the integer $N_0$ controls the spatial scale at which energy is injected into the system, and the real number B specifies the Fourier-space width of $\hat{\psi}$ at time $t = 0$. The initial condition IC2 is like IC1 but, in addition, it has a finite initial condensate population $N_0^i = |\hat{\psi}(\mathbf{k}=0, t)|^2 (\Delta k)^2$ at time $t = 0$. Note the study of [41] uses a hyper-viscosity term $\nu(-\nabla^2)^n \psi$, which is absent in our study; such hyperviscosity terms can
modify energy spectra in important ways, as has been discussed in the context of turbulence in the Navier–Stokes equation in [50, 63].

We obtain the initial condition IC3 by solving the 2D, stochastic, Ginzburg–Landau equation (SGLE), which follows from the free-energy functional

$$ \mathcal{F} = \int_{\mathcal{A}} d^2x \left( |\nabla \psi|^2 - \mu |\psi|^2 + \frac{1}{2} g |\psi|^4 \right), \qquad (22) $$

where $\mu$ is the chemical potential$^4$. The SGLE is

$$ \frac{\partial \psi}{\partial t} = -\frac{\delta \mathcal{F}}{\delta \psi^*} + \zeta(\mathbf{x}, t), \qquad (23) $$

where $\zeta$ is a zero-mean, Gaussian white noise with

$$ \langle \zeta(\mathbf{x}, t) \zeta^*(\mathbf{x}', t') \rangle = D \delta(\mathbf{x} - \mathbf{x}') \delta(t - t'), \qquad (24) $$

where $D = 2T$, in accordance with the fluctuation–dissipation theorem [64], $T$ is the temperature and $\delta$ the Dirac delta function. Finally, the SGLE (23) becomes

$$ \frac{\partial \psi}{\partial t} = \nabla^2 \psi - \mu \psi + g |\psi|^2 \psi + \zeta, \qquad (25) $$

which we solve along with the following ad hoc equation

$$ \frac{d\mu}{dt} = -\frac{\nu_N}{\mathcal{A}}(N - N_{\text{av}}), \qquad (26) $$

to control the number of particles $N$; the parameter $N_{\text{av}}$ controls the mean value of $N$; and $\nu_N$ governs the rate at which the SGLE equilibrates. We solve the SGLE by using a pseudospectral method, similar to the one described above for the 2D GP equation, with periodic boundary conditions in space, an implicit-Euler scheme, with time step $\Delta t$, for time marching and the method of reference [65] (see page 25 of this reference).

The motivation for choosing the parameters for our runs is to explore the initial-condition dependence of our results by varying the energy, the value $k_0$ of the wave number at which the initial energy is concentrated, the spread $\sigma$ of the energy about $k_0$ at time $t = 0$, the interaction strength $g$ and the presence of the initial condensate density $N_0^i$. For our systematic study, we have performed numerous runs; these are listed in table 1; the parameters for these runs have been chosen to highlight one type of behaviour or another; e.g. the time dependence of the self-truncation wave number $k_c$ (see below) changes dramatically as we change the parameters of our runs (see table 2); this exploration of parameter space is necessary in order to gain a comprehensive understanding of the dynamics of the 2D, Fourier-truncated, GP equation.

Note that we introduce the initial conditions IC2 and IC3 to obtain the relevant behaviour in a shorter time span than is possible with IC1. The initial condition IC1 has zero initial condensate population $N_0^i$ and, therefore, several vortices; its dynamics involves a build up of the condensate population; this takes a long time. The initial condition IC2 is similar to IC1 but with an initial condensate population $N_0^i > 0$; thus, condensate build-up is bypassed. In the same manner, IC3 allows us to study the late stages of the self-truncation regime.

<footnote>$^4$ Recall that the SGLE can be thought of as an imaginary-time GP equation with external, additive noise (see e.g. [38]).</footnote>
<table><caption>**Table 1.** Parameters for our DNS runs A1–A13, B1–B2 and C1–C6: $N_c^2$ is the number of collocation points, $k_0$ is the energy-injection scale, $\sigma$ is the Fourier-space width of $\hat{\psi}$ at $t=0$; $g$ is the effective interaction strength; $N_0^i$ is the initial condensate population; $D$ and $k_c^{in}$ are, respectively, the variance of the white noise and the initial value of the truncation wave number, which we use in the initial conditions of type IC3; $E$ is the total energy; we use a square simulation domain of area $A=L^2$; we choose $L=32$.</caption><thead><tr><th></th><th>N<sub>c</sub></th><th>k<sub>0</sub> (×Δk)</th><th>σ (×Δk)</th><th>g</th><th>N<sub>0</sub><sup>i</sup></th><th>√D (×10<sup>-3</sup>)</th><th>k<sub>c</sub><sup>in</sup></th><th>E</th></tr></thead><tbody><tr><td>A1</td><td>1024</td><td>5</td><td>2</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>2.120</td></tr><tr><td>A2</td><td>1024</td><td>5</td><td>2</td><td>2000</td><td>—</td><td>—</td><td>—</td><td>3.045</td></tr><tr><td>A3</td><td>1024</td><td>5</td><td>2</td><td>5000</td><td>—</td><td>—</td><td>—</td><td>5.82</td></tr><tr><td>A4</td><td>1024</td><td>35</td><td>5</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>49.69</td></tr><tr><td>A5</td><td>512</td><td>5</td><td>2</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>2.15</td></tr><tr><td>A6</td><td>256</td><td>5</td><td>2</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>2.07</td></tr><tr><td>A7</td><td>128</td><td>5</td><td>2</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>2.1</td></tr><tr><td>A8</td><td>64</td><td>5</td><td>2</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>2.2</td></tr><tr><td>A9</td><td>256</td><td>5</td><td>2</td><td>2000</td><td>—</td><td>—</td><td>—</td><td>2.94</td></tr><tr><td>A10</td><td>256</td><td>5</td><td>2</td><td>5000</td><td>—</td><td>—</td><td>—</td><td>5.57</td></tr><tr><td>A11</td><td>256</td><td>15</td><td>2</td><td>1000</td><td>—</td><td>—</td><td>—</td><td>9.86</td></tr><tr><td>A12</td><td>256</td><td>15</td><td>2</td><td>2000</td><td>—</td><td>—</td><td>—</td><td>10.82</td></tr><tr><td>A13</td><td>256</td><td>15</td><td>2</td><td>5000</td><td>—</td><td>—</td><td>—</td><td>13.68</td></tr><tr><td>B1</td><td>128</td><td>5</td><td>1</td><td>10000</td><td>0.95</td><td>—</td><td>—</td><td>5.44</td></tr><tr><td>B2</td><td>128</td><td>5</td><td>1</td><td>1000</td><td>0.95</td><td>—</td><td>—</td><td>0.59</td></tr><tr><td>C1</td><td>256</td><td>—</td><td>—</td><td>5000</td><td>—</td><td>8</td><td>6</td><td>2.536</td></tr><tr><td>C2</td><td>256</td><td>—</td><td>—</td><td>1000</td><td>—</td><td>8</td><td>6</td><td>0.583</td></tr><tr><td>C3</td><td>256</td><td>—</td><td>—</td><td>1000</td><td>—</td><td>10</td><td>6</td><td>0.637</td></tr><tr><td>C4</td><td>256</td><td>—</td><td>—</td><td>1000</td><td>—</td><td>8</td><td>9</td><td>0.7</td></tr><tr><td>C5</td><td>256</td><td>—</td><td>—</td><td>1000</td><td>—</td><td>8</td><td>15</td><td>1.085</td></tr><tr><td>C6</td><td>256</td><td>—</td><td>—</td><td>1000</td><td>—</td><td>8</td><td>20</td><td>1.557</td></tr></tbody></table>

### 3. Results

We first present the time evolution of the different energies, the PDFs of the velocity components and the population $N_0$ in the zero-wave-number mode. We then give a detailed statistical characterization of the temporal evolution of the Fourier-truncated, 2D GP equation in the four regimes mentioned in the introduction (section 1).

#### 3.1. Evolution of energies, velocity probability distribution functions and the zero-wave-number population

We show the early stages of the time evolution of the energies $E_{kin}^i$, $E_{kin}^c$, $E_{int}$ and $E_q$, from our DNS runs A1–A4, B1 and C6 in figure 1. The runs A1–A4 use initial conditions of type IC1, in which $E_{kin}^i$ is a significant fraction of the total initial energy; the runs B1 and C6 start with initial configurations of types IC2 and IC3, respectively, in which $E_{kin}^i$ is negligibly small at $t=0$. The transient nature of the early stages of the dynamical evolution of the dissipationless, unforced,
<table><caption>**Table 2.** Summary of the self-truncation results from our DNS runs A1–A4, B2 and C1–C6: *E* is the total energy; *k*<sub>max</sub> = 2π*N*<sub>c</sub>/2*L*; *ξ* = *L*/√*g* is the healing length; *k*<sub>c</sub><sup>i</sup> and *k*<sub>c</sub><sup>f</sup> are the initial and final values of *k*<sub>c</sub> (averaged over a few time steps); *α*<sub>1</sub> is the slope obtained from the log-log (base 10) plot of *k*<sub>c</sub> versus *t* and *α*<sub>2</sub> = 1/(1 − *χ*), where *χ* is the slope obtained from the log-log (base 10) plot of *dk*<sub>c</sub>/dt versus *k*<sub>c</sub>/*k*<sub>max</sub>.</caption><thead><tr><th></th><th><i>E</i></th><th><i>k</i><sub>max</sub></th><th>ξ</th><th>ξ<i>k</i><sub>max</sub></th><th><i>k</i><sub>c</sub><sup>i</sup></th><th><i>k</i><sub>c</sub><sup>f</sup></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th></tr></thead><tbody><tr><td>A1</td><td>2.120</td><td>100.53</td><td>1.01</td><td>101.73</td><td>4.52</td><td>12.42</td><td>0.28</td><td>0.26</td></tr><tr><td>A2</td><td>3.045</td><td>100.53</td><td>0.72</td><td>71.9</td><td>5.39</td><td>18.72</td><td>0.28</td><td>0.28</td></tr><tr><td>A3</td><td>5.82</td><td>100.53</td><td>0.45</td><td>45.49</td><td>7.11</td><td>31.3</td><td>0.29</td><td>0.27</td></tr><tr><td>A4</td><td>49.69</td><td>100.53</td><td>1.01</td><td>101.73</td><td>17.31</td><td>30.53</td><td>0.2</td><td>0.21</td></tr><tr><td>B2</td><td>0.589</td><td>12.57</td><td>1.01</td><td>12.72</td><td>2.23</td><td>9.23</td><td>0.24</td><td>0.25</td></tr><tr><td>C1</td><td>2.536</td><td>25.13</td><td>0.45</td><td>11.37</td><td>7.08</td><td>19.91</td><td>0.22</td><td>0.22</td></tr><tr><td>C2</td><td>0.583</td><td>25.13</td><td>1.01</td><td>25.43</td><td>6.15</td><td>8.90</td><td>0.12</td><td>0.14</td></tr><tr><td>C3</td><td>0.637</td><td>25.13</td><td>1.01</td><td>25.43</td><td>6.18</td><td>10.05</td><td>0.14</td><td>0.15</td></tr><tr><td>C4</td><td>0.6999</td><td>25.13</td><td>1.01</td><td>25.43</td><td>9.05</td><td>11.07</td><td>0.09</td><td>—</td></tr><tr><td>C5</td><td>1.085</td><td>25.13</td><td>1.01</td><td>25.43</td><td>15.09</td><td>16.08</td><td>0.04</td><td>—</td></tr><tr><td>C6</td><td>1.557</td><td>25.13</td><td>1.01</td><td>25.43</td><td>20.17</td><td>20.87</td><td>0.02</td><td>—</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Plots versus time *t* of the four components of the total energy *E*<sup>i</sup><sub>kin</sub>, *E*<sup>c</sup><sub>kin</sub>, *E*<sub>int</sub> and *E*<sub>q</sub>, during the initial stages of evolution, from our DNS runs (a) A1, (b) A2, (c) A3, (d) A4, (e) B1 and (f) C6 (see table 1).</figcaption></figure>

2D GP equation is evident from figure 1, in which we observe a rapid conversion of *E*<sup>i</sup><sub>kin</sub> into the other three components, with a significant fraction being transferred to *E*<sup>c</sup><sub>kin</sub>; moreover, the transient stage depends on the initial conditions, as we describe below. Figures 1(a)–(c)
show comparisons of the temporal evolution of the energies, from the runs A1–A3; we observe, in particular, that the conversion of $E_{\text{kin}}^i$ into the other energy components is accelerated as $g$ increases from 1000 to 5000 (cf [42]); and there is a corresponding acceleration in the approach to thermalization. The time evolution of the incompressible kinetic energy follows that of the total number of vortices $N_v$; see e.g. [66], and figure F1, which we have included in the supplementary material (available from stacks.iop.org/NJP/15/113025/mmedia). Moreover, the larger the value of $E_{\text{kin}}^i$, the larger is the time required for thermalization, as we can see by comparing figures 1(a) and (d), for the runs A1 and A4, respectively; the run A4 starts with a high value of $E_{\text{kin}}^i$ ($t=0$) because of a large number of vortices and anti-vortices, so it takes a long time to thermalize; indeed, if the spatial resolution of our DNS is very high, the computational cost of achieving a statistically steady state is prohibitively high for initial conditions A1–A4. In contrast, the runs B1 and C6 have negligibly small values of $E_{\text{kin}}^i$ ($t=0$) to begin with (figures 1(e) and (f), respectively); and $E_{\text{kin}}^i(t)$ remains close to zero throughout the dynamical evolution here. For run B1, both $E_{\text{kin}}^c$ and $E_q$ start from values close to zero, grow at the cost of $E_{\text{int}}$ and finally saturate to small, statistically steady values. For run C6, there are hardly any vortices in the initial configuration, so the energies start fluctuating about their statistically steady values very rapidly.

In figure 2 we plot, at three instants of time, the PDFs of $v_x$ and $v_y$, the Cartesian components of the velocity, for our DNS runs A1, B1 and C6, which correspond, respectively, to initial conditions of types IC1, IC2 and IC3. For the run A1, these PDFs, in figures 2(a)–(c), show a cross-over from a distribution with power-law tails to one that is Gaussian; the right and left tails of the PDFs in figure 2(a) can be fitted to the form $\sim v_i^{-\gamma}$, with $\gamma \simeq 3.2$, and $i=x$ or $y$ (we show fits only for $i=x$). Such power-law tails in velocity-component PDFs have been seen in experiments [67] and some numerical studies [39, 45, 68, 69]. However, it has not been noted hitherto that, for turbulence in the Fourier-truncated, 2D GP equation with low-energy initial conditions, such PDFs evolve, as $t$ increases, from PDFs with power-law tails (figure 2(a) for run A1), to ones with a Gaussian form near the mean, followed by broad tails (figure 2(b) for run A1), and then to more-or-less Gaussian PDFs (figure 2(c) for run A1), but with tails that can be fitted to an exponential form. This evolution towards Gaussian PDFs is associated with the annihilation of vortices and anti-vortices. The video S1 in the supplementary material shows the temporal evolution of this PDF in the left panel and the spatiotemporal evolution of the pseudocolour plot of the vorticity in the right panel. The analogues of figures 2(a)–(c) for runs B1 and C1, both of which have a negligibly small value of $E_{\text{kin}}^i$ at $t=0$, are given, respectively, in figures 2(d)–(f) and (g)–(i).

To calculate the velocity PDFs, during the various stages of the evolution of the system, we obtain the velocity at every grid point in our simulation domain; thus, we make these numerical measurements at a length scale that is always less than the inter-vortex separation. The power-law tails, which we observe in the velocity PDFs for our DNS run A1, arise because of the singular nature of quantum vortices [12]; and the cross-over from such power-law tails, in the initial stages of evolution, to the more-or-less Gaussian PDFs, in the partially thermalized state, arises because of the depletion of the vortex density with time. Thus, our results complement those in [70] insofar as the cross-over from power law to Gaussian tails occurs as our system evolves in time and not as we change the length scale of our measurement as in [70].

We turn now to the time evolution of the population $N_0(t)$, in the $k=0$ mode [36, 40, 71], and its dependence on the initial conditions. In figure 3(a) we plot $N_0$ versus $t$ for the runs A1–A4 (red, blue, green and brown curves, respectively), which use initial configurations of type IC1;
<figure><img src="image_3.png" /></figure>

**Figure 2.** Semilog (base 10) plots of the PDFs of the *x* (red circles) and *y* (green squares) components of the velocity from our DNS runs: (a)–(c) A1, (d)–(f) B1 and (g)–(i) C6, corresponding to each of the three types of initial conditions IC1, IC2 and IC3, respectively. The complete time evolution of the PDFs (a)–(c) for the run A1 is illustrated in the top-left panel of the video S1; see the supplementary material available from stacks.iop.org/NJP/15/113025/mmedia. The blue-dashed lines in (b)–(i) indicate fits to Gaussian PDFs; the dashed lines in (a) indicate power-law fits to the left (blue-dashed line) and right (orange-dashed line) tails of the PDFs (see text).

these figures show that $N_0(t)$ increases with $t$, on average, and depends on $E$, $g$, $k_0$ and $\sigma$. For the runs A1 and A2 (red and blue curves in figure 3(a)), $N_0(t)$ approaches a saturation value for the time scales probed by our simulations; figure 3(a) also shows that, as we increase $g$ (red, blue and green lines in figure 3(a)), the fluctuations in $N_0$ are enhanced and its large $t$ value, which it seems to approach asymptotically, diminishes. By comparing the runs A1 and A4 (red and brown lines in figures 3(a)), we see that the latter has a higher value of $E$ than the former, because both $k_0$ and $\sigma$ are smaller for A1 than for A4; thus, $N_0(t)$ grows more slowly in A4 than in A1; and, after equal lengths of simulation time, its value in A4 is nearly an order of magnitude
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Plots versus time *t* of the population *N*₀, in the zero-wave-number mode, from our DNS runs (a) A1–A4 (initial condition of type IC1), (b) B1 (initial condition of type IC2), (c) C6 (initial condition of type IC3) and (d) A1 and A5–A8, for five values for the number of collocation points *N*<sub>c</sub><sup>2</sup>, namely, 1024<sup>2</sup>, 512<sup>2</sup>, 256<sup>2</sup>, 128<sup>2</sup> and 64<sup>2</sup>.</figcaption></figure>

lower than in A1; the former shows large fluctuations in *N*₀(*t*) and no sign of saturation. The run B1 (figure 3(e)) uses an initial configuration of type IC2, with a large value of *N*₀(*t* = 0) = 0.95; in this case, after a period of initial transients, *N*₀(*t*) → 0.98 over our simulation time. The run C6 (figure 3(f)) uses an initial condition of type IC3; here *N*₀(*t*) fluctuates slightly but remains close to its initial value (cf [40, 71]).

To study the dependence of $N_0(t)$ on the number of collocation points $N_c^2$, we evolve the initial configuration of A1 for $N_c = 512$ (run A5), 256 (run A6), 128 (run A7) and 64 (run A8). Figure 3(g) shows plots of $N_0(t)$ versus $t$ for these five runs; clearly, the initial evolution of $N_0(t)$ depends significantly on $N_c$; however, the large $t$ values of $N_0(t)$, on the time scales of our runs, are comparable ($\simeq 0.9$) for the runs with $N_c = 128$ (run A7), 256 (run A6) and 1024 (run A1). In contrast, the saturation value for the run with $N_c = 64$ (run A8) is $\simeq 0.8$. For the run A5 ($N_c = 512$), $N_0(t)$ shows large fluctuations and no sign of saturation over the time scale that we have covered; this suggests that $N_0(t)$ also depends on the realization of the random phases $\Theta(k_x, k_y)$ in (21). These plots of $N_0(t)$ illustrate that complete thermalization proceeds very slowly for $N_0$; in the completely thermalized state of the Fourier-truncated, 2D GP system, $N_0$ must vanish in the thermodynamic limit by virtue of the Hohenberg–Mermin–Wagner theorem [58, 59]; however, it is not easy to realize this limit in finite-size systems and with the
<figure><img src="image_5.png" /></figure>

**Figure 4.** Log-log (base 10) plots of the spectra $E_{\text{kin}}^{i}(k)$ from our DNS runs (a) A1, (b) A2, (c) A3 and (d) A4 at different times $t$ (indicated by curves of different colours); a $k^{-3}$ power law is shown by orange-dashed lines. The complete time evolution of the spectra in (a)–(d) is illustrated in the video S2 (see the supplementary material, available from stacks.iop.org/NJP/15/113025/mmedia).

limited run times that are dictated by computational resources. We discuss these issues again in section 3.4 and also refer the reader to [71, 72].

## 3.2. Initial transients and the onset of thermalization

The initial stages of the evolution of energy spectra for the Fourier-truncated, 2D GP equations are qualitatively different for initial conditions of types IC1–IC3. The first type begins with a sizeable incompressible kinetic energy spectrum $E_{\text{kin}}^{i}(k)$; and the initial transients are associated with the annihilation and creation of vortex–antivortex pairs, the associated depletion of $E_{\text{kin}}^{i}(k)$, and the growth of the other energy components [41]. In contrast, runs with initial conditions of types IC2 and IC3 start with a very small incompressible-energy component, therefore, even the early stages of their dynamical evolution are akin to the late stages of the dynamical evolution with initial conditions of type IC1. In figures 4(a)–(d) we show the time evolution of the spectra $E_{\text{kin}}^{i}(k)$, for the runs A1–A4, to ascertain the presence of scaling behaviour, if any. We find that, in the low-$k$ region, $E_{\text{kin}}^{i}(k)$ lacks a well-defined scaling region (unlike in [42]); indeed, this region depends on the initial configuration, changes continuously with time, and, in particular,
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Log-log (base 10) plots of the spectra $E_{\text{kin}}^{\text{c}}(k)$ from our DNS runs (a)–(c) A1, (d)–(f) A4 and (g)–(i) B1 at different times $t$ (indicated by curves of different colours); a $k$ power law is shown by orange-dashed lines.</figcaption></figure>

a $k^{-5/3}$ scaling region is tenable (a) over a range of wave numbers that is very tiny and (b) over a fleetingly short interval of time (around $t=50$ for the run A1). At large wave numbers, $E_{\text{kin}}^{\text{i}}(k) \sim k^{-3}$, during the initial stages of evolution, because of the presence of the vortices [44]; this power-law form holds over the same time scales for which the PDF $P(v_x/\sigma_{v_x}) \sim v_x^{-\gamma}$ (figures 2(a) and (b)).

It is useful to define a wave number $k_{hd} = 2\pi/\delta$ based on the average vortex separation $\delta = \lambda^{-1/2}$, where $\lambda = N_v/A$ is the vortex density. We calculate $k_{hd}$ and plot it versus time $t$, for the DNS runs A1–A3, in figure F2 in the supplementary material (available from stacks.iop.org/NJP/15/113025/mmedia). In our simulations, only a small number of modes have $k \le k_{hd}$, especially for our DNS runs A1–A3 (see figure 4 in which we have plotted $E_{\text{kin}}^{\text{i}}(k)$). In our study, given the system sizes and initial conditions we use, the system evolves towards states in which the vortex density is low.

The initial transients described above are followed by a regime in which the energy and occupation-number spectra begin to show power-law-scaling behaviour, but the power-law exponent and the extent of the scaling region change with time and depend on the initial conditions; we regard this as the onset of thermalization, which is shown in figures 5 and 6,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Log-log (base 10) plots of the spectra $E_{\text{kin}}^{\text{c}}(k)$ from our DNS runs (a)–(c) A7 and (d)–(f) B2 at different times $t$ (indicated by curves of different colours); a $k$ power law is shown by orange-dashed lines.</figcaption></figure>

where we illustrate the time evolution of $E_{\text{kin}}^{\text{c}}$. Figure 5(a) shows $E_{\text{kin}}^{\text{c}}(k)$ for the run A1; we begin to see a power-law region here with $E_{\text{kin}}^{\text{c}}(k) \sim k$, on the low-$k$ side of the peak after which the spectrum falls steeply. Similar $E_{\text{kin}}^{\text{c}}(k) \sim k$ behaviour starts to emerge in the region $k \lesssim k_{\text{max}}$ for the run B1 (figure 5(g)). In this onset-of-thermalization regime, we also see the development of the following power laws: $E_{\text{int}}(k) + E_q(k) \sim k$ (figure 7) and $n(k) \sim 1/k$ (figure 8).

### 3.3. Partial thermalization and self-truncation

#### 3.3.1. Partial thermalization.
In the third stage of the dynamical evolution of the 2D, Fourier-truncated, GP equation, which we refer to as the partial-thermalization stage, well-defined, power-law-scaling behaviour appears in energy and occupation-number spectra, with exponents that are independent of the initial conditions as illustrated by the compressible-kinetic-energy spectra in figures 5(b), (c) (e), (f) and 6(b) for initial conditions of type IC1, and figures 5(h) and 6(e), (f), for initial conditions of type IC2. It is important to distinguish between (I) spectra that fall steeply at large values of $k$, e.g. the spectra in figures 5(b), (c) (e), (f) and 6(e), (f), and (II) spectra that increase all the way to $k_{\text{max}}$, e.g. the spectra in figures 5(h) and 6(b), (c). In case (I), we have spectral convergence to the 2D GP partial differential equation (PDE); in case (II), the effects of Fourier truncation are so pronounced that our truncated 2D GP system does not provide a good representation of the 2D GP PDE. As we show below, case (I) can be further subdivided into (A) a subclass in which the maximum, at $k = k_c$ in $E_{\text{kin}}(k) = (E_{\text{kin}}^{\text{c}}(k) + E_{\text{kin}}^{\text{i}}(k))$, referred to as the self-truncation wave number [52], moves out to $k_{\text{max}}$ as a power of $t$ and (B) a subclass in which $k_c$ moves out to $k_{\text{max}}$ at a rate that is slower than a power of $t$.

Figures 5(g)–(i), from the run B1, show how $E_{\text{kin}}^{\text{c}}(k)$ evolves as the spectral convergence to the GP PDE is lost in case (II); note that the scaling region with $E_{\text{kin}}^{\text{c}} \sim k$ sets in at high
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Log-log (base 10) plots of the spectra $E_{\text{int}}(k) + E_q(k)$ from our DNS runs (a)–(c) A1, (d)–(f) A4 and (g)–(i) B1 at different times $t$ (indicated by curves of different colours); a $k$ power law is shown by orange-dashed lines.</figcaption></figure>

wave numbers close to $k_{\text{max}}$ and then extends to the low-wave-number regime. For case (IA), analogous plots of $E_{\text{kin}}^c(k)$ are given in, e.g., figures 6(a)–(c). We give plots for case (IB) in the next subsection, where we study in detail the time dependence of $k_c$. Illustrative plots of the spectra ($E_i(k) + E_q(k)$) and $n(k)$ in this regime of partial thermalization are given in figures 7 and 8, respectively.

### 3.3.2. Self-truncation.
We now present a detailed characterization of the partial-thermalization regime, when energy spectra display self-truncation at wave numbers beyond $k_c(t)$, which can be defined as follows:

$$k_c = \sqrt{\frac{2 \int_0^{k_{\text{max}}} k^2 E_{\text{kin}}(k) dk}{\int_0^{k_{\text{max}}} E_{\text{kin}}(k) dk}}, \qquad (27)$$

as the system approaches complete thermalization, $k_c(t) \to k_{\text{max}}$. In particular, we explore how the scaling ranges in energy spectra grow with $t$ for different values of $g$, with the initial
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

configuration and number of collocation points $N_c$ held fixed. For an initial condition of type IC1, with $k_0 = 5\Delta k$, $\sigma = 2\Delta k$ and $N_c = 256$, we obtain the time evolution of energy spectra for $g = 1000$ (run A6), 2000 (run A9) and 5000 (run A10) in figures 9(a)–(c), respectively, and their video analogues (video S3 (panel V2) in the supplementary material, available from stacks.iop.org/NJP/15/113025/mmedia). The larger the value of $g$, the more rapid is the thermalization, and the consequent loss of spectral convergence, as we can see by comparing the sky-blue (run A10), green (run A9) and purple (run A6) spectra in figures 9(a)–(c); run A6 loses spectral convergence around $t = 2500$. We obtain the same qualitative $g$ dependence, with $k_0 = 15\Delta k$, $\sigma = 2\Delta k$ and $N_c = 256$, for $g = 1000, 2000$ and $5000$, i.e. runs A11, A12 and A13, respectively, for which energy spectra are portrayed in figures 9(d)–(f) and video S3 (panel V3) in the supplementary material.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

In figures 9(g)–(i) we explore the $N_c$ dependence of the self-truncation of energy spectra, for initial conditions, with $k_0 = 5\Delta k$, $\sigma = 2\Delta k$ and $g = 1000$, and five different values of $N_c$, namely, $N_c = 1024$ (run A1), 512 (run A5), 256 (run A6), 128 (run A7) and 64 (run A8). We find, not surprisingly, that the lower the value of $N_c$, the more rapidly the system loses spectral convergence.

Note the dual nature of solutions to the truncated GPE: in the early part of the dynamical evolution of this system, which lasts only as long as spectral convergence is ensured, the solutions approximate numerically genuine solutions of the original PDE (i.e. the untruncated GPE). In subsequent evolution, when spectral convergence is lost, the truncated system evolves to a truncation-dependent thermodynamic equilibrium, which we call ‘complete thermalization’. These two distinct aspects might appear, at first sight, to be mutually exclusive. Indeed, in the first part of the evolution, the truncation wave number $k_{max}$ has no effect (or

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Log-log (base 10) plots of the spectra $E_{\text{kin}}(k)$ from our DNS runs (initial conditions of type IC3): (a) C1, (b) C2, (c) C3, (d) C4, (e) C5 and (f) C6.</figcaption></figure>

a vanishingly small one) on the solution, whereas, in the next part, the very existence of the equilibrium needs the influence of $k_{\text{max}}$. However, figure 9 shows that, in a regime that we call ‘partial thermalization’, the system ‘self-truncates’ with a physical cutoff at momentum $k_c < k_{\text{max}}$. In this new regime the system is in a state that both approximates a solution to the original PDE and can also be thought of as a thermalized state with a slowly growing ‘self-truncation’ at wave number $k_c$. We show below that the self-truncation regime can last a long time when $\xi k_c$ is large.

Initial conditions of type IC2 lead to energy spectra whose time evolution, and their dependence on $g$ and $N_c$, is similar to those that are obtained from initial conditions of type IC1.

With initial conditions of types IC1 and IC2, we cannot control the initial value $k_c(t=0) = k_c^{\text{in}}$ easily. However, initial conditions of type IC3, which we obtain from the SGLE, allow us to control $k_c^{\text{in}}$ and start, therefore, with initial spectra that display partial thermalization for $k < k_c^{\text{in}}$ [52] and a sharp fall thereafter. In figure 10, we show the time evolution of $E_{\text{kin}}^c(k)$ for such initial conditions from runs C1–C6. For different representative values of $k_c^{\text{in}}$, $g$ and $D$, we now study the time evolution of $k_c(t)$, which characterizes the growth of the partially thermalized scaling region. Here too, as with initial conditions of types IC1 and IC2, if all other parameters like $k_c^{\text{in}} = 6.0$ and $D$ are held fixed, the speed of thermalization increases with $g$ (cf figure 10(a) for the run C1, with $g = 5000$, and figure 10(b) for the run C2, with $g = 1000$). For these runs C1–C6, the growth of the energy spectra, in the region $k > k_c^{\text{in}}$, starts with the smoothing of the sharp cutoff at $k_c^{\text{in}}$; the higher the value of $k_c^{\text{in}}$, the slower is this growth (cf figures 10(b), (d), (e) and (f) for runs C2, C4, C5 and C6, respectively). By contrast, an increase in $D$ (or $T$) in the SGLE accelerates this growth (cf figures 10(b) and (c) for runs C2 and C3, respectively).
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Plots of (a) the self-truncation wave-number $k_c(t)$ versus time $t$ and (b) $dk_c/dt$ versus $k_c/k_{max}$ from our DNS runs A1–A4, B2 and C1–C6.</figcaption></figure>

The growth of $k_c(t)$ with $t$, illustrated in figure 11(a), can be fitted to the form $k_c(t) \sim t^\alpha$; however, as we show below, $\alpha$ depends on the initial condition. We obtain the exponent $\alpha$ either from slopes of log-log plots of (i) $k_c(t)$ versus $t$ or (ii) $dk_c/dt$ versus $k_c/k_{max}$; we denote the values from procedures (i) and (ii) as $\alpha_1$ and $\alpha_2$, respectively. Note that in (ii) we have a parametric plot [38, 52], shown in figure 11(b); this yields a straight-line scaling regime with slope $\chi$ and $\alpha_2 = 1/(1-\chi)$. The values of $\alpha_1$ and $\alpha_2$, listed in table 2, show that $\alpha_1 \approx \alpha_2$; the discrepancy between these two values for $\alpha$ is a convenient measure of the errors of our estimates. For runs C4, C5 and C6, we cannot obtain $\alpha_2$ reliably; the small values of $\alpha_1$ for these runs indicate very slow growth of $k_c(t)$; indeed, in runs C5 and C6, a case can be made for a logarithmic growth of $k_c(t)$ with $t$.

### 3.4. Complete thermalization

The partially thermalized stage of the dynamical evolution of the 2D, Fourier-truncated, GP equation may either gradually become completely thermalized, in which state a power-law scaling region is present in the entire energy and the occupation number spectra, or remain self-truncated with logarithmic growth. In figures 5(g)–(i) and 6(a)–(c), we show the compressible kinetic energy spectra $E_{\text{kin}}^c$ for the runs B1 and A7, where $E_{\text{kin}}^c$ shows power-law scaling over the entire wave-number range, from $k = 2\pi/L$ up to $k_{\text{max}}$, towards the end of the respective simulations; a naive fit is consistent with $E_{\text{kin}}^c(k) \sim k$ (but see below).

#### 3.4.1. Correlation functions and the Berezinskii–Kosterlitz–Thouless (BKT) transition.
A uniform, 2D, interacting Bose gas exhibits a BKT phase at low energies (temperatures in the canonical ensemble). Thus, the completely thermalized state of the 2D, Fourier-truncated, GP equation should yield a BKT phase [53, 55], with the correlation function $c(r) \sim r^{-\eta}$, at energies $E < E_{\text{BKT}}$; and $c(r)$ should decay exponentially with $r$ if $E > E_{\text{BKT}}$. We show this explicitly now by using initial conditions of type IC1 with $N_c = 64$ and 128, and $g = 1000$; we obtain different energies by changing $k_0$ and $\sigma$ (runs D1–D13 and E1–E12 in table 3).

In figure 12, we present plots of the correlation functions $c(r)$. To illustrate the BKT transition clearly, we present log-log plots of $c(r)$ versus $r$, for $E < E_{\text{BKT}}$, in figures 12(a) and (d), where the straight lines indicate power-law regimes; and, for $E > E_{\text{BKT}}$, we use semi-log
<table><caption>**Table 3.** List of parameters for our complete-thermalization DNS runs D1–D13 ($N_c^2 = 128^2$) and E1–E12 ($N_c^2 = 64^2$): $N_c^2$ is the number of collocation points; $k_0$ is the energy-injection scale; $\sigma$ is Fourier-space width of $\psi$ at $t = 0$; $E$ is the total energy; and $\eta$ is the exponent of the correlation function $c(r) \sim r^{-\eta}$ for $E < E_{BKT}$. $g = 1000$ for all the DNS runs and they have been performed on a square simulation domain of area $A = L^2$, with $L = 32$.</caption><thead><tr><th rowspan="2">N<sub>c</sub> = 128</th><th colspan="2">k<sub>0</sub></th><th colspan="2">σ</th><th colspan="2">k<sub>0</sub></th><th colspan="2">σ</th></tr><tr><th>(×Δk)</th><th>(×Δk)</th><th>E</th><th>η</th><th>N<sub>c</sub> = 64</th><th>(×Δk)</th><th>(×Δk)</th><th>E</th><th>η</th></tr></thead><tbody><tr><td>D1</td><td>5</td><td>2</td><td>2.1</td><td>0.008</td><td>E1</td><td>0</td><td>2</td><td>1.12</td><td>0.012</td></tr><tr><td>D2</td><td>10</td><td>2</td><td>5.05</td><td>0.024</td><td>E2</td><td>3</td><td>2</td><td>1.64</td><td>0.025</td></tr><tr><td>D3</td><td>12</td><td>2</td><td>6.74</td><td>0.034</td><td>E3</td><td>5</td><td>2</td><td>2.2</td><td>0.040</td></tr><tr><td>D4</td><td>14</td><td>2</td><td>8.74</td><td>0.047</td><td>E4</td><td>8</td><td>2</td><td>3.68</td><td>0.083</td></tr><tr><td>D5</td><td>16</td><td>2</td><td>11.05</td><td>0.080</td><td>E5</td><td>10</td><td>2</td><td>5.04</td><td>0.164</td></tr><tr><td>D6</td><td>18</td><td>2</td><td>13.68</td><td>0.111</td><td>E6</td><td>11</td><td>2</td><td>5.84</td><td>0.255</td></tr><tr><td>D7</td><td>20</td><td>2</td><td>16.62</td><td>0.181</td><td>E7</td><td>12</td><td>2</td><td>6.75</td><td></td></tr><tr><td>D8</td><td>21</td><td>2.5</td><td>18.34</td><td>0.239</td><td>E8</td><td>13</td><td>2</td><td>7.74</td><td></td></tr><tr><td>D9</td><td>24</td><td>3</td><td>23.75</td><td></td><td>E9</td><td>14</td><td>2</td><td>8.78</td><td></td></tr><tr><td>D10</td><td>25</td><td>2</td><td>25.3</td><td></td><td>E10</td><td>15</td><td>2</td><td>9.88</td><td></td></tr><tr><td>D11</td><td>26</td><td>2</td><td>27.27</td><td></td><td>E11</td><td>16</td><td>2</td><td>11.05</td><td></td></tr><tr><td>D12</td><td>28</td><td>2</td><td>31.44</td><td></td><td>E12</td><td>17</td><td>2</td><td>12.32</td><td></td></tr><tr><td>D13</td><td>30</td><td>2</td><td>35.9</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

plots, as in figures 12(b) and (e), where the straight lines signify an exponential decay of $c(r)$ with $r$. Given the resolution of our DNS runs, we find that, in a small energy range in the vicinity of $E_{BKT}$, we cannot fit power-law or exponential forms satisfactorily; this leads to an uncertainty in our estimate for $E_{BKT}$. Apart from this uncertainty, the behaviour of $c(r)$, in the regime of complete thermalization, is in accord with our expectations for the BKT phase; in particular, the exponent $\eta$ (see equation (18)) depends on $E$ for $E < E_{BKT}$ as shown in figures 12(c) and (f). Our values for $\eta$, for the runs with $E < E_{BKT}$ and with $N_c = 64$ and 128, are listed in table 3. Note that $E_{BKT} \simeq 6$ ($N_c = 64$) and $E_{BKT} \simeq 19$ ($N_c = 128$), i.e. $E_{BKT}$ depends on $N_c$, the number of collocation points; we show analytically below how a low-temperature analysis can be used to understand this dependence of $E_{BKT}$ on $N_c$. In the completely thermalized state of the Fourier-truncated, 2D GP system, $N_0$ must vanish in the thermodynamic limit by virtue of the Hohenberg–Mermin–Wagner theorem [58, 59] and $n(k) \sim k^{-1+\eta}$; it is not easy to realize this limit in finite-size systems and with the limited run times that are dictated by computational resources (see the plots of $N_0$ in figure 3); however, finite-size scaling can be used to extract the exponent $\eta$ from the $k=0$ part of $n(k)$, as shown in reference [72]; similarly, $E_{kin}^c(k)$ should also show a power-law form with an exponent that depends on $\eta$, but this is difficult to realize in numerical calculations with limited spatial resolutions and run lengths.

### 3.4.2. Analytical estimation of the energy of the BKT transition.
The energy of a pure condensate of a uniform, weakly interacting, 2D Bose gas, which is described by the GP equation (1), is $E_0 = g/(2A)$. We define the energy of our system to be $E = E_0(1 + \delta E)$; this energy $E$ is fixed by the initial condition; and $\delta E$ measures the relative amount by which $E$
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Plots of $c(r)$ versus $r$ for different energies in the complete-thermalization regime, for $N_c^2 = 128^2$ ((a) and (b)) and $N_c^2 = 64^2$ ((d) and (e)). (a) and (d) Log-log (base 10) plots of $c(r)$ versus $r$ for different energies $E < E_{KT}$; the slopes of the linear parts of these plots yield the exponent $\eta$ (table 3); (b) and (e) semi-log (base 10) plots of $c(r)$ versus $r$ for different energies $E > E_{KT}$; (c) ($N_c^2 = 128^2$) and (f) ($N_c^2 = 64^2$) show plots of $\eta$ and $N_0$ versus $E$ (on the time scales of our runs $N_0$ is nonzero; see the text for a detailed discussion).</figcaption></figure>

exceeds $E_0$. As we show in appendix B, the $N_c$ dependence of the energy $E_{BKT}$, at which the BKT transition occurs, can be obtained approximately as follows. We begin with

$$ \delta\mathcal{E}_{BKT} = \delta\tilde{\mathcal{E}}_{BKT} \frac{8}{\log\left(\pi^2 N_c^2 \left(1 + \frac{\pi^2 N_c^2}{2g}\right)\right)}, \qquad (28) $$

where $\delta\tilde{\mathcal{E}}_{BKT}$, the estimate for the BKT transition energy that follows from an energy-entropy argument (see (20) in the appendix and [53]), is

$$ \delta\tilde{\mathcal{E}}_{BKT} = \frac{\pi^2 N_c^2}{2g} = \frac{\xi^2 k_{\max}^2}{2}, \qquad (29) $$

whence we obtain

$$ \delta\mathcal{E}_{BKT} = \frac{4 k_{\max}^2 \xi^2}{\log\left(k_{\max}^2 A\left(1 + \frac{k_{\max}^2 \xi^2}{2}\right)\right)}. \qquad (30) $$

We can now write

$$ E_{BKT} = E_0 \left( 1 + \frac{4\pi^2 N_c^2}{g \log(\pi^2 N_c^2 \left(1 + \frac{\pi^2 N_c^2}{2g}\right))} \right), \qquad (31) $$

by using this expression, we can determine the ratio $E_{BKT}(N_c^a)/E_{BKT}(N_c^b)$ for runs with two different values, $N_c^a$ and $N_c^b$, for the number of collocation points; we can also obtain this ratio from our DNS, by determining the value of $E$ at which the exponent $\eta$ becomes 1/4.
<table><caption>**Table 4.** The values of $E_0$, $\delta\tilde{\mathcal{E}}_{BKT}$ (see (29)), $\delta\mathcal{E}_{BKT}$ (see (30)), $E_{BKT}^A$ (see (31)) and $E_{BKT}^{\text{DNS}}$ from our DNS runs D1–D13 ($N_c = 64$) and E1–E12 ($N_c = 64$). $E_0$ is the ground state energy of a pure condensate of a uniform, interacting, 2D Bose gas and $E_{BKT}^{\text{DNS}}$ is BKT-transition energy determined using our DNS runs.</caption><thead><tr><th>N<sub>c</sub></th><th>E<sub>0</sub></th><th>δε̃<sub>BKT</sub></th><th>δε<sub>BKT</sub></th><th>E<sup>A</sup><sub>BKT</sub></th><th>E<sup>DNS</sup><sub>BKT</sub></th></tr></thead><tbody><tr><td>64</td><td>0.488</td><td>20.21</td><td>11.84</td><td>6.27</td><td>5.84</td></tr><tr><td>128</td><td>0.488</td><td>80.85</td><td>39.44</td><td>19.75</td><td>18.34</td></tr></tbody></table>

In table 4, we compare $E_{BKT}(N_c)$ for $N_c = 64$ and 128; our analytical approximation (31) yields $E_{BKT}^{128}/E_{BKT}^{64} \simeq 3.15$; this is in excellent agreement with the value $\simeq 3.14$ that we obtain for this ratio from our DNS results.

The thermalized state in the run A1 is in the BKT phase, because its total energy $E < E_{BKT}^{1024} = 818.7$ (31); thus, the system should be devoid of any free vortices, so the power-law tails in the velocity PDFs should disappear; this is indeed what we find. By contrast, for the thermalized states with energy $E > E_{BKT}$, free vortices and antivortices are present; so the velocity PDFs should show power-law tails, in such states, as we show explicitly in figure F3 in the supplementary material (available from stacks.iop.org/NJP/15/113025/mmedia) for the run D10, which has $E > E_{BKT}$ in the thermalized state.

## 4. Conclusions

We have carried out an extensive study of the statistical properties of the dissipationless, unforced, 2D, Fourier-truncated, GP equation. Our study has been designed specifically to study and identify the universal features, if any, of the turbulent evolution of the solutions of this equation, by undertaking a systematic DNS. In our study, we have used statistical measures such as velocity-component PDFs and energy and occupation-number spectra, for a large number of initial conditions. To the best of our knowledge, such a comprehensive study of the Fourier-truncated, 2D GP equation has not been attempted hitherto.

Our comprehensive study of the Fourier-truncated, 2D GP equation, which makes use of the three types of initial conditions (section 2.2) and a wide range of parameters (tables 1 and 3), allows us to systematize the dynamical evolution of this system into four different regimes, with qualitatively different statistical properties. This demarcation of the evolution into different regimes has not been systematized in earlier studies, which have concentrated only on one or two of these regimes. For example, the study of White et al [39] has investigated states with a significant number of vortex–anitvortex pairs and obtained for them PDFs of velocity components that have power-law tails of the type shown in figure 2. Small et al [47], Foster et al [55] and Damle et al [72] have investigated the BKT nature of the thermalized state. Wave-turbulence studies [32, 41, 73] have focused on power-law regions in energy and occupation-number spectra of the type we find in our third regime. The DNS studies in [41–45, 48, 74] have considered the time evolution of spectra and PDFs for the Fourier-truncated, 2D GP equation; in some cases, these studies introduce dissipation or hyperviscosity and forcing; they have also reported different power laws in spectra [42, 43, 45]. A recent theoretical and numerical study [48] has studied power-law regimes, associated with metastable transient states in the 2D
GP system; it suggests that these power-law regimes occur because of non-thermal fixed points. Our work suggests that, for the class of initial conditions which we have considered, at least in the dissipationless, unforced, Fourier-truncated, 2D GP equation, the only robust power laws in spectra are the ones we have reported above; all other apparent power laws occur either (a) for very special initial conditions [44] or (b) last for fleetingly small intervals of time and extend over very small ranges of k.

To recapitulate, we find that, in the first dynamical evolution regime of the Fourier-truncated, 2D GP equation, there are initial-condition-dependent transients. In the second regime the energy and the occupation-number spectra start to develop power-law scaling regions, but the power-law exponent and the extent of the scaling region change with time and are influenced by the initial conditions. In the third regime, of partial thermalization, we find $E_{\text{kin}}^c(k)$ and $E_{\text{int}}(k) + E_q(k) \sim k$, and $n(k) \sim 1/k$, for $k < k_c(t)$ and, for $k > k_c$, we find an initial-condition-dependent self-truncation regime, in which the spectra drop rapidly; the self-truncation wave number $k_c(t)$ grows either as $t^\alpha$ or logarithmically for different initial conditions (table 2). In the fourth, complete-thermalization regime, power-law forms of correlation functions and spectra, for $E < E_{\text{BKT}}$, are consistent with their nontrivial BKT forms; however, considerable care must be exercised, as explained in section 3.4.1 and [47, 55, 72], to distinguish these nontrivial power laws from their wave-turbulence analogues [32, 41, 73].

We have calculated a variety of order-p structure functions in our study; a full analysis of these structure functions lies outside the scope of this paper. However, we comment briefly on the use of the extended-self-similarity (ESS) procedure [75–77] here. Recall that, in the ESS procedure, log-log plots of the order-p structure functions $S_p(r)$ versus the third-order structure function $S_3(r)$ extend the scaling region, and their slopes yield estimates for the multiscaling exponent ratios $\zeta_p/\zeta_3$ that are better than the estimates of $\zeta_p$ obtained from slopes of log-log plots of $S_p(r)$ versus r; this procedure is especially valuable if $\zeta_3 = 1$, as it is in 3D fluid turbulence, by virtue of the von Kármán–Howarth relation. We have calculated structure functions [78] of the $\psi$ [79], velocity and vorticity fields in the 2D, Fourier-truncated, GP equation. The ESS procedure works here insofar as it extends the range over which scaling occurs; however, the exponents that follow from such plots evolve in time in a manner that mirrors the evolution of the spectral exponents that we have described above. A full elucidation of ESS and multiscaling in the 2D, Fourier-truncated, GP system will be presented elsewhere [80].

## Acknowledgments

We thank CSIR, DST and UGC (India) for financial support, and SERC (IISc) for computational resources. We acknowledge support from the Indo-French Centre for Applied Mathematics (UMI, IFCAM).

## Appendix A

The GP equation, which describes the dynamical evolution of the wave function $\psi(\mathbf{x}, t)$ of a weakly interacting, 2D Bose gas at low temperatures, is

$$i\hbar \frac{\partial \psi(\mathbf{x}, t)}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi(\mathbf{x}, t) + g_{2D} |\psi|^2 \psi(\mathbf{x}, t), \quad (A.1)$$
where $g_{2D}$ is the effective interaction strength. As we have mentioned earlier (see (2) and (3)), the GP equation conserves the energy, given by the Hamiltonian

$$H = \int_{\mathcal{A}} d^2x \left( \frac{\hbar^2}{2m} |\nabla \psi|^2 + \frac{g_{2d}}{2} |\psi|^4 \right), \qquad (A.2)$$

and the total number of particles $N = \int_{\mathcal{A}} |\psi|^2 d^2x$. To obtain (1), we first divide (A.1) by $\hbar$ and define $g = g_{2D}/\hbar$; we then set $\hbar/2m = 1$, with $m = 1$, so that $|\psi|^2$ is the same as $\rho$; this is tantamount to using units with $\hbar = 2$. The energy of the system, as expressed in equation (2), is obtained by dividing the equation (A.2) by $\hbar = 2$. A comparison with the experimental values can be made by noting that the healing energy $E_h = \hbar^2/(2m\xi^2)$; in our units this is simply $E_h = 1/\xi^2$; $\xi$ depends on $g$ through equation (16); and $g$ is related to $g_{2D}$ as mentioned above. The interaction strength $g_{2D}$ depends, inter alia, on the scattering length of the 2D interaction potential and the size of the BEC-cloud; for more details about the actual form of $g_{2D}$ in 2D see [57]. The wave number $k = n_i \Delta k$, where $\Delta k = 2\pi/L$ and $n_i = 0, 1, 2, \dots$; $L$ is the length of the side of a square simulation domain.

## Appendix B

The BKT transition is best studied by using the renormalization group [53]; here, we restrict ourselves to the heuristic, energy-entropy argument to obtain a rough estimate of the BKT transition temperature $T_{BKT}$. In the XY model, this transition is studied by using the Hamiltonian

$$H_{XY} = -J \sum_{\langle i, j \rangle} \cos(\theta_i - \theta_j), \qquad (B.1)$$

where $\langle i, j \rangle$ denotes nearest-neighbour pairs of sites, on a 2D square lattice, $J$ is the nearest-neighbour exchange coupling and $(\theta_i - \theta_j)$ is the angle between the nearest-neighbour, XY spins on sites $i$ and $j$. In the continuum limit, the above Hamiltonian becomes, to lowest order in spatial gradients,

$$H_{XY} = \frac{J}{2} \int d^2x (\nabla\theta(x))^2. \qquad (B.2)$$

By comparing (B.2) with the kinetic-energy term in (A.1), we find that

$$J = \frac{|\langle\psi\rangle|^2 \hbar^2}{m} = \frac{\rho \Gamma^2}{(2\pi)^2}, \qquad (B.3)$$

where $\Gamma$ denotes the Onsager–Feynman quantum of velocity circulation $\Gamma = 4\pi\hbar/2m = h/m$. A rough estimate for the BKT transition temperature $T_{BKT}$ is given below:

$$\tilde{T}_{BKT} = \frac{\pi J}{2k_B} = \frac{\pi |\langle\psi\rangle|^2 \hbar^2}{2mk_B} = \frac{\rho \Gamma^2}{8\pi k_B}, \qquad (B.4)$$

where $\tilde{T}_{BKT}$ denotes the estimate for $T_{BKT}$ that follows from an energy-entropy argument [53]. For $T < T_{BKT}$, the phase correlation function $c(r)$ (see (17)) and the angle-integrated spectrum $\hat{c}(k)$, which follows from a Fourier transform of $c(r)$, scale as

$$c(r) \sim (a/r)^{\frac{T}{4\pi k_B}} \qquad (B.5)$$
and

$$ \hat{c}(k) \sim k^{-1+\frac{T}{4T_{\text{BKT}}}}, \qquad (B.6) $$

respectively. Above $T_{\text{BKT}}$ the correlation length

$$ \ell = \frac{\int k^{-1} E(k) \, dk}{\int E(k) \, dk} \qquad (B.7) $$

is finite; and, as $T \to T_{\text{BKT}}$, it displays the essential singularity

$$ \ell \sim \exp(b(T_{\text{BKT}}/(T - T_{\text{BKT}}))^{1/2}). \qquad (B.8) $$

We now develop an analytical framework, which is valid at low-temperatures $T \ll T_{\text{BKT}}$, that can be used to test some of the results of our DNS runs in the region of complete thermalization. We first calculate equilibrium thermodynamic functions for a weakly interacting, 2D Bose gas, in the grand-canonical ensemble; we then obtain their analogues in the microcanonical ensemble. In the grand-canonical ensemble the probability of a given state is

$$ \mathbb{P} = \frac{1}{\Xi} e^{-\beta(H - \mu N)}, \qquad (B.9) $$

where $\Xi$ is the grand partition function, $\beta$ the inverse temperature, $\mu$ the chemical potential and $N$ the number of bosons. The grand-canonical potential is

$$ \Omega = -\beta^{-1} \log(\Xi), \qquad (B.10) $$

and the mean energy $E$, entropy $S$ and $N$ are

$$ N = -\frac{\partial \Omega}{\partial \mu}, \qquad (B.11a) $$

$$ S = \beta^2 \frac{\partial \Omega}{\partial \beta}, \qquad (B.11b) $$

$$ E = \frac{\partial \Omega}{\partial \beta} + \mu N = \frac{S}{\beta} + \mu N. \qquad (B.11c) $$

We adapt to 2D the 3D study of [52], expand $\psi$ in terms of Fourier modes $A_k$ and obtain $\Omega$ as the sum of the saddle-point part $\Omega_{sp}$ and $\Omega_Q$, the deviations from the saddle point that are quadratic in $A_k$. We write $\Omega = \Omega_{sp} + \Omega_Q$, where $\Omega_{sp} = -A\mu^2/2g$ and

$$ \Omega_Q = - \int_0^{p_{\text{max}}} dp \frac{\left( p A \log \left( \frac{2m}{\beta \sqrt{p^4 + 4mp^2\mu}} \right) \right)}{2\pi \beta \hbar^2}. \qquad (B.12) $$

We can also calculate the condensate depletion $\delta N$, where the particle number $N = N_0 + \delta N$ and $N_0$ is the number of particles in the $k=0$ mode, as follows:

$$ \delta N = \int_0^{p_{\text{max}}} dp \frac{mpA \left( p^{-2} + \frac{1}{p^2 + 4m\mu} \right)}{2\pi \beta \hbar^2}. \qquad (B.13) $$

The integrals in (B.12) and (B.13) can be performed analytically, but, in contrast to the 3D case where the primitives are zero at $p=0$, the 2D primitive for $\Omega_{ph}$ is finite at $p=0$ and for $\delta N$ it is infra-red (IR) divergent. By subtracting the IR finite and divergent terms from $\Omega_Q$ and $\delta N$,
respectively, we obtain the following expressions, in 2D, in the thermodynamic limit $A \to \infty$:

$$ \Omega = - \frac{\mu^2 A}{2g} - \frac{p_{\max}^2 A}{4\pi \beta \hbar^2} + \frac{m\mu A \log \left(1 + \frac{p_{\max}^2}{4m\mu}\right)}{2\pi \beta \hbar^2} - \frac{p_{\max}^2 A \log \left(\frac{2m}{\beta\sqrt{p_{\max}^4 + 4m\mu p_{\max}^2}}\right)}{4\pi \beta \hbar^2} \quad (B.14) $$

and

$$ \delta N = \frac{m A \left( \log \left( 1 + \frac{p_{\max}^2}{4m\mu} \right) + \log \left( \frac{p_{\max}^2 A}{\hbar^2} \right) \right)}{4\pi \beta \hbar^2} \quad (B.15) $$

By using the thermodynamic relations (B.11), we obtain

$$ N = \frac{\mu A}{g} - \frac{m A \log \left( 1 + \frac{p_{\max}^2}{4m\mu} \right)}{2\pi \beta \hbar^2} \quad (B.16) $$

and

$$ E = \frac{\mu^2 A}{2g} + \frac{p_{\max}^2 A}{4\pi \beta \hbar^2} - \frac{m \mu A \log \left( 1 + \frac{p_{\max}^2}{4m\mu} \right)}{2\pi \beta \hbar^2} \quad (B.17) $$

We next determine the chemical potential $\mu$, which fixes the total density $\rho = mN/A$ at a given value, by solving the equation

$$ \rho - \frac{m\mu}{g} + \frac{m^2 \log \left( 1 + \frac{p_{\max}^2}{4m\mu} \right)}{2\pi \beta \hbar^2} = 0 \quad (B.18) $$

at $\beta = \infty$, i.e. zero temperature (subscript 0), we obtain

$$ \mu_0 = \frac{g \rho}{m} \quad (B.19) $$

to order $\beta^{-1}$ we obtain

$$ \mu = \mu_0 + \delta\mu, \quad (B.20) $$

where

$$ \delta\mu = \frac{mg(4g\rho^2 + \rho p_{\max}^2) \log \left(1 + \frac{p_{\max}^2}{4g\rho}\right)}{m^2 p_{\max}^2 + 2\pi\beta\hbar^2 \rho p_{\max}^2 + 8\pi\beta\hbar^2 g\rho^2} \quad (B.21) $$

We insert $\mu$ from (B.20) into (B.15), define the change in density $\delta\rho = m\delta N/A$, use the energy $E$ from (B.17) and then expand to order $\beta^{-1}$ to obtain

$$ \delta\rho = \frac{m^2 \left( \log \left( 1 + \frac{p_{\max}^2}{4g\rho} \right) + \log \left( \frac{p_{\max}^2 A}{\hbar^2} \right) \right)}{4\pi\beta\hbar^2} \quad (B.22) $$

and

$$ E = \frac{g\rho^2 A}{2m^2} + \frac{p_{\max}^2 A}{4\pi\beta\hbar^2} \quad (B.23) $$
By using (B.4) and $\rho = m|\langle\psi\rangle|^2$, we obtain

$$ \tilde{\beta}_{\text{BKT}} = \frac{1}{k_{\text{B}} \tilde{T}_{\text{BKT}}} = \frac{2m^2}{\pi \rho \hbar^2}, \qquad (B.24) $$

which we can use along with (B.22) to relate the condensate relative depletion $\delta\rho/\rho$ to $\beta/\tilde{\beta}_{\text{BKT}}$, where $\beta = 1/(k_{\text{B}}T)$ and $k_{\text{B}}$ is the Boltzmann constant, as given below:

$$ \frac{\delta\rho}{\rho} = \frac{\tilde{\beta}_{\text{BKT}}}{8\beta} \log \left( \frac{p_{\text{max}}^2 \left( 1 + \frac{p_{\text{max}}^2}{4g\rho} \right) A}{\hbar^2} \right). \qquad (B.25) $$

We use this low-temperature result (B.25) to estimate the inverse-temperature scale $\beta_{\text{BKT}}$, at which the depletion of the $k=0$ condensate mode becomes significant for a finite-size system with $N_c^2$ collocation points (which fixes the maximum momentum $p_{\text{max}}$); in particular, we can solve (B.25), for $\delta\rho/\rho = 1$, to obtain

$$ \frac{\beta_{\text{BKT}}}{\tilde{\beta}_{\text{BKT}}} = \frac{1}{8} \log \left( \frac{p_{\text{max}}^2 \left( 1 + \frac{p_{\text{max}}^2}{4g\rho} \right) A}{\hbar^2} \right). \qquad (B.26) $$

By making the replacements that correspond to defining $\hbar$, $m$ and $g$ in terms of $c$ and $\xi$, as in [52], $p_{\text{max}} \rightarrow \hbar k_{\text{max}}$, $\hbar \rightarrow \sqrt{2cm\xi}$ and $g \rightarrow c^2m^2/\rho$, we can rewrite (B.26) as

$$ \frac{\beta_{\text{BKT}}}{\tilde{\beta}_{\text{BKT}}} = \frac{1}{8} \log \left( k_{\text{max}}^2 A \left( 1 + \frac{k_{\text{max}}^2 \xi^2}{2} \right) \right). \qquad (B.27) $$

Our DNS runs, which use initial conditions of types IC1 and IC2, give the dynamical evolutions of the Fourier-truncated, 2D GP equation, which is a Hamiltonian system. The energy $E$, particle number $N$ and area $A$ are conserved in this evolution, so our calculation can be viewed as a simulation of this Hamiltonian system in the microcanonical ensemble, which yields, eventually, the fully thermalized state that we have described above. Therefore, we now transform the results, which we have obtained in the previous subsection, into their counterparts in the microcanonical ensemble. In the low-temperature limit, (B.23) yields

$$ \beta = \frac{m^2 p_{\text{max}}^2 A}{2\pi \hbar^2 (2m^2 E - g\rho^2 A)}. \qquad (B.28) $$

The energy of a pure condensate is

$$ E_0 = \lim_{\beta \to \infty} E = \frac{g \rho^2 A}{2 m^2} \qquad (B.29) $$

and the energy and the inverse temperature $\beta$ (B.28) can be related as follows:

$$ E = E_0(1 + \delta\mathcal{E}), \qquad (B.30) $$

where $\delta\mathcal{E}$ is the relative increase of energy above $E_0$, and

$$ \beta = \frac{m^2 p_{\text{max}}^2}{2\pi \hbar^2 g \rho^2 \delta\mathcal{E}}. \qquad (B.31) $$

If we now substitute $\beta = \beta_{\text{BKT}}$ by using (B.26), we obtain, in terms of $c$, $\xi$ and $\rho$ (see text just below (B.26))

$$ E_0 = \frac{c^2 \rho A}{2}, \qquad (B.32) $$
$$ \delta\tilde{\mathcal{E}}_{\text{BKT}} = \frac{k_{\text{max}}^2 \xi^2}{2} \qquad (\text{B.33}) $$

and

$$ \delta\mathcal{E}_{\text{BKT}} = \frac{4k_{\text{max}}^2 \xi^2}{\log\left(k_{\text{max}}^2 A \left(1 + \frac{k_{\text{max}}^2 \xi^2}{2}\right)\right)} \qquad (\text{B.34}) $$

All the energies mentioned in the main paper are dimensionless; thus, to convert the energies given in this appendix to dimensionless forms, we divide them by $\hbar$. Hence, the energy of a pure condensate is obtained, in the dimensionless form, by dividing (B.29) by $\hbar$, which gives

$$ E_0 = \frac{g}{2A} = \frac{1}{2} \frac{g}{L^2}. \qquad (\text{B.35}) $$

## References

[1] Feynman R P 1955 Application of quantum mechanics to liquid helium *Progress in Low Temperature Physics* vol I ed C J Gorter (Amsterdam: North-Holland) p 17

[2] Vinen W F 1957 Mutual friction in a heat current in liquid helium II: I. Experiments on steady heat currents *Proc. R. Soc. Lond. A* **240** 114–27

[3] Vinen W F 1957 Mutual friction in a heat current in liquid helium II: II. Experiments on transient effects *Proc. R. Soc. Lond. A* **240** 128–43

[4] Vinen W F 1957 Mutual friction in a heat current in liquid helium II: III. Theory of the mutual friction *Proc. R. Soc. Lond. A* **242** 493–515

[5] Tilley D R and Tilley J 1990 *Superfluidity and Superconductivity* 3rd edn (Bristol: IOP Publishing)

[6] Donnelly R J 1991 *Quantized Vortices in Helium II* vol 3 (Cambridge: Cambridge University Press)

[7] Barenghi C F, Donnelly R J and Vinen W F (ed) 2001 *Quantized Vortex Dynamics and Superfluid Turbulence* (*Lecture Notes in Physics* vol 571) (Berlin: Springer)

[8] Vinen W F and Niemela J J 2002 Quantum turbulence *J. Low Temp. Phys.* **128** 167–231

[9] Niemela J J 2008 Reconnecting to superfluid turbulence *Physics* **1** 26

[10] Procaccia I and Sreenivasan K R 2008 The state of the art in hydrodynamic turbulence: past successes and future challenges *Physica D: Nonlinear Phenom.* **237** 2167–83 (Euler equations: 250 years on—*Proc. Int. Conf.*)

[11] Tsubota M 2008 Quantum turbulence *J. Phys. Soc. Japan* **77** 111006

[12] Paoletti M S and Lathrop D P 2011 Quantum turbulence *Annu. Rev. Condens. Matter Phys.* **2** 213–34

[13] Skrbek L and Sreenivasan K R 2012 Developed quantum turbulence and its decay *Phys. Fluids* **24** 011301

[14] Maurer J and Tabeling P 1998 Local investigation of superfluid turbulence *Europhys. Lett.* **43** 29

[15] Smith M R, Donnelly R J, Goldenfeld N and Vinen W F 1993 Decay of vorticity in homogeneous turbulence *Phys. Rev. Lett.* **71** 2583–6

[16] Stalp S R, Skrbek L and Donnelly R J 1999 Decay of grid turbulence in a finite channel *Phys. Rev. Lett.* **82** 4831–4

[17] Skrbek L, Niemela J J and Donnelly R J 2000 Four regimes of decaying grid turbulence in a finite channel *Phys. Rev. Lett.* **85** 2973–6

[18] Skrbek L and Stalp S R 2000 On the decay of homogeneous isotropic turbulence *Phys. Fluids* **12** 1997–2019

[19] Stalp S R, Niemela J J, Vinen W F and Donnelly R J 2002 Dissipation of grid turbulence in helium II *Phys. Fluids* **14** 1377–9

[20] Fisher S N, Hale A J, Guénault A M and Pickett G R 2001 Generation and detection of quantum turbulence in superfluid $^3$He – B *Phys. Rev. Lett.* **86** 244–7

[21] Henn E A L, Seman J A, Roati G, Magalhães K M F and Bagnato V S 2009 Emergence of turbulence in an oscillating Bose–Einstein condensate *Phys. Rev. Lett.* **103** 045301
[22] Neely T W, Bradley A S, Samson E C, Rooney S J, Wright E M, Law K J H, Carretero-González R, Kevrekidis P G, Davis M J and Anderson B P 2012 Characteristics of two-dimensional quantum turbulence in a compressible superfluid arXiv:1204.1102v2

[23] Sun C, Jia S, Barsi C, Rica S, Picozzi A and Fleischer J W 2012 Observation of the kinetic condensation of classical waves *Nature Phys.* **8** 470–4

[24] Roche P-E, Barenghi C F and Leveque E 2009 Quantum turbulence at finite temperature: the two-fluids cascade *Europhys. Lett.* **87** 54006

[25] Salort J, Roche P-E and Lévêque E 2011 Mesoscale equipartition of kinetic energy in quantum turbulence
*Europhys. Lett.* **94** 24001

[26] Schwarz K W 1985 Three-dimensional vortex dynamics in superfluid $^4$He: line–line and line–boundary interactions *Phys. Rev. B* **31** 5782–804

[27] Schwarz K W 1988 Three-dimensional vortex dynamics in superfluid $^4$He: homogeneous superfluid turbulence *Phys. Rev. B* **38** 2398–417

[28] Barenghi C F and Samuels D C 1999 Self-consistent decay of superfluid turbulence *Phys. Rev. B* **60** 1252–60

[29] Adachi H, Fujiyama S and Tsubota M 2010 Steady-state counterflow quantum turbulence: simulation of vortex filaments using the full Biot–Savart law *Phys. Rev. B* **81** 104511

[30] Kagan Yu M, Svistunov B V and Shlyapnikov G V 1992 Kinetics of Bose condensation in an interacting Bose gas *Sov. Phys.—JETP* **75** 387

[31] Zakharov V E and Nazarenko S V 2005 Dynamics of the Bose–Einstein condensation *Physica D* **201** 203–11

[32] Svistunov B V 1991 Highly nonequilibrium Bose condensation in a weakly interacting gas *J. Moscow Phys. Soc.* **1** 373

[33] Lvov Yu, Nazarenko S V and West R J 2003 Wave turbulence in Bose–Einstein condensates *Physica D* **184** 333–51

[34] Koplik J and Levine H 1993 Vortex reconnection in superfluid helium *Phys. Rev. Lett.* **71** 1375–8

[35] Nore C, Abid M and Brachet M E 1997 Kolmogorov turbulence in low-temperature superflows *Phys. Rev. Lett.* **78** 3896–9

[36] Berloff N G and Svistunov B V 2002 Scenario of strongly nonequilibrated Bose–Einstein condensation *Phys. Rev. A* **66** 013603

[37] Kobayashi M and Tsubota M 2005 Kolmogorov spectrum of superfluid turbulence: numerical analysis of the Gross–Pitaevskii equation with a small-scale dissipation *Phys. Rev. Lett.* **94** 065302

[38] Krstulovic G and Brachet M 2011 Dispersive bottleneck delaying thermalization of turbulent Bose–Einstein condensates *Phys. Rev. Lett.* **106** 115303

[39] White A C, Barenghi C F, Proukakis N P, Youd A J and Wacks D H 2010 Nonclassical velocity statistics in a turbulent atomic Bose–Einstein condensate *Phys. Rev. Lett.* **104** 075301

[40] Connaughton C, Josserand C, Picozzi A, Pomeau Y and Rica S 2005 Condensation of classical nonlinear waves *Phys. Rev. Lett.* **95** 263901

[41] Nazarenko S and Onorato M 2007 Freely decaying turbulence and Bose–Einstein condensation in Gross–Pitaevskii model *J. Low Temp Phys.* **146** 31–46

[42] Numasato R, Tsubota M and L'vov V S 2010 Direct energy cascade in two-dimensional compressible quantum turbulence *Phys. Rev. A* **81** 063630

[43] Nowak B, Sexty D and Gasenzer T 2011 Superfluid turbulence: nonthermal fixed point in an ultracold Bose gas *Phys. Rev. B* **84** 020506

[44] Bradley A S and Anderson B P 2012 Energy spectra of vortex distributions in two-dimensional quantum turbulence *Phys. Rev. X* **2** 041001

[45] Nowak B, Schole J, Sexty D and Gasenzer T 2012 Nonthermal fixed points, vortex statistics and superfluid turbulence in an ultracold Bose gas *Phys. Rev. A* **85** 043627

[46] Reeves M T, Anderson B P and Bradley A S 2012 Classical and quantum regimes of two-dimensional turbulence in trapped Bose–Einstein condensates *Phys. Rev. A* **86** 053621

[47] Small E, Pugatch R and Silberberg Y 2011 Berezinskii–Kosterlitz–Thouless crossover in a photonic lattice
*Phys. Rev. A* **83** 013806
[48] Schole J, Nowak B and Gasenzer T 2012 Critical dynamics of a two-dimensional superfluid near a nonthermal fixed point *Phys. Rev. A* **86** 013624

[49] Cichowlas C, Bonaïti P, Debbasch F and Brachet M 2005 Effective dissipation and turbulence in spectrally truncated Euler flows *Phys. Rev. Lett.* **95** 264502

[50] Frisch U, Kurien S, Pandit R, Pauls W, Ray S S, Wirth A and Zhu J-Z 2008 Hyperviscosity, Galerkin truncation and bottlenecks in turbulence *Phys. Rev. Lett.* **101** 144501

[51] Ray S S, Frisch U, Nazarenko S and Matsumoto T 2011 Resonance phenomenon for the Galerkin-truncated Burgers and Euler equations *Phys. Rev. E* **84** 016301

[52] Krstulovic G and Brachet M 2011 Energy cascade with small-scale thermalization, counterflow metastability and anomalous velocity of vortex rings in Fourier-truncated Gross-Pitaevskii equation *Phys. Rev. E* **83** 066311

[53] Kogut J B 1979 An introduction to lattice gauge theory and spin systems *Rev. Mod. Phys.* **51** 659-713

[54] Chaikin P M and Lubensky T C 1995 *Principles of Condensed Matter Physics* (Cambridge: Cambridge University Press)

[55] Foster C J, Blakie P B and Davis M J 2010 Vortex pairing in two-dimensional Bose gases *Phys. Rev. A* **81** 023623

[56] Pethick C J and Smith H 2001 *Bose-Einstein Condensation in Dilute Gases* (Cambridge: Cambridge University Press)

[57] Posazhennikova A 2006 Colloquium: weakly interacting, dilute Bose gases in 2D *Rev. Mod. Phys.* **78** 1111-34

[58] Mermin N D and Wagner H 1966 Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models *Phys. Rev. Lett.* **17** 1133-6

[59] Hohenberg P C 1967 Existence of long-range order in one and two dimensions *Phys. Rev.* **158** 383-6

[60] Berezinskii V L 1971 Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group: I. Classical systems *Sov. Phys.—JETP* **32** 493

[61] Kosterlitz J M and Thouless D J 1973 Ordering, metastability and phase transitions in two-dimensional systems *J. Phys. C: Solid State Phys.* **6** 1181

[62] Minguzzi A, Succi S, Toschi F, Tosi M P and Vignolo P 2004 Numerical methods for atomic quantum gases with applications to Bose-Einstein condensates and to ultracold fermions *Phys. Rep.* **395** 223-355

[63] Frisch U, Ray S S, Sahoo G, Banerjee D and Pandit R 2013 Real-space manifestations of bottlenecks in turbulence spectra *Phys. Rev. Lett.* **110** 064501

[64] Kubo R 1966 The fluctuation-dissipation theorem *Rep. Prog. Phys.* **29** 255

[65] Miguel M and Toral R 2000 Stochastic effects in physical systems *Instabilities and Nonequilibrium Structures VI* (Nonlinear Phenomena and Complex Systems vol 5) ed E Tirapegui, J Martínez and R Tiemann (Dordrecht: Springer) pp 35-127

[66] Abid M, Huepe C, Metens S, Nore C, Pham C T, Tuckerman L S and Brachet M E 2003 Gross-Pitaevskii dynamics of Bose-Einstein condensates and superfluid turbulence *Fluid Dyn. Res.* **33** 509

[67] Paoletti M S, Fisher M E, Sreenivasan K R and Lathrop D P 2008 Velocity statistics distinguish quantum turbulence from classical turbulence *Phys. Rev. Lett.* **101** 154501

[68] Adachi H and Tsubota M 2011 Numerical study of velocity statistics in steady counterflow quantum turbulence *Phys. Rev. B* **83** 132503

[69] Davide P, Nazarenko S and Onorato M 2012 Sustained turbulence in the three-dimensional Gross-Pitaevskii model *Physica D* **241** 304-14 (special issue on small scale turbulence)

[70] Baggaley A W and Barenghi C F 2011 Quantum turbulent velocity statistics and quasi-classical limit *Phys. Rev. E* **84** 067301

[71] Düring G, Picozzi A and Rica S 2009 Breakdown of weak-turbulence and nonlinear wave condensation *Physica D* **238** 1524-49

[72] Damle K, Majumdar S N and Sachdev S 1996 Phase ordering kinetics of the Bose gas *Phys. Rev. A* **54** 5037-41

[73] Nazarenko S and Onorato M 2006 Wave turbulence and vortices in Bose-Einstein condensation *Physica D* **219** 1-12
[74] Matthew T R, Billam T P, Anderson B P and Bradley A S 2013 Inverse energy cascade in forced 2D quantum turbulence *Phys. Rev. Lett.* **110** 104501

[75] Briscolini M, Santangelo P, Succi S and Benzi R 1994 Extended self-similarity in the numerical simulation of three-dimensional homogeneous flows *Phys. Rev. E* **50** R1745–7

[76] Pandit R, Perlekar P and Ray S S 2009 Statistical properties of turbulence: an overview *Pramana* **73** 157–91

[77] Chakraborty S, Frisch U and Ray S S 2010 Extended self-similarity works for the Burgers equation and why *J. Fluid Mech.* **649** 275

[78] Frisch U 1996 *Turbulence* (Cambridge: Cambridge University Press)

[79] Dyachenko A and Falkovich G 1996 Condensate turbulence in two dimensions *Phys. Rev. E* **54** 5095–9

[80] Shukla V, Brachet M and Pandit R 2013 in preparation