# Variational calculations for anisotropic solitons in dipolar Bose-Einstein condensates

###### Abstract

We present variational calculations using a Gaussian trial function to calculate the ground state of the Gross-Pitaevskii equation and to describe the dynamics of the quasi-two-dimensional solitons in dipolar Bose-Einstein condensates. Furthermore we extend the ansatz to a linear superposition of Gaussians improving the results for the ground state to exact agreement with numerical grid calculations using imaginary time and split-operator method. We are able to give boundaries for the scattering length at which stable solitons may be observed in an experiment. By dynamical calculations with coupled Gaussians we are able to describe the rather complex behavior of the thermally excited solitons. The discovery of dynamically stabilized solitons indicates the existence of such BECs at experimentally accessible temperatures.

###### pacs:

03.75.Lm, 05.30.Jp, 05.45.-a## I Introduction

Since the prediction and the experimental realization of Bose-Einstein condensates (BECs) the field of cold atomic gases has been subject of multiple theoretical and experimental investigations. BECs with long-ranged interaction which can experimentally be realized by the condensation of atoms with a magnetic dipole moment such as are of special interest Griesmaier05a (); Beaufils08 (); Lahaye09 (). To stabilize such condensates, in general optical traps are applied in all three spatial directions. However, it has been shown by Tikhonenkov et al. Tikhonenkov08a () that stable quasi-2d solitons are possible where a trap is applied only in one direction perpendicular to the axis of the aligned atomic dipoles.

Solitons are a nonlinear effect which arises from the dispersion and the nonlinearity cancelling out each other. The solitons suffer from two kinds of instabilities: First, strong attractive particle interactions can cause the collapse of the condensate and second, when the interactions are only weakly attractive or even repulsive the BEC can dissolve in the directions where the external trap is open.

A prerequisite for the experimental realization of the solitons is a detailed theoretical investigation of the parameter ranges where the condensate is stable. The relevant parameters are the trap frequency, the scattering length of the contact interaction, and the excitation energy which allows for an estimation of the temperature range where the solitons can exist. The detailed analysis of the stationary states and the dynamics of solitons based on extended variational calculations is the objective of this article.

In the mean-field approximation the BEC with particle number is described by the extended Gross-Pitaevskii equation (GPE). Introducing “natural” units for mass , action , length , energy , frequency , abbreviating and making use of the scaling properties of the GPE leads to the scaled extended time-dependent GPE in “natural” units

(1) |

where the tilde is omitted. The dipole moments of the atoms are aligned in direction by an external magnetic field, and is the angle between and the magnetic field axis. In the trap geometry assumed here only a trap in the -direction perpendicular to the magnetic field is present, i.e. .

In Ref. Tikhonenkov08a () the GPE (1) has been solved approximately using a variational approach with a Gaussian type orbital and numerically exact by simulations on a grid. The grid calculations are numerically quite expensive. For a detailed analysis of the parameter space of the quasi-2d solitons we therefore introduce and employ an extended variational method based on coupled Gaussian functions, which has already turned out, for dipolar BECs with an axisymmetric three-dimensional trap, to be a full-fledged alternative to grid simulations Rau10 (); Rau10a (); Rau10b (). In this article the stable ground state of the condensate is computed by imaginary time evolution of an initial wave function, and it is shown that typically three to six coupled Gaussians are sufficient to obtain fully converged results. For a given trap frequency stable solitons exist in a finite range of the scattering length , outside that range the condensate collapses or dissolves.

The dynamics of energetically excited solitons is studied by solving the equations of motions for the variational parameters in real time. The mean-field energy of the ground state is typically only slightly below the energy threshold where the soliton can dissolve. However, the investigation of the dynamics reveals the existence of dynamically stabilized solitons at energies far above that threshold, indicating the possible experimental realization of solitons at temperature or even higher in the limit of the GPE. The transition temperature of a chromium BEC is experimentally given by Griesmaier05a (). Therefore the soliton can be dynamically stabilized in all such BECs.

The paper is organized as follows: In Sec. II we investigate the solitons with the ansatz of a single Gaussian wave packet. The appealing simplicity of this model is that both the stationary states and the dynamics of the solitons can be obtained from the Hamiltonian of a pseudo particle moving in a three-dimensional potential. In Sec. III the variational ansatz will be extended to a linear superposition of Gaussian wave packets (GWPs), and the TDVP will be applied to GWPs. The imaginary time evolution method will be used in Sec. IV for the evaluation of the ground state, and the results will be compared to the calculations with one Gaussian as well as with numerical grid calculations with the split-operator method. The analysis of the real time dynamics allows us to estimate the stability of excited solitons at finite temperatures. Concluding remarks are given in Sec. V.

## Ii Variational approach with a single Gaussian

Although the variational ansatz with a single Gaussian function cannot provide quantitatively correct results the simple model already allows us to gain deep insight into the physics of the quasi-2d solitons.

The ansatz with a single Gaussian as a trial function in the TDVP reads

(2) |

with the complex variational parameters

(3a) | ||||

(3b) |

The parameters determine the real half-widths of the Gaussian function, describes the amplitude used for normalization of the wave function, and is a global phase. The variational ansatz (2) of course strongly simplifies the problem. Nevertheless the investigation of this model yields physical insight and the results for the ground state are, as will be shown, qualitatively correct.

The mean-field energy and the chemical potential are

(4) | ||||

(5) |

respectively. The expectation values with the wave function given in Eq. (2) read

(6) |

for the kinetic term,

(7) |

for the trapping potential,

(8) |

for the scattering potential, and

(9) |

for the dipolar potential with and . The term

(10) |

denotes an elliptic integral of second kind in “Carlson’s form” which can be evaluated with a fast and stable approximation algorithm Carlson94a (); Carlson95a () more efficiently than the integral representation given in Tikhonenkov08a ().

### ii.1 Hamiltonian form of the equations of motion

The ground state of the GPE can in principle be obtained by minimizing the mean-field energy in Eq. (4). Alternatively, the TDVP can be applied to derive equations of motion for the variational parameters in Eq. (3). The equations of motion can be used to investigate the dynamics of the soliton, and the stationary states of the system can be found by searching for the fixed points of these equations. The stable fixed point with the lowest mean-field energy denotes the ground state.

For the ansatz of a single Gaussian the system can be transformed to the descriptive form of a Hamiltonian system by the coordinate transformation

(11) |

The Hamiltonian in the coordinates

(12) | ||||

has the conventional form . It can be shown that Hamilton’s equations

(13) |

lead to the same equations of motion as the TDVP applied to the variational parameters in Eq. (3).

The potential is visualized in Fig. 1.

The isopotential surfaces close to the stable fixed point have ellipsoidal form. At the trap energy the potential gets open. The unstable fixed point at smaller values of than the local minimum is given by the saddle point where the surrounding isosurfaces have hyperbolic form. In Fig. 1(c) the potential at a smaller scattering length close to the bifurcation point is shown. The local minimum and the saddle point approach one another and coincide at the critical scattering length .

### ii.2 Stationary states and linear stability

The fixed points of the system can now be easily calculated by a nonlinear root search for and in Hamilton’s equations (13). In Fig. 2 the mean-field energy and the chemical potential of the stable ground state and of an excited unstable state are shown as a function of the scattering length.

The two branches emerge in a tangent bifurcation at the critical scattering length . In Fig. 2(c) the width parameters corresponding to the ground state show, that with increasing scattering length the soliton gets very broad and finally dissolves.

To analyze the stability of the solution of the GPE the eigenvalues of the six-dimensional Jacobi matrix

(14) |

or equivalently in the canonical coordinates the eigenvalues of the three-dimensional real symmetric Hesse matrix of the potential are calculated. The results are shown in Fig. 3.

The eigenvalues appear in pairs of different sign and are all pure imaginary for the stable state. The linear stability for different trap frequencies shows qualitatively similar behavior.

### ii.3 Dynamics with a frozen Gaussian

Hamilton’s equations (13) describe the dynamics of the soliton with the potential . The systematic investigation and visualization of the dynamics of a Hamiltonian system with three degrees of freedom is a nontrivial task. However, the force in the direction is dominated by the strong harmonic trap potential. If therefore is sufficiently large at least for the ground state of the condensate takes nearly the value of the harmonic oscillator ground state. As a further simplification of the problem we therefore restrict the dynamics to the plane given by the condition

(15) |

which corresponds to a frozen Gaussian ansatz in the direction. For the plane is marked in Fig. 1(b).

The dynamics of the wave function (2) with fixed parameters , is described in the canonical coordinates and by the two-dimensional potential

(16) | |||

which is illustrated in Fig. 4.

To describe the dynamics of this two-dimensional system it is convenient to analyze Poincaré surfaces of section (PSOS). An adequate choice of the PSOS shown in Fig. 4 is using the rotated coordinates and momenta

(17a) | ||||

(17b) |

with and the crossing condition . For constant energy different initial conditions are integrated. The border of the allowed energy range is given by

(18) |

In Fig. 5(a) a PSOS at an energy close to the ground state is shown. The motion is completely regular and an elliptic fixed point of the Poincaré map is present, which belongs to the antisymmetric periodic oscillation of the condensate. The symmetric oscillation is not visible in the PSOS for the surface of section being almost parallel to it. Increasing the energy towards , the turning points of the periodic oscillations in Fig. 4 move to larger values of and (for the symmetric oscillation only one turning point, for the antisymmetric both). For they lie at infinity, i.e. the soliton dissolves. In the Poincaré map the resonant tori decay according to the Poincaré-Birkhoff therorem building the same number of elliptic and hyperbolic fixed points in the Poincaré map. The closer to the bifurcation point the scattering length is, the lower the energy is, where this takes place. This can be seen especially in Fig. 5. Further increasing of the energy leads to regions of chaotic oscillations while the elliptic fixed point moves outwards to larger values of .

In the picture with a frozen Gaussian the solitons dissolve at energy in agreement with Tikhonenkov08a (); Santos09 (). If this were always true it would imply that in an experiment with realistic parameters (e.g., a condensate with 10000 particles at and ) the solitons must be cooled down to temperatures of about because the energy gap between the ground state and the threshold is very small (for a more detailed discussion see Sec. IV.1). However, a dynamical stabilization of the solitons at energies is possible as will be discussed for an ansatz with a single Gaussian in the next Section II.4 and for coupled Gaussians in Sec. IV.2.

### ii.4 Three-dimensional dynamics

In the frozen Gaussian approximation any excitation energy of the soliton must be completely deposited in the motion, which leads to the dissolving of the soliton at the low threshold . In this section we demonstrate that for the fully three-dimensional dynamics with the ansatz (2) a large amount of excitation energy can be stored in the motion dominated by the one-dimensional harmonic trap.

An example trajectory for trap frequency , scattering length , and with mean-field energy far above the threshold at is presented in Fig. 6.

In (a) and (b) the expectation values

(19) |

are drawn as functions of time. The slow oscillations in and generate the Lissajous type motion of the quasi-2d soliton visualized in Fig. 6(c). The fast oscillation in results dominantly from the external harmonic trap. The calculations with a single Gaussian thus indicate that a rather highly excited soliton may be dynamically stabilized and does not dissolve in contrast to the calculations with a frozen Gaussian Tikhonenkov08a (); Santos09 (). However, it must be clarified (in Sec. IV.2) whether the dynamical stabilization is also possible with coupled Gaussians.

## Iii Ansatz with coupled Gaussians

Though the ansatz made in Sec. II offers a descriptive analysis of the soliton, comparison with calculations where the GPE is solved numerically on a grid show that the results for the ground state only hold qualitatively. As it will be shown, the results for the dynamics of the soliton are qualitatively different as well. To gain quantitatively correct results the ansatz for the trial function has to be modified. As shown in Rau10a (); Rau10b () the condensate wavefunction can be well described by a superposition of Gaussians. We therefore apply the TDVP to coupled Gaussian wave packets McLachlan1964a (); Heller75a (); Fab09a (); Fab09b (). The ground state will be determined by the imaginary time evolution (ITE) with the norm conservation embedded in the TDVP as a constraint.

In this section we elaborate the theory and derive the formulae for the computations with coupled Gaussians. The results for the stationary states and the dynamics of the solitons are presented in Sec. IV.

### iii.1 Time-dependent variational principle for GWPs

In this article the TDVP provides the basis for the solution of the time-dependent GPE. Its application to GWPs was originally introduced by Heller heller:4979 (); Heller81a () for the description of atomic and molecular quantum dynamics and later on applied to BEC Rau10a (); Rau10b (). For the convenience of the reader we here shortly recapitulate this method.

The quantity

(20) |

is to be minimized with respect to . Afterwards it is set . The wave function shall be parameterized by the variational parameters . The variation of in Eq. (20) carries over to the variation of and and results in

(21) |

for the wave function being constant. Because the variational parameters are complex quantities both brackets in Eq. (21) have to vanish separately. This yields

(22) |

which can be written shortly as

(23) |

where is an hermitian positive definite matrix. The linear system (23) has to be solved for every time step to integrate the equations of motion .

We now choose a linear superposition of Gaussians as the trial function

(24) |

where are complex diagonal matrices of dimension and denotes the relative phases and the amplitude of the Gaussians, respectively. The Gaussians are fixed at the origin. The time evolution can be considered as the motion in an effective time-dependent harmonic potential

(25) |

The dynamics of the GWPs are now determined by the TDVP, which variationally fits the effective time-dependent harmonic potential coefficients to the underlying potential. Splitting the Hamiltonian and operating on the trial wave function (24) yields

(26) |

The parameters are scalars and the matrices are diagonal matrices. The equations of motion then read

(27a) | ||||

(27b) |

The coefficients and have to be calculated from the TDVP

(28) |

(29) |

Combining the coefficients and in a complex vector the set of equations (29) can be written as

(30) |

The positive definite hermitian matrix includes in contrast to Eq. (23) the kinetic operator. Eq. (30) has to be solved for every time step and the resulting vector has then to be inserted in Eq. (27).

In the numerical integration of the equations of motion (27) (e.g. with a standard algorithm like Runge-Kutta), the quadratic term can lead to the numerical overflow heller:4979 (). It is therefore reasonable to split up the matrices into two matrices and according to . For the case of diagonal matrices , and are diagonal, too. Hence the splitting can be done component-by-component. This leads to

(31) |

where and are taken as initial values.

### iii.2 Constraints in the TDVP

The most important reason for the introduction of constraints is the conservation of the norm in the imaginary time evolution (cf. Sec. III.4). Apart from that inequality constraints can be used to avoid matrix singularities which arise from overcrowding the basis set Fabcic08a (). In principle constraints can be introduced combined in the -dimensional vector . The constraints are embedded then by a set of Lagrangian multiplicators

(32) |

where is a real matrix. To find the minimum of

(33) |

has to be calculated which yields the system of equations

(34) |

with

In the special case of Gaussian wave packets the equations of motion (split into real and imaginary part and combining the variational parameters to can be written in the form

(35) |

where consists of the terms linear in and the constant terms are absorbed in the vector . The constraints combined in yield

(36) |

The set of linear equations (34) then reads

(37) |

with

(38) |

Solving this system of linear equations and inserting the result vector in Eq. (27) enables one to integrate the equations of motion.

### iii.3 Energy functional

To gain the mean-field energy of the soliton the time-dependent variational parameters of the GWPs are used in the evaluation of the integrals in the energy functional. The energy functional for coupled Gaussians reads

(39) |

With the integrals in (39) yield

(40) | |||

(41) |

(42) |

(43) |

with the abbreviations

and

The additional integrals in the TDVP , and are obtained easily from the integrals in the energy functional.

### iii.4 Imaginary time evolution

In principle the ground state of the system of coupled Gaussians could be obtained by a nonlinear root search in the equations of motion. Since complex initial values are required, it is difficult to find the fixed points this way for a large number of Gaussians . The imaginary time evolution represents an alternative for the determination of the ground state.

In the linear Schrödinger equation the transformation leads to an exponential damping of the excited states and converges to the ground state for . The method can be applied to the nonlinear GPE as well. In imaginary time the norm is not conserved any more. This is crucial in the nonlinear GPE because the decay of the norm has the meaning of particle losses. To avoid the norm decay we implement the norm conservation as a constraint

(44) |

in the TDVP as elaborated above. Supposing a general complex rotation of the time with the equations of motion can be written in the general form

(45) |

with

where the abreviative notation in can be understood as follows: Every entry is an submatrix . If all elements of , otherwise . The vectors , , and contain the real and imaginary parts of the coefficients of the effective potential. The entries in are vectors with the index running from to . One special case is the real time evolution , and the other special case the imaginary time evolution , . The gradient vector in Eq. (36) is given by

(46a) | ||||

(46b) | ||||

(46c) | ||||

(46d) |

with which one obtains and .

With this method the norm is conserved during the numerical integration for rather long times. However for very long times a small drift in the norm is present due to the numerical error of the integration. The ITE of a system of 10 coupled Gaussians is shown in Fig. 7.

If the relative accuracy of the integrater is chosen coarse, the mean-field energy does not decay monotonously but for long times it converges to the same value as an integration with a fine relative accuracy thus a coarse accuracy can be chosen for the sake of time. The ITE is very robust to the initial choice of the wave function. This also holds for large numbers of Gaussians coupled, where the nonlinear root search often fails. However, only the stable ground state is accessible with this method.

## Iv Results with coupled Gaussians

Applying the variational ansatz with coupled Gaussians to the wave functions of dipolar condensates we are now able to obtain quantitatively correct ranges of the parameters where stable quasi-2d solitons can exist. For an experimental realization of solitons with chromium atoms realistic trap frequencies are typically a few hundred Hertz for condensates with about to particles. These values roughly correspond to an interval of . We therefore focus our calculations on two values for the scaled trap frequency, viz. and .

### iv.1 Stationary ground state

In Figs. 8 and 9 the results for the ground state of the soliton are shown at trap frequencies and , respectively.

In Figs. 8(a) and 9(a) the mean-field energy as a function of the scattering length is plotted for up to 6 coupled Gaussians. For comparison the calculation using a nonlinear root search for one Gaussian as discussed in Sec. II.2 is displayed, as well. The mean-field energy and also the chemical potential [see Figs. 8(b) and 9(b)] converge with increasing number of Gaussians . The detailed convergence properties of the calculations with increasing number of Gaussians are illustrated in Fig. 10, where calculations for up to 12 Gaussians are shown at constant scattering length .

The ground state energy converges fast with the number of Gaussians and the corrections above a number of about 5 Gaussians can be neglected. The small energy differences between the variational and the grid calculations in Fig. 9(a) may be due to the finite grid size, which is limited by an acceptable computation time.

In the calculations with a single Gaussian a stable and an unstable state are created in a tangent bifurcation at a scattering length . Using coupled Gaussians the critical scattering length where the condensate collapses is shifted to higher values, e.g., from to at trap frequency and from to at trap frequency . The ITE can only provide the stable ground state but no unstable states which are certainly also involved in the bifurcation. For dipolar condensates in an axisymmetric trap a complicated bifurcation scenario has been revealed Rau10 (); Rau10b (), and it will be an interesting future task to study bifurcations of the soliton states in more detail.

The corrections to the calculations with a single Gaussian decrease with growing scattering lengths. This can be understood when looking at the spatial size of the soliton [cf. Fig. 2(c)]. For large expansion of the condensate the nonlinear interaction terms in the GPE are small and thus the mean-field energy is dominated by the ground state energy of the harmonic trap in the direction. Above a certain threshold value for the scattering length the soliton is no longer bound but dissolves. The thresholds at, e.g., for and for obtained with a single Gaussion (see Fig. 2) are only very slightly shifted to higher values when using coupled Gaussians.

The extended variational approach allows us to compute for each trap frequency an accurate lower and upper critical value of the scattering length, i.e., a range where stable solitons in dipolar BECs can exist. However, in an experiment thermal excitations of the soliton may cause the dissolving of the soliton, and therefore it is also necessary to determine an energy or temperature limit for the existence of stable solitons. As a dissolving condensate without any interaction must have at least the zero point energy of the harmonic trap in the direction the difference between this threshold and the mean-field energy of the ground state,

(47) |

can be taken for a very conservative estimation of the temperature at which the soliton should survive thermal excitations. In Figs. 8(c) and 9(c) that temperature is plotted as a function of the scattering length for two condensates with and particles. With this conservative estimate a soliton at, e.g., trap frequency and scattering length must be cooled down to about . However, solitons may exist at much higher temperatures when they are dynamically stabilized as discussed in Sec. II.4 for the ansatz with a single Gaussian. We now extend that discussion to the ansatz with coupled Gaussians.

### iv.2 Dynamics with coupled Gaussians

The ansatz with a single Gaussian used to describe the dynamics of the solitons in Sec. II.4 and especially the further simplification with a frozen Gaussian in Sec. II.3 are basically mathematical model systems, and the results obtained with these models have to be seen from an academic point of view. For a realistic description it is necessary to investigate the dynamics of the solitons with the extended ansatz (24) of coupled Gaussians which has shown to significantly improve the results for the ground state in Sec. IV.1. Of special interest is to search for the existence of dynamically stabilized solitons at energies above the threshold . To this aim the equations of motion (27) are integrated with the initial wave function chosen appropriately to yield the desired mean-field energy.

The condensate wave function (24) with coupled Gaussians is parametrized by real time-dependent variational parameters. For the visualization of the wave function we use the reduced set of parameters

(48) |

which are related to the extension of the soliton in the three spatial dimensions and can be directly compared with the results for a single Gaussian using the coordinates in Eq. (19).

For energies close to the ground state energy periodic oscillations can be found. The parameters of the turning points of the periodic oscillations show a rather complicated behavior when the mean-field energy is increased. The different periodic oscillations vanish in multiple bifurcations close to the ground state energy. Above no periodic oscillations were found. However increasing the energy above yields oscillations of the soliton which do not destroy it. Applying a frozen Gaussian approximation yields stable oscillations only slightly above this limit. But allowing oscillations in all spatial directions stabilizes the condensate and gives rise to an oscillating soliton as shown in Fig. 11 for a trajectory at trap frequency , scattering length , and mean-field energy .