Numerical method for the stochastic projected Gross-Pitaevskii equation

Numerical method for the stochastic projected Gross-Pitaevskii equation


We present a method for solving the stochastic projected Gross-Pitaevskii equation (SPGPE) for a three-dimensional Bose gas in a harmonic-oscillator trapping potential. The SPGPE contains the challenge of both accurately evolving all modes in the low energy classical region of the system, and evaluating terms from the number-conserving scattering reservoir process. We give an accurate and efficient procedure for evaluating the scattering terms using a Hermite-polynomial based spectral-Galerkin representation, which allows us to precisely implement the low energy mode restriction. Stochastic integration is performed using the weak semi-implicit Euler method. We extensively characterize the accuracy of our method, finding a faster than expected rate of stochastic convergence. Physical consistency of the algorithm is demonstrated by considering thermalization of initially random states.



Providing a quantitative description of non-equilibrium dynamics of Bose-Einstein condensates (BECs) at finite temperature is an ongoing challenge [1]. Classical field methods have been a popular tool to describe finite temperature BECs, utilizing the tractability of the Gross-Pitaevskii equation (GPE) to simulate the dynamics of many highly occupied modes of the system [3]. These highly occupied coherent modes are treated with a classical field approximation, enabling the dynamics of these modes to be treated nonperturbatively. These methods have been used to treat of both equilibrium and non-equilibrium systems in a range of finite temperature systems [5].

Introducing coupling between the coherent region and the remaining incoherent reservoir results in a stochastic GPE, which includes damping and noise terms from the reservoir interaction. Such a description has been derived microscopically [16], allowing for these methods to provide an ab initio description of non-equilibrium dynamics. There are now numerous examples in the literature of the application of stochastic GPEs to a range of systems, including vortex decay [20], soliton decay [23], defect formation across phase transitions [25], spinor condensates [30], polariton condensates [34], equilibrium properties [35], and low dimensional systems [36].

Currently, most applications of stochastic GPEs have only included growth processes, where collisions between two incoherent region atoms lead to a change in population of the coherent region. The damping term that arises from this process is similar to that of the damped GPE [46], and is relatively simple to implement numerically. In the stochastic projected GPE (SPGPE) of Gardiner et al. there are scattering reservoir processes, where a collision between coherent and incoherent atoms results in energy change with no population transfer. These scattering terms have been recently implemented in Ref. [51], showing a dominant effect on highly non-equilibirum dynamics.

In this paper, we present a numerical method for evolving the SPGPE, including the scattering terms. Numerically solving the SPGPE involves two major technical challenges: (i) All moderately occupied coherent modes play an important role in finite temperature non-equilibrium dynamics. Thus all modes beneath a well-chosen cutoff must be propagated accurately. (ii) The deterministic term arising from the reservoir interaction is non-local, while the noise is multiplicative and spatially correlated. Thus accurately and efficiently implementing the scattering terms is a challenge. Previous work has shown how a spectral approach [52] can be used within the PGPE formalism to precisely implement the energy cutoff, and provide a means for accurate evolution of all low energy coherent modes for a Bose gas with both contact [54] and dipolar interactions [55]. In this work we extend these methods for the PGPE, providing a numerical method for evolving the full SPGPE that evaluates the scattering terms, while still propagating all coherent modes accurately. We thus give a complete spectral-Galerkin method for the SPGPE.

This paper is organized as follows: In the rest of this section we briefly review the SPGPE formalism. In Sec. ? we outline our spectral approach, and outline the equations for the mode amplitudes we need to evaluate to solve the SPGPE. In Sec. ? we present our algorithm for evaluating the scattering terms using a Gauss-Hermite quadrature approach. In Sec. ? we characterize the accuracy of our approach, with results demonstrating the convergence with quadrature grid size and evolution time step size. We present an example of the usage of the SPGPE, showing the evolution of a random initial state to equilibrium.

1.2SPGPE theory

We briefly outline the key formalism of the SPGPE relevant to our work (see Refs. [16] for full details). Here we consider a system of bosons confined in a three-dimensional trapping potential, described by the dimensionless single-particle Hamiltonian

where , and

is the dominant contribution to the single particle Hamiltonian and defines our c-field basis, any time dependent perturbation potential is included in , and is the relative trap frequency in each direction . Note that we work in dimensionless units throughout this paper, using harmonic oscillator units of length , energy , and time , where is the particle mass, and is a chosen reference frequency.

In our c-field approach we divide the modes of our system into two subspaces according to their occupation, which we refer to as the coherent (C) and incoherent (I) regions [6]. The low energy C region contains highly occupied modes, which are dominated by classical fluctuations. In this case, we may use a classical field approximation to treat these highly occupied modes [1]. The remaining I region contains sparsely occupied modes, and plays the role of a thermal reservoir. The SPGPE is a stochastic equation of motion for the C-field, which accounts for the reservoir interactions from the I region. Treating the I region semi-classically and assuming spatially constant reservoir interaction rates [25], the evolution equation is given in Stratonovich (S) form by [51]


where the evolution of the c-field is formally restricted to the C region by the projector operator

where are eigenstates of the single-particle Hamiltonian satisfying , and the summation includes all modes the C region defined as

where is the single-particle energy cutoff defining the C region.

The first term of the SPGPE (Equation 2) describes the Hamiltonian evolution, where

is the Hamiltonian evolution operator for the C region, where the dimensionless nonlinearity constant is , where is the s-wave scattering length. On its own Eq. (Equation 2) is the PGPE, a nonlinear Schröndiger equation for the C-field as an isolated microcanonical system [1].

The remaining terms in the SPGPE account for distinct reservoir interactions.

Growth reservoir interaction

Eq. ( ?) describes the growth reservoir interaction, where two I-region atoms collide leading to population growth of the C-region. The rate of this process is set by [25], and the Gaussian complex noise has the non-vanishing correlation


is a delta function in the C region.

Scattering reservoir interaction

Eq. ( ?) describes the number conserving scattering reservoir interaction where energy and momentum are transferred between the C and I regions, without population transfer. This process is described by the effective potential

where , the C-field current is


The real scattering noise has the non-zero correlation

where the correlation is anti-diagonal in -space so can be easily sampled numerically [see Sec. ?].

Scattering SPGPE

By neglecting the growth terms in the SPGPE (Equation 1), we arrive at the scattering SPGPE

In this work we are concerned with developing a method for evaluating the terms involved with the scattering processes. Since including the growth terms is a relatively simple extension of the PGPE algorithm [25], we will consider only the scattering SPGPE only from here onwards.

2Formal algorithm

2.1Spectral representation

We use the single-particle eigenstates as our spectral basis, so we write the c-field as

where are time dependent complex amplitudes and represents all quantum numbers required to specify a single-particle state. This choice is convenient because it allows us to efficiently implement the projection by restricting the spectral modes [as indicated in Eq. (Equation 10)] to the set indicated in Eq. (Equation 3) defining the C region.

2.2Mode evolution

Spectral-Galerkin formulation

We exploit that is diagonal in the spectral basis and use a Galerkin approach [52], where we project the scattering SPGPE (Equation 9) onto the spectral basis. This leads to a system of equations for the evolution of the amplitudes, i.e.


are the nonlinear matrix elements of the two-body interaction, scattering effective potential, and scattering noise terms respectively. We introduce the functions later [see Sec. ?], and is the standard real Wiener process satisfying

There are two main steps in solving this equation: (i) time-evolution to step this equation forward in time [see Sec. ?]; and (ii) evaluating the non-linear matrix elements (Equation 12)-( ?) at each time step [see Sec. ?].

Stochastic time evolution algortihm

Because the noise associated with the scattering reservoir interaction is multiplicative, we use the weak vector semi-implicit Euler algorithm [56] to evolve our stochastic equations forward in time. As this algorithm is extensively discussed in the literature we briefly review the algorithm here. Equation (Equation 11) is of the general form

where , we use the notation to represent the dependence of matrix elements on the full field (), and the noise matrix elements depends on the Wiener process . The solution is propagated to a set of discrete times , where is the step size, and we denote that solution at time as . Using this solution, the solution at the next time-step is computed as , where


Note formally , however in practice we sample as a real Gaussian distributed random variable of variance [c.f. Eq. ( ?)].

3Evaluation of the scattering SPGPE matrix elements

Here we give a full description of our algorithm to efficiently and accurately evaluate the scattering SPGPE matrix elements in the harmonic oscillator basis. This expands on the brief overview of our method presented in Ref. [51].

3.1Harmonic-oscillator state properties

We firstly discuss some important properties of our single-particle basis used in this work.


The eigenstates of the basis Hamiltonian () are separable into one dimensional basis states, that is,

where are eigenstates of the dimensionless 1D harmonic-oscillator Hamiltonian

where the harmonic-oscillator states take the form

where the normalization constant is , and is a Hermite polynomial of degree , defined by the recurrence relation

with and . The eigenvalue is the single-particle energy given by , where is a non-negative integer (we use Greek subscripts to denote 1D eigenstates). The C region is defined by the region containing all modes below the single-particle energy cutoff

so that within the C region there exist distinct 1D eigenstates in each direction, and total 3D basis states in the C region.

For this work we will consider a spherically symmetric system (), to avoid cumbersome notation involving . However the work presented in this paper can be easily generalized for any trap anisotropy, by retaining the general form of Eq. (Equation 17). Thus we drop any reference to , and consider a system with modes in the C region in each direction.

Step operators

Step operators allow us to represent certain operators exactly in the spectral basis. The step operators are defined as

We then find that the matrix representation of the step operators in the spectral basis is

and similarly,

For our purposes, this allows us to differentiate in the spectral basis exactly, since

This 1D procedure is applied in the same manner to include and . It is important that we can apply the derivative operator exactly in the spectral basis, since we need to evaluate the c-field current (Equation 6) to evaluate the scattering effective potential term.

For a further discussion on the use of step operators to calculate observables, see Ref. [54].

3.2Two-body interaction term

Here we briefly review our scheme for numerically integrating the nonlinear contact interaction term (Equation 12) using Gauss-Hermite quadrature, which is described in complete detail in Ref. [54]. This algorithm calculates the nonlinear interaction term exactly for a harmonically trapped system, and is required to integrate the scattering SPGPE (Equation 9).

Firstly, using (Equation 10) and the form of the basis states (Equation 17), we can write the c-field as


is a polynomial of maximum degree in each independent coordinate due to the energy cutoff.

The interaction term (Equation 12) is forth order in the field, so we can write it in the form


is a polynomial of maximum degree in each coordinate. We evaluate (Equation 23) using Gauss-Hermite quadrature. The general form of the point quadrature rule is

where is a Gaussian weight function, and and are the quadrature weights and roots. The Gauss-Hermite quadrature is exact if is a polynomial of maximum degree . Since the exponential in (Equation 23) takes the form of the appropriate weight function for Gauss-Hermite quadrature, we can evaluate the interaction matrix element exactly by

using a three-dimensional spatial grid with points in each direction, where the and are the roots and weights of the one-dimensional Gauss-Hermite quadrature with weight function .

3.3Scattering effective potential term

To compute the effective potential matrix elements ( ?), we adapt the scheme developed for evaluating the dipolar interaction term in the PGPE [55]. Here the calculation involves the Fourier transform of the current (Equation 6), whose momentum space form is then Fourier transformed to form the effective potential (Equation 5). As we are using a nonuniform quadrature grid to represent the c-field, we follow Ref. [55] in using an auxiliary harmonic-oscillator basis to perform the Fourier transforms.

Firstly we calculate the c-field current density (Equation 6) by

where the index represents each spatial dimension with corresponding unit vector , and

where Eq. (Equation 28) can be calculated exactly using step operators [see section ?]. Now using similar arguments to section ?, each component of the current can be written in the form

where is polynomial of maximum degree in each coordinate (due to the projector).

Following Ref. [55], we introduce a set of auxiliary harmonic-oscillator states,

where the exponential argument is chosen to match Eq. (Equation 29). These states are eigenstates of the harmonic oscillator Hamiltonian, with a trapping potential twice as tight as that defining the spectral basis (Equation 16). Noting this, the auxiliary states are related to the spectral basis modes by .

Due to the same exponential factor, we can represent each component of the current [i.e. (Equation 28)] as

where is a set of auxiliary basis coefficients for each component of . Since the auxiliary oscillator states form an orthonormal basis, we can evaluate these basis coefficients by


is a polynomial of maximum degree in each spatial dimension. We can use the same Gauss-Hermite quadrature to integrate (Equation 31) as we used for the nonlinear interaction term (Equation 23), as the weight function and maximum polynomial degree are both . Thus Eq. (Equation 31) can be calculated exactly for each component of the current, by

where the quadrature roots , and weights , are the same as used in (Equation 26).

To evaluate the effective potential (Equation 5), we must Fourier transform each component of the current. Since the harmonic-oscillator states are eigenstates of the Fourier transform operator,

we can efficiently evaluate the Fourier transform with knowledge of auxiliary basis amplitudes of the current by

Thus (Equation 35) will represent the Fourier transform of the current on the quadrature grid exactly. We can now evaluate the scattering effective potential term (Equation 5) by evaluating

and computing the inverse Fourier transform of Eq. (Equation 36). To compute the inverse Fourier transform, we expand in the auxiliary basis


We integrate (Equation 38) using a Gauss-Hermite quadrature, via


Here the quadrature roots and weights are found from the weight function . Since in general cannot be represented exactly in the oscillator basis, our quadrature rule (Equation 38), and thus Eq. (Equation 37), are approximations. Hence the number of -grid quadrature points required is somewhat arbitrary [see Sec. ?]. We investigate how the accuracy of the scattering matrix element depends on the number of quadrature points in Sec. ?.

From the form of Eq. (Equation 37), we see that the effective potential matrix elements we require ( ?) are of the form


where represents the c-field and is given by Eq. (Equation 22). Although our expression for is approximate, once in this form, takes the form of a polynomial of degree in each coordinate. Thus we can integrate Eq. (Equation 41) exactly using the same Gauss-Hermite quadrature sum as for Eq. (Equation 26), i.e. using a three-dimensional spatial grid with points in each direction, where the and are the roots and weights corresponding to the weight function .

3.4Scattering noise term

Correlation function

A noise with the necessary correlations (Equation 8) can be constructed as

Here we have made use of the Fourier transformed modes . The phase factors arising in Fourier transforming the modes generates the antidiagonal delta-function we require, since

where we have used Eq. (Equation 4).

Computing the matrix elements

To implement the noise in our Galerkin approach, we need to evaluate the matrix elements corresponding to Eq. ( ?). This means we need to project the noise in position space onto our spectral basis. Since numerically we form the correct scattering noise correlation in Fourier space [see Eq. (Equation 42)], we need to compute

where the are real Gaussian distributed random variables [see Eq. (Equation 13)-( ?)], and where


where the are the Fourier transformed basis states. Thus to evaluate the scattering noise matrix elements ( ?), we must apply to the Gaussian distributed random variables. We now outline how we evaluate the scattering noise matrix elements.

Firstly, we can represent Eq. (Equation 45) in the form


In order to compute the inverse Fourier transform of Eq. (Equation 46) to form , we use a similar process as in Sec. ?. Here we do not need to use the auxiliary basis states necessary for computing , since the exponential term in Eq. (Equation 46) matches that of the standard basis states. Thus, we can calculate using the expansion of in the standard spectral basis


Eq. (Equation 48) is approximate [c.f. Eq. (Equation 37)] since we can not represent as a polynomial. We evaluate the integral using a Gauss-Hermite quadrature


where the quadrature roots and weights correspond to the weight function . Since our Gauss-Hermite quadrature is approximate, the number of points is again somewhat arbitrary here. See Sec. ? for further discussion.

Finally, since the matrix elements we require ( ?) are third order in the field, they can be cast in the form


where represents the c-field and is given by Eq. (Equation 22), so is a polynomial of degree in each coordinate. Thus we can evaluate Eq. (Equation 51) using a Gauss-Hermite quadrature

This quadrature is exact using a three-dimensional grid with quadrature points in each coordinate, where the grid points and quadrature weights correspond to the weight function .

3.5Summary of the algorithm

Here we give a summary of our algorithm for calculating the scattering SPGPE matrix elements. We split our summary into two sections, first dealing with the deterministic terms followed by the noise term.

Fourth order field matrix elements - and

For each Euler step, the procedure for calculating these matrix elements is as follows:

  • We firstly transform the c-field from the spectral to spatial representation via

    where , are the quadrature points associated with the weight function , and are the precomputed transformation matrices. These are the one-dimensional basis states evaluated at the quadrature grid points

  • Differentiate the field

    where the derivative operators are given by Eq. (Equation 20), and the are the identity operators.

  • Transform the differentiated fields to position space

    so we can then form the position space current current via Eq. (Equation 27), where we have

  • The weighted position space current is constructed for each component of

    where are the quadrature weights associated with the weight function .

  • We compute the Fourier transform of each component of the current by

    Here the precomputed transformation matrix is

    where the -space quadrature grid is based on the weight function . Eq. (Equation 56) combines both steps of the Fourier transform into one operation [, and ].

  • In Fourier space we form the projection of the current onto . We then premultiply with the quadrature weights to form

  • Inverse transforming takes us back to position space. and computes the effective potential :

  • We combine the two-body interaction and scattering effective potential terms into a single integrand (since the quadrature sum has the same weight function for both cases),

  • Inverse transforming completes the integration of Eq. (Equation 12) and Eq. ( ?), giving the required matrix elements

Third order field matrix elements - Scattering noise

To implement the noise, we generate at the start of each Euler step. Then for each iteration of Eq. (Equation 14), we perform the following steps:

  • We transform the random variables to Fourier space via

    where the precomputed transformation matrices are the Fourier transformed basis states

    where the tilde notation implies the -space quadrature grid is based on the weight function .

  • We form the appropriate correlation function in -space, and then compute the inverse Fourier transform:


    where the hat notation implies the position space quadrature grid is based on the weight function .

  • We now need to transform the c-field to position space, onto the grid based on the three-field weight function:

    where now the precomputed transformation matrices are

  • We form the appropriate integrand of ,

  • The inverse transform of completes the required integration

4Accuracy of algorithm

In this section, we discuss how we can control the accuracy of our algorithm, and quantify the accuracy of our approach of both spatial and temporal integration using various measures.

4.1-space quadrature grid

Both the deterministic and noise terms involve -space integrals which are evaluated as approximate quadrature sums [evaluating Eq. (Equation 57), and Eq. (Equation 60) respectively]. To control the accuracy of these steps, we can vary the size of the appropriate -space quadrature grid.

For a given , the transform to momentum space is exactly invertible if the quadrature grid sizes are chosen by , , where . We choose an even number of -grid points to avoid a quadrature point at where both the deterministic and noise terms are singular. To investigate the effect of the -grid on the accuracy of our method, we vary the number of quadrature points used by

where is the total number of quadrature grid points, and is the number of quadrature points added to the reference value . Eq. (Equation 61) refers to the number of quadrature points used to evaluate the deterministic term, i.e. the quadrature points are based on the weight function [see Eq. (Equation 40)]. Eq. ( ?) refers to the quadrature grid for evaluating the noise term, i.e. the quadrature points are based on the weight function [see Eq. (Equation 52)]. For the remainder of this section we associate and with the quadrature grid for the deterministic and noise terms respectively.

4.2Effective potential term convergence for a breathing mode

A Gaussian wave function undergoing radial breathing provides a rare example where we can analytically calculate for a case where . Here we consider the Gaussian wave function

which corresponds to a radial breathing oscillation. This wave function has a non-zero current density, given by

from which we find [using Eq. (Equation 5)] the radially symmetric potential


with . The potential is most significant at , where

To quantify the accuracy of our numerical evaluation of , we use the relative error measure

In Figure 1 we plot as a function of , for the Gaussian wave function [Eq. (Equation 62)], with and . We show the relative error for evaluated with and , corresponding to and respectively. In both cases rapidly converges with