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 and partial thermalization takes place for modes with ; the self-truncation wave-number depends on the initial conditions and it grows either as a power of or as . 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.
pacs:47.27.Ak, 47.27.Gs, 47.37.+q, 67.25.dk, 67.25.dj
The elucidation of the nature of superfluid turbulence, which began with the pioneering studies of Feynman  and of Vinen and Hall [2, 3, 4, 5, 6], has continued to engage the attention of experimentalists, theoreticians, and numerical simulators [7, 8, 9, 10, 11, 12, 13]. Experimental systems, in which such turbulence is studied, include the bosonic superfluid He, its fermionic counterpart He, and Bose-Einstein condensates (BECs) of cold atoms in traps and their optical analogues; for representative studies we refer the reader to [14, 15, 16, 17, 18, 19, 20, 21, 22, 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 (NLS) equations [30, 31]. These models have been studied by a combination of theoretical methods, such as wave-turbulence theory [32, 30, 33, 31], and numerical simulations [34, 35, 36, 37, 38, 39, 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 past few years [41, 42, 43, 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, though not all, parts of our results are contained in earlier simulations [41, 42, 43, 45, 44, 46, 47]. 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 [48, 49, 50]; a similar study for the 3D GP equation has been carried out by Krstulovic and Brachet [38, 51].
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 power-law 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 initial-condition-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 [52, 53, 54, 47], 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 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 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 of a weakly interacting 2D Bose gas at low temperatures, is
and the total number of particles
where is the area of our 2D, periodic, computational domain of side . From (1) we obtain the continuity equation
where is interpreted as the particle density and the velocity is
We can use the Madelung transformation , where is the phase of , to write , whence we get 
where the kinetic, interaction, and quantum-pressure energies are defined, respectively, as
We separate the compressible (supercript ) and the incompressible (superscript ) parts of the kinetic energy by making use of the decomposition
where and , whence we obtain the following:
The spectra for these energies are defined as follows:
furthermore, we define an occupation-number spectrum via
here we denote the Fourier transform of by ; and, for notational convenience, we do not show explicitly the dependence of these spectra on time . 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).
where the sound velocity is 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 
for temperatures below the transition temperature (or energy in the microcanonical ensemble),
where and the critical exponent for ; and at . The BKT phase shows bound vortex-antivortex pairs; these unbind above , so
in the disordered phase, with the correlation length.
2.2 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 on a square simulation domain of side with 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 , 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 . Thus, the maximum wave number , where , 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 , , and , always normalized to correspond to a total number of particles (3) . The first of these is best represented in Fourier space as follows:
where , are random numbers distributed uniformly on the interval ; and , where the integer controls the spatial scale at which energy is injected into the system, and the real number specifies the Fourier-space width of at time . The initial condition is like but, in addition, it has a finite initial condensate population at time .
We obtain the initial condition by solving the 2D, stochastic, Ginzburg-Landau equation (SGLE), which follows from the free-energy functional
where is the chemical potential222Recall that the SGLE can be thought of as an imaginary-time GP equation with external, additive noise (see, e.g. reference ). The SGLE is
where is a zero-mean, Gaussian white noise with
which we solve along with the following, ad-hoc equation
to control the number of particles ; the parameter controls the mean value of ; and 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 , for time marching and the method of reference  (see page 25 of this reference).
We first present the time evolution of the different energies, the probability distribution functions (PDFs) of the velocity components, and the population 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 PDFs, and the zero-wave-number population
We show the early stages of the time evolution of the energies , , , and , from our DNS runs -, , and in figure 1. The runs - use initial conditions of type , in which is a significant fraction of the total initial energy; the runs and start with initial configurations of type and , respectively, in which is negligibly small at . 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 into the other three components, with a significant fraction being transferred to ; 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 -; we observe, in particular, that the conversion of into the other energy components is accelerated as increases from to (cf. ); and there is a corresponding acceleration in the approach to thermalization. Moreover, the larger the value of the larger is the time required for thermalization, as we can see by comparing figures 1 (a) and (d), for the runs and , respectively; the run starts with a high value of 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 -. In contrast, the runs and have negligibly small values of to begin with (figures 1 (e) and (f), respectively); and remains close to zero throughout the dynamical evolution here. For run , both and start from values close to zero, grow at the cost of , and finally saturate to small, statistically steady values. For run , 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 and , the Cartesian components of the velocity, for our DNS runs , , and , which correspond, respectively, to initial conditions of types , , and . For the run , 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 , with , and or (we show fits only for ). Such power-law tails in velocity-component PDFs have been seen in experiments  and some numerical studies [39, 65, 66, 45]. 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 increases, from PDFs with power-law tails (figure 2 (a) for run ), to ones with a Gaussian form near the mean, followed by broad tails (figure 2 (b) for run ), and then to more-or-less Gaussian PDFs (figure 2 (c) for run ), 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 and , both of which have a negligibly small value of at , are given, respectively, in figures 2(d)-(f) and figures 2(g)-(i).
We turn now to the time evolution of the population , in the mode [36, 40, 67], and its dependence on the initial conditions. In figure 3 (a) we plot versus for the runs - (red, blue, green, and brown curves, respectively), which use initial configurations of type ; these figures show that increases with , on average, and depends on , , , and . For the runs and (red and blue curves in figure 3 (a)), approaches a saturation value for the time scales probed by our simulations; figure 3 (a) also shows that, as we increase (red, blue, and green lines in figure 3 (a)), the fluctuations in are enhanced and its large- value, which it seems to approach asymptotically, diminishes. By comparing the runs and (red and brown lines in figures 3 (a)), we see that the latter has a higher value of than the former, because both and are smaller for than for ; thus, grows more slowly in than in ; and, after an equal amount of simulation time, its value in is nearly an order of magnitude lower than in ; the former shows large fluctuations in and no sign of saturation. The run (figure 3 (e)) uses an initial configuration of type , with a large value of ; in this case, after a period of initial transients, over our simulation time. The run (figure 3 (f)) uses an initial condition of type ; here fluctuates slightly but remains close to its initial value (cf. [40, 67]).
To study the dependence of on the number of collocation points , we evolve the initial configuration of for (run ), (run ), (run ), and (run ). Figure 3 (g) shows plots of versus for these five runs; clearly, the initial evolution of depends significantly on ; however, the large- values of , on the time scales of our runs, are comparable () for the runs wth (run ), (run ), and (run ). In contrast, the saturation value for the run with (run ) is . For the run (), shows large fluctuations and no sign of saturation over the time scale that we have covered; this suggests that also depends on the realisation of the random phases in (21). These plots of illustrate that complete thermalization proceeds very slowly for ; in the completely thermalized state of the Fourier-truncated, 2D, GP system, 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 [68, 67].
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 , , and . The first type begins with a sizeable incompressible kinetic energy spectrum ; and the initial transients are associated with the annihilation and creation of vortex-antivortex pairs, the associated depletion of , and the growth of the other energy components . In contrast, runs with initial conditions of types and 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 . In figures 4 (a)-(d) we show the time evolution of the spectra , for the runs , , , and , to ascertain the presence of scaling behaviour, if any. We find that, in the low- region, lacks a well-defined scaling region (unlike in ); indeed, this region depends on the initial configuration, changes continuously with time, and, in particular, a 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 for the run ). At large wave numbers, , during the initial stages of evolution, because of the presence of the vortices ; this power-law form holds over the same time scales for which the PDF (figures 2 (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 . Figure 5 (a) shows for the run ; we begin to see a power-law region here with , on the low- side of the peak after which the spectrum falls steeply. A similar behaviour starts to emerge in the region for the run (figure 5 (g)). In this onset-of-thermalization regime, we also see the development of the following power laws: (figure 7) and (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 , and figures 5 (h) and 6 (e), and (f), for initial conditions of type . It is important to distinguish between (I) spectra that fall steeply at large values of , 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 , 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 in , referred to as the self-truncation wave number , moves out to as a power of and (B) a subclass in which moves out to at a rate that is slower than a power of .
Figures 5 (g)-(i), from the run , show how evolves as the spectral convergence to the GP PDE is lost in case (II); note that the scaling region with sets in at high wave numbers close to and then extends to the low-wave-number regime. For case (IA) analogous plots of 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 . Illustrative plots of the spectra and in this regime of partial thermalization are given in figures 7 and 8, respectively.
We now present a detailed characterization of the partial-thermalization regime, when energy spectra display self-truncation at wave-numbers beyond , which can be defined as follows:
as the system approaches complete thermalization, . In particular, we explore how the scaling ranges in energy spectra grow with for different values of , with the initial configuration and number of collocation points held fixed. For an initial condition of type , with , , and , we obtain the time evolution of energy spectra for (run ), (run ), and (run ) 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 , the more rapid is the thermalization, and the consequent loss of spectral convergence, as we can see by comparing the sky-blue (run ), green (run ), and purple (run ) spectra in figures 9 (a)-(c); run loses spectral convergence around . We obtain the same qualitative dependence, with , , and , for , , and , i.e., runs , , and , 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 dependence of the self-truncation of energy spectra, for initial conditions, with , , and , and five different values of , namely, (run ), (run ), (run ), (run ), and (run ). We find, not surprisingly, that the lower the value of the more rapidly does the system lose spectral convergence.
Initial conditions of type lead to energy spectra whose time evolution, and their dependence on and , is similar to those that are obtained from intial conditions of type .
With initial conditions of types and , we cannot control the initial value easily. However, initial conditions of type , which we obtain from the SGLE, allow us to control and start, therefore, with initial spectra that display partial thermalization for  and a sharp fall thereafter. In figure 10 we show the time evolution of for such initial conditions from runs -. For different representative values of , , and , we now study the time evolution of , which characterizes the growth of the partially thermalized scaling region. Here too, as with initial conditions of types and , if all other parameters like and are held fixed, the speed of thermalization increases with (cf. figure 10 (a) for the run , with , and figure 10 (b) for the run , with ). For these runs -, the growth of the energy spectra, in the region , starts with the smoothening of the sharp cut-off at ; the higher the value of , the slower is this growth (cf. figures 10 (b), (d), (e), and (f) for runs , , , and , respectively). By contrast, an increase in (or ) in the SGLE, accelerates this growth (cf. figures 10 (b) and (c) for runs and , respectively).
The growth of with , illustrated in figure 11 (a), can be fit to the form ; however, as we show below, depends on the initial condition. We obtain the exponent either from slopes of log-log plots of (i) versus or (ii) versus ; we denote the values from procedures (i) and (ii) as and , 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 . The values of and , listed in table 2, show that ; the discrepancy between these two values for is a convenient measure of the errors of our estimates. For runs , , and , we cannot obtain reliably; the small values of for these runs indicate very slow growth of ; indeed, in runs and , a case can be made for a logarithmic growth of with .
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 for the runs and , where shows power-law scaling over the entire wave number range, from up to , towards the end of the respective simulations; a naïve fit is consistent with (but see below).
3.4.1 Correlation functions and the 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 [52, 54], with the correlation function , at energies ; and should decay exponentially with if . We show this explicitly now by using initial conditions of type with and and ; we obtain different energies by changing and (runs D1-D13 and E1-E12 in table 3).
In figure 12, we present plots of the correlation functions . To illustrate the BKT transition clearly, we present log-log plots of versus , for , in figures 12 (a) and (d), where the straight lines indicate power-law regimes; and, for , we use semi-log plots, as in figures 12 (b) and (e), where the straight lines signify an exponential decay of with . Given the resolution of our DNS runs, we find that, in a small energy range in the vicinity of , we cannot fit power-law or exponential forms satisfactorily; this leads to an uncertainty in our estimate for . Aside from this uncertainty, the behavior of , 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 for as shown in figures 12 (c) and (f). Our values for , for the runs with and with and , are listed in table 3. Note that and , i.e., depends on , the number of collocation points; we show analytically below how a low-temperature analysis can be used to understand this dependence of on . In the completely thermalized state of the Fourier-truncated, 2D, GP system, must vanish in the thermodynamic limit by virtue of the Hohenberg-Mermin-Wagner theorem [57, 58] and ; 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 in figure 3); however, finite-size scaling can be used to extract the exponent from the part of as shown in reference ; similarly, 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.
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 . We define the energy of our system to be ; this energy is fixed by the initial condition; and measures the relative amount by which exceeds . As we show in the B, the dependence of the energy , at which the BKT transition occurs, can be obtained approximately as follows. We begin with
whence we obtain
We can now write
by using this expression we can determine the ratio for runs with two different values, and , for the number of collocation points; we can also obtain this ratio from our DNS, by determining the value of at which the exponent becomes . In Table 4 we compare for and ; our analytical approximation (31) yields ; this is in excellent agreement with the value that we obtain for this ratio from our DNS results.