A variational approach to Bogoliubov excitations and dynamics of dipolar BECs

# A variational approach to Bogoliubov excitations and dynamics of dipolar Bose-Einstein condensates

Manuel Kreibich, Jörg Main and Günter Wunner 1. Institut für Theoretische Physik, Universität Stuttgart, 70550 Stuttgart, Germany
###### Abstract

We investigate the stability properties and the dynamics of Bose-Einstein condensates with axial symmetry, especially with dipolar long-range interaction, using both simulations on grids and variational calculations. We present an extended variational ansatz which is applicable for axial symmetry and show that this ansatz can reproduce the lowest eigenfrequencies of the Bogoliubov spectrum, and also the corresponding eigenfunctions. Our variational ansatz is capable of describing the roton instability of pancake-shaped dipolar condensates for arbitrary angular momenta. After investigating the linear regime we apply the ansatz to determine the dynamics and show how the angular collapse is correctly described within the variational framework.

###### pacs:
03.75.Kk, 67.85.De, 05.45.-a

## 1 Introduction

As an alternative to the direct numerical solution of the Gross-Pitaevskii equation (GPE) and the Bogoliubov-de Gennes equations (BDGE), which describe the ground states and excitations of Bose-Einstein condensates [1], Rau et al[2, 3] proposed a variational approach with coupled Gaussians to calculate stationary solutions and the stability of a condensate using the methods of nonlinear dynamics. Specifically, they used the time-dependent variational principle (TDVP) to map the GPE to a nonlinear dynamical system for the variational parameters whose fixed points correspond to ground states of the GPE. The stability can be analysed by means of the linearised equations of motion, i. e. in terms of the eigenvalues of the Jacobian. However, it remained unclear whether there is a direct relationship between the eigenvalues of the BDGE and those of the Jacobian. In a recent work [4] we analysed this problem for spherically symmetric systems and showed that a variational ansatz with coupled Gaussians has limitations in that it can only describe excitations with an angular momentum . We present an extended variational ansatz and show that there is indeed a good agreement between the lowest eigenvalues of the BDGE and the Jacobian for arbitrary angular momenta, with or without additional long-range interaction (LRI).

A variational ansatz with a single Gaussian [5] cannot reproduce the biconcave structure of dipolar condensates revealed in certain parameter regions [6, 7]. Using the variational ansatz of coupled Gaussians, however, it was possible to obtain non-Gaussians shapes [2, 3, 8] and to reproduce the wave function and the energy of structured ground states. Additionally, the stability was analysed which revealed that the collapse breaks the axial symmetry within the variational approach. However, the ansatz used in Ref. [8] is only capable of describing modes with angular momentum and (see discussion in [4]). The angular collapse with symmetry, calculated numerically on a grid in [9], is therefore not accessible to this variational ansatz of coupled Gaussians.

Several extensions of the Gaussian ansatz have been proposed in the literature [10, 11, 12], but they allow for no systematic inclusion of arbitrary angular momenta. It is the purpose of this work to extend the variational ansatz of coupled Gaussians in such a way that an instability with an arbitrary projection of the angular momentum on the axis can be described for systems with axial symmetry. We present an extended variational ansatz including angular exponentials , which are eigenfunctions of the component of the angular momentum operator. This enables us to describe modes with arbitrary angular momentum within the variational framework, which will be demonstrated by applying the ansatz to condensates with contact interaction only, and with dipole-dipole interaction (DDI).

In Sec. 2 we discuss the stability and excitations of BECs. We first present the method of numerically solving the BDGE, which we need for a comparison with the variational method. The extended variational ansatz is described and applied to calculate the eigenfrequencies of elementary excitations. We compare the results with those obtained from the full-numerical solution. In Sec. 3 we apply the variational ansatz to calculate the time evolution of a dipolar BEC and demonstrate that the ansatz can describe the angular collapse. In Sec. 4 we draw conclusions and give an outlook on future work.

## 2 Stability and excitations

The dime-dependence of a BEC is described in the mean-field approximation by the time-dependent GPE [1]

 iℏ∂ψ∂t(r,t)=[−ℏ22MΔ+M2ω2ρ(ρ2+λ2z2)+NΦint(r,t)]ψ(r,t), (1)

where is the particle number,

 Φint(r,t)=∫d3r′Vint(r−r′)∣∣ψ(r′,t)∣∣2 (2)

the mean field produced by the inter-particle interaction, the wave function of the condensate normalised to unity, and the mass of the atom species. The harmonic external trap is given by the trapping frequencies and in the and direction, and can be characterised by the trap aspect ratio . We measure all lengths in units of the harmonic oscillator length , all frequencies in units of , and the time in units of . A stationary solution of Eq. (1) has the time dependence , where is the chemical potential of the condensate.

The interaction potential is given by [13]

 Vint(r−r′)=4πaℏ2Mδ(r−r′)+μ0μ2m4π1−3(z−z′)2(r−r′)2|r−r′|3. (3)

The quantity is the s-wave scattering length and the magnetic moment of the atom species and is, e. g., ( the Bohr magneton) for \ce^52Cr [14]. The dipoles are assumed to be aligned along the axis by an external magnetic field. The interactions can be characterised by the dimensionless parameters and [15]. In this article, we investigate BECs with and without LRI. For the condensate with DDI we use values of and , since for these parameters structured ground states appear [6, 3]. For the BEC without LRI we use a trap aspect ratio of , which was used in experiments [16]. The scattering length can be tuned in a wide range via Feshbach resonances [17]. Thus, we use this quantity as a free parameter in our calculations.

### 2.1 Full-numerical treatment

Besides the variational approach we also perform grid calculations to directly solve the GPE and BDGE, which allows us to compare both methods. Before solving the BDGE, one first needs to find the ground state of a BEC. In this work, we use the imaginary time evolution (ITE) for this task. We evolve an initial wave function on a grid in imaginary time defined by , and as time evolves the wave function converges to the ground state. The time evolution on the grid is done via the split operator method [18], in which the time evolution operator in imaginary time is given by

 e−^HΔτ/ℏ=e−12^TΔτ/ℏe−^VΔτ/ℏe−12^TΔτ/ℏ+O(Δτ3), (4)

where and represent the kinetic and the potential terms of the mean-field Hamiltonian, respectively. The action of the kinetic part is trivial in Fourier space, whereas the action of the potential part is trivial in real space. Since we are dealing with a cylindrically symmetric system, we can use the Fourier-Hankel algorithm of [15] to compute the necessary Fourier transforms for a time step and the dipolar integrals via the convolution theorem.

To derive the BDGE, one starts from the usual ansatz for the perturbation of a stationary state with chemical potential [1]

 ψ(r,t)=[ψ0(r)+λ(u(r)e−iωt+v∗(r)eiω∗t)]e−iμt/ℏ, (5)

where is the frequency and the amplitude of the perturbation (). After inserting this ansatz into the time-dependent GPE (1), while neglecting terms of second order in , and collecting terms evolving in time with and , respectively, one obtains the BDGE

 ℏωu(r)= [H0−μ+∫d3r′Vint(r−r′)∣∣ψ0(r′)∣∣2]u(r) +∫d3r′Vint(r−r′)ψ∗0(r′)u(r′)ψ0(r) +∫d3r′Vint(r−r′)ψ0(r′)v(r′)ψ0(r), (6a) −ℏωv(r)= [H0−μ+∫d3r′Vint(r−r′)∣∣ψ0(r′)∣∣2]v(r) +∫d3r′Vint(r−r′)ψ∗0(r′)u(r′)ψ∗0(r) +∫d3r′Vint(r−r′)ψ0(r′)v(r′)ψ∗0(r), (6b)

where . Due to the cylindrical symmetry, we can make an ansatz for the solutions

 u(r)=eimϕu(ρ,z), v(r)=eimϕv(ρ,z), (7)

where is the usual quantum number of the projection of the angular momentum operator on the axis. Since the Laplacian in cylindrical coordinates applied to the ansatz (7) yields

 Δeimϕu(ρ,z)=(Δρ,z−m2ρ2)eimϕu(ρ,z), (8)

which only depends on , there is a degeneracy of the eigenmodes for . Thus, we can restrict ourselves to solutions with . Another symmetry of the system is the reflection at the plane (). Therefore, the solutions can be described by the parity under this reflection. We denote the set of quantum numbers with . There are always solutions of the BDGE for and the trapping frequencies, which represent the gauge mode and the centre of mass oscillations [1, 4].

We solve the linear, nonlocal BDGE (6) by calculating the matrix elements of the right-hand side in coordinate representation by means of the Fourier-Hankel algorithm [15]. The lowest eigenvalues of the resulting high-dimensional matrix are then calculated using the Arnoldi iteration [19], implemented in the Fortran ARPACK routines [20].

The solutions of the BDGE (6) yield the frequencies and the shape of elementary oscillations of the condensates. But it may happen that the Bogoliubov spectrum of the ground state of a dipolar BEC contains imaginary frequencies [9], which correspond to dynamical instabilities. A small perturbation of the ground state then leads to an exponential decay.

We calculated the ground states and Bogoliubov spectra for a BEC with an attractive short-range interaction (SRI) only () and a trap aspect ration of , and for a dipolar BEC with and , for which it was shown that the structured ground state changes its stability in a pitchfork bifurcation [8] at the scattering length 111In [8] a Gaussian ansatz with Gaussians was used. The scattering length, where the bifurcation takes place, however, converges slowly with increasing number of Gaussians and may be different in our calculations.. The structured ground states, where the maximum of particle density lies away from the centre, appear in isolated regions in parameter space with [6]. For our choice of parameters, and , the ground state assumes a biconcave structure for the scattering lengths .

For both systems ground states only exist for scattering lengths larger than the critical scattering length, . For the condensate with SRI, the critical scattering length is given by . BECs at small negative scattering lengths may exist, since the kinetic energy compensates the attractive interaction and stabilises the condensate. The critical scattering length was probed experimentally in a condensate of \ce^85Rb, and significant deviations from the results of the GPE were measured [21], which may be explained by taking into account higher-order nonlinear effects [22]. The critical scattering length of a dipolar BEC considerably depends on the trap geometry, in our case (, ) its value is given by . A good agreement between experiment and the results of the GPE was found [23].

Fig. 1 shows the Bogoliubov spectra. Since both systems have already been considered in the literature (see, e. g., [24] for the BEC with SRI, and [15] for the dipolar BEC), we only give a short review. The spectrum of the BEC without LRI (Fig. 1a) can be understood by taking the limit , which yields the exact eigenvalues of the cylindrically symmetric harmonic oscillator

 ωnρ,nz,mωρ=2nρ+|m|+λnz, (9)

where are the quantum numbers for excitations along the and direction, respectively. When decreasing the scattering length towards the critical value, below which no stationary solutions exist anymore, the frequencies vary slightly, which means that in this regime the dynamics is dominated by the external trap, and the interaction acts as a quasi perturbation. Just slightly above the critical point the lowest mode with drops to zero, which marks the global collapse of the condensate. The oscillations of the centre of mass, which constantly lie at the trapping frequencies, are represented by the lowest odd mode () and the lowest even mode ().

The Bogoliubov spectrum of the dipolar condensate (Fig. 1b) indicates a stable condensate when the scattering length is far from the critical point. However, for decreasing scattering length, not only the mode with but also frequencies for higher angular momenta decrease and drop to zero. At these points, those modes become unstable, as the eigenfrequencies are imaginary. Fig. 2 shows the real and imaginary parts of the lowest eigenfrequencies for the angular momenta . The modes with and are the first ones which turn unstable, and thus mark the local collapse of the condensate. This means, that collapse occurs via local density-fluctuations in contrast to a global collapse, where the condensate collapse to the trap centre [9]. In Sec. 3 we discuss the dynamics of such a local collapse, with the result that an extended variational ansatz is necessary to describe the collapse. This is a unique feature of dipolar BECs and is related to the biconcave structure of the ground state occurring for specific trapping geometries [8, 9]. These instabilities are referred to as rotons, or, more specifically, as discrete “angular rotons”, since these excitations assume a shape with azimuthal nodal surfaces [6].

Worth mentioning also is the avoided crossing in the lowest modes for the angular momenta , whereas for there is a regular crossing since is an exact solution for every scattering length, which does not allow for an avoided crossing.

### 2.2 Variational approach

The grid calculations are very accurate, if grid size and resolution are chosen carefully, but computationally expensive. The variational method provides an alternative in that the computation time reduces significantly compared to grid calculations. With this motivation, Rau et al[2, 3, 8] used the variational ansatz

 ψ=NG∑k=1ei(Akxx2+Akyy2+Akzz2+γk) (10)

to calculate stationary solutions and excitations of dipolar condensates. Here, is the number of Gaussians, , and describe the widths and the amplitude of each Gaussian. The quantities are the variational parameters of this ansatz. In a recent article [4] we showed that the variational ansatz with coupled Gaussians can only describe modes with an angular momentum of , . For that reason we proposed an extended variational ansatz with coupled Gaussians and spherical harmonics, which can describe modes with arbitrary angular momenta in spherically symmetric systems. We found a good agreement of the lowest eigenmodes with those obtained by grid calculations.

For the systems with axial symmetry considered here, this ansatz is not applicable, since the angular momentum is no longer a good quantum number and therefore the angular part of the excitations is not described by spherical harmonics. However, the component of the angular momentum, represented by the quantum number and the eigenfunctions , is still conserved. Using this knowledge we construct an extended variational ansatz, which is specific to cylindrically symmetric systems, and which is given by

 ψ=NG∑k=1⎛⎝1+∑m≠0∑p=0,1dkm,pρ|m|zpeimϕ⎞⎠e−Akρρ2−Akzz2−pkzz−γk. (11)

The complex variational parameters and describe the width in the and direction, is the displacement along the axis, and describes the amplitude of each Gaussian. With these parameters alone, the wave function does not depend on the angular coordinate , thus the parameters are responsible only for any excitation. To account for higher angular momenta and arbitrary parity , we include the sums over and in the brackets ( belongs to , and to ). The amplitude of each of these excitations is given by the variational parameters .

To derive the equations of motion for the variational parameters, we follow the procedure in [2, 4] and apply the Dirac-Frankel-McLachlan TDVP [25, 26]. It states that the norm of the difference between the left- and the right-hand side of the Schrödinger or GPE

 I=||iχ(t)−^Hψ(t)||2 (12)

must be minimised for a given variational ansatz with respect to the variational parameters . The variation of , while is fixed, and the identification of with leads to the necessary condition [27]

 K˙z=−ih, (13)

where the matrix and the vector are defined by

 Kij=⟨∂ψ∂zi∣∣∣∂ψ∂zj⟩, hi=⟨∂ψ∂zi∣∣∣^Hψ⟩. (14)

The remaining task is to evaluate all integrals appearing in Eq. (14). Technical details are presented in A. Since the ground state in the GPE corresponds to a fixed point of the equations of motion (13), we can find the ground state of a system by a nonlinear root search, for which we require

 ˙zi={iμfor zi≡γk,0else, (15)

where is the chemical potential.

Stability and excitations in the framework of the TDVP are given by the linearised equations of motion, for which one first needs to rewrite the complex equations of motion (13) into real ones by defining the new vector of variational parameters as . The time-dependence of a small perturbation is then given by [2]

 δ˙~zi=∑j∂˙~zi∂~zj∣∣∣~z=~z0δ~zj≡∑jJijδ~zj, (16)

with the Jacobian evaluated at the fixed point . These linearised equations of motion are as usual solved by a harmonic time-dependence of the perturbation

 δ~z(t)=eiωtδ~z(0), (17)

which leads to the eigenvalue problem

 iωδ~z(0)=Jδ~z(0). (18)

A real belongs to an oscillatory motion, as for real eigenvalues of the BDGE. The symmetry of the system, which leads to the quantum numbers and in the Bogoliubov spectrum, gives rise to a block structure of the Jacobian which is further explained in B.

The eigenvalues of the Jacobian and the BDGE can be directly compared. But is there also a relationship between the eigenvectors of the Jacobian and the eigenfunctions of the BDGE and ? For a given eigenvector , a wave function can be constructed by taking the difference of the excitation with variational parameters and the wave function of the stationary solution with the parameters ,

 δψ(~z0,δ~z)≡ψ(~z0+δ~z)−ψ(~z0). (19)

Since we are dealing with an excitation in the linear regime, one can Taylor expand the first summand in Eq. (19) to first order in . This yields the function [2]

 δψ=δ~z⋅∂ψ∂~z∣∣∣~z=~z0, (20)

which is the scalar product of the eigenvector and the gradient of the wave function with respect to the variational parameters evaluated at the stationary solution. Thus, every eigenvector belongs uniquely to a function .

As the Jacobian is real and a stable mode is represented by an imaginary eigenvalue , for every eigenvalue with there is also an eigenvalue , therefore we can construct corresponding functions and which evolve in time as and , respectively. They can be combined to the perturbation

 δψ=δψ−e−iωt+δψ+eiωt. (21)

Comparing this with the Bogoliubov ansatz from Eq. (5) leads to the relationship

 u↔δψ−, v∗↔δψ+, (22)

which links the eigenvectors of the Jacobian with the eigenfunctions of the BDGE.

### 2.3 Results for the Bec with scattering interaction

We now apply the variational ansatz (11) to a condensate with an attractive SRI. Fig. 3 shows the results, where coupled Gaussians have been used (solid lines), compared to the full-numerical calculations (dotted lines). For the sake of clarity we distinguish between even modes in Fig. 3a and odd modes in Fig. 3b. One recognises that the lowest mode of each angular momentum and parity can perfectly be reproduced by the variational ansatz. Interestingly, as with our previous ansatz for spherically symmetric systems [4], the frequencies of the centre of mass oscillations (lowest odd at , lowest even at ) agree, within numerical accuracy, with an arbitrary number of coupled Gaussians, even with . For the second lowest modes we also find a very good agreement as long as the scattering length is not close to the critical value. For the case , even more eigenvalues of the Jacobian are close to the numerical results.

To demonstrate the power of our extended ansatz, we calculated eigenmodes for angular momenta up to (see Fig. 4) for a fixed scattering length . Even for these angular momenta the lowest modes agree well. For the second excitations and lower angular momenta we find also a good agreement, however, the deviations of the variational and numerical results become larger for increasing angular momenta. Nevertheless, there is an obvious relationship between the eigenfrequencies of the Bogoliubov spectrum and the eigenvalues of the Jacobian, where our extended variational ansatz (11) – specific to cylindrically symmetric systems – has been used. Thus, the variational ansatz is a valid alternative to simulations on grids if one is interested only in the lowest modes.

### 2.4 Results for the dipolar Bec

We now want to include the DDI in the calculations and especially verify if the extended variational ansatz can describe the stability change for higher angular momenta . We first investigate the eigenmodes in the stable regime with a scattering length . Fig. 5 shows the comparison of both spectra for the angular momenta . With coupled Gaussians, the lowest modes of each angular momentum and parity (solid lines) are converged to their corresponding numerical frequencies (dotted lines). As in the case of the BEC without LRI, the centre of mass oscillations agree within numerical accuracy. For the modes with even parity, there is a good agreement even for the second lowest modes, where the differences become larger for rising scattering length and angular momentum. In the case of odd parity, the deviations for the second lowest modes become quite large for , when the scattering length is large.

We have also investigated eigenmodes with an angular momentum up to in the stable regime. Fig. 6 shows the results for a fixed scattering length . The variational ansatz can reproduce the lowest modes even for these angular momenta and, since the scattering length is not too large, also the second lowest modes. For the case , , we find the lowest three modes converged, and a good agreement of the fourth lowest mode. These investigations show that the variational ansatz (11) can also reproduce the full-numerical Bogoliubov spectrum, if only the lowest modes are compared.

We now investigate the interesting case of a small negative scattering length, where the Bogoliubov spectrum shows unstable roton modes for different angular momenta (see Fig. 2). We focus here on the case , which cannot be described by the ansatz of coupled Gaussians Eq. (10), hence the extended variational ansatz is necessary. Fig. 7 shows the lowest modes with obtained using different numbers of Gaussians (solid and dashed lines). For comparison the full-numerical results are also plotted (dotted line). We find that with an increasing number of coupled Gaussians the accuracy of the variational calculations increases gradually.

Other angular momenta show a similar behaviour, the convergence is in general better for lower angular momenta. Thus we can conclude, that the variational ansatz (11) can describe, though not quantitatively perfect, the dynamical instabilities of the Bogoliubov spectrum for angular momenta .

It is known that at the critical scattering length below which no stationary solution exists anymore an unstable collectively excited state emerges in a tangent bifurcation [3] in addition to the ground state (see the inset in Fig. 8). For condensates without DDI, this excited state is unstable with respect to an excitation [3]. The variational approach offers the possibility to investigate this state of a dipolar BEC, which is not accessible to the full-numerical ITE method, and its collapse mechanism. Fig. 8 shows the imaginary frequencies of the excited state obtained with Gaussians. As both states merge at the critical scattering length the same holds for the eigenfrequencies, and the excited state shares the instabilities of the ground state at the critical scattering length, except the additional instability. In contrast to the case with SRI only or with monopolar LRI [4] the excited state of the dipolar BEC for the trapping frequencies is unlikely to collapse with symmetry due to the additional instabilities for , which have a higher imaginary frequency and thus a larger collapse rate.

The excitation in Fig. 8 changes its stability with increasing scattering length, in a similar way as the unstable modes of the ground state (see Fig. 7). Note, however, that this does not mean a global change of the stability due to the other unstable modes. The excited state stays unstable for all scattering lengths considered in our calculations.

### 2.5 Shape of Bogoliubov excitations

After linking the eigenvalues of the BDGE and the Jacobian of the time-dependent variational approach, we will check if there is also a direct correspondence between the eigenfunctions and eigenvectors of these two approaches. In Sec. 2.2 we discussed how to construct the functions and from the eigenvectors and how they are related to the functions and . For a better comparison of both methods, we use the usual normalisation of the Bogoliubov functions [1], and assume the same for the variational method, . We point out that in calculating there is no ambiguity, except for a global phase, which we choose in such a way that the functions are real. Fig. 9 shows the comparison of both methods for some modes. The first two functions with have one and two radial nodes, respectively, because the lowest mode with represents the gauge mode, which is no physical excitation. For other angular momenta, the -th excitation in direction has radial nodes. The asymptotic behaviour of all functions in Fig. 9 for is . The nodal structure for the functions is different. All for have one node, and for the other angular momenta zero nodes. For higher excitation numbers and angular momenta, the amplitude of becomes smaller and smaller, indicating that in these limits the coupling between and in the BDGE can be neglected, which is in accordance with previous calculations in the literature [28, 1, 15].

For the second excitation with we find small deviations, which is to be expected, since there are also deviations in the eigenfrequencies (cf. Fig. 6). By comparing more eigenmodes we can conclude: The better the agreement of the eigenfunctions, the better is also the accordance of the frequencies. For higher modes, where there is no quantitative agreement between the eigenfrequencies, only the asymptotic behaviour and the nodal structure can be described by the variational ansatz, even though there is at least a qualitative agreement in those aspects.

The comparison of the eigenfunctions of the BDGE with those obtained from the eigenvectors of the Jacobian reveals a very good agreement. This is remarkable, since both functions originate from quite different approaches of stability. The extended variational ansatz (11) cannot only quantitatively describe the frequencies of the lowest modes, but also the shape of the oscillations, even for angular momenta, which were not accessible by the variational approach without the extension.

## 3 Dynamics

The extended variational ansatz (11) was proposed to find a direct relationship between the eigenvalues of the Jacobian and the eigenfrequencies of the BDGE. We have derived the equations of motion using the TDVP to calculate the Jacobian. The equations of motion (13) are not only valid in the vicinity of a fixed point, but for all variational parameters for which the integrals in the equations of motion converge. Wilson et al[9] simulated the dynamical evolution of a collapse for a dipolar condensate with biconcave structure by reducing the scattering length below the critical value, which yields a symmetry breaking collapse (angular collapse) with symmetry. Rau et al[3] simulated this scenario within the variational framework with an ansatz of coupled Gaussians. However, as mentioned in Sec. 2.2, this ansatz cannot describe wave functions that exhibit the symmetry. Motivated by these papers we will now integrate the equations of motion numerically with appropriate initial conditions to verify, whether or not the extended variational ansatz is also applicable in calculating the dynamics of dipolar condensates and especially to describe the angular collapse for arbitrary angular momenta, which is caused by the roton instability (cf. Secs. 2.1 and 2.4). To do so we only address the collapse of ground states with biconcave structure.

As a first example, we consider the system with Gaussians and calculate the dynamics of an induced collapse by lowering the scattering length below the critical value. As an initial state for , we take the stable ground state for a scattering length of and add a small perturbation for a specific angular momentum by changing the variational parameter to a small value larger than zero, e.g. , which is just a small perturbation of the ground state. We let the condensate evolve until , after which the scattering length is linearly ramped down until to a value of , which is below the critical scattering length. The condensate then evolves at a fixed scattering length, until it collapses to a point, from where on no further integration is possible. We calculated this scenario for an initial and excitation, which cannot be described by the variational ansatz without the extension and thus is an advance as compared to previous calculations of a collapsing condensate within the variational framework [3].

Fig. 10 shows the resulting density of the collapsing condensate for both cases. For an initial excitation, the collapsing cloud exhibits the clear shape with an symmetry, the angular dependency is proportional to . This result shows that the extended variational ansatz is not only able to describe the instability in the vicinity of the ground state, but also the nonlinear dynamics given by the equations of motion (13). If the initial state is perturbed by an excitation, the resulting density distribution shows a superposition of several angular momenta. There are altogether four density peaks, of which two are significantly larger than the other two. The appearance of four peaks corresponds to an symmetry, but the distribution of the particle density clearly belongs to . Of the two large peaks, one is slightly higher than the other one which indicates an symmetry. The calculation shows that the different angular momenta do not evolve freely, but instead are coupled and influence each other. This is due to the nonlinearity of the GPE (1) which carries over to the equations of motion for the variational parameters.

There are many excitations one can add to the ground state as a perturbation. Above, we added specific eigenmodes to investigate the behaviour of the extended variational ansatz when considering the collapse dynamics. We now want to utilise a more physical motivated initial state. In Ref. [9] the authors added modes with weights determined by the Bose-Einstein distribution, which leads to

 ψ0(ρ,z)+∑n,m√nn,mNe2πiαn,m[un,m(ρ,z)eimϕ+v∗n,m(ρ,z)e−imϕ], (23)

where is the ground state. The weight of each excitation, where is the angular momentum and summarises the other quantum numbers, is given by

 nn,m=(eℏωn,m/kBT−1)−1. (24)

Here, is the energy of each excitation, the Boltzmann constant, and the temperature. The quantities in Eq. (23) determine the initial phase of each excitation and are random numbers between and . The overall factor , with the particle number, is necessary, since the wave function is normalised to unity instead of .

We now can make use of the relationship between the eigenvectors of the Jacobian and the eigenfunctions of the BDGE (Sec. 2.5) to prepare an initial state similar to that given by Eq. (23) within the variational framework, where instead of the functions, the eigenvectors of the variational parameters have to be added to the ground state vector. For the following calculations we simulate a condensate consisting of \ce^52Cr atoms at a temperature of . We start with an initial wave function according to Eq. (23) and a scattering length of , which is ramped down to a value of over a time of .

We performed several calculations with this scenario, for which the only differences are the initial phases . As the scattering length decreases the particle density establishes its biconcave structure, which is a typical feature of dipolar BECs [9, 8, 6, 7]. If the scattering length further decreases, the critical value is reached and the atomic cloud starts to collapse to the region with the highest density, which is on a ring around the centre. Due to the initial excitations with angular momentum , which are now unstable, the symmetry is broken and several density peaks emerge. Fig. 11 shows the density distribution of four different calculations shortly before collapse. Each condensate has three peaks located almost equidistantly on a ring. But the distribution along the ring differs: In Figs. 11 (a)-(b) there are two large peaks and a smaller one, while in Figs. 11 (c)-(d) it is the other way round. Additionally the initial phases determine the overall rotation of the collapsed condensate in the -plane. We can conclude that both, the and mode, are responsible for the shape of the collapsed cloud, which corresponds to the numerical solution of the BDGE (see Fig. 2), in which both modes turn unstable almost simultaneously.

A question, which naturally arises, is how the number of Gaussians influences the dynamics calculation. We performed similar calculations with and Gaussians. With Gaussians we find a different behaviour in that a collapse with symmetry is favoured. However, with Gaussians almost the same results as in the case with Gaussians are obtained. The only difference is the time taken before the condensate collapses, which is shorter with Gaussians. This can be explained by the fact that the critical scattering length is shifted to a higher value and thus the critical point is reached earlier during the simulation. The eigenfrequencies with our Gaussian ansatz do not fully agree with the full-numerical results (cf. Sec. 2.4), however, we expect the variational results to not being completely converged with or Gaussians.

## 4 Conclusion and outlook

We presented an extended variational ansatz, which is specific to cylindrically symmetric systems and which can describe arbitrary angular momenta, and derived the equations of motion for the variational parameters using the TDVP. The eigenfrequencies of two different systems (without and with DDI) were calculated using the variational ansatz and compared with the direct solutions of the BDGE. We found a good agreement of the lowest-lying modes and could confirm the direct relationship of both methods. This extended our investigations of spherically symmetric systems [4] and we could show that this agreement is of generic type, independent of the symmetry and the interactions.

The roton-instability of a pancake-shaped dipolar BEC was investigated using the extended ansatz. We could find instabilities for arbitrary values of the angular momenta, which is a great progress in comparison with the variational ansatz with coupled Gaussians without the extension. With more and more Gaussians we obtained gradually more accurate results within the variational framework. Thus, our variational method is a valid approximation to grid calculations.

We further constructed functions of the eigenvectors of the Jacobian and found that these functions agree very well with the solutions of the BDGE, as long as the corresponding eigenvalues agree. There is not only a direct link between the eigenfrequencies, but also between the eigenvectors and the eigenfunctions.

After considering the dynamics in the linear regime we directly integrated the nonlinear equations of motion and showed that our variational ansatz, assuming that the initial conditions and the scenario of the time-dependent scattering length are appropriately chosen, is capable of describing the angular collapse with arbitrary angular momenta. We used the relationship between the eigenvectors of the Jacobian and the eigenfunctions of the BDGE to construct an initial state which corresponds to a condensate in thermal equilibrium within the mean field approximation, and simulated a realistic experiment of the angular collapse. We found that the random initial phase distribution strongly influences the shape of the collapsed atomic cloud, which makes the resulting shape of an experiment irreproducible.

Altogether it was shown that the extended variational ansatz is a significant progress compared to coupled Gaussians, and as yet unaccessible results could be obtained. Our ansatz is a valid alternative to calculate eigenfrequencies and shapes of elementary excitations, if the lowest eigenmodes of stable condensates are considered. The roton-instabilities can be obtained, but we could not find a good quantitative agreement, leaving our ansatz an approximation.

In future applications of the extended variational ansatz it would be possible to include the particle loss in the dynamics calculation, as the losses are necessary to accurately simulate a collapse scenario of a dipolar BEC [9, 29]. Furthermore, our ansatz is applicable to dipolar BECs in a one-dimensional optical lattice [30, 31, 32] to calculate the ground state, stability, excitations and the dynamics. The direct relationship between the eigenvectors of the Jacobian and the eigenfunctions of the BDGE allows several applications of the variational approach, e. g. to describe a BEC at finite temperature, where a self-consistent solution of the GPE and BDGE is necessary [33].

This work was supported by Deutsche Forschungsgemeinschaft.

## Appendix A Calculation of the integrals for the variational ansatz

This appendix is dedicated to the integrals necessary for the matrix and the vector , which are defined by

 Kij =⟨∂ψ∂zi∣∣∣∂ψ∂zj⟩=∫d3r(∂ψ∂zi(r))∗∂ψ∂zj(r), (25a) hi =⟨∂ψ∂zi∣∣∣^Hψ⟩=∫d3r(∂ψ∂zi(r))∗^Hψ(r), (25b)

to set up the equations of motion (13) of the variational ansatz (11) via the time-dependent variational principle. We write the variational ansatz in the form

 ψ=NG∑k=1∑m∑p=0,1dkm,pρ|m|zpeimϕe−Akρρ2−Akzz2−pkzz−γk, (26)

where and are no variational parameters, but constants. We further introduce the abbreviations

 ⟨zi∣∣zj⟩≡⟨∂ψ∂zi∣∣∣∂ψ∂zj⟩, ⟨zi|Xψ⟩≡⟨∂ψ∂zi∣∣∣Xψ⟩ (27)

for the matrix elements, where is an operator of the mean field Hamiltonian with the kinetic part, the external potential and the scattering and dipolar mean field potentials.

We do not give a detailed way on how to calculate all necessary integrals, nor do we present all results. We rather want to sketch the calculations and give results of some examples. For all details we refer to the master thesis of one of the authors [34].

### a.1 Integrals of the K-matrix

Before giving the results, we define the abbreviations

 Aklρ≡Akρ+(Alρ)∗, Aklz≡Akz+(Alz)∗, pklz≡pkz+(plz)∗, γkl≡γk+(γl)∗ (28)

and the basic integrals

 Ikl(m,n) ≡∞∫0dρρ2m−1e−Aklρρ2∞∫−∞dzzne−Aklzz2−pklzze−γkl=e−γklΓ(m)2(Aklρ)mIklz(n). (29)

The integral can be identified with the gamma function by a simple substitution. The integral is a Gaussian integral for , higher orders can be obtained by the recursion formula

 Iklz(n+1)=−∂∂pklzIklz(n). (30)

The integrals of the matrix then yield

 ⟨dlm2,p2∣∣dkm1,p1⟩ =2πδm1m2Ikl(|m1|+1,p1+p2), (31a) ⟨dlm2,p2∣∣Akρ⟩ =−2π∑p1Ikl(|m2|+2,p1+p2), (31b) ⟨dlm2,p2∣∣Akz⟩ =−2π∑p1Ikl(|m2|+1,p1+p2+2), (31c) ⟨dlm2,p2∣∣pkz⟩ =−2π∑p1Ikl(|m2|+1,p1+p2+1), (31d) ⟨dlm2,p2∣∣γk⟩ =−2π∑p1Ikl(|m2|+1,p1+p2). (31e)

The other integrals can similarly be expressed by Eq. (29).

### a.2 Integrals of the kinetic term

For the integrals of the kinetic term the action of the Laplacian on the variational ansatz is necessary. This is most easily evaluated in cylindrical coordinates. The integrals can then be expressed by Eq. (29). For the parameters, e. g., one obtains

 (32)

### a.3 Integrals of the trapping potential

The trapping potential is of the form . Thus, the corresponding integrals are easy to evaluate. As an example we obtain

 (33)

### a.4 Integrals of the scattering interaction

The integrals for the interaction terms require some more attention. First an expression for the modulus squared wave function is needed. This enters the integrals for the interaction potentials. The scattering interaction is easier to calculate, since this kind of interaction is local, in contrast to the DDI.

We further need new abbreviations, which depend on four indices and which are defined by

 Aijklρ≡Aijρ+Aklρ, Aijklz≡Aijz+Aklz, pijklz≡pijz+pklz, γijkl≡γij+γkl. (34)

The new basic integral, , can be obtained from the former one, Eq. (29), by substituting all indices with . The scattering integral then reads

 (35)

### a.5 Integrals of the dipolar interaction

In order to calculate the integrals of the dipolar interaction, we first need the expression for the mean field potential . This potential can be rewritten using Fourier transforms and the convolution theorem, which yields [35]

 Φd(r)=F−1[~Vd(k)⋅~n(k)], (36)

where the functions with tilde designate the Fourier transforms of the dipole potential and the particle density. The transform of the dipole potential is well-known to be [35]

 ~Vd(k)/ℏωρ=4π3(3k2zk2−1)D/N. (37)

Our next task is to calculate . Writing the square of the absolute value of the wave function with the given variational ansatz, one notices that the angular dependency is of the form with some . In [15] it was noted that the three-dimensional Fourier transform of such a function can be expressed as a one-dimensional Hankel transform for the direction and a remaining one-dimensional Fourier transform in the direction. The Hankel transform brings in the Bessel function, and the integral can be worked out by Taylor expanding the Bessel function. The Fourier transform is a Gaussians integral with linear terms and variable powers of , for which an analytical expression can be found in established tables of integrals [36, Eq. 3.462.2].

In the original integral, , the integral over is similar to the Fourier transform of and can easily be expressed in the same way, which yields the intermediate result

 (38)

where we defined

 σm≡(signm)m, μm1,m2≡{0for signm1≠signm2,min(|m1|,|m2|)otherwise, (39)

and used the Pochhammer symbol , which is defined by

 (a)n≡a(a−1)(a−2)⋯(a−n+1)=Γ(a+1)Γ(a−n+1). (40)

The remaining three dimensional integral reads

 (41)

with the integrals

 ~Iijz(p) ≡∞∫−∞dzzpe−Aijzz2−(pijz+ikz)z, (42a) ~Ikl−z(p) ≡∞∫−∞dzzpe−Aklzz2−(pklz−ikz)z, (42b)

which are evaluated in Ref. [36, Eq. 3.462.2]. The Fourier transform of the dipole potential consists of two parts, . The integral of containing is equivalent to the integrals of the scattering interaction222The Fourier transform of is , and thus proportional to . This means that the whole integral of this part is equivalent to the integrals of the scattering interaction., and requires no further attention. The part of with , which we denote by , can be rewritten to

 I1k/ℏωρ =4π2(D/N)δm1+m3,m2+m4n!(αijklρ)n ×∞∫1dtt−(n+1)∞∫−∞dkzk2ze−(t−1)αijklρk2z~Ikl−z(p1+p2)~Iijz(p3+p4). (43)

The integral is of Gaussian type, and similar to and . The remaining integral can be expressed by the hypergeometric Function [37], assuming all parameters are zero. In the general case, the integral can be efficiently evaluated using a Taylor series.

## Appendix B Structure of the Jacobian

The solutions of the BDGE (6) can be classified by the quantum numbers and , which express the symmetries of the system. Applying the TDVP and linearizing the equations of motion yields an alternative description of linear oscillations and stability, in which the Jacobian plays a central role. The Jacobian, similar to the BDGE, also express the symmetries of the system and the quantum numbers, in that it assumes a block structure, which can in general be written as

 (44)

where is the submatrix belonging to the angular momenta and the parity . This block structure allows one to set up and diagonalise the Jacobian for each and individually, analogously to solving the BDGE. The only difference is that in the Jacobian the angular momenta and are coupled, which is due to the coupling of those angular momenta in the integrals necessary for the equations of motion (cf. A).

A simple example, which shows that the coupling may not be neglected, is the submatrix of the , Jacobian obtained with Gaussian, for a dipolar BEC with and

 J+1±1=⎛⎜ ⎜ ⎜⎝0−2.59360−2.39312.59360−2.393100−2.39310−2.5936−2.393102.59360⎞⎟ ⎟ ⎟⎠ωρ, (45)

where the variational parameters are ordered as for a scattering length of . The eigenvalues of the matrix (45) are . The eigenvectors of the eigenvalue can then be combined to an eigenvector belonging to and another eigenvector for , the same holds for the eigenvalue . In general, for arbitrary number of Gaussians and angular momenta, by taking appropriate linear combinations of the eigenvectors, eigenvectors belonging to a specific angular momentum or can be obtained.

## References

• [1] L. P. Pitaevskii and S. Stringari. Bose-Einstein Condensation. Oxford University Press, 2003.
• [2] S. Rau, J. Main, and G. Wunner. Variational methods with coupled Gaussian functions for Bose-Einstein condensates with long-range interactions. I. General concept. Phys. Rev. A, 82:023610, 2010.
• [3] S. Rau, J. Main, H. Cartarius, P. Köberle, and G. Wunner. Variational methods with coupled Gaussian functions for Bose-Einstein condensates with long-range interactions. II. Applications. Phys. Rev. A, 82:023611, 2010.
• [4] M. Kreibich, J. Main, and G. Wunner. Relation between the eigenfrequencies of Bogoliubov excitations of Bose-Einstein condensates and the eigenvalues of the Jacobian in a time-dependent variational approach. Phys. Rev. A, 86:013608, 2012.
• [5] V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller. Dynamics of Bose-Einstein condensates: Variational solutions of the Gross-Pitaevskii equations. Phys. Rev. A, 56:1424–1432, 1997.
• [6] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn. Radial and Angular Rotons in Trapped Dipolar Gases. Phys. Rev. Lett., 98:030406, 2007.
• [7] O. Dutta and P. Meystre. Ground-state structure and stability of dipolar condensates in anisotropic traps. Phys. Rev. A, 75:053604, 2007.
• [8] S. Rau, J. Main, P. Köberle, and G. Wunner. Pitchfork bifurcations in blood-cell-shaped dipolar Bose-Einstein condensates. Phys. Rev. A, 81:031605(R), 2010.
• [9] R. M. Wilson, S. Ronen, and J. L. Bohn. Angular collapse of dipolar Bose-Einstein condensates. Phys. Rev. A, 80:023614, 2009.
• [10] D. Buccoliero, A. S. Desyatnikov, W. Krolikowski, and Y. S. Kivshar. Laguerre and Hermite Soliton Clusters in Nonlocal Nonlinear Media. Phys. Rev. Lett., 98:053901, 2007.
• [11] D. Buccoliero and A. S. Desyatnikov. Quasi-periodic transformations of nonlocal spatial solitons. Opt. Express, 17:9608–9613, 2009.
• [12] F. Maucher, S. Skupin, M. Shen, and W. Krolikowski. Rotating three-dimensional solitons in Bose-Einstein condensates with gravitylike attractive nonlocal interaction. Phys. Rev. A, 81:063617, 2010.
• [13] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos. Observation of Dipole-Dipole Interaction in a Degenerate Quantum Gas. Phys. Rev. Lett., 95:150406, 2005.
• [14] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau. Bose-Einstein Condensation of Chromium. Phys. Rev. Lett., 94:160401, 2005.
• [15] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn. Bogoliubov modes of a dipolar condensate in a cylindrical trap. Phys. Rev. A, 74:013623, 2006.
• [16] D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell. Collective Excitations of a Bose-Einstein Condensate in a Dilute Gas. Phys. Rev. Lett., 77:420–423, 1996.
• [17] S. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle. Observation of Feshbach resonances in a Bose-Einstein condensate. Nature, 392:151–154, 1998.
• [18] M. D. Feit, J. A. Fleck Jr., and A. Steiger. Solution of the Schrödinger equation by a spectral method. J. Comput. Phys., 47:412–433, 1982.
• [19] W. E. Arnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue problem. Q. Appl. Math., 9:17–29, 1951.
• [20] R. B. Lehoucq, D. C. Sorensen, and C. Yang. ARPACK USERS GUIDE: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods, 1997.
• [21] J. L. Roberts, N. R. Claussen, S. L. Cornish, E. A. Donley, E. A. Cornell, and C. E. Wieman. Controlled Collapse of a Bose-Einstein Condensate. Phys. Rev. Lett., 86:4211–4214, 2001.
• [22] A. Gammal, T. Frederico, and L. Tomio. Critical number of atoms for attractive Bose-Einstein condensates with cylindrically symmetrical traps. Phys. Rev. A, 64:055602, 2001.
• [23] T. Koch, T. Lahaye, J. Metz, B. Frohlich, A. Griesmaier, and T. Pfau. Stabilization of a purely dipolar quantum gas against collapse. Nat. Phys., 5:218–222, 2008.
• [24] M. Edwards, P. A. Ruprecht, K. Burnett, R. J. Dodd, and C. W. Clark. Collective Excitations of Atomic Bose-Einstein Condensates. Phys. Rev. Lett., 77:1671–1674, 1996.
• [25] A. D. McLachlan. A variational solution of the time-dependent Schrodinger equation. Mol. Phys., 8:39–44, 1964.
• [26] P. A. M. Dirac. Note on Exchange Phenomena in the Thomas Atom. Math. Proc. Cambridge, 26:376–385, 1930.
• [27] H. Cartarius, T. Fabčič, J. Main, and G. Wunner. Dynamics and stability of Bose-Einstein condensates with attractive interaction. Phys. Rev. A, 78:013615, 2008.
• [28] F. Dalfovo, S. Giorgini, M. Guilleumas, L. Pitaevskii, and S. Stringari. Collective and single-particle excitations of a trapped Bose gas. Phys. Rev. A, 56:3840–3845, 1997.
• [29] T. Lahaye, J. Metz, B. Fröhlich, T. Koch, M. Meister, A. Griesmaier, T. Pfau, H. Saito, Y. Kawaguchi, and M. Ueda. -Wave Collapse and Explosion of a Dipolar Bose-Einstein Condensate. Phys. Rev. Lett., 101:080401, 2008.
• [30] P. Köberle and G. Wunner. Phonon instability and self-organized structures in multilayer stacks of confined dipolar Bose-Einstein condensates in optical lattices. Phys. Rev. A, 80:063601, 2009.
• [31] A. Junginger, J. Main, and G. Wunner. Variational calculations on multilayer stacks of dipolar Bose-Einstein condensates. Phys. Rev. A, 82:023602, 2010.
• [32] R. M. Wilson and J. L. Bohn. Emergent structure in a dipolar Bose gas in a one-dimensional lattice. Phys. Rev. A, 83:023623, 2011.
• [33] N. P. Proukakis and B. Jackson. Finite-temperature models of Bose–Einstein condensation. J. Phys. B, 41:203002, 2008.
• [34] M. Kreibich. Diploma thesis, Universität Stuttgart, 2011, available at http://itp1.uni-stuttgart.de/publikationen/abschlussarbeiten/kreibich_diplom_2011.pdf.
• [35] K. Góral and L. Santos. Ground state and elementary excitations of single and binary Bose-Einstein condensates of trapped dipolar gases. Phys. Rev. A, 66:023613, 2002.
• [36] I. S. Gradshteyn and I. M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, New York and London, 4th edition, 1965.
• [37] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions. United States Department of Commerce, 10th edition, 1972.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

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