The Stochastic Gross-Pitaevskii Equation and some Applications

Abstract

The stochastic Gross-Pitaevskii equation represents a versatile approach for studying the
dynamics of trapped degenerate ultracold Bose gases in the presence of large phase and density fluctuations. Following a brief review of the original formulation of Stoof, which highlights the benefits of this approach and its relation to alternative theories, we present selected applications for the dynamics of effectively one-dimensional systems, and briefly discuss the generalization to two-dimensional systems, highlighting certain potential pitfalls in their numerical implementations.

The ultracold atomic samples now routinely produced in many laboratories worldwide, provide excellent means by which to investigate many-body quantum systems on a macroscopic scale. From a theoretical perspective, these systems are appealing as they support an array of rich dynamics, in many cases readily amenable to experimental observation. With this close relation between theory and experiment (1); (2), many such experimental features have been subject to analytical and numerical investigations: much attention has been paid to numerical analyses based upon the Gross-Pitaevskii equation (GPE) in particular (3), which has been shown to successfully describe condensate dynamics in the limit when the many-body wavefunction represents the condensate alone, thus corresponding to the low-temperature regime.

Most experiments on trapped ultracold gases feature, to a greater or lesser extent, a depletion of the condensate due to the effect of finite temperature. While in three dimensional systems, finite temperature effects play an increasingly important role as the critical temperature is approached, for low dimensional systems, thermal effects harbour greater importance still, due to the existence of enhanced phase fluctuations within systems under tight confinement in one or two directions
(4); (5); (6). Moreover, the dynamics of the thermal component can play a significant role, for example in the case of an ultracold gas under strong external perturbations (7); (8); (9); (10).

A number of finite-temperature techniques have been developed to incorporate the effects of temperature into a description of trapped Bose gases, as recently summarised in (11); (12); (13); (14), which feature a detailed comparison of the relative weaknesses and merits of such approaches.
In brief, various theories have been put forward and applied to the experiments, as discussed below: In mean-field type approaches, the condensate is treated as a classical object, separated from the effect of the thermal cloud; various such perturbative approaches have been formulated to second order in the effective interaction strength
(15); (16); (17); (18), with the common feature of obtaining a set of coupled dynamical equations for the condensate and the non-condensate, which interact both via mean fields and by particle exchange.
The condensate obeys an appropriately generalised GPE, which,
in the most commonly adapted version following from early work of Kirkpatrick and Dorfman (19); (20); (21), is coupled to a quantum Boltzmann equation for the thermal cloud; this model has been extensively used by Zaremba, Nikuni, Griffin and co-workers, who introduced the terminology ‘ZNG’ (15).
Although such approaches work remarkably well in predicting a large range of experimentally-relevant phenomena, they do not fully describe regimes where fluctuations in the condensate phase are large, as is appropriate for example in low-dimensional systems
(4); (5); (6).
To incorporate such effects, generalised zero temperature (25) and finite temperature (26) modified methods were developed, to describe the static properties of such systems in an ab initio manner.

However, in order to go beyond mean field effects in a dynamical manner, various alternative approaches have been developed and implemented. In first instance, the time-dependent mean field approach mentioned above was reformulated in a somewhat more justified number-conserving manner (27), following from early work (28); (29); (30); although this has led to good agreement with experiments in at least one case (31), a full dynamical model of such theories is still lacking.

Classical field approaches are based upon the observation that (32); (33), as an equation for a classical field, the GPE should be able to describe the classical aspects of an ultracold trapped Bose gas. This includes all highly occupied modes, for which it is possible to approximate the field operators for the creation/annihilation of particles, by c-number amplitudes (34). Consistency requires only the highly occupied modes to be propagated this way, which can be achieved by introducing a projector to the GPE (35). Temperature effects can be incorporated through the initial conditions, which are typically set to some non-equilibrium distribution, and subsequently allowed to equilibrate at some temperature(34), via the (projected-)GPE, under the restriction of fixed particle number and energy.

A variant of this technique is known as the Truncated Wigner method, originally developed in the context of quantum optics and first used for Bose gases in (36); a key element of this method is the approximate incorporation of quantum effects,
through the addition of a prescribed amount of quantum noise to the wavefunction initial conditions.
The numerical procedure which results, is to propagate a set of stochastic matter fields, which serve to sample the Wigner distribution under this approximation, via an ‘appropriate’ equation of motion for the Wigner distribution. The terminology ‘truncated’ Wigner arises from the fact that the propagation is carried out in an approximate manner, via the (projected) GPE.
A limitation of this approach is its inability to accurately describe thermal effects, although certain remedies of limited applicability have been proposed (37). In the latter two approaches, a temperature can only be assigned at the end of the simulations, with the thermalisation arising respectively due to ergodicity (34), while simultaneously exhibiting an unphysical heating.

To gain better control over the temperature to which the system equilibrates, requires detailed consideration of the coupling between such modes, and the higher lying, thermal region in our system. Within such schemes, an appropriate generalisation has been formulated in (38). Earlier, rather distinct, yet analogous in spirit,
means of achieving this was provided independently by Stoof
(39); (40); (41); (42); (43), and by Gardiner, Davis and co-authors (44); (45); (46); (47); (54), using very different theoretical approaches.

This led to the generalisation of the GPE to a Langevin equation (41), often referred to as the Stochastic-Gross-Pitaevskii equation (SGPE) (47), a terminology that we also adopt here; these two independent formulations have many formal similarities and are closely related, up to some subtle difference discussed in (12); (47); (45). The SGPE describes the stochastic evolution of the condensate plus highly occupied, low-lying thermal modes, under the influence of coherent and incoherent interactions with the particles occupying higher energy modes.

Given the diverse range of approaches currently available, our aim in this manuscript is to highlight the key elements of this latter approach, paying particular attention on some of its appealing features via somewhat idealised experimental scenarios. We start our discussion with a brief overview of one formulation of this approach, which has already attracted a relatively large literature
(41); (42); (49); (50); (51); (52); (53), with related work discussed in (47); (54); (55); (56); (57).

ii.1 The Stochastic Gross-Pitaevskii equation

The stochastic Gross-Pitaevskii equation is a non-equilibrium microscopic theory to describe the evolution of an ensemble of ultracold atoms in contact with a thermal cloud. Using the Keldysh non-equilibrium formalism, Stoof derived a Fokker-Planck equation for the the system in (40); such an equation is equivalent to the following Langevin equation (41),

iℏ∂Φ(r;t)∂t=[−ℏ2∇22m+Vext(r)−μ−iR(r;t)+g|Φ(r;t)|2]Φ(r;t)+η(r;t),

(1)

where here Φ(r;t) is the order parameter for the lowest, coherent system modes
and the higher modes are assumed to obey their own coupled dynamics, being in general governed by a quantum Boltzmann equation as in (15); (16); (12).
In addition to the usual Gross-Pitaevskii terms for the kinetic energy, trapping potential (Vext(r)) and nonlinear interaction term (where g=4πℏ2a3d/m, parameterised by the s-wave scattering length a3d), this equation features two additional contributions: firstly, R(r;t) denotes a dissipative term which represents the coupling between the condensate and a thermal reservoir at a temperature T, and accounts for the transfer of particles between the high and low energy modes of the system; the contribution η(r;t) is a noise term which accounts for the random nature of incoherent collisions within the system and is somewhat analogous to the random particle jitter in Brownian-type motion (see, e.g. (42)).
Finally, μ denotes an effective chemical potential.
The above noise contribution is characterized by Gaussian correlations of the form

⟨η∗(r,t)η(r′,t′)⟩=i(ℏ2/2)ΣK(r,t)δ(t−t′)δ(r−r′),

(2)

where ⟨…⟩ denotes averaging over different realisations of the noise (and it is generally understood that the results of simulations will undergo appropriate averaging, although single-runs can also offer qualitative results as discussed later).
The amplitude of the noise is in general position and time dependent, and is governed by the so-called Keldysh self-energy ΣK(r,t). For a thermal cloud which is sufficiently close to equilibrium, this takes the form (41)

ΣK(r,t)

=

−i(4πℏ)g2∫dp2(2πℏ)3∫dp3(2πℏ)3∫dp4(2πℏ)3

(3)

×

(2πℏ)3δ(p2−p3−p4)δ(εc+ε2−ε3−ε4)

×

[(1+N2)N3N4+N2(1+N3)(1+N4)].

Although the occupation numbers, Ni, and therefore also energies εi are in general time-dependent (their values governed by a quantum Boltzmann equation and self-consistent Hartree-Fock energies (41); (12)), current numerical implementations only focus on near-equilibrium situations, for which the thermal cloud is treated as static; in this case, Ni can be approximated by the usual Bose-Einstein distributions Ni=[exp(β(εi−μ)))−1]−1, where β=1/kBT and pi labels the momentum of a particle in the ith single particle energy level.

As one might expect, the strength of the dissipative effects arising from dynamical particle exchange between the two subsystems is ‘balanced’ by the magnitude of fluctuations present, a relation encapsulated by the fluctuation-dissipation relation for the system, upon noting explicitly the dependence of these quantities on the system energy εc, namely
iR(r;εc)=−(1/2)ℏΣK(r;εc)[1+2N(εc)]−1.
In order to make the solution of (1) numerically tractable, in first instance one can allow the system to equilibrate to a classical distribution by noting that for large occupation numbers, [1+2N(εc)]−1≈1/2β(εc−μ).
Since the condensate energy is still an operator of the form εc=(−ℏ2∇2/2m+Vext(r)+g|Φ(r;t)|2),
the Langevin equation takes the somewhat simplified form

The very first numerical implementation of this SGPE in the context of ultracold Bose gases was undertaken in 2001 by Stoof and Bijlsma (41) to model the reversible growth of a condensate, as seen in the MIT experiment of Stamper-Kurn et al. (58): in this experiment, a dimple potential was introduced in a reversible manner to the centre of a harmonic trap, thereby inducing localised phase space compression, and leading to a (periodic) crossing through the phase transition.
Good qualitative agreement was demonstrated between the numerical results (41) and the observed lagging behind of the condensate growth, relative to the addition of the dimple trap, found in the experiment. The method, which is also amenable to analytic variational calculations of stochastic dynamics (42); (51),
has since been applied (in collaboration with one of us, NPP) to numerical studies of coherence
(53); (49), atom laser (50), and atom chip dynamics (52).
It should be noted that related work has been recently undertaken by Gardiner, Davis and co-workers (46); (47); (54); (55); (56); (57).

The aim of this paper is to briefly review the main concepts and illustrate the versatility of the SGPE formalism for describing systems in which fluctuations are of great importance; for numerical convenience, most of our discussion is restricted to effectively one-dimensional systems, for which we also consider the reduction of the SGPE to simpler theories, which incorporate finite temperature effects, to differing levels of approximation. The extension of this approach to two-dimensional systems is also briefly reviewed. Although, by construction, the results of such simulations are best interpreted after numerical averaging over the different noisy initial conditions, we shall see that important information is also contained in single runs, as also noted in (56); (57).

iii.1 Spatiotemporal Evolution of Coherence in a One-dimensional Bose Gas

Non-equilibrium Coherence

Low dimensional systems exhibit a richer phase diagram than their three-dimensional counterparts: in the case of a purely 1d homogeneous system, long wavelength fluctuations in the phase (59), prevent the onset of Bose-Einstein condensation, even down to zero temperature. For a trapped, one-dimensional gas, however, the long-wavelength physics is altered due to the low energy cut-off imposed by the harmonic trapping potential, and the quantum degenerate regime is instead described by two characteristic temperatures, Td and Tϕ, below which density and phase fluctuations, respectively, become suppressed (25). The intermediate regime (60), is found to be one in which density fluctuations no longer dominate, yet fluctuations in the phase remain; this is the regime in which the system is said to contain a quasi-condensate (59).

Under tight transverse confinement, the description of the dynamics of highly elongated three-dimensional systems effectively reduces to consideration of a one-dimensional equation, as dynamics in the transverse direction becomes restricted to the lowest mode. As this is the case in many current experiments with atom chips (5), it is crucial to fully understand and characterize such fluctuations.
The SGPE is well suited to investigations of phase-fluctuating, or quasi-condensates, as fluctuations about the mean field are implicitly retained within the order parameter, thus giving access to information on coherence, available through the evaluation of various temporal and spatial numerical auto-correlation functions.

Figure 1: Snapshots of the density and first-order correlation function during the growth of a one-dimensional 23Na quasi-condensate. The simulations were performed for an effective chemical potential μ=30Eosc, where Eosc=ℏωz is the harmonic oscillator energy, with ωz=2π×30Hz at a temperature T=100nK for a final quasi-condensate atom number N≈20000.
The effective 1D interaction strength g1d is obtained by averaging over transverse Gaussian profiles of width l⊥=√ℏ/mω⊥, where ω⊥=2π×120Hz. For the chosen somewhat artificial parameters, the phase degeneracy temperature is relatively high, slightly in excess of 100nK.
Displayed snapshots were taken at t=0.05teq, t=0.15teq, t=0.3teq and t=0.45teq, where teq denotes an approximate time for the system to reach dynamical equilibrium.

A study of the lowest order correlation functions of the gas reveals crucial information about the system. In particular, information about the phase of the system is contained in the first order correlation function, with the next higher order correlation containing crucial details about density fluctuations. Such quantities can be directly measured in interferometric ‘juggling’ experiments (22); (23); (24), and can also be easily obtained numerically within this formalism. More specifically, the renormalised first order correlation function at time t can be obtained in ‘symmetric form’ about a point z0, via the relation

We start our discussion by the basic
demonstration of the stochastic method, investigating the nonequilibrium evolution of the density and first-order correlation function during the growth of a one-dimensional quasi-condensate from a static thermal cloud. Fig. 1 shows snapshots of the time evolution of both the density (top row) and coherence (bottom row) towards their equilibrium values, for a gas in contact with a heat bath at a fixed temperature.
As the atomic system is cooled and heads towards the new equilibrium, the spatial extent of the correlation function increases as phase coherence is established, with the final shape attained dictated by the temperature, for a given particle species, number and trap geometry.
In this and all subsequent figures on coherence, distances are scaled to the effective quasi-condensate spatial extent RTF(T) at equilibrium. This quantity takes account of the varying system size at different temperatures due to quasi-condensate depletion, and reduces at T=0 to the ‘usual’ Thomas-Fermi solution RTF(T=0)=√2μ/mω2z. There are two ways to obtain RTF(T), as discussed in (53): firstly, by considering the modified low-dimensional theory of Andersen et al.(26), or equivalently by using the curvature of the density profiles to obtain a corresponding quasi-condensate size.
The particular example considered here corresponds to the case where the coherence length of the system is comparable to the system size, with T≈0.9Tϕ.
The observed growth in the densities is a direct consequence of the transfer of particles from the heat bath to the system, and the noise is crucial in order to seed the growth process (40); (41).

For temperatures such that T>Tϕ, the first-order correlation function has been found experimentally to exhibit exponential decay, whereas for T<Tϕ the corresponding spatial dependence becomes Gaussian (22); (23); (24).
Thus, one possible way to quantify the degree of coherence and extract a uniformly-varying coherence length is based on a simple fitting function of the form (53), f(z)=exp{−[(z/Lcoh)+ζ(z/Lcoh)2]},
in which ζ is a parameter characterizing the crossover from exponential to Gaussian behaviour.

A comparison between the spatial coherence at the trap centre (i.e. for z0=0) for early and late times, together with the results of the fitting procedure, is shown in Fig. 2 for two temperatures, corresponding roughly to the case of a ’pure’ condensate (left images) and a gas at the crossover from ‘pure’ to ‘quasi’-condensation (right images).

Figure 2: Coherence growth at the trap centre for the temperatures T=0.2Tϕ and T=0.9Tϕ. The upper row shows the coherence evaluated at a time early in the growth, t=0.15teq, for each temperature, whereas the plots in the lower row display the coherence at the equilibrium time, t=teq.
Solid (noisy) curves correspond to numerical results obtained after suitable averaging over a few thousand runs, whereas dashed lines indicate the ‘optimum’ fit using the function f(z) given in the text. Other parameters are as in Fig. 1.

As is evident from this figure (right images), for
the higher temperature results, corresponding to the parameters of Fig. 1, the equilibrium correlation function decays more rapidly than for the lower temperature limit.
Although the spatial extent of the system differs between these two cases (RTF(0.2Tϕ)=7.4lz vs. RTF(0.9Tϕ)=6.6lz), the appropriately scaled correlation function plotted in the top row of Fig. 2 exhibit very similar behaviour, unlike the case at their respective equilibria (bottom row).
Nonetheless, the chosen fitting function seems to appropriately capture the main features; importantly, this enables us to study in a universal manner the growth of coherence as the system is equilibrating, being driven by the heat bath.

Numerous experiments (61); (62); (63); (64) have been performed to study the relation between condensate growth and evolution of coherence, and a detailed comparison to such experiments is pending
(but see also (66); (65) and references therein). For our present illustrative purposes, we simply highlight that the growth of coherence appears to exhibit similar behaviour to typical quasi-condensate number growth curves, although these generally follow such curves with an additional small lagging time; most notably, once the norm has reached its steady-state value, the coherence is still growing, albeit at a very slow rate.

With this means of quantifying the degree of coherence across a sample, it is then interesting to assess how Lcoh varies with time during the cooling process.
The analysis of these results for the two different temperatures in the context of this somewhat idealised scenario is plotted in Fig. 3.
As these simulations are based on a fixed value of μ, the final atom numbers are slightly increased in the high temperature case (top right image); note also the significantly increased relaxation time and observation of the initial spontaneous initiation process prior to the domination of the bosonic stimulation effect (top left).

Figure 3:
Quasi-condensate atom number growth curves for a system equilibrating in contact with a prescribed heat bath (top row), versus corresponding temporal evolution of the coherence length Lcoh, obtained to approximately 10% accuracy from the fitting function defined in the text (bottom row). Note the differing y-axis scales between the coherence growth for two temperatures, with results scaled to the temperature dependent Thomas-Fermi radius RTF(T).

Figure 4: Evolution of the coherence (from left to right within each image) during the growth to equilibrium with the heat bath for a temperature T=0.9Tϕ. The correlation function is evaluated at z0=0.7RTF(T) (lower left plot) and z0=1.4RTF(T) (lower right plot), with z0 indicated by the circles on the equilibrium density profile for the chosen temperature (top plot).

Dependence of Equilibrium Coherence on Position

We also investigate the position dependence of the correlation function at equilibrium, by comparing the evolution of g(1) at different positions z0, for a fixed temperature corresponding to the high temperature case presented in the above results. The results of this investigation are shown in Fig. 4 for z=0.7RTF and z=1.4RTF for various times during the growth to equilibrium.
The final profile attained at this temperature, together with the Thomas-Fermi profile for the same number of atoms, is also shown, with the positions considered indicated.
The coherence within the Thomas-Fermi radius decreases uniformly but at a relatively slow rate as soon as one probes off-centred regions, with the effective coherence length decreasing rapidly as soon as one moves beyond the Thomas-Fermi radius (right images).

Potential ’Identification’ of a Quasi-condensate

An appealing feature of the stochastic approach that is often not fully appreciated, is that it actually generates total density profiles without needing to resort to a special mean-field treatment of the condensate mode. This means that, as in the experiments, one does not have an automatic ’visualisation’ of the condensate in the system. So, in order to analyse these results, one may consider applying bimodal fits precisely as done in the experiments. However, an alternative definition has been proposed (67), in which the quasi-condensate density is identified as
nQC(z)=√⟨|Φ(z)|2⟩2−⟨|Φ(z)|4⟩.
This approach has been implemented in (53),
with the density profile contribution for the thermal cloud over low-lying modes determined self-consistently via the relation
nT(z)=⟨|Φ(z)|2⟩−nQC(z),
within the region z<RTF(T), and given simply by nT(z)≡n(z) outside of the Thomas-Fermi radius.
The identification of these quantities enables spatial plots of the temperature dependent quasi-condensate density profiles to be obtained in an ab initio manner, as shown for different temperatures in Fig. 5.
Such a definition was used more recently to distinguish between a condensate and a quasi-condensate (56); (48)

Figure 5: Quasi-condensate (solid) and thermal (dashed) density profiles obtained from the SGPE by manipulation of density-density correlations. Plotted images correspond to the following temperatures T/Tϕ= 0.1, 0.9, 1.8, and 2.9 (top left to bottom right).

Having discussed some basic elements of finite temperature trapped, quasi-one-dimensional Bose gases via the SGPE, we now discuss certain limiting cases of the above theory which are commonly used in the literature, and which
shed more light on the benefits and versatility of the SGPE.

iii.2 Reduction to Alternative Theories

We have already argued that the two important features of the SGPE are that it includes the correct damping due to its dynamical coupling to the heat bath, and that it extends beyond mean field effects, by including a stochastic noise term. Removing either (or both) of these elements leads to simpler theories (12), which can nonetheless be useful in certain regimes, as discussed below.

Removing Contact with the Heat Bath

In first instance, we can consider what happens when the coupling to the heat bath is abruptly (and perhaps somewhat artificially) removed. We recall that
upon beginning a simulation with the SGPE, the effect of the iR(r;t) term in (1) is to pump particles from the heat bath, and populate the low lying modes to their equilibrium occupations, in accordance with the classical equilibrium for a given set of trap frequencies, chemical potential and temperature. On reaching a dynamical equilibrium, the scattering rate ∝iR(r;t) becomes on average equal to zero, and there is no net particle exchange with the heat bath, although this continuous exchange can cause additional damping. However, any externally imposed perturbations of the system about this set of system parameters, would typically lead to the introduction or loss of particles, upon further evolution with (1).

In general, the SGPE is coupled to a quantum Boltzmann equation for the population of the higher-lying modes, and these two equations should be solved self-consistently. If this is not done, then modifying for example the trap geometry will lead to a change in the norm of the SGPE. If we are only considering the evolution of the low-lying modes, in order to keep the atom number within our system fixed, one could consider ‘fixing’ the norm in the numerics by removing the coupling to the heat bath. A better alternative might perhaps be to change the effective chemical potential in the SGPE to the final value corresponding to the same atom number in the new trap, although the observed dynamics will typically depend on how quickly the chemical potential is adjusted to its final value.

One procedure is thus to propagate initially the full stochastic solution in order to obtain the desired equilibrium state, before switching off the noise and damping terms of (1), and subsequently propagating the noisy initial condition via the usual GPE which enforces particle number conservation. This corresponds physically to the removal of contact with the heat bath upon reaching equilibrium.

As a method, this procedure has strong parallels with the Truncated Wigner method, and a comparative study of the two is under way currently (68).
The main difference arises in the fact that in the usual truncated Wigner implementations (see, e.g. (36); (69); (70) or (12) for a more extensive list), only quantum (as opposed to thermal) noise is typically added in the initial conditions (but see also (37); (71)), whereas the SGPE in general contains both quantum and thermal fluctuations (40). Although existing numerical implementations of the SGPE focus on the thermal effects (due to the ‘classical’ approximation mentioned earlier), the important advantage of this approach is that
there is an explicit fluctuation-dissipation relation which guarantees convergence to the correct classical equilibrium at any pre-determined temperature; in other words, the initial thermal state for the unperturbed system can be correctly assigned.
The validity of conventional truncated Wigner implementations is on the other hand limited to very low temperatures (37) and occupation numbers for which quantum fluctuations have a dominant role - see (68) for a more detailed discussion.

Figure 6: Suppression of shockwaves induced upon the introduction of a deep dimple trap to a one-dimensional ultracold Bose gas of around 3500 23Na atoms, as manifested in the temporal evolution (in units of the inverse trap frequency) of the corresponding damped central density oscillations obtained after numerical averaging, and scaled to the peak density for the given system. (Details of other parameters used in these simulations can be found in (52).) Contact with the heat bath was removed prior to the introduction of the dimple trap to maintain number conservation.

This finite temperature approach based on the removal of the heat bath once the correct equilibrium has been obtained via the SGPE
was used to probe the dynamics of an initially equilibrated thermal stationary state following the addition of a deep Gaussian dimple trap (52).
A large and sudden perturbation is found to lead to the appearance of shock waves in the system, and in the averaged profiles typically studied, these manifest themselves in terms of an initial increase in the central density, followed by periodic damped decrease-increase cycles, as shown in Fig.6; this can be interpreted as direct evidence of the temperature-induced suppression of the initial shockwave amplitude induced by the addition of the new trap.
It should be noted here that the introduction of the correct thermal noise in the initial conditions is crucial to the generation of the damping which would not be present in simulations of the ordinary GPE starting from a stationary condensed state.

A Microscopically-based Dissipative GPE

An even simpler limiting case of the SGPE is based on
removing only the noise term of (1); then, the SGPE reduces to a ‘dissipative GPE’, as noted for example in (54); (72). The form of the damping is given by the self energy, and as such is spatially dependent, consistent with the spatially inhomogeneous background of the harmonic trap.
Such a model has been used to study vortex and soliton dynamics leading to a number of interesting predictions. However, the origin of such models can be traced to rather crude approximations in the more advanced microscopic models, such as the SGPE, with most such publications justifying this equation phenomenologically, and actually using a constant value for the damping γ(73); (74); (75); (76).

Figure 7: Soliton dynamics within a one-dimensional gas of N≈2000023Na atoms, for the zero temperature (GPE) case, plus at temperatures of 100nK and 200nK. The top row plots show density snapshots at three times during the dynamics, t=24ω−1z, t=37.2ω−1z and, t=60ω−1z; shown in each figure are the results for all three temperature cases at that time (zero temperature - blue, dot-dashed line; 100nK - green, dashed line; 200nK - red, solid line). The spacetime plots are ordered downwards in order of increasing temperature, for the highest temperature case the soliton can be seen to decay rather rapidly.

The dissipative equation which results from the SGPE in general takes, within the ‘classical approximation’, the form

iℏ∂Φ(r;t)∂t=(1−iγ(r,T,t))[−ℏ2∇22m+Vext(r)−μ+g|Φ(r;t)|2]Φ(r;t),

(5)

where the damping term (for a system made up of a fixed total number of atoms) is given by

γ(r,T,t)=iβ4ℏΣK(r,t),

(6)

although our discussion so far has additionally ignored the time-dependence for numerical simplicity.

It is very easy to see that
the resulting evolution is no longer norm conserving, and it therefore becomes nessessary to renormalize the wavefunction as each time step to maintain a constant particle number.

For illustrative purposes, we apply this approach here to model the oscillations of an ‘ideally generated’ dark soliton within a one-dimensional condensate.
The dynamics for three identical initial soliton solutions are shown at different temperature in Fig. 7 (top images), with each image showing the position of the dark soliton for three different temperatures, corresponding to a different evolution time.
To best illustrate the different temperature-dependent decay rates, we give corresponding
space-time plots of the density along the z-axis
with such damping already observed in experiments (77).

Although consideration of the 1D case enables us to bring out the main physical features of the SGPE model, comparison to experiments requires us to consider its properties in higher-dimensional systems, and we conclude our presentation by a simple illustration in two-dimensional geometries, and some related comments.

As in one dimension, for a uniform two dimensional Bose gas, thermal fluctuations remain even at very low temperatures, preventing the onset of Bose-Einstein condensation
(78); (79). For interacting trapped two dimensional gases, however, quasi-condensation is possible for low temperatures (80), with a Berezinskii-Kosterlitz-Thouless phase transition to a superfluid state associated with the the pairing of vortices of opposite circulation beneath some critical temperature. Recently, this has spawned much interest in the physics of two-dimensional trapped ultracold systems, particularly in assessing the nature of the phase transition at the critical point (6); (81); (82); (83).

iv.1 Spontaneous Processes - Vortex Nucleation

Phase transitions are strongly associated with such topological defects, the appearance of which are predicted by the Kibble-Zurek mechanism
(84); (85).
An appealing feature of the stochastic approach is that in including fluctuations about the mean field, excitations in the form of spontaneous vortices can be seen during the growth process, as investigated recently using an alternative formulation of the SGPE (57).
Our discussion so far has been restricted to analysing predictions of the SGPE based on numerical averaging over different initial noisy conditions. However, the monitoring of the evolution within a single run can also provide important information on the system parameters, which may otherwise be lost through averaging - e.g. information related to spontaneous events, such as topological defect formation (86).

Within the numerical SGPE applied here, condensate formation from the static thermal cloud is initiated for a positive value of the chemical potential in (II.1). A rapid quench can be simulated by an instantaneous change in the chemical potential, leading to a rapid condensate growth. Within simulations, the appearance of a (random) number of vortices is common, although their details vary from run to run. The appearance of these actually becomes ‘washed out’ following the averaging process, as, due to the random nature of the nucleation of vortices, the positions at which they appear varies between different realisations.
Such studies should be contrasted to studies based on the GPE, in which a vortex is actually ‘artificially’ imprinted at a predetermined position; when modelling the ensuing dynamics, one often introduces damping by resorting to the dissipative GPE (???), with the damping term, γ(x,y), which is actually proportional to the two-dimensional self-energy at a particular temperature, typically treated as a constant.

Figure 8: Nonequilibrium density and phase snapshots for the growth of a two dimensional harmonically trapped quasi-condensate in the x-y plane, for N≈6000023Na atoms, dimensionless coupling constant g2d=√8πa3d/lz≈0.042 and trap frequencies ωx=ωy=2π×200Hz and ωz=2π×4kHz. The single run density profiles (top row) show the spontaneous appearance of vortices, as apparent also in the single run phase plots (middle row). Time increases to the right in the plots; at the centre of the single run images two vortices can be seen to come together and annihilate. The bottom row profiles are density profiles obtained following averaging over 100 single runs, with snap-shots taken at the same times as for the single runs. The results show a far smoother density profile is obtained upon averaging, though vortices are now no longer visable.

In density plots, vortices appear as dark regions, corresponding to positions at which singularities in the phase of the system can be seen. An example of the results from a typical simulation based on the SGPE is shown in Fig. 8. In the density and phase snapshots shown, two spontaneously generated vortices (top left) can be seen to come together (middle) and then annihilate (right), illustrating the spontaneous dynamics captured within simulations. The averaged density profiles included show the result of averaging over 100 runs, and the smooth density plots which result.

iv.2 Numerical Issues - Dependence on Momentum Cutoff

Figure 9: Dependence of the norm on the numerical grid spacing for the two-dimensional SGPE. The fit shows the dependence of the particle number on the high momentum cut-off, Λ, introduced by the lattice for fixed chemical potential, temperature and trap frequencies. This is inversely related to the grid spacing Δz, against which the results are plotted. The divergence is found to be ∝log(Λ) in the two dimensional case, as shown by the line indicating the fit to a function ∝log(π/Δz).

The numerical means of evolving the system discussed so far make use of the
macroscopic occupation of the low lying modes, assuming these states to be well
described by a classical field. In two dimensions, known divergences arise,
associated with large momenta, or equivalently small grid spacings, as the
continuum limit is approached on the lattice. This is due to the high momentum cutoff, Λ, related to the numerical grid spacing, Δz, through the relation Λ=π/Δz.

We present here preliminary results of an investigation into this dependence for the SGPE as a ‘caution’ to the reader. Here we make no adjustments to parameters to compensate for grid dependence and find the expected logarithmic divergence in the momentum cut-off, with the system norm (or number of particles in the classical region) increasing with decreasing spatial discretization Δz, as shown in Fig. 9.

In principle such divergences can be removed from the problem through the appropriate renormalization of physical quantities from their original (‘bare’) values. For example in homogeneous systems, the lattice effects on temperature or density can be taken into account in an a posteriori manner, as described in (87). An alternative method is to introduce counter-terms to the system Hamiltonian, which has been successfully applied to Langevin equations previously in, for example, (89); (88); (90); (91).

The stochastic Gross-Pitaevskii is emerging as a key approach for the equilibrium and dynamical properties of degenerate Bose gases, and is particularly beneficial for describing regimes of large phase and density fluctuations, such as near (or at) the phase transition, or for weakly-interacting low-dimensional Bose gases. Single run results contain important (at least qualitative) information on spontaneous processes, such as defect formation characteristics, whereas averaged profiles are most suitable for determination of correlation functions, quasi-condensate fractions and smooth density profiles. Although these stochastic simulations have obvious advantages, their existing numerical implementation have not yet explored such theories to their full potential, with current models typically assuming a static (near-equilibrium) thermal cloud, and additionally discarding quantum fluctuations in simulations. Finally, the important issue of the cut-off dependence, already highlighted by other researchers should not be overlooked, particularly for dimensions higher than one, and various techniques are available for renormalising the observed profiles to their actual values, thereby eliminating (at least to leading order) such dependences. We believe that use of the stochastic Gross-Pitaevskii equation, which is also coupled to a quantum Boltzmann equation for the self-consistent dynamical treatment of the higher-lying thermal modes, should be able to describe essentially all main features of weakly-interacting Bose gases currently pursued experimentally, and therefore look forward to the exciting times which lie ahead.

Acknowledgments

We are grateful to Henk Stoof for early collaboration and discussion on these ideas and to the UK EPSRC for funding.