Source: https://www.groundai.com/project/truncated-wigner-method-for-bose-gases/
Timestamp: 2019-04-24 13:47:04+00:00

Document:
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.
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.
where in 1d the interaction strength g→g1d=2ℏω⊥a, the s-wave scattering length a, and V 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 ∫dx|Φ(x)|2=N+m/2, where m denotes the number of modes in the initial state and N the total number of atoms.
Unlike in the usual GPE, here Φ(x,t) 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) .
with ⟨^b†q^bq⟩=nBE(ϵq)=[exp(ϵq/kBT)−1]−1 denoting the ideal Bose-Einstein distribution. At T=0 we have nBE(ϵq)=0.
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 (βq,β∗q′) in the place of the quantum operators (^bq,^b†q′) in Eq. (3).
where N′s fluctuates in each realization with the ensemble average ⟨N′s⟩e=N′=∑q|vq|2 at T=0. In Eq. (6) we have transformed the symmetric ordering of the Wigner representation to quantum expectation values of normally-ordered operators by subtracting ⟨β∗qβq⟩e=1/2 from each mode. The ground-state atom number is then obtained from the fixed total atom number N, so that in each stochastic realization Nc=N−N′s and we set β0=√Nc+1/2. The ensemble average of the ground-state population is obtained from ⟨Nc⟩e=N−⟨N′s⟩e=N−N′.
Here Φ(x) is a stochastic representation of the full field operator for the atoms.
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.
where the ground-state density n0,W(x)=(Nqc+1/2)|ϕ0(x)|2.
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) .
Calculation of ⟨^n2j⟩ then, however, results in double integrals over the sites that can be computationally slow when performed over a large number of realizations.
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 uj(x),vj(x) 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.
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 T=0, 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.
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 Vϕ to half of the atom cloud, for time τ, so that in the corresponding classical case the light sheet imprints a phase jump of ϕc=Vϕτ/ℏ at x=0, preparing a dark soliton. Classically the imprinted soliton oscillates in a harmonic trap at the frequency ω/√2 Busch and Anglin (2000) with the initial velocity v/c=cos(ϕc/2), depending on ϕc and the speed of sound c. The soliton is stationary (dark) for ϕc=π, 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 N′/N. At T=0 we keep the nonlinearity Ng1d fixed, but adjust the ratio g1d/N. This is tantamount to varying the effective interaction strength γint=mg1d/ℏ2n Kheruntsyan et al. (2003) . We can also study the effects of thermal depletion by varying T.
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 |Φ|2 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 ⟨^x⟩ and its uncertainty Δx=√⟨^x2⟩−⟨^x⟩2 Martin and Ruostekoski (2010, 2010) . Our results are summarized in Fig. 1.
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 87Rb 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 s(0)=48ER to s(τ)=96ER. 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 T=0, 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.
Figure 2: The numerical solution of the lowest two HFB modes in a six-site optical lattice showing (a) u1(x) (dotted) and v1(x) (solid); (b) u2(x) (dotted) and v2(x) (solid) at s=24 and T≃5.5nK; (c) Relative atom number (or spin) squeezing at different lattice height between two central nearest-neighbor sites ξpq=[Δ(^np−^nq)]2(np+nq)/(4npnq). The different data sets correspond to temperatures (curves from top to bottom) T≃5.5nK, T≃4.5nK, T≃4.0nK, and T=0. The harmonic trap frequency is ω=2π×21Hz, the atom number N≃5000, and the lattice spacing d≃5.7μm.
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.
Figure 3: Differences between initial noise generation in TWA. Dark soliton dynamics in a 1d harmonic trap showing (a-b) the Wigner density |Ψw(x,t)|2 for individual stochastic realizations of TWA for ϕc=2, T=0, and N=50 (parameters as explained in IV.1). In (a) the initial state is generated within the Bogoliubov approximation and in (b) using the quasi-condensate formalism; (c) The quantum mechanical expectation values for the soliton position ⟨^x⟩ (solid lines) and its uncertainty δx (shaded regions). The lighter curve with larger oscillation amplitude corresponds to the Bogoliubov case and the darker one the quasi-condensate case.
We acknowledge financial support from EPSRC and Leverhulme Trust.
C. Gardiner and P. Zoller, Quantum Noise (Springer, Berlin, 1999).
A. Sinatra, C. Lobo, and Y. Castin, Journal of Physics B: Atomic, Molecular and Optical Physics, 35, 3599 (2002).
P. D. Drummond and A. D. Hardman, EPL (Europhysics Letters), 21, 279 (1993).
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).
L. Isella and J. Ruostekoski, Phys. Rev. A, 72, 011601 (2005).
L. Isella and J. Ruostekoski, Phys. Rev. A, 74, 063625 (2006).
P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, Advances in Physics, 57, 363 (2008).
A. Polkovnikov, Annals of Physics, 325, 1790 (2010), ISSN 0003-4916.
A. D. Martin and J. Ruostekoski, Phys. Rev. Lett., 104, 194102 (2010a).
J. Javanainen and J. Ruostekoski, arXiv:1104.0820 (2011).
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).
T. Kinoshita, T. Wenger, and D. S. Weiss, Nature, 440, 900 (2006).
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).
J. Mun, P. Medley, G. K. Campbell, L. G. Marcassa, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett., 99, 150604 (2007).
C. Gross, J. Estève, M. K. Oberthaler, A. D. Martin, and J. Ruostekoski, Phys. Rev. A, 84, 011609 (2011).
J. Ruostekoski and L. Isella, Phys. Rev. Lett., 95, 110403 (2005).
U. Shrestha, J. Javanainen, and J. Ruostekoski, Phys. Rev. A, 79, 043617 (2009).
A. D. Martin and J. Ruostekoski, New Journal of Physics, 12, 055018 (2010b).
M. Olshanii, Phys. Rev. Lett., 81, 938 (1998).
C. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases (Cambridge University Press, Cambridge, 2002).
D. A. W. Hutchinson, R. J. Dodd, and K. Burnett, Phys. Rev. Lett., 81, 2198 (1998).
N. P. Proukakis, S. A. Morgan, S. Choi, and K. Burnett, Phys. Rev. A, 58, 2435 (1998).
K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond, and G. V. Shlyapnikov, Phys. Rev. Lett., 91, 040403 (2003).
C. Mora and Y. Castin, Phys. Rev. A, 67, 053615 (2003).
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).
J. Javanainen and J. Ruostekoski, Journal of Physics A: Mathematical and General, 39, L179 (2006).
W. T. V. W. H. Press, S. A. Teukolsky and B. P. Flannery, Numerical Recipes in C++ (Cambridge University Press, Cambridge, 2002).
P. D. Drummond and C. W. Gardiner, Journal of Physics A: Mathematical and General, 13, 2353 (1980).
I. Carusotto, Y. Castin, and J. Dalibard, PHYSICAL REVIEW A, 63 (2001).
P. Drummond, P. Deuar, and K. Kheruntsyan, PHYSICAL REVIEW LETTERS, 92 (2004).
L. I. Plimak, M. K. Olsen, M. Fleischhauer, and M. J. Collett, Europhys. Lett., 56, 372 (2001).
S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G. V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett., 83, 5198 (1999).
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.
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).
Z. Dutton, M. Budde, C. Slowe, and L. V. Hau, 293, 663 (2001).
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).
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).
Y. S. Kivshar and B. Luther-Davies, Physics Reports, 298, 81 (1998), ISSN 0370-1573.
T. Busch and J. R. Anglin, Phys. Rev. Lett., 84, 2298 (2000).
J. Dziarmaga, Z. P. Karkuszewski, and K. Sacha, Journal of Physics B: Atomic, Molecular and Optical Physics, 36, 1217 (2003).
R. V. Mishmash and L. D. Carr, PHYSICAL REVIEW LETTERS, 103 (2009).
J. Esteve, C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler, Nature, 455, 1216 (2008).

References: V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V.