Variational approaches to dipolar Bose-Einstein
condensates
D I S S E R T A T I O N
Von der Fakult¨at Mathematik und Physik der Universit¨at Stuttgart zur Erlangung der W¨urde eines Doktors der Naturwissenschaften (Dr. rer. nat.) genehmigte Abhandlung
Vorgelegt von
R¨
udiger Matthias Fortanier
aus Stuttgart
Hauptberichter:
Prof. Dr. J¨
org Main
Mitberichter:
Prof. Dr. Hans Peter B¨
uchler
Tag der m¨
undlichen Pr¨
ufung: 17.04.2014
1. Institut f¨
ur Theoretische Physik der Universit¨
at Stuttgart
2014
Contents
Abstract 7 Inhaltsangabe 9 1. Introduction 13 1.1. Preface . . . 13 1.2. Outline . . . 16 2. Bose-Einstein condensation 19 2.1. The ideal Bose-gas . . . 192.2. Fundamental description of ultra-cold atomic gases . . . 20
2.3. The Gross-Pitaevskii equation . . . 22
2.3.1. Contact interaction . . . 22
2.3.2. Dipole-dipole interaction . . . 23
2.3.3. Dipolar Gross-Pitaevskii equation . . . 24
2.3.4. Properties and consequences . . . 25
2.4. Discussion of common solution methods . . . 27
2.4.1. Stationary states and real-time evolution . . . 27
2.4.2. Implementation of the grid calculations . . . 29
2.4.3. Stability analysis and excitation spectrum . . . 29
3. The time-dependent variational principle 31 3.1. The McLachlan variational principle . . . 32
3.1.1. Gaussian wave packets . . . 34
3.2. Analytical integrals for the equations of motion . . . 39
3.2.1. Norm integral and K-matrix . . . 40
3.2.2. Computation of the r-vector . . . 41
3.3. Dipolar potential . . . 42
3.3.1. Numerical evaluation of Jklij σ . . . 51
3.3.2. Cases with analytic solutions for Jklij σ . . . 53
3.3.3. Optimizations . . . 54
3.4. Techniques relying on the EOM . . . 56
3.4.1. Dynamical simulations . . . 56
3.4.3. Stability of stationary states . . . 59
3.5. Implementation . . . 59
3.5.1. Graphical user interface . . . 59
4. Two-dimensional solitons 63 4.1. Natural units . . . 63
4.2. Important symmetries of the Gross-Pitaevskii equation . . . 64
4.2.1. Galilei invariance . . . 64
4.2.2. Parity . . . 65
4.3. Stationary states . . . 65
4.4. Stability and dynamics . . . 67
4.4.1. Soliton without translational and rotational degrees of freedom 67 4.4.2. Soliton with rotational degrees of freedom . . . 68
4.5. Collisions of solitons . . . 69
4.5.1. Procedure . . . 70
4.5.2. Collisions with k = 10 . . . 71
4.5.3. Collisions of solitons with various initial momentum . . . 73
4.5.4. Collisions including finite angular momentum . . . 76
5. Dipolar condensates in a triple-well potential 81 5.1. External potential and units . . . 82
5.2. Method . . . 84
5.3. Results for the repulsive configuration . . . 85
5.3.1. Cut at N add = 0.6 . . . 87
5.3.2. Cut at N add = 0.2 . . . 90
6. Condensates in PT-symmetric potentials 97 6.1. PT symmetry . . . 97
6.2. BECs in a PT-symmetric double well without dipolar interaction . 99 6.3. Dipolar BECs in a PT-symmetric double-well potential . . . 100
6.3.1. Results for the repulsive configuration . . . 103
6.3.2. Results for the attractive configuration . . . 108
6.3.3. Comparison with the non-dipolar case . . . 108
7. Conclusion 111 7.1. Outlook . . . 113
A. Supplementary material for LSE construction 115 A.1. Supplementary material for derivatives of I0 . . . 115
A.2. Supplementary material for the dipolar integrals . . . 116
Contents
A.2.2. Quadratures . . . 118 A.2.3. Recursion formula for the terms of the Taylor series . . . 120 A.2.4. Pad´e approximation and -Wynn algorithm . . . 123 A.2.5. Additional dipolar integrals reducible to elliptic integrals . . 123 A.3. Error function and related functions . . . 126 A.4. Dipoles parallel to strong confinement . . . 127
A.4.1. Taylor expansion of the integrand for polarization parallel to the strong confinement . . . 127 A.4.2. Elliptic integrals for polarization axis parallel to the strong
confinement . . . 128 A.5. Gaussian integrals over subspace . . . 130 A.6. Excerpt from the source code . . . 131
Bibliography 135
Zusammenfassung in deutscher Sprache 147
Lebenslauf 155
Abstract
In this work several aspects of dipolar Bose-Einstein condensates are investigated. These condensates have attracted much interest since their experimental realiza-tion and offer the chance to explore the novel, fascinating, and sometimes unex-pected physics originating from the dipole-dipole interaction. For the theoreti-cal description the time-dependent variational principle with an ansatz of coupled Gaussian wave packets is applied to the analysis of this form of matter. In order to describe physical systems, where the shape of the condensate wave function cannot be approximated adequately by origin-centered Gaussian functions, translational and rotational degrees of freedom of the particular wave packets are included. This provides not only a full-fledged alternative method to full-numerical grid cal-culations that are demanding in computation time, used for the determination of ground states and investigation of the real-time dynamics, but also allows for the calculation of excited states and the examination of stability properties and offers the possibility for a deeper physical understanding of dipolar condensates.
The variational method is applied to three different physical systems. In the first one anisotropic quasi-two-dimensional solitons are investigated. Ground states and small excitations are presented and the full flexibility of the approach is utilized in the simulation of the collision of two solitons. Thereby simulations of collisions are performed, where initially the solitons are in the repelling side-by-side con-figuration and move towards each other with a specific momentum. It is shown that the results are in excellent agreement with full-numerical grid calculations. Collisions with various initial parameters such as the initial momentum and im-pact factor are performed, allowing to classify the qualitative interactions of two solitons and illuminating the effects of the solitons’ two-dimensional nature on the collision process.
In the second system dipolar Bose-Einstein condensates in a triple-well potential are investigated. This constitutes a well-suited model system for periodic optical potentials with important contributions of the non-local and anisotropic dipole-dipole interaction, which show a variety of effects such as self-organization and formation of patterns. Here, a macroscopic sample of dipolar bosons in the mean-field limit is addressed. The analysis goes beyond the calculation of ground states and clarifies the role of excited and metastable states in such systems. In particular, the formation of phases is found to originate from the interplay of several states with distinct stability properties. As some of the phases are formed by metastable
states, special attention is paid to the characteristics of phase transitions in real time and the dynamical stabilization of the condensate. The dynamical simulations that support these results, constitute a way for an experimental observation of these findings.
The third physical application is an open quantum system. Here, an external potential models gain and loss of particles and features parity with simultane-ous time-reversal symmetry. For these dipolar condensates real eigenvalues of the
nonlinear Hamiltonian are found, preserving thePT symmetry of the external
po-tential, as well as states that break this symmetry. A detailed stability analysis of these, in general complex states, is performed. The comparison with a similar system featuring a non-dipolar condensate demonstrates the role of the dipolar interaction for novel effects. Furthermore, the real-time dynamics is investigated, which reveals, among other effects, long-living oscillations of the two wells’ popu-lation.
Inhaltsangabe
In dieser Arbeit werden verschiedene Aspekte dipolarer Bose-Einstein-Kondensate untersucht. Diese Kondensate lenken seit ihrer experimentellen Realisierung großes Interesse auf sich und bieten die M¨oglichkeit, die neue, faszinierende und manchmal unerwartete Physik zu erforschen, deren Ursprung die Dipol-Dipol-Wechselwirkung darstellt.
Zur theoretischen Beschreibung wird das zeitabh¨angige Variationsprinzip mit einem Ansatz aus gekoppelten Gaußschen Wellenpaketen auf die Untersuchung dieser Art der Materie angewandt. Um physikalische Systeme zu beschreiben, de-ren Wellenfunktion nicht ad¨aquat durch im Urprung zentrierte Gaußfunktionen angen¨ahert werden kann, werden Translations- und Rotationsfreiheitsgrade der einzelnen Wellenpakete ber¨ucksichtigt. Dies erm¨oglicht nicht nur eine vollwerti-ge Alternative zu rechenaufw¨andivollwerti-gen voll-numerischen Gitterrechnunvollwerti-gen, die zur Bestimmung von Grundzust¨anden und zur Untersuchung der Realzeitdynamik ver-wendet werden, sondern erlaubt auch die Berechnung angeregter Zust¨ande, sowie die Betrachtung der Stabilit¨atseigenschaften und er¨offnet damit die M¨oglichkeit f¨ur ein tieferes physikalisches Verst¨andnis dipolarer Kondensate.
Das Variationsverfahren wird auf drei unterschiedliche physikalische Systeme an-gewandt. Im ersten werden quasi-zweidimensionale Solitonen untersucht. Grund-zust¨ande, sowie kleine Anregungen werden dargestellt, wobei die ganze Flexibilit¨at des Verfahrens in der Simulation von zwei kollidierenden Solitonen ausgenutzt wird. Dabei werden Simulationen von Kollisionen durchgef¨uhrt, bei welchen sich die So-litonen zu Anfang nebeneinander in abstoßender Position befinden und sich mit ei-nem bestimmten Anfangsimpuls aufeinander zubewegen. Es wird gezeigt, dass sich
die Ergebnisse in exzellenter ¨Ubereinstimmung mit voll-numerischen
Gitterrech-nungen befinden. Kollisionen werden mit verschiedenen Anfangsparametern, wie beispielsweise Anfangsimpuls und Streuparameter durchgef¨uhrt, was eine Klassifi-zierung der qualitativen Wechselwirkung zweier Solitonen zul¨asst und die Effekte auf den Kollisionsprozess, die sich aus der zweidimensionalen Natur der Solitonen ergeben, beleuchtet.
Im zweiten System werden dipolare Bose-Einstein-Kondensate in einem Dreimul-denpotential untersucht. Dieses stellt ein geeignetes Modellsystem f¨ur periodische optische Potentiale mit wesentlichen Beitr¨agen der nichtlokalen und anisotropen Dipol-Dipol-Wechselwirkung, die eine Vielzahl von Effekten wie Selbstorganisation und Strukturbildung zeigt, dar. Hier wird sich mit einem makroskopischen System
dipolarer Bosonen im Mean-Field-Grenzfall befasst. Die Untersuchung geht ¨uber die Berechnung von Grundzust¨anden hinaus und kl¨art die Rolle angeregter und metastabiler Zust¨ande in solchen Systemen. Insbesondere zeigt sich, dass die Bil-dung von Phasen aus dem Wechselspiel verschiedener Zust¨ande mit unterschied-lichen Stabilit¨atseigenschaften entspringen. Da manche Phasen von metastabilen Zust¨anden gebildet werden, wird der Charakteristik der Phasen¨uberg¨ange in Real-zeit und der dynamischen Stabilisierung des Kondensates besondere Aufmerksam-keit geschenkt. Dynamische Simulationen, welche diese Ergebnisse belegen, weisen eine M¨oglichkeit f¨ur experimentelle Beobachtungen dieser Resultate auf.
Die dritte physikalische Anwendung stellt ein offenes Quantensystem dar. Hier modelliert ein externes Potential Zunahme und Verlust an Teilchen und weist Parit¨ats- mit gleichzeitiger Zeitumkehrsymmetrie auf. F¨ur solche dipolaren Kon-densate k¨onnen sowohl reelle Eigenwerte des nichtlinearen Hamiltonian gefunden
werden, welche die PT Symmetrie des externen Potentials erhalten, als auch
Zu-st¨ande, die diese brechen. Eine detaillierte Stabilit¨atsanalyse der, im Allgemeinen komplexen, Zust¨ande wird vorgenommen. Der Vergleich mit dem gleichen System mit einem nicht-dipolaren Kondensat zeigt die Rolle der dipolaren Wechselwirkung bei neuen Effekten auf. Des Weiteren wird die Realzeitdynamik untersucht, welche unter anderem langlebige Oszillationen der Besetzungen der Mulden offenbart.
List of abbreviations
Abbreviations in the text
BEC Bose-Einstein condensate
BdGE Bogoliubov-de Gennes equations
DDI dipole-dipole interaction
EOM equations of motion
FFT fast Fourier transform
GPE Gross-Pitaevskii equation
GWP Gaussian wave packet
ITE imaginary-time evolution
LSE linear system of equations
MQST macroscopic quantum self trapping
NLSE nonlinear Schr¨odinger equation
ODLRO off-diagonal long-range order
PCA principle component analysis
TW triple well
Physical constants
Symbol denotation value [1]
aB Bohr radius 5.2917721092(17)× 10−11m
~ reduced Planck constant 1.054571726(47)× 10−34 Js
c speed of light 299 792 458 ms−1
µ0 vacuum permeability 4π× 10−7 TmA−1
0 dielectric constant 1/µ0c2
List of used Symbols
Some of the symbols are reused contextually which is noted at the particular point.
Symbol denotation value / comment
¯
Ψ condensate wave function
Ψ scaled condensate wave function Ψ = N1/2Ψ use in calc. of¯
integrals and TDVP as well
m atomic mass
Emf mean-field energy
µ chemical potential
a (scaled) scattering length
Γ amplitude of PT-sym. potential
N particle number or contextual usage
number of GWPs
z variational parameters contains A, p, q, γ
g Gaussian wave packet
F, F−1 Fourier and
inverse Fourier transform
Φ(z) complex error function
erfc(z) complementary error function
w(z) Faddeeva function
RD(x, y, z) elliptic integral of
1. Introduction
1.1. Preface
If bosonic particles are cooled down below a critical temperature Tc (usually some
hundred nK), a quantum phase transition takes place. This quantum phenomenon has been predicted in 1924 by Satyendranath Bose and Albert Einstein [2–4] and is known as Bose-Einstein condensation. The phase transition does not occur due to interactions of the particles but is a result of the statistics only. It took more than 70 years to reach the low temperatures which are necessary for the experimental observation. Finally, in 1995 experimental realizations of Bose-Einstein conden-sates (BECs) consisting of 87Rb [5,6] and sodium [7] atoms could be reported. In
2001 Cornell, Wiemann and Ketterle received the Nobel Prize in Physics for this remarkable achievement which has stimulated numerous further experiments and theoretical investigations in the field of cold gases.
The effect of Bose-Einstein condensation is closely related to the quantum effects of superfluidity and superconductivity and the theoretical description is based on similar ideas. Those effects have been experimentally observed much earlier (Liq-uefaction of4He in 1908, discovery of superconductivity in 1911, and description of
superfluidity in 1938) and they are key discoveries in the field of low-temperature physics. However, the development of an appropriate microscopic theory for con-ventional superconductors took a rather long time, whereas a microscopic descrip-tion of superfluidity with its fascinating effects still constitutes a challenging part of ongoing research. Few years after the first experiments with BECs, atomic
de-generate Fermi gases could be de-generated [8–11] which opened the way to Fermi
superfluidity and to a theoretical description in the context of the
Bardeen-Cooper-Schrieffer (BCS) theory [12] and the BEC-BCS crossover. The fundamentals of
these topics are nowadays explained in numerous text books, e.g. [13,14].
Bose-Einstein condensation is a macroscopic quantum phenomenon enabling one to observe quantum mechanical effects not only on an atomic or subatomic but macroscopic scale. This provides for a high degree of controllability in experiments with BECs. Most of the parameters in a BEC, such as temperature, particle num-ber, and strength of interactions are accessible to a selective variation. The latter
is tunable by Feshbach resonances [15, 16] even to negative values resulting in
possess-ing a large magnetic moment generates BECs with strong dipolar interaction. In this field the condensation of 52Cr [17], 164Dy [18], and 168Er [19] have been
re-ported. Furthermore, there has been fast progress towards the condensation of
polar molecules with electric dipole moments [20], where the dipole-dipole
inter-action is even more dominant.
The dipole-dipole interaction (DDI) behaves completely different to the short-range contact interaction. Fascinating effects appear even in classical systems,
like the Rosensweig instability of a ferrofluid [21]. The novel, interesting, and
sometimes unexpected physics originate from the unique nature of the DDI: It is long-ranged, non-local, nonlinear and anisotropic. Dipolar BECs show a variety of new effects such as an elongation of the atomic cloud, structured ground states, the roton-maxon spectrum, different collapse dynamics, solitons, pattern formation, new quantum phases, and more. A review on the recent developments in the field of dipolar BECs has been given by Lahaye et al. [22].
There are several theoretical models for the description of Bose-Einstein con-densates and in particular, of dipolar BECs. In general, these concon-densates are described by the many-particle Schr¨odinger equation (SE). In principle, this equa-tion can be solved exactly, e.g. within statistical error bounds by quantum Monte
Carlo (QMC) methods [23, 24]. However, in many cases the application of such
numerically demanding methods is not necessary or even impractical. In optical lattices the Bose-Hubbard model can be applied. Alternatively, it is often sufficient to solve the system in a mean-field approximation that leads to the Gross-Pitaevskii equation (GPE). This is a nonlinear Schr¨odinger equation (NLSE). The numer-ically exact solution of the GPE by means of representing the wave function on a grid and carrying out the time evolution e.g. by the split operator method is a challenging task and an alternative method is desired. The long-range nature of the DDI and the request for a flexible calculation of excited states enhance this demand.
In this work we investigate the physics of dipolar condensates in three scenarios. All of them have in common the DDI, but differ in the properties of the external potential. In the first system that we address, the particles are confined by an external trap only in one direction, whereas the potential is open in the other two. Thereby, it is a unique feature of dipolar BECs that the macroscopic wave function is not dissolving in the two open directions for an adequate choice of the scattering length and initial conditions. This wave function with constant shape, called a quasi-2d soliton, has been theoretically predicted [25], yet the experimental proof of such bright matter-waves has still to be accomplished. After a simple
description in [25] of the soliton ground state with a single Gaussian function
placed at the origin (supported by numerical grid calculations) more advanced
1.1. Preface
of Gaussian wave packets (GWPs), centered at the origin, were used to represent the wave function. This method has successfully been applied to dipolar BECs [26–30] and proven to be a full-fledged alternative to numerical grid calculations. A highly interesting field to explore is opened by interactions of several solitons. However, in the works given above relying on GWPs, the GWPs are centered at the origin and a variational description of soliton interactions is impossible with this ansatz. Numerical grid calculations [31–33] have demonstrated the prospects of a detailed analysis in situations such as the collision of solitons. For this reasons an essential part of this work is devoted to the development of an alternative method on the basis of the time-dependent variational principle (TDVP) and Gaussian wave packet propagation including additional translational and rotational degrees of freedom. With this method we will be able to simulate collisions of solitons with various initial conditions. The results have been published in references [26, 34].
In its original form the variational approach with Gaussian functions has been introduced by Heller [35–37] for the computation of atomic and molecular spectra and has successfully been applied to the nonlinear GPE as well [27,28,38,39]. We employ a linear superposition of GWPs as an ansatz for the TDVP and, on that basis, develop a numerical tool to evolve the wave function in real and imaginary time, far beyond the linear regime and the stationary states. As mentioned above we go beyond the ans¨atze of [26–30, 39] and include rotational and translational degrees of freedom.
In the second system we turn our attention to the interplay of long- and short-ranged interactions in dipolar BECs which is the origin of multiple interesting effects, such as self-organization and pattern formation. In the literature these ef-fects are often investigated in optical lattices [22,40, 41]. For non-dipolar conden-sates one of the most simple model systems is the double-well potential [42, 43]. Although dipolar BECs in double-well potentials have been subject of studies [44,45], the effects of the DDI’s long-range nature can be investigated thoroughly only in a system, where we are able to distinguish between on-site effects of the short-ranged interaction and the site effects caused by the long-ranged inter-action. The triple-well (TW) potential constitutes such a system. In the inves-tigation of a mesoscopic sample of dipolar bosons in a TW potential in reference
[46] the occurrence of several phases has been shown. With increasing particle
number one observes distinct phase transitions for certain paths in the parameter space of short- and long-ranged interaction. These calculations on the basis of the Bose-Hubbard model have been extended to a macroscopic description within
the GPE by Peter et al. [47], demonstrating that some of the phases are unstable
with respect to the collapse. Yet, what is missing is a fundamental understanding
of the mechanisms behind the phase transitions. As pointed out in reference [40]
necessary to consider not only the ground but also excited states. The method we apply to the collision of quasi-2d solitons provides for an easy access of the excited states and is hence suitable for the study of the phase transitions. We therefore use this method to calculate ground and excited states in the triple-well system and support our results with dynamical simulations. The results for the triple-well system have been published in reference [48].
The third system of our studies is an open quantum system, a dipolar BEC
in a complex double-well potential providingPT symmetry (parity reflection and
time-reversal) of the Hamiltonian. Bender and Boettcher have shown [49,50] for a
linearPT-symmetric Hamiltonian that for a certain form of the complex potential
real eigenvalues exist even for non-vanishing imaginary part of the potential. The
research domain of PT symmetry was driven forward by the experimental
obser-vation ofPT symmetry in optical systems [51–53]. Graefe et al. [54–56] transferred
the concept of PT symmetry to BECs, where the imaginary part of the potential
has the meaning of coherent in- and out-coupling of particles. These works are essentially based on a description involving a two-mode approximation and have
been carried on to a more realistic model in the framework of the GPE and a
PT-symmetric double-delta [57] and double-well potential [58]. To transfer these works to dipolar BECs we make use of the techniques used in the analysis of the triple-well potential and extend the method for the description of complex ground states. This enables us to calculate the stationary states, investigate their stability prop-erties, examine the real-time evolution and offers the possibility to compare with
the non-dipolar system. Results regarding the PT-symmetric system have been
published for non-dipolar BECs in references [58–60]. Results for dipolar BECs in
the PT-symmetric double-well potential are prepared for publication [61].
1.2. Outline
This work is organized as follows. In chapter 2 we will briefly introduce
Bose-Einstein condensation. We concentrate on the mean-field description of this phe-nomenon as all our calculations are based on this approximation. In this context
we will be concerned with the Gross-Pitaevskii equation in section2.3 and discuss
the two important inter-particle interactions contained in this equation and their consequences. Special unit systems for several physical systems will be discussed in the specific sections, where the results are presented (see below). Furthermore, we give a short review on the commonly-used techniques for the solution of the GPE, which we will use to compare our results with.
In chapter 3 we introduce the time-dependent variational principle and the
method used in this thesis that relies on it. After a short discussion of the basic
1.2. Outline
apply an ansatz of coupled Gaussian wave packets in section 3.1.1. The most
im-portant results in this section are the equations of motion (EOM) for the variational parameters. The computation of these requires the construction and solution of a linear system of equations. The construction, involving the computation of several
integrals, is the topic of sections 3.2 and 3.3. Thereby it is a major advantage
of our ansatz that most of these integrals can be calculated analytically. These
analytical integrals are discussed in section 3.2. The only remaining integral, the
dipolar integral, is dealt with in section 3.3. We address the numerical
implemen-tation including the compuimplemen-tation of the dipolar integral, special cases, where the evaluation can be simplified, and optimizations in the corresponding subsections of this section. With the EOM available we are able to analyze dipolar BECs by the computation of ground and excited states via either a nonlinear root search or an imaginary-time evolution, by the investigation of the stability with techniques of nonlinear dynamics, and by the real-time evolution of the EOM. These
meth-ods are presented in section 3.4. In section 3.5 the implementation is discussed,
involving a presentation of the graphical user interface.
The physics of dipolar condensates in three different scenarios is investigated in
chapters 4, 5, and 6. We discuss quasi-2d solitons in chapter 4. An appropriate
unit system is introduced in section4.1 and symmetries of the GPE that are
rele-vant for the investigations in the following sections are discussed. In section4.3the stationary states are treated and a short review of previous results is given. Fur-thermore, in section 4.4the stability of states is investigated with respect to small deviations that preserve the origin-centered, non-rotated form of the wave function and such that require the newly introduced translational and rotational degrees
of freedom. In section 4.5 we present the collision of two solitons. Additionally,
we test various sets of initial parameters, compare the results to numerically exact calculations and discuss the dependency of the solution on the initial momentum.
In chapter 5 dipolar BECs in a triple-well potential are investigated. First, the
configuration of the external potential along with an appropriate unit system is
introduced. Then, in section5.2, the methods which provide the stationary states
and their stability are introduced. Detailed results are presented in two separate
sections 5.3.1and 5.3.2, each handling one distinct cut in the parameter space.
We study a dipolar BEC in aPT-symmetric double-well potential in chapter 6.
After a short introduction toPT-symmetric systems and a review on a non-dipolar
BEC in such a system in section6.2we introduce the configurations (repulsive and
attractive) of the external and dipolar potential and present the used methods in
section 6.3. Thereafter we address the repulsive configuration in section6.3.1 and
the attractive configuration in section 6.3.2 and compare with the non-dipolar
The results of this thesis are summarized in chapter 7, where further prospects are shown as well. For a summary in German see the corresponding section on
2. Bose-Einstein condensation
In this chapter we give an overview on the different descriptions of the phenomenon Bose-Einstein condensation. We start from thermodynamic arguments regarding a bosonic system without interactions to show the evidence of the phase tran-sition. Then, we briefly present the many-particle description of BECs and the
approximations that lead to the GPE, discussed in section 2.3.
2.1. The ideal Bose-gas
One of the key ideas in the description of cold atomic gases is that the particles are considered as indistinguishable identical particles. This idea, originally used by Bose for an alternative derivation of the black-body radiation formula, can be applied to bosons featuring integer spin and fermions featuring half-integer spin. Both kinds of particles obey different quantum statistics, i.e. bosons obey Bose-Einstein statistics and fermions Fermi-Dirac statistics. In the former one there is no restriction on the number of occupied single-particle states, whereas for the latter statistics only one particle can occupy a single-particle state. Here we will focus on bosonic particles, although fermionic BECs, where two fermions together show bosonic statistical properties (by the formation of molecules or a binding
mechanism similar to Cooper-pairs) have been successfully generated as well [63].
For details of the following see references [13,14].
By the use of the corresponding combinatorics it can be shown that the average number of bosonic particles of energy kis given by the Bose-Einstein distribution
¯ nk=
1
eβ(k−µ)− 1, (2.1)
where β = (kBT)−1and µ is the chemical potential which can be interpreted via the
first law of thermodynamics. The sum Pkn¯k yields the total number of particles
N. In the thermodynamic limit
this sum becomes an integral and we obtain the particle density n (T, µ) as a function of temperature and the chemical potential
n = N V = 1 (2π)3 Z d3k 1 eβ(k−µ)− 1. (2.3)
This is an implicit equation for the chemical potential so we can conclude that if the temperature is lowered and n is kept constant then µ increases and eventually
becomes zero at a certain temperature Tc. With the dispersion relation of an
ideal Bose-gas and µ = 0 one obtains after a short derivation [14] the transition
temperature Tc for Bose-Einstein condensation of an ideal Bose-gas
kBTc = 2π~ 2 m n ζ(3/2) 2/3 ≈ 2π~ 2 m n 2.612 2/3 , (2.4)
where ζ(x) is the Riemann zeta function and ζ(3/2) ≈ 2.612. At the critical
temperature Tc the heat capacity cV = ∂u/∂T with the average internal energy
per particle u has the form cV = 15 4 kB ζ(5/2) ζ(3/2) T Tc 3/2 for T ≤ Tc, (2.5a) cV ≈ 3 2kB " 1 + ζ(3/2) 27/2 Tc T 3/2# for T ≥ Tc, (2.5b)
where we have used a first-order approximation of the generalized zeta function gν(z) =P∞l=1(zl/lν) (c.f. reference [64]). Therefore, cV has a cusp (see figure 2.1)
which implies that this indeed is a thermodynamic phase transition. Furthermore, we find from equation (2.3) for T < Tc the relation of the condensate fraction of
the ideal bose gas
N0 N = 1− T Tc 3/2 , (2.6)
where N0 is the number of particles in the k = 0 state. This shows that if
T Tc almost all particles are part of the condensate what justifies the
mean-field approximation in section 2.3, where T = 0 is assumed.
2.2. Fundamental description of ultra-cold atomic
gases
A multitude of interesting BEC-properties originates from the presence of inter-actions. In this section we give a short view on the correct description of an
2.2. Fundamental description of ultra-cold atomic gases 3 2kB 0 Tc cV T cV(T ≤ Tc) cV(T ≥ Tc)
Figure 2.1.: Heat capacity of the ideal Bose gas as a function of temperature for T ≤ Tc and T ≥ Tc, corresponding to equations (2.5). The small discontinuity at Tc
originates from the approximation used in equation (2.5b).
interacting BEC and show the approximations made to handle the problem in the mean-field limit.
For a system of N bosons the symmetrized wave function is given by
ψ(x1, . . . ,xi, . . . ,xj, . . . , xN) = ψ (x1, . . . ,xj, . . . ,xi, . . . , xN) . (2.7)
This wave function obeys the many-particle Schr¨odinger equation. If the system exhibits Bose-Einstein condensation, the de Broglie wavelength of the particles
overlap and particles at different positions x and x0 become indistinguishable.
This means, that the process of removing a particle at position x0 and inserting it
at x, has a coherent quantum amplitude and phase. We can measure this property by the definition of the one-particle density matrix
ρ(x− x0) = Dψˆ†(x) ˆψ(x0)E , (2.8)
where ˆψ†(x) and ˆψ(x0) are the creation and annihilation field operators, respec-tively. If ρ (x− x0) becomes constant for large separations |x − x0| the system
possesses off-diagonal long-range order (ODLRO). This is equivalent to a large eigenvalue of ρ (x− x0), proportional to the number of particles. By the use of the
coherent state representation of the field operators ˆ ψ(x) = √1 V X k ˆ akeikx, (2.9)
where ˆakdenotes the annihilation operator of bosons with wave vector k, we obtain
ρ(x− x0) = N0 V + 2 (2π)3 Z d3ke−ik(x−x0)f(k) , (2.10)
where f (k) is a smooth function (for details see e.g. reference [13]). The second term in equation (2.10) vanishes for large |x − x0| and we can rewrite the density
matrix as D
ˆ
ψ†(x) ˆψ(x0)E→Dψˆ†(x)E Dψˆ(x0)E , for |x − x0| → ∞ . (2.11) At this point one defines the condensate wave function or order parameter
¯
Ψ (x) =Dψˆ(x)E , (2.12)
a quantity appearing in the thermodynamic limit describing the ODLRO in the condensate. This represents the mean-field approximation used in the derivation
of the Gross-Pitaevskii discussed in section 2.3. The condensate wave function is
applicable not only to the stationary states but to non-equilibrium situations, as well. The condensate density can then be expressed by
ρ(x, t) = ¯Ψ∗(x, t) ¯Ψ (x, t) . (2.13)
2.3. The Gross-Pitaevskii equation
In this section we will introduce the central equation of this thesis, the extended time-dependent GPE. Two important nonlinear terms, viz. the short-range contact interaction and the long-range dipole interaction, account for the properties which characterize this equation. In the following we introduce the description of dipolar condensates in the low-temperature limit and discuss the particular interaction terms in the GPE.
At sufficiently low temperatures all particles in a BEC can be assumed to occupy the ground state. In this case a product ansatz can be used to variationally solve the many-particle Schr¨odinger equation in the mean-field limit. Therefore, the energy functional of the SE, with the product ansatz inserted, is minimized under the constraint of a fixed norm with a Lagrange parameter µ that appears to be the chemical potential.
2.3.1. Contact interaction
We omit the dipolar interaction of the particles here and address it in section 2.3.2. As the temperature is low, the short-range interaction of the particles can be described by the s-wave scattering if interactions are sufficiently small. This is experimentally realized by generating a dilute gas.
Describing the scattering properties of the particles by the s-wave scattering yields a short-ranged scattering potential, where we use here the following defini-tion of short- and long-ranged interacdefini-tion. If both, the differential and the total
2.3. The Gross-Pitaevskii equation
cross section converge and show no singular behavior, the scattering potential is considered as short-ranged. It can be shown that this is the case if the interaction
potential decays faster than 1/r3 (with r =|x|). For example, the Lennard-Jones
potential, often used to model the interactions between neutral atoms, is short-ranged as its slowest decay contributions are proportional to 1/r6.
The contact interactions of the condensate particles yield a nonlinear potential within the GPE. This short-range scattering potential
Vc= 4π~
2N
m ac|Ψ (x, t)|
2
, (2.14)
is isotropic and its strength is given by the scattering length ac. This quantity is
adjustable by Feshbach resonances [15,65] and we will treat acas a free parameter
in our calculations. We introduced the single-particle wave function Ψ = N1/2Ψ¯
that is normalized to unity hΨ| Ψi = 1 and thereby obtain an additional factor N
in front of the interaction terms (see also for the dipolar interaction below). The properties of the second nonlinear potential, the long-range dipole-dipole interac-tion, will be discussed in section 2.3.2.
Replacing the Lagrange parameter with µ→ i~∂tyields a nonlinear Schr¨odinger
equation, known as the time-dependent Gross-Pitaevskii equation (GPE) − ~ 2 2m∆ + Vext(x) + 4π~ 2ac mN|Ψ (x, t)| 2 Ψ (x, t) = i~d dtΨ (x, t) . (2.15)
This equation, originally derived by E. P. Gross and L. P. Pitaevskii [66, 67] for superfluids and Bose gases, is one special case of nonlinear Schr¨odinger equations which appear in the description of a number of physical applications, such as nonlinear optics, photonic crystals or water waves.
2.3.2. Dipole-dipole interaction
In this work, we deal with particles that possess a large dipole moment. The interaction of two dipoles d1 and d2 that we assume to be polarized in one direction
can be described by
Udd(x, x0) = cdd
1− 3 cos2θ
|x − x0|3 , (2.16)
where θ is the angle between the polarization axis of the dipoles and the vector
x− x0 and cdd denotes the strength of the DDI which is given by
cdd = ( d2 4π0 (electric dipoles) , µ0d2 4π (magnetic dipoles) , (2.17)
Figure 2.2.: Visualization of the kernel(2.16) of the DDI potential. On the left-hand side a 3d-visualization is shown, on the right-hand side a 2d-cut is plotted, where the ordinate is parallel to the polarization axis of the dipoles, the abscissa is one of the other perpendicular directions. Blue regions mark the repulsive, red regions the attractive character of the DDI.
with the dielectric constant 0, the vacuum permeability µ0, and d = |d1| = |d2|.
Two dipoles on top of each other experience an attractive force, whereas two dipoles beside one another repel each other. The form of the DDI potential (2.16) is visualized in figure2.2. When the potential (2.16) is included in the original SE the variation described above yields the following interaction term in the GPE
Vd = cddN
Z
d3x0|Ψ (x0, t)|2 1− 3 cos
2θ
|x − x0|3 , (2.18)
where cdd is a factor depending on the actual form of the dipoles, the unit system,
and the particle number. This potential is not only nonlinear depending on |Ψ|2,
but also nonlocal as it is a functional of the whole position space. Furthermore, the dipole potential is of anisotropic nature.
The DDI is long-ranged as it decays with 1/r3. The long-range nature of the DDI
is the reason for many new effects in dipolar BECs and is often competing with the short-range scattering interaction (2.14). Moreover, some effects like pattern formation originate from the interplay of long- and short-ranged interactions.
2.3.3. Dipolar Gross-Pitaevskii equation
Including the DDI term in the GPE (2.15) yields the extended time-dependent Gross-Pitaevskii equation
2.3. The Gross-Pitaevskii equation i~d dtΨ (x, t) = " − ~ 2 2m∆ + Vext+ 4π~2 m acN|Ψ (x, t)| 2 + cddN Z d3x01− 3 cos 2θ |x − x0|3 |Ψ (x 0, t)|2 # Ψ (x, t) , (2.19)
where m is the particle mass, Vextis the external trapping potential, N the number
of particles, ac the scattering length, and cdd the dipole strength. A derivation of
equation (2.19) can be found e.g. in reference [14]. We will deal with different
types of external potentials, however, the mathematical form of equation (2.19) will not change throughout this work.
The GPE (2.19) can be brought to dimensionless form in several ways. The appropriate choice of the unit system depends on the system we want to investigate. In particular, we have to consider the most important length for the choice of the length scale. The chosen unit system will be presented and discussed in the corresponding chapters.
The mean-field energy is defined by Emf = Ψ − ~ 2 2m∆ Ψ +hΨ| Vext|Ψi + 1 2 hΨ| Vc|Ψi + hΨ| Vd|Ψi , (2.20)
whereas the chemical potential
µ=hΨ| H |Ψi = Ψ − ~ 2 2m∆ Ψ
+hΨ| Vext|Ψi + hΨ| Vc|Ψi + hΨ| Vd|Ψi ,
(2.21) is equivalent to the expectation value of the Hamiltonian of the GPE, c.f. reference [33].
2.3.4. Properties and consequences
In the derivation of the GPE (see e.g. [14]) the many-particle wave function is
approximated by a product of one-particle wave functions. Thereby all quantum correlations have been neglected as well as any finite-temperature effects. With the GPE (2.19) we have obtained an effective one-particle nonlinear Schr¨odinger equation. In the picture of the many-particle Schr¨odinger equation the condensate fraction is equal to unity (all particles are in the ground state) and the one-particle density matrix shows perfect off-diagonal long-range order (ODLRO).
It has been shown that despite these restrictions the results obtained on the basis of the GPE are in very good agreement with numerically exact solutions
of the many-particle Schr¨odinger equation [23], particularly for dilute gases. The corrections made by methods like quantum Monte-Carlo simulations (see e.g.
ref-erence [23, 24]) grow larger if interactions cannot be treated as weak anymore.
An example is superfluid helium, where strong interactions lead to non-negligible quantum fluctuations.
The nonlinearity of the GPE has a couple of important consequences. One of the most important properties of the usual Schr¨odinger equation is the fact that the Hamiltonian is a linear operator. For this reason it is possible to construct a vector space, the Hilbert space, on the eigenstates of the Hamiltonian. In the GPE the Hamiltonian is no longer a linear operator and thus, strictly speaking, no vector space can be constructed by the eigenfunctions of the Hamiltonian. Therefore, the superposition principle is not maintained for all possible states in the manifold of the eigenstates of the Hamiltonian. The function Ψ(x, t) in the GPE (2.19) cannot be interpreted as a quantum mechanical wave function anymore but represents the macroscopic wave function, or the order parameter, of the system as introduced
in section 2.2. Although Ψ is not a quantum mechanical wave function, most
properties are similar to that case. For example, the square modulus |Ψ(x, t)|2
can be interpreted as the condensate density.
One direct consequence of the nonlinearity is that the application of the time evolution operator U (t, t0) = e−iH(t−t0) on an initial state Ψi does not necessarily
lead to a state Ψf which can be expanded in the same basis of eigenfunctions of H
as Ψi.
Further examples of effects which arise from the nonlinearity are:
• Collapse of the wave function. Obviously, in equation (2.19) an infinite
amount of energy can be stored in the scattering potential when the
scat-tering length is negative and the density |Ψ|2 grows. In case of a global
collapse only one density peak is observed, whereas during a local collapse the condensate density diverges at several positions simultaneously.
• Bifurcations and exceptional points. In nonlinear systems not only the energy eigenvalues may degenerate, but the wave functions of different states may coalesce, too. This happens at the bifurcation points, which have been shown to be exceptional points [68, 69].
• Solitary waves. In the linear Schr¨odinger equation a localized wave packet will disperse, if it is not confined by an external trap. However, in nonlinear quantum mechanics the nonlinearity and dispersion may cancel out each other. In this case, the wave function has a constant shape and it may move with a finite velocity maintaining its shape (Galilean invariance of the GPE [70]).
2.4. Discussion of common solution methods
2.4. Discussion of common solution methods
In this section we will give a short overview on the commonly used methods in finding the stationary states, evolving the wave function in real-time, investigating the stability of the stationary states, and in calculating the excitation spectrum of
the GPE. We will use those methods for comparison purposes in chapters4 and5
and it is therefore helpful to know about their capabilities, limits, and restrictions, which we address in this section.
2.4.1. Stationary states and real-time evolution
One of the most common methods to solve nonlinear wave equations and particu-larly the GPE is the solution on a grid. The initial wave function Ψi is discretized
on an – in general 3d – grid. The real-time or imaginary-time evolution (one
method to compute the ground state; see section 3.4.2) of Ψi is performed by the
application of the time evolution operator U = e−iHt in real or imaginary time,
respectively. This application is carried out for small time steps by splitting U
symmetrically with the use of the Baker-Campell-Hausdorff formula [71]
U(∆t) = e−iH∆t= e−i(T +V )∆t ≈ e−i1
2T ∆te−iV ∆te−i 1
2T ∆t, (2.22)
where the potential V = Vext+ Vc+ Vd consists of the harmonic external potential,
the scattering potential and the DDI potential, respectively. One projects the action of this approximated time evolution operator on the basis of the position operator and makes use of the possibility to insert R dν|νi hν| = 1:
Ψ(x, t + ∆t) =hx| U(∆t) |Ψi
= Z
d3p0d3x0d3phx| e−ip22∆t|p0i
hp0| e−iV (x)∆t|x0i hx0| e−ip22∆t|pi hp| Ψi
=√1
2π9 Z
d3p0d3x0d3peixp0
e−ip022 ∆t
e−ip0x0e−iV (x0)∆teix0pe−ip22∆tΨ(p).˜ (2.23)
The structure of (2.23) suggests the following algorithmic procedure:
• Fourier transform of Ψ(x) in order to obtain ˜Ψ(p).
• Inverse Fourier transform to real space. • Multiply by e−iV (x)∆t.
• Fourier transform to momentum space. • Multiply by e−ip2
2 ∆t.
• Inverse Fourier transform to real space.
The scattering potential Vc and the DDI potential Vd have to be calculated for
each time step. The latter can be evaluated by means of the convolution theorem, which results in two more Fourier transforms:
Φdd(x) = 4π 3 F −1 3k2 z k2 − 1 F{|Ψ(x)|2} . (2.24)
Here k and kz denote the absolute value of the momentum and the momentum
in z-direction, respectively. Altogether, we have to perform six Fourier transforms for each time step. Note that the first and last Fourier transforms described in the algorithmic procedure of the time evolution are only necessary if one is interested in physical quantities, e.g. the mean-field energy, whose evaluation requires the wave function in coordinate space. We conclude that this numerical scheme consists of a series of Fourier and inverse Fourier transforms. These are usually performed by a fast Fourier transform (FFT).
The advantage of such methods is that only little knowledge about the wave function and its evolution in time is required. However, some information about the wave function is necessary here as well: The whole wave function and the potentials have to be mapped on the grid. This results in conditions for the size of the grid and its resolution.
In case of the dipolar GPE the product of dipolar interaction with the wave
function is decaying with 1/r3. This means that contributions from relatively
large distances have to be taken into account. For an adequate resolution of the wave function the grid has to be very fine. These conditions result in grids large in number of cells. To perform three-dimensional calculations the typical grid
size is 512× 128 × 512 for quasi-2d solitons [72]. The whole numerical scheme
involving FFTs of such large grids renders it very demanding in computation time
and hardware resources. The implementation for graphics processing units [34]
reduces the problem, however, the computation time is still rather large and special hardware is necessary.
Another drawback of grid calculations are the difficulties in finding excited states. The imaginary-time evolution in its basic form does not yield excited states. In some cases, where the different states are separated in parameter space,
2.4. Discussion of common solution methods
the excited states are accessible by the ITE [40]. This is a consequence of the
nonlinearity of the GPE (remember the ITE is globally convergent in the linear Schr¨odinger equation). In such systems the convergence of the ITE can be very sensitive to the exact initial conditions (cf. [40]) and stochastic approaches, where the initial conditions are varied, once again struggle with the large amount of necessary calculation time.
2.4.2. Implementation of the grid calculations
We use grid calculations based on the numerical scheme presented above to com-pare our results presented in chapters4and5. K¨oberle [33,72] provided an imple-mentation for the shared-memory computation on multi-core processors. However,
for our purposes a grid size of 512× 128 × 512 grid points is necessary in order
to resolve all details in the dynamical evolution of the wave function. Since the three-dimensional Fourier transform is numerically very demanding, Zajec [34,73] implemented this scheme for graphics processing units (GPUs) using CUDA, en-abling a very high degree of parallelization. By the implementation for the Tesla C2070 computing architecture Zajec improved the performance of the algorithm (including the well-known FFTW library) by a factor of about 80 in comparison to the corresponding C algorithm∗ [72]. All grid calculations used for comparison
purposes in this thesis have been implemented and performed by Zajec.
2.4.3. Stability analysis and excitation spectrum
The stability of the stationary states obtained by grid calculations can be in-vestigated by considering small perturbations and by linearizing the perturbation amplitude in the GPE. Therefore, one applies the usual ansatz for the perturbation of the stationary state Ψ0 (cf. [74])
Ψ (x, t) =Ψ0(x) + λ u (x) e−iωt + v∗(x) eiω ∗t
e−iµt/~, (2.25) with the frequency ω and the amplitude λ of the perturbation. After inserting (2.25) into equation (2.19) and neglecting terms of second order in λ, one yields a system of coupled equation for u (x) and v (x), the Bogoliubov-de Gennes equations (BdGE). By solving this system of equations (usually on a grid as well) one is able to calculate the excitation spectrum of the stationary solution.
In chapter 3 we will discuss an alternative method on the basis of the
time-dependent variational principle, for which Kreibich et al. [74] have shown the cor-respondence between the excitation spectrum obtained from the eigenvalues of the BdGE and the eigenvalues of the Jacobian, built up on the EOM in the framework of the TDVP.
3. The time-dependent variational
principle
We already touched on the fact in the introduction of this thesis that an alterna-tive method to full-numerical grid calculations is desired for a number of reasons. To face the issues of numerical effort, accessibility of excited states, and flexibility, apparently the number of parameters has to be reduced. This can be achieved by parameterizing the wave function. In this chapter we will discuss our parametriza-tion in terms of Gaussian wave packets and show, how to derive and compute equations of motion for the variational parameters. Although this derivation is crucial in the development of an efficient tool that is able to describe dipolar BECs, it is a somehow technical step to take and the reader interested in the
physical results only may therefore proceed to chapter 4, directly.
GWPs are a natural choice for BECs confined by an external trap and a linear superposition has proven to be a full-fledged alternative to full-numerical grid calculations in the description of BEC ground states [26–28, 30, 39]. The works given above restrict the GWPs to the origin and do not allow any rotations of the wave function. Yet, these additional translational and rotational degrees of freedom are crucial in the description of the dynamics of dipolar condensates, in particular for condensates in external potentials which are more complicated
than a harmonic trap (multi-wells, open potentials,PT-symmetric potentials with
gain and loss term, etc.). For small perturbations of the stationary solutions the condensate wave function may be described without these degrees of freedom as long as the stationary solution exhibits the appropriate symmetries. If dynamics far beyond this linear regime are investigated, an advanced variational ansatz has to be applied. We therefore use the ansatz
Ψ (x, t) = N X k=1 e− (x−qk(t))TAk(t)(x−qk(t))−i(pk(t))T(x−qk(t))+γk(t) ≡ N X k=1 gk, (3.1)
where N is the number of Gaussians∗, Ak are 3× 3 complex symmetric matrices,
pk and qk are three-dimensional real vectors and γk are complex numbers. The
symbol T denotes the transposition. We omit the time argument of the variational parameters from here on and gather them in the vector zk = Ak, pk, qk, γk (if
no upper index k is supplied, z refers to the whole set of variational parameters).
The matrices Ak represent the widths of the Gaussians where the imaginary part
describes a position-dependent phase. The vectors pk and qk are the momentum
and the center of the particular GWP, respectively. The real part of γk determines
the amplitude Ak ∼ e− Re γk
of the GWP and its imaginary part is a phase factor of the Gaussian. Rotational degrees of freedom are introduced in the off-diagonal
elements of Ak, translational degrees of freedom by the vectors pk and qk. As
introduced in equation (3.1) we will use the abbreviation gk for the kth GWP.
In this chapter we use a dimensionless, particle-number scaled form of the GPE, i.e. id dtΨ(r, t) = " −∆ + Vext+ 8πa|Ψ(r, t)|2 | {z } Vc + Z d3r01− 3 cos 2θ |r − r0|3 |Ψ(r 0, t)|2 | {z } Vd # Ψ(r, t) , (3.2)
with all quantities, units and scaling laws given in section 4.1, as we use this
equation for calculations regarding quasi-2d solitons in chapter4.
3.1. The McLachlan variational principle
The application of the time-dependent variational principle (TDVP) offers the possibility to decrease the amount of degrees of freedom and still obtain a valuable wave function. The idea is to use a parametrized wave function Ψ (z(t)), where z(t) are the time-dependent variational parameters. Of course, we have in mind the ansatz (3.1) consisting of GWPs, yet we will consider no specific parametrization here and apply the explicit ansatz in section3.1.1. The variational solution of the GPE is a subset of the manifold of the exact solutions representing the macroscopic wave function. The sub-manifold of the variational solution is spanned by the variational parameters z(t). For a detailed description of the TDVP see references [35, 37, 38].
Several formulations of the TDVP exist but they all yield the same equations of
motion (EOM) ˙z(t) for the variational parameters. An overview is given in [38].
In this work we use the descriptive formulation of McLachlan [62] and recapitulate the derivation of the equations of motion. In the original form of the McLachlan variational principle the norm of the deviation of the left- and right-hand side of the Schr¨odinger equation is minimized with respect to the trial function. To minimize the quantity
3.1. The McLachlan variational principle
I =||iφ(t) − HΨ(t)||2 != min , (3.3)
φ is varied and set φ ≡ ˙Ψ afterwards. The wave function Ψ(t) is supposed to be
given at time t. The fact that the Hamiltonian H is nonlinear in the GPE does not imply any difficulties in the further calculation, because only φ is varied.
If we evaluate the quantity I in equation (3.3) and substitute Ψ(t) = Ψ(z(t)) we yield I = ∂Ψ ∂z ˙z ∂Ψ ∂z ˙z − i HΨ ∂Ψ∂z ˙z + i ∂Ψ ∂z ˙z HΨ +hHΨ| HΨi . (3.4)
In the minimization of (3.4) the variation affects only φ and carries over to the parameters z. As the wave function Ψ in equation (3.3) is not varied, we can use δz = 0. We obtain for the variation δφ expressed in terms of the variational parameters δ|φi = δ ˙Ψ (z, ˙z)E= ∂ ˙Ψ ∂zδz + + ∂ ˙Ψ ∂ ˙zδ˙z + = ∂ ∂˙z ∂Ψ ∂z ˙z δ˙z = ∂Ψ∂zδ˙z , (3.5) and for the variation of I in equation (3.3)
δI =Dδ ˙Ψ ˙ΨE+DΨ˙ δ ˙ΨE−Diδ ˙Ψ HΨE−DHΨ iδ ˙ΨE = ∂Ψ ∂zδ˙z Ψ + iHΨ˙ + ˙ Ψ + iHΨ ∂Ψ ∂zδ˙z . (3.6)
For complex z this yields the condition ∂Ψ ∂z iΨ˙ − HΨ = 0 , (3.7)
because the variations δ ˙z can be chosen arbitrarily and δ ˙z and δ ˙z∗are independent
and must vanish separately. Equation (3.7) can be rewritten in the form ∂Ψ ∂z ∂Ψ∂z ˙z =−i ∂Ψ ∂z HΨ , (3.8)
which may be abbreviated to the short form ¯
where ¯K is an Hermitian positive definite matrix. The linear system of equations (LSE) (3.9) has to be solved at given time t to obtain the equations of motion ˙z = f (z) at this point in time. To evolve the wave function in time the LSE has to be solved numerically for every time step of an numerical integrator like e.g. the Runge-Kutta algorithm. In some very simple cases (e.g. a single GWP, see section
3.1.1for GWPs) the LSE (3.9) can be solved analytically, however, the integration
of the EOM has still to be performed numerically, except for some very simple cases of h in equation (3.9) [38].
3.1.1. Gaussian wave packets
The application of the TDVP to wave packet dynamics has been developed by E. J. Heller [35–37]. We will use the ansatz (3.1), a linear superposition of Gaussian wave packets, however, a wide range of parametrizations of the wave function can be used in the TDVP. Actually, the only requirement is that the wave function is an analytic function of the parameters with respect to complex differentiation (cf. [38]).
Insertion of the ansatz (3.1) in equation (3.7) immediately shows that the sum over all GWPs can be pulled outside the bra- and the ket-vectors. The same holds for equation (3.8) and it is therefore reasonable to apply all operators on a single GWP and sum over the Gaussians afterwards. Splitting the Hamiltonian in the kinetic and potential part via H = T + V yields for a single GWP
i ˙gk− T gk = " − 2 Tr Ak+ qk4 Ak2− i ˙Akqk+ 4ipkAkqk− pk2+ qkp˙k + pkq˙k− 2iqkAkq˙k− i ˙γk − qk4 Ak2− i ˙Akx− 4ipkAkx− ˙pkx+ 2i ˙qkAkx + x4 Ak2− i ˙Akx # gk. (3.10)
Sorting by powers of x suggests the introduction of an effective time-dependent potential
Veffk = v0k+ v1kx+ xV2kx , (3.11)
with the coefficients vk0 =−2 Tr Ak+ qkV 2qk+ 4ipkAkqk− pk 2 + qkp˙k+ pkq˙k− 2iqkAkq˙k− i ˙γk, (3.12a) vk1 =− ˙pk+ 2iAk q˙k− 2pk− 2V2qk, (3.12b) V2k = 4 Ak2− i ˙Ak, (3.12c)
3.1. The McLachlan variational principle
where Vk
2 is a 3× 3 complex matrix that possesses the same symmetry as Ak,
v1k is a complex three-dimensional vector, and vk
0 is a complex number. Equation
(3.12b) can be split into real and imaginary parts
Re v1k=− ˙pk− 2 Im Akq˙k− 2pk− 2 Re Vk 2 qk
Im v1k= 2 Re Ak q˙k− 2pk− 2 Im V2kqk. (3.13)
Solving equations (3.12) for the time-derivatives yields the equations of motion for the variational parameters
˙ Ak=−4i Ak2+ iV2k, (3.14a) ˙ pk=− Re vk1 − 2 Im Ak q˙k− 2pk− 2 Re V2kqk, (3.14b) ˙ qk= 2pk+ 1 2 Re A k−1 Im vk 1 + 2 Im V2kqk , (3.14c) ˙γk= 2i Tr Ak− iqkV2kqk+ 4pkAkqk+ i pk2− iqkp˙k− ipkq˙k− 2qkAkq˙k+ ivk0. (3.14d) Although all quantities carry the index k the dynamics of the GWPs is usually coupled. The reason for this is that in order to compute the EOM (3.14), at every time step the LSE
Kv = r (3.15)
must be solved, where the matrix K and the vector r are discussed below, and
the vector v = v1
0. . . v0N, v11. . . vN1 , V21. . . V2N
T
contains the coefficients of the effective potential (3.11) that have to be inserted into equations (3.14). The vector ron the right-hand side of equation (3.15) contains the “physical potentials” and its calculation will constitute a major task of the subsequent sections. Equation (3.15) is similar to equation (3.9), however, note that in equation (3.15) the kinetic part of the Hamiltonian is included in the matrix K and (3.15) is obtained by inserting the result of equation (3.10) into equation (3.7) which then reads
* αmβngl N X k=1 v0k+ v1kx+ xV2kx− V (x) gk + = 0 , (3.16)
with l = 1, . . . , N ; α, β = x, y, z, and m + n = 0, 1, 2. If the “physical” potentials
V(x) are pulled on the right-hand side of equation (3.16), we obtain the LSE
(3.15). For details of the derivation see reference [38].
In this work we always choose the condensates to be confined strongly in one
translations in this direction or rotations about the perpendicular axes. Without loss of generality we choose the y-direction as the direction of strong confinement.
The following form of the matrices Ak and vectors pk, qk accounts for these
demands, as Ak
xz provides deformations of the condensate wave function in
xz-direction Ak= A k xx 0 Akxz 0 Akyy 0 Akxz 0 Akzz , (3.17a) pky = qk y = 0 . (3.17b)
With equations (3.17) the LSE (3.15) can be simplified and we obtain for the coefficients of v in equation (3.15) v1k= v1,xk 0 vk1,z , and Vk 2 = Vxxk 0 Vxzk 0 Vyyk 0 Vxzk 0 Vzzk , (3.18)
which reduces the dimension of the LSE (3.15) by 3N . In this reduced form the
matrix K on the l.h.s. of equation (3.15) has dimension 7N× 7N and reads
K = (1)lk (x)lk (z)lk (x2) lk (y2)lk (z2)lk (xz)lk (x)lk (x2) lk (xz)lk (x3)lk (xy2)lk (xz2)lk (x2z)lk (z)lk (xz)lk (z2) lk (x2z)lk (y2z)lk (z3)lk (xz2)lk (x2) lk (x3)lk (zx2)lk (x4)lk (x2y2)lk (x2z2)lk (x3z)lk (y2) lk (xy2)lk (y2z)lk (x2y2)lk (y4)lk (y2z2)lk (xy2z)lk (z2) lk (xz2)lk (z3)lk (x2z2)lk (z2y2)lk (z4)lk (xz3)lk (xz)lk (x2z) lk (xz2)lk (x3z)lk (xy2z)lk (xz3)lk (x2z2)lk , (3.19)
where all entries (α)lk denote N × N matrices according to
(α)lk = gl=1 α gk=1 · · · gl=1 α gk=N ... ... gl=N α gk=1 · · · gl=N α gk=N . (3.20)
With equation (3.18) the reduced solution vector v reads
v = v0k v1,xk vk1,z V2,xxk V2,yyk V2,zzk V2,xzk , (3.21)
3.1. The McLachlan variational principle
where each entry (for k = 1, . . . , N ) is a vector of length N . The vector r on the r.h.s. of equation (3.15) contains, as said before, the external and interaction potentials in the following (reduced) form
r = N X k=1 gl V (x) gk gl xV (x) gk gl zV (x) gk gl x2V(x) gk gl y2V(x) gk gl z2V(x) gk gl xzV (x) gk , (3.22)
where l = 1, . . . , N and V (x) = Vext + Vc+ Vd from the GPE (3.2). The form of
the external potential Vext depends on the actual system, whereas the contact and
dipole interaction potential are given in equations (2.14) and (2.16), respectively. Transformation of matrix A to CB-matrices
In the numerical integration of the EOM (3.14) the terms ∝ Ak2 can lead to
numerical difficulties [36]. It is therefore reasonable to introduce two auxiliary
matrices Ck and Bk. With
Ak ≡ Bk(Ck)−1, (3.23)
the l.h.s. of the equations of motion (3.14a) for the width matrices can be written as
˙
Ak= ˙Bk Ck−1− Bk Ck−2C˙k, (3.24)
where Ck and Bk are 3× 3 complex matrices. There is some freedom of choice of
the matrices Ckand Bkto fulfill equation (3.23) and furthermore, the CB -matrices
do not possess the same symmetry properties as the matrices Ak. Therefore the
transformation increases the number of (real) parameters from 2(4+1)+2+2 = 14 per GWP to 2(5 + 5 + 1) + 2 + 2 = 26 per GWP. Omitting the index k we obtain from the equations (3.24) and (3.14a)
B−1BC˙ −1− C−2C˙ =−4iB−1A2+ iB−1V2 (3.25)
⇒ B−1B˙ − C−2CC˙ =−4iC−1B+ iB−1V2C . (3.26)
By comparison of the l.h.s. and the r.h.s. we yield the equations of motion for Ck
˙
Bk= iV2kCk, (3.27a)
˙
Ck= 4iBk. (3.27b)
Equations (3.27) can then be integrated instead of equation (3.14a) with the initial conditions
Ck(0) = 1 and Bk(0) = Ak(0) . (3.28)
This splitting of Akcopes with most of the instabilities in the numerical integration.
However, note that the matrices Bk(t) and Ck(t) do not become stationary at a
fixed point Ak(t) = const. of the EOM (3.14) [30].
Conservation laws
Although the TDVP yields an approximate solution of the GPE, it can be shown that the norm and the mean-field energy (2.20) are conserved exactly, as long as they are conserved in the GPE. The method in this work aims for the calculation of stationary states and dynamical simulations of the GPE. The conservation of
norm and energy in the TDVP have been shown in [38, 75] for the linear SE, and
we will shortly recapitulate the findings here.
It has been shown in reference [75] (for a linear Hamiltonian Hl) that the
expecta-tion valuehΨ| A |Ψi of any operator A, where Ψ is the variational wave function, is conserved within the McLachlan variational principle if AΨ lies within the tangent
space, spanned by the derivative of|Ψi with respect to the variational parameters
and [Hl, A] = 0. Thus, we may write [75]
d
dthΨ| A |Ψi = 2 Re hAΨ| ˙Ψi = 2 Re hAΨ| − iHl|Ψi = hΨ| − i [Hl, A]|Ai = 0 . (3.29) As the identity operator commutes with the Hamiltonian, norm conservation is obviously fulfilled.
This still holds for the GPE and, furthermore the mean-field energy is conserved, i.e. d dtEmf = d dt Ψ − ∆ + Vext+ 1 2(Vc+ Vd) Ψ = 0 . (3.30)
Note that this only applies to the mean-field energy, which is the expectation value of the many-body Hamiltonian, and is in general not valid for the chemical potential, as the GPE-Hamiltonian depends on the square modulus of the
3.2. Analytical integrals for the equations of motion
3.2. Analytical integrals for the equations of motion
The setup and solution of the LSE (3.15) is the central step for the computation of the EOM (3.14) and (3.27). We now turn to the construction of the LSE (3.15) and therefore derive the integrals which appear in that. One of the major advantages in using Gaussian wave packets as an ansatz is that almost all integrals can be calculated analytically. The only exception is the dipolar integral which we discuss in section3.3. In fact, all integrals, except the dipolar integral, can be expressed in terms of simple Gaussian integrals. For the calculation of the Gaussian integrals and their derivatives (see below) it is convenient to rewrite the ansatz (3.1) to the equivalent form Ψ = N X k=1 e− xTA˜kx+(p˜k)Tx+˜γk , (3.31)with the complex symmetric 3 × 3 matrices ˜Ak, the three-dimensional complex
vector ˜pk, and the complex number ˜γk. The relation between the variational
parameters in equations (3.1) and (3.31) is ˜ Ak = Ak, (3.32a) Im ˜pk =−2 Im Akqk− pk, (3.32b) Re ˜pk =−2 Re Akqk, (3.32c) ˜ γk = qkT Akqk+ ipk+ γk, (3.32d) and qk=−1 2 Re ˜Ak−1Re ˜pk, (3.33a) pk=− Im ˜pk− 2 Im ˜Akqk. (3.33b)
The set of parameters can be transformed easily between the two representations. For the calculation of the integrals we use the form given in equation (3.31) and omit the tilde from here on.
The fundamental three-dimensional Gaussian integral, which is the starting point for the following derivations reads with (3.31)
I0(A, p, γ) ≡ Z d3xe−(xTAx+pTx+γ) = r π3 det Ae 1 4p TA−1p −γ, (3.34)
