also observe increases in the peak phase space density of several orders of mag-nitude, when an attractive dimple potential was ramped up adiabatically. It is notable that without detrimental effects, the position of the dimple potential could be chosen slightly off-centre in the quadrupole potential. This was done in order to yield a larger reduction in the Majorana loss rate at the zero mag-netic field centre region. This technique may thus become a practical method to achieve BEC phase space densities with particle species, which have sub-optimal collision properties. A Gaussian trap dimple is shown at a far off-centre position in figure 4.14.
−6 −8
−2 −4 2 0
6 4 8
−10
−5 0
5 10
−4
−2 0 2 4 6 8 10 12
Figure 4.14: Linear quadrupole trap potential in two spatial dimensions with a “dimple potential” depression (harmonic approximation) at an offset (x,y) position from the central potential minimum.
4.6 Results
Oxygen 17O2 was identified as a molecular species, which can potentially be trapped in magnetic traps due to its favourable Zeeman energy levels (shown in figure 3.1 b on page 49). We find, however, that the feasibility of evapora-tive cooling of oxygen strongly depends on the lowest temperatures and initial population sizes, which can be achieved by trap loading procedures, such as buffer gas cooling, as a starting point for the evaporative cooling process. The lower the initial temperature is, the larger the ratio between elastic and inelas-tic collisions will be. Additionally, for oxygen trapping, magneinelas-tic traps need to
be very shallow due to the adverse influence of strong magnetic fields on the inelastic collision rate, as explained in section 2.2.4 (page 34).
This practically rules out magnetic traps for trapping at high temperatures in the millikelvin range, since oxygen with its relatively small magnetic mo-ment of µ= 2µB requires fairly strong magnetic field gradients for trapping at such high temperatures. For a practical maximum magnetic field magnitude of 53 G, as pointed out in [54], the magnetic trap needs to be very shallow with a depth of merely 7 mK, and it would naturally need to be loaded at tem-peratures significantly below this figure. Presently no trap loading mechanism of this performance exists for molecular oxygen. Purely optical trapping and evaporation from optical traps may prove to be more feasible.
In our computer simulations of evaporative cooling of oxygen we have made optimistic assumptions about the starting parameters for evaporative cooling, and we have shown the feasibility of the evaporative cooling process in the presence of the adverse effects of Majorana spin flip particle loss. Furthermore we have shown how reduced dimensionality of the cooling process decreases the cooling efficiency and may make it infeasible, particularly if sufficiently low starting temperatures for the evaporative cooling runs cannot be achieved in the trap loading procedures.
We have investigated traps combining magnetic and optical potentials, and we have shown that an optical “dimple” potential can lead to significant in-creases in peak phase space density. The combination of magnetic and optical trapping methods has proven to work for problematic particle species, such as caesium, which could not be cooled to quantum degeneracy in purely magnetic traps. Thus we are confident that BEC phase space densities may in the fu-ture also be reached in oxygen experiments using such combined methods. One major obstacle for O2 cooling, however, is the large inelastic collision rate at temperatures above 10 mK.
In the future, further computational work is planned on cooling using com-bined magnetic and optical methods. Sympathetic cooling using multiple par-ticle species may turn out to be a viable route to BEC of molecular oxygen.
Our simulation program has been designed for use with ensembles of multiple particle species and can also be applied to a wide range of other dynamical particle ensemble problems in addition to magnetic trapping and cooling.
Chapter 5
Simulations of Bose-Einstein Condensates
Bose-Einstein Condensates (BEC) in dilute atomic quantum gases are coher-ent macroscopic matter-wavefunctions. Experimcoher-entally realised BEC clouds typically consist of 105-107 atoms and have sizes in the µm range. Despite the extremely small physical size of the BEC, accurate simulation of such clouds is a big challenge for computational physics. Computers have signifi-cantly evolved over the last decade and are now capable of making accurate fully three-dimensional numerical simulations of condensate clouds possible on normal workstation computers using the Gross-Pitaevskii Equation (GPE).
While 1D and even 2D simulations are possible using limited hardware, even using slow high-level programming languages like Matlab [120] or Octave [121], 3D simulations are a class of its own. They require a lot of memory and CPU time. High level languages are inefficient and cannot be used fir this purpose because of the overhead associated with them. Most importantly, such long run-ning simulations also mandate a consistent operating environment, which is able to run for many days or even months without overheating or crashing. These requirements strongly favour well-chosen computer hardware, the Linux operat-ing system, compiled programmoperat-ing languages and carefully debugged software.
The author would like to explicitly mention the C++ memory debugging utility valgrind [122], which is an invaluable tool for hunting down memory leakages and sloppy programming. Large and complex programs would not be possible without this or similar tools.
As an example for memory demands, consider a 1283 spherically symmetric discrete numerical grid of complex double precision floating point numbers to hold the complex wavefunction. At 16 Bytes for a complex double, we need 32 Megabytes for a single wavefunction. In order to do numerical work, we need at least three copies of this wavefunction in memory (see section 5.3.4 and appendix E), increasing the memory required to close to 100MB. Note that 128 gridpoints in one dimension cannot be considered to be a very good resolution and that a factor of two improvement in grid resolution increases memory de-mand by a factor of 23. Computer memory is very cheap today and thus the bottleneck for most numerical simulations is the computational power of the
105
computer CPU and the limited speed at which data can be read and written from and to main memory, rather than its actual size. One conceptual ansatz to improving the situation is parallelisation. Parallelisation involves identifying program components which can be split into several smaller chunks that can be processed by different computers or CPUs simultaneously. After separate computation, the constituent parts are then put back together again to form the solution or an intermediate solution. This approach works well for problems like weather forecasts on large scale parallel so-called “super computers”. In our GPE simulations we also use parallelisation across the system’s CPUs for Fourier transforms.
In the following sections we will outline the theory behind the numerical sim-ulations and describe the implementation of our simulation program before we discuss the problems we have investigated and applied the simulation program to. Much of the work described in this chapter has been done in collaboration with experimental BEC groups in Oxford and Konstanz.
5.1 Bose-Einstein Condensation
In this section we outline some of the basic theoretical concepts required to understand the simulation work presented in the subsequent sections, since a more detailed introduction is beyond the scope of this thesis. The interested reader will find a large body of literature on the subject. A good starting point can be found in [123].
5.1.1 Mean field theory and the GP Equation
A theoretical description of interacting ultra-cold bosonic gases starts with the many-body Hamiltonian forN bosons trapped in an external potentialVext(r, t)
Hˆ = where ˆΨ(r, t) is the Bose field operator and V(r−r0) is the particle interaction potential.
Particle interactions have a tremendous influence on the physics of a BEC.
In liquid helium, where first observations of superfluidity were made and the occurrence of BEC was suggested in 1938 [124], particle interactions are strong due to its liquid nature. Today, we know that BEC in liquid He is highly depleted, which means that there is a significant amount of excitations into states other than the lowest bosonic energy state, and that only a fraction of approximately 10% of the particles are condensed into the BEC.
In contrast to the situation for liquid helium, for dilute gases such as al-kali BEC in magnetic traps the interactions are very weak and the interaction potential is determined by simple s-wave scattering only (see section 2.1.4).
For this case we may replace the interaction potential with the simple binary collision effective interaction potential
V(r−r0) =U0δ(r−r0), (5.2)
5.1. BOSE-EINSTEIN CONDENSATION 107 where U0 = 4π~2a/m, containing the s-wave scattering length a. A simple description of the BEC based on the assumption of a mean field Ψ(r, t), which makes up the full Bose field ˆΨ(r, t) in combination with a small and negligible amount of excitations ˆδ(r, t) [125]
Ψ(r, t) = Ψ(r, t) + ˆˆ δ(r, t). (5.3) Ψ(r, t) =hΨ(r, t)ˆ i is the mean field of the Bose field operator and represents a complex wavefunction, the dynamics of which are described by the Schr¨odinger equation.
With the Heisenberg equation of motion for the Bose field operator i~∂
∂tΨ = [ ˆˆ Ψ,H],ˆ (5.4)
and replacing the Bose quantum field operator by a classical mean field wave-function ˆΨ → Ψ, we get the time-dependent so-called Gross-Pitaevskii (GP) equation Note that we use the normalisation convention R
|Ψ|2dr= 1, while some other authors may use the BEC particle numberN0in the normalisationR
|Ψ|2dr=N0. In such cases the Bose quantum field operator ˆΨ would be replaced by a clas-sical field Ψ as in ˆΨ → √
N0Ψ. The GP equation thus describes the atomic field in a classical approximation neglecting quantum field fluctuations. Due to particle interactions, a non-vanishing quantum depletion of the bosonic ground state is not avoidable. Beyond the mean field treatment in terms of the GP equation, BECs with small excitations can be described by the Bogoliubov-deGennes equations, which are obtained substituting equation (5.3) into the equation of motion (5.4) for the full Bose quantum field operator [126].
A more thorough derivation of the GP equation can be found, for example, in [15]. Therein, Castin points out the interesting view that a pure BEC at T=0 is a classical state of the atomic quantum field in the same way as a laser is a classical state of the electromagnetic quantum field.
Eigenstates of the trapping potential exhibit no dynamic spatial evolution.
In such cases the wavefunction Ψ(r, t) can be separated into a part of spatial and a second part of temporal dependence.
Ψ(r, t) = Ψ(r)e−iµt/~ (5.6)
Substituting this into the time-dependent GP equation leads to the time-inde-pendent GP equation. µis the chemical potential of the wavefunction.
µΨ(r) =
An excellent theoretical review of the dynamics of Bose-Einstein condensation can be found in [127].
The Thomas-Fermi approximation
In cases of strong nonlinearity as in large condensate clouds, the kinetic energy term in the time-independent GP equation can be neglected in comparison with the remaining energetic contributions. This leads to an approximation of the ground state cloud shape, which has the form of an inverted parabola for a harmonic external trap potential Vtrap
|Ψ|2 = 1
U0(µ−Vtrap) (5.8)
In practice, only the immediate cloud borders deviate from this approximation, as can be seen in the gradual improvement of a Gaussian initial state guess in figure 5.2.
5.1.2 Irrotational flow and vortices
The GP equation describes a superfluid, which cannot exhibit the common rotational flow known from and observed in normal fluids due to its nature as a complex wavefunction ψ = |ψ|eiS, where S = argψ is the wavefunction phase. The phase is multi-valued in the sense that the phase S as the complex argument is specified only on the range argψ∈(−π, π], where it is determined plus or minus an integer number multiple of 2π.
With a condensate density ofn(r, t) =|ψ(r, t)|2, a BEC at zero temperature T=0 fulfills the hydrodynamic continuity and force equations of superfluids.
The continuity equation applying to the condensate is
∂n
∂t +∇ ·(nv) = 0. (5.9)
The force equation, with g characterising the interaction strength of the parti-cles with mass min an external potentialVext, is
m∂v
where we have neglected the kinetic pressure term 2m~√2 n∇2√
nfrom inside the brackets. This is a reasonable approximation to make for large particle numbers and strong interactions. Under such conditions the condensate density is very smooth and homogeneous in a cloud centre region making this term negligible.
Deviations due to a breakdown of this approximation must be expected at the cloud borders.
The velocity fieldv(r, t) is related to the wavefunction phaseS by [128]
v(r, t) = ~
m∇S(r, t). (5.11)
For our choice of computational units, the velocity fieldv(r, t) of a BEC is given by v= 2∇S. With equations (5.9) and (5.10) this leads to the result that the flow is irrotational.– The curl, or the vorticity, of the velocity field vanishes:
∇ ×v= 0 (5.12)
5.2. THE GPESIM PROGRAM 109