Truncated Wigner Method for Bose Gases

# Truncated Wigner Method for Bose Gases

## Abstract

We discuss stochastic phase-space methods within the truncated Wigner approximation and show explicitly that they can be used to solve non-equilibrium dynamics of bosonic atoms in one-dimensional traps. We consider systems both with and without an optical lattice, and address different approximations in the stochastic synthesization of quantum statistical correlations of the initial atomic field. We also present a numerically efficient projection method for analyzing correlation functions of the simulation results, and demonstrate physical examples of non-equilibrium quantum dynamics of solitons and atom number squeezing in optical lattices.

## I Introduction

In stochastic phase-space methods based on sampling classical probability distributions the common approach is to unravel the evolution dynamics into stochastic trajectories each of which obey the classical mean-field dynamics, with or without additional dissipative coupling terms to environment. Each trajectory is a representative of a probabilistic initial state distribution that is numerically generated by a Monte Carlo type of sampling. The probability distribution is selected in order to synthesize as closely as necessary, e.g., the thermal distribution or quantum statistical correlations of the initial state. The phase-space representation that most accurately reproduces classical mean-field dynamics is the Wigner representation because of the ‘correct amount’ of quantum noise in the initial state Gardiner and Zoller (1999); Sinatra et al. (2002). It has become common to call mean-field dynamical simulations together with sampling of quantum noise as the Truncated Wigner Approximation (TWA) Drummond and Hardman (1993); Steel et al. (1998); Sinatra et al. (2002); Isella and Ruostekoski (2005, 2006); Blakie et al. (2008); Polkovnikov (2010); Martin and Ruostekoski (2010).

In this Chapter, we consider the unitary evolution in the TWA formalism in ultra-cold atom systems without additional dissipative coupling to an environment (for TWA employed in an open system, see for instance Ref. Javanainen and Ruostekoski (2011)) that is particularly suitable for analyzing dissipative non-equilibrium quantum dynamics in 1d systems. The TWA approach is able to represent systems with a large number of degrees of freedom using a stochastic representation of the atomic field operator. In this method, dissipative dynamics emerge from a microscopic treatment of the unitary quantum evolution, due to energy dissipation within the large phase-space without any additional explicit damping terms in the Hamiltonian. Quantum and thermal fluctuations of the atoms are included in the initial state, and the resulting quantum statistical correlations of the initial state may be accurately synthesized for different quantum states in the Wigner representation.

Atomic systems with enhanced quantum fluctuations that may be modeled with TWA can, for instance, be prepared in tightly-confined cigar-shaped atom traps, where the strong transverse confinement suppresses density fluctuations along the radial direction of the trap (see, e.g., Ref. Tolra et al. (2004); Kinoshita et al. (2006); Fertig et al. (2005); Mun et al. (2007)). Quantum effects may be further strengthened by reducing the kinetic energy of the atoms by means of applying an optical lattice potential along the axial direction Fertig et al. (2005); Mun et al. (2007).

In this Chapter we briefly summarise the TWA method, before addressing different approximations in synthesization of quantum and thermal noise in the initial state. We start with a simple uniform system and phonon excitations within the Bogoliubov approximation. These are extended to non-uniform systems and situations where the back-action of the excited-state correlations on the ground-state atoms is included in a self-consistent manner. A particular problem of analyzing non-equilibrium quantum dynamics in TWA, related to symmetric operator-ordering of the Wigner distributions is addressed by providing a numerically practical projection technique.

Applications of the approaches presented here include the studies of dissipative non-equilibrium systems both with and without an optical lattice, such as the fragmentation of a BEC by ramping on an optical lattice Isella and Ruostekoski (2005, 2006); Gross et al. (2011), dissipative atom transport Ruostekoski and Isella (2005), dynamically unstable lattice dynamics Shrestha et al. (2009), and dark solitons Martin and Ruostekoski (2010, 2010), some of which are discussed in Sec. IV.

## Ii Methodology

We describe non-equilibrium quantum dynamics of bosonic atoms by the evolution that can be considered with a good approximation to be unitary. Quantum and thermal noise enter the equations of motion for an atomic field only through the stochastic initial state. The dynamics is governed by the Gross-Pitaevskii Equation (GPE),

 iℏ∂∂tΦ=(−ℏ2∇22m+Vext+g|Φ|2)Φ, (1)

where in 1d the interaction strength , the -wave scattering length , and is the trapping potential Olshanii (1998). We concentrate on 1d dynamics in tightly-confined atom traps. In higher dimensions, the unitary evolution may typically be replaced by a model with an explicit low-momentum cut-off Sinatra et al. (2002). Other realizations of TWA have, e.g., treated the system and environment separately with a coupling between the two, that results in an explicit dynamical noise term in each time step Drummond and Hardman (1993), or by including a continuous quantum measurement process Javanainen and Ruostekoski (2011). Here, depending on the physical problem, we could also include additional terms in Eq. (1), e.g., atom losses via collisions and spontaneous emission that would generate also dynamical noise terms for each time step. In Eq. (1) we have explicitly included the atom number in the nonlinear coefficient, so we use the normalization , where denotes the number of modes in the initial state and the total number of atoms.

Unlike in the usual GPE, here should be considered as a stochastic phase-space representation of the full field operator describing the time evolution of the ensemble of Wigner distributed wavefunctions. The time evolution is unraveled into stochastic trajectories, where the initial state of each realization for the classical field is stochastically sampled in order to synthesize the quantum statistical correlation functions for the initial state.

Since for the unitary evolution all the noise is incorporated in the initial state, it is especially important that the quantum mechanical correlation functions for the initial state of the atomic field operator are synthesized as accurately as practical for each particular physical problem. Here we follow our basic formalism of Refs. Isella and Ruostekoski (2005, 2006); Shrestha et al. (2009); Martin and Ruostekoski (2010, 2010); Gross et al. (2011).

### ii.1 Initial State Generation in the TWA

#### Uniform system

In the case of a weakly interacting bosonic gas in a uniform space at the simplest approach to model quantum fluctuations of the atoms, if we are not interested in the conservation of the total atom number, is the Bogoliubov approximation. In the Bogoliubov theory we calculate the linearized fluctuations of the ground state (or a stationary GPE solution) in which case the back-action of the excited-state atoms on the ground state is ignored Pethick and Smith (2002). We write the decomposition

 ^Ψ(x)=ϕ0(x)^b0+^ψ′(x), (2)

where the total number of ground state atoms . The fluctuation part for the excited states can be written in terms of quasi-particle operators and as

 ^ψ′(x,t)=1√L∑q≠0(uq^bqeiqx−v∗q^b†qe−iqx). (3)

The normal mode frequencies and the quasi-particle amplitudes and can be solved straightforwardly Pethick and Smith (2002) and the number of excited-state atoms in the Bogoliubov theory is

 N′=∑q(|uq|2+|vq|2)nBE(ϵq)+∑q|vq|2, (4)

with denoting the ideal Bose-Einstein distribution. At we have .

In order to construct the initial state for the atoms in the TWA evolution we replace the quantum field operators by the classical fields by using complex stochastic variables in the place of the quantum operators in Eq. (3).

In the Bogoliubov theory, the operators form a set of ideal harmonic oscillators and at , (for ) are obtained by sampling the corresponding Wigner distribution function Gardiner and Zoller (1999)

 W(βq,β∗q)=2πexp[−2|βq|2]. (5)

In this case each unoccupied excitation mode is uncorrelated with Gaussian-distributed noise. The expectation value specifies the width of the distribution and represents vacuum noise, resulting from the symmetric ordering of the expectation values in the Wigner representation. The noise is distributed in space according to the plane waves, with a constant density. In the absence of any correlations between the modes, the vacuum noise in the uniform space could be replaced by uncorrelated Gaussian noise on evenly-spaced numerical grid points. However, if we do not want to allow the total atom number to fluctuate between different trajectories (conserved atom number), the simplest modification to the Bogoliubov expansion is to fix the total atom number in each stochastic realization. This introduces long-wavelength correlations between the ground-state mode and the excited-state phonon modes, so that already in this simple example there exist non-trivial spatial noise correlations Martin and Ruostekoski (2010, 2010). For each stochastic realization the number of excited-state atoms satisfies

 N′s=∑q(|uq|2+|vq|2)(β∗qβq−12)+∑q|vq|2. (6)

where fluctuates in each realization with the ensemble average at . In Eq. (6) we have transformed the symmetric ordering of the Wigner representation to quantum expectation values of normally-ordered operators by subtracting from each mode. The ground-state atom number is then obtained from the fixed total atom number , so that in each stochastic realization and we set . The ensemble average of the ground-state population is obtained from .

At we replace Eq. (5) by Gardiner and Zoller (1999)

 W(βq,β∗q)=2πtanh(ϵq/2kBT)exp[−2|βq|2tanh(ϵq/2kBT)]. (7)

The Wigner function is Gaussian-distributed with width . The formula introduces thermal populations of each quasi-particle mode and generates more complex spatial noise correlations. After the noise generation, the initial state for stochastic evolution at time may be written as

 Φ(x)=ϕ0(x)β0+1√L∑q≠0(uqβqeiqx−v∗qβq∗e−iqx). (8)

Here is a stochastic representation of the full field operator for the atoms.

#### Trapped Gases

Placing atoms in a non-uniform potential results in a spatially-varying initial noise distribution even at . For a combined harmonic trap and optical lattice we write the external potential as , where is the lattice photon recoil energy and is the lattice period. The Bogoliubov equations now become spatially dependent and need to be solved numerically Pethick and Smith (2002). In Eq. (3) we replace and , where the index refers to the mode number. In the lowest order approximation the quasi-particle mode functions and are obtained in the Bogoliubov theory. In several cases of interest where the multi-mode structure of the excitations become important, the Bogoliubov approximation is insufficient due to the large contribution of the quadratic fluctuation terms. One consequently needs to use a higher-order theory in which case the ground-state and the excited-state populations are solved self-consistently. One such candidate is the gapless Hartree-Fock-Bogoliubov (HFB) formalism: this is similar to the usual HFB approach, but constructed in such a manner that there is no gap in its excitation spectrum at zero momentum. The coupled equations for the ground state and excitations thus take the general form Hutchinson et al. (1998); Proukakis et al. (1998)

 (^L−U\scriptsize c¯Nc|ϕ0|2)ϕ0=0 (9) ^Luj−U\scriptsize c¯Ncϕ20vj=ϵjuj, ^Lvj−U\scriptsize c¯Ncϕ∗20uj=−ϵjvj. (10)

where and () are restricted to the subspace orthogonal to . In order to express these in a form amenable to stochastic simulations, we have used the notation . Here

 ^L≡−ℏ22m∂2∂x2+V(x)+2U\scriptsize c¯Nc|ϕ0|2+2U\scriptsize e% n′(x)−μ, (11)

and is the chemical potential. This general notation contains numerous theories as sub-cases for appropriate choices of the interaction strength of a condensed atom with another condensed atom () or a thermal atom (): the Bogoliubov approximation is obtained by setting , in Eq. (II.1.2). Setting and yields the gapless HFB theory (the G1 version in Ref. Proukakis et al. (1998)); here is the depleted density, and the anomalous pair correlation, both of which introduce back-action of the excitations on the ground-state. Hence, Eqs. (9) and (II.1.2) must be solved iteratively until the solutions converge. In the non-uniform case the number of excited-state atoms is given by

 N′=∫dx∑j[(|uj(x)|2+|vj(x)|2)nBE(ϵj)+|vj(x)|2], (12)

and the total atom number may be fixed in each realization as in the uniform case Martin and Ruostekoski (2010, 2010). The gapless HFB theory was introduced as a stochastic sampling technique for TWA simulations in Ref. Gross et al. (2011) to model reduced atom number fluctuations, fragmentation and spin-squeezing in optical lattice systems.

#### Quasicondensate Description

In tightly-confined 1d traps, the phase fluctuations may be enhanced compared to those obtained using the standard Bogoliubov theory Kheruntsyan et al. (2003). A more accurate description can be calculated using quasi-condensate formalism Mora and Castin (2003) that can be particularly important, e.g., to phase kinks Martin and Ruostekoski (2010, 2010). In the quasi-condensate description we write the field operator as

 ^Ψ(x)=√n0(x)+δ^n(x)exp[i^θ(x)]. (13)

The density and phase operators are written as (for )

 ^θ(x) = −i2√n0(x)∑j(θj(x)^bj−θ∗j(x)^b†j), (14) δ^n(x) = √n0(x)∑j(δnj(x)^bj+δn∗j(x)^b†j), (15)

where and are given in terms of the solutions to the Bogoliubov equations (see the previous section). This results in a stochastic Wigner representation of phase and density operators. The stochastic initial state at for the time evolution then reads Martin and Ruostekoski (2010, 2010)

 Φ(x)=√n0,W(x)+δnW(x)exp(iθW(x)), (16)

where the ground-state density .

#### Relaxation

We may also consider an ideal, non-interacting BEC as an initial state for the TWA simulations, but before the actual time evolution, we can continuously turn up the nonlinear interactions between the atoms. If the process is slow enough and relaxes to the ground state, we may be able to produce the stochastic initial state of the interacting system. Although this may simplify the calculations, in practice the technique in a closed system does not necessarily converge to the correct interacting state Isella and Ruostekoski (2006). More complex models with open systems, kinetic equations and time-dependent noise can help the relaxation process at finite temperatures Cockburn et al. (2010).

### ii.2 Wigner Representation and Symmetric Ordering

The Wigner distribution returns symmetrically-ordered expectation values of any stochastic representations of quantum operators. In particular, the expectation values of the full multi-mode Wigner fields in the TWA simulations of the time-dynamics are symmetrically ordered with respect to every mode. In general, this can significantly complicate the analysis of the numerical results when quantum fluctuations are important Isella and Ruostekoski (2006). A numerically practical transformation of the symmetrically-ordered expectation values to the normally-ordered expectation values of physical observables can be done using projection techniques (see Refs. Isella and Ruostekoski (2005, 2006); Gross et al. (2011)). In the presence of an optical lattice, a natural approach is to project the stochastic field on to the several lowest mode functions of the individual lattice sites. We denote the annihilation operator for the atoms in the th vibrational mode of the site as . We write the corresponding stochastic amplitude as which can numerically be obtained from

 ai,j(t)=∫ithwelldx[φi,j(x,t)]∗Φ(x,t), (17)

where is the stochastic field and is the th vibrational mode function of the site . The integration is performed over the th site and the normally ordered quantum expectation values for the site populations reads

 ⟨^ni⟩=∑j⟨^a†i,j^ai,j⟩=∑j[⟨a∗i,jai,j⟩e−1/2], (18)

Fluctuations are calculated using analogous transformations. For the on-site fluctuations of the atom number in the th site we obtain

 (Δni)2 = ⟨^n2i⟩−⟨^ni⟩2 (19) = ∑j,k[⟨a∗i,jai,ja∗i,kai,k⟩e−⟨a∗i,jai,j⟩e⟨a∗i,kai,k⟩e−δjk/4].

Similarly, the relative atom number fluctuations between the sites and are obtained from

 [Δ(^np−^nq)]2 = ∑i,k[⟨(a∗p,iap,i−a∗q,iaq,i)(a∗p,kap,k−a∗q,kaq,k)⟩e (20) − ⟨a∗p,iap,i−a∗q,iaq,i⟩e⟨a∗p,kap,k−a∗q,kaq,k⟩e−δik/2].

Alternatively, we could have, for instance, written

 ⟨^nj⟩=∫jdx⟨^Ψ†(x)^Ψ(x)⟩ = ∫jdx⟨Φ∗(x)Φ(x)⟩e (21) − 12∫jdx∑i(|ui(x)|2−|vi(x)|2).

Calculation of then, however, results in double integrals over the sites that can be computationally slow when performed over a large number of realizations.

### ii.3 Numerical Implementation

In the numerical implementation the initial state fluctuations are solved by first finding a (stable) stationary state or a local energetic minimum in GPE. If the system is assumed to be initially in thermal equilibrium, we find the corresponding ground state by imaginary time evolution of GPE, e.g., by using the nonlinear split-step Fourier methods Javanainen and Ruostekoski (2006). Using the ground state solution we then diagonalize the linearized equations for the quasi-particle excitations in order to obtain the eigenfunctions and the corresponding eigenenergies. In the case of self-consistent HFB method, the excitations and the ground state are solved iteratively until the solutions converge Hutchinson et al. (1998).

The time evolution of the ensemble of Wigner distributed wavefunctions is unraveled into stochastic trajectories, where the initial state of each realization for the classical stochastic field is generated using the quasiparticle mode functions and amplitudes Isella and Ruostekoski (2006). The complex, Gaussian-distributed stochastic mode amplitudes are sampled using the Box-Muller algorithm W. H. Press and Flannery (2002). During the time evolution we simulate some physical process that describes the changing of the equilibrium configuration, e.g., displacement of atoms from the trap center Ruostekoski and Isella (2005) or turning up of the optical lattice potential Isella and Ruostekoski (2005, 2006). The integration of the time dynamics is also performed using the nonlinear split-step methods Javanainen and Ruostekoski (2006), typically on a spatial grid of a few thousand grid points. In several cases the sufficient convergence is obtained after 600-1000 realizations.

In order to transform the symmetrically-ordered expectation values of the Wigner representation into normally-ordered expectation values, we numerically introduce an orthonormal basis, e.g., in each lattice site. The stochastic field is then at different times projected onto this basis and the desired expectation values are evaluated using the transformations for each projected mode function, as described in the previous section.

## Iii Validity Issues

In 2d and 3d the TWA can have implementation problems. Firstly, the atom cloud can heat during time evolution due to rapid nonlinear dynamics between the vacuum modes Sinatra et al. (2002). Secondly, physical observables can diverge as a function of the number of modes (or equivalently energy cut-off or grid spacing). Importantly for the present discussions, 1d systems are more robust to these effects. In particular, TWA has been successful in describing superfluid dynamics in the presence of considerable quantum fluctuations in 1d systems, even though it is clearly insufficient, e.g., in a Mott-insulator regime and at very low atom numbers. For instance, TWA simulations Ruostekoski and Isella (2005) were qualitatively able to produce the experimentally observed damping rate of center-of-mass oscillations of bosonic atomic cloud in a very shallow, strongly confined 1D optical lattice, corresponding to the dissipative atom transport experiments of Ref. Fertig et al. (2005) in which case atom numbers approximately down to 70-80 were used in an elongated trap of a very large phase space.

The accuracy of the initial state noise generation can be a crucial limitation, especially in simulations involving very short time dynamics. The spatial distribution of phonon excitations in trapped systems can result in very rapid noise variation where, e.g., phase fluctuations dominate near the edges of the atom cloud Isella and Ruostekoski (2006). For dark soliton dynamics, the differences in the soliton trajectories between the cases in which the noise was generated within the quasi-condensate description and in the Bogoliubov theory are notable Martin and Ruostekoski (2010), indicating that the soliton imprinting process and dynamics are sensitive to enhanced phase fluctuations of the quasi-condensate description Mora and Castin (2003). Evaluating phonon modes in the linearized Bogoliubov approximation may also become inaccurate, compared to self-consistent HFB methods even at , as demonstrated in the case coupled condensates in a few-site lattice system Gross et al. (2011).

Stochastic phase-space methods based on sampling classical probability distributions are necessarily approximate, unless the problem is reformulated, e.g., by doubling the phase space and considering and as independent fields. Such positive-P Drummond and Gardiner (1980); Gardiner and Zoller (1999); Carusotto et al. (2001); Drummond et al. (2004) or positive-Wigner L. I. Plimak et al. (2001) methods can in principle provide exact solutions but frequently run into numerical problems due to rapidly growing sampling errors.

## Iv Applications

### iv.1 Dark Solitons

Dark solitons have been actively studied in BECs Burger et al. (1999); Denschlag et al. (); Anderson et al. (2001); Dutton et al. (2001); Weller et al. (2008); Stellmer et al. (2008), and in nonlinear optics Kivshar and Luther-Davies (1998). Although there exist numerous studies of classical solitons, the quantum properties of dark solitons are much less known. Numerical TWA simulations are suitable for the studies of the creation and non-equilibrium quantum dynamics of solitons in 1d traps. We consider the experimental imprinting method Burger et al. (1999); Denschlag et al. (), where a soliton is generated by applying a ‘light-sheet potential’, of value to half of the atom cloud, for time , so that in the corresponding classical case the light sheet imprints a phase jump of at , preparing a dark soliton. Classically the imprinted soliton oscillates in a harmonic trap at the frequency Busch and Anglin (2000) with the initial velocity , depending on and the speed of sound . The soliton is stationary (dark) for , with zero density at the kink. Other phase jumps produce moving (grey) solitons, with non-vanishing densities at the phase kink.

In TWA simulations we generate the initial state using the quasi-condensate formalism and vary the ground-state depletion . At we keep the nonlinearity fixed, but adjust the ratio . This is tantamount to varying the effective interaction strength Kheruntsyan et al. (2003). We can also study the effects of thermal depletion by varying .

In the presence of noise, soliton trajectories in the TWA fluctuate between different realizations due to quantum and thermal fluctuations Dziarmaga et al. (2003); Mishmash and Carr (2009); Martin and Ruostekoski (2010, 2010); Cockburn et al. (2010). Individual stochastic realizations of in a harmonic trap represent possible experimental observations of single runs Isella and Ruostekoski (2006); Martin and Ruostekoski (2010, 2010). In the TWA we can ensemble average hundreds of stochastic realizations in order to obtain quantum statistical correlations of the soliton dynamics. We numerically track the position of the kink at different times in individual realizations and calculate the quantum mechanical expectation values for the soliton position and its uncertainty Martin and Ruostekoski (2010, 2010). Our results are summarized in Fig. 1.

### iv.2 Atom Number Squeezing

In this section we consider an example of a TWA calculation of spin and relative atom number squeezing due to turning up of an optical lattice. Unlike in Refs. Isella and Ruostekoski (2005, 2006) where a BEC fragmentation in TWA was investigated by a lattice with a large number of small sites, we simulate a six-site system, analogous to the recent experimental observations of spin and relative atom number squeezing Esteve et al. (2008); Gross et al. (2011) as well as reduced on-site atom number fluctuations and long-range correlations Gross et al. (2011) between coupled condensates. Bose-condensed Rb atoms are confined to a cigar-shaped optical dipole trap where an optical lattice is applied along the axial direction. The lattice potential is slowly turned up from to . Due to large individual lattice sites the multi-mode structure of the fluctuations is important, and the atom number fluctuations are evaluated by using the projection technique on to several modes in each site, as explained in the previous section. The Bogoliubov approximation is not accurate due to phonon-phonon interactions, indicating a significant contribution of higher-order terms to atom number fluctuations even at , and the initial state is calculated using the HFB method Gross et al. (2011). The TWA simulation results demonstrated a qualitative agreement with the experimental observations, although the experiment was not performed in a tightly-confined 1d trap with completely suppressed radial density oscillations. The spatially non-uniform distribution of quantum and thermal fluctuations is clearly seen in Fig. 2. The lowest HFB modes dominantly occupy the outer regions of the atom cloud with significantly enhanced atom number and phase fluctuations in those sites. Such fluctuations could not be represented, e.g., by a uniform stochastic noise sampling.

## V Comparisons

The accuracy of the initial state noise generation can contribute to the simulation results. For dark soliton dynamics, the differences in the soliton trajectories between the cases in which the noise was generated within the quasi-condensate description and in the Bogoliubov theory are notable Martin and Ruostekoski (2010). Some examples are illustrated in Fig. 3. As previously noted, also evaluating phonon modes in the linearized Bogoliubov approximation may become inaccurate, compared to self-consistent HFB methods.

## Acknowledgments

We acknowledge financial support from EPSRC and Leverhulme Trust.

### References

1. C. Gardiner and P. Zoller, Quantum Noise (Springer, Berlin, 1999).
2. A. Sinatra, C. Lobo,  and Y. Castin, Journal of Physics B: Atomic, Molecular and Optical Physics, 35, 3599 (2002).
3. P. D. Drummond and A. D. Hardman, EPL (Europhysics Letters), 21, 279 (1993).
4. M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls,  and R. Graham, Phys. Rev. A, 58, 4824 (1998).
5. L. Isella and J. Ruostekoski, Phys. Rev. A, 72, 011601 (2005).
6. L. Isella and J. Ruostekoski, Phys. Rev. A, 74, 063625 (2006).
7. P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh,  and C. W. Gardiner, Advances in Physics, 57, 363 (2008).
8. A. Polkovnikov, Annals of Physics, 325, 1790 (2010), ISSN 0003-4916.
9. A. D. Martin and J. Ruostekoski, Phys. Rev. Lett., 104, 194102 (2010a).
10. J. Javanainen and J. Ruostekoski, arXiv:1104.0820 (2011).
11. B. L. Tolra, K. M. O’Hara, J. H. Huckans, W. D. Phillips, S. L. Rolston,  and J. V. Porto, Phys. Rev. Lett., 92, 190401 (2004).
12. T. Kinoshita, T. Wenger,  and D. S. Weiss, Nature, 440, 900 (2006).
13. C. D. Fertig, K. M. O’Hara, J. H. Huckans, S. L. Rolston, W. D. Phillips,  and J. V. Porto, Phys. Rev. Lett., 94, 120403 (2005).
14. J. Mun, P. Medley, G. K. Campbell, L. G. Marcassa, D. E. Pritchard,  and W. Ketterle, Phys. Rev. Lett., 99, 150604 (2007).
15. C. Gross, J. Estève, M. K. Oberthaler, A. D. Martin,  and J. Ruostekoski, Phys. Rev. A, 84, 011609 (2011).
16. J. Ruostekoski and L. Isella, Phys. Rev. Lett., 95, 110403 (2005).
17. U. Shrestha, J. Javanainen,  and J. Ruostekoski, Phys. Rev. A, 79, 043617 (2009).
18. A. D. Martin and J. Ruostekoski, New Journal of Physics, 12, 055018 (2010b).
19. M. Olshanii, Phys. Rev. Lett., 81, 938 (1998).
20. C. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases (Cambridge University Press, Cambridge, 2002).
21. D. A. W. Hutchinson, R. J. Dodd,  and K. Burnett, Phys. Rev. Lett., 81, 2198 (1998).
22. N. P. Proukakis, S. A. Morgan, S. Choi,  and K. Burnett, Phys. Rev. A, 58, 2435 (1998).
23. K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond,  and G. V. Shlyapnikov, Phys. Rev. Lett., 91, 040403 (2003).
24. C. Mora and Y. Castin, Phys. Rev. A, 67, 053615 (2003).
25. S. P. Cockburn, H. E. Nistazakis, T. P. Horikis, P. G. Kevrekidis, N. P. Proukakis,  and D. J. Frantzeskakis, Phys. Rev. Lett., 104, 174101 (2010).
26. J. Javanainen and J. Ruostekoski, Journal of Physics A: Mathematical and General, 39, L179 (2006).
27. W. T. V. W. H. Press, S. A. Teukolsky and B. P. Flannery, Numerical Recipes in C++ (Cambridge University Press, Cambridge, 2002).
28. P. D. Drummond and C. W. Gardiner, Journal of Physics A: Mathematical and General, 13, 2353 (1980).
29. I. Carusotto, Y. Castin,  and J. Dalibard, PHYSICAL REVIEW A, 63 (2001).
30. P. Drummond, P. Deuar,  and K. Kheruntsyan, PHYSICAL REVIEW LETTERS, 92 (2004).
31. L. I. Plimak, M. K. Olsen, M. Fleischhauer,  and M. J. Collett, Europhys. Lett., 56, 372 (2001).
32. S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G. V. Shlyapnikov,  and M. Lewenstein, Phys. Rev. Lett., 83, 5198 (1999).
33. J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, S. L. Rolston, B. I. Schneider,  and W. D. Phillips, Science, 287, 97.
34. B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark,  and E. A. Cornell, Phys. Rev. Lett., 86, 2926 (2001).
35. Z. Dutton, M. Budde, C. Slowe,  and L. V. Hau, 293, 663 (2001).
36. A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, M. K. Oberthaler, D. J. Frantzeskakis, G. Theocharis,  and P. G. Kevrekidis, Phys. Rev. Lett., 101, 130401 (2008).
37. S. Stellmer, C. Becker, P. Soltan-Panahi, E.-M. Richter, S. Dörscher, M. Baumert, J. Kronjäger, K. Bongs,  and K. Sengstock, Phys. Rev. Lett., 101, 120406 (2008).
38. Y. S. Kivshar and B. Luther-Davies, Physics Reports, 298, 81 (1998), ISSN 0370-1573.
39. T. Busch and J. R. Anglin, Phys. Rev. Lett., 84, 2298 (2000).
40. J. Dziarmaga, Z. P. Karkuszewski,  and K. Sacha, Journal of Physics B: Atomic, Molecular and Optical Physics, 36, 1217 (2003).
41. R. V. Mishmash and L. D. Carr, PHYSICAL REVIEW LETTERS, 103 (2009).
42. J. Esteve, C. Gross, A. Weller, S. Giovanazzi,  and M. K. Oberthaler, Nature, 455, 1216 (2008).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters