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

We undertake a systematic, direct numerical simulation (DNS) 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. First, 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 depended 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 forms.

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: (1) The first is the region of initial transients; this depends on the initial conditions. (2) 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 behaviours, but the powerlaw exponent and the extents of the scaling regions change with time and depend on the initial conditions. (3) In the third regime, which we call the region of partial thermalization, these spectra show clear, power-law, scaling behaviours, with a power that is independent of the initial conditions, and, at large wave vectors, an initialcondition-dependent, self-truncation regime, where spectra drop rapidly; (4) 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,[52][53][54], 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,51].
The remaining part of this paper is organised 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 subsections 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 energy and the occupation-number spectra for different initial conditions); (c) a similar statistical characterization, as in subsection 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 the details of some analytical calculations are presented in Appendix A and Appendix B, respectively.

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.

The Gross-Pitaevskii Equation
The GP equation, which describes the dynamical evolution of the wave function ψ of a weakly interacting 2D Bose gas at low temperatures, is ψ(x, t) is a complex, classical field and g is the effective interaction strength [55,56]. This equation conserves the total energy and the total number of particles where A = L 2 is the area of our 2D, periodic, computational domain of side L. From (1) we obtain the continuity equation where ρ = |ψ| 2 is interpreted as the particle density and the velocity is We can use the Madelung transformation ψ( where the kinetic, interaction, and quantum-pressure energies are defined, respectively, as We separate the compressible (supercript c) and the incompressible (superscript i) parts of the kinetic energy by making use of the decomposition where ∇ · (ρ 1/2 v) i = 0 and ∇ × (ρ 1/2 v) c = 0, whence we obtain the following: The spectra for these energies are defined as follows: and furthermore, we define an occupation-number spectrum n(k) via here we denote the Fourier transform of A(x) by 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. [48,49] for studies of the Fourier-or Galerkin-truncated Euler equation). The Bogoluibov dispersion relation ω B (k) is obtained by linearizing (1) around a constant ψ. For a total number of particles (3) N = 1, it is where the sound velocity is c = √ 2g L and the coherence length is 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, Berezenskii-Kosterlitz-Thouless (BKT) phase [57][58][59][60], which shows quasi-long-range order with an algebraic decay of the spatial correlation function [52] for temperatures T below the transition temperature T BKT (or energy E BKT in the microcanonical ensemble), where r ≡ |r| and the critical exponent η < 0.25 for T < T BKT ; and η = 0.25 at T = T BKT [53]. The BKT phase shows bound vortex-antivortex pairs; these unbind above T BKT , so c(r) ∼ e −r/ℓ , in the disordered phase, with ℓ the correlation length.

Numerical Methods and Initial Conditions
To perform a systematic, pseudospectral, direct numerical simulation (DNS) of the spatiotemporal evolution of the 2D, Fourier-truncated, GP equation, we have developed a parallel, MPI code in which we discretize ψ(x, t) on a square simulation domain of side L = 32 with N 2 c 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 ∆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 [61]. Thus, the maximum wave number k max = (N c /2)∆k, where ∆k = 2π/L, and 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.
To initiate turbulence in the 2D GP equation we use three types of initial conditions IC1 [41], IC2, and IC3 [51], always normalized to correspond to a total number of particles (3) N = 1. The first of these is best represented in Fourier space as follows: where k = k 2 x + k 2 y , Θ(k x , k y ) are random numbers distributed uniformly on the interval 0, 2π ; k 0 = N 0 ∆k and σ = B∆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 ψ at time t = 0. The initial condition IC2 is like IC1 but, in addition, it has a finite initial condensate population We obtain the initial condition IC3 by solving the 2D, stochastic, Ginzburg-Landau equation (SGLE), which follows from the free-energy functional where µ is the chemical potential §. The SGLE is where ζ is a zero-mean, Gaussian white noise with where D = 2T , in accordance with the fluctuation-dissipation theorem [62], T is the temperature, and δ the Dirac delta function. Finally, the SGLE (23) becomes which we solve along with the following, ad-hoc equation to control the number of particles N; the parameter N av controls the mean value of N; and ν 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 ∆t, for time marching and the method of reference [63] (see page 25 of this reference). § Recall that the SGLE can be thought of as an imaginary-time GP equation with external, additive noise (see, e.g. reference [38]) Table 1. Parameters for our DNS runs A1-A13, B1-B2, and C1-C6: N 2 c is the number of collocation points, k 0 is the energy-injection scale, σ is the Fourierspace width of ψ at t = 0; g is the effective interaction strength; N i 0 is the initial condensate population; D and k in c 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.

Results
We first present the time evolution of the different energies, the probability distribution functions (PDFs) of the velocity components, and the population N 0 in the zero-wavenumber 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).

Evolution of energies, velocity PDFs, and the zero-wave-number population
We show the early stages of the time evolution of the energies E i kin , E c kin , 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 i kin is a significant fraction of the total initial energy; the runs B1 and C6 start with initial configurations of type IC2 and IC3, respectively, in which E i kin is negligibly small at t = 0. The transient nature of the early stages of the dynamical evolution of the dissipationless, unforced, 2D, GP equation is evident from figure 1, in which we observe a rapid conversion of E i kin into the other three components, with a significant fraction being transferred to E c kin ; 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 i kin 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. Moreover, the larger the value of E i kin 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 i kin (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 intial conditions A1-A4. In contrast, the runs B1 and C6 have negligibly small values of E i kin (t = 0) to begin with (figures 1 (e) and (f), respectively); and E i kin (t) remains close to zero throughout the dynamical evolution here. For run B1, both E c kin and E q start from values close to zero, grow at the cost of E 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 crossover 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 fit to the form ∼ v −γ i , with γ ≃ 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 [64] and some numerical studies [39,45,65,66]. 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 fit 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 pseudocolor 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 i kin at t = 0, are given, respectively, in figures 2(d)-(f) and figures 2(g)-(i).
We turn now to the time evolution of the population N 0 (t), in the k = 0 mode [36,40,67], 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; these figures show that N 0 (t) increases with t, on average, and depends on E, g, k 0 , and σ. 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 σ are smaller for A1 than for A4; thus, N 0 (t) grows more slowly in A4 than in A1; and, after an equal amount of simulation time, its value in A4 is nearly an order of magnitude lower than in A1; the former shows large fluctuations in N 0 (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 0 (t = 0) = 0.95; in this case, after a period of initial transients, N 0 (t) → 0.98 over our simulation time. The run C6 (figure 3 (f)) uses an initial condition of type IC3; here N 0 (t) fluctuates slightly but remains close to its initial value (cf. [40,67]).
To study the dependence of N 0 (t) on the number of collocation points N 2 c , we evolve the initial configuration of A1 for N c = 512 (run A5), N c = 256 (run A6), N c = 128 (run A7), and N c = 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 (≃ 0.9) for the runs wth N c = 128 (run A7), N c = 256 (run A6), and N c = 1024 (run A1). In contrast, the saturation value for the run with N c = 64 (run A8) is ≃ 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 realisation of the random phases  Θ(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 [57,58]; however, it is not easy to realize this limit in finite-size systems and with the limited run times that are dictated by computational resources.
We discuss these issues again in section 3.4 and also refer the reader to [67,68].

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, IC2, and IC3. The first type begins with a sizeable incompressible kinetic energy spectrum E i kin (k); and the initial transients are associated with the annihilation and creation of vortexantivortex pairs, the associated depletion of E i kin (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 i kin (k), for the runs A1, A2, A3, and A4, to ascertain the presence of scaling behaviour, if any. We find that, in the low-k region, E i kin (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, 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 i kin (k) ∼ 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  (a)-(b)).
The initial transients described above are followed by a regime in which the energy and occupation-number spectra begin to show power-law-scaling behaviours, 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, where we illustrate the time evolution of E c kin . Figure 5 (a) shows E c kin (k) for the run A1; we begin to see a power-law region here with E c kin (k) ∼ k, on the low-k side of the peak after which the spectrum falls steeply. A similar E c kin (k) ∼ k behaviour starts to emerge in the region k k 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 int (k) + E q (k) ∼ k (figure 7) and n(k) ∼ 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 behaviours appear 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), and (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) and (f), and (II) spectra that increase all the way to k max , e.g., the spectra in figures 5 (h) and 6 (b) and (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 kin (k) = (E c kin (k) + E i kin (k)), referred to as the self-truncation wave number [51], moves out to k max as a power of t and (B) a subclass in which k c moves out to k max at a rate that is slower than a power of t. convergence to the GP PDE is lost in case (II); note that the scaling region with E c kin ∼ k sets in at high wave numbers close to k max and then extends to the low-wavenumber regime. For case (IA) analogous plots of E c kin (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.

Self-truncation
We now present a detailed characterization of the partialthermalization regime, when energy spectra display self-truncation at wave-numbers beyond k c (t), which can be defined as follows: as the system approaches complete thermalization, k c (t) → k max . In particular, we explore how the scaling ranges in energy spectra grow with t for different values of g, with the initial configuration and number of collocation points N c held fixed. For an initial condition of type IC1, with k 0 = 5∆k, σ = 2∆k, and N c = 256, we obtain the time evolution of energy spectra for g = 1000 (run A6), g = 2000 (run A9), and g = 5000 (run A10) in figures 9 (a), (b), and (c), respectively, and their video analogues (Videos S3 (panel V2) in the Supplementary Material). 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∆k, σ = 2∆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. 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∆k, σ = 2∆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 does the system lose spectral convergence.
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 intial conditions of type IC1.
With initial conditions of types IC1 and IC2, we cannot control the initial value k c (t = 0) ≡ k in c easily. However, initial conditions of type IC3, which we obtain from the SGLE, allow us to control k in c and start, therefore, with initial spectra that display partial thermalization for k < k in c [51] and a sharp fall thereafter. In figure 10 we show the time evolution of E c kin (k) for such initial conditions from runs C1-C6. For different representative values of k in c , 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 in c = 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 in c , starts with the smoothening of the sharp cut-off at k in c ; the higher the value of k in c , 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).
The growth of k c (t) with t, illustrated in figure 11 (a), can be fit to the form k c (t) ∼ t α ; however, as we show below, α depends on the initial condition. We obtain the exponent α 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 α 1 and α 2 , respectively. Note that in (ii) we have a parametric plot [38,51], shown in figure 11 (b); this yields a straight-line scaling regime with slope χ and α 2 = 1/ (1 − χ). The values of α 1 and α 2 , listed in table 2, show that α 1 ≃ α 2 ; the discrepancy between these two values for α is a convenient measure of the errors of our estimates. For runs C4, C5, and C6, we cannot obtain α 2 reliably; the small values of α 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.

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  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.
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 c kin for the runs B1 and A7, where E c kin shows power-law scaling over the entire wave number range, from k = 2π/L up to k max , towards the end of the respective simulations; a naïve fit is consistent with E c kin (k) ∼ k (but see below).  Table 2. Summary of the self-truncation results from our DNS runs A1-A4, B2, and C1-C6: E is the total energy; k max = 2πN c /2L; ξ = L/ √ g is the healing length; k i c and k f c are the initial and final values of k c (averaged over a few time steps); α 1 is the slope obtained from the log-log (base 10) plot of k c versus t and α 2 = 1/(1 − χ), where χ is the slope obtained from the log-log (base 10) plot of dk c /dt versus k c /k max .

Correlation functions and the BKT transition
should yield a BKT phase [52,54], with the correlation function c(r) ∼ r −η , at energies E < E BKT ; and c(r) should decay exponentially with r if E > E BKT . We show this explicitly now by using initial conditions of type IC1 with N c = 64 and N c = 128 and g = 1000; we obtain different energies by changing k 0 and σ (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 BKT , in figures 12 (a) and (d), where the straight lines indicate power-law regimes; and, for E > E BKT , we use semi-log 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 . Aside from this uncertainty, the behavior of c(r), in the regime of complete thermalization, is in accord with our expectations for the BKT phase; in particular, the exponent η (see equation (18)) depends on E for E < E BKT as shown in figures 12 (c) and (f). Our values for η, for the runs with E < E BKT and with N c = 64 and N c = 128, are listed in table 3. Note that E BKT ≃ 6(N c = 64) and E BKT ≃ 19(N c = 128), i.e., E BKT depends on N c , the number of collocation points; we show analytically below how a lowtemperature 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 [57,58] and n(k) ∼ k −1+η ; 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 η from the  Table 3. List of parameters for our complete-thermalization DNS runs D1-D13 (N 2 c = 128 2 ) and E1-E12 (N 2 c = 64 2 ): N 2 c is the number of collocation points; k 0 is the energy-injection scale; σ is Fourier-space width of ψ at t = 0; E is the total energy; and η is the exponent of the correlation function c(r) ∼ r −η 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. k = 0 part of n(k) as shown in reference [68]; similarly, E c kin (k) should also show a power-law form with an exponent that depends on η, but this is difficult to realize in numerical calculations with limited spatial resolutions and run lengths.

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 + δE); this energy E is fixed by the initial condition; and δE measures the relative amount by which E exceeds E 0 . As we show in the 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 where δẼ BKT , the estimate for the BKT transition energy that follows from an energyentropy argument (see (20) in the Appendix and [52]), is whence we obtain    Table 4. The values of E 0 , δẼ BKT (see (29)), δE (see (30)), E A BKT (see (31)), and E DNS BKT 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 DNS BKT is BKT-transition energy determined using our DNS runs.
We can now write by using this expression we can determine the ratio E BKT (N a c )/E BKT (N b c ) for runs with two different values, N a c and N b c , 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 η becomes 1/4. In Table 4 we compare E BKT (N c ) for N c = 64 and N c = 128; our analytical approximation (31) yields E 128 BKT /E 64 BKT ≃ 3.15; this is in excellent agreement with the value ≃ 3.14 that we obtain for this ratio from our DNS results.

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 occupationnumber 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 demarkation 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 reference [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. References [47,54,68] have investigated the BKT nature of the thermalized state. Wave-turbulence studies [32,41,69] have focussed on powerlaw regions in energy and occupation-number spectra of the type we find in our third regime. The DNS studies in [41][42][43][44][45]70] 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]. Our work suggests that, at least in the dissipationless, unforced, Fourier-truncated, 2D, GP equation, the only robust power laws in spectra are the 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 Fouriertruncated, 2D, GP equation, there are initial-condition-dependent transients. In the second regime the energy and the occupation-number spectra start to develop powerlaw 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 c kin (k) and E int (k) + E q (k) ∼ k, and n(k) ∼ 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 α or logarthimically for different intial conditions (table 2). In the fourth, complete-thermalization regime, power-law forms of correlation functions and spectra, for E < E BKT , are consistent with their nontrivial BKT forms; however, considerable care must be exercised, as explained in section 3.4.1 and [47,54,68], to distinguish these nontrivial power laws from their wave-turbulence analogs [32,41,69].
hereT BKT denotes the estimate for T BKT that follows from an energy-entropy argument [52]. For T < T BKT , the phase correlation function c(r) (see (17)) and the angle-integrated spectrumĉ(k), which follows from a Fourier tranform of c(r), scale as respectively. Above T BKT the correlation length is finite; and, as T → T BKT , it displays the essential singularity We now develop an analytical framework, which is valid at low-temperatures T ≪ T 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 where Ξ is the grand partition function, β the inverse temperature, µ the chemical potential, and N the number of bosons. The grand-canonical potential is We adapt to 2D the 3D study of Ref. [51], expand ψ in terms of Fourier modes A k , and obtain Ω as the sum of the saddle-point part Ω sp and Ω Q , the deviations from the saddle point that are quadratic in A k . We write Ω = Ω sp + Ω Q , where Ω sp = −Aµ 2 /2g and We can also calculate the condensate depletion δN, where the particle number N = N 0 + δN and N 0 is the number of particles in the k = 0 mode, as follows: The integrals in the (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 Ω ph is finite at p = 0 and for δN it is infra-red (I.R.) divergent. By subtracting the I.R. finite and divergent terms from Ω Q and δN, respectively, we get the following expressions, in 2D, in the thermodynamic limit A → ∞: We next determine the chemical potential µ, which fixes the total density ρ = mN/A at a given value, by solving the equation We use this low-temperature result (B.25) to estimate the inverse-temperature scale β BKT , at which the depletion of the k = 0 condensate mode becomes significant for a finite-size system with N 2 c collocation points (which fixes the maximum momentum p max ); in particular, we can solve (B.25), for δρ/ρ = 1, to obtain By making the replacements that correspond do defining , m, and g in terms of c and ξ, as in [51], p max → k max , → √ 2cmξ, and g → c 2 m 2 /ρ, we can rewrite (B.26) as Our DNS runs, which use initial conditions of type 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 . (B.28) The energy of a pure condensate is E 0 = lim β→∞ E = g ρ 2 A 2 m 2 ; (B.29) and the energy and the inverse temperature β (B.28) can be related as follows: where δE is the relative increase of energy above E 0 , and  .

(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 . Hence, the energy of a pure condensate is obtained, in the dimensionless form, by dividing (B.29) by , which gives (B.35)