Bifurcations, order, and chaos in the Bose-Einstein condensation of dipolar gases

We apply a variational technique to solve the time-dependent Gross-Pitaevskii equation for Bose-Einstein condensates in which an additional dipole-dipole interaction between the atoms is present with the goal of modelling the dynamics of such condensates. We show that universal stability thresholds for the collapse of the condensates correspond to bifurcation points where always two stationary solutions of the Gross-Pitaevskii equation disappear in a tangent bifurcation, one dynamically stable and the other unstable. We point out that the thresholds also correspond to"exceptional points,"i.e. branching singularities of the Hamiltonian. We analyse the dynamics of excited condensate wave functions via Poincare surfaces of section for the condensate parameters and find both regular and chaotic motion, corresponding to (quasi-) periodically oscillating and irregularly fluctuating condensates, respectively. Stable islands are found to persist up to energies well above the saddle point of the mean-field energy, alongside with collapsing modes. The results are applicable when the shape of the condensate is axisymmetric.

degenerate quantum gases with adjustable long-range and short-range interactions. The achievement of Bose-Einstein condensation in a gas of chromium atoms [14], with a large dipole moment, has in fact opened the way to promising experiments on dipolar gases [15], which could show a wealth of novel phenomena [16][17][18][19]. In particular, the experimental observation of the collapse of dipolar quantum gases has also been reported [20] which occurs when the contact interaction is reduced, for a given particle number, below some critical value using a Feshbach resonance.
In this experimental situation it is most timely and appropriate to extend the investigations of the effects of the nonlinearity of the Gross-Pitaevskii equation to dipolar quantum gases, and this is the goal of the present paper. To model these effects we will pursue, for the sake of simplicity, a variational ansatz. We do this in the spirit of Refs. [3,4,21] where also variational techniques were applied to model the dynamics of dilute ultracold atom clouds in the Bose-Einstein condensed phase by solving the Gross-Pitaevskii equation without dipole-dipole interaction. In fact, quite recently Parker et al. [22] have pointed out that in dipolar Bose-Einstein condensates the Gaussian variational method gives excellent agreement with the full numerical solutions of the Gross-Pitaevskii equation in wide ranges of the physical parameters. For oblate dipolar Bose-Einstein condensates the Gaussian approximation appears to agree only for weak dipolar interactions. The approximation used may not be valid in the limit of interest and clearly further exact studies based on exact solutions of the dipolar Gross-Pitaevskii equation need to be carried out to verify the results. Full numerical quantum calculations, for condensates with only the contact interaction [3,4,21] present, or with an additional attractive gravity-like 1/r interaction [23][24][25][26], have confirmed that properties of the solutions of the Gross-Pitaevskii equation found in the variational calculations are recovered in the quantum calculations. We should therefore expect that a simple variational approach will also capture essential features of the dynamics of condensate wave functions of dipolar gases.
We treat the problem in the "atomic" units provided by the magnetic dipole-dipole interaction, i.e., we measure lengths in units of the "dipole length" energies in units of E d =h 2 /(2ma 2 d ), frequencies in units of ω d = E d /h and time in units ofh/E d , respectively. In (1), α is the fine-structure constant, a 0 the Bohr radius, m/m e the ratio of the atom and electron mass, and µ the magnetic moment. For 52 Cr, with µ = 6µ B (µ B the Bohr magneton), one has a d = 91 a 0 , E d = 1.7 × 10 −8 eV, and ω d /2π = 4.2 × 10 6 Hz.
In these "atomic" units, the Hartree equation of the ground state of a system of N identical bosons in an external trapping potential V trap = m(ω 2 r r 2 + ω 2 z z 2 )/2, all in the same single-particle orbital ψ, interacting via the contact interaction and the magnetic dipole-dipole interaction, assumes the dimensionless form Here, ε is the chemical potential, a is the scattering length, and γ r,z = ω r,z /(2ω d ) are the dimensionless trap frequencies. Setting Ψ = √ N ψ, one recovers the familiar form of the time-independent Gross-Pitaevskii equation for dipolar quantum gases.
In the above "atomic" units, the mean-field Hamiltonian in (2) obeys a scaling law with respect to the number N of atoms: Letψ(r) be a solution of the (formal) one-boson problem for a given scattering length a/a d and trap frequenciesγ r,z , then ψ(r) := N −3/2ψ (r), with r = Nr, solves the N-boson problem for the same scattering length a/a d , but with trap frequencies γ r,z =γ r,z /N 2 and chemical potential ε =ε/N 2 . Note that the particle number scaling leaves the aspect ratio λ = γ z /γ r invariant. Thus, the physical properties of Bose condensates of dipolar quantum gases quite generally only depend on the value of the scattering length a/a d and the particle number scaled trap frequencies N 2 γ r,z or, alternatively, on the aspect ratio λ and the scaled geometric mean of the trap frequencies N 2γ = N 2 γ 2/3 r γ 1/3 z . We note that the dimensionless parameter D introduced by Dutta and Meystre [27] and Ronen et al. [28] to measure the effective strength of the dipole interaction in trap frequency units is related to our scaling parameters by D = (N 2 γ r /2) 1/2 .
As an application of the particle number scaling law we emphasise that the experimental results reported by Koch et al. [20] for the stabilisation of dipolar chromium quantum gases for particle numbers ∼ 20.000 and a trap frequencyω/2π = 720 Hz correspond to a value of the scaled trap frequency N 2γ = 3.4 × 10 4 . Therefore they directly carry over to any pairs of particle numbers and mean trap frequencies with this value of the scaled trap frequency and the same aspect ratios! To study the nonlinearity effects of the time-independent Gross-Pitaevskii equation for dipolar gases we adopt the familiar variational ansatz of a (normalised) Gaussian type orbital (e. g. [9,18,20]) and exploit the time-dependent variational principle for the mean-field energy to determine the width parameters A r and A z . It is well known that solutions can be found only in certain parameter ranges and that at critical values the condensate collapses.
What -to the best of our knowledge -for dipolar quantum gases has gone unnoticed before is that at these stability thresholds actually two solutions of the extended Gross-Pitaevskii equation (2) disappear. The situation is depicted in figure 1 where the chemical potential is plotted as a function of the scattering length for different values of the trap aspect ratio λ and the fixed value of the scaled geometric mean of the trap frequency of N 2γ = 3.4 × 10 4 . It can be seen that, as the scattering length is increased, . Stability thresholds, i.e. the critical scattering lengths a crit /a d , as a function of the scaled trap parameter N 2γ and the aspect ratio λ. In the limit N 2γ → ∞, λ → 0, one has a crit /a d = 1/6.
for every value of λ two solutions are born in a tangent bifurcation, one corresponding to the ground state and the other to a collectively excited state. An inspection of the mean-field energies shows that the excited state corresponds to the branch of the chemical potential which diverges for a/a d → 1/6. The behaviour of the stability thresholds over a wide region of the parameter space spanned by N 2γ and λ is shown in figure 2. Large values of N 2γ correspond to the regime where the dipole-dipole interaction is dominant, as in the experiments of Koch et al. [20].
For N 2γ < 1 the contact interaction prevails, and the dipole-dipole interaction is only a small perturbation.
In what follows we shall investigate the structure of the bifurcation, the stability and the dynamics of dipolar condensates. The first aspect we want to highlight is that dipolar Bose-Einstein condensates near the stability thresholds are experimental realizations of physical systems with "exceptional points", a phenomenon hitherto observed only in open quantum systems, described by non-Hermitian Hamiltonians ( [29,30], see also [31] and references therein). Exceptional points are positions in the parameter space where both the energy eigenvalues and the wave functions of (usually) two eigenstates pass through a branch point singularity as functions of the parameters and, consequently, are identical at the critical set of parameters. While in linear Schrödinger equations two eigenstates can only become identical if the Hamiltonian is non-Hermitian, two coalescing eigenstates can also occur for a nonlinear Schrödinger equation as an effect of its nonlinearity.
This indeed the case for the stationary solutions of the Gross-Pitaevskii equation with dipole-dipole interaction. At the stability thresholds the energy eigenvalues and the corresponding wave functions of the two bifurcating states pass through a branch point singularity of the Hamiltonian, the energies N 2 ε and the corresponding wave functions are identical. The way to reveal the branch point singularity structure is to continue the scattering length a/a d into the complex plane and to check a well-known property of exceptional points: If in traversing a full circle around the critical value a crit /a d at the point of bifurcation, a/a d = (a crit /a d ) + ̺ e iϕ , ϕ = 0 . . . 2π, a permutation of the two solutions occurs, an exceptional point is located within the circular area [29].
For the parameters of the Koch et al. [20] experiment we have performed such an analysis, and for the case of an almost purely dipolar quantum gas the results are shown in figure 3. As one goes around a small circle in the complex plane with the critical scattering length at the centre, the two solutions are permuted, i.e. we have an exceptional point. It must be said, however, that the usual way of experimentally proving the occurrence of an exceptional point in open quantum systems, changing two real physical parameters to traverse a circle in the complex energy plane, cannot be applied here. In the case of the Gross-Pitaevskii equation it is by continuing one real parameter, the scattering length, into the complex plane that the nature of the stability thresholds as exceptional points is revealed.
We now analyse the dynamics of the condensate wave functions. To do so we start from the time-dependent Gross-Pitaevskii equation, which is (2) with the replacement ε → i d dt , and generalise the ansatz (5) to a time-dependent Gaussian type orbital with complex width parameters A r (t) and A z (t) and complex normalisation factor A(t) (cf. [21]). We can then apply the time-dependent variational principle iψ(t)−Hψ(t) 2 ! = min. [32,33] to derive a system of ordinary nonlinear differential equations for the time evolution of the real and imaginary parts of the variational parameters A r (t) and A z (t). These can be shown to be equivalent to canonical equations of motion belonging to a two-dimensional nonintegrable autonomous Hamiltonian system. For the time evolution one can therefore expect all the features familiar from nonlinear dynamics studies of such systems, including a transition to chaos. For brevity, the details of these calculations are relegated to a subsequent paper. Here we shall concentrate on the results. We first investigate the stability of the two independent solutions of the stationary Gross-Pitaevskii equation which emerge from the bifurcations. To this end the four equations of motion for the real and imaginary parts of A r and A z are linearised around the values corresponding to the stationary states. In this way for each of the two states we obtain four eigensolutions ψ (lin) ∼ e κt of the linearised system of equations, with eigenvalues κ. For the ground state all eigenvalues turn out to be purely imaginary, proving that the state is indeed dynamically stable. For the collectively excited state one also finds a positive real eigenvalue, and hence the state is dynamically unstable. Similar behaviour was found by Huepe et al. [3,4] in their study of the stability of bifurcating solutions with only an attractive contact interaction present.
The system of nonlinear first-order differential equations also serves to investigate the time evolution of any initial state of the condensate by following the corresponding trajectories in the four-dimensional configuration space spanned by the coordinates of the real and imaginary parts of A r and A z . Since the total mean-field energy is a constant of motion the trajectories are restricted to three-dimensional hyperplanes, and their behaviour can most conveniently be visualised by two-dimensional Poincaré surfaces of section defined by requiring one of the coordinates to assume a fixed value.
We consider Poincaré surfaces of section defined by the condition that the imaginary part of A z (t) is zero. Each time the trajectory crosses the plane Im(A z ) = 0, the real and imaginary parts of A r (t) = A r r (t)+iA i r (t) are recorded. In figure 4 surfaces of section are plotted for seven different, increasing, values of the mean-field energy. Again the physical parameters of the experiment of Koch et al. [20] are adopted, and the scattering length is fixed to a/a d = 0.1, away from its critical value. At these parameters, the variational mean-field energy of the ground state is NE gs = 4.24 × 10 5 and represents the local minimum on the two-dimensional mean-field energy landscape, plotted as a function of the (real) width parameters. The variational energy of the second, unstable, stationary . Poincaré surfaces of section of the condensate wave functions represented by their width parameters for seven different mean-field energies at the scaled trap frequency N 2γ = 3.4 × 10 4 , aspect ratio λ = 6, and the scattering length a/a d = 0.1. The surfaces of section labelled 1 to 7 correspond to the increasing values of the meanfield energy N E = 4.24 × 10 5 , 5.00 × 10 5 , 6.00 × 10 5 , 9.00 × 10 5 , 1.50 × 10 6 , 3.00 × 10 6 , 6.00 × 10 6 , respectively. state at these experimental parameters is NE es = 6.24×10 5 , it corresponds to the saddle point on the mean-field energy surface. Between these two energy values the motion on the trajectories is bound, while for energies above the saddle point energy the motion on the trajectories can become unbound: once the saddle point is traversed by a trajectory A r (t), A z (t), the parameters run to infinity, A r r (t), A r z (t) → ∞, meaning a shrinking of the quantum state to vanishing width, i.e. a collapse of the condensate takes place.
The surface of section labelled 1 in figure 4 belongs to the energy of the ground state. At this energy the kinematically allowed region for the crossing points of the trajectories is confined to a single stable stationary point. At the next higher energy (surface of section labelled 2) the kinematically accessible region in configuration space has grown, and the initially stationary state has evolved into a periodic orbit (fixed point in the surface of section), corresponding to a state of the condensate whose motion is periodic. The oscillations of the width parameters A r (t) and A z (t) represent oscillatory stretchings of the condensate along the r and z directions. The stable periodic orbit in the surface of section is surrounded by elliptical, quasi-periodic orbits, representing quasi-periodic oscillations of the condensate. The surface of section 3, at the next higher energy, reveals that bifurcations have occurred, creating new stable and unstable periodic states, manifested by the emergence of additional elliptical islands and separatrices in the surface of section.
The surface of section labelled 4 is the first in figure 4 with an energy value above the saddle point energy. Now chaotic orbits have appeared which surround the stable regions. In contrast to the (quasi-) periodic stretching oscillations of the condensate within the elliptical islands, the chaotic motion of the parameters describes a condensate which does not yet collapse but whose widths fluctuate irregularly.
One might imagine that well above the saddle point energy stable condensate wave functions no longer exist. However, in the surfaces of section labelled 5, 6, and 7 regular islands are still clearly visible. These stable islands are surrounded by chaotic trajectories. Since ergodic motion along these trajectories comes close to every point in the configuration space, the chaotic motion sooner or later leads to a crossing of the saddle point and then to the collapse of the condensate wave functions. It can be seen that with growing energy above the saddle point the sizes of the stable regions gradually shrink. The reason why the kinematically allowed regions surrounding the stable islands are hardly recognisable any more in these surfaces of section is that high above the saddle point energy the chaotic motion becomes more and more unbound, which means that the trajectories cross the Poincaré surfaces of section only a few times, if ever, before they escape to infinity and collapse takes place.
It must be emphasised, however, that stable islands do persist even far above the saddle point energy, implying the existence of quasi-periodically oscillating nondecaying modes of the condensate wave functions.
The prediction of such modes of energetically excited solutions of the timedependent Gross-Pitaevskii equation for cold dipolar quantum gases is a result of our analysis. It would certainly be an intriguing and challenging task to examine whether it is possible to prepare excited states of dipolar quantum gases of this type in both the regular as well as in the chaotic regions, to distinguish between the two different dynamics, and to access the stable regions high above the energy of the unstable stationary state. One way of creating the collectively excited states one might imagine is to prepare the condensate in the ground state, and then to non-adiabatically reduce the trap frequencies. Clearly experimental investigations along these lines are strongly encouraged.
One might ask, however, whether the Gross-Pitaevskii equation underlying our calculations is adequate at all to describe complex dynamics of this type in real dipolar quantum gases. In particular in the chaotic regime local density maxima might occur for which losses by two-body or three-body collisions would have to be taken into account. However, by virtue of the scaling law (4) parameter ranges can always be found where the particle densities remain small even in these regimes and the Gross-Pitaevskii equation is applicable.
In this paper we have used a simple variational ansatz (5) to model the dynamics of axisymmetric solutions of the Gross-Pitaevskii equation of dipolar gases. The advantage of the variational technique is that the analysis of the nonlinearity effects becomes especially transparent. Of course the results must be checked by accurate numerical simulations. Such simulations, and the extensions to fully three-dimensional and structured condensate wave functions [27,34], are under way. Comparisons of variational and accurate numerical results [24][25][26] in an alternative system with a long-range interaction, viz. Bose condensates with an attractive gravity-like 1/r interaction [23], have shown that the nonlinear dynamical properties found in the variational calculation are recovered in the numerical calculations, and that the quantum behaviour may be even richer. As a common feature of the propagation of quantum wave functions, the real-time evolution of arbitrary condensate wave functions of dipolar gases will exhibit complicated fluctuations. An interpretation of their full dynamics, and in particular the search for periodic, quasiperiodic or chaotic structures, therefore will hardly be possible without the guidance of the variational results obtained in this paper.