Scrambled Mean Field Approach to the Quantum Dynamics of Degenerate Bose Gases

TL;DR

This paper introduces a new method for modeling the quantum dynamics of degenerate Bose gases that extends beyond traditional mean field approaches, capturing phenomena like coherence collapse with simplified equations.

Contribution

It develops a scrambled mean field approach using a multimode coherent-state ansatz that accurately models quantum effects beyond the Gross--Pitaevskii equation.

Findings

Captures coherence collapse in Bose gases

Provides equations simpler than traditional many-body methods

Extends modeling capabilities beyond mean field approximation

Abstract

We present a novel approach to modeling dynamics of trapped, degenerate, weakly interacting Bose gases beyond the mean field limit. We transform a many-body problem to the interaction representation with respect to a suitably chosen part of the Hamiltonian and only then apply a multimode coherent-state ansatz. The obtained equations are almost as simple as the Gross--Pitaevskii equation, but our approach captures essential features of the quantum dynamics such as the collapse of coherence.

Figures

Equations

$$i\hslash \frac{\partial }{\partial t}|\tilde{\Psi }⟩=\hat{\tilde{V}}(t)|\tilde{\Psi }⟩,$$

$$\hat{\tilde{V}}(t)={\hat{U}}_{\mathrm{r}}^{-1}(t)\hat{V}{\hat{U}}_{\mathrm{r}}(t).$$

$${\hat{U}}_{\mathrm{r}}^{-1}(t){\hat{b}}_j{\hat{U}}_{\mathrm{r}}(t)={\hat{\mathcal{A}}}_j{\hat{b}}_j,{\hat{U}}_{\mathrm{r}}^{-1}(t){\hat{b}}_j^{†}{\hat{U}}_{\mathrm{r}}(t)={\hat{b}}_j^{†}{\hat{\mathcal{A}}}_j^{†},$$

$$\hat{\cal A}_{j}=\exp\bigg{(}-i\omega_{j}t-i\sum_{j^{\prime}}g_{jj^{\prime}}\hat{b}^{\dagger}_{j^{\prime}}\hat{b}_{j^{\prime}}t\bigg{)}.~{}$$

$$|\tilde{\Psi }⟩=\prod _j\left[e^{-|{\psi }_j{|}^2/2}\sum _{n_j=0}^{\infty }\frac{({\psi }_j{\hat{b}}_j^{†}{)}^{n_j}}{n_j!}\right]|0⟩,$$

$$i\hslash \frac{\partial }{\partial t}{\psi }_j=\frac{\partial }{\partial {\psi }_j^{\ast }}⟨\tilde{\Psi }|\hat{\tilde{V}}(t)|\tilde{\Psi }⟩.$$

$$⟨\tilde{\Psi }|{\hat{\tilde{V}}}_{j_1^{\prime }\dots j_{l^{\prime }}^{\prime }}^{j_1\dots j_l}|\tilde{\Psi }⟩$$

$${\nu }_{j_1^{\prime }\dots j_{l^{\prime }}^{\prime }}^{j_1\dots j_l}$$

$$G_{j_1^{\prime }\dots j_{l^{\prime }}^{\prime };j}^{j_1\dots j_l}$$

$$\hat{H}$$

$${\hat{b}}_x$$

$${\hat{b}}_x^{†}$$

$${\hat{b}}_y$$

$${\hat{b}}_y^{†}$$

$${\omega }_{\mathrm{J}}^2(1-s^2)f_j(s)-{\omega }_{\Vert }^2\frac{d}{ds}\left[(1-s^2)\frac{d{f}_j(s)}{ds}\right]={\omega }_j^2{f}_j(s),$$

$${\hat{b}}_j$$

$${\hat{b}}_j^{†}$$

$$\hat{H}$$

$$g_{j{j}^{\prime }}=-\frac{U_{\mathrm{c}}}{4\hslash (1+{\delta }_{j{j}^{\prime }})}\frac{{\omega }_{\mathrm{J}}^2}{{\omega }_j{\omega }_j^{\prime }}\sum _{k=1}^{\mathcal{M}}{w}_k(1-s_k^2)f_j^2(s_k)f_{j^{\prime }}^2(s_k).$$

Full text

The mean-field approximation has become a powerful tool for modeling dynamics of degenerate gases of weakly interacting bosonic atoms Dalfovo ; ufn98 ; Leggett-rev . In this approximation, a Bose–Einstein condensate (BEC) or, in a case of low-dimensional geometry, a quasicondensate is described by a complex-valued classical field subject to the time-dependent Gross-Pitaevskii equation (GPE). Thermal and even quantum zero-point fluctuations can be incorporated into this classical-field picture within the truncated Wigner approximation (TWA) TWA1 via initial conditions. Unfortunately, the TWA with quantum noise provides physically meaningful solutions on rather a restricted time scale TWA2 . An alternative to the mean-field calculations is given by the quantum Boltzmann equation that can be derived by the standard techniques of the non-equilibrium quantum field theory Gasenzer1 ; Werner ; Gasenzer2 . Buchhold and Diehl DiehlEPJD derived kinetic equations not only for populations, but also for anomalous correlations in phonon modes in one dimension (1D). However, the quantum field theory methods are developed for a bulk medium in the thermodynamic limit, where the excitation spectrum is continuous and on-shell self-energies have therefore a non-zero imaginary part. In experiment, the finite size of trapped ultracold atomic clouds makes their excitation spectra discrete with frequencies of different excitation modes being well resolvable Jin96 ; Mewes96 . The multiconfigurational time-dependent Hartree method for bosons (MCTDHB) Alon08 is suitable for numerical modeling the non-mean-field dynamics of finite-size systems. However, it seems that the MCTDHB well describes the experimental data only in situations where the number of involved configurations remains small because of limitations specific for a given process, such as parametric excitation of a Bose–Einstein condensate (BEC) Hulet19 , and remains otherwise suitable mainly for few-body problems. A recently developed truncated conformal space approach Takacs2 ; Takacs1 works well at relatively low excitation energies of a bosonic system, as numerical diagonalization methods in general do.

Experiments with ultracold bosonic gases exhibiting effects beyond the mean field include, first of all, observations of the collapse and revival dynamics in optical lattices Bloch1 ; Bloch2 . Moreover, redistribution of atomic population between lattice sites accompanying this phenomenon has been detected in a recent experiment Zhou . The available theory phase-diff ; Imamoglu-theor ; Kuklov does not account for multimode aspects of the problem and remains a matter-wave analog of the well-known Jaynes–Cummings model in quantum optics cum-collapse ; eber-collapse .

The multimode approach that we present here can be called a scrambled mean field method. It bears certain similarities to the rotated Hartree method Cederbaum87 , but is remarkably simpler. Consider a Hamiltonian $\hat{H}=\hat{H}_{\mathrm{h}}+\hat{H}_{\mathrm{a}}$ that describes a BEC or a low-dimensional quasicondensate in collective variables Dalfovo ; ufn98 ; Leggett-rev . Its harmonic part, $\hat{H}_{\mathrm{h}}$, can be written after diagonalization as $\hat{H}_{\mathrm{h}}=\hbar\sum_{j}\omega_{j}\hat{b}_{j}^{\dagger}\hat{b}_{j}$, where $\hat{b}_{j}^{\dagger}$ and $\hat{b}_{j}$ are the creation and annihilation operators, respectively, of excitation quanta in the $j$th elementary mode with the fundamental frequency $\omega_{j}$ ($\hat{b}_{j}^{\dagger}$ and $\hat{b}_{j}$ obey the standard bosonic commutation rules). We denote the eigenstates of $\hat{H}_{\mathrm{h}}$, i.e., Fock states of elementary excitations, by $|\mathbf{n}\rangle=|\{n_{1},\,n_{2},\,n_{3},\,\dots\}\rangle$, where $\hat{b}_{j}^{\dagger}\hat{b}_{j}|\mathbf{n}\rangle=n_{j}|\mathbf{n}\rangle$. The anharmonic part $\hat{H}_{\mathrm{a}}$ of the Hamiltonian describes interaction between elementary excitations. This interaction is assumed to be small in order to make elementary excitations well defined. Usually $\hat{H}_{\mathrm{a}}$ can be expanded in Taylor series in $\hat{b}_{j}^{\dagger}$ and $\hat{b}_{j}$ beginning from a cubic term in a general case. The first-order perturbative correction to the energy of the state $|\mathbf{n}\rangle$ is given by the matrix element $\langle\mathbf{n}|\hat{H}_{\mathrm{a}}|\mathbf{n}\rangle$. The lowest order term that contributes to $\langle\mathbf{n}|\hat{H}_{\mathrm{a}}|\mathbf{n}\rangle$ is a quartic one, more precisely, its diagonal part $\hat{H}_{\mathrm{qd}}=\frac{\hbar}{2}\sum_{j,j^{\prime}}g_{jj^{\prime}}\hat{b}^{\dagger}_{j}\hat{b}^{\dagger}_{j^{\prime}}\hat{b}_{j^{\prime}}\hat{b}_{j}$, so that $\langle\mathbf{n}|\hat{H}_{\mathrm{a}}|\mathbf{n}\rangle=\frac{\hbar}{2}\sum_{j,j^{\prime}}g_{jj^{\prime}}n_{j}(n_{j^{\prime}}-\delta_{jj^{\prime}})+\Delta E_{\mathbf{n}}$, where $\Delta E_{\mathbf{n}}$ is the contribution of higher-order terms and $g_{jj^{\prime}}=g_{j^{\prime}j}$.

Now we need to introduce a unitary operator that induces quantum correlations between the modes. In contrast to Ref. Cederbaum87 , it contains no time-dependent parameters to be determined from the variational principle, but we derive instead its form from the perturbation theory considerations. We rearrange the terms in the Hamiltonian as $\hat{H}=\hat{H}_{0}+\hat{V}$, where $\hat{H}_{0}=\hat{H}_{\mathrm{h}}+\hat{H}_{\mathrm{qd}}$ and $\hat{V}=\hat{H}_{\mathrm{a}}-\hat{H}_{\mathrm{qd}}$. The evolution of the wave function $|\Psi\rangle$ os the system is governed by the Schrödinger equation $i\hbar\frac{\partial}{\partial t}|\Psi\rangle=\hat{H}|\Psi\rangle$. Our first step is to introduce the interaction representation with respect to the diagonalizable anharmonic Hamiltonian $\hat{H}_{0}$. This is done by an unitary transformation $|\Psi\rangle=\hat{U}_{\mathrm{r}}(t)|\tilde{\Psi}\rangle$, where $\hat{U}_{\mathrm{r}}(t)=\exp(-i\hat{H}_{0}t/\hbar)$. This transformation does not mix different Fock states of elementary excitations, but induces correlations between modes via the energy shift for elementary excitations depending on the quantum state of all the other modes and hence “scrambles” $|\tilde{\Psi}\rangle$. After this transformation the Schrödinger equation reads as

[TABLE]

where

[TABLE]

The field operators for the modes transform as

[TABLE]

[TABLE]

Eqs. (3, 4) allow us to write $\hat{\tilde{V}}$ in explicit form.

We assume that initially, at $t=0$, the state of the system is a product of coherent states (normalized to unity eigenstates of the respective annihilation operators Glauber ) for each mode. This type of initial conditions is also assumed in the mean field theory. Now we make the variational ansatz in the coherent state form not for $|\Psi\rangle$, but for the wave function in the interaction representation:

[TABLE]

where $|\mathbf{0}\rangle$ is the vacuum state for all the modes. By minimizing the action ${\cal S}=\int dt\,\langle\tilde{\Psi}|\left(i\hbar\frac{\partial}{\partial t}-\hat{\tilde{V}}\right)|\tilde{\Psi}\rangle$ we obtain the evolution equations for the complex functions $\psi_{j}(t)$ parametrizing the coherent states in Eq. (5):

[TABLE]

To give a recipe for the calculation of $\langle\tilde{\Psi}|\hat{\tilde{V}}|\tilde{\Psi}\rangle$, we assume that $\hat{V}$ can be expanded in Taylor series in creation and annihilation operators and consider a term $\hat{V}^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}}}=h^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}}}\prod_{\kappa=1}^{l}\hat{b}_{j_{\kappa}}^{\dagger}\prod_{\kappa^{\prime}=1}^{l^{\prime}}\hat{b}_{j^{\prime}_{\kappa^{\prime}}}$, where $h^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}}}$ is a constant. A straightforward calculation based on elementary properties of coherent states yields

[TABLE]

where

[TABLE]

The number of the system modes to be taken into account is to be determined in each particular case from physics considerations. A good guidance can be obtained from the thermalization argument. We calculate the mean energy $E=\langle\Psi|\hat{H}|\Psi\rangle$ at $t=0$, assume that the system equilibrates at $t\rightarrow\infty$, and calculate the temperature that corresponds to the internal energy $E$ (if the total number of elementary excitation is conserved, we need to determine also the chemical potential). The number of modes with the mean number of quanta larger than 1 at thermal equilibrium will give an estimation for the minimal number of modes to be considered.

After solving Eq. (6) with the initial conditions $\psi_{j}(0)=\psi_{0j}$, we can the find quantum-mechanical expectation value $\langle\hat{\cal O}\rangle\equiv\langle\Psi|\hat{\cal O}|\Psi\rangle$ for any observable $\hat{\cal O}$ as a function of time. For example, $\langle\hat{b}_{j}+\hat{b}_{j}^{\dagger}\rangle=\psi_{j}\exp\{-i\omega_{j}t-\sum_{j^{\prime}}|\psi_{j^{\prime}}|^{2}[1-\exp(ig_{jj^{\prime}}t)]\}+\mathrm{c.c.}\,$.

We test our method on the Hamiltonian of a two-dimensional (2D) harmonic oscillator with a quartic perturbation:

[TABLE]

The co-ordinates $x$ and $y$ in Eq. (10) are dimensionless. Since this Hamiltonian describes only two modes, many its eigenstates $|\Phi_{\lambda}\rangle$ and respective eigenenergies $\epsilon_{\lambda}$ can be found numerically with a high precision up to pretty high values of $\epsilon_{\lambda}$. In Fig. 1 we show the quantum mechanical mean value and variance of $x$ obtained by solving Eq. (6) in compartison with the results directly following from the expansion $|\Psi(t)\rangle=\sum_{\lambda}e^{-i\epsilon_{\lambda}t/\hbar}|\Phi_{\lambda}\rangle\langle\Phi_{\lambda}|\Psi(0)\rangle$. The same data for $y$ as well as the covariance $\mathrm{Cov}\{xy\}=\langle xy\rangle-\langle x\rangle\langle y\rangle$ are shown in Supplemental Material SM . Our numerical method reproduces the behavior of the expectation values and second-order correlations quite well. As the initially coherent wave packet disperses in the anharmonic potential, its regular motion characterized by oscillating expectation values of the co-ordinates is damped and the co-ordinate variances reach their asymptotic values. We tested numerical energy conservation for our method and found $\langle\Psi|\hat{H}|\Psi\rangle$ not exhibiting a systematic drift and deviating from its initial value by 1% at maximum SM .

We choose the phase difference dynamics in an extended bosonic Josephson junction as the first application of our method to a system with a nontrivially large number of modes. We consider ultracold gas of bosonic atoms in a trap consisting of two tunnel-coupled atomic waveguides, to be referred as the left and right waveguides. The longitudinal trapping is harmonic with the fundamental frequency $\omega_{\|}$. In order to simplify the overview of the example, we make a few approximations. We assume that the number of atoms is small enough to neglect the dependence of the radial width of the atomic cloud on the local density Salasnich , but large enough to provide an inverted parabolic Thomas–Fermi longitudinal profile of the 1D atomic density. Also we consider only antisymmetric modes (out-of-phase motion in the left and right waveguides) and neglect their coupling to the symmetric (in-phase) modes, thus reducing our problem to an ultracold-atom implementation of the sine-Gordon model, but, in contrast to Ref. polk-sg , the mean 1D density profile is in our case non-uniform and proportional to $1-s^{2}$, where $s$ is the dimensionless (scaled to the Thomas–Fermi equilibrium radius $R_{\mathrm{TF}}$) longitudinal co-ordinate. The Hamiltonian reads then

[TABLE]

Here $\hat{\phi}$ is the operator of the local phase difference between the left and right waveguides, $\hat{\rho}$ is the operator of the conjugate (density-difference) variable. The latter is made dimensionless by scaling to $R_{\mathrm{TF}}^{-1}$, since we use the dimensionless co-ordinate $s$, so that the commutation relation $[\hat{\rho}(s),\hat{\phi}(s^{\prime})]=i\delta(s-s^{\prime})$ holds. The charging energy $U_{\mathrm{c}}$ is positive, since we consider repulsive interactions between atoms. In repulsively interacting quantum gases at low energies the kinetic energy is dominated by phase fluctuations Dalfovo ; ufn98 , therefore we omit a term proportional to $[\partial\hat{\rho}/(\partial s)]^{2}$ in Eq. (11) from the very beginning. The Josephson oscillation frequency corresponding to the peak mean density (at $s=0$) is denoted by $\omega_{\mathrm{J}}$.

In Fig. 2 the results of modeling of the mean and variance of the global phase difference $\bar{\phi}=\frac{1}{2}\int_{-1}^{1}ds\,\hat{\phi}$ are presented (we drop the operator hat above $\bar{\phi}$ to keep notation simple; this observable can be measured by standard experimental techniques Pigneur ). Regular Josephson oscillations decay and the quantum uncertainty of $\bar{\phi}$ becomes large compared to its zero-point level, but still within a range corresponding to high visibility of the integrated interference picture. Their damping of oscillations is not as fast and perfect as in the experiment Pigneur , perhaps, because the model Hamiltonian (11) designed to demonstrate the proof of principle is too simplified.

The numerical method to solve Eq. (6) is overviewed in Supplemental Materials SM . Its most non-trivial part is related to the evaluation of exponential factors appearing in the r.h.s. of Eq. (7). On first glance one may get an impression that, e.g., estimation of a quartic interaction Hamiltonian for ${\cal M}$ modes requires independent calculation of ${\cal O}({\cal M}^{4})$ different terms. However, this “curse of dimensionality” curde can be circumvented in a very efficient way. At times $g_{jj^{\prime}}t\lesssim 1$, which are sufficiently long to observe the collapse of coherence, we can make two simplifications. Firstly, we can neglect the non-commutativity of $\hat{b}_{j}$ and $\hat{\cal A}_{j^{\prime}}$ and, hence, set $\nu^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}}}\approx\sum_{\kappa=1}^{l}\omega_{l}-\sum_{\kappa^{\prime}=1}^{l^{\prime}}\omega_{l^{\prime}}$ [see Eqs. (3,  4,  8)]. Secondly, we can write $\exp[-|\psi_{j}|^{2}(1-e^{iG^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}};j}t})]\approx\exp[i|\psi_{j}|^{2}G^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}};j}t-\frac{1}{2}(|\psi_{j}|G^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}};j}t)^{2}]$, employ the identity $\exp[-\frac{1}{2}(|\psi_{j}|G^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}};j}t)^{2}]=\int_{-\infty}^{\infty}d\sigma\,\exp(i\sigma|\psi_{j}|G^{j_{1}\dots j_{l}}_{j^{\prime}_{1}\dots j^{\prime}_{l^{\prime}};j}t-\frac{1}{2}\sigma^{2})/\sqrt{2\pi}$, and replace the integration by numerical averaging over a normally distributed pseudorandom parameter $\sigma$. These two approximations radically reduce rank of tensors used in numerical evaluation of the r.h.s. of Eq. (6).

To summarize, we developed a novel approach to numerical simulation of the dynamics of finite-size bosonic systems beyond the mean field approximation. Our method is free from the curse of dimensionality and designed to evaluate time scales of the multimode quantum dynamics manifested through the collapse of coherence as well as expectation values and correlations of simple observables. The description of scattering of quanta into initially empty modes remains beyond its scope. The main advantage of our method compared to the multiconfiguration approaches Cederbaum90 ; TC2008 is its remarkable simplicity and numerical efficiency in terms of the computational time and resources; it can be applied to obtain qualitative estimations on such long time scales that the number of configurations needed for the solution by standard methods becomes impractically large. Our method may be used not only in physics of trapped ultracold atomic gases, but also in other fields, such as molecular and chemical physics.

The author thanks C. Lévêque, N. J. Mauser, J.-F. Mennemann, J. Schmiedmayer, and H.-P. Stimming for helpful discussion. This work is supported by the Wiener Wissenschafts- und Technologie Fonds (WWTF) via Grant No. MA16-066 (SEQUEX).