Efficient description of Bose-Einstein condensates in time-dependent rotating traps

TL;DR

This paper introduces an efficient analytical and numerical approach to model Bose-Einstein condensates in time-dependent rotating traps, aiding high-precision quantum sensing applications.

Contribution

It presents a novel combined analytical and numerical method for accurately predicting condensate dynamics in complex, time-dependent rotating traps.

Findings

Analytical expressions for spatial and momentum distributions.

Comparison with previous methods shows improved efficiency.

Results are directly applicable to experimental quantum sensors.

Abstract

Quantum sensors based on matter-wave interferometry are promising candidates for high-precision gravimetry and inertial sensing in space. The favorable source for the coherent matter waves in these devices are Bose-Einstein condensates. A reliable prediction of their dynamics, which is governed by the Gross-Pitaevskii equation, requires suitable analytical and numerical methods which take into account the center-of-mass motion of the condensate, its rotation and its spatial expansion by many orders of magnitude. In this chapter, we present an efficient way to study their dynamics in time-dependent rotating traps that meet this objective. Both, an approximate analytical solution for condensates in the Thomas-Fermi regime and dedicated numerical simulations on a variable adapted grid are discussed. We contrast and relate our approach to previous alternative methods and provide further results, such as analytical expressions for the one- and two-dimensional spatial density distributions and the momentum distribution in the long-time limit that are of immediate interest to experimentalists working in this field of research.

Tables

Table 1. Table 1: List of the grid settings, physical trap frequencies and 2D-interaction strength in terms of the simulation parameters.

number of discretization steps	$N_x=N_y=2^7$
physical length of the grid	$L_x=L_y=20a_x$
trap frequency along the $x$-axis	${\omega }_x$
trap frequency along the $y$-axis	${\omega }_y=ϵ{\omega }_x$
2D-interaction strength	$g_{2D}N=\overline{g}N\hslash {\omega }_x{a}_x^2$

Equations

$$\text{rmi}\hslash \frac{\partial \psi }{\partial t}=\left[-\frac{{\hslash }^2}{2m}{\nabla }_{\text{bi}x}^2+V(t,\text{bi}x)+g|\psi (t,\text{bi}x){|}^2\right]\psi (t,\text{bi}x),$$

$${\int }_{\text{mathds}{R}^d}|\psi (t,\text{bi}x){|}^2{\mathrm{d}}^{d}x=N.$$

$$V(t,\text{bi}x)=V(t,\mathbf{\rho }(t))-\text{bi}F(t)[\text{bi}x-\mathbf{\rho }(t)]+\frac{m}{2}[\text{bi}x-\mathbf{\rho }(t){]}^{\mathrm{T}}{\Omega }^2(t)[\text{bi}x-\mathbf{\rho }(t)].$$

$$\psi (t,\text{bi}x)={\text{rme}}^{\frac{\text{rmi}}{\hslash }{\mathcal{S}}_1(t)}{\text{rme}}^{\frac{\text{rmi}}{\hslash }\left[\mathcal{P}(t)\text{bi}x-\mathcal{R}(t)\widehat{\text{bi}}p\right]}{\psi }_{\mathrm{D}}(t,\text{bi}x)$$

$$\frac{\mathrm{d}\mathcal{R}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}$$

$$\mathcal{L}(\mathcal{R},\dot{\mathcal{R}},t)=\frac{m}{2}{\dot{\mathcal{R}}}^2-V(t,\mathbf{\rho })+\text{bi}F(\mathcal{R}-\mathbf{\rho })-\frac{m}{2}(\mathcal{R}-\mathbf{\rho }{)}^{\mathrm{T}}{\Omega }^2(\mathcal{R}-\mathbf{\rho }).$$

$${\mathcal{S}}_k(t)={\int }_0^t\mathcal{L}(\mathcal{R},\dot{\mathcal{R}},t^{\prime })\mathrm{d}{t}^{\prime }-\frac{k}{2}\left[\mathcal{R}(t)\mathcal{P}(t)-\mathcal{R}(0)\mathcal{P}(0)\right],$$

$$\text{rmi}\hslash \frac{\partial {\psi }_{\mathrm{D}}}{\partial t}=\left[-\frac{{\hslash }^2}{2m}{\nabla }_{\text{bi}x}^2+\frac{m}{2}\text{bi}{x}^{\mathrm{T}}{\Omega }^2(t)\text{bi}x+g|{\psi }_{\mathrm{D}}(t,\text{bi}x){|}^2\right]{\psi }_{\mathrm{D}}(t,\text{bi}x)$$

$$t$$

$$\text{bi}x$$

$${\psi }_{\mathrm{D}}(t,\text{bi}x)=\frac{1}{\sqrt{\det \Lambda (\tau )}}{\text{rme}}^{\frac{\text{rmi}}{\hslash }\left[{\mathbf{\xi }}^{\mathrm{T}}A(\tau )\mathbf{\xi }-\beta (\tau )\right]}{\psi }_{\mathrm{\Lambda }}(\tau ,\mathbf{\xi })$$

$$\beta (\tau )\equiv {\int }_0^{\tau }\frac{\mu }{\det \Lambda ({\tau }^{\prime })}\mathrm{d}{\tau }^{\prime }$$

$$A(\tau )\equiv \frac{m}{2}{\Lambda }^{\mathrm{T}}(\tau )\frac{\mathrm{d}\Lambda }{\mathrm{d}\tau }.$$

$$\text{rmi}\hslash \frac{\partial {\psi }_{\mathrm{\Lambda }}}{\partial \tau }=$$

$${\Lambda }^{\mathrm{T}}(\tau )\left[\frac{{\mathrm{d}}^2\Lambda }{\mathrm{d}{\tau }^2}+{\Omega }^2(\tau )\Lambda (\tau )\right]=\frac{{\Omega }^2(0)}{\det \Lambda (\tau )}.$$

$$\Lambda (0)=\text{mathds}1\text{mbox}and{\frac{\mathrm{d}\Lambda }{\mathrm{d}\tau }|}_{\tau =0}=0$$

$${\Lambda }^{\mathrm{T}}(\tau )\frac{\mathrm{d}\Lambda }{\mathrm{d}\tau }=\frac{\mathrm{d}{\Lambda }^{\mathrm{T}}}{\mathrm{d}\tau }\Lambda (\tau ),$$

$$\psi (t,\text{bi}x)=\frac{1}{\sqrt{\det \Lambda (t)}}{\text{rme}}^{\text{rmi}\Phi (t,\text{bi}x)}{\psi }_{\mathrm{\Lambda }}\left(t,{\Lambda }^{-1}(t)\left[\text{bi}x-\mathcal{R}(t)\right]\right),$$

$$\Phi (t,\text{bi}x)\equiv \frac{1}{\hslash }\left\{{\mathcal{S}}_1(t)-\beta (t)+\mathcal{P}(t)\left[\text{bi}x-\frac{1}{2}\mathcal{R}(t)\right]+\frac{m}{2}{\left[\text{bi}x-\mathcal{R}(t)\right]}^{\mathrm{T}}C(t)\left[\text{bi}x-\mathcal{R}(t)\right]\right\}$$

$$C(t)\equiv \frac{\mathrm{d}\Lambda }{\mathrm{d}t}{\Lambda }^{-1}(t).$$

$${\psi }_{\mathrm{\Lambda }}(t,\text{bi}x)=\sqrt{\det \Lambda (t)}{\text{rme}}^{-\text{rmi}\Phi (t,\Lambda (t)\text{bi}x+\mathcal{R}(t))}\psi \left(t,\Lambda (t)\text{bi}x+\mathcal{R}(t)\right).$$

$${\Omega }^2(t)\equiv \mathrm{diag}\left[{\omega }_1^2(t),\dots ,{\omega }_d^2(t)\right],$$

$$\frac{{\mathrm{d}}^2{\lambda }_i}{\mathrm{d}{\tau }^2}+{\omega }_i^2(\tau ){\lambda }_i(\tau )=\frac{{\omega }_i^2(0)}{{\lambda }_i(\tau ){\prod }_{k=1}^d{\lambda }_k(\tau )},$$

$$\lambda_{i}(0)=1\quad\mbox{and}\quad\frac{\mathrm{d}\lambda_{i}}{\mathrm{d}\tau}\Big{|}_{\tau=0}=0\,.$$

$${\nabla }_{\text{bi}x}V(t,\text{bi}x){|}_{\text{bi}x=\mathbf{\rho }(t)}=0$$

$$\mathcal{R}(0)$$

$$\mathcal{P}(0)$$

$$m\ddot{\mathbf{\rho }}(t)=-{\nabla }_{\text{bi}x}V(t,\text{bi}x){|}_{\text{bi}x=\mathbf{\rho }(t)}$$

$$\mathcal{R}(t)$$

$$\mathcal{P}(t)$$

$$\mu \varphi (\text{bi}x)=\left\{-\frac{{\hslash }^2}{2m}{\nabla }_{\text{bi}x}^2+\frac{m}{2}[\text{bi}x-\mathcal{R}(0){]}^{\mathrm{T}}{\Omega }^2(0)[\text{bi}x-\mathcal{R}(0)]+g|\varphi (\text{bi}x){|}^2\right\}\varphi (\text{bi}x),$$

Full text

In this chapter, we present an efficient way to study their dynamics in time-dependent rotating traps that meet this objective. Both, an approximate analytical solution for condensates in the Thomas-Fermi regime and dedicated numerical simulations on a variable adapted grid are discussed. We contrast and relate our approach to previous alternative methods and provide further results, such as analytical expressions for the one- and two-dimensional spatial density distributions and the momentum distribution in the long-time limit that are of immediate interest to experimentalists working in this field of research.

Keywords: Bose-Einstein condensate, Gross-Pitaevskii equation, Thomas-Fermi approximation, scaling approach, time-dependent rotating trap, numerical simulation, Hamiltonian formalism, integrated density distribution

Date: July 7th, 2017

Preprint submitted to Elsevier

1 Introduction

Since the first creation of a Bose-Einstein condensate (BEC) in the mid 1990s [Anderson95, Davis95], the field of ultra-cold quantum gases [DGPS99, Giorgini2008] has enormously developed and BECs are nowadays commonly used in a broad variety of applications. They include the generation of vortices [Fetter2009], the exploration of different quantum phase transitions [Bloch2008] as well as inertial sensors that are realized with the help of matter-wave interferometry [Berman1997, Cronin2009, Tino2014]. In many of these phenomena the Gross-Pitaevskii (GP) equation [Gross61, Gross63, Pitaevskii61] provides a reliable theoretical description of the BEC dynamics. Analytic solutions of the GP equation can be derived within the Thomas-Fermi (TF) approximation when combined with the so-called scaling approach [kagan96, castin96, kagan97], where most of the BEC dynamics is described by an appropriate time-dependent coordinate transformation. In this chapter we generalize this scaling approach and apply it directly to the GP equation to obtain an efficient analytic description of the dynamics of a BEC in time-dependent rotating traps.

Our approach was motivated and used by seminal experiments realized within the QUANTUS collaboration [VanZoest2010, Muentinga2013], which has successfully performed matter-wave interferometry with BECs in microgravity at the drop tower in Bremen (ZARM). These pioneering experiments study and manipulate BECs that expand freely for several seconds and hence evolve into matter waves of macroscopic dimensions. In order to reliably predict and describe the outcome of such experiments, it is essential to have dedicated analytical and numerical tools at hand which take into account translational and rotational motions of the BEC as well as the fact that the spatial size of the condensate changes by many orders of magnitude during its free expansion. To account for these effects, the scaling approach has been generalized within the hydrodynamic framework [storey00, Edwards2002] by employing an affine transformation, where a translation first eliminates the center-of-mass motion and a subsequent linear map absorbs most of the remaining dynamics of the BEC. Here we apply this affine transformation directly to the GP equation resulting in an efficient description of the time evolution of a BEC in a time-dependent quadratic potential including slow rotations. Our approach has a straightforward application to matter-wave interferometry [Roura2014], facilitates an efficient numerical computation of the condensate wave function and provides valuable analytical insights into the dynamics of BECs.

Our chapter is organized as follows. In section 2 we first present the affine transformation of the GP equation, where we displace the wave function to eliminate the center-of-mass motion and then apply a linear transformation of the coordinates to account for the inner dynamics of the BEC. The combination of the affine approach and the time-dependent TF approximation results in an approximate, but astonishingly accurate solution for the dynamics of a BEC subject to a time-dependent rotating trap. Moreover, we derive expressions for the one- and two-dimensional integrated density distributions, which are experimentally accessible through time-of-flight pictures.

To verify the accuracy of our efficient description of the BEC dynamics, we perform in section 3 full numerical simulations of the GP equation for a purely rotating trap as well as for the free expansion of an initially rotating BEC. We find an excellent agreement with our approximate analytical solution and show that the affine transformation can be used to improve the performance of numerical simulations by solving the transformed GP equation rather than the original one. Indeed, since the affine transformation itself does not contain any approximation, our technique does not alter the accuracy of the numerical simulations but speeds up the computation.

In section 4 we establish a connection between the dynamics of the time-dependent affine transformation and the corresponding Hamiltonian formalism triggered by [kagan97]. In this context, we discuss two constants of motion of the affine transformation matrix, which we relate to the conservation of the total energy and the angular momentum of a BEC.

Since almost every experiment dealing with BECs makes use of time-of-flight pictures, it is essential to have a thorough understanding of the dynamics of a BEC during its free expansion. Indeed, due to their mean-field interaction, BECs posses a more complex time-of-flight dynamics compared with the simple ballistic expansion known from non-interacting quantum gases. Therefore, in section 5 we compare the free time evolution of the GP equation with the Schrödinger equation and derive a relation between the long-time behavior of the momentum distribution of a freely expanding BEC and its initial spatial density distribution.

In order to keep our consideration self-contained but focused on the central ideas, we have included further details and calculations in three appendices. In A we sketch the derivation of the affinely transformed GP equation and show how to calculate the integrated density distributions of a BEC within the TF regime. The important steps to determine the total energy of a BEC and its angular momentum within the time-dependent TF approximation are presented in B. Finally, in C we discuss the behavior of the affine transformation matrix for a freely expanding BEC that is released from an isotropic harmonic trap.

2 Efficient description of the time evolution of a Bose-Einstein condensate

The present section introduces a natural generalization of the scaling approach [kagan96, castin96, kagan97] to the case of rotating harmonic traps. In contrast to the hydrodynamical approach [storey00], we carry out this generalization directly to the GP equation. We apply these results to a BEC within the TF regime and obtain an approximate, analytical description of its time evolution, the so-called time-dependent TF approximation. Moreover, we then derive simple expressions for the integrated density distributions of a BEC within the TF regime which establish a direct link to experimental time-of-flight observations. We conclude this section with a brief discussion of an alternative, but equivalent description [Edwards2002] of the internal dynamics of a BEC. Further details on the results presented in this section are collected in A.

2.1 Affine transformation of the Gross-Pitaevskii equation

The affine transformation of the GP equation is realized in two steps [PhD_ME], namely (i) we eliminate the center-of-mass motion using Kohn’s theorem [Kohn1961, Dobson1994, Bialynicki-Birula2002], and (ii) we introduce a linear mapping of the coordinates in order to account for the main contributions to the internal dynamics of a BEC. Moreover, we show how this linear mapping is reduced to the well-known scaling approach [kagan96, castin96, kagan97] in the case of a non-rotating harmonic potential.

2.1.1 Basic setting

We describe the dynamics of a BEC in an inertial frame of reference by the macroscopic wave function $\psi(t,\bi{x})$ which satisfies the GP equation

[TABLE]

where the position vector $\bi{x}$ is considered here as an element of the $d$-dimensional vector space $\mathds{R}^{d}$ with $d=1,2,{\rm or}\,3$. In this way, our results are applicable to the general three-dimensional case [cornell796, ketterle896, Stamper-Kurn1998, DGPS99, Ketterle1999, Hechenblaikner2002, Edwards2002, Clancy2007] as well as to the cases of BECs confined in one or two dimensions [Petrov2000, Plaja2002, Bongs2001, Goerlitz2001]. Throughout this article the macroscopic wave function $\psi(t,\bi{x})$ is normalized to the number of particles $N$ in the condensate, that is

[TABLE]

The general potential $V(t,\bi{x})$ in Eq. (1) describes the interaction of an atom of mass $m$ with the external fields that correspond for example to a magneto-optical trap and Earth’s gravity. Moreover, we assume that the atoms interact with each other via a repulsive contact interaction, leading to the non-linear term $g\,|\psi(t,\bi{x})|^{2}$ in Eq. (1) with the positive coupling constant $g$.

In many experimental situations a semiclassical treatment of the underlying quantum mechanical dynamics allows an accurate and transparent interpretation of the measurement results. Here, we apply this approach to study the BEC dynamics around a classical trajectory $\bm{\rho}(t)$ and expand the potential $V(t,\bi{x})$ into a Taylor series up to second order around this trajectory

[TABLE]

Any anharmonicity of the potential can be neglected, as long as the size of the BEC remains sufficiently small within the vicinity of the trajectory $\bm{\rho}(t)$. The latter can be associated with the center-of-mass motion of the atomic cloud or the minimum of the external potential $V(t,\bi{x})$, as discussed in more detail in section 2.1.6.

Each term in Eq. (3) has a clear meaning. Indeed, the zeroth order term $V(t,\bm{\rho}(t))$ represents the value of the potential along the trajectory $\bm{\rho}(t)$. The second term ${-\bi{F}(t)[\bi{x}-\bm{\rho}(t)]}$ corresponds to the force $\bi{F}(t)$ acting on the atoms at the point $\bi{x}=\bm{\rho}(t)$. The third term $\frac{1}{2}m[\bi{x}-\bm{\rho}(t)]^{\!{\rm T}}\Omega^{2}(t)[\bi{x}-\bm{\rho}(t)]$ is the purely quadratic trapping potential represented in terms of the symmetric positive definite matrix ${\Omega^{2}(t)}$. Its eigenvalues coincide with the squared trap frequencies $\omega_{i}^{2}(t)$ along the principal axes of the harmonic trap.

2.1.2 Elimination of the center-of-mass motion

We incorporate the center-of-mass motion of the BEC in a straightforward way and thereby eliminate all $\bm{\rho}$-dependent terms in the potential (3) by making use of the transformation

[TABLE]

from the original wave function $\psi(t,\bi{x})$ to the new one $\psi_{\mathrm{D}}(t,\bi{x})$. Note that $\hat{\bi}{p}\equiv-\rmi\hbar\bm{\nabla}_{\bi{x}}$ represents the momentum operator in the position representation.

The time-dependent vectors $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$ in Eq. (4) describe the time evolution of the center of mass of the condensate, see section 4.2.2 for more details, and obey the classical equations of motion

[TABLE]

which are the Hamilton equations corresponding to the Lagrangian function

[TABLE]

The gobal phase $\mathcal{S}_{1}(t)$ in Eq. (4) depends on the classical action via the generalized definition

[TABLE]

where the additional integer $k\in\mathds{Z}$ has been introduced for later purposes. As outlined in A.1, we arrive at the GP equation

[TABLE]

for the transformed wave function $\psi_{\mathrm{D}}(t,\bi{x})$ by inserting Eq. (4) into Eq. (1) and taking advantage of Eqs. (2.1.2)-(7) for $k=1$.

The decoupling of the center-of-mass motion of a BEC is possible as long as the potential is at most quadratic [Dobson1994, Bialynicki-Birula2002, Nandi07]. It no longer holds true for anharmonic potentials for which a nontrivial coupling of the center-of-mass motion and the inner dynamics of the condensate exists [dum98].

2.1.3 Linear transformation as a natural generalization of the

scaling approach

After having eliminated the center-of-mass motion, we perform a time-dependent linear transformation to account for the main internal dynamics of the BEC. For this purpose, we also introduce a new time coordinate $\tau$ such that the linear mapping between the “original coordinates” $(t,\bi{x})$ and the new ones $(\tau,\bm{\xi})$ reads

[TABLE]

with $\Lambda(\tau)$ being an arbitrary, time-dependent matrix. In what follows, we refer to $\Lambda(\tau)$ as “adaptive matrix” and to $(\tau,\bm{\xi})$ as “adapted coordinates”.

In analogy to Eq. (4), the coordinate transformation (9) goes hand in hand with the transformation

[TABLE]

to the so-called affinely transformed wave function $\psi_{\mathrm{\Lambda}}(t,\bi{x})$. The scalar phase $\beta(\tau)$ and the symmetric matrix $A(\tau)$ introduced in this transformation depend on the adaptive matrix $\Lambda(\tau)$ via

[TABLE]

and

[TABLE]

The constant $\mu$ that appears in the definition (11) represents the chemical potential associated with the ground state of the initial BEC at the time $\tau=0$. In A.2 we outline the derivation of the affinely transformed GP equation

[TABLE]

for the wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$, where we have made use of the shorthand notation $\Lambda^{-{\rm T}}\equiv(\Lambda^{-1})^{\rm T}$ and the requirement that the matrix $\Lambda(\tau)$ obeys the nonlinear matrix differential equation

[TABLE]

The additional assumption that at $\bi{x}=\bm{\xi}$ and $t=0$ we have $\psi_{\mathrm{D}}(0,\bi{x})=\psi_{\mathrm{\Lambda}}(0,\bi{x})$ yields the initial conditions

[TABLE]

for the nonlinear matrix differential equation (14). As discussed in A.3, the symmetry of the matrix $A(\tau)$, Eq. (12), is connected to the so-called irrotationality condition

[TABLE]

that gives rise to $d(d-1)/2$ constants of motion of the matrix differential equation (14).

So far, we have made no approximation in deriving the affinely transformed GP equation (13). Although it looks much more complicated than Eq. (8), the solutions of Eq. (13) show almost no time dependence within the TF regime. This observation enables us to establish an efficient description of the BEC dynamics, as discussed in subsection 2.2.

2.1.4 The affinely transformed wave function

By combining the two transformations given by Eqs. (4), (9) and (10), we obtain the following relation between the original and the affinely transformed wave function

[TABLE]

where we have introduced the local phase

[TABLE]

and the time-dependent, symmetric matrix

[TABLE]

Since the matrix $C(t)$ determines the time evolution of the quadratic phase term, we refer to it as “quadratic phase matrix”. For later purposes we also note, that by inverting Eq. (17), we find the affinely transformed wave function $\psi_{\mathrm{\Lambda}}(t,\bi{x})$ in terms of $\psi(t,\bi{x})$ via

[TABLE]

In summary, the affine transformation of the macroscopic wave function is realized by Eq. (17) and leads to the affinely transformed GP equation (13). In order to determine $\psi(t,\bi{x})$, one has to solve the time-dependent partial differential Eq. (13) together with the ordinary differential equations (2.1.2) and (14) for the center-of-mass variables $\bm{\mathcal{R}}(t)$, $\bm{\mathcal{P}}(t)$ and the adaptive matrix $\Lambda(t)$, respectively. Throughout this article, we refer to this generalization of the standard scaling method [kagan96, castin96, kagan97] as the “affine approach”.

2.1.5 Non-rotating trap: connection to the scaling approach

We continue by pointing out the connection between our affine approach and the scaling method introduced by [kagan96, castin96, kagan97]. We consider a non-rotating trap and assume for simplicity that the principal axes of the harmonic potential coincide with the coordinate axes. This can always be achieved by an appropriate choice of coordinates for which the trap matrix $\Omega^{2}(t)$ in Eq. (3) possesses the diagonal form

[TABLE]

with $\omega_{i}(t)$ being the trapping frequency along the $x_{i}$ direction. Accordingly, we assume that the adaptive matrix $\Lambda(\tau)$ is of diagonal form and substitute $\Lambda(\tau)\equiv{\rm diag}\left[\lambda_{1}(\tau),\ldots,\lambda_{d}(\tau)\right]$ into Eq. (14) which yields the coupled nonlinear differential equations

[TABLE]

where $i\in 1,\dots,d$. The functions $\lambda_{i}(\tau)$ characterize the time evolution of the condensate in terms of an individual scaling along the three principal axes of the potential. From Eq. (15), we obtain for the corresponding initial conditions

[TABLE]

Clearly, the irrotationality condition, Eq. (16), is automatically satisfied for a diagonal adaptive matrix $\Lambda(\tau)$. For the case $d=3$, Eqs. (22) and (23) precisely rephrase the well-known equations derived for the scaling factors in Ref. [castin96, kagan97].

2.1.6 Interpretation of the classical trajectory

During the preparation of the BEC by laser and evaporative cooling at $t\leq 0$, we denote by $\bm{\rho}(t)$ the position of the minimum of the full external potential $V(t,\bi{x})$ which acts on the atomic cloud. For this reason, the classical trajectory $\bm{\rho}(t)$ is found as a solution of the equation

[TABLE]

for all times $t\leq 0$. Since the trapping potential $V(t,\bi{x})$ is located in an Earth-bound laboratory, a capsule freely falling in a drop tower or a satellite in space, inertial effects due to the local acceleration and rotation of the comoving frame of reference attached to the trapping potential $V(t,\bi{x})$ along $\bm{\rho}(t)$ do occur in general. However, here we assume that these inertial effects can be neglected for all times $t\leq 0$ due to the dominating influence of the trapping potential and the repulsive interaction of the atoms.

As pointed out in section 2.1.2, the time-dependent vectors $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$ are associated with the center-of-mass motion of the atomic cloud. This interpretation implies that the initial conditions $\bm{\mathcal{R}}(0)$ and $\bm{\mathcal{P}}(0)$ for the classical equations of motion (2.1.2) are directly linked to the preparation of the BEC in the external potential $V(t,\bi{x})$ according to

[TABLE]

In other words, the center of mass of the initial BEC at $t=0$ is supposed to be at rest within the comoving frame of reference attached to the external potential $V(t,\bi{x})$ at the position $\bm{\rho}(0)$.

After the initial preparation phase, the dynamical evolution of the BEC for times $t>0$ can be analyzed based on two different interpretations of the classical trajectory $\bm{\rho}(t)$. In case one is interested in the relative motion of the center of mass of the BEC with respect to the trap minimum, the trajectory $\bm{\rho}(t)$ can still be associated with the minimum of the external potential $V(t,\bi{x})$ and Eq. (24) has to be satisfied by $\bm{\rho}(t)$ also for $t>0$.

In terms of a semiclassical approach, one can likewise associate $\bm{\rho}(t)$ for all $t>0$ with the center-of-mass motion of the atomic cloud itself. In this case, the trajectory $\bm{\rho}(t)$ is determined as the solution of the classical equation of motion

[TABLE]

for $t>0$ and describes the semiclassical motion of the center of mass of the BEC in the full external potential $V(t,\bi{x})$. Thus, the time-dependent vectors $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$ follow directly from $\bm{\rho}(t)$ via

[TABLE]

whereas the classical equations of motion, Eq. (2.1.2), mathematically coincide with Eq. (26).

2.2 The macroscopic wave function in the time-dependent Thomas-Fermi approximation

We start this subsection by recalling the TF approximation [Fetter1998, pethick02] for the ground state of a BEC in a harmonic trap. We then apply this method also to the case of a time-dependent rotating trap in order to provide an efficient description of the dynamics of the macroscopic wave function of a BEC within the TF regime [kagan96, castin96, kagan97, storey00, DGPS99].

In addition, we present a simple relationship between the spatial density distribution of a BEC in three dimensions and its corresponding integrated density distributions in one and two dimensions, which allows us to fully characterize the BEC dynamics within the TF regime by a sequence of three mutually orthogonal time-of-flight pictures. We conclude this subsection by briefly discussing an alternative description of the internal dynamics of a rotating BEC within the TF regime [Edwards2002].

2.2.1 Thomas-Fermi approximation for the initial ground state

Without loss of generality, we assume that $V(0,\bm{\rho}(0))=0$ at the initial time $t=0$. Furthermore, we recall that $\bi{F}(0)=0$ for the harmonic potential (3) due to the validity of Eq. (24). Taking also Eq. (2.1.6) into account, the ground state $\phi(\bi{x})$ of a BEC at $t=0$ is defined as the solution of the stationary GP equation

[TABLE]

with $\mu$ being the chemical potential of the ground state.

The TF regime is characterized by a dominant contribution of the atomic interactions to the total energy of the BEC [Fetter1998, pethick02]. In this case, the kinetic energy term in the stationary GP equation (28) only accounts for a negligible contribution to the total energy in comparison with those arising from the harmonic potential and the atom-atom interactions. For this reason, within the TF approximation one simply neglects the kinetic energy term which gives rise to the approximate ground state of the BEC

[TABLE]

Here, the chemical potential $\mu_{\mathrm{TF}}$ within the TF approximation

[TABLE]

is determined by the normalization condition (2) for the wave function $\phi_{\mathrm{TF}}(\bi{x})$, $\mathrm{\Gamma}(x)$ denotes the Gamma function and

[TABLE]

stands for the positive part of $x$ to the power of $\alpha$, $\alpha>0$, with $\Theta(x)$ being the Heavyside function. The additional phase $\rmi\bm{\mathcal{P}}(0)\,\bi{x}/\hbar$ in Eq. (29) is due to the non-vanishing momentum $\bm{\mathcal{P}}(0)$, Eq. (2.1.6), of the initial BEC as discussed in section 2.1.6.

The ground state within the TF approximation $\phi_{\mathrm{\Lambda,TF}}(\bm{\xi})$ rewritten in terms of the adapted coordinates $\bm{\xi}$ follows directly from Eqs. (20) and (29) according to

[TABLE]

Besides the initial conditions for the adaptive matrix $\Lambda(t)$, Eq. (15), we have used the defining Eqs. (7) for $k=1$, (11) and (19) to find $\Phi(0,\bi{x})=\bm{\mathcal{P}}(0)\left[\bi{x}-{\textstyle\frac{1}{2}}\bm{\mathcal{R}}(0)\right]/\hbar$ for the local phase (18) at $t=0$ in the derivation of Eq. (32).

2.2.2 Time-dependent Thomas-Fermi approximation for the dynamical evolution

Within the time-dependent TF approximation [castin96, kagan97], one neglects the kinetic energy term in the affinely transformed GP equation (13) for all times $\tau\geq 0$ in order to arrive at the approximate, ordinary differential equation

[TABLE]

As we will show below, the solution $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ of Eq. (33) is time-independent and therefore agrees with the approximate ground state $\phi_{\mathrm{\Lambda,TF}}(\bm{\xi})$, Eq. (32), for all times $\tau\geq 0$

[TABLE]

This observation is the most important result of the time-dependent TF approximation and it manifests its full strength by providing an efficient and accurate description of the BEC dynamics within the TF regime.

In order to prove the validity of Eq. (34), we first multiply the ordinary differential equation (33) by $\psi^{*}_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ and subtract the imaginary part of the resulting equation to find

[TABLE]

Thus, the absolute value of the affinely transformed wave function does not change in time and we conclude that

[TABLE]

We infer from the Eqs. (32) and (36), that the square bracket on the right hand side of the ordinary differential equation (33) vanishes, giving rise to $\partial\psi_{\mathrm{\Lambda}}/\partial\tau=0$ which finally justifies the validity of Eq. (34).

To find the above mentioned efficient description of the BEC dynamics, we take advantage of Eq. (17) with the function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ being determined by combining the Eqs. (34) and (32). As a result, we obtain the central expression

[TABLE]

that characterizes the time evolution of a BEC within the time-dependent TF approximation. In Eq. (37) we have introduced the new local phase

[TABLE]

where the generalized action $\mathcal{S}_{2}(t)$ is given by Eq. (7) for $k=2$. Moreover, we have defined the so-called “TF matrix”

[TABLE]

which is positive definite.

According to Eq. (37), the boundary at which the spatial density distribution of the BEC vanishes corresponds to an ellipse if $d=2$ or to an ellipsoid if $d=3$. The TF matrix (39) defines the orientation and semi-principal axes of this ellipse or ellipsoid whereby its points $\bi{x}\in\mathds{R}^{d}$ satisfy the condition ${[\bi{x}-\bm{\mathcal{R}}(t)]^{{\rm T}}\Sigma^{-1}(t)[\bi{x}-\bm{\mathcal{R}}(t)]=2\mu_{\mathrm{TF}}/m}$. The eigenvalues $\sigma_{i}$ of the TF matrix $\Sigma(t)$ determine the lengths of the semi-principal axes, which coincide with the individual TF radii $r_{i}$ of the BEC, via $r_{i}=\sqrt{2\mu_{\mathrm{TF}}\,\sigma_{i}/m}$.

Fig. 1 illustrates the time evolution of such an ellipsoid in three dimensions. The initial density distribution of the ground state is centered around $\bm{\mathcal{R}}(0)$ with its principal axes being parallel to the coordinate axes, as shown in Fig. 1(a). Changes in the orientation and position of the quadratic potential (3) give rise to a center-of-mass motion $\bm{\mathcal{R}}(t)$ and a transformation of the elliptic contour of the BEC density distribution via the time evolution of the TF matrix $\Sigma(t)$, as depicted in Fig. 1(b).

We summarize this section by pointing out once more that the expression given by Eq. (37) provides a valuable description of the time evolution of a rotating BEC within the TF regime. Having the time-dependent TF approximation in mind, one might be tempted to believe that Eq. (37) is only valid for times shortly after $t=0$. However, Eq. (37) yields very accurate results also for freely expanding BECs, since in this case the adaptive matrix $\Lambda(\tau)$ grows linearly in time, as discussed in section 5.1, and therefore strongly reduces the influence of the kinetic energy term in Eq. (13), the so-called quantum pressure.

2.2.3 Integrated density distributions for time-of-flight pictures

In practice, absorption imaging techniques [Reinaudi2007] are the prevalent method to study the dynamics of a three-dimensional BEC. Identically prepared BECs are illuminated by a laser field at different times of flight and their shadows are detected by a CCD camera resulting in a sequence of two-dimensional pictures of the spatial density distribution of a BEC. For this reason, these techniques provide only partial information about the orientation and size of the original three-dimensional density distribution.

In this subsection we discuss the one- and two-dimensional density distributions that result from an originally three-dimensional density profile of a BEC within the TF regime. In addition, we sketch the determination of the underlying three-dimensional density distribution from three mutually orthogonal time-of-flight pictures of a single BEC. The high efficiency of our approach is due to the knowledge of the density profile of the macroscopic wave function, Eq. (37). In case one does not have any a priori knowledge about the wave function under consideration, a reconstruction is still possible with the help of the Radon transformation [Gindikin1994, Schleich2001].

We start with the density distribution of a three-dimensional BEC within the TF regime

[TABLE]

which directly follows from Eq. (37). The two-dimensional integrated density distribution in the $x$-$y$ plane

[TABLE]

can be evaluated as discussed in A.4 and reads

[TABLE]

Here, we have introduced the reduced vectors $\bar{\bi}{x}\equiv(x,y)^{\rm T}$ and $\bar{\bm{\mathcal{R}}}\equiv(\mathcal{R}_{1},\mathcal{R}_{2})^{\rm T}$ together with the symmetric submatrix

[TABLE]

which is obtained from the TF matrix $\Sigma$, Eq. (39), by removing the row $\Sigma_{3i}$ and the column $\Sigma_{i3}$ that correspond to the $z$-coordinate. In a similar way, the two-dimensional integrated density distributions $n^{\mathrm{2D}}_{\mathrm{TF}}(t,y,z)$ and $n^{\mathrm{2D}}_{\mathrm{TF}}(t,x,z)$ depend on the two-dimensional submatrices in which the rows and columns corresponding to the $x$- and $y$-coordinate have been eliminated, respectively.

For the one-dimensional integrated density distribution along the $x$-axis

[TABLE]

we obtain

[TABLE]

The one-dimensional integrated density distributions $n^{\mathrm{(1D)}}_{\mathrm{TF}}(t,y)$ and $n^{\mathrm{(1D)}}_{\mathrm{TF}}(t,z)$ possess a similar structure. A general derivation of the previous results for arbitrary dimensions is presented in A.4.

To reconstruct the unknown three-dimensional density distribution of a BEC within the TF regime, Eq. (40), from three mutually orthogonal time-of-flight pictures, one needs to fit the two-dimensional integrated density distributions $n^{\mathrm{(2D)}}_{\mathrm{TF}}(t,x,y)$, $n^{\mathrm{(2D)}}_{\mathrm{TF}}(t,y,z)$ and $n^{\mathrm{(2D)}}_{\mathrm{TF}}(t,x,z)$ to the corresponding experimental data and thereby determine the submatrix (43) and its two other counterparts. These three submatrices provide us with all elements of the TF matrix $\Sigma$, where each diagonal element of $\Sigma$ is found twice. This redundancy can be used to check the quality of the data obtained from the three mutually orthogonal time-of-flight pictures.

2.2.4 Alternative description of rotating and vortex-free condensates

Our description of the internal dynamics of a BEC within the TF regime rests upon the solution $\Lambda(t)$ of the nonlinear matrix differential equation (14) with the corresponding initial conditions, Eq. (15). Based on the adaptive matrix $\Lambda(t)$, we then determine the TF matrix $\Sigma(t)$ and the quadratic phase matrix $C(t)$ defined by Eqs. (39) and (19), respectively. The two symmetric matrices $\Sigma(t)$ and $C(t)$ characterize the time evolution of the absolute value and the phase of the macroscopic wave function $\psi(t,\bi{x})$, Eq. (37). Whereas the TF matrix $\Sigma(t)$ plays a crucial role in the context of time-of-flight pictures, the quadratic phase matrix $C(t)$ is of great importance for BEC based interferometry [Torii2000, Simsarian2000, Debs2011, Altin2013, Muentinga2013].

One might wonder whether there is a way to study the BEC dynamics directly in terms of the experimentally more accessible matrices $\Sigma^{-1}(t)$ and $C(t)$ without a link to the adaptive matrix $\Lambda(t)$ at all. In fact, such an alternative description is possible and involves the solution of the two coupled first order nonlinear matrix differential equations

[TABLE]

with the corresponding initial conditions

[TABLE]

These first order differential equations can be derived by taking the time derivative on both sides of Eqs. (19) and (39) and using the matrix differential equation (14), the two identities $(d\Sigma^{-1}/dt)=-\Sigma^{-1}(d\Sigma/dt)\Sigma^{-1}$ and ${\det{\Lambda(t)}=\sqrt{\det{\Omega^{2\!}(0)}/\det{\Sigma^{-1}(t)}}}$ as well as the symmetry of the matrix $C(t)$.

At present, only a simplified version [Edwards2002] of the system of equations (45) and (46) has been used to study the dynamics of a freely expanding BEC that was initially prepared in a rotating anisotropic harmonic trap [Hechenblaikner2002]. The trap was instantly turned off at the time $t=T_{\mathrm{off}}$. In fact, by inserting the trap matrix $\Omega^{2\!}(t)=0$ for all $t\geq T_{\mathrm{off}}$ into both Eqs. (45) and (46), it can be shown that the resulting matrix differential equations are mathematically equivalent to the equations of motion (12) and (13) derived in Ref. [Edwards2002]. However, the approach presented in this article facilitates an analysis of the time evolution of BECs within the TF regime that goes far beyond the scope considered in Refs. [Edwards2002, Hechenblaikner2002].

Next, we briefly outline how to obtain $\Lambda(t)$ from the quadratic phase matrix $C(t)$ and the TF matrix $\Sigma(t)$. When the time evolution of $C(t)$ is known, the adaptive matrix $\Lambda(t)$ follows from Eq. (19) by simple integration of the first order matrix differential equation

[TABLE]

with the initial condition $\Lambda(0)=\mathds{1}$. In contrast, the knowledge of the TF matrix $\Sigma(t)$ does not suffice to fully determine the adaptive matrix $\Lambda(t)$. From Eq. (39), the latter can only be identified up to an arbitrary orthogonal matrix $U(t)$ via the relation

[TABLE]

where the square root of $\Sigma(t)$ is defined in terms of its spectral decomposition. The specification of the orthogonal matrix $U(t)$ involves again the quadratic phase matrix $C(t)$.

There are two main reasons that make our approach based on the adaptive matrix $\Lambda(t)$ more suitable for the characterization of BEC dynamics within the TF regime than the alternative description just presented. First, the second order matrix differential equation (14) for $\Lambda(t)$ can be recast in the form of Hamilton’s equations, as shown in subsection 4.1. This fact allows us to apply the full mathematical machinery available for Hamiltonian mechanics to study the time evolution of $\Lambda(t)$. Second, for studying the density distribution of a BEC near its surface at the scale of the corresponding healing length, the TF ground state given by Eq. (29) provides no longer an adequate solution of Eq. (28) and the macroscopic wave function Eq. (37) does not suffice for this purpose. Hence, the TF matrix $\Sigma(t)$, which rests upon the validity of Eq. (37), loses its immediate physical significance. In contrast, the affine approach based on the adaptive matrix $\Lambda(t)$ provides a valuable tool to study scenarios that are beyond the scope of the TF approximation with the help of dedicated numerical simulations, as discussed in the next section.

3 Application of the affine approach to numerical simulations

In this section we show how numerical simulations of the BEC dynamics can benefit from the affine approach and we likewise apply this approach to quantify the accuracy of the time-dependent TF approximation for several different scenarios. For simplicity, we present here the study of the dynamics of rotating BECs on a 2D grid. We emphasize that full 3D simulations of BECs in the quasi-2D regime [petrov-2000, salasnich-2002, mateo-2008] have also been performed, which verify the validity of our 2D results.

3.1 Efficient simulation of the time evolution of a BEC

Solving numerically the time-dependent GP equation (1) in the case of a freely expanding BEC in the original coordinates $(t,\bi{x})$, one immediately faces the problem that the size of the condensate grows by several orders of magnitude after switching off the trap. To deal with this problem, a large and well resolved grid is required, leading to a huge increase of the computational costs of the numerical simulations, especially for two or three dimensional cases.

Here we present an alternative approach to overcome this problem. Namely we are solving the affinely transformed GP equation (13) in the adapted coordinates $(\tau,\bm{\xi})$, rather than the GP equation (1) in the original coordinates $(t,\bi{x})$. Since the external and internal dynamics of the BEC are handled by the time dependence of the center-of-mass position $\bm{\mathcal{R}}(t)$ and the adaptive matrix $\Lambda(\tau)$, respectively, only the time evolution of the affinely transformed wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ is left and can be computed very efficiently [PhD_ME]. As a result, in order to obtain the macroscopic wave function $\psi(t,\bi{x})$ in the original coordinates it is sufficient to apply the transformation Eq. (17) once at the end of the simulation. Moreover, our method is not limited to the case of free expansion, but can enhance the efficiency of numerical simulations for various experimental scenarios such as rotating traps or delta-kick collimation [ammann-1997].

3.2 Quantifying the accuracy of the time-dependent Thomas-Fermi approximation

Since the time-dependent TF approximation plays a key role in deriving the analytic expression, Eq. (37), for the wave function $\psi(t,\bi{x})$ of a BEC, it is worthwhile to have a closer look at its accuracy. The time-dependent TF approximation is based on assumption (34) meaning that the wave function in adapted coordinates $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ remains approximately in its initial value $\psi_{\mathrm{\Lambda}}(0,\bm{\xi})$ and does not undergo any time evolution at all. With our numerical simulations we have quantified the validity of this assumption for two different scenarios involving rotating traps.

We take advantage of the Bures distance to quantify how much the numerically obtained state $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ differs from the initial state $\psi_{\mathrm{\Lambda}}(0,\bm{\xi})$. In the following, we first recall the definition of the Bures distance and its properties and then describe the 2D model of a rotating BEC. Finally, we discuss the results of the simulations and characterize the parameter regime where it is safe to apply the TF approximation.

3.2.1 Bures distance

The Bures distance [Bures1969] of two macroscopic wave functions $\psi_{1}$ and $\psi_{2}$ is given by

[TABLE]

where we have used the normalized condition given by Eq. (2). First we note, that one obtains $B(\psi_{1},\psi_{2})=0$ if $\psi_{1}=\psi_{2}$ and $B(\psi_{1},\psi_{2})=\sqrt{2}$ if $\psi_{1}$ and $\psi_{2}$ are orthogonal. Second, for two $d$-dimensional spherically symmetric Gaussian wave packets

[TABLE]

with different widths $\sigma_{1}$ and $\sigma_{2}$, the Bures distance Eq. (48) reads

[TABLE]

In the case of two-dimensional Gaussian wave packets with the relative difference $\delta_{\sigma}\equiv(\sigma_{2}-\sigma_{1})/\sigma_{1}$ in their widths, Eq. (49) gives rise to

[TABLE]

Thus, a small relative change $\delta_{\sigma}$ in the width of a two dimensional Gaussian wave packet results in a small Bures distance, that is $B_{\mathrm{G}}^{(2)}(\delta_{\sigma})\cong\delta_{\sigma}\ll 1$. Despite the fact that the macroscopic wave function of a BEC is not typically Gaussian, this estimate can help us to evaluate the time dependence of the remaining inner dynamics which is not included in the affine approach.

3.2.2 2D model of a Bose-Einstein condensate in a rotating trap

The numerical simulations discussed in this section have all been performed with a 2D model based on realistic values that are accessible in state-of-the-art experiments. Our system consists of atoms of mass $m$ which are condensed and harmonically trapped along the $x$- and $y$-direction with the frequencies $\omega_{x}$ and $\omega_{y}$, respectively. The time evolution of this system is described by the affinely transformed GP equation (13) with $d=2$. The number of discretization steps in time and the physical length of the grid are chosen as a trade-off between accuracy and computational time.

The free parameters of our system are $(i)$ the anisotropy factor $\epsilon\equiv\omega_{y}/\omega_{x}$, $(ii)$ the strength $\bar{g}N$ of the zero-range interaction between the atoms, and $(iii)$ the final rotation rate $\dot{\varphi}_{\mathrm{end}}$ of the trap. All quantities are measured with respect to the chosen time scale $1/\omega_{x}$ and the length scale $a_{x}\equiv\sqrt{\hbar/(m\omega_{x})}$, respectively. Table 1 displays the parameters for the numerical simulations as well as the relations between the physical quantities and the free parameters. For convenience we refer to the dimensionless quantity $\bar{g}N=g_{\mathrm{2D}}N/(\hbar\omega_{x}a_{x}^{2})$ as the interaction strength as it combines the 2D coupling constant $g_{\mathrm{2D}}$ with the number of particles $N$.

It is worth emphasizing that our 2D model can be experimentally realized as a 3D disk-shaped BEC in a highly anisotropic trap [Goerlitz2001, burger-2001, rychtarik-2004], with the confinement length $a_{z}\equiv\sqrt{\hbar/(m\omega_{z})}$ along the $z$-direction being much smaller than those along the $x$- and $y$-direction [petrov-2000, salasnich-2002, mateo-2008], where $\omega_{z}$ is the trap frequency along the $z$-axis. As a result, the 2D coupling constant $g_{\mathrm{2D}}$ is determined by the 3D one $g_{\mathrm{3D}}$ as $g_{\mathrm{2D}}=g_{\mathrm{3D}}/(\sqrt{2\pi}a_{z})$, where $g_{\mathrm{3D}}=4\pi\hbar^{2}a_{\mathrm{s}}/m$ and $a_{\mathrm{s}}$ is the $s$-wave scattering length. Using these relations, the results based on the full 3D simulations have been shown to be in good agreement with the results of the 2D simulations.

3.2.3 Time evolution within a rotating trap

The starting point for all simulations is the ground-state wave function $\psi_{\mathrm{\Lambda}}(0,\bm{\xi})$ of the BEC obtained by imaginary time propagation of Eq. (13) with $\Lambda(\tau)=\mathds{1}$ in combination with the Newton method. For the time evolution itself, the affinely transformed GP equation (13) is solved with an implicit Adams-Bashforth-Moulton multi-step algorithm [book-numerical-recipes-2007].

The first scenario in which the accuracy of the time-dependent TF approximation is checked numerically is the rotating trap. In this case the time-dependent trap matrix ${\Omega^{2}(\tau)=O(\tau)D\,O^{\rm T}(\tau)}$ contains the diagonal matrix

[TABLE]

with the time-independent trap frequencies and an orthogonal matrix

[TABLE]

describing a rotation around the $z$-axis with the time-dependent angle $\varphi(\tau)$.

Triggered by experiments [Hechenblaikner2002, Muentinga2013], the rotation rate $\dot{\varphi}(\tau)$ of the trap has been chosen to have a sigmoidal ramp for $0\leq\tau\leq T_{\mathrm{end}}$, as depicted in Fig. 2, and to keep equal to the rotation rate $\dot{\varphi}_{\mathrm{end}}$ for $\tau\geq T_{\mathrm{end}}$. We assume that the characteristic time $T_{\mathrm{end}}$ is large enough, that is $\omega_{x}T_{\mathrm{end}}\gg 1$, to have an adiabatically slow increase of the rotation rate in order to avoid collective excitations of the BEC. In addition, the final rotation rate $\dot{\varphi}_{\mathrm{end}}$ is kept significantly smaller than the trap frequency $\omega_{x}$ in order to avoid the generation of vortices. For the rest of this section we introduce the dimensionless characteristic time $\bar{T}_{\mathrm{end}}\equiv\omega_{x}T_{\mathrm{end}}/(2\pi)$ and use $\bar{T}_{\mathrm{end}}=15$ for all numerical simulations.

Fig. 3 shows the time dependence of the Bures distance $B(\tau)=B\left(\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi}),\psi_{\mathrm{\Lambda}}(0,\bm{\xi})\right)$, Eq. (48), between the time evolved state $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$, being the numerical solution of Eq. (13), and the initial state $\psi_{\mathrm{\Lambda}}(0,\bm{\xi})$ in adapted coordinates. The Bures distance $B(\tau)$ exhibits oscillations over the whole simulation time $\tau$ and the local maxima of this oscillation grow as long as the rotation rate increases until they reach their maximal value $B_{\mathrm{rot}}$ at the characteristic time $\bar{T}_{\mathrm{end}}$ and stay constant from that point on. Since the magnitude of the Bures distance is very small, $B_{\mathrm{rot}}<0.07$, the time evolved state does not differ significantly from the initial state.

For the parameters used in Fig. 3 the Bures distance oscillates with the frequency $\omega_{\mathrm{Bures}}=1.8\,\dot{\varphi}_{\mathrm{end}}$ for $\tau>T_{\mathrm{end}}$. Since the Bures distance, Eq. (48), only measures the absolute value of the overlap between two wave functions, the time evolved state $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ actually undergoes collective oscillations with the frequency $0.5\,\omega_{\mathrm{Bures}}$. In general, the frequency of the collective oscillations caused by the rotation of the trap depends on the rotation rate $\dot{\varphi}(\tau)$ as well as the anisotropy factor $\epsilon$ of the trap.

In Fig. 4 we contrast the dynamics of a rotating BEC in original and adapted coordinates while the whole time evolution is shown in movie 1, which is available in the online supplementary material. Fig. 4a depicts the two-dimensional (non-integrated) density distributions $|\psi(t,\bi{x})|^{2}$ of a BEC in the original coordinates at different evolution times, whereas the first column of Fig. 4b displays the corresponding density distributions $|\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})|^{2}$ in the adapted coordinates. In contrast to the density distribution $|\psi(t,\bi{x})|^{2}$, which follows the clockwise rotation induced by the rotating trap, the density distribution $|\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})|^{2}$ in adapted coordinates does not show any visible changes. To highlight the residual dynamics, we present in the second column of Fig. 4b the difference $||\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})|-|\psi_{\mathrm{\Lambda}}(0,\bm{\xi})||^{2}$ of the absolute values of the wave functions, which depends only on local effects of the density distributions and not on their phases. This difference has non-zero values only at the edge of the BEC forming a ring-like structure, which rotates clockwise in accordance with the density distribution $|\psi(t,\bi{x})|^{2}$ presented in Fig. 4a. However, the magnitude of this residual dynamics is very small compared to the magnitude of the density distribution $|\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})|^{2}$ and it oscillates in accordance with the Bures distance shown in Fig. 3.

In summary, we have proven that the time-dependent TF approximation applied to the affinely transformed GP equation (13) provides an accurate description of the dynamics of a rotating BEC. Only at the very edge of the BEC a minor amount of residual dynamics occurs that goes beyond the time-dependent TF approximation. In subsection 3.2.5 we study how the magnitude of this residual dynamics depends on the different setup parameters.

3.2.4 Free expansion after switching off a rotating trap

As the second scenario to verify the time-dependent TF approximation, we have studied the subsequent free expansion of an initially rotating BEC. Here, the setup is the same as in the first scenario with the exception that the trap is switched off at $\tau=T_{\mathrm{off}}>T_{\mathrm{end}}$, that is shortly after the final rotation rate $\dot{\varphi}_{\mathrm{end}}$ is reached. Since the timing of the switch off has a strong influence on the final value $B_{\mathrm{free}}$ of the Bures distance, we have varied $T_{\mathrm{off}}$ within the first period of the Bures distance after $T_{\mathrm{end}}$ for each set of parameters to find the maximal value of $B_{\mathrm{free}}$.

The time dependence of the Bures distance for this case is depicted in Fig. 5. During the slow initial ramp, shown in Fig. 2, the Bures distance displays the collective oscillations discussed in subsection 3.2.3. After switching off the trap at $\tau=T_{\mathrm{off}}$, these oscillations stop and the Bures distance approaches its final value $B_{\mathrm{free}}$. This behavior is a consequence of the fact that for a freely expanding BEC the adaptive matrix grows linearly for large times, as discussed in section 5.1, and thus the right-hand side of the affinely transformed GP equation (13) vanishes asymptotically.

The two-dimensional (non-integrated) density distributions in the original coordinates $|\psi(t,\bi{x})|^{2}$, Fig. 6a, as well as in the adapted coordinates $|\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})|^{2}$, Fig. 6b, are displayed at different times indicated by the black dots in Fig. 5. The entire time evolution of these density distributions is presented in movie 2, which is available in the online supplementary material. Whereas the size of the BEC in the original coordinates quickly grows after the trap is switched off, the rotation of the BEC stops completely after a further rotation of about $\pi/2$, in accordance with the irrotationality condition, as well as the energy and momentum conservation. It is worth mentioning that our results are in good agreement with the theoretical consideration [Edwards2002] and the experimental data [Hechenblaikner2002].

However, in adapted coordinates, Fig. 6b, neither an increase of the size of the BEC nor a significant rotation of the density distribution can be observed. Only by looking at the residual dynamics visualized again by the difference $||\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})|-|\psi_{\mathrm{\Lambda}}(0,\bm{\xi})||^{2}$, minor changes of $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ become visible at the edge of the condensate. Thus, the time-dependent TF approximation used for the affine approach is an excellent tool to predict the free expansion of an initially rotating BEC as well.

3.2.5 Dependence of the Bures distance on the setup parameters

In extension of the two previous scenarios, we now provide a broader study of the accuracy of the time-dependent TF approximation. In particular, we have solved numerically Eq. (13) for a rotating BEC with and without subsequent free expansion for many different combinations of the interaction strength $\bar{g}N$, the anisotropy factor $\epsilon$ and the final rotation rate $\dot{\varphi}_{\mathrm{end}}$, thereby using $\bar{T}_{\mathrm{end}}=15$ as before. The results for the Bures distances $B_{\mathrm{rot}}$ and $B_{\mathrm{free}}$ are presented in Fig. 7.

As displayed in Fig. 7a the Bures distance $B_{\mathrm{rot}}$ decreases with an increase of the interaction strength $\bar{g}N$ and increases for larger values of the anisotropy factor $\epsilon$ and the maximal rotation rate $\dot{\varphi}_{\mathrm{end}}$. These observations can be explained as follows.

(i) With an increase of $\bar{g}N$, the interaction energy term in Eq. (13) becomes more dominant compared to the kinetic energy term, resulting in an improved accuracy of the time-dependent TF approximation. For this reason, the Bures distance $B_{\mathrm{rot}}$ decreases as the interaction strength $\bar{g}N$ grows.

(ii) Since a spherically symmetric trap with $\epsilon=1$ cannot transfer any angular momentum to a BEC, the wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ will remain in the initial state $\psi_{\mathrm{\Lambda}}(0,\bm{\xi})$, giving rise to $B_{\mathrm{rot}}=0$. However, for $\epsilon>1$ the rotation of the trap affects the dynamics of the BEC. The Bures distance $B_{\mathrm{rot}}$ increases with an increase of $\epsilon$, since the kinetic energy term in Eq. (13) depends on the spatial derivatives of the wave function and therefore increases for a stronger confinement along one axis. Thus, the time-dependent TF approximation works better for a system with small anisotropies.

(iii) For larger values of the maximal rotation rate $\dot{\varphi}_{\mathrm{end}}$, the total energy transfered from the trap to the BEC increases, raising the chance to excite collective modes of the condensate. However, these modes are not taken into account within the framework of the TF approximation. As a result, the amount of residual dynamics grows with the increase of the maximal rotation rate $\dot{\varphi}_{\mathrm{end}}$. The quadrupole mode is well known to be excited when $\dot{\varphi}_{\mathrm{end}}$ approaches the critical value $\omega_{\mathrm{c}}=0.71\,\omega_{x}$ [madison-2001, hodby-2001] and the time-dependent TF approximation breaks down completely in this regime. We have verified this fact with our numerical simulations and for the results presented here, we have always made sure to stay well below the threshold for vortex creation.

Since, the dependencies of $B_{\mathrm{rot}}$ and $B_{\mathrm{free}}$ on the simulation parameters are very similar to each other, as shown in Fig. 7, we conclude here that the subsequent free expansion of an initially rotating BEC does not substantially amplify the residual time evolution in the adapted coordinates.

In summary, we find that the accuracy of the time-dependent TF approximation depends on various parameters with the most important one being the interaction strength $\bar{g}N$. For values $\bar{g}N\geq 100$, the transformed wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ is practically constant in time and the time-dependent TF approximation holds true. In most experimental setups, one typically aims at a large number of particles in a BEC to improve the signal-to-noise ratio. For this reason, the TF approximation provides a valuable tool to describe the dynamics of a large class of rotating BECs.

4 Constants of motion

In the preceding section we have confirmed the accuracy of our method based on the time-dependent TF approximation by dedicated numerical simulations of the affinely transformed GP equation (13). However, the reliability of such numerical simulations themselves must also be guaranteed. One possible way to test it is based on simulations of suitable scenarios that possess certain underlying symmetries which by Noether’s theorem give rise to specific constants of motion. In this chapter, we focus on the time-translational invariance and the rotational symmetry and discuss the corresponding energy and angular momentum conservation for both, the matrix differential equation (14) and the GP equation (1).

In this section, we first show how the matrix differential equation (14) can be rewritten in terms of a Hamiltonian formalism. Second, we analyze the center-of-mass motion of the BEC with regard to a generalized version of the Ehrenfest theorem [Bodurov1998] and third, we present expressions for specific constants of motion of the GP equation (1) valid within the time-dependent TF approximation.

4.1 Hamiltonian formalism for the matrix differential equation

Here we establish a connection between the matrix differential equation (14) and the corresponding Hamiltonian formalism. For this purpose, we first perform a transformation of the adaptive matrix and the time coordinate to bring the matrix differential equation (14) into a specific form that does no longer contain the initial trap matrix $\Omega^{2}(0)$. We call the resulting equation the canonical form of the matrix differential equation. We then introduce the Hamiltonian and verify that the corresponding Hamilton equations of motion are equivalent to the canonical form of the matrix differential equation. With the help of this Hamiltonian we finally obtain two constants of motion for the matrix differential equation, which we later on relate to the conservation of the total energy and the angular momentum of the BEC.

4.1.1 Canonical form of the matrix differential equation

We start by defining the orthogonal matrix $O$ and the diagonal matrix $D$ via the diagonalization of the initial trap matrix $\Omega^{2}(0)$ according to

[TABLE]

Moreover, we introduce a new time scale with the help of the quantity

[TABLE]

Based on these new quantities, we now consider the transformation

[TABLE]

to the new adaptive matrix $\tilde{\Lambda}$ and the dimensionless time $\tilde{t}$. The canonical form of the matrix differential equation is simply obtained by inserting the transformation (4.1.1) into the matrix differential equation (14) thereby taking into account the relation $\det(\Lambda)=\det(\tilde{\Lambda})$. It describes the time evolution of the new adaptive matrix $\tilde{\Lambda}(\tilde{t})$ and reads

[TABLE]

where we have introduced the transformed trap matrix

[TABLE]

The initial conditions (15) accordingly transform to

[TABLE]

whereas the irrotationality condition (16) preserves its form under the transformation (4.1.1)

[TABLE]

Next we show that the canonical form of the matrix differential equation (54) can be embedded in a Hamiltonian formalism.

4.1.2 Hamiltonian and equations of motion

We start by introducing the momenta

[TABLE]

as conjugate variables to the matrix elements $\tilde{\Lambda}_{\alpha\beta}$ together with the Hamiltonian

[TABLE]

The Hamilton equations of motion

[TABLE]

are mathematically equivalent to the canonical form of the matrix differential equation (54). Here, we have introduced the corresponding Poisson bracket

[TABLE]

of two functions $f$ and $g$ that depend on the matrix elements $\tilde{\Lambda}_{\alpha\beta}$ and $\tilde{\Pi}_{\alpha\beta}$ and made use of the relation

[TABLE]

When we accordingly define the “momenta” $\Pi\equiv\mathrm{d}\Lambda/\mathrm{d}t$, Eq. (4.1.1) implies the following transformation of the conjugate variables

[TABLE]

We emphasize that Eq. (4.1.2) is not a canonical transformation and that a corresponding Hamiltonian formalism for the original matrices $\Lambda$ and $\Pi$ does not exist, as can be seen by inserting the transformation (4.1.2) into the expression for the Poisson-bracket (61). As a result, we arrive at

[TABLE]

which clearly shows that the structure of the Poisson bracket (61) is not preserved by the transformation (4.1.2). Thus, there exist no corresponding Hamilton equations for the original matrix differential equation (14).

4.1.3 Energy conservation

When we consider a time-independent trapping potential with the constant trap matrix $\tilde{\Omega}^{2}(\tilde{t})=\tilde{\Omega}_{\ast}^{2}$, the energy is a constant of motion given by the Hamiltonian $H(\tilde{\Lambda},\tilde{\Pi})$, Eq. (59). In this spirit, we call the quantity

[TABLE]

the total energy associated with the matrix differential equation in canonical form also for a time-dependent trap matrix $\tilde{\Omega}^{2}(\tilde{t})$.

Using Eq. (58) together with the inverse of the transformation (4.1.1) in Eq. (65), the total energy associated with the matrix differential equation $E_{\Lambda}(t)=\tilde{E}(\tilde{\Lambda},\tilde{\Pi};\tilde{t})$ reads in terms of the original adaptive matrix

[TABLE]

where $\Omega^{2}(t)=\alpha^{2}\,O\,\tilde{\Omega}^{2}(\tilde{t})\,O^{\rm T}$ and $\Omega^{-2}(0)=O\,D^{-1}\,O^{\rm T}$. The fact that the total energy (66) is indeed a constant of motion if $\Omega^{2}(t)=\Omega_{\ast}^{2}$ can be easily verified by taking the time derive of Eq, (66) on both sides and inserting the matrix differential equation (14) on the right hand side of the resulting equation.

4.1.4 Angular momentum conservation

For an isotropic harmonic potential with corresponding trap matrix $\tilde{\Omega}^{2}(\tilde{t})=\tilde{\omega}^{2}(\tilde{t})\,\mathds{1}$, the angular momentum of a BEC is preserved. In order to see how this conservation law is related to the matrix differential equation (14), we study the quantity

[TABLE]

which we denote as angular momentum matrix associated with the matrix differential equation in canonical form. The angular momentum matrix (67) can be found with the help of Noethers theorem and the irrotationality condition (57) which reads in terms of both conjugate variables $\tilde{\Lambda}^{\rm T}\,\tilde{\Pi}-\tilde{\Pi}^{\rm T}\,\tilde{\Lambda}=0$.

We now prove that this angular momentum matrix is preserved for an isotropic trapping potential. For this purpose, we take the total derivative of $\tilde{L}$ with respect to the time $\tilde{t}$ and obtain

[TABLE]

where the Poisson bracket $\{H,\tilde{L}\}_{(\tilde{\Lambda},\tilde{\Pi})}$ is defined by Eq. (61). Since the angular momentum matrix $\tilde{L}$ does not explicitly depend on the time $\tilde{t}$, the partial derivative $\partial\tilde{L}/\partial\tilde{t}$ vanishes. With the Hamiltonian (59), the evaluation of the Poisson bracket finally yields

[TABLE]

Thus, the angular momentum matrix $\tilde{L}$ is a constant of motion if the matrix $\tilde{\Lambda}\,\tilde{\Lambda}^{\rm T}\,\tilde{\Omega}^{2}(\bar{t})$ is symmetric, which is indeed the case for an isotropic trap with $\tilde{\Omega}^{2}(\tilde{t})=\tilde{\omega}^{2}(\tilde{t})\,\mathds{1}$.

When we insert the transformation (4.1.2) into the definition (67) and use the transformation law $L_{\Lambda}(t)=O\,\tilde{L}(\tilde{\Lambda},\tilde{\Pi})\,O^{\rm T}$, we obtain the angular momentum matrix in terms of the original adaptive matrix

[TABLE]

We emphasize that the angular momentum matrix is antisymmetric and explicitly depends on the initial trapping potential.

4.2 Constants of motion of the Gross-Pitaevskii equation

Our efficient description of the time evolution of a rotating BEC within the TF regime, Eq. (37) rests upon two different, but related approximations: (i) the initial macroscopic wave function $\psi\left(0,\bi{x}\right)$ is assumed to be given by the ground state in the TF approximation, Eq. (29), and (ii) the kinetic energy term in the affinely transformed GP equation (13) is supposed to be negligible in comparison to the potential and interaction energies involved. According to this second assumption known as the time-dependent TF approximation, the affinely transformed wave function, being the solution of Eq. (13), does not display any essential dynamics

[TABLE]

In section 2.2.2 we used a combination of both approximations, (i) and (ii), to derive the efficient description of the BEC dynamics within the TF regime given by Eq. (37).

However, for BECs that are not deep within the TF regime and whose ground state differs significantly from Eq. (29), the assumption (71) can still provide very accurate results for some scenarios such as freely expanding BECs. For this reason, we first seek an efficient description that does not involve approximation (i), but includes the exact ground state $\psi(0,\bi{x})$ of a BEC or its numerically determined approximation. Based on this slightly more general description, we derive approximate expressions for the total energy and the angular momentum of a BEC that allow further tests of assumption (71) by analyzing possible time dependencies of the constants of motion that should not occur due to the chosen symmetry of the GP equation. Eventually, we also include approximation (i) and relate the resulting expressions for the total energy and the angular momentum based on the TF ground state (29) to the corresponding constants of motion associated with the matrix differential equation.

4.2.1 Effective description of the dynamics based solely on the time-dependent Thomas-Fermi approximation

In order to determine the slightly generalized efficient description of the BEC dynamics, we first relate the initial state in adapted coordinates $\psi_{\mathrm{\Lambda}}(0,\bm{\xi})$ with the original ground state $\psi(0,\bi{x})$ by evaluating Eq. (20) at $t=0$ giving rise to

[TABLE]

Next, we obtain $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ simply by combining Eqs. (71) and (72). We insert the result into Eq. (17) and thereby use the substitution $\bm{\xi}=\Lambda^{-1}(t)\left[\bi{x}-\bm{\mathcal{R}}(t)\right]$. Hence, we find as efficient description of the BEC dynamics that is solely based on the time-dependent TF approximation

[TABLE]

where we have introduced the phase

[TABLE]

Note that by insertion of $\psi(0,\bi{x})=\phi_{\mathrm{TF}}(\bi{x})$ together with Eq. (29) into Eq. (73), we consequentially arrive at the expression (37).

4.2.2 Generalized Ehrenfest Theorem for a Bose-Einstein condensate in a harmonic trap

According to the generalized Ehrenfest-Theorem [Bodurov1998], the expectation values

[TABLE]

of the position $\hat{\bi}{x}=\bi{x}$ and the momentum $\hat{\bi}{p}\equiv-\rmi\hbar\bm{\nabla}_{\bi{x}}$ operators satisfy the classical equations of motion

[TABLE]

The expectation values are evaluated with the full solution $\psi(t,\bi{x})$ of the GP equation (1) for the harmonic potential, Eq. (3). Since the Eqs. (4.2.2) coincide with the Eqs. (2.1.2) for the time-dependent vectors $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$, we can identify the latter ones with the center-of-mass position $\left\langle\hat{\bi}{x}\right\rangle_{\psi(t)}$ and the average momentum $\left\langle\hat{\bi}{p}\right\rangle_{\psi(t)}$ by simply choosing the initial conditions according to

[TABLE]

Here, we have introduced the shorthand notation of the exact initial state $\psi_{0}\equiv\psi(0,\bi{x})$. We emphasize that the parameters $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$ introduced in section 2.1.2, have strictly speaking no physical meaning unless their initial conditions are fixed in accordance with Eq. (78).

Next, we discuss if the above assignment

[TABLE]

is also consistent with our efficient description, Eq. (73), based on the time-dependent TF approximation. For this purpose, we insert Eq. (73) into the Eqs. (75) and (76) and make use of the substitution

[TABLE]

to arrive at the expectation values

[TABLE]

Thus, using $\bm{\mathcal{R}}(0)$ and $\bm{\mathcal{P}}(0)$ according to Eq. (78), the last two expressions simply reduce to Eqs. (79), which proves that our efficient description based on the time-dependent TF approximation, Eq. (73), satisfies the generalized Ehrenfest theorem.

4.2.3 Energy of a Bose-Einstein condensate

Here we use the approximate BEC wave function $\psi(t,\bi{x})$ given by Eq. (73) to derive an explicit expression for the total energy

[TABLE]

of the BEC within the time-dependent TF approximation, where the kinetic, the potential and the interaction energy are defined by

[TABLE]

In B.1 we evaluate each of these energy terms and obtain as the energy per particle

[TABLE]

where we have assumed that the initial wave function $\psi(0,\bi{x})$ is given by

[TABLE]

Here, $\phi_{\mathrm{D}}(\bi{x})$ is the real-valued ground state of the BEC in the comoving frame of the trap as discussed in section 2.1.6 and the phase $\rmi\bm{\mathcal{P}}(0)\,\bi{x}/\hbar$ accounts for the initial momentum $\bm{\mathcal{P}}(0)$ of the BEC in the inertial frame of reference $\bi{x}$, analogously to the TF ground state (29). The first two terms in Eq. (4.2.3) display the kinetic and the potential energy of the center of mass of the BEC, respectively, whereby we have used $V(t,\bm{\mathcal{R}}(t))$ as shorthand notation for the harmonic potential (3). The third term corresponds to the quantum pressure and its expectation value is evaluated with respect to the absolute value $|\psi(0,\bi{x})|$ of the initial wave function. The forth term describes the mean-field interaction energy. The last term is determined by the explicit form of the initial wave function $\psi(0,\bi{x})$ and represent the inner dynamics of the condensate associated with the deformation of the BEC which corresponds to the evolution of the adaptive matrix $\Lambda(t)$.

When the initial wave function $\psi(0,\bi{x})$ is given by the TF ground state $\phi_{\mathrm{TF}}(\bi{x})$, Eq. (29), the energy per particle reads

[TABLE]

as shown in B.2. Here $E_{\Lambda}(t)$ denotes the total energy associated with the matrix differential equation given by Eq. (66) whereas $\mu_{\mathrm{TF}}$ is the chemical potential within the TF approximation, Eq. (30).

By applying Eq. (15) to the energy per particle (89) and assuming without loss of generality that $V(0,\bm{\mathcal{R}}(0))=0$, the initial energy per particle for a BEC in the TF ground state $\phi_{\mathrm{TF}}(\bi{x})$ reads

[TABLE]

and is determined solely by the initial momentum $\bm{\mathcal{P}}(0)$ and the chemical potential $\mu_{\mathrm{TF}}$.

4.2.4 Angular momentum of a Bose-Einstein condensate

In analogy to the previous subsection, we now discuss the angular momentum of a BEC within the time-dependent TF approximation. Indeed, as derived in B.3 the angular momentum operator with respect to the center-of-mass coordinates $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$ is defined by

[TABLE]

In B.4 we show that the expectation value of the angular momentum operator, Eq. (91), evaluated with respect to the approximate wave function (73), is given by the relation

[TABLE]

where we have used the initial conditions (78). Here, the first vector of the cross product represents the position, while the second vector corresponds to the momentum associated with the inner dynamics of the BEC. For the sake of a concise mathematical notation, we map the angular momentum operator to the so-called angular momentum matrix operator via $\hat{L}_{kl}=\epsilon_{jkl}\hat{L}_{j}$, where $\varepsilon_{jkl}$ denotes the Levi-Civita symbol. Applied to Eq. (92) this mapping yields the following expression for the matrix operator

[TABLE]

This matrix notation of the angular momentum is especially useful if the initial wave function $\psi(0,\bi{x})$ is approximated by the TF ground state $\phi_{\mathrm{TF}}(\bi{x})$, Eq. (29), since in that case Eq. (4.2.4) reduces to

[TABLE]

as discussed in B.5 with the angular momentum matrix $L_{\Lambda}(t)$ given by Eq. (70).

In conclusion, the deep connection between the constants of motion of the GP equation, Eq. (1), and the matrix differential equation, Eq. (14), illustrated by Eqs. (89) and (94), once more highlights the benefit of the affine approach for the description of the BEC dynamics. Not only grants the solution of Eq. (14) valuable insights into the time evolution of a BEC, but it is also relevant for numerical applications.

5 Special properties of free expanding Bose-Einstein condensates

Experiments dealing with quantum gases usually involve time-of-flight pictures to deduce the state of the system under study. It is well known, and we will briefly recall this statement later, that in the long-time limit the spatial density distribution obtained by time-of-flight measurements is determined by the initial momentum distribution of the quantum gas if the free time evolution is described by the Schrödinger equation. In other words, a non-interacting quantum gas undergoes a ballistic expansion if it is released from the trapping potential. However, it is not obvious that this statement holds true for free expanding BECs governed by the non-linear GP equation. In fact, for BECs a reverse relation applies [PhD_ME], namely that the momentum distribution in the long-time limit is given by the initial spatial density distribution. This relation is not only relevant for the interpretation of time-of-flight pictures of BECs, but it also plays a role for matter-wave interferometry [Kleinert2015], where the visibility of the phase oscillations quickly decreases when the momentum distribution of the atomic clouds exceed the corresponding line width of the laser pulses.

This section starts with the discussion of the analytic solutions of the matrix differential equation (14) for freely expanding BECs in the long-time limit that start from isotropic initial traps. We then recall the long-time behavior for the Schrödinger equation before we turn to the GP equation.

5.1 Analytic solutions of the matrix differential equation

For a pure free expansion, that is $\Omega^{2}(t)=0$ for all times $t>0$, the matrix differential equation (14) possesses analytical solutions that can be used to determine the macroscopic wave function of a BEC in the long-time limit. In order to derive these solutions we assume without loss of generality that the coordinate system $\bi{x}$ is aligned with the principal axes of the trap at the initial time $t=0$. As a consequence, the adaptive matrix $\Lambda(t)$ is diagonal during the free expansion, that is $\Lambda(t)={\rm diag}\left[\lambda_{1}(t),\ldots,\lambda_{d}(t)\right]$ for $t>0$, and the scaling parameters $\lambda_{i}$ with $i=1,2,...,d$ are determined by the system of non-linear differential equations

[TABLE]

which are a special case of equations (14) or (22) for free evolution. Here the frequencies $\omega_{i}(0)$ denote the initial trap frequencies along the principal axes. The initial conditions (15) reduce to

[TABLE]

accordingly.

5.1.1 Long-time limit

Since the scaling parameters $\lambda_{i}(t)$ determine the size of the condensate, they grow as the condensate expands freely. Hence, the right hand side of Eq. (95) vanishes in the long-time limit and the scaling parameters approach the linear dependence

[TABLE]

as $t\rightarrow\infty$. In general, the constants $a_{i}$ and $b_{i}$ are determined numerically. However, in the case of an isotropic initial trap Eq. (95) can be solved analytically also yielding the values for $a_{i}$ and $b_{i}$ as shown in the next subsection.

5.1.2 Isotropic case

For an isotropic initial trap, that is $\omega_{i}(0)=\omega_{0}$, the scaling parameters $\lambda_{i}(t)$ are given by the single scaling parameter $\lambda(t)$ for all $i=1,2,...,d$ and the system of non-linear differential equations, Eq. (95), decouples. We show in C that the function $\lambda(t)$ is a solution of the integral equation

[TABLE]

Depending on the dimensionality $d$, Eq. (98) gives rise to explicit, $d=2$, or implicit, $d=1$ and $d=3$, solutions for the scaling parameter $\lambda(t)$, which are discussed in the next paragraphs.

The case $d=1$

The left-hand side of Eq. (98) gives rise to the transcendental equation

[TABLE]

which can also be written as

[TABLE]

It is clear that for long times, that is $\omega_{0}t\gg 1$, the function $\lambda\cong\sqrt{2}\,\omega_{0}t$, since the denominator in the right-hand side of Eq. (99) approaches unity as $\lambda\rightarrow\infty$. Taking into account the terms in the denominator scaling as $\lambda^{-1}$ and $\ln(\lambda)\lambda^{-1}$, we finally arrive at the approximate solution

[TABLE]

for $\omega_{0}t\gg 1$. Here we have neglected all terms, that vanish in the long-time limit.

The case $d=2$

In the case of two dimensions, the integral in the left-hand side of Eq. (98) yields $\sqrt{\lambda^{2}(t)-1}$, leading to the exact solution

[TABLE]

valid for all $t\geq 0$.

The case $d=3$

For a three dimensional BEC the integral in the left-hand side of Eq. (98) can be expressed in terms of the Gaussian hypergeometric function ${}_{2}F_{1}(a,b,c;z)$ leading to the implicit equation

[TABLE]

for $\lambda(t)$, where $\Gamma(x)$ denotes the Gamma function. Since the hypergeometric function ${}_{2}F_{1}(-1/3,1/2,2/3;1/\lambda^{3})$ approaches unity as $\lambda\rightarrow\infty$, the solution $\lambda(t)$ of Eq. (102) in the long time limit, that is $\omega_{0}t\gg 1$, is given by

[TABLE]

When we compare the relations (100), (101), and (103) with Eq. (97), we find that the coefficient $b_{i}=\sqrt{2/d}\,\omega_{0}$.

5.2 Connection between the long-time behavior of the wave function and the initial state

Here we show how the long-time behavior of the macroscopic wave function of a BEC can be estimated in the case of free expansion, although the dynamics is governed by the non-linear GP equation. Since the approximate analytic result for BECs contrasts the long-time solution of non-interacting wave functions, we first have a look at the Schrödinger equation before discussing the GP equation.

5.2.1 Schrödinger equation

The Schrödinger equation describing the free evolution of the non-interacting wave function $\psi(t,\bi{x})$ in a $d$-dimensional space reads

[TABLE]

and has the formal solution

[TABLE]

in terms of the initial wave function $\psi(0,\bi{x})$. When we use the Fourier transformation

[TABLE]

of the wave function $\psi(0,\bi{x})$, we can cast Eq. (105) into the form

[TABLE]

For long times, that is in the limit $t\rightarrow\infty$, the exponential function inside the integral on the right-hand side of Eq. 107 oscillates very rapidly and only the values of $\bi{p}$ in a narrow region around $\bi{p}=m\,\bi{x}/t$ contribute to the integral. Hence, when we assume that the momentum distribution $\tilde{\psi}(0,\bi{p})$ is a smooth function of $\bi{p}$ in this region, we can evaluate $\tilde{\psi}(0,\bi{p})$ at $\bi{p}=m\,\bi{x}/t$ and perform the remaining integral, leading to the relation

[TABLE]

Consequently, as is well known, the spatial density distribution in the long time limit

[TABLE]

is determined by the shape of the initial momentum distribution $\tilde{\psi}(0,\bi{p})$.

5.2.2 Gross-Pitaevskii equation

For free expanding BECs a similar relation can be established based on the affine approach. In order to compare the results derived for the Schrödinger equation and the GP equation, respectively, it is suitable to consider the wave function $\psi_{\mathrm{D}}(t,\bi{x})$ that describes the inner dynamics of the BEC where the center-of-mass motion has already been eliminated. As shown in section 2.1, the time evolution of $\psi_{\mathrm{D}}(t,\bi{x})$ is governed by the transformed GP equation (8) which is related to the affinely transformed wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ by Eq. (10). If the initial wave function $\psi_{\mathrm{D}}(0,\bi{x})$ coincides with the TF ground state, Eq. (32), the time-dependent TF approximation (34) holds true during the free expansion and we arrive at the solution

[TABLE]

where we have used the definition (12) of the symmetric matrix $A(t)$.

In the long-time limit and for free expansion the adaptive matrix $\Lambda(t)$ grows linearly, Eq. (97), that is $\Lambda(t)\cong B\,t$ as $t\to\infty$ with the diagonal matrix $B={\rm diag}(b_{1},\ldots,b_{d})$. When we apply this limit to expression (110), we obtain

[TABLE]

As a result, the spatial density distribution in the long time limit

[TABLE]

is determined by the initial spatial density distribution and the time-dependent scaling. Obviously, Eq. (112) is a direct consequence of the time-dependent TF approximation on which this derivation is based and stands in contrast to relation (109) derived for the Schrödinger equation.

In the next step we obtain a similar relation for the momentum distribution of a free expanding BEC in the long-time limit. In order to do so, we insert Eq. (111) into Eq. (106) for the Fourier transform and arrive at the expression

[TABLE]

where we have introduced the new integration variable $\bi{x}^{\prime}\equiv\bi{x}-\bi{p}t/m$. For long times $t\rightarrow\infty$, the argument of the initial wave function $\psi_{\mathrm{D}}$ gets practically independent of $\bi{x}^{\prime}$, which allows to perform the remaining integration and to obtain

[TABLE]

Hence, in the case of the GP equation for free expansion, the momentum distribution in the long-time limit

[TABLE]

is given by the initial spatial density distribution, while for the Schrödinger equation it is vice versa, Eq. (109). Moreover, according to Eq. (115) the momentum distribution $|\tilde{\psi}_{\mathrm{D}}(t,\bi{p})|^{2}$ of a free expanding BEC becomes time-independent in the long-time limit. However, we emphasize that for the derivation of Eq. (115) we have required the validity of the time-dependent TF approximation during the free expansion. As a consequence, for a proper interpretation of time-of-flight pictures of BECs the interactions play an important role, especially deep within the TF regime.

In addition, we point out that the methods presented in this section have been successfully applied by [Roura2014] to study the influence of gravity gradients on the properties of an atom interferometer.

6 Conclusion and Outlook

In this chapter we have introduced a natural generalization of the scaling approach for time-dependent rotating traps, the so-called affine approach. In contrast to the hydrodynamical approach, we have carried out this generalization directly to the GP equation by first eliminating the center-of-mass motion and by subsequently applying a linear mapping to account for the main part of the remaining internal dynamics of a BEC. With the help of this affine approach we have established an efficient analytic description of the BEC dynamics based on two main assumptions: (i) the external potential in the vicinity of the BEC is at most quadratic and (ii) the time-dependent TF approximation holds true during the whole time evolution, which manifests itself in the fact that the density distribution of the condensate preserves its parabolic shape.

We have verified the accuracy of this approach by performing full numerical simulations of the GP equation for two scenarios: (i) a pure rotation of the trap, and (ii) a free expansion of an initially rotating BEC. In both cases we have found an excellent agreement of the numerical results with our approximate analytical solution for many different combinations of the interaction strength, the anisotropy factor and the final rotation rate of the trap. Moreover, the affine transformation can be successfully used to improve the performance of numerical simulations for a variety of experimentally relevant scenarios by solving the transformed GP equation rather than the original one.

By employing an additional transformation, we have introduced the canonical form of the matrix differential equation that governs the time evolution of the adaptive matrix and presented its corresponding Hamiltonian formalism. We have furthermore analyzed the relation between the constants of motion of this matrix differential equation and the constants of motion of the GP equation in connection with the conservation of the total energy and the angular momentum of a BEC. For both cases, we have derived explicit expressions valid within the TF regime.

In agreement with our affine approach, the spatial density distribution of a freely expanding BEC in the long-time limit does not display the initial momentum distribution in general, unlike a freely expanding atomic cloud that evolves according to the Schrödinger equation. Instead, it resembles the affinely transformed initial spatial density distribution. However, there exists a converse relation between the spatial and momentum distributions for a BEC within the TF regime, that is the momentum distribution of a free expanding BEC in the long-time limit is indeed given by its initial spatial density distribution. These relations are a direct consequence of the nonlinearity of the GP equation and stand in contrast to the well-known ballistic expansion that wave functions governed by the free Schrödinger equation show. Hence, in order to properly interpret time-of-flight pictures of BECs, the mean-field interaction has to be taken into account. This observation is of special importance for the TF regime where the interactions play a dominant role.

In summary, scaling solutions like the affine approach presented in this chapter are valuable tools for studying the dynamics of strongly interacting quantum gases. Not only do they provide efficient ways to accurately describe their time evolution, but they also allow to unravel and separate the different layers of their dynamics. Similar approaches have been successfully used in a multitude of interesting physical phenomena, such as (i) the Gaussian wave packet, which keeps its form in an at most harmonic potential, (ii) the Airy wave packet [Berry1979, Kajari2010] and (iii) the Bessel beam [Durnin1987, DurninEberly1987], which freely propagate without dispersion.

We emphasize that there are a few scenarios in which one of the two above-mentioned assumptions for our efficient analytic description of the BEC dynamics is not fulfilled. Indeed, the external potential could be of a form incompatible with the harmonic expansion, Eq. (3), as is the case for (i) realistic trapping potentials created by optical dipole traps, optical lattices and magnetic chip traps, where the anharmonic terms play a crucial role especially in the context of delta-kick collimation techniques [ammann-1997, Zeller2016], and (ii) trapping potentials which are not differentiable due to sharp corners like linear potential wells (the irrigation canal) [Gallas1995, Bestle1995] or gravitational traps [Aminoff1993, Nesvizhevsky2002]. Moreover, the TF approximation is well-known to break down when the interaction strength approaches very small values, which can be caused by a low number of atoms, or a decreasing coupling constant. In order to accurately describe BECs in this regime one can either expand the theoretical model beyond the TF approximation [Jamison2011] or resort to numerical methods based on the affinely transformed GP equation. For all these scenarios the affine approach can still offer valuable insights by dedicated numerical simulations that take advantage of the affine transformation to substantially reduce the computational costs of the simulations [PhD_ME].

Finally, it is noteworthy that the affine approach can even be applied to multi-species BECs to study their expansion dynamics both from an analytical and numerical perspective. However, a detailed discussion of this method is beyond the scope of this chapter and will be addressed in a future publication [Meister2017].

Acknowledgments

We would like to thank H. Ahlers, M. Buser, A. Friedrich, N. Gaaloul, E. Giese, N. Harshman, W. Herr, J. Jenewein, S. Kleinert, M. Krutzik, W. Lewoczko-Adamczyk, A. Roura, E. Sadurni, S. T. Seidel, E. M. Rasel, C. Ufrecht, V. Yakovlev, W. Zeller and T. van Zoest for many fruitful discussions and helpful suggestions. M. A. E. thanks the Alexander von Humboldt Stiftung and W. P. S. is grateful to the Hagler Institute for Advanced Study at Texas A&M University for a Faculty Fellowship. This project was generously supported by the German Space Agency (DLR) with funds provided by the Federal Ministry for Economic Affairs and Energy (BMWi) under the grant numbers 50WM0346, 50WM0837, 50WM1136 and 50WM1556.

Appendix A Affine transformation

In this appendix we outline the derivation of the affinely transformed GP equation (13). This is done by first eliminating the center-of-mass motion of the condensate and second, by performing a linear transformation of the coordinates together with a unitary transformation of the wave function. Moreover, we show how to determine the integrated density distributions of a BEC within the TF regime.

A.1 Center-of-mass motion

Our starting point is the GP equation (1)

[TABLE]

with the explicit time-dependent harmonic potential given by Eq. (3). We then rewrite the transformation (4) by making use of the displacement operator in position representation

[TABLE]

and obtain

[TABLE]

The displacement operator satisfies the relations

[TABLE]

for any smooth wave function $\psi_{\mathrm{D}}(t,\bi{x})$ and phase-space operator $f(\bi{x},\hat{\bi}{p})$. Taking advantage of Eqs. (118) - (120), as well as the fact that

[TABLE]

we arrive at the transformed GP equation

[TABLE]

for the wave function $\psi_{\mathrm{D}}$.

In order to eliminate all terms in Eq. (A.1) that depend linearly on $\bi{x}$ and $\hat{\bi}{p}$, as well as arbitrarily on $\bm{\rho}$, we require the parameters $\bm{\mathcal{R}}(t)$ and $\bm{\mathcal{P}}(t)$ to satisfy the classical equations of motion (2.1.2) and the phase $\mathcal{S}_{1}(t)$ to be a solution of the first order differential equation

[TABLE]

Here, we have used the first equation in Eq. (2.1.2) and the definition of the Lagrange function $\mathcal{L}(\bm{\mathcal{R}},\dot{\bm{\mathcal{R}}},t)$, Eq. (6). As a result, we have simplified the transformed GP equation (A.1) to Eq. (8).

A.2 Linear transformation

We continue the derivation of the affinely transformed GP equation by recalling Eq. (9) for the linear transformation of the spatial coordinates and the substitution of the time variable. This mapping implies that the corresponding partial derivatives in Eq. (8) transform as

[TABLE]

As mentioned in subsection 2.1.3, the linear mapping (9) is accompanied by the unitary transformation of the wave function (10) which guarantees that the affinely transformed wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$ is again normalized according to Eq. (2). Thus, by using the Eqs. (9), (10) and (A.2), as well as the identity

[TABLE]

we arrive at the GP equation for the affinely transformed wave function $\psi_{\mathrm{\Lambda}}(\tau,\bm{\xi})$

[TABLE]

Now we first choose $A(\tau)$ to be given by Eq. (12) which serves the purpose of eliminating the last term in Eq. (125) as well as the part that contains the trace in the second term on the right-hand side. Next we insert $\beta(\tau)$ as defined by Eq. (11). In order to simplify the third term on the right-hand side of Eq. (125), we therein replace $A(\tau)$ according to Eq. (12) and require the adaptive matrix $\Lambda(\tau)$ to satisfy the matrix differential equation (14). Finally, by taking advantage of the irrotationality condition (16) whose validity is a direct consequence of Eqs. (14) and (15) as discussed in the next section, we achieve our goal to simplify Eq. (125) to the affinely transformed GP equation (13).

A.3 Irrotationality condition

Within our approach the adaptive matrix $\Lambda(\tau)$ automatically satisfies the irrotationality condition, Eq. (16), since the time derivative of the auxiliary matrix function

[TABLE]

vanishes for all ${\tau\geq 0}$ provided $\Lambda(\tau)$ is a solution of the matrix differential equation (14) with the initial conditions (15). For this reason, the antisymmetric matrix $Z(\tau)=Z(0)=0$ represents $d(d-1)/2$ constants of motion of the second order matrix differential equation (14) and the matrix $A(\tau)$, Eq. (12), is indeed symmetric.

However, as pointed out by [storey00], the irrotationality condition (16) can likewise be derived by using the notion of classical particle trajectories that constitue a velocity field $\bi{v}(t,\bi{x})$, which is irrotational, namely $\bm{\nabla}\times\bi{v}(t,\bi{x})=0$ at any point in space. We briefly sketch this derivation within our approach. Let $\bi{x}=\bi{x}(\bi{x}_{0},t)$ be the classical trajectory of a particle within a BEC, which depends on the initial position $\bi{x}_{0}$ and the time $t$. As we have shown in section 2.2.2 a BEC does not show any time evolution in the adapted coordinates if the time-dependent TF approximation holds true. Hence, the trajectory in adapted coordinates is given by $\bm{\xi}(t)=\bm{\xi}_{0}=\bi{x}_{0}$, keeping in mind that the original and the adapted coordinates coincide at $t=0$. Since the transformation between the two coordinate systems is defined by the linear mapping $\bm{\xi}=\Lambda^{-1}(t)\left[\bi{x}-\bm{\mathcal{R}}(t)\right]$ with the center-of-mass position $\bm{\mathcal{R}}(t)$ and the adaptive matrix $\Lambda(t)$, the classical trajectory is given by

[TABLE]

or rather

[TABLE]

Thus, the velocity field

[TABLE]

is irrotational, that is $\bm{\nabla}\times\bi{v}(t,\bi{x})=0$, if the matrix $\left(\mathrm{d}\Lambda/\mathrm{d}t\right)\Lambda^{-1}$ is symmetric

[TABLE]

which is equivalent to the irrotationality condition (16).

A.4 Integration of the density distribution of a Bose-Einstein condensate within the Thomas-Fermi regime

In this appendix, we discuss the integration of the density distribution $\left|\psi_{\mathrm{TF}}(t,\bi{x})\right|^{2}$ of a $d$-dimensional BEC within the TF regime over $n$ dimensions with $n<d$. We denote by $\bi{x}_{[2]}\in\mathds{R}^{n}$ the vector consisting of the components of $\bi{x}\in\mathds{R}^{d}$ over which the integration is performed. Accordingly, the vector $\bi{x}_{[1]}\in\mathds{R}^{d-n}$ is composed of the remaining coordinates. Thus, the integrated density distribution is defined as

[TABLE]

Without loss of generality, we can rearrange the coordinates in the original density distribution, Eq. (40), such that $\bi{x}^{\rm T}=(\bi{x}_{[1]}^{\rm T},\bi{x}_{[2]}^{\rm T})$, which furthermore translates into ${\bm{\mathcal{R}}^{\rm T}=(\bm{\mathcal{R}}_{[1]}^{\rm T},\!\bm{\mathcal{R}}_{[2]}^{\rm T})}$ for the center-of-mass position and

[TABLE]

for the TF matrix defined by Eq. (39). By using the decomposition of the matrix $\Sigma$ into the Schur complements [HornJohnson1990]

[TABLE]

the inverse of the TF matrix can be written as

[TABLE]

By inserting this expression into Eq. (40), the quadratic form $(\bi{x}-\bm{\mathcal{R}})^{{\rm T}}\Sigma^{-1}(\bi{x}-\bm{\mathcal{R}})$ splits into two terms: $(i)$ a term involving only the free coordinates $\bi{x}_{[1]}$ and $(ii)$ a term depending on both, the coordinates $\bi{x}_{[2]}$ over which the integration is performed and the free coordinates $\bi{x}_{[1]}$. In view of the integration in Eq. (130), we perform a change of the integration variables from $\bi{x}_{[2]}$ to $\bi{y}$ via

[TABLE]

Note, that the inverse square root of the matrix $(\Sigma_{[22]}-\Sigma_{[12]}^{\,{\rm T}}\Sigma_{[11]}^{-1}\Sigma_{[12]})$ is defined in terms of its spectral decomposition. As a result, we obtain for the integrated density distribution (130)

[TABLE]

where we have introduced the auxiliary function

[TABLE]

We perform the integration in Eq. (132) with the help of the $n$-dimensional spherical coordinates

[TABLE]

where $S_{n-1}$ denotes the $(n-1)$-dimensional spherical surface. Moreover, by applying the determinant on both sides of the decomposition (131), we find with the help of Eq. (39) the useful relation

[TABLE]

Taking advantage of the latter, the integrated density distribution Eq. (132) finally reads

[TABLE]

For $n=1$ and $n=2$, this result reduces to Eqs. (42) and (44), respectively.

Appendix B Energy and angular momentum of a Bose-Einstein condensate

This appendix presents the important steps to determine the different contributions to the total energy of a BEC and its angular momentum within the time-dependent TF approximation.

B.1 Energy terms in the time-dependent Thomas-Fermi approximation

We now evaluate the energy terms, Eqs. (84) - (86), by using the approximate wave function, Eq. (73), valid in the time-dependent TF approximation. With the help of Eq. (80) we can rewrite the expectation values in Eqs. (84) - (86) in terms of the expectation value with respect to the initial state $\psi(0,\bi{x})$. Thus, for the kinetic energy, Eq. (84), we obtain

[TABLE]

In the same manner the potential energy given by Eqs. (85) and (3) reads

[TABLE]

Likewise, the dependence on time of the interaction energy Eq. (86) is given by

[TABLE]

We can further simplify Eqs. (134) and (135) for the kinetic and the potential energies by taking into account the initial conditions (78) and by using the identity

[TABLE]

which is valid for any symmetric matrix $M$. Hence, Eqs. (134) and (135) reduce to

[TABLE]

and

[TABLE]

respectively, where $V(t,\bm{\mathcal{R}}(t))$ denotes the harmonic potential (3).

In addition, by making use of the explicit form of the initial wave function, Eq. (88), we can evaluate the second and last term of Eq. (138). We start by inserting Eq. (88) into the expectation value of the second term of Eq. (138) and arrive at the relation

[TABLE]

where the expectation values are now evaluated with respect to the absolute value of the initial state. Since the BEC is initially at rest in the comoving frame of the trap, we obtain $\left\langle\hat{\bi}{p}\right\rangle_{|\psi_{0}|}=0$ and thus the second term of Eq. (138) fully reads

[TABLE]

Next, we apply Eq. (88), to the last term of Eq. (138) and arrive at the expression

[TABLE]

which vanishes if we take into account the initial conditions, Eq. (78), for the center-of-mass position $\bm{\mathcal{R}}(0)$ of the BEC. Thus, by inserting Eqs. (141) and (B.1) into Eq. (138), the kinetic energy is finally given by

[TABLE]

In conclusion, by summing up the contributions given by Eqs. (136), (139) and (143), we arrive at Eq. (4.2.3) for the energy per particle of a BEC within the time-dependent TF approximation.

B.2 Energy terms for the Thomas-Fermi ground state

In this appendix we use the TF ground state $\phi_{\mathrm{TF}}(\bi{x})$, Eq. (29), to derive explicit expressions for the energy of a BEC given by Eq. (4.2.3). Since we have neglected the quantum pressure term in the derivation of the TF ground state, Eq. (29), we consequently neglect it also for the calculation of the energy within the TF approximation, that is

[TABLE]

Next we use Eq. (29), the transformation $\bm{\xi}=\bi{x}-\bm{\mathcal{R}}(0)$, and the relation (30) to rewrite the expectation value

[TABLE]

in the form

[TABLE]

which is valid in $d$ dimensions. Similarly, by applying Eq. (29) to Eq. (86) we obtain

[TABLE]

Finally, when we insert Eqs. (144) - (146) into Eq. (4.2.3), we arrive at the expression for the total energy of a BEC, Eq. (89), valid within the TF regime.

B.3 Angular momentum operator in original coordinates

Here we establish a relation for the angular momentum operator with respect to the center of mass in original coordinates. We start with the displaced wave function $\psi_{\mathrm{D}}(t,\bi{x})$, which is the solution of the GP equation (8) in the center-of-mass coordinates. In these coordinates the angular momentum operator $\hat{\bi}{L}_{\mathrm{D}}$ is defined in the usual way by

[TABLE]

With the help of the inverse of transformation (4) we can rewrite the expectation value of the angular momentum operator (147) as follows

[TABLE]

where \hat{D}\big{(}\bm{\mathcal{R}}(t),\bm{\mathcal{P}}(t)\big{)} is the displacement operator which fulfills the relations (119) and (120). Hence, the angular momentum operator in original coordinates reads

[TABLE]

B.4 Angular momentum in the time-dependent Thomas-Fermi approximation

We now evaluate the expectation value, Eq. (148), of the angular momentum operator defined by Eq. (149) with respect to the approximate wave function, Eq. (73), valid in the time-dependent TF approximation. By making use of the coordinate transformation (80) we obtain the relation

[TABLE]

Taking into account the initial condition, Eq. (78), for the center-of-mass position $\bm{\mathcal{R}}(0)$, the second term in Eq. (150) vanishes and we arrive at Eq. (92).

B.5 Angular momentum for the Thomas-Fermi ground state

When we consider the TF ground state $\phi_{\mathrm{TF}}(\bi{x})$, Eq. (29), as the initial state $\psi(0,\bi{x})$ for the expectation value of the angular momentum matrix operator, Eq. (4.2.4), the last two terms vanish, that is

[TABLE]

This results from the fact that the contributions of these two terms cancel each other when the momentum operator is applied to $|\phi_{\mathrm{TF}}(\bi{x})|$. In addition, the contributions due to the phase $\rmi\bm{\mathcal{P}}(0)\,\bi{x}/\hbar$ of the TF ground state, Eq. (29), vanish in consequence of the initial condition, Eq. (78), for $\bm{\mathcal{R}}(0)$. What remains is the contribution of the first term in Eq. (4.2.4), which can be evaluated with the help of Eq. (145). Hence, we obtain the expression

[TABLE]

for the expectation value of the angular momentum matrix operator within the TF regime. By applying the angular momentum matrix associated with the matrix differential equation as given by Eq. (70) to (152), we obtain Eq. (94).

Appendix C Solutions of the matrix differential equation for isotropic traps

In this appendix we find the solution of the non-linear differential equation (95) in the implicit form given by Eq. (98) valid for isotropic initial traps. In the case of such a trap, that is $\omega_{i}(0)=\omega_{0}$, all scaling parameters $\lambda_{i}(t)$ are identical, $\lambda_{i}(t)=\lambda(t)$, and Eq. (95) then reads

[TABLE]

with the initial conditions

[TABLE]

resulting from Eq. (96).

By considering Eq. (153) as the classical Newton equation for the “coordinate” $\lambda$, we can integrate this equation with the integral of motion

[TABLE]

which plays the role of “energy” corresponding to Eq. (153). Here we have used the initial conditions given by Eq. (154). After solving Eq. (155) with respect to $d\lambda/dt$ and using again Eq. (154), we obtain the implicit solution

[TABLE]

of Eq. (153) for the function $\lambda(t)$.