Dynamical Behaviour of Radio Frequency Output Couplers for a Bose-Einstein Condensate Jan Max Walter Kr¨ uger March 2000
A thesis submitted for the degree of Master of Science in Physics of the University of Otago, Dunedin New Zealand
Abstract In the ultra-low temperature regime, Bose-Einstein condensates are well described by the Gross-Pitaevskii (GP) equation for the single particle wavefunction Ψ(r, t). In this thesis we use the GP equation to study the behaviour of two-component Bose-Einstein condensates in two spatial dimensions, coupled by an rf field. These coupled equations represent a model for an atom laser output coupler, and we study aspects of the coupler behaviour in two distinct scenarios. In the first, we simulate continuous-wave output in a harmonic trap under the influence of gravity. We study the effects of spatially selective coupling, and show how it depends on the usual mechanisms of bandwidth for coupling of two quantum states, including field strength, detuning and pulse length. We investigate the effects of noisy coupling fields on the coherence of the outcoupled populations. In the second scenario, we simulate the full dynamical behaviour of the Time-Orbiting Potential (TOP) trap in regimes of strong and weak coupling interaction. We examine the evolution in the two condensate states under the influence of the rotating cw coupling field, and we discover rich vortex dynamics. We present a simple model of vortex formation and give an explanation for vortex formation and temporal evolution in the TOP trap. Furthermore, we identify two distinct regimes in which we can interpret the mechanisms of vortex formation.
i
Acknowledgements First of all, I would like to thank my excellent supervisor Robert J. Ballagh. He was always very helpful and encouraging and despite all his other commitments, he always had enough time for me. Thanks also to my lovely girlfriend Katharine J. Webb, who kept me sane by dragging me out into the outback more than once. She believes that physicists are the only people in the world who say, “Oh, nooo! It’s Friday again!”. I am also very much indebted to all other members of the theory group who had to put up with my questions. Blair Blakie, Ben Caradoc-Davies, Adam Norrie, Andreas Penckwitt and Terrence Scott. They all were very helpful and enlightened me on many things. Furthermore, they shared my fascination for everything GNU and Linux, thanks to which no Micros˜1 software is needed at all. Credits also go to Malcolm Fraser, our group’s cluefull system administrator and my later flatmate. Thanks especially to Andreas for challenging my German and for being a brother in arms, struggling to apply fading German language abilities to physics. Often, we remained speechless and were forced to fall back to English. I wish to thank the experimentalists for helpful input and ideas, especially Andrew Wilson, Jos Martin and Callum McKenzie. Thanks also to the rest of the BEC gang, Reece Geursen, Nick Thomas and Nicola van Leuwen. Thanks go to my New Zealand friends and flatmates inside and outside the physics department. Sadly, many of them left Dunedin at some stage during the two years of my stay, after completing their degrees. I gratefully acknowledge financial support by the University of Otago 1999 Study Award, which paid for my tuition fees and living expenses, and I conclude the credits with a big “Thank You!” to my brother Hannes, sister Anna and my parents and relatives who supported me and had to put up with me going so far away for such a long time. However, thanks to the internet, I was never a lot further away than the next keyboard. My apologies to numerous people who I forgot to mention here.
Contents 1 Introduction 1.1 Bose-Einstein Condensation . . . . . . . . . . . . . . . . . . . . 1.2 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Background 2.1 Mean Field theory and the GP equation . . . . . . . . . . 2.1.1 Thomas-Fermi Approximation . . . . . . . . . . . . 2.1.2 Condensate Energy . . . . . . . . . . . . . . . . . . 2.1.3 Energy conservation . . . . . . . . . . . . . . . . . 2.2 The two-state Atom . . . . . . . . . . . . . . . . . . . . . 2.2.1 Pseudospin vector description of a two-state Atom . 2.2.2 The Rotating Wave Approximation . . . . . . . . . 2.2.3 Rabi-solution . . . . . . . . . . . . . . . . . . . . . 2.2.4 Coupling Pulses . . . . . . . . . . . . . . . . . . . . 3 Phase Evolution 3.1 Dressed States . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Shift of Energy Levels due to Coupling . . . . . . . . . . 3.2.1 Using detuned lasers to shape trapping potentials 3.3 Phase Evolution in a driven Two State System . . . . . . 3.3.1 Phase evolution for arbitrary initial inversions . . 4 Numerical Implementation 4.1 Dimensionless Units . . . . 4.2 Interaction Picture . . . . 4.3 Fourier Transform Method 4.4 Runge-Kutta Method . . .
. . . .
. . . . ii
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . . .
. . . .
. . . . . . . . .
. . . . .
. . . .
. . . . . . . . .
. . . . .
. . . .
1 1 3
. . . . . . . . .
5 5 10 11 12 13 13 15 16 18
. . . . .
20 20 22 25 26 29
. . . .
30 30 31 32 33
CONTENTS
iii
4.5
Numerical Limits . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.5.1 Limitations imposed by grid resolution . . . . . . . . . . 34 4.5.2 Limited simulation stepsize . . . . . . . . . . . . . . . . 36
5 RF 5.1 5.2 5.3 5.4 5.5
. . . . . . . . . . . . .
38 38 39 43 43 49 49 51 55 59 61 61 61 63
. . . . . . . . . . . . . . .
68 69 70 71 73 75 76 77 80 82 83 84 91 96 102 104
Output Coupler in 2D Output Couplers . . . . . . . . . . . . . . . . Spatial Effects in Coupling . . . . . . . . . . . Influence of gravity on trapped BEC . . . . . Localized Output Coupling . . . . . . . . . . . Coupling Bandwidth . . . . . . . . . . . . . . 5.5.1 Power Broadening . . . . . . . . . . . . 5.5.2 Broadening through finite Pulse Times 5.5.3 Broadening using Noisy Signals . . . . 5.6 Limits of the Simulation . . . . . . . . . . . . 5.7 Absorber . . . . . . . . . . . . . . . . . . . . . 5.7.1 Spatially filtering the Wavefunction . . 5.7.2 Complex Potential . . . . . . . . . . . 5.8 Validation of simulations . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
6 TOP Trap Simulations 6.1 The Quadrupole Trap . . . . . . . . . . . . . . . . 6.2 Time-Orbiting Potential (TOP) Magnetic Trap . . 6.2.1 Calculating the Time Average . . . . . . . . 6.3 Numerical Model of a TOP Trap . . . . . . . . . . 6.4 Spatial coupling . . . . . . . . . . . . . . . . . . . . 6.4.1 Coupling Bandwidth . . . . . . . . . . . . . 6.4.2 Localized evolution of population and phase 6.4.3 Force on a condensate in a TOP trap . . . . 6.5 Vortices . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Vortex dynamics in a TOP trap . . . . . . . . . . . 6.6.1 Vortices in the strong coupling regime . . . 6.6.2 Model of vortex formation . . . . . . . . . . 6.6.3 Vortices in the weak coupling regime . . . . 6.6.4 Angular momentum . . . . . . . . . . . . . 6.7 Discussion . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . . . .
CONTENTS 7 Conclusion
iv 105
A Adiabatic Elimination 108 A.1 Two-State System in Rotating Wave Approximation . . . . . . . 108 A.2 Adiabatic Elimination . . . . . . . . . . . . . . . . . . . . . . . 110 A.3 Accuracy of the approximation . . . . . . . . . . . . . . . . . . 111 B Simulation and Experiment
113
Chapter 1 Introduction 1.1
Bose-Einstein Condensation
Quantum Mechanics classifies all particles either as fermions or as bosons. Systems of fermions are described by antisymmetric wavefunctions and have half-integral spin, while systems of bosons are described by symmetric wavefunctions and have integral spin. Although the difference between a fermionic and a bosonic gas is small at high temperatures, both families of particles follow different quantum statistics with striking differences in the low temperature regime. The Pauli exclusion principle prohibits two fermions to occupy the same single-particle quantum state of a system, while bosons have no such fundamental restriction. In fact, a single-particle quantum state of a bosonic system can be occupied by a macroscopic number of particles. At very low temperatures, more and more particles in a bosonic system will transfer into the singleparticle state with the lowest energy, leading to the build-up of a macroscopic population. Macroscopic population further increases scattering into the same state, and one can refer to this condensation phenomenon as a phase transition into a fifth state of matter. In such a Bose-Einstein condensate (BEC), the deBroglie wavelength of the cold bosons reaches the mean particle separation. The bosons lose their individual identity and the whole assembly starts acting like, and can be described by, a single coherent macroscopic wavefunction. Einstein predicted BEC as a consequence of Bose statistics in 1924. He generalized the statistical methods used by Bose to describe photons to the case 1
CHAPTER 1. INTRODUCTION
2
of arbitrary indistinguishable particles. In 1938, London attributed the phenomenon of superfluidity in liquid helium 4 He to Bose-Einstein condensation. It has subsequently been established that because of the strong interactions between particles in a liquid, only 10% of liquid 4 He are condensed while 90% are depleted into quantum states of higher energies. Experimental realization of BEC in weakly interacting atomic systems requires advanced techniques in cooling dilute atomic clouds to ultra-low temperatures. Diluteness is an important requirement, because the negligible probability of three-particle collisions prevents condensate depletion by classical condensation into liquids or solids. A theoretical description of BEC was formulated in 1961, when the GrossPitaevskii (GP) equation was derived. This equation describes the single particle state for a weakly interacting BEC, and incorporates the interparticle interaction in an approximation using the condensate’s mean field. The full system state for N particles is then assumed to be the simple product of N such single particle wavefunctions, neglecting particle correlations, depletion and finite temperature effects. The GP equation could not be applied to describe BEC in liquid 4 He because of the large interaction and depletion. In the early 1990s, advances in technologies for magnetic trapping, and for laser cooling atoms, paved the road to achieving BEC experimentally, and in 1995, Anderson et al. [2] achieved BEC in a weakly-interacting dilute gas of magnetically trapped 87 Rb atoms. At the present time, a number of groups have experimentally created BEC in dilute vapours of alkali metals and in hydrogen, including one at Otago University, using 87 Rb. This thesis presents a study of solutions of the GP equation, a nonlinear partial differential equation, for a variety of scenarios related to current experimental concerns. Numerical methods are used to solve the two component GP equations, representing a condensate with two hyperfine internal states, coupled by classical radio frequency fields. This system can be used to model an output coupler [30], and this is the major focus of our work. Output coupling of BEC led to experiments trying to realize an “atom laser” in [8], the notion of which has been strictly defined in [39]. Atom lasers are promising to revolutionize atom optics and are subject of huge interest in current research. The first theoretical work in this area was done in one spatial dimension by Ballagh et. al. in [3]. Subsequently, the two-component equations have also
CHAPTER 1. INTRODUCTION
3
been solved in two spatial dimensions. In [14], Eschmann et. al. investigated the formation of Ramsey fringes in pulsed two-state coupling in a two dimensional condensate system. The two dimensional treatment allowed the effects of gravity to be modelled by Zhang in [40], Wallis et. al. [35] and by Jackson in [22], where the behaviour of out-coupled pulses, and also the flow through a potential constriction, have been studied. The two-component 2D version of the GP equation also allows us to model some detailed dynamics of the TOP trap. This implementation of magnetic trapping is an important current tool in experimental BEC, and in most cases it is assumed to provide a static harmonic magnetic trapping field. However, in fact only its time-averaged field is harmonic, and we show in chapter 6 that the dynamical behaviour of the trapping field, in conjunction with the radio frequency coupling field, can lead to some interesting dynamical behaviour of the condensate, including the formation of vortices.
1.2
Thesis Overview
Chapter 2 introduces the reader to the general background necessary to understand the simulations which were carried out and described in this work. It outlines the description of a Bose-Einstein condensate (BEC) as a quantum mechanical single particle state wavefunction obeying the GP equation, and it introduces the physics of coupled two-state systems. Couplers permit output coupling of population density from a trapped condensate into a fine coherent atom beam. This configuration represents a basic element of the so called atom lasers [21] and is thus of considerable importance. In Chapter 3, we review the mechanism of phase encoding of an atomic wavefunction by potentials and coupling interactions, and give detailed and to our knowledge novel interpretation of the mechanisms. Manipulation of condensate phase is an essential element for the manipulation and control of coherent atomic beams. In normal optics, this phase manipulation is achieved with well understood devices such as lenses, but a new method is required for atom optics. In chapter 4, we outline and discuss the numerical implementation of the physical problem, including the methods used to process and simulate the phys-
CHAPTER 1. INTRODUCTION
4
ical system in a computer. We discuss issues of accuracy and computational limitations of the simulated condensate systems. In chapter 5, the results of output coupling simulations under gravity are presented. Spatial coupling effects are used to create shaped output pulses and beams. We discuss bandwidth effects of radio frequency coupling, investigate the influence of a noisy rf coupling field on spatial coherence of an outcoupled “beam”, and we propose and investigate absorber methods that alleviate problems arising from condensate population, which reaches the edges of the computational grid. Chapter 6 presents the model used to describe and simulate BEC in a dynamical TOP trap in two spatial dimensions. TOP traps are based on the simple quadrupole trap with a magnetic field zero in its centre and a linear increase in all radial directions. By means of a rotating bias field, the point of the magnetic field zero is rotated in the xy plane, yielding a dynamic magnetic trapping field that time-averages to a harmonic field. We present results of two-state coupling simulations and show that under certain conditions vortices form. We also present a simple model explaining the formation of vortices in two-state coupling in TOP traps, which we have developed, and we identify two simple parameter regimes in which the mechanics of vortex formation can be analyzed. In chapter 7 we present the conclusions from this thesis.
This thesis is also available electronically from the author’s internet pages at http://hubble.physik.uni-konstanz.de/jkrueger/ or http://max.yi.org/jkrueger/ The author’s email address is
[email protected] (July 2000.)
Chapter 2 Background In this chapter we outline the theoretical framework the work in this thesis is based upon. We outline the mathematical description of Bose-Einstein condensates as macroscopic single particle wavefunctions, following a mean field approach. More comprehensive discussions of the subject are in [18][31]. In the second part of this chapter we introduce the concept of two-state coupling and outline the model of the two-state atom [1] which is the basis of our simulations of two-component BEC.
2.1
Mean Field theory and the GP equation
In an ultra-cold atomic cloud undergoing Bose-Einstein condensation, the majority of the atoms condense into the same single particle quantum state and lose their individuality. Since any given atom is not aware of the individual behaviour of its neighbouring atoms in the condensate, the interaction of the cloud with any single atom can be approximated by the cloud’s mean field, and the whole ensemble can be described by the same single particle wavefunction. Simple mean field descriptions have widely been used to theoretically describe and simulate BEC. Although effects like depletion of the condensate by higher order interactions between the condensate atoms, and the existence and interaction with a thermal fraction, are neglected in these descriptions, they give satisfactory results in the ultra-low temperature limit. In the following, a derivation of the Gross-Pitaevskii (GP) equation will be outlined. ˆ bosons trapped by an external field The hamiltonian of a gas consisting of N 5
CHAPTER 2. BACKGROUND
6
Vtrap can be described quantum mechanically using the boson field operator ˆ t). Assuming solely pairwise collisional interaction in the very dilute gas, Ψ(r, ˆ is the hamiltonian operator H 2 2 Z h ¯ ∇r 3 ˆ† ˆ t) ˆ + Vtrap (r, t) Ψ(r, H = d rΨ (r, t) − 2m Z Z 1 ˆ † (r, t)Ψ ˆ † (r0 , t)V (r − r0 )Ψ(r ˆ 0 , t)Ψ(r, ˆ t) + d3 rd3 r0 Ψ (2.1) 2 V (r − r0 ) represents an interatomic potential characteristic for species and internal states of the atoms. This full potential is commonly approximated heuristically by a simplified binary collision pseudo-potential V (r − r0 ) = U0 δ(r − r0 ),
(2.2)
treating binary collisions as hard-sphere collisions. The effective interaction strength U0 is related to the s-wave scattering length a by 4π¯h2 a , (2.3) U0 = m where m is the atomic mass. In this thesis we only consider the dilute nearzero temperature regime with a λdB , i.e. scattering length a much less h2 , excluding the treatment of length than the deBroglie wavelength λdB = 2π¯ mkT scales shorter than the scattering length. There is a considerable body of literature about this. A very comprehensive discussion of interaction effects can be found in [31]. Morgan showed how divergences of the above pseudopotential approximation in the finite temperature regime can be mitigated by making the approximation on higher order terms of the scattering T matrix series, which V is the first term of. There are a number of different ways of deriving the GP equation, perhaps the most conceptually consistent of which can be found in [18]. Other discussions are also in [31] [34]. In the following, we will give a somewhat naive treatment which nevertheless has the virtue of indicating in a simple way how the nonlinear term in the GP equation arises. ˆ t) satisfy the following commutation relaThe boson field operators Ψ(r, tions: h i h i 0 † † 0 ˆ ˆ ˆ ˆ Ψ(r, t), Ψ(r , t) = Ψ (r, t), Ψ (r , t) = 0 (2.4) h i ˆ t), Ψ ˆ † (r0 , t) = δ(r − r0 ) Ψ(r, (2.5)
CHAPTER 2. BACKGROUND
7
ˆ t) can be From these relations, the Heisenberg equation of motion for Ψ(r, calculated and one obtains i¯h
i h ˆ t) ∂ Ψ(r, ˆ t), H ˆ (2.6) = Ψ(r, ∂t 2 2 h ¯ ∇r † ˆ ˆ ˆ t) = − + Vtrap (r, t) + U0 Ψ (r, t)Ψ(r, t) Ψ(r, 2m
ˆ t) can in general be written as a sum over all parThe field operators Ψ(r, ticipating single-particle wave functions Ψi (r, t) and the corresponding boson creation and annihilation operators. ˆ t) = Ψ(r,
X
Ψi (r, t)ˆ ai
(2.7)
i
The boson creation and annihilation operators obey the commutation rules h i † a ˆi , a ˆj = δij h i [ˆ ai , a ˆj ] = a ˆ†i , a ˆ†j = 0 (2.8)
and their operation on the system state in Fock space (tensor product of all single state Hilbert spaces) is defined as a ˆ†i |n0 , n1 , . . . , ni , . . . i = a ˆi |n0 , n1 , . . . , ni , . . . i =
√
ni + 1 |n0 , n1 , . . . , ni + 1, . . . i √ ni |n0 , n1 , . . . , ni − 1, . . . i.
(2.9)
where the ni denote the bosonic populations of the particle states. Since the main characteristic of BEC is that most participating particles condense into the lowest single particle quantum state, it is possible to separate out the condensate part Ψ0 (r, t)ˆ a0 of the generalised mean field operator. With a total number of particles N , the population n0 of the lowest state is macroscopic such that N − n0 N . In this case (with N0 ≡ n0 ), there is no significant physical difference between states with N0 and N0 + 1 so that operators a ˆ0 and a ˆ†0 (defined above) in the nonlinear mean field part of the √ generalised version of equation (2.1) can be replaced by N0 . This is well known as the Bogoliubov approximation. Using the Bogoliubov approximaˆ t) is written as a sum of its expectation value tion, the field operator Ψ(r, √ ˆ t) representing the ˆ t)i = N0 Ψ0 (r, t) and an operator φ(r, Ψ(r, t) = hΨ(r,
CHAPTER 2. BACKGROUND
8
remaining populations in thermal states, which can be considered vanishingly small in the zero temperature BEC regime. ˆ t) ˆ t) = Ψ(r, t) + φ(r, Ψ(r,
(2.10)
ˆ †Ψ ˆΨ ˆ term The decomposition (2.10) leads to the following expression for the Ψ in (2.1) [34] ˆ †Ψ ˆΨ ˆ = |Ψ|2 Ψ + 2|Ψ|2 φˆ + 2Ψφˆ† φˆ + Ψ2 φˆ† + Ψ∗ φˆφˆ + φˆ† φˆφˆ Ψ
(2.11)
Only the first term of this result is used in the GP equation. The next four terms represent interactions between the condensate component and excited states, which are neglected. The last term represents more complex interactions between particles, which are also neglected. Thus, it should be kept in mind that the GP equation we are about to obtain is an approximation to a more complete quantum field theoretical description, with validity restricted to applications of dilute gases without significant depletion of the condensed component at temperatures of about T ≈ 0. A conceptual difficulty worth noting, and arising from this simple form of the mean field approach, is the inconsistency of a non-zero average of the field operator a ˆ 0 with the conservation of particle number. This is sometimes referred to as the issue of broken symmetry. Again, the interested reader is referred to a comprehensive resolution of these difficulties by Gardiner [18]. By substituting the decomposition (2.10), within the approximation of (2.11) as outlined above, into equation (2.1), and normalising the condensate R wavefunction Ψ(r, t) to dr|Ψ(r, t)|2 = N0 , one obtains the Gross-Pitaevskii equation 2 2 h ¯ ∇r ∂Ψ(r, t) 2 = − + Vtrap (r, t) + U0 |Ψ(r, t)| Ψ(r, t) (2.12) i¯h ∂t 2m ˆ t) have As indicated above, all terms containing the perturbation operator φ(r, been neglected in (2.12). Strictly speaking, it is therefore only valid in the zero temperature case when all atoms are in the condensate component. However, even in that case, the GPE is an approximation, since even at T = 0, a certain fraction of the atoms is still depleted from the ground state due to interactions. In 4 He, the depletion of the condensed component is about 90%, while in condensed alkali gases such as Rb and Cs, the depletion is in the order of a mere 0.5%, justifying application of the GP equation in the latter case.
CHAPTER 2. BACKGROUND
9
Coupled two-component GP equations In this work we consider coupled two-component BEC, as have previously been investigated computationally in one dimension in [3] and [14]. We are going to analyse and describe properties of such systems in two dimensions in detail at a later time in this thesis. However, in context with outlining the derivation of the single component GP equation, we want to state at this point how the extension to two components is done. The GP equations for a coupled two state system under the influence of gravity, as used in chapter 5, introduce the concept of two coupled states (see section 2.2) and include the potential of gravity in addition to the harmonic trapping potential. The coupled equations for the two components, which are denoted by the appropriate subscripts, are: r2 Ω ∂ψ1 = i∇2 ψ1 −i ψ1 +i ψ2 −iC |ψ1 |2 +w|ψ2 |2 ψ1 −iGyψ1 (2.13) ∂t 4 2 2 ∂ψ2 r Ω = i∇2 ψ2 −ik ψ1 +i ψ1 −i∆ψ2 −iC |ψ2 |2 +w|ψ1 |2 ψ2 −iGyψ2 (2.14) ∂t 4 2
The units are given in a dimensionless computational basis, which is described in detail in section 4.1. We assume a harmonic trapping potential with an untrapped, equally trapped or stronger confined second state (adjustable through k = [0, 1, 2]), and Ω is the Rabi frequency providing the coupling interaction between the two states. ∆ is the detuning in the coupling process, and w is the ratio of intra– and intercomponent scattering lengths. It is considered to be w = 1, as we will see. C represents the coefficient of the nonlinear term, and G accomodates the potential of gravity, which acts along the positive y axis in the representation we have chosen in our simulations. The time-independent GP equation In certain cases, i.e. for eigenstates of a harmonic trap, the wavefunction Ψ(r, t) can be separated into parts of spatial and time dependence µt
Ψ(r, t) = Ψ(r)e−i h¯
(2.15)
with eigenvalue µ representing the chemical potential of the system at zero temperature. Substitution of 2.15 into the time-dependent GP equation leads
CHAPTER 2. BACKGROUND
10
to the time independent GP equation 2 2 h ¯ ∇r 2 + Vtrap (r) + U0 |Ψ(r)| Ψ(r). µΨ(r) = − 2m
(2.16)
While this thesis concentrates on the time-dependent GP equation, the steady state case of equation (2.16) is needed for analysis purposes and for computation of initial (eigen-)states. A comprehensive discussion of the timeindependent GP equation and the generation of eigenstates in harmonic traps can be found in [7].
2.1.1
Thomas-Fermi Approximation
The time independent GP equation 2.16, with nonlinearity C, and for a harmonic trapping potential r2 Vtrap = , (2.17) 4 can be simplified in the so-called “Thomas-Fermi Approximation”: For large nonlinearities C, and thus for high particle numbers and high self-energy as in real experimental BEC, the kinetic energy term ∇2 Ψ becomes small compared to the high self-energy and can be neglected. We then get µ=
r2 + C|Ψ|2 , 4
(2.18)
which gives a simple analytic solution for the wavefunction. The population
|Ψ| 2
Harmonic 1D Eigenstate in TF Limit, C=1000
x=2 µ
-30
-20
-10
0 x
10
20
30
Figure 2.1: Shape of a harmonic eigenstate in 1D within the TF limit. C=1000.
CHAPTER 2. BACKGROUND
11
density |Ψ|2 for this case of a harmonic trap is plotted in figure 2.1 and goes √ to zero at x = 2 µ. Consequently the normalisation condition in 1D is √
2 Z
µ
|Ψ|2 dx ≡ 1,
√
(2.19)
−2 µ
which leads to the result for the chemical potential µ
µ=
3C 8
23
.
(2.20)
In two dimensions, the normalisation condition is √
2 Z2π Z 0
µ
|Ψ|2 r dr dφ ≡ 1,
(2.21)
0
which results in µ=
C 2π
12
,
(2.22)
and for completeness, in 3D we have normalisation condition √
2 Zπ Z2π Z 0
0
µ
|Ψ|2 r2 sin(θ) dr dθ dφ ≡ 1,
(2.23)
25
(2.24)
0
resulting in µ=
15C 64π
.
It can be seen that the C values for given chemical potentials increase with dimensionality.
2.1.2
Condensate Energy
The total energy of a condensate is the sum of three components: • kinetic energy
Ekin = −
• potential energy • interaction energy
Epot =
R
R
V
V
Ψ∗ ∇2 ΨdV Ψ∗ Vtrap ΨdV
Eself = 12 C
R
V
|Ψ|4 dV
CHAPTER 2. BACKGROUND
12
The kinetic energy arises from the Laplacian term in the Hamiltonian (2.1), and gives rise to the spatial movement of the condensate atoms. In the ThomasFermi approximation this term is zero. Potential energy originates from the confining trap and, in some of the cases we consider, also from gravity. The interaction energy, also called the selfenergy, originates from collisional forces between the atoms, and in multicomponent cases can be further subdivided into intra- and intercomponent interaction energies. In the cases we consider, the interaction energies are always positive. The condensate energy expectation value for a single component BEC is Z 1 4 ∗ 2 ∗ (2.25) E = hΨ|H|Ψi = −Ψ ∇ Ψ + Ψ V Ψ + C|Ψ| dV. 2 V It is worth remarking on the difference in the factor of 21 that accompanies the self energy term in the Hamiltonian (2.1) and in the GPE (2.12). While the chemical potential µ in the GP equation (2.16) represents the energy one adds to the system by adding one more atom, the Hamiltonian represents the average energy of an atom within the condensate. In the next section it will be shown that total energy is conserved in systems without interaction of an external coupling field. This is also given for twostate systems. Section 5.8 outlines, how energy components are calculated in numerical simulations.
2.1.3
Energy conservation
The GP equation is unitary and preserves population. It also conserves energy as we now show. The time derivative of equation (2.25) is Z dE ∂Ψ∗ 2 ∂Ψ ∂Ψ∗ ∂Ψ = − ∇ Ψ − Ψ ∗ ∇2 + V Ψ + Ψ∗ V dt ∂t ∂t ∂t ∂t V (2.26) ∗ ∂Ψ ∗ 2 ∂Ψ ∗2 +C Ψ Ψ +C ΨΨ dV. ∂t ∂t For energy conservation this expression must vanish. To simplify the above, Green’s second identity [9] can be used. Z Z 2 2 u∇ v − v∇ u dV = (u∇v − v∇u) · da (2.27) V
S
CHAPTER 2. BACKGROUND
13
When u and v are wavefunctions with zero density for r → ∞, the surface integrals on the RHS of Green’s second identity vanish, leading to the following equality: Z Z 2 v∇2 u dV. (2.28) u∇ v dV = V
V
term Using the above equation allows us to rearrange the order of the Ψ∗ ∇2 ∂Ψ ∂t in (2.26) before substituting the time dependent GP equation and its complex conjugate into it. With this, it is straightforward to show that equation (2.26) becomes dE d = hΨ|H|Ψi = 0. (2.29) dt dt Thus energy is conserved in single component simulations. Because the classical coupling fields we use to model output coupling in chapters 5 and 6 are not conserving energy, we do not expect energy conservation in our coupled two-state simulations.
2.2
The two-state Atom
In this section, a review of the physical description of coupled two-state systems will be given. The significance of the model becomes evident when we take away kinetic energy and nonlinear self-interaction terms from the coupled twocomponent GP equations, as we have done e.g. for localized analysis and for debugging purposes. In this case, the behaviour of each microscopic point in a condensate accurately follows the coupled two-state model, and thus it must be considered a fundamental principle of our two-component coupling simulations.
2.2.1
Pseudospin vector description of a two-state Atom
Bloch’s spin vector formalism [1] provides a convenient way of visualizing the effects of coupling between the two states of a spin 12 system (or a general two-state system) such as an atom in an electromagnetic field. If the system state is a linear combination of two states as follows, |Ψ(t)i = a(t)|+i + b(t)|−i,
(2.30)
CHAPTER 2. BACKGROUND
14
it proves useful to define the following quantities u = ab∗ + a∗ b v = −i(ab∗ − a∗ b)
(2.31)
w = a∗ a − b∗ b.
For a spin 21 system these can be immediately recognised as the spin components sx , sy and sz . More generally, for an arbitrary two-state system, we interpret these as the in-phase and out of phase components of the magnetic dipole, and the inversion. From population conservation in the two states, i.e. |a|2 + |b|2 = 1, it can be shown that probabilities are conserved over time, and that vector s(t) = (u, v, w) traces out curves on the Bloch unity sphere as it evolves. u2 (t) + v 2 (t) + w 2 (t) = (|a(t)|2 + |b(t)|2 )2 = 1. (2.32) The system state |Ψi (2.30) evolves as i¯h
∂ |Ψ(t)i = H|Ψ(t)i ∂t
(2.33)
with Hamiltonian H =h ¯
"
0 m·B cos(ωt) h ¯
m·B h ¯
cos(ωt) ω0
#
.
(2.34)
The off-diagonal coupling terms represent the magnetic dipole interaction of dipole moment m of the atomic state and the oscillating magnetic component B(t) = B0 cos(ωt) of the coupling field, while h ¯ ω0 is the energy difference between the two states. (A discussion purely in terms of the two states can be found in appendix A.1). With the above hamiltonian, the components of s˙ (t) become u(t) ˙ =
−ω0 v(t),
2mB(t) w(t), h ¯ − 2mB(t) v(t). h ¯
v(t) ˙ = ω0 u(t) + w(t) ˙ =
(2.35)
Feynman, Vernon and Hellwarth [15] showed that this has a geometrical explanation. Defining a vector 2m · B (2.36) Ω= − , 0, ω0 , h ¯
CHAPTER 2. BACKGROUND
15
the equations for components u, v and w of s(t) obey d s = Ω × s, (2.37) dt which is exactly the equation for a magnetic dipole (s) precessing in a torque (Ω).
2.2.2
The Rotating Wave Approximation
For all practical purposes the transition energy h ¯ ω0 is orders of magnitude larger than the magnetic dipole interaction energy m · B, even for strong coupling fields. Thus, Ω has its largest component along the z axis, with a small component along the x axis, oscillating very fast with the frequency ω0 of the interacting field B. It is convenient to decompose the oscillation along x into two circular components, so that Ω can be written as Ω = Ω+ (t) + Ω− (t) + Ω0
(2.38)
where Ω0 =
(0, 0, ω0 )
cos ωt, − 2mB sin ωt, 0) Ω+ = (− 2mB h ¯ h ¯ Ω− =
(− 2mB cos ωt, 2mB sin ωt, 0) h ¯ h ¯
(2.39) (2.40) (2.41)
Ω+ represents a vector, rotating clockwise at the driving frequency ω whereas Ω− rotates counter-clockwise. In the absence of coupling, the torque vector causes precession of s about the z axis at a rate ω0 . We can simplify the description by moving into a frame of reference which is co-rotating with the precession. Assuming conditions near resonance, i.e. ω ≈ ω0 , Ω+ rotates with (small) frequency ω−ω0 in this frame and accounts for a steady and cumulative influence on the vector of the system state. Vector Ω− is counter-rotating with ω + ω0 , which is a comparatively large frequency, and thus it does not exhibit a cumulative effect on the system state. It reverses itself far too quickly and has a negligible time-average. The Rotating Wave Approximation (RWA) neglects Ω− . After dropping the Ω− term, the torque vector Ω becomes Ω0 within the rotating wave approximation, and has the components 2m · B 0 Ω = − , 0, ω0 − ω = (−Ω0 , 0, ∆) . (2.42) h ¯
CHAPTER 2. BACKGROUND Now vector equation
becomes
d 0 s = Ω0 × s0 dt
u 0 −∆ 0 u d 0 Ω0 v v = ∆ dt w 0 −Ω0 0 w
16
(2.43)
(2.44)
and after transforming s into its rotating frame equivalent with Ω0 = 2mB h ¯ 0 0 0 0 s = (u , v , w ) using transformation
u0 cos ωt sin ωt 0 u 0 v = − sin ωt cos ωt 0 v . 0 0 1 w0 w
(2.45)
Above, we specified the Rabi-frequency Ω in terms of the interaction between magnetic moment m and the coupling magnetic field B. The Rabifrequency Ω is the rate which coherent transitions between the two levels of the system are induced at, when the system is on resonance, i.e. ∆ = 0. We will call this Ω0 from now on. For non-zero detunings ∆, the Rabi frequency as the oscillation frequency between the two states becomes detuning dependent, and we get q (2.46) Ω(∆) = Ω20 + ∆2 .
In the rotating wave approximation we made the assumption that the magnitude of the magnetic dipole interaction is small compared with the transition energy. In the case of optical transitions and electric dipole interaction, caution is indicated only for the strongest interacting laser fields [38]. For magnetic dipole interactions, transition energy and dipole interaction both are orders of magnitude smaller. Thus the rotating wave approximation can be considered to be valid for all magnetic interaction fields that can practically be generated in a laboratory.
2.2.3
Rabi-solution
Some well known solutions exist for the steady state of (2.44), for the case of constant ∆ and Ω0 . The “Rabi solution” can be obtained by two consecutive rotations of the “torque” vector Ω. The first rotation is about the v axis to
CHAPTER 2. BACKGROUND
17
align Ω with the u axis. (In the general case with detuning, ∆ is in the u − v plane.) w
∆
Ω
χ
v
Ω0 u
Figure 2.2: Torque vector Ω is constant in the rotating frame and lies in the u-w plane. The angle Ω is rotated by is tan χ = Ω∆0 as is illustrated in figure 2.2. This rotation is decribed by 0 cos θ 0 sin θ u u 0 (2.47) 0 1 0 v . v = 0 w − sin θ 0 cos θ w
After this rotation, the time evolution of the Bloch vector (2.44) simplifies to 0 u 0 0 0 u d 0 (2.48) 0 Ω(∆) v . v = 0 dt w 0 −Ω(∆) 0 w0
And vector ρ0 is now counter-rotated by an angle of −Ω(∆)t, thus moving into a coordinate frame where s00 = (u00 , v 00 , w00 ) is stationary. The relation between the two rotated vectors s0 and s00 is: 0 00 u 1 0 0 u 0 00 (2.49) v = 0 cos(Ω(∆)t) − sin(Ω(∆)t) v w0
0 sin(Ω(∆)t)
cos(Ω(∆)t)
w00
Now, knowing the initial value of the Bloch vector s = (u0 , v0 , w0 ), an analytic solution for all times can be obtained by applying all three rotations
CHAPTER 2. BACKGROUND
18
consequently in the order discussed above: Rotate s until Ω is parallel with the u axis, rotate depending on the length of the application of the “torque” and reverse the first rotation. The resulting matrix multiplication leads to an analytical expression for the temporal evolution of the inversion: ∆ · Ω0 (1 − cos(Ω(∆)t)) Ω(∆)2 Ω0 −v0 sin(Ω(∆)t) Ω(∆) ∆2 + Ω20 cos(Ω(∆)t) +w0 Ω(∆)2
w(t, ∆) = −u0
(2.50)
An interesting special case of this general solution is that of an initial condition u0 = v0 = 0 and w0 = −1 (ground state): w0 (t, ∆) = −1 +
t 2Ω20 sin2 (Ω(∆) ) 2 Ω(∆) 2
(2.51)
From an initial inversion w0 = −1, the inversion oscillates with a frequency depending on the h i effective Rabi frequency Ω(∆), with a maximum inversion Ω20 of 2 Ω2 +∆2 − 1 , which also depends on the detuning. 0 This is illustrated for four different detunings in figure 2.3. The illustration also shows the effect of gradual “dephasing” of Rabi oscillations with different detunings, with an initial simultaneous population transfer across a wide range of detunings during an initial time t ≈ π/2. In the case of zero detuning, the previous expression simplifies to w(t; 0) = − cos(Ω0 t), describing simple resonant Rabi-cycling with complete population transfer.
2.2.4
Coupling Pulses
As has been shown in the previous section, that if one starts with a system in its ground state, so that w0 = −1 and u0 = v0 = 0, complete population inversion can be achieved with a so called “π-pulse”. After a time δt such that Ωδt = π,
(2.52)
where Ω is the effective Rabi frequency. The solutions of equation (2.51) show that, for Ω = Ω0 , w = 1 after a coupling pulse like this. For non-zero Ω20 detunings, the inversion after a “π-pulse” is w = 2 Ω2 +∆ 2 − 1. For other times 0
CHAPTER 2. BACKGROUND
19
1 0.8 0.6
Inversion
0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
0
2π
4π
6π
t
Figure 2.3: Rabi solutions for Ω0 =1 showing the effects of different detunings ∆. The highest curve represents resonance. The curves with lower peaks are detuned from resonance by 0.2, 1 and 2 times the on-resonance Rabi frequency, respectively. δt the quantity Ω0 (t2 − t1 ) gives the overall angle the pseudospin vector s was rotated by due to the pulse, to calculate the corresponding inversion. In the more general case of time dependent Rabi frequencies due to time dependent fields, the relevant property of the pulse is its envelope area A(t) =
Zt
Ω(t0 )dt0 = θ(t)
(2.53)
−∞
Pulses of f π, where f is some fraction between 0 and 1, have been used in a huge variety of applications, such as atomic clocks, Ramsey fringes [14], and more recently in measuring phase coherence of BECs [20].
Chapter 3 Phase Evolution Control over phase is the key to all wave optics. In classical optics, controlling the phase of waves permits engineering beam shapes. Lenses of media with different phase velocities are used to impose specific phase signatures on beams as they pass through to make them converge or diverge. In atom optics, phase has the same importance, but also has a specific interpretation: Phase gradients correspond to velocities. By imposing specific phase profiles on matter waves, they can be made evolve in similar ways as their classical optical analogs, exhibiting effects like interference and diffraction. Because of the central importance for atom optics, in this chapter we will describe the evolution of phase in coupled two-state systems in detail.
3.1
Dressed States
In Atom Optics, coupled systems may be described in the familiar “dressed state” picture [12]. It is not only the energy levels of the atomic system that undergo a quantized treatment, but also the interacting radiation field. Thus, atoms and field are included in the system hamiltonian and diagonalisation yields the dressed states. Looking at a two state system interacting with a strong radiation field like that of a laser, the description will be in terms of states that describe the atomic states of |ui (“upper”) and |li (“lower”) at the same time as the Fock states of the radiation field. While there is an arbitrary number of general
20
CHAPTER 3. PHASE EVOLUTION
21
atom-field states for all numbers N |ui|N i,
∀N
|li|N i,
∀N,
(3.1)
the following atom-field states are chosen as a two state basis that constitute a closed system, within some well defined approximation. |Ii = |ui|N − 1i
(3.2)
|IIi = |li|N i In state |IIi the system is in the upper atomic state |ui after having absorbed one photon from the N photon radiation field. The unperturbed energies of these two atom-field states are: H0 |Ii = h ¯ (ω0 − ωL ) = −¯h∆
(3.3)
H0 |IIi = 0. Here, we introduced the unperturbed hamiltonian H0 of the atoms plus radiation, and the detuning ∆ as the difference between the transition frequency and the laser frequency coupling the two states. The energy is arbitrarily set to zero for state |IIi and H0 is diagonal in the {|Ii, |IIi} basis. When the photon energy is higher than the transition energy, ωL > ω0 (positive detuning), state |IIi is energetically higher than state |Ii. For negative detunings the situation is reversed. Due to the coupling radiation field of frequency ωL , a perturbation is introduced into the system and the states of equation (3.2) are no longer eigenstates. Under the rotating wave approximation (compare section 2.2.2), diagonalisation of the perturbed hamiltonian (3.9) connecting states |Ii and |IIi yields the eigenstates |N +i and |N −i. These so-called “dressed states” [11] are linear combinations of the atom-field states (3.2). |N +i = α|Ii + β|IIi
(3.4)
|N −i = β|Ii − α|IIi Eigenenergies of the dressed states are E+ and E− , also found by diagonalising the perturbed hamiltonian.
CHAPTER 3. PHASE EVOLUTION
3.2
22
Shift of Energy Levels due to Coupling
Without coupling, the eigenenergies of the states |Ii and |IIi are in general H0 |Ii = E1 |Ii
(3.5)
H0 |IIi = E2 |IIi. Coupling between the two states introduces a perturbation W so that the hamiltonian of the coupled system becomes H = H0 + W.
(3.6)
This perturbation shifts the energy levels of the original states and the eigenequations become H|N +i = E+ |N +i
(3.7)
H|N −i = E− |N −i. Assuming a purely non-diagonal coupling matrix (W), with Rabi frequency Ω 0 characterising the coupling strength, where ∗ = W12 = W21
h ¯ Ω0 , 2
the matrix representing H in the {|Ii, |IIi} basis is written # " E1 W12 (H) = W21 E2
(3.8)
(3.9)
and one finds the eigenenergies 1 (E1 + E2 ) + 2 1 (E1 + E2 ) − = 2
E+ = E−
1p (E1 − E2 )2 + 4|W12 |2 2 1p (E1 − E2 )2 + 4|W12 |2 . 2
(3.10)
The relations between the {|Ii, |IIi} basis and the new eigenvector basis {|N +i, |N −i} are |N +i = α|Ii + β|IIi
(3.11)
|N −i = β|Ii − α|IIi |Ii = α|N +i + β|N −i |IIi = β|N +i − α|N −i
(3.12)
CHAPTER 3. PHASE EVOLUTION
23
where α and β are: α=
r
Ω0 + ∆ , 2Ω0
β=
r
Ω0 − ∆ 2Ω0
(3.13)
Here, the generalized Rabi frequency Ω0 , which depends on the detuning h ¯∆ = E1 − E2 , is defined as q 1p 0 (¯h∆)2 + (2|W12 |)2 . (3.14) Ω = ∆2 + Ω20 = h ¯ E
E+ E1
Ωo /2 Em
’ Ω/2
∆
E2 E-
Figure 3.1: Variation of the eigenenergies E+ and E− of the coupled system with respect to the detuning ∆. In the absence of coupling interaction, the energies E1 and E2 cross at the origin. Under the effect of non-diagonal coupling, the energy levels repel each other, leading to the “avoided crossing”. Using Em = 21 (E1 + E2 ) and the above expressions, the eigenenergies become h ¯ Ω0 2 h ¯ Ω0 = Em − . 2
E+ = E m + E−
(3.15)
An illustration visualising the effect of detuning in a coupled system is found in figure 3.1. Without the coupling W , the eigenenergies cross at the
CHAPTER 3. PHASE EVOLUTION
24
origin for zero detuning ∆ = 0. The coupling interaction between the two states separates the energy levels for all detunings, leading to the “avoided crossing” in the illustration. For ∆ = 0, the shift is ±|W12 | for the respective states, separating the levels by the energy h ¯ Ω0 . The asymptotic nature of this level shift at large detuning is easily found 12 |: from a power series expansion [11] of Eq. (3.10) in terms of | W h ¯∆ ! W12 2 h ¯∆ − ... . 1 + 2 E± = E m ± (3.16) 2 h ¯∆
Figure 3.2 illustrates the situation of the level shifts for a fixed detuning ∆. The “bare” atom-field states are separated by the energy h ¯ ∆, the dressed states Unperturbed levels
Levels perturbed by coupling
, h (Ω−∆)/2
|l>|N> h∆ |u>|N-1>
, hΩ , - h (Ω−∆)/2
Figure 3.2: Shift of dressed state energy levels in the presence of a detuned coupling field. (Diagram drawn for the case ∆ > 0). are separated by h ¯ Ω0 . Since equations (3.10) describe a symmetric level shift, 0 for all values of the magnitude of the shifts of the respective levels is h ¯ Ω −∆ 2 ∆. Notice that for the bare atom-field states, state |li|N i has higher energy than |ui|N − 1i for positive detuning ∆, while for negative detunings the order is reversed.
CHAPTER 3. PHASE EVOLUTION
Corresponding atom-field state |ui|N − 1i |li|N i
25
Level shift of dressed state for detuning ∆>0 ∆<0 0 0 Ω +|∆| −¯h Ω −|∆| −¯ h 02 02 Ω −|∆| +¯h +¯h Ω +|∆| 2 2
Coupling always leads to new eigenstates which are energetically separated further than the bare atom-field states.
3.2.1
Using detuned lasers to shape trapping potentials
It is worth noting that the potentials created by far detuned laser beams used for creating spatially controlled perturbations of the potentials confining condensate clouds [10] [26], are based on exactly the same physical principle of the level shifts described above, expressed in Eq. (3.16). A far blue detuned Ω2 of the ground state (see beam, for example, causes a positive energy shift 4∆ appendix A) and thus behaves as a repulsive potential. A red detuned beam behaves as an attractive potential. Fundamentally, the potential, or “force”, arises by the mechanism of coherent scattering of photons from one mode (or k vector) to another, and one can see that for a net exchange of momentum between the atom and the radiation field, a spatially shaped, or non zero intensity gradient is required. This so-called dipole force is discussed in detail in [12].
CHAPTER 3. PHASE EVOLUTION
3.3
26
Phase Evolution in a driven Two State System
The temporal evolution of a system state can be easily expressed in terms of the eigenstates and eigenenergies. For the case of the initial population being fully in state |Ii, temporal evolution can be expressed as follows: e
−iHt h ¯
|Ii = e
−iHt h ¯
= αe
[α|N +i + β|N −i]
−iE+ t h ¯
|N +i + βe
−iE− t h ¯
(3.17)
|N −i
Now expressing the eigenstates |N +i and |N −i in terms of the bare atom-field states, equation (3.17) becomes = e
i∆t 2
≡
0
cos Ω2 t −
i∆ Ω0
0 sin Ω2 t |Ii −
0 Ω i sin Ω2 t |IIi Ω0
a(t)|Ii + b(t)|IIi
(3.18)
The magnitudes of the complex coefficients a(t) and b(t) relate to the populations in the respective states |Ii and |IIi, while the phase of the states is evaluated from the complex arguments of the coefficients. Conservation of total population is expressed as |a(t)|2 + |b(t)|2 = 1.
(3.19)
a(t) = |a(t)|eiΦ|Ii (t)
(3.20)
b(t) = |b(t)|eiΦ|IIi (t) ,
(3.21)
Now if we write
and it is easy to see that Φ|IIi evolves as ∆2 t, but with a π phase jump every period T = 2π , as illustrated in figures 3.3 and 3.4. The phase jumps occur when Ω0 the upper state population goes to zero. Averaging over an integral number of periods, we find Ω0 + ∆ t. (3.22) Φ|IIi (t) = 2 For state |Ii, it is easy to show that ∆t ∆ Ω0 t −1 Φ|Ii (t) = − tan tan . (3.23) 2 Ω0 2
CHAPTER 3. PHASE EVOLUTION
27
This solution is smooth (e.g. see figures 3.3 and 3.4). Over an integer number of periods it averages to Ω0 − ∆ t. (3.24) Φ|Ii (t) = − 2
a)
20 15
Phase angle [rad]
10 5 0 -5 -10 -15
b)
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
0.5
1
1.5
2
2.5 Time [t/T]
3
3.5
4
4.5
5
1
Relative populations
0.8
0.6
0.4
0.2
0
Figure 3.3: a) phase evolution b) population evolution for states |Ii (solid lines) and |IIi (dashed lines), driven by a cw field for the case Ω0 = 1, ∆ = 0.1, initial inversion w = −1. Time scales in terms of Ω2π0 . For large detunings ∆, i.e. ∆2 Ω2 , this phase slope is approximated by Ω2 t. This result is in agreement with the phase encoding that can be ΦI (t) = − 4∆ calculated directly for the far detuned regime, when the two-state system can be approximated as a single-state (ground state) system using the so-called “adiabatic elimination”. (We have outlined this in Appendix A.) Equations (3.22) and (3.24) describe the averaged phase evolution accurately for ∆ > 0 as well as for ∆ < 0. The equivalence of the gradients of
CHAPTER 3. PHASE EVOLUTION
a)
28
30 25
Phase angle [rad]
20 15 10 5 0 -5
b)
0
1
2
0
1
2
3
4
5
3
4
5
1
Relative populations
0.8
0.6
0.4
0.2
0
Time [t/T]
Figure 3.4: a) phase evolution b) population evolution for states |Ii (solid lines) and |IIi (dashed lines), driven by a cw field for the case Ω0 = 30, ∆ = 100, initial inversion w = −1. Time scales in terms of Ω2π0 . phase evolution with the shift of energy levels previously discussed in section 3.2 will be interpreted in the following section. Note that we have chosen the energy of atom-field state |Ii to be E1 = 0, which sets the origin of the energy scale. All level shifts and, as we will see, all phase gradients are measured relative to this energy level.
CHAPTER 3. PHASE EVOLUTION
3.3.1
29
Phase evolution for arbitrary initial inversions
Two interesting initial conditions to look at are the superpositions of atomfield states |Ii and |IIi that constitute the eigenstates |N +i and |N −i of the coupled system (equations 3.11). In these cases, the initial values for a0 (t) and b0 (t) are α and β (3.13) respectively for |N +i, and β and −α for |N −i. In the |N +i case, the gradients of phase evolution for both physical states |Ii and |IIi are constant and both are equal to the average described in equation (3.22). In the |N −i case, the gradients are constant and equal to (3.24) for both physical components |Ii and |IIi. Furthermore, the populations in states |Ii and |IIi remain constant during temporal evolution of the system. The above two cases show us that we can interprete the phase evolution for arbitrary initial populations a0 and b0 by expressing it as the evolution of a superposition of the system’s two eigenstates with relative components c and d as follows: e−
iHt h ¯
[c|N +i + d|N −i] = ce−
iE+ t h ¯
|N +i + de−
iE− t h ¯
|N −i
(3.25)
Expressing the eigenstates in terms of the physical states (3.11), this becomes i h i h iE+ t iE− t iE− t iE+ t (3.26) = cαe− h¯ + dβe− h¯ |Ii + cβe− h¯ − dαe− h¯ |IIi.
For pure eigenstates, either c or d go to zero and the phase evolves with either gradient E+ or E− , as observed, and populations in |Ii and |IIi remain constant. For all other combinations, the coefficients of Ii and |IIi are sums of complex numbers of different magnitudes (cα, cβ, dα, dβ) and differently evolving phases. The result is an oscillating population and an oscillating phase gradient, which averages to the value of the dominant component (cα, cβ,. . . ). This interpretation also gives an explanation for the phase jumps, observed when populations in state |Ii or |IIi go through zero: Addition of the evolving complex components of the state’s complex coefficient causes the sum to go through the origin of the complex plane, causing a π phase skip. If the magnitude of any of the two components is slightly larger than the other, phase evolution of the complex sum misses the origin in the complex plane and causes a smooth but very rapid phase change of nearly π, as can be seen, for example, for state |Ii in figure 3.3a.
Chapter 4 Numerical Implementation This section describes the implementation of the method that was used to solve numerically the nonlinear GP equations. The Fourier Transform method, which we outline below, was chosen because it is a well developed algorithm that was available within the group [16].
4.1
Dimensionless Units
To overcome problems of machine precision and numerical errors, which arise through truncations and roundings, it is worthwhile in a numerical implementation of a physical problem to transform the equations into “computational units”. For reasons of simplicity we will discuss the single component GP equation only, since it is trivial to extend this description to the coupled two component case. The GP equation, for a trapped BEC in two spatial dimensions, and in a gravitational field is ¯ h ¯2 2 ¯ 1 ∂Ψ ¯ − N U0 |Ψ| ¯ 2Ψ ¯ − mg y¯Ψ ¯ i¯h =− ∇ Ψ + m(ωx2 x¯2 + ωy2 y¯2 )Ψ ¯ ∂t 2m 2
(4.1)
Equation (4.1) is written in SI units (marked with bars on the symbols) and allows for nonisotropic traps by introducing different trap frequencies in x and y directions, such that ωy = λωx . ¯ into dimensionless (i.e. computational) quantities We transform t¯, x¯, y¯, Ψ
30
CHAPTER 4. NUMERICAL IMPLEMENTATION
31
by the transformations x¯ =
x0 x
y¯ =
x0 y
t¯ = t0 t ¯ = Ψ0 Ψ Ψ
(4.2)
in which t0 , x0 , y0 and Ψ0 carry dimensions. Substituting equations (4.2) into equation (4.1), we find that by setting t0 ≡ ω1x , equation (4.1) reduces to i ∂Ψ = i∇2 Ψ − (x2 + λy 2 )Ψ + iC|Ψ|2 Ψ + iGyΨ ∂t 4
(4.3)
where x0 ≡ C≡ G≡
q
h ¯ 2mωx
N U0 h ¯ ωx x20 mg x. h ¯ ωx 0
(4.4)
Notice that Ψ0 has the dimension x10 , as can be inferred from its normalisation equation Z∞ |Ψ|2 dx dy ≡ 1. (4.5) −∞
In the above definition of C, U0 has the dimensions of an area, so that C is nondimensional. The two-dimensional nonlinearity used in our simulations can be derived from the three-dimensional value using an integration method as used in [41], or by matching density peaks of corresponding 3D- and 2Dwavefunctions, as outlined in [29].
4.2
Interaction Picture
It is convenient to separate the right hand side of equation (4.3) into the part involving spatial derivatives, and the remainder. We write ∂Ψ ˆ +N ˆΨ = DΨ ∂t where ˆ = i∇2 , D
ˆ = − i (x2 + λy 2 ) + iC|Ψ|2 + iGy N 4
(4.6)
(4.7)
CHAPTER 4. NUMERICAL IMPLEMENTATION
32
This allows the Laplacian operator to be handled separately from effects of the nonlinear term in the GP equation. Since it is easy to account for the kinetic energy term (i∇2 ) alone, we transform into an interaction picture as follows ˆ
ΨI (τ, x, y) = e−D(τ −t) Ψ(τ, x, y),
(4.8)
where τ is the time at which the interaction picture separates from the normal picture. The GP equation (4.6) now becomes ∂ΨI (τ, x, y) ˆ −t) ˆ ˆ I (τ, x, y)ΨI (τ, x, y) = e−D(τ N Ψ(τ, x, y), (4.9) =N ∂t which is formally similar to an ordinary differential equation, and thus can be solved by Runge Kutta methods. The only remaining problem is the evaluation of the right hand side of equation (4.9), which is described in the following section.
4.3
Fourier Transform Method
For the sake of simplicity, the following description applies to the one-dimensional case. The extension to two spatial dimensions is relatively straightforward conceptually, but computationally intensive. The use of Fourier transforms (FTs) simplifies the evaluation of the right ˆ is diagonalised in this reprehand side of equation (4.9), since the operator D sentation. ˆ (τ, x)Ψ(τ, x) (as in [33]) Defining the Fourier transform of N θ(τ, f ) =
Z∞
ˆ (τ, x)Ψ(τ, x)ei2πf x dx, N
(4.10)
−∞
then the inverse transform will be Z∞
ˆ (τ, x)Ψ(τ, x) = N
θ(τ, f )e−i2πf x df.
(4.11)
−∞
It is now easy to calculate the right hand side of equation (4.9) as follows: R∞
ˆ
ˆ (τ, x)Ψ(τ, x) = e−D(τ −t) N
2
θ(τ, f )e−i(τ −t)∇ e−i2πf x df
−∞
=
R∞ h
−∞
θ(τ, f )e
−i(τ −t)(−i2πf )2
i
e−i2πf x df
(4.12)
CHAPTER 4. NUMERICAL IMPLEMENTATION
33
To summarize the whole calculation, if one defines ˜ f ) = θ(τ, f )e−i(τ −t)(−i2πf )2 , θ(τ,
(4.13)
the calculation of (4.9) can be described as follows: • Calculate θ(τ, f ) by simple FT. 2 ˜ • Multiply each coefficient with e−i(τ −t¯)(−i2πf ) to obtain θ.
ˆ ˆ Ψ. • Carry out inverse FT on θ˜ to obtain e−D(τ −t) N
4.4
Runge-Kutta Method
The Fourier method described in the previous section allowed us to conveniently evaluate the function f (ΨI , t) representing the derivative in equation (4.9). ∂ I Ψ (t) = f (ΨI , t) (4.14) ∂t The Runge-Kutta (RK) method provides a well tested “workhorse” for integrating differential equations, and a 4th order implementation to equation (4.9) is convenient. For each timestep the solution is advanced by the RK4 algorithm, the derivatives are evaluated four times in different places; once at the initial point, twice at the center between the step points and once at a trial endpoint. The “true” endpoint, the solution of the timestep, is then calculated as a weighted sum of the four derivative evaluations as follows. Defining [33] f (ΨI , tn )
k1 = k2 =
f (ΨI + k1 ∆t , tn + 2
∆t ) 2
k3 =
f (ΨI + k2 ∆t , tn + 2
∆t ) 2
(4.15)
k4 = f (ΨI + k3 ∆t, tn + ∆t), the weighted sum is k1 k2 k3 k4 + + + Ψ (tn+1 ) = Ψ (tn + ∆t) = Ψ (tn ) + ∆t. 6 3 3 6 I
I
I
(4.16)
An advantage of using the interaction picture is that by chosing the origin of the picture (τ ) at each time step to be at the centre of the step, the factor
CHAPTER 4. NUMERICAL IMPLEMENTATION
34
ˆ
of e−D(τ −t) disappears from the derivative, so that FFTs are saved in two of the four evaluations of the differential necessary in each timestep.
4.5
Numerical Limits
Our computer simulations are limited, as usual, by the available resources of CPU power, memory and data storage. The main impact on our work is on the choice of grid- and stepsize. We discuss the implications of these limitations in the following sections.
4.5.1
Limitations imposed by grid resolution
The velocity of a coherent BEC superfluid is a linear function of its phase gradient (see section 6.5). In our computational units, the following relation holds: v = 2 ∇S (4.17) where S is the phase of a complex wavefunction of the form |Ψ|eiS . To simulate 40 30
v
10
velocity
20
0 -10 -20 -30 0
5
10
15
20
Phase Gradient
25
30
35
40
∆
-40
S
Figure 4.1: Centre of mass condensate velocity vs. phase gradient. Aliasing occurs at the critical velocity limit v=19.95. Solid line for a Gaussian eigenstate and dotted line for a plane wave with zero bandwidth. In the plane wave case the velocity changes sign exactly at the limit.
CHAPTER 4. NUMERICAL IMPLEMENTATION
35
a BEC cloud moving through computational space, we apply a phase gradient ∇S across a BEC eigenstate and propagate it in time. Now, if 20 spatial units of computational space are stretched out over 128 gridpoints, the spatial resolution is ∆x = 20/127 = 0.1575 . . . , and the phase change ∆S from one gridpoint to the next must be less than π for a positive velocity, since a phase greater than exactly π can be ambiguously interpreted. Additionally, one must consider the problem of aliasing. The Fourier spectrum associated with the grid 1 ) must be large enough to accommodate the Fourier components (kmax = ∆x of the BEC. Figure 4.1 the centre of mass velocity of the simulated BEC cloud is plotted vs. the applied phase gradient (solid line). Since the BEC cloud has a near-Gaussian Fourier spectrum (Gaussian for a C=0 eigenstate), aliasing occurs in a region of a phase gradient between 18 and 22, while a plane wave with vanishing bandwidth exhibits a sharp aliasing “transition” at a phase gradient of exactly π/0.1575 · · · ≈ 19.95 (dotted line). This demonstrates that simulations of dynamic processes, i.e. output coupling with acceleration under gravity, require minimum grid resolutions, depending on the expected maximum velocities, in order to yield physical results. To allow high velocities in output coupling simulations, i.e. velocities of up to v=80, we used spatial resolutions of ∆x = 20/255 = 0.078 . . .
CHAPTER 4. NUMERICAL IMPLEMENTATION
4.5.2
36
Limited simulation stepsize
It is desirable to use simulation timesteps as large as possible, to limit computational time. The RK based algorithm described in sections 4.3 and 4.4 has been shown to be accurate to 4th order in step size [16], and thus ensures accurate solutions with larger step size than in a standard split step method. Nevertheless, the step size cannot be chosen too large. The temporal phase gradient of an eigenstate (compare equation 2.15), is the chemical potential µ as defined in equation (2.16). Forqtwo-dimensional C , as shown in wavefunctions, the chemical potential is evaluated as µ = 2π section 2.1.1. Thus, the aliasing limit for a C=200 eigenstate in 2D is reached for stepsizes greater than ∆t = π/µ = 0.54. In figure 4.2, the centre of mass (c.o.m.) velocity of a BEC eigenstate with a constant phase gradient of 10 and a corresponding velocity of vx = 20 is plotted vs. temporal step size; the initial state was propagated for a single time step, after which the c.o.m. velocity was calculated. With increasing stepsize, the aliasing limit is reached at ∆t = 0.54 and the measured c.o.m. velocity is reversed, i.e. phase rotations larger than π within one timestep are interpreted as negative phase rotations. With further increasing stepsize, measured velocity reverses several times while tiny differences in the condensate’s phase gradient (due to adiabatic expansion in absence of a trapping potential) become increasingly important, leading to a “washing out”, and causing the average c.o.m. velocity eventually to go to v = 0. When practically determining the step size for dynamic numerical simulations, we must also consider effects of changing wavefunction amplitudes, which can cause large temporal phase gradients. Furthermore, although the RK method works very well, it cannot accurately follow large gradients close to the aliasing limit. Thus, for all practical purposes, the step sizes we chose were many orders of magnitude smaller than the step size at the aliasing limit in the example described above. For example, while the aliasing limit for a C=200 2D eigenstate is at 12 steps per trap cycle Ttrap = 2π, we typically used step sizes as small as 1000 or 2000 steps per Ttrap .
CHAPTER 4. NUMERICAL IMPLEMENTATION
37
20
15
cloud c.o.m. velocity
10
5
0
−5
−10
−1
10
0
1
10
10
2
10
temporal step size
Figure 4.2: Illustration of simulation failure for too large stepsizes. The temporal phase gradient µ (chem. potential) of a C=200 condensate 2D eigenstate reaches the aliasing limit for a step size of ∆t=0.54.
Chapter 5 RF Output Coupler in 2D This chapter describes simulations of continuous wave (cw) radio frequency (rf) output coupling in a two-state system in two dimensions. While cw output coupling has been simulated in one spatial dimension [3], a key difference in the present work is the addition of gravity. The setup of a cw output coupler for BEC has been called an “atom laser” and has been realized experimentally [8]. While the term “atom laser” appears to be an appropriate description for a “beam” of matter with properties similar to optical laser beams, most importantly with spatial coherence, the notion is not universally agreed upon and still a matter of controversy. A strict definition has been proposed by Wiseman [39]. The effects of gravity open up a range of new behaviour, but also give rise to certain computational problems. For example, gravity provides an easy way to get rid of outcoupled untrapped populations by just letting it fall out of the trap, thus forming a “beam”. However, the output “beam” will now inevitably reach the edge of the computational grid in the consequent acceleration. Furthermore, increased velocities in the simulation demand higher computational grid resolutions and smaller temporal stepsizes.
5.1
Output Couplers
To obtain Bose-Einstein condensation, atoms are trapped magnetically, utilizing the Zeeman shift of their hyperfine states. Different spin states, i.e. MF = 1, 0, −1, have different trap properties since the force, which the trap38
CHAPTER 5. RF OUTPUT COUPLER IN 2D
39
ping field exercises on the atoms, is F ∝ MF
∂B(r) , ∂r
(5.1)
where B is the magnetic field and r the distance from the field zero, i.e. the trap centre. While, for example, atoms with MF = 1 are trapped, atoms with MF = 0 are indifferent to a magnetic field, i.e. they are untrapped, and atoms with MF = −1 are “anti-trapped” and repelled from the trap. The energy difference between trapped and untrapped states (we will use MF = 1 and MF = 0 for this purpose from now on) is in the radio frequency (rf) regime, so that two basic strategies for an output coupler are viable: Population can be transferred and deposited into an untrapped state by simple radio frequency coupling, and left there to evolve from the same spatial region (rf coupler), or it can be transferred and actively ejected from the spatial region of coupling (Raman couplers). We will only discuss in this work the simple rf coupler case, as has been experientally realised in [30] and [8].
5.2
Spatial Effects in Coupling
A coupled two-state system is described by the following set of GP equations r2 Ω ∂ψ1 = i∇2 ψ1 − i ψ1 + i ψ2 − iC |ψ1 |2 + w|ψ2 |2 ψ1 ∂t 4 2 2 ∂ψ2 r Ω = i∇2 ψ2 − ik ψ1 + i ψ1 − i∆ψ2 − iC |ψ2 |2 + w|ψ1 |2 ψ2 . ∂t 4 2
(5.2) (5.3)
Ω is the Rabi frequency, i.e. the coupling strength, ∆ is the detuning of the coupling field for the centre of the trap, and k is the relative trap parameter of the second component. For an untrapped second state, k is equal to zero. Other physical systems with different magnetic momenta, e.g. in the MF = [2, . . . , −2] manifold, can be simulated by adjusting this parameter appropriately. w is the relative inter-state collision strength, the ratio between the inter-component interaction and the self interaction. The magnitude of w has been a matter of discussion and was incorrectly set at w = 2 in reference [3]. Recent experimental results in 87 Rb vapour [19] show scattering length ratios of a1 : a12 : a2 :: 1.03 : 1 : 0.97, where a1 and a2 are the scattering
CHAPTER 5. RF OUTPUT COUPLER IN 2D
40
lengths within the first and second component, and a12 is the intercomponent scattering length. Since all three values are very close to each other and approximately equal to one, using a single value of w = 1 for all three scattering lengths for the numerical simulations in this work seems to be justified. In our simulations of output coupling, only the first state is trapped by a harmonic potential, while the second state is untrapped. Thus the detuning between the two components is not uniform at all spatial positions in the cloud. In the absence of the ∇2 term, i.e. within the TF approximation, we can define an effective detuning ∆eff 5.1): 1 ∆eff (r) = −C(w − 1) |Ψ2 (r, t)|2 − |Ψ1 (r, t)|2 + ∆ − r2 . 4
(5.4)
which is obtained by taking the difference of the diagonal terms of (5.2) and
E Trapping potential first state
1 2 r 4 Untrapped second state
r
Figure 5.1: Spatial dependence of energy difference between trapped and untrapped states. (w=1) (5.3). ∆eff reduces to ∆ at the centre of the trap, when w=1, and increases radially, as we illustrate in figure 5.1. As we will show in detail in section 5.5, the coupling strength depends on the detuning. We define a spatially dependent Rabi frequency: q Ω(r) = ∆2eff (r) + Ω20 . (5.5)
The spatial dependence of the effective detuning (and hence the generalized Rabi frequency) lead to “spatial effects”: Coupling is no longer homogeneous
CHAPTER 5. RF OUTPUT COUPLER IN 2D
41
throughout the cloud and is more effective at certain spatial locations. As shown in [3], for spatial effects to be negligible, i.e. for homogeneous population transfer, we require that Ω(r) is approximately constant over the whole cloud and thus approximately equal to the trap centre cycling frequency Ω0 . This requirement can be expressed as Ω(r) ≈ 1. Ω0
(5.6)
Using the Thomas-Fermi expression (see section 2.1.1) 1 C|Ψ|2 = µ − r2 , 4
(5.7)
we see that for w = 2 as in [3], ∆eff is reduced to an initial value (for vanishing population in the second state) of 1 ∆eff (r) ≈ µ + ∆ − r2 . 2
(5.8)
For w = 1, we have no net contribution of self energies, and we get 1 ∆eff (r) ≈ ∆ − r2 . (5.9) 4 √ √ Evaluating (5.6) in the range −2 µ < r < 2 µ, where the wavefunction is nonzero within the Thomas-Fermi limit, the condition for spatially independent Rabi cycling (i.e. no spatial effects of coupling) for both values of w (w=1,2) becomes ∆2 + Ω20 µ(µ + 2|∆|). (5.10) For the remainder of this thesis, we will always assume a value of w = 1. Since inequality 5.10 is not easy to interpret, we illustrate it for a 2D C=200 eigenstate in figure 5.2. Curves are plotted for an inequality factor of 10, and also for a factor 1. This provides a rough impression of the shape of the parameter regime, in which spatially dependent coupling can be expected, but it should be interpreted with appropriate care. Spatial effects will become more significant with smaller ∆ and Ω. Spatially homogeneous coupling can occur for two distinct reasons: Either large bare Rabi frequencies or large detuning ∆, dominating the spatial variation in the effective detuning ∆eff .
CHAPTER 5. RF OUTPUT COUPLER IN 2D
42
60 50
Coupling Ω0
40 30 20 10 0 -150
-100
-50
0
50
Detuning ∆
100
150
Figure 5.2: Illustration of regimes in which spatial effects are expected. Parameter µ=5.64 represents a 2D C=200 eigenstate. Plotted are functions ∆2 + Ω20 = µ2 + 2µ|∆| (solid line) and ∆2 + Ω20 = 10(µ2 + 2µ|∆|) (dashed line). Equation (5.10) is valid for parameters under the dashed curve.
−3
3
x 10
−3
x 10 3
a)
b)
2.5
2.5 2
2
1.5
|Ψ|
|Ψ|
2
2
1 0.5
1.5
1
0 10 0.5
5 10
y
0
5 0
−5 −5 −10
−10
x
0 −8
−6
−4
−2
0
2
4
6
8
x
Figure 5.3: Two dimensional normalised C=200 eigenstate. a) Spatial probability density |Ψ|2 . b) Profile along y=0 (solid line) and Thomas-Fermi approximation (dashed line).
CHAPTER 5. RF OUTPUT COUPLER IN 2D
5.3
43
Influence of gravity on trapped BEC
Gravity adds another potential to the coupled GP equations. This potential is linear and vertically increasing, and it leads to an offset of the condensate’s equilibrium position from the center of the magnetic trapping potential. This gives rise to new behaviour of output couplers. The coupled nonlinear Gross-Pitaevskii equations for a two-level system under the influence of gravity are (in computational units): Ω ∂ψ1 = i∇2 ψ1 −iVtrap ψ1 +i ψ2 −iC |ψ1 |2 +w|ψ2 |2 ψ1 −iGyψ1 (5.11) ∂t 2 ∂ψ2 Ω 2 = i∇ ψ2 −ikVtrap ψ2 +i ψ1 −i∆ψ2 −iC |ψ2 |2 +w|ψ1 |2 ψ2 −iGyψ2 (5.12) ∂t 2
G is the acceleration of gravity scaled into dimensionless units in equation (4.4). All other symbols have been discussed in the previous section. The equilibrium offset of the ground state condensate centre of mass, relative to the center of the magnetic trapping potential is easily calculated. We have the vertical component of the trapping potential Vtrap,v = − 41 y 2 and the linear vertical potential of gravity Vgrav = −Gy. The condensate equilibrium position is where the joint potential is at a minimum, i.e. where the forces resulting from the harmonic trapping potential and gravity are equal. This is easily found to be at y = 2G. In our simulations of output coupling, we use nonlinear two-dimensional initial states with values of C=200, unless stated differently. We plot in figure 5.3a a true 2D wave function, and in figure 5.3b, we compare a profile of this eigenstate with the corresponding TF approximation. We can see that the latter is a very accurate approximation, except at its edges.
5.4
Localized Output Coupling
While gravity shifts the equilibrium position of the condensate, the spatial location of the coupling is still determined by the effective detuning. (Notice that gravity effects both components equally and thus does not alter the expressions for ∆eff of equation 5.4.) Output coupling occurs along the surface of a sphere (a circle in the 2D representation) centred at the trap origin,
CHAPTER 5. RF OUTPUT COUPLER IN 2D
44
i.e. along a magnetic equipotential surface, and this is no longer symmetric about the condensate due to its offset of y = 2G. This is illustrated in figure 5.4. Population coupled into the untrapped second state is released from this “cut” along the magnetic equipotential line to fall down under the influence of gravity. y magnetic equipotential lines of trapping field
x
g
condensate cloud
Figure 5.4: Equlibrium position of condensate cloud in a harmonic trapping potential under the influence of gravity. The position of the “cutting” magnetic equipotential line can be determined by setting the effective detuning in equation (5.4) to ∆eff = 0. For w = 1 this becomes independent of populations and nonlinearity, and the solution is √ r = 2 ∆.
(5.13)
For arbitrary values of w, the cutting position depends on the component population densities according to r=2
p
∆ − C(w − 1)[|ψ2 |2 − |ψ1 |2 ].
(5.14)
In figure 5.5, the relationship between spatial position of coupling and detuning is illustrated. (Data acquired in simulations in the TF regime, i.e.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
45
0 2 4
y
6 8 10 12 14 16
0
5
10
15
20
25 detuning
30
35
40
45
50
∆
Figure 5.5: Spatial region of maximum coupling in a profile along the x=0 √ axis. Function y = 2 ∆ (solid line). Simulation data points for w=1 (crosses) and data points for w=2 (stars). Condensate equilibrium position at y=10 for G=-5. ∇2 term ignored.) For w=1, all data points obtained from simulation-runs lie on the theoretical curve corresponding to Eq. (5.13). For w=2 the situation is slightly different. Since coupling is relatively weak, population in state ψ 2 remains small and equation 5.14 can be approximated by r=2
p
∆ + C|ψ1 |2 .
(5.15)
Higher self-energy increases the effective detuning between the two components at a certain spatial position (5.4). Since coupling occurs along a line of constant detuning ∆, this moves the coupling region to a spatial position where an increased trap potential compensates the shift, and coupling occurs at a different spatial position, i.e. at increased y. This shift is dependent on population densities, and thus, the w=2 data points edge closer to the line, which represents w=1, with increasing distance y from the condensate centre at density maximum at y=10 in figure 5.5. Figure 5.6 shows the temporal evolution of an output coupled population stream in the weak coupling limit. The initial state (a) is left almost unchanged during the time of the simulation, while a region of the cloud is subject to Rabi
CHAPTER 5. RF OUTPUT COUPLER IN 2D
−20
−20 a)
−10
46
−2
b) −10
t = 0.42412
10
t = 0.42412 −3
10 0
0
10
10
20
20
30
30
10
40
40
10
50
50
−4
10
y
−5
10
−6
−7
−8
10 60 −20
−10
0
10
20
−20
60 −20
0
10
20
−20 c)
−10
−10
10−2
d) −10
t = 1.3666
t = 2.5447
10−3 0
0
10
10
20
20
30
30
10−6
40
40
10−7
50
50
y
10−4
60 −20
−10
0 x
10
20
60 −20
10−5
10−8
−10
0 x
10
20
Figure 5.6: Population densities in condensate states |1i and |2i for cw output coupling in a gravitational field. (a) state |1i, (b) state |2i at the same time t as in (a). A temporal sequence of the evolution of state |2i is shown in (b), (c) and (d), plotted in a logarithmic density scale. (Simulation parameters: Ω0 =0.2, ∆=40, G = −5, C = 200, 256 by 512 grid.)
CHAPTER 5. RF OUTPUT COUPLER IN 2D
−20
y
b) −10
t = 0.42412
0
0
10
10
20
20
30
30
40
40
50
50
60 −20
−10
0
10
20
−20
y
π/2
0
−π/2
−π −10
0
10
20
π
−20 d) −10
t = 1.3666
0
0
10
10
20
20
30
30
40
40
50
50
60 −20
t = 0.42412
60 −20
c) −10
π
−20 a)
−10
47
−10
0 x
10
20
60 −20
t = 2.5447
π/2
0
−π/2
−π −10
0 x
10
20
Figure 5.7: Phase of condensate states |1i and |2i for the same simulation as in figure 5.6. Phase plotted for population densities larger than 10−8 .
CHAPTER 5. RF OUTPUT COUPLER IN 2D
48
cycling with the second state, shown in (b) at the same time t as (a), and in (c) and (d) at later times. (At a Rabi fequency of Ω=0.1, a full Rabi cycle takes ∆t=20π, thus the small effect on the first state.) In figure 5.7 we plot the phase of the two wavefunctions for the same case as in figure 5.6. Gravity accelerates the outcoupled population downwards, leading to an increasing phase gradient. At the same time, self energy and interaction energy (mutual repulsion between the states) cause the untrapped cloud to expand and accelerate radially. Overlaying these two gradients leads to the “interference” fringes most clearly seen in (d). A magnification of the interesting region is shown in figure (5.8). It has been verified that this is not an artefact of printer resolution by taking the population in question and evolving it after subtracting the phase gradient caused by gravity. Note that in figures (5.7) and (5.8) we have plotted phase only for population densities greater than 10−8 in front of a uniform background.
π
25
30
π/2
y
35
40
0
45
−π/2 50
55 −15
−10
−5
0 x
5
10
15
−π
Figure 5.8: Magnification of figure 5.7d. Population falling under gravity (strong vertical phase gradient) while undergoing free expansion (small radial phase gradient) leads to an “interference like” phase plot.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
5.5
49
Coupling Bandwidth
In the previous sections we have shown how localized spatial coupling can be accomplished by controlling the coupling field strength and detuning under appropriate trapping conditions. We have also seen that the coupling bandwidth plays an important role in output coupling simulations by determining the spatial width of the coupling region, i.e. very high bandwidths lead to homogeneous coupling through power broadening. For a precise engineering of output coupled “beams”, controlled small coupling bandwidths are essential.
5.5.1
Power Broadening
The natural linewidth of our transitions depends on upper state lifetime τ , leading to a Lorentzian shaped spectrum of |E(ω)|2 ∝
1 ∆2 +
γ2 4
(5.16)
for detunings ∆ = ω − ω0 and with γ = 1/τ . Compared with optical dipole transitions, the spontaneous decay rate for magnetic dipole transitions between different MF states is negligibly small, even on much larger timescales than that of our simulation. Furthermore, collisional coupling at low densities is also negligible, and thus power broadening is the main broadening mechanism. It is well known from the Rabi solutions (e.g. see chapter 2.2.3 or [25]) that the maximum upper state population in a coupled system undergoing Rabi cycling with a detuning ∆ is pmax =
Ω20 Ω20 = ∆2 + Ω20 Ω0 2
(5.17)
This leads to a Lorentzian shaped curve of maximum population vs. detuning, with a HWHM value ∆ωH given by the Rabi frequency ∆ωH = Ω0 .
(5.18)
In the output coupler described in the previous sections, the detuning ∆ in (5.17) is replaced by the spatially dependent effective detuning ∆eff (eq. 5.4). √ The position of the centre of coupling was found to be r = 2 ∆. Using the HWHM condition ∆eff = ∆ωH = Ω0 , (5.19)
CHAPTER 5. RF OUTPUT COUPLER IN 2D
50
4 Spatial FWHM of coupling region
3.5 3 2.5 2 1.5 1 0.5 0
0
2 4 6 8 10 12 Distance of centre of coupling region from trapping field zero
14
Figure 5.9: Width of coupling region plotted vs. distance of the resonant √ coupling region from the trapping magnetic field zero (r=2 ∆). Widths calculated from equation (5.20) for a purely power broadened transition. Curves are plotted for Ω=0.2, 0.4, 1, 2 (solid, dot-dashed, dashed and dotted lines respectively). we calculate the full spatial width of the (FWHM) coupling region ∆y as: p p ∆+Ω0 − ∆−Ω0 for ∆ > Ω0 , (5.20) ∆y = 2 p for − Ω0 < ∆ < Ω0 , ∆y = 2 ∆+Ω0 and ∆y = 0
for ∆ ≤ −Ω0 .
Figure 5.9 illustrates how the power broadening converts into a region of spatial coupling for various choices of the centre of the coupling.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
5.5.2
51
Broadening through finite Pulse Times −20
−20
a) −10
−2
b) −10
t = 2.4976
10
t = 2.4976 −3
10 0
0
10
10
20
20
30
30
10
40
40
10
50
50
−4
10
y
−5
10
−6
−7
−8
10 60 −20
−10
0 x
10
20
60 −20
−10
0 x
10
20
Figure 5.10: Population density in condensate states |1i (a) and |2i (b) on a logarithmic scale, illustrating coupling linewidth (gradually approaching the power broadened limit) for long pulse-time at t = 2.4976. Simulation in the TF limit. Ω=0.2, ∆=40, G=-5, C = 200. In this and the previous sections, we develop a concept for interpreting results of dynamic coupling bandwidths, which we find in our simulations. Note that these are not models on which our simulation code is based, since the GP equation takes care of all this. We are merely interpreting the results, solutions of the GP equation. The most important contribution to coupling linewidths in our simulations originates from the pulse-like nature of our coupling field. It has a finite pulsetime τpulse starting at t=0, leading to a spectral FWHM width of Γpulse =
2π τpulse
.
(5.21)
Neglecting the insignificantly small natural linewidth, power broadening and pulse spectral width combine to a linewidth Γ which is calculated from a convolution of the spectra of the two separate broadening mechanisms. For all
CHAPTER 5. RF OUTPUT COUPLER IN 2D
52
practical purposes in this work, it can be approximated by Γ ≈ Γpower + Γpulse .
(5.22)
(Note that the GP equation provides an accurate description of the broadening phenomena without approximations in any case.) Pulse spectral width is the reason why broad coupling over the whole cloud is observed at the early times of a coupling simulation. As the coupling (pulse-) time increases, bandwidth decreases significantly, ultimately approaching the limit of a purely power broadened signal. Figure 5.10 illustrates a very narrow bandwidth, coupling an arc of high population density from the first into the second component. Pulse time in this simulation is approaching half a trap period (which is Ttrap /2=π) and the simulation is performed in the TF regime to prevent the second state component to fall down. In figure 5.11 we illustrate the temporal evolution of the measured spatial bandwidth. It is evident that finite pulse time constitutes the main broadening mechanism on a timescale of full trap cycles. The effect of the bandwidth due to transient effects is also dramatically illustrated in the case of negative detuning. In this case, ∆eff is not resonant at any spatial point, and thus we do not expect long term coupling. Figure 5.12 shows a temporal sequence of the evolution in the second (untrapped) state: Initially, after turning on the coupling cw field, coupling occurs throughout the whole cloud despite the negative effective detuning. With increasing interaction time, the bandwidth of the coupling field decreases significantly, leading to a very narrow energy range of coupling, which is not resonant at any position in the cloud. The system eventually settles into a state of constant low coupling bandwidth where rapid, small amplitude Rabi cycling occurs between the two states without significant output coupling. The part of the condensate that was output coupled in the initial high bandwidth period falls down under the influence of gravity forming a low density crescent shaped pulse as its selfenergy causes it to expand (reminiscent of the original MIT output coupler [30]). Repetitive crescent shaped pulses like the above are obtained when a resonant coupling field is operated in pulse mode [22]. Due to inter-state repulsion with the falling outcoupled population, the parent cloud, i.e. state |1i component, experiences a small “recoil” directed upwards. Because of this, the parent cloud undergoes small amplitude oscillation in its trap. However,
CHAPTER 5. RF OUTPUT COUPLER IN 2D
53
25
20
∆y
15
10
5
0
0.5
1
1.5
2
2.5
t
Figure 5.11: Temporal evolution of the spatial FWHM coupling bandwidth ∆y of the simulation in figure 5.10 (solid line). The dashed line represents the bandwidth component caused by finite pulse time (5.21) and the dotdashed line at ∆y=0.063 represents the small bandwidth contribution of power broadening at a coupling strength of Ω=0.2 at a detuning of ∆=40. in the low coupling strength regimes of our simulation (i.e. low Rabi frequencies and small outcoupled densities), this effect, which is described in [22], was too small to measure. In regimes of stronger coupling, we did not obtain any conclusive results. Increasing energy selectivity of the coupling field (i.e. bandwidth narrows), as time progresses is also evident in Fig. 5.6 on page 46, where it leads to fringes of shrinking width as population undergoes quick Rabi oscillations between the two states, with a spatially dependent Ω(r), Eq. (5.5), in (b) and (c). This effect combined with vertical acceleration of outcoupled population also leads to the density “pulse”-like fluctuations in the output coupled stream, which are visible in (d). These “pulses” are an effect of the transient due to turning on the coupling field, (i.e. high coupling bandwidth, or low energy selectivity) and subside as time progresses.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
−20
y
10−5
−20 a)
−10
54
b) −10
t = 0.1885
0
0
10
10
20
20
30
30
40
40
t = 0.6597 −6
10
10−7
−8
10 50
50
60 −20
−10
0
10
20
−20
60 −20
y
0
10
20
10−5
−20 c)
−10
−10
d) −10
t = 1.3666
0
0
10
10
20
20
30
30
40
40
t = 2.3091 −6
10
10−7
−8
10 50
60 −20
50
−10
0 x
10
20
60 −20
−10
0 x
10
20
Figure 5.12: Output coupling with negative detuning. Pulse output is due to temporal broadening at early times. Temporal evolution sequence of probability density of component |2i in figures (a)-(d). (Parameters: Ω=0.2, ∆=-20, C = 200, 256 by 512 grid)
CHAPTER 5. RF OUTPUT COUPLER IN 2D
5.5.3
55
Broadening using Noisy Signals
We have considered the effect of noisy cw coupling fields on the output coupling process using a simple stochastic model. The noise was modelled as small time dependent fluctuations of the driving frequency, implemented as random phase jumps between temporal simulation steps. We use a complex amplitude z(t) = z0 eiω(t)t ,
(5.23)
where the frequency exhibits small time dependent fluctuations α(t) around the mean ω0 . ω(t) = ω0 + α(t) (5.24) Formally, this can be represented by the Stratonovich equation [17] i h p dz = i ω0 dt + 2γ dW (t) z
(5.25)
where γ is a measure of the strength of the noise on our complex amplitude and W (t) is the Wiener function. The solution to (5.25) is z(t) = z0 ei[ω0 t+
√ 2γW (t)]
.
The time dependent frequency fluctuation α(t) becomes √ 2γ α(t) = W (t) t
(5.26)
(5.27)
and (5.26) gives a frequency spectrum S(∆) =
I γ , π ∆2 + γ 2
(5.28)
with a HWHM bandwidth of ∆ω = γ. Randomness in the simulation was represented as a Wiener function W (t) with expectation value hW (t)i = 0 and variance hW (t)2 i = t. Numerical simulation of this process was implemented using the MATLAB function randn(), which is a timer seeded Gaussian pseudo random number generator (PRNG), with variance v = σ 2 = 1. By multiplying the output numbers from randn() √ with ∆t, where ∆t is the simulation timestep, we obtain a Wiener random process function with variance hW (t)2 i = t.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
−20
10−4
−20 a)
−10
56
b) −10
t = 0.7697
t = 1.2881 −5
10 0
0
10
10 −6
y
10 20
20
30
30
40
40
50
50
60 −20
−10
0
10
20
−20
10−8
60 −20
−10
0
10
20
10−4
−20 c)
−10
−7
10
d) −10
t = 1.8064
t = 2.3405
10−5 0
0
10
10
y
10−6 20
20
30
30
40
40
50
50
60 −20
−10
0 x
10
20
60 −20
10−7
10−8
−10
0 x
10
20
Figure 5.13: Output coupling with a noisy cw field. (a)-(d) gives a temporal sequence of the probability density of untrapped component |2i. (Parameters: Ω=0.2, ∆=10, C = 200, γ=10, 256 by 512 grid)
CHAPTER 5. RF OUTPUT COUPLER IN 2D
57
In output coupling simulations, noise destroys the coherence of the outcoupled beam. Noise leads to density fluctuations and stochastic phase evolution along the output “beam” (see Fig. 5.13) and thus causes a loss of the properties that are characteristic for “atom lasers” [39]. To demonstrate this in more detail, we compared the self correlation functions Γ(y) = Ψ∗2 (y 0 )Ψ2 (y) (5.29) of the outcoupled state at time t = 2.3405 of a simulation without noise, to various realizations of a simulation at the same noise level (e.g. one of the realizations is in Fig. 5.13d). In each case we selected a profile at x = 0. Figure 5.14a shows the function Γ(y) for a simulation with noiseless rf coupling field, with the fixed value y 0 = 16.6 (vertical mark). We note that y 0 was chosen somewhat arbitrarily, but any y 0 within the outcoupled population yields similar results, a smooth almost uniform correlation along the beam. The situation is different for a simulation with a noisy coupling field, as shown in figure 5.14b. The magnitude of the correlation function exhibits large fluctuations which are also apparent in a density profile (e.g. see Fig. 5.13d). The loss of coherence, which is not illustrated very well by the magnitude function, becomes evident by averaging Γ(y) over a number of simulations with the same noise level. All realisations have different random dephasings which eventually average to zero when enough individual runs are included. Figure 5.14c shows an average of 25 realisations of the same simulation. The width of the peak at y 0 can be regarded as the coherence length of the “beam” as we will see shortly. Figure 5.14d shows the same data but using y 0 = 33.855. The correlation peak decreases as the population density decreases and widens because of the spread due to gravitational acceleration. A rough estimate of the coherence lengths can be made as follows: A noise level of γ = 10 has a bandwidth of ∆ω = 10, as found in equation (5.28). This bandwidth relates to a spectral time τpulse = 2π/γ = 0.628. Now, using simple Newtonian mechanics, we can calculate the time t taken by the outcoupled population to fall from the coupling region to position y 0 . We expect the beam to be coherent between times t + τpulse /2 and t − τpulse /2, which is easily converted to a corresponding spatial distance along the beam.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
58
−5
2
x 10
a) |Γ|
1
0 1
0 −5 x 10
10
20
30
40
50
60
b)
|Γ| 0.5
0 1
0 −5 x 10
10
20
30
40
50
60
c) 〈|Γ|〉 0.5
0 1
0 −5 x 10
10
20
30
40
50
60
d)
〈|Γ|〉 0.5
0
0
10
20
30 y
40
50
60
Figure 5.14: Spatial coherence analysis of outcoupled “beams”; (a): No noise present. (b): Single realization of simulation at noise level γ = 10. (c): Average over 25 different realisations of (b). (d): same as (c) but with different y 0 . The positions of the y 0 values used are marked by vertical solid lines. Profiles at x = 0 of simulations with Ω0 = 0.2, ∆ = 10, C = 200, γ = 10 in a 256 by 512 grid at t = 2.3405.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
59
For figure 5.14c, the calculation yields ∆y ≈ 6.3 and for figure 5.14d we get ∆y ≈ 8.4, which agrees with the graphs.
5.6
Limits of the Simulation
There are significant difficulties in simulating output couplers compared to simulating BEC in a trap. First, the output coupled beam, falling under gravity, can reach the critical speed where reflections due to numerical aliasing occur, as described in section 4.5. This can be resolved by increasing the grid resolution, and typically we used grids of 256 (horizontal) by 512 (vertical). The ultimate restriction on the size of grid that can be used is set by the available physical memory of the machine. The grid matrix must fit into memory, since many matrix operations and Fourier transforms are performed on it in every timestep of the simulation. A typical (i.e. 256 by 512 grid) simulation running in MATLAB on a Pentium Pro 200 CPU with 128MB RAM took approximately eight hours to propagate through a full trap cycle Ttrap . We illustrate in figure 5.15 how aliasing in a coarse grid leads to the reversal of the outcoupled “beam”. The critical velocity of v=40 is reached at around y=35. Consequently, the reversed beam interferes with the initial beam. The other, more difficult problem was that of output coupled material hitting the edges of our simulation grid and reflecting back– clearly an unphysical scenario. In an ideal case, material “falling” out of the grid should just disappear, without causing any reflections or refractions that can interact with the initial populations. We have mitigated the effects of these reflections by two separate strategies. The simplest was to monitor population density in the grid border regions as the simulation progressed, stopping the simulation run when a certain limit was reached. Typically, we stopped the simulations when the combined relative density of both states in a 5% rim around the grid exceeded 0.001. The other strategy was to add a thin absorbing layer at the edge of the grid, as we discuss in the next section.
CHAPTER 5. RF OUTPUT COUPLER IN 2D
−20
−20 a)
−10
60
−2
b) −10
t = 1.3666
10
t = 1.3666 −3
10 0
0
10
10
20
20
30
30
10
40
40
10
50
50
−4
10
y
−5
10
−6
−7
−8
10 60 −20
−10
0
10
20
−20
60 −20
0
10
20
−20 c)
−10
−10
10−2
d) −10
t = 2.0735
t = 2.7803
10−3 0
0
10
10
20
20
30
30
10−6
40
40
10−7
50
50
y
10−4
60 −20
−10
0 x
10
20
60 −20
10−5
10−8
−10
0 x
10
20
Figure 5.15: Reflection due to velocity aliasing on spatial resolution limit. Figure a) shows the first component, mostly unchanged in the coupling process. Figures b), c) and d) show the evolution of the outcoupled population in the second component, reflecting at the Nyquist limit and interfering with the initial beam. Figures a) and b) show the two states at the same time t. (Parameters: Ω=0.2, ∆=40, C = 200, 128 by 256 grid)
CHAPTER 5. RF OUTPUT COUPLER IN 2D
5.7
61
Absorber
Our simulation should model as closely as possible the real physical situation, where the outcoupled beam rapidly leaves the region of interaction. We have attempted to achieve this by using a thin absorbing shell at the edge of the computational grid, which ideally should remove all condensate population that enters the shell, and cause no reflections and no interactions with the incoming beam.
5.7.1
Spatially filtering the Wavefunction
The most straight forward implementation of an absorber just spatially filters the wavefunction in every simulation timestep, multiplying it with what essentially is a unity matrix that smoothly goes to zero at the borders. This is computationally very inexpensive and flexible. However, it is not really an absorber, but an ongoing reinitialization of the simulation. A smooth function that has proven to be very well suited is a sin1/8 (rise) and a cos1/8 (falloff) in a shell at the edge of the grid. An illustration for a 128 point grid and a shell width of 20% of the grid can be found in figure 5.16a. For larger grid sizes, the relative width of the shell can be reduced to 10% or even 5%, as long as about 20 to 30 grid points are used for the edge functions.
5.7.2
Complex Potential
An alternative approach to the absorber problem is to use a complex component added to the trapping potential in the edge regions, where absorption is required. The same considerations as above, concerning the width of the shell in the edge regions, apply in this situation. Considering only the trap part of the time dependent GP equation ∂ψ = ...V ψ..., (5.30) ∂t we see that for a complex potential V = −ik, equation (5.30) becomes i
∂ψ = · · · − kψ . . . ∂t
(5.31)
ψ(t) = ψ0 e−kt .
(5.32)
which has the solution
CHAPTER 5. RF OUTPUT COUPLER IN 2D
62
This describes amplitude, damping with absorption coefficient k. By adding a complex component to the potential in the grid edge regions, we design a spatially dependent absorption coefficient. To avoid reflections, the complex 1.2 1
factor
0.8
a)
0.6 0.4 0.2 0
−10
−8
−6
−4
−2
0
2
4
6
8
10
−6
−4
−2
0 grid position
2
4
6
8
10
complex potential cmpt.
0 −20
b)
−40 −60 −80 −100 −120
−10
−8
Figure 5.16: Profiles of matrices used in the two different absorber approaches. Figure a): shaping method with a sin1/8 rise at the left edge and a cos1/8 falloff at the right edge of the grid. Figure b): complex component of the trapping potential with a (− sin1/8 · exp) falloff to the right edge. Opposite edges are mirror images, and the functions in both approaches affect 20% of the grid on each edge in this illustration. potential component should increase as smoothly as possible from zero at the inner border of our absorbing “shell” and reach high magnitudes at the grid edge to prevent tunneling “transmissions” due to periodic boundary effects. We used an exponential function to obtain a steep increase and multiplied it with −(1 − cos)1/8 (and − sin1/8 at the opposite edge) to smooth the transition to zero at the inner edges. (Sine and cosine periods were adjusted to obtain an increase/decrease between 0 and 1 over the appropriate intervals.) This method is mathematically more elegant and justifies the notion of an absorber better than the previous approach. We carried out a comparison of the performances of the two different methods. An isotropic two dimensional condensate wavefunction of Gaussian ra-
CHAPTER 5. RF OUTPUT COUPLER IN 2D
63
dial profile, i.e. C = 0, was placed in the centre of a 1282 computational grid (spatial extents of the grid x, y = [−10, 10]), with an absorbing shell of 20% of the grid, and was given a centre of mass velocity by imposing a linear phase gradient on the wavefunction. The critical velocity, where grid aliasing reverses the cloud’s motion, is v = 20 at this grid size. On the velocity interval v = [6, 14], well under the aliasing limit, yet with reasonably fast cloud movement, we measured the residual overall population after the cloud has completely passed through the absorbing region. The results are plotted in figure 5.17. We found that the complex potential method performed better for higher velocities while the shaping method performs well for low velocities. For low velocities, i.e less than v = 10, mostly reflections lead to residual populations. For cloud velocities higher than v = 10, the main cause of residual population was tunneling “transmission” to the opposite side of the grid (due to periodic boundary effects). Note that the performance analysis discussed above represents a very specific case and that adjustments to the width of the absorbing shell, to the gridsize and to the functions used to model the absorption coefficients leave a wide margin for customizations. In general, however, the complex potential method appears to be superior to the shaping method for higher velocities, while the reverse is true for lower velocities. Since the outcoupled streams in our cw rf coupling simulations accelerated to near critical grid velocities while falling under gravity towards the edge of the grid, we have chosen to use the complex potential method in our simulations.
5.8
Validation of simulations
Since coding errors and conceptual mistakes are easily made, simulation code must be tested to make sure that simulations are physically valid. We can test simulation code on simple scenarios with predictable outcome. A very simple and predictable simulation is, for example, that of an atomic cloud undergoing trap oscillations. The cloud oscillates with the the trap frequency around the trap equilibrium and, because of the absence of external interactions, the total energy and the normalisation, the particle number, are conserved [32].
CHAPTER 5. RF OUTPUT COUPLER IN 2D
64
−3
10
−4
residual pupoluation
10
−5
10
−6
10
−7
10
−8
10
6
7
8
9
10 velocity
11
12
13
14
Figure 5.17: Residual relative population after a moving condensate cloud has completely passed through the absorbing or spatially filtering region at the grid edge. The solid line is the result from the shaping method and the dashed line is from the complex potential method. The diagrams in figure 5.18 illustrate the simulation of a cloud at initial position y=0, oscillating about the equilibrium position at y=10 (vertical offset due to gravity). Using Ψ for Ψ(x, y, t), the values monitored are the vertical centre of mass position Z∞ Z∞ |Ψ|2 y dx dy, hy(t)i =
(5.33)
−∞−∞
the kinetic energy, Ekin (t) = −
Z∞ Z∞
Ψ∗ ∇2 Ψ dx dy,
(5.34)
−∞−∞
and the potential energy Epot (t) =
Z∞ Z∞
−∞−∞
1 2 2 r |Ψ| dx dy. 4
(5.35)
CHAPTER 5. RF OUTPUT COUPLER IN 2D
65
The internal energy, or self energy, of the cloud is calculated as Eself =
Z∞ Z∞
1 C|Ψ|4 dx dy. 2
(5.36)
−∞−∞
a)
0
centre of mass position y
Self energy does not change over time in the simulation described above, since the cloud retains its shape undergoing oscillations in the trap.
5
15
20
b) energy [arbitr. units]
trap equilibrium
10
0
2
4
6
8
10
12
14
1
0
-1 0
2
4
6
8
10
12
14
t
Figure 5.18: (a) Single component condensate cloud in vertical oscillating around trap equilibrium at y=10. While the cloud oscillates with ωtrap =1, (b) kinetic energy Ekin (solid line) and potential energy Epot (dotted line) oscillate with ω=2. Total energy Etot (dashed line) remains constant. Self energy Eself is also constant (but not plotted). In general, conserved quantities are calculated and monitored as a routine operation since these are very sensitive indicators of simulation validity. Normalisation of the wavefunction, for example, is the single most important indicator for code errors and misconfigurations. In coupled two-state 2D simulations, normalisation over the two present components is defined as the spatial
CHAPTER 5. RF OUTPUT COUPLER IN 2D
66
integral over the sum of the two states Z∞ Z∞
−∞−∞
|Ψ1 |2 + |Ψ2 |2 dx dy ≡ 1,
(5.37)
representing bosonic particle conservation over all participating components. Generally, normalisation fluctuates slightly during the course of a simulation because of rounding errors. In valid simulations of 256 by 512 point grids, these fluctuations were within the order of 10−5 or 10−4 . For lower grid resolutions, these fluctuations were higher, starting to become intolerable for 64 2 grids. Although we used large grids for “production” simulations, quick runs on 1282 grids could be used as rough indicators or for rapid probing of parameter regimes. This posed another validity test: Changing the grid size, or, as another option, decreasing the time step size, should not change the fundamental physical outcome of any simulation run. To evaluate the magnitude of numerical errors, we can stop our simulations after a certain time and then reverse the time steps. Running the simulation backwards to t = 0 should finally reconstitute the initial state, since the GP equation is unitary. All deviations are then due to numerical errors. Output coupling and TOP simulations (presented in chapter 6) reconstituted the initial state with a difference between initial state Ψi and final state Ψf Z∞ Z∞
−∞ −∞
|Ψf |2 − |Ψi |2 dx dy
(5.38)
of less than 10−4 or 10−3 , depending on how far forward it was run initially and what grid and time step size was used. Time reversal of simulations which included damping at grid edges, vastly amplified numerical errors and can thus not be considered representative for the precision of the simulation code. For debugging purposes and in order to illustrate specific interpretations, it was very convenient to disable spatial diffusion by omitting the (∇2 Ψ)terms in the coupled GP equations, thus effectively making the Thomas-Fermi approximation which has been described earlier. For example, figure 5.10 on page 51, which illustrates the effect of finite pulse time broadening, was generated with diffusion disabled. With population densities staying in place under this provision, spatial effects in coupling can be investigated in greater
CHAPTER 5. RF OUTPUT COUPLER IN 2D
67
detail and evolution at individual spatial points can be checked against analytic solutions.
Chapter 6 TOP Trap Simulations One of the most commonly used magnetic traps in BEC is the Time-averaged Orbiting Potential (TOP) trap, and in particular it is the trap used in the Otago BEC experiment. This ingenious configuration is based on the simple quadrupole trap, consisting of two Helmholtz coils on the z axis, generating magnetic fields of opposite polarity with a magnetic field zero in the trap centre between the coils. The magnetic field increases linearly in all spatial directions and thus atoms in low field seeking atomic states are trapped. However, atoms undergo spin flips into untrapped states in the zero-field region and can thus escape from the simple quadrupole trap. Cornell et. al. proposed to introduce a magnetic bias field of constant magnitude in the xy plane. This bias field moves the magnetic field zero away from the trap centre. By rotating the bias field, which is generated by two sets of Helmholtz coils on the x and y axis respectively, the magnetic field zero can be moved on a circle in the xy plane around the centre of the trap. This rotation can be done so quickly that the atoms in the trap centre do not have time to re-centre at the shifted magnetic field zero and thus never enter regions of near-zero magnetic field magnitudes, where trap loss can occur. The physics and the numerical implementation of quadrupole and TOP traps will be described in detail in the following sections. We simulate two-state coupling in the TOP trap in two different parameter regimes. We describe and analyse the dynamic vortex behaviour which we discovered, and we propose a model describing the spatial coupling effects and the formation of vortices.
68
CHAPTER 6. TOP TRAP SIMULATIONS
6.1
69
The Quadrupole Trap z
B1
Pair of antiHelmholtz coils
y
x Cloud of trapped atoms around magnetic field zero point
B2
Figure 6.1: Schematic illustration of a simple quadrupole trap. The Helmholtz coils positioned on the z axis carry anti-polar currents and create anti-polar fields B1 and B2 . The cloud of trapped atoms rests near the magnetic field zero point at the origin. A quadrupole trap consists of two Helmholtz coils centred on the z axis with equal currents in opposite directions (see figure 6.1), creating a magnetic field zero in the centre between the coils and a linearly increasing magnetic field amplitude in all spatial directions. For small displacements from the centre of the trap, i.e. less than the radius of the Helmholz coils and less than the separation of the two, the magnetic field can be described as B = Bx x ˆ + By y ˆ + Bz ˆ z = (Bq0 x)ˆ x + (Bq0 y)ˆ y + (2Bq0 z)ˆ z
(6.1)
where Bq0 , the gradient of the magnetic field in radial or (x,y) directions, is constant, and x, y and z are displacements from the origin. Note that the field gradient in the axial direction is twice as large as in the radial, leading to a stronger trap confinement along the z axis (and ultimately to “pancake shaped” condensates).
CHAPTER 6. TOP TRAP SIMULATIONS
70
The potential energy of a magnetic dipole µ in a magnetic field B is U = −µ · B. Under the assumption that the atoms’ magnetic moments always align with the magnetic field lines and follow any small changes in B adiabatically, this simplifies to U = −µB. The trapping force is proportional to the local field gradient, and given by F = −∇U
(6.2)
Thus the magnitude of the trapping force e.g. in the z direction is ∂B ∂B ∂U =µ = g F MF µ B . (6.3) ∂z ∂z ∂z 87 Rb F =2 states, as used in the Otago experiment, have a Land´e-factor of gF = 21 , so that they are low field seeking for magnetic quantizations MF =[1,2], magnetically untrapped with MF =0 and high field seeking, and thus expelled from the trap, for magnetic quantizations of MF =[-2,-1]. Trapping low field seeking states works as long as the Larmor frequency µB ωLarmor = gF MF B (6.4) h ¯ of the atoms is larger than the rate of change of the magnetic field due to the atom’s movement in the inhomogeneous magnetic trapping field. This is required so that the rapid precession of the atom’s magnetic moment around the magnetic field lines allows the atom’s mean magnetic moment to follow adiabatically the changes in local field direction. At the trap centre, the magnetic field B and the Larmor frequency go to zero, allowing trapped atoms to undergo Majorana spin-flip transitions into untrapped states [24]. This process constitutes the trap “hole” mentioned before, the major drawback of the quadrupole trap. F =
6.2
Time-Orbiting Potential (TOP) Magnetic Trap
The “hole” in the bottom of the quadrupole trap can be plugged by rotating the homogeneous bias field of two pairs of Helmholtz coils in the x-y plane. This rotating field, Btop = Bbias [cos(ωtop t)ˆ x + sin(ωtop t)ˆ y] ,
(6.5)
CHAPTER 6. TOP TRAP SIMULATIONS
71
constitutes the so called “TOP” field. The bias field moves the magnetic field zero point away from the origin. Because of the geometry of the bias field coils centred on x and y axes, the resultant field is simply the original field displaced sideways by an amount Bbias in the xy plane. The rotation of the bias field Bbias causes the displaced quadrupole, and in particular the point of zero field, to rotate around the origin with an angular velocity ωT OP at a radius r¯0 =
p
x¯2 + y¯2 =
Bbias . Bq0
(6.6)
The B field zero traces out a trajectory which is commonly known as “the circle of death”. Actually, there is a small spatial region about the zero field point with relatively high probability of population loss into untrapped states, and this is more accurately called the “Doughnut of Doom” [28]. Neglecting influences of gravity, the trapped atomic cloud resides at the origin, well within radius r0 , as we will demonstrate in the following section. Adding the quadrupole field and the rotating bias, we get the following time dependent magnetic field, which we will occasionally refer to as the “rotating TOP field” [4]: B = (Bq0 x + Bbias cos(ωtop t))ˆ x + (Bq0 y + Bbias sin(ωtop t))ˆ y + (2Bq0 z)ˆ z
(6.7)
A schematic illustration of the TOP field is shown in figure 6.2. We show the trajectory of the magnetic field minimum around the origin, where the trapped cloud rests, and an instantaneous projection of the magnetic field magnitude in the xy plane. The height of the cone, which increases linearly in x and y directions from the minimum point, over the xy plane characterises the field magnitude. Note that this simplified figure is only illustrating the field magnitude in the xy plane, the projection we chose for our numerical simulations.
6.2.1
Calculating the Time Average
The reason why TOP traps are extremely useful and in widespread use in BEC is that the orbiting potential time-averages to a harmonic potential, as we now show. Starting with equation (6.7) and integrating over one TOP rotation
CHAPTER 6. TOP TRAP SIMULATIONS
72
z
Atom cloud in trap centre y
Bquad
x
path of magnetic field B=0 region "circle of death" in x-y plane
Figure 6.2: Schematic illustration of the field in the xy plane of a TOP trap. Instantaneous view of the displaced quadrupole field. (Topographical view of magnetic field magnitude in the xy plane in a 2D projection along the z-axis.) p θ = [0, 2π], with r = x2 + y 2 , we obtain Z 1 2π q 2 0 2 2 dθ r Bq + Bbias + 2Bq0 Bbias (x cos θ + y sin θ) + 4z 2 Bq0 2 Baverage = 2π 0 (6.8) This can be rearranged in terms of the radius of the circle of death (c.o.d.) r0 = Bbias /Bq0 , and we get Z 2π s Bq0 q 2 2r0 (x cos θ + y sin θ) dθ 1 + Baverage = r0 + r2 + 4z 2 . (6.9) 2π r02 + r2 0 We can make a binomial expansion of Eq. (6.8), which is valid for regions well within the c.o.d., which is where most atoms reside. The zeroth order yields a constant term and first order integrates to zero, and thus we need to include the second order term to get a result that varies with position. We obtain r2 r2 + 4z 2 0 Baverage ≈ Bq r0 1 + 1− 2 2r02 4r0 4 2 2 r + 4z 2 r2 r + 8z 0 − = B q r0 + + ... (6.10) 4r0 8r03
CHAPTER 6. TOP TRAP SIMULATIONS
73
When terms in r0−3 and smaller are neglected, we are left with a time-averaged harmonic potential in all three dimensions, where the trap “spring” constant in z direction is eight times larger than that in the xy plane. This leads to “pancake” shaped, or elliptic confining potentials, and condensates that are flattened in the z direction because of the stronger confinement, and circular in their xy plane projection. The constant offset field Bq0 r0 does not affect the trapping properties (since the trapping force is proportional to the field gradient), but remedies the zero field hole of the simple quadrupole arrangement. The usual form taken for the time averaged potential is (from 6.10) 0
V¯ = µBaverage
Bq2 1 2 (¯ r + 8¯ z 2 ). ≈ µBbias + µ Bbias 4
(6.11)
Here, we can see that the trap’s “spring constant” is given by the ratio of (B q0 )2 to Bbias . Thus, by increasing the magnetic field gradient Bq0 or by decreasing the bias field Bbias , we can increase trap stiffness. In dimensionless units, we get 0 µBbias µx20 Bq2 1 2 V¯ (r + 8z 2 ). (6.12) ≈ + V = h ¯ ωr h ¯ ωr h ¯ ωr Bbias 4 A very important parameter is the angular velocity ωT OP of the rotating TOP field. ωT OP must be chosen large enough so that the trapped atoms do not have enough time to recentre in the shifted field and small enough so that the atomic magnetic momenta will remain adiabatically aligned with the field direction. ωLarmor ωT OP ωtrap (6.13)
6.3
Numerical Model of a TOP Trap
Our interest in this chapter is to dynamically simulate condensate evolution in the presence of the rotating field in a TOP trap. We explicitly take into account the time evolution of the magnetic field and determine its effect on the condensate dynamics. However, we are limited by computational resources to consider only a two dimensional model of the trap. We have chosen to consider the behaviour in the z = 0 plane. In SI units, the nonlinear GP equations for a coupled two component condensate with trapped first and untrapped second state in a “dynamic” TOP
CHAPTER 6. TOP TRAP SIMULATIONS
74
trap are (symbols with bars are in SI units, and we use ωT ≡ ωT OP for clarity) ¯1 p h ¯2 2 ¯ ∂Ψ ¯ 1 (6.14) =− ∇ Ψ1 +µBq0 (¯ x − r¯0 cos(¯ ωT t¯))2 + (¯ y − r¯0 sin(¯ ωT t¯))2 Ψ i¯h ∂ t¯ 2m h ¯Ω ¯ ¯ 2 ¯ 2 ¯ +i Ψ 2 − N U0 |Ψ1 | + w|Ψ2 | Ψ1 2
¯2 h ¯Ω ¯ h ¯2 2 ¯ ∂Ψ ¯ 2 (6.15) ¯ 2 − N U 0 |Ψ ¯ 2 |2 + w|Ψ ¯ 1 |2 Ψ =− ∇ Ψ2 + i Ψ h∆Ψ 1 − i¯ ∂ t¯ 2m 2 In our dimensionless units, the coupled GP equations are p ∂Ψ1 = i∇2 Ψ1 − iκ (x − r0 cos(ωT t))2 + (y − r0 sin(ωT t))2 Ψ1 (6.16) ∂t Ω +i Ψ2 + iC |Ψ1 |2 + w|Ψ2 |2 Ψ1 2 i¯h
Ω ∂Ψ2 = i∇2 Ψ2 + i Ψ1 − i∆Ψ2 + iC |Ψ1 |2 + w|Ψ2 |2 Ψ2 (6.17) ∂t 2 where we define the dimensionless parameter r µBq0 µBq0 x0 h ¯ = , (6.18) κ≡ h ¯ ωr h ¯ ωr 2mωr and where ωT ≡ ωT OP has been scaled with the trap frequency ωr , and ωT = ω ¯T . ωr The rotating TOP field averages to a harmonic field (Eq. 6.12) with a 0
coefficient of
µx20 Bq2 , 4¯ hωr Bbias
so that we need to set 0
µx20 Bq2 = 1, h ¯ ωr Bbias
(6.19)
in order to obtain a coefficient of 14 for the time-average of this trap. We are motivated to obtain a time-averaged trapping potential of 1 U = r2 , (6.20) 4 by the fact that we would be able to re-use in our TOP simulations the initial eigenstates previously generated for our harmonic simulations in the previous chapters. We can calculate Bq0 from (6.19) and then insert it into κ (6.18). Using (6.6), this yields that κ needs to be equal to r, the dimensionless radius of the “circle of death”, r¯0 κ= = r. (6.21) x0 With this definition, Eq. (6.11) describes an offset harmonic potential with the desired coefficient 14 .
CHAPTER 6. TOP TRAP SIMULATIONS
6.4
75
Spatial coupling
In a TOP trap, the dynamically evolving trapping field leads to a behaviour in the presence of coupling fields that is qualitatively different from the behaviour expected in a harmonic trap.
y
Region of strongest coupling Coupling field tuned to resonance at the centre of the cloud
x
Magnetic equipotential lines of trapping TOP field "Circle of death", movement of TOP field zero point
Figure 6.3: Motion of the instantaneous coupling surface during the TOP rotation. (Illustration not to scale.) For Rabi frequencies Ω larger than the frequency ωT OP , the coupling interaction occurs on a smaller timescale than the TOP field rotation period TT OP , and thus the resonant instantaneous coupling surface is a surface of constant TOP field magnitude and not a surface of constant time averaged (harmonic) field magnitude, as would be the case in a genuine harmonic trap. Thus, an approximate picture for the coupling interaction is a “static” B field, centred on the instantaneous field zero, as illustrated in figure 6.3. For TOP rotation frequencies close to, or larger than the Rabi frequency, this interpretation becomes problematic, and in fact, the TOP trap becomes more similar to a simple harmonic trap. In all of our TOP trap simulations involving coupling fields, the centre of the condensate is tuned to resonance, i.e. the effective detuning at the cloud
CHAPTER 6. TOP TRAP SIMULATIONS
76
centre is always zero during the TOP field rotation. Note that the curved coupling surface, and also the region of maximum coupling marked in the above illustration, is symmetric only about a radial line towards the point of the magnetic field zero on the circle of death (i.e. about the x axis in the illustration). In our simulations, the curvature of the coupling region is not negligible, leading to asymmetries and to transfer of angular momenta due to the coupling process, as we will show later in this chapter. At this point it is also worth noting that the above discussion is only to give an interpretation of the situation observed. The full GP equation treatment in our simulations describes the physical system correctly in any case.
6.4.1
Coupling Bandwidth
Coupling bandwidths in TOP traps are very similar to the bandwidths in harmonic traps, particularly in terms of spectral broadening and power broadening. The discussion of these contributions in section 5.5 applies here without modifications. A new feature of TOP traps, however, is the modulation of effective detuning at a given spatial point caused by temporal changes of the magnetic trapping field. While the interacting rf field is tuned to resonance with the condensate cloud centre at all times, the resonant surface (see figure 6.3) rotates with an angular velocity of ωT OP , so that all positions in the condensate, with exception of the centre, experience an oscillating detuning. We note that a different approach to the interpretation of the effects of oscillating detuning has been made by the experimental BEC group at Otago University. In [27], the changing detuning is simply treated as another linewidth component Bq0 , (6.22) Γtransit = 2πrωT OP Γpulse where Γpulse is the linewidth due to finite coupling pulse time as defined in (5.21). For our purposes, however, it is more appropriate to interpret it as a modulated detuning, the effects of which we will investigate in detail in the following section.
CHAPTER 6. TOP TRAP SIMULATIONS
6.4.2
77
Localized evolution of population and phase
The TF approximation is valid for times t Ttrap , i.e. for timescales of a few TOP field rotations. (Where Ttrap is the harmonic trap oscillation period.) Within this timescale, diffusion caused by non-zero kinetic energy is negligible and grid points can be treated as if they were independent of each other. In other words, as long as we only look at a few TOP rotations, simulations in the TF limit do not differ significantly from simulations including the kinetic energy term. Under the Thomas-Fermi (TF) approximation of zero kinetic energy, the set of coupled GP equations simplifies to a set of linear differential equations, and the evolution of population and phase at any position in our two-dimensional computational grid can be described independently using the simple two-state model outlined in section 2.2. This enables a significant simplification of the full equations, for the purposes of interpreting the mechanisms at work. An important process in systems with time-dependent detunings is that of adiabatic inversion of population, which we will outline in the following paragraphs. Adiabatic inversion and adiabatic following In the geometric representation of the system state within the rotating wave approximation, as outlined in section 2.2, the state vector s = (u, v, w) (Eq. 2.31) precesses around the coupling “torque” vector Ω= (−Ω0 , 0, ∆) (Eq. 2.37). By using a time dependent detuning ∆(t) = −µt, it is possible to adiabatically move, and ultimately invert, the state vector from population component w = 1 to w = −1 [37]. Starting at t0 = −∞, we reach inversion w = 0 at time t = 0. This is illustrated in figure 6.4. Ω starts out almost parallel to the 3 axis, moves through (−Ω0 , 0, 0) at t = 0, and ends up parallel to the negative 3 axis, and the state vector follows it along, as it precesses around it. For an adiabatic following description to be valid, we require ∂∆(t) Ω, ∂t
(6.23)
i.e., the precession frequency must be larger than the chirp rate of the detuning sweep. As long as this requirement is met, the inversion process is qualitatively independent of the strength of the coupling interaction Ω0 .
CHAPTER 6. TOP TRAP SIMULATIONS
−µ t
78
s Ω −Ω0
1
2
3
Figure 6.4: Adiabatic population inversion by means of a detuning sweep from ∆ = ∞ to ∆ = −∞. While requirement (6.23) must be met for a truly adiabatic inversion process, inversions with residual population in the originating state can still be achieved when the chirp rate is larger than the precession rate by factors of up to around 20 (as illustrated in figure 6.5b). The magnitude of residual populations in these cases can depend on the specific phase relation between population and coupling interaction during the sweep, and the notion of “inversion” may not justified in all cases. If the chirp rate is even greater than that, the population does not pass the equilibrium point and may even partly couple back into the originating state in the second half of the chirp. Figure 6.5c illustrates a case where the sweep rate is 40 times as large as the precession frequency, well out of the adiabaticity condition. The population equilibrium point between the two states is reached only after resonance at t = 0.5, and we note that for an oscillating detuning in a TOP trap, when sweep rates in this order of magnitude are reached, population may not even cross the equlibrium point before it “swings” back due to a reversed detuning sweep. When we apply this model to the situation in the TOP trap, we also need to understand that in most positions close to the cloud centre we are not looking at a complete inversion process since the detuning only oscillates around values
CHAPTER 6. TOP TRAP SIMULATIONS
a)
79
1
0.5
0
b)
0
1
2
3
4
5
6
7
8
9
10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1
0.5
0
c)
1
0.5
0
t Figure 6.5: Adiabatic population inversion for Ω0 = 40 and sweep rates of (a) µ = −20, (b) µ = −800 and (c) µ = −1600. Resonance (∆ = 0) is reached at t = 5 in (a) and at t = 0.5 in (b) and (c). close to zero. Thus, the behaviour in this regime is dominantly the usual Rabi oscillations with a small modulation of the peaks due to the oscillating detuning. For points off the centre of the condensate at r = (x, y), the detuning oscillates with the TOP field rotation and the time-dependent effective detuning is p 0 2 2 ∆eff (t) = Bq r0 − (x − r0 cos(ωT OP t)) + (y − r0 sin(ωT OP t)) , (6.24)
where r0 is the radius of the c.o.d. For displacements r from the condensate centre, small compared to r0 , i.e. for all positions within the condensate, the effective detuning takes the simple form ∆eff (t) = −r0 r cos(ωT OP t)
(6.25)
CHAPTER 6. TOP TRAP SIMULATIONS
80
where Bq0 = r0 is the magnetic field gradient (Eq. 6.21). This yields an approximate maximum rate of change of the detuning of ∂ ∆eff (t) = r0 rωT OP . ∂t
(6.26)
For ωT OP = 60 and r0 = 20, we get a maximum rate of change in the detuning at the Thomas-Fermi “cloud border” (i.e. at r ≈ 7) of ∂t ∆ef f,max = 8400. Compared with the coupling frequency we used in the strong coupling case (where Ω0 ωT OP ), Ω0 = 350, there is a factor of 24 difference, placing the strong coupling regime into the region where “pseudo-adiabatic” inversion, or pseudo-adiabatic following, can be expected (figure 6.5b). In the weak coupling regime, i.e. where Ω ≈ ωT OP , this ratio is an order of magnitude larger, i.e. 140, so that only regions close to the condensate centre exhibit a behaviour, qualitatively comparable to the strong coupling regime with pseudoadiabatic following. At radii of larger than r ≈ 3, the coupling process cannot be as easily interpreted as in the strong coupling case, where pseudo-adiabatic following and Rabi cycling are dominant. We will return to the case, where pseudo-adiabatic following fails, in section 6.6.3, and we define the pseudoadiabatic following regime to be where 1 ∂ ∆eff (t) Ω0 . 100 ∂t
(6.27)
Specifically for positions r in a TOP trap, with Eq. (6.26) this becomes 1 rr0 ωT OP Ω0 . 100
(6.28)
Note that in the full GP equation treatment, the changing effective detuning is fully accounted for in form of a complete dynamic description of the TOP potential. The above discussions only concerns the interpretation of observations in our simulations, and the numerical solutions of the GP equations are accurate in all cases.
6.4.3
Force on a condensate in a TOP trap
To understand the coupling effects in a TOP trap, it is important to look at the evolution of phase and population on a small time scale, i.e. on the time scale characteristic of the coupling, which is less than one TOP field rotation. When
CHAPTER 6. TOP TRAP SIMULATIONS
81
integrated over one TOP period TT OP , the rotating linear field averages to a stationary harmonic field. The situation is more complicated on shorter time scales. For an eigenstate in a harmonic trap, phase remains flat and uniform, evolving slowly depending on the eigenvalue of the condensate. In a TOP trap, however, the cloud keeps developing a phase gradient towards the field zero and develops a velocity in response to the force, and a phase gradient. Since the field zero moves, this gradient also rotates. As we will show, it rotates with half of the TOP field rotation frequency and reaches its maximum after half a TOP field rotation. Since condensate velocity is proportional to the phase gradient (equation 4.17), the trap field rotation causes the condensate, which is starting out from a harmonic eigenstate at the trap centre, to gradually start “sloshing around” in our simulated trap as its centre of mass moves in the direction y of the time average of the phase gradient. In order to understand this behaviour, we consider the following simple treatment in terms of classical mechanics. The force F towards the magnetic field minimum of the rotating TOP field in the condensate in the trapped state can be expressed as F ∝ cos (ωT OP t) x + sin (ωT OP t) y.
(6.29)
The condensate velocity due to this force is proportional (and in dimensionless units equal to) v∝
Zt
Fdt0 =
0
1 ωT OP
sin (ωT OP t) x +
2 ωT OP
sin2
ω
T OP
2
t y.
(6.30)
This describes a simple oscillation along the x axis, and an oscillation with half the TOP frequency in the positive y direction. From (6.30), the evolution of the condensate displacement ∆r(t) can be calculated. We get ∆r =
Zt 0
vdt0 =
1 ωT2 OP
(1 − cos (ωT OP t)) x +
1 ωT2 OP
(ωT OP t − sin (ωT OP t)) y
(6.31) This calculation shows that the condensate centre of mass oscillates with the TOP frequency around a point slightly offset in positive x. The magnitude and
CHAPTER 6. TOP TRAP SIMULATIONS
82
offset of this oscillation decreases with the inverse square of the TOP frequency. A small drift occurs in positive y direction. The y offset increases with t/ω T OP due to this drift. This becomes more evident for longer simulation times, but it is negligible for simulations of only a few (i.e. 10) TOP field rotations.
6.5
Vortices
In order to understand the nature of superfluids like BEC, and the occurrence of vortices, we consider the change of population density in a volume V and derive the probability current density vector ∗ Z Z ∂Ψ ∂Ψ ∂ 2 ∗ |Ψ| dr = Ψ + Ψ dr ∂t V ∂t ∂t V Z Z i¯h ∗ ∗ = ∇·[Ψ (∇Ψ)−(∇Ψ )Ψ]dr = − ∇·jdr, (6.32) 2m V V where the second line is obtained by using the time-dependent GP equation h ¯2 2 ∂ 2 (6.33) ∇ + V + U0 |Ψ| Ψ i¯h Ψ = − ∂t 2m and its complex conjugate. For notational convenience we use Ψ for Ψ(r, t). Equation (6.32) allows us to interpret j as j(r, t) =
h ¯ Ψ∗ (∇Ψ)−(∇Ψ∗ )Ψ . 2mi
(6.34)
It is possible to write a condensate wavefunction Ψ = |Ψ0 |eiS where S is a real quantity, and can be interpreted as the phase, which leads to the following expression for the velocity [13]: v(r, t) = |Ψ0 |−2 j(r, t) =
h ¯ ∇S(r, t) m
(6.35)
From this result follows ∇ × v = 0.
(6.36)
This means that BECs are irrotational superfluids, since the curl of the condensate velocity vanishes everywhere where v is defined. The only way the superfluid can exhibit rotational flow is by having points where phase, and thus velocity, are undefined, i.e. a zero of the wavefunction. A structure in
CHAPTER 6. TOP TRAP SIMULATIONS
83
which the phase undergoes a positive or negative 2π winding is called a vortex (and in three spatial dimensions it becomes a vortex line). In stirred homogeneous superfluids like liquid helium, the flow is quantized, and one vortex forms for every unit of angular momentum [36]. Clouds of magnetically trapped alkali BECs are not of homogeneous density and thus there is no strict quantization of the superfluid flow. While a single centre ˆ z values occur while a single vortex vortex represents a unity flow, fractional L travels from the low population density cloud border to the high density centre [10].
6.6
Vortex dynamics in a TOP trap
The external coupling field connects trapped condensate state |1i and untrapped condensate state |2i. This creates a situation where population of state |1i with a specific phase gradient, which has developed due to the trapping potential, is temporarily transferred into and “stored” in state |2i. In this state, in the absence of a trapping field, which causes a radial phase gradient to develop on the time scale of Ttrap in the trapped state, the phase profile remains largely unchanged, evolving only slowly due to the slight energy shift provided by the coupling interaction as outlined in chapter 3. Half a Rabi cycle later, this “stored” population is transferred back from |2i and superposed onto the residual population in initial state |1i, which, being “exposed” to the rotating magnetic TOP field, has in the meantime evolved into a slightly different, rotated, phase profile. This is a situation in which phase singularities, “vortices”, can occur, as we will describe in detail in our model of vortex formation in section 6.6.2. In our simulations we can destinguish between two regimes with different vortex dynamics • Strong coupling, i.e. Ω ωT OP and fully in the pseudo-adiabatic following regime. • Weak coupling, i.e. Ω ≈ ωT OP and only partially (condensate centre) in the pseudo-adiabatic following regime.
CHAPTER 6. TOP TRAP SIMULATIONS
84
Note that these are the same two regimes that we distinguished before, in section 6.4.2. Since the medium of choice for the presentation of our simulations is a computer generated and played movie, it is difficult to illustrate some of the complex behaviour, which we see in many simulations, by means of limited sequences of static images. With this in mind, we hope the reader can excuse our somewhat excessive use of plots in the following sections.
6.6.1
Vortices in the strong coupling regime
In the strong coupling case, large Rabi frequencies and pseudo-adiabatic following within the entire condensate cloud lead to significant population transfer between the two states at all spatial positions on time scales much smaller than the TOP field rotation frequency. In this case, the picture of spatially dependent coupling, which we developed in section 6.4, is applicable and the TOP field rotates through only a small angle during one complete Rabi cycle. Also, the coupling bandwidth is large due to power broadening, and coupling occurs on a wide band across the whole condensate cloud. This band, however, has a Lorentzian profile across the cloud, so that the spatial coupling profile changes as the TOP field rotates around. This also leads to a slightly non-uniform density distribution of |1i compared with |2i, because of a changed coupling profile between two maxima of the Rabi cycle. We keep this in mind, because it is important for vortex formation by two-state coupling, as we will elaborate later. In this regime, strong coupling leads to effects that occur periodically with the Rabi frequency, making it the simplest to describe and interpret. In figure 6.6, we present a temporal sequence of plots of density and phase evolution of coupled condensates showing vortex formation. The left hand column shows the population density on a logarithmic scale, and the right hand column shows the condensate phase. The number N on the images indicates the number of TOP field rotations that have occured since the beginning of the simulation. Fig 6.7 and on page 89 and Fig. 6.9 show the evolution of total population in each state during this process. Two solid vertical lines in 6.7 mark the position of the beginning and the end of the sequence shown in 6.6, i.e. the beginning and end of the period between vortex creation and annihilation.
CHAPTER 6. TOP TRAP SIMULATIONS
85
The sequence in Fig. 6.6 displays a short fraction of a Rabi cycle, several Rabi cycles into the simulation. We see from Fig. 6.7, that at the first image of 6.6, the total population density in state |1i is decreasing towards zero, due to coupling into component |2i. At this time, coupling starts to transfer population back from |2i into state |1i, starting at the low density edge region of the condensate on the axis, on which we will see vortices travelling in a short time later. This coupling-in appears to start at the condensate cloud edge before the total population density minimum in the concerning state is reached. An explanation for the observed “inward sweep” of population minima and onsetting in-coupling is the delay between the start of in-coupling and the time when the in-coupling wavefunction actually starts to increase the population density. Initially, the in-coupling wavefunction may decrease the total amplitude at certain locations (like the diameter along which the vortices travel) in a complex addition to the residual wavefunction. This causes a transient decrease in local population density, which depends on the relative phases, and also the wavefunction amplitude profiles, of the two coupled components. Population from untrapped state |2i is now (in Fig.6.6b, c) coupled back, in the presence of a magnetic trapping profile that differs slightly from the one it was coupled out by, because the coupling region has rotated by a small angle due to the TOP field rotation. Thus, the phase profile of state |2i “lags” (has a gradient of a slightly different direction and magnitude) behind the one of the residual population in |1i, onto which it is now being added in the shape of the new (slightly rotated) spatial coupling region. Two vortices, minima in the density profile, appear at the positions along a cloud diameter, where the addition of complex wavefunctions yields a zero result. The two vortices are of opposite polarity and they travel inwards during only small fractions of a TOP cycle. In the condensate centre, they meet and annihilate shortly after the total population density goes through its minimum (in Fig. 6.7) at approximately N = 0.77. Summarizing, we can say that what we are observing is a radial inward sweep of superimposing one complex wavefunction profile onto another with a different phase and amplitude profile. This will be explained in more detail in our simplified model of vortex formation in section 6.6.2. In figure 6.7, the first vertical mark also represents the minimum in local
CHAPTER 6. TOP TRAP SIMULATIONS
86
population density at the cloud border, on the track, which the vortex travels along, and the minimum of the curve of local Rabi cycling Ωeff . It is important to note, that the vortices appear (with a short delay depending on the residual wavefunction amplitude) when peripheral coupling starts to couple population into the concerned condensate component. This happens alternatingly in both states |1i and |2i, with periodicity of the Rabi frequency Ω0 at the condensate centre, with half a Rabi cycle difference between the two states. At a Rabi frequency of Ω = 350 and a TOP rotation frequency of ωT OP = 60, this vortex scenario occurs almost 6 times in both states during one TOP rotation cycle. Vortices occur for the first time during the first minimum in state |2i (initially, at N = 0, the full population is in |1i). Surprisingly, no vortices appear during the first minimum of state |1i. This may be due to high initial bandwidth with very homogeneous coupling preventing a spatially different density profile to develop between the two coupled states immediately after turning on the coupling interaction at N = 0. The diameter, along which the vortices travel inwards at certain fractions of the Rabi cycle, rotates with the direction of the phase gradient of state |1i, i.e. with ωT OP /2, and there appears to be a 90◦ angle difference between the orientation of these diameters in states |1i and |2i. For simulations with finite C values (shown graphs represent simulations with C = 0 for reasons of simplicity), the same effects can be observed. However, the condensate clouds are slightly bigger. Also, the presence or absence of the kinetic energy term in the GP equation does not qualitatively change the simulations. This can be understood easily in terms of separated time scales. Diffusion occurs only on a much larger time scale than the coupling in the strong coupling regime. A situation which appears qualitatively different to the one described above occurs when the phase profiles of |1i and |2i, which are added in the coupling process, are nearly “flat”, i.e. do not have a gradient. A situation like this occurs during completion of (or in the early parts of) a full TOP rotation, i.e. around N = 1. This case is illustrated for state |2i in figures 6.8, using the same representation as previously in figures 6.6. Evolution of total population density is plotted in figure 6.9. The sequence shows that for tiny phase gradients across the condensate at the time of a specific Rabi cycle, the
CHAPTER 6. TOP TRAP SIMULATIONS
87
two vortices of opposite polarity “spread out” into oval or stretched density minima, connecting to form a complete “vortex-ring” within the condensate cloud. This “ring” contracts as the vortex pair travels inwards. The vortex pair annihilates in the same way as described before for the case of Fig. 6.6 and the total population density in Fig. 6.9 evolves as expected and similar to Fig. 6.7. The vortex ring travels into the condensate centre at the same velocity as any other typical vortex pair in the same coupling scenario. We conclude that the velocity, which the vortices travel at, depends on the Rabi frequency and not on the phase gradient across the condensate cloud, because the occurance of vortices is a coupling phenomenon as we will show in detail in the next section.
CHAPTER 6. TOP TRAP SIMULATIONS −10
π
−10
a1) t = 0.077493 N = 0.74
−2
10
10−3
88
a2) t = 0.077493 N = 0.74
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
π/2 10−4
10−5
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
0
5
10
π
−10
b2)
b1) t = 0.079063 N = 0.755
10−2
10−3
−π −5
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.079063 N = 0.755
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
c1)
−2
10
−π −5
c2)
t = 0.080634 N = 0.77
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.080634 N = 0.77
π/2 10−4
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
d1)
−2
10
−π −5
d2) t = 0.081158
t = 0.081158 N = 0.775
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
N = 0.775
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−π −5
0
5
10
Figure 6.6: Temporal evolution of condensate state |1i in the large coupling regime, during a small fraction of a Rabi cycle, after N = 0.74 TOP field rotations. Left hand column, population probability density. Right hand column, phase. The bar on the left gives the density scale and that on the right gives the phase scale. (Parameters: 2562 grid, Ω0 = 350, ωT OP = 60, r0 = 20, C = 0, coupling on resonance at cloud centre.)
CHAPTER 6. TOP TRAP SIMULATIONS
89
1
0.9
0.8
0.7
|Ψ|
2
0.6
0.5
0.4
0.3
0.2
0.1
0
0.68
0.7
0.72
0.74
0.76
N=f
TOP
0.78
t
0.8
0.82
0.84
0.86
Figure 6.7: Temporal evolution of the total population densities for figures 6.6. Solid lines represent state |1i and dashed lines represent state |2i. Vortices appear before the minimum in total density is reached (first solid vertical mark), and annihilate shortly after the minimum (second solid vertical mark).
CHAPTER 6. TOP TRAP SIMULATIONS −10
10−3
π
−10
a1) t = 0.10472
−2
10
90
a2)
N=1
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.10472 N = 1
π/2 10−4
10−5
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
b1) t = 0.10577
10−2
−π −5
b2)
N = 1.01
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.10577 N = 1.01
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
c1)
−2
10
−π −5
c2)
t = 0.10734 N = 1.025
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.10734 N = 1.025
π/2 10−4
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
d1)
−2
10
−π −5
d2)
t = 0.10838 N = 1.035
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.10838 N = 1.035
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−π −5
0
5
10
Figure 6.8: Temporal sequence of condensate state |2i in the large coupling regime, during a small fraction of a Rabi cycle after N = 1 TOP rotation. Left hand column, population probability density. Right hand column, phase. The bar on the left gives the density scale and that on the right gives the phase scale. (Parameters: 2562 grid, Ω0 = 350, ωT OP = 60, r0 = 20, C = 0, coupling on resonance at cloud centre.)
CHAPTER 6. TOP TRAP SIMULATIONS
91
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0.94
0.96
0.98
1
1.02
N=f
TOP
1.04
t
1.06
1.08
1.1
1.12
Figure 6.9: Temporal evolution of the total population densities for figures 6.8. Solid lines represent state |1i and dashed lines represent state |2i. Vortices appear before the minimum in total density is reached (first solid vertical mark), and annihilate shortly after the minimum (second solid vertical mark).
6.6.2
Model of vortex formation
In this section we explain the mechanism of vortex formation that was observed in our TOP simulations of the previous section. From the previous section we see that the basic condition for vortex formation is the following: We have an existing population in, for example, state |1i, with a given (near uniform) non-zero phase gradient. A second component in state |2i, with a different spatial density profile and a (near uniform) phase gradient, which points in a different direction than that of the existing population in |1i, is added to the existing population in state |1i. We can understand how this leads to vortices by means of the following simple model. Let us consider adding a wavefunction Ψ1 with a simple density hump (an inverted parabola) to a homogeneous (i.e. constant) wavefunction Ψ2 , where the two wavefunctions have some relative velocity. Without loss of generality we can assume that Ψ1 has a constant phase gradient along the x axis (fig
CHAPTER 6. TOP TRAP SIMULATIONS
|Ψ1|
b)
20
20
15
15
|Ψ|
|Ψ|
a)
92
10
|Ψ2 |
10
5
5
0 5
0 5 5 0
y
0 −5
−5
x
5 0
y
0 −5
−5
x
Figure 6.10: Simplified model of vortex formation. Plotted are probability amplitudes |Ψ| of (a) residual population with linear phase gradient along x axis, and (b) homogeneous probability amplitude of “in-coupled” population with “flat” phase. 6.10a) and the homogeneous wavefunction Ψ2 (fig 6.10b) has a constant phase (i.e. a zero gradient). The superposition of Ψs = Ψ1 + Ψ2 at any spatial point is obtained by adding the complex values of the fields at that point. In order to obtain a zero result, it is necessary that the magnitudes of the two added wavefunctions are equal, i.e. |Ψ1 | = |Ψ2 |, and in our model this occurs along a circle centred in the origin. For a zero result, furthermore it is also necessary that the phase of Ψ1 relative to Ψ2 be π, and we can easily see that due to its phase gradient, Ψ1 has the appropriate π phase difference along lines of constant x, of which there may be zero, one, two, three etc., intersecting the ring which is defining equal densities. The simplest case, which is adequate for our discussion, is when there is only one line of this appropriate phase difference intersecting the ring, and we see immediately that there will be two zeros of net amplitude, one at each end of this line. Furthermore, the superposition of Ψ1 and Ψ2 in the region about each zero leads to a 2π phase circulation as we now demonstrate. (Example with three vortex pairs shown in fig. 6.11.) We can visualize the situation with the help of Fig. 6.12, where we represent the superposition of the complex fields at a number of points on a small circle about the zero net amplitude of Ψ1 and Ψ2 . In each of Fig. 6.12 a, b, c, the solid arc represents (part of) the circle
CHAPTER 6. TOP TRAP SIMULATIONS
93
in the xy plane where |Ψ1 | = |Ψ2 |. The cross ‘x’ marks the point on the arc where the phases differ by π, and hence result in zero net field Ψs . On the dashed circle, we have shown at a number of points the complex value of of the field, by means of a vector with length representing amplitude and angle representing phase. Fig. 6.12a gives the field Ψ1 . We see that the length of the vectors increases as we move vertically up the page, i.e. as we move closer to the centre of the inverted parabola (at the origin of the xy coordinate system). The phase is constant along lines of constant x. In Fig. 6.12b, we represent the amplitude of the field Ψ2 . This is homogeneous, and all the vectors are of the same length and point in the same direction. Finally, in Fig. 6.12c, we represent the sum Ψs of Ψ1 and Ψ2 . It is clear that as we move around the small circle, the phase of Ψs goes through 2π, i.e. we have a 2π phase circulation and thus the structure is a vortex. For a case where wavefunction Ψ2 increases magnitude over time, i.e. due to population transfer by coupling, the circle |Ψ1 | = |Ψ2 | shrinks, i.e. moves inwards, as “higher” regions of the inverted parabola Ψ1 intersect with rising Ψ2 . This causes the vortices in Ψs to move inwards and ultimately to disappear (at the latest, in case the π phase difference occurs along x = 0) when Ψ1 reaches the level of Ψ2 at the top of the inverted parabola.
CHAPTER 6. TOP TRAP SIMULATIONS
−6
π
−6
b)
a)
10−2
94
−4
−4
−2
−2
0
0
2
2
4
4
−3
10
π/2
−5
10
10−6
y
10−4
0
−π/2
−7
10
10−8 6 −6
−4
−2
0
2
4
6
6 −6
−4
−2
x
0
2
4
6
−π
x
c) 7 log(|Ψ|2)
6 5 4 3 2 1 5 0 5
y
−5
0 −5
x
Figure 6.11: Superposition of the two different wavefunctions of figure 6.10. Three pairs of vortices of opposite signs appear. (a) shows the population density on a logarithmic scale, (b) shows the phase. The solid arrow points to the vortex studied closer in figure 6.12. In (c) a meshplot of the population density on a logarithmic scale is shown. The six vortices are the singularity points (→ −∞) on the logarithmic density scale.
CHAPTER 6. TOP TRAP SIMULATIONS To population maximum at parabola centre
−α
95
Ψ1
α
β
−β
a) −α
y
α
|Ψ1 |=|Ψ2 |
x
Ψ2
b)
Ψs
c)
Figure 6.12: Vortex formation in Ψs (c) by superposition of complex fields (a) Ψ1 , an inverted parabolic amplitude with constant phase gradient along x, and (b) Ψ2 , homogeneous population with zero phase, in a xy plane projection. Length of solid arrows symbolizes complex wavefunction amplitude, direction symbolizes relative phase angle. The vortex point is located in the centre of the dashed circle. (c) Ψs = Ψ1 + Ψ2 exhibits a 2π phase rotation around the vortex point.
CHAPTER 6. TOP TRAP SIMULATIONS
6.6.3
96
Vortices in the weak coupling regime
In the case of weak coupling interaction between the two states, i.e. with Rabi frequencies in the order of Ω ≈ ωT OP , we are in a regime, where pseudoadiabatic following (p.a.f.) does not occur through the entire condensate cloud. In this case, where p.a.f. is restricted to centre regions of the condensate, we have observed a more complex dynamical behaviour. In this regime, power broadening is not the main broadening mechanism of interest. Finite pulse time leads to significant initial broadening, and the width of the interacting region exhibits a significant and ongoing temporal decrease during a few TOP field rotations. This has previously been illustrated in figure 5.11. Furthermore, within time scales of full TOP rotations, the population in state |2i remains untrapped for longer than in the previous case, and thus the simple spatial coupling model of figure 6.3 can not be applied, since the time scales are not easily separated. For the sake of clarity, we have stopped our simulations after two TOP field rotations. All important effects can be illustrated within this period. In contrast to the strong coupling regime (where p.a.f. does occur through the entire condensate), simulations in the weak coupling regime show vortices and features of spatial coupling effects in both condensate components at the same time, in different spatial regions. We observe the creation of one vortexpair in each state per Rabi cycle. The vortex pair does not travel to the cloud centre, however, but rotates around the centre. One pair of vortices of opposite polarity orbits the cloud centre in each condensate state, positioned on rotating cloud diameters with a π/2 angle difference. Figures 6.13 (state |1i) and 6.14 (state |2i) show the dynamics of each of the states during the first two TOP field rotations at rotation numbers N = (a) 0.325, (b) 0.65, (c) 1, and (d) 2. To illustrate that vortices occur in both states at the same time in different spatial positions, we repeat the density profiles of figures 6.13 and 6.14 side by side in Fig. 6.15. To rule out effects of expansion in the untrapped state, we simulated a C = 0 state in the TF limit, i.e. simulation without the diffusion term in the GP equation. While the situation for states with finite C values is qualitatively the same as in the figures shown, the diffusion term adds a significant expansion (and thus a radial phase gradient) to state |2i. This gradient couples over into
CHAPTER 6. TOP TRAP SIMULATIONS
97
state |1i, so that eventually both states expand. In frames (a) of figures 6.13 and 6.14, we can see the first vortex pair (of opposite polarity) forming in |2i as population with a quickly increasing and rotating phase gradient is transferred into this state. This is essentially the same mechanism of vortex formation, by superposition of two different density profiles of different spatial phase, that leads to vortices in the large coupling regime. At N = 0.5, a situation between frames (a) and (b) which is not shown, the vortex pair in |2i lies on a vertical condensate diameter parallel to x = 0 and coupling has carved out a horizontal line of density minimum in state |1i while the TOP field has travelled half way around. This density minimum in |1i evolves into a vortex pair of opposite polarity, due to the same mechanism as before as population is transferred back into this state after half a TOP field rotation. After a full TOP field rotation and one complete Rabi cycle (at resonance in the cloud centre), instead of a full reversion resulting in complete population transfer back into state |1i, as one would expect in a simple harmonic trap, state |1i is left with a vortex pair on a vertical axis while a ring and a horizontal line of density minimum at the end of the completed Rabi cycle appear in state |2i; shown in frames (c). A second full TOP field rotation superposes a repetition of the above description with a simple evolution of the state as found at N = 1 (after one complete TOP rotation). New vortex pairs, created in the second TOP rotation period, appear closer to the cloud centre. In further TOP field rotations, the spatial separation of the vortex pairs, and thus the radius of vortex “rings”, decreases. The simulation state after two complete TOP field rotations is shown in frames (d). The dynamics in the weak coupling regime can be interpreted using the same approach as to the strong coupling regime, using the model of vortex formation for the strong coupling regime as presented in the previous section, and keeping in mind two key differences. • Coupling width decreases significantly during a few TOP field rotations, and the power broadened bandwidth is small compared with the strong coupling case, so that coupling the spatial coupling width with significant population transfer rapidly decreases to a slim rotating band across the condensate centre.
CHAPTER 6. TOP TRAP SIMULATIONS
98
• The simple spatial model in figure 6.3 and the pseudo-adiabatic following model (in section 6.4.2) fail for all but the immediate condensate centre regions. Hence, we only get complete population transfer by Rabi cycling in a small centre region of the condensate, the radius of which decreases with increasing interaction time. This explains the decreasing radius of the vortex occurances and ring features, which enclose the region where the simple models are still valid. The validity is limited to the centre regions of the condensate, where the temporal oscillation of the detuning is small enough (Eq. 6.28) compared with the Rabi frequency. In respect to the strong coupling regime, this TOP trap scenario can be considered a transition to a regime where ωT OP Ω. For increasingly large TOP rotation frequencies ωT OP , the behaviour of the TOP trap becomes more and more similar to that of a simple harmonic trap.
CHAPTER 6. TOP TRAP SIMULATIONS −10
π
−10
a1)
−2
10
10−3
99
a2)
t = 0.051051 N = 0.325
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.051051 N = 0.325
π/2 10−4
10−5
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
0
5
10
π
−10
b1)
10−2
10−3
−π −5
t = 0.1021
b2)
N = 0.65
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.1021
N = 0.65
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
c1)
−2
10
−π −5
t = 0.15708
c2)
N=1
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.15708
N=1
π/2 10−4
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
d1)
−2
10
−π −5
t = 0.31416
d2)
N=2
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.31416 N = 2
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−π −5
0
5
10
Figure 6.13: Temporal sequence of dynamics in condensate state |1i in the weak coupling regime during two complete TOP field rotations. Left hand column, population probability density. Right hand column, phase. The bar on the left gives the density scale and that on the right gives the phase scale. (Parameters: Ω0 = 40, ωT OP = 40, r0 = 20, C = 0, 2562 grid.)
CHAPTER 6. TOP TRAP SIMULATIONS −10
10−3
π
−10
a1)
−2
10
100
t = 0.051051
a2)
N = 0.325
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.051051
N = 0.325
π/2 10−4
10−5
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
b1)
10−2
−π −5
b2)
t = 0.1021 N = 0.65
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.1021
N = 0.65
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
10−8
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
c1)
−2
10
−π −5
t = 0.15708
c2)
N=1
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.15708
N=1
π/2 10−4
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−10
10−3
0
5
10
π
−10
d1)
−2
10
−π −5
t = 0.31416
d2)
N=2
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.31416
N=2
π/2 −4
10
−5
10
−6
10
0
−π/2
10−7
−8
10
8
8
10 −10
10 −10
−5
0
5
10
−π −5
0
5
10
Figure 6.14: Temporal sequence of dynamics in condensate state |2i in the weak coupling regime during two complete TOP field rotations. Left hand column, population probability density. Right hand column, phase. The bar on the left gives the density scale and that on the right gives the phase scale. (Parameters: Ω0 = 40, ωT OP = 40, r0 = 20, C = 0, 2562 grid.)
CHAPTER 6. TOP TRAP SIMULATIONS −10
101
−10
a1)
a1)
t = 0.051051 N = 0.325
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
8
8
t = 0.051051
N = 0.325
−2
10
10−3
10−4
10−5
−6
10
10−7
10 −10
−5
0
5
10
−10
10−8
10 −10
−5
0
5
10
−10
b1)
t = 0.1021
b1)
N = 0.65
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.1021 N = 0.65
10−2
10−3
10−4
−5
10
10−6
10−7
8
8
10 −10
10 −10
−5
0
5
10
−10
10−8
−5
0
5
10
−10
c1)
t = 0.15708
c1)
N=1
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.15708
N=1
−2
10
10−3
10−4
−5
10
−6
10
10−7
8
8
10 −10
10 −10
−5
0
5
10
−10
−8
10
−5
0
5
10
−10
d1)
t = 0.31416
d1)
N=2
−8
−8
−6
−6
−4
−4
−2
−2
0
0
2
2
4
4
6
6
t = 0.31416
N=2
−2
10
10−3
−4
10
−5
10
−6
10
10−7
8
8
10 −10
10 −10
−5
0
5
10
−8
10
−5
0
5
10
Figure 6.15: Rearrangement of figures 6.13 and 6.14 for a more convenient comparison. A temporal sequence of dynamics in states |1i (left column) and |2i (right column) in the weak coupling regime. Plotted are population densities on a logarithmic scale.
CHAPTER 6. TOP TRAP SIMULATIONS
6.6.4
Angular momentum
0.06
0.04
b)
a) 0.02
0.02
0
Lz
0.04
0
−0.02
−0.02
−0.04
−0.04
−0.06
−0.06
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
−0.08
−4
1
102
0
4
0.4
0.6
0.8
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
1
1.2
1.4
1.6
1.8
2
x 10
c)
0.8
0.2 −4
x 10
d) 3
0.6
0.4 2
Lz
0.2
0
1
−0.2 0 −0.4
−0.6 −1 −0.8
−1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
−2
0
N
N
Figure 6.16: Temporal evolution of angular momenta in simulation of figures 6.13 and 6.14. For figure 6.13: (a) state |1i, (b) state |2i. For figure 6.6: (c) state |1i and (d) state |2i. The solid lines represent total angular momentum ˆ z i, dashed lines represent angular momentum of the cloud centre of mass hL ˆ z,origin i, and dot-dashed lines represent the around the coordinate origin hL ˆ z,com i. difference, “residual” angular momentum about cloud centre of mass hL A mirror-symmetry, in all of the TOP trap simulations shown, about an axis through the condensate components can be found in all images of the temporal sequences. We have symmetry along the x axis at integer number of TOP field rotations, and vortices only appear in form of pairs of opposite polarity. Looking at the development of angular momentum of the two components undergoing coupling, we find, however, that this symmetry is not perfect. We calculated global angular momentum as ˆ z i = hxˆ hL py − y pˆx i
(6.37)
CHAPTER 6. TOP TRAP SIMULATIONS
103
and the angular momentum arising from centre of mass condensate motion around the trap origin as hLz,origin i = hxihˆ py i − hyihˆ px i.
(6.38)
While total angular momentum and angular momentum around the coordinate origin were expected to oscillate, because of the coupling process and the slight “sloshing” of the condensate in the TOP trap as described in section 6.4.3, the compensated centre of mass angular momentum of the cloud, Lz,com = hLz i − hLz,origin i,
(6.39)
was expected to vanish because of said symmetry. This was not the case, however, and a typical result for the temporal evolution of the angular momenta is shown in Fig. 6.16. These results can be explained by the asymmetry of the coupling region, illustrated in Fig. 6.3. This asymmetry, and the time-dependent width of the coupling region, account for the residual angular momenta about the cloud centre of mass. Figures 6.16 illustrate the differences between the two coupling regimes, which we have discussed. In the strong coupling case, i.e. in frames (c) and (d), the rapid periodicity with Ω0 almost entirely reverses population transfers in every Rabi cycle, while the spatial coupling width decreases on time scales of ωT OP , due to increasing interaction time. This explains the fact that all angular momentum components vanish after complete TOP rotations at N = 1 and N = 2. In the weak coupling regime, i.e. Fig. 6.16, frames (a) and (b), population transfers occur on time scales of the TOP field rotation, so that the spatial bandwidth and the coupling region change significantly before transfers can be reversed by Rabi cycling. This leaves behind asymmetric population density profiles, which account for the angular momenta. Note that the angular momenta in the strong and weak coupling regimes differ by two orders of magnitude.
CHAPTER 6. TOP TRAP SIMULATIONS
6.7
104
Discussion
Our models for vortex formation in a spatially dependent coupling has enabled us to explain the formation of vortices, which we have seen in our TOP trap simulations in the strong coupling regime. This regime is also used by the experimental BEC group, as discussed in Appendix B. In the weak coupling regime, we have also seen vortex formation, but this is more difficult to interpret and it remains a problem of interest. However, we understand the ways in which our models fail, and we can interpret the weak coupling regime as a transition from the strong coupling regime to a behaviour more similar to that of a static and purely harmonic trap, as we approach Ω0 ωT OP . Investigation of TOP trap dynamics showed that the TOP trap is far more complex than expected and not very well understood, although it is widely used by experimental groups. Problems with TOP traps, like unexpected heating as described in [27], may be attributable to false assumptions of similarity with harmonic traps. Within the TOP trap condition (6.13) ωLarmor ωT OP ωtrap , these problems can possibly be avoided or mitigated by moving the experimental setup to smaller coupling interactions and larger TOP frequencies.
Chapter 7 Conclusion The Gross-Pitaevskii (GP) equation provides a relatively simple but surprisingly accurate description of dilute Bose-Einstein condensates in the ultra-low temperature regime. In this thesis we have applied the GP equation in numerical simulations in two spatial dimensions to two qualitatively different dynamic scenarios. In both situations we simulate a two-component condensate system, corresponding to two different hyperfine states for the condensate atoms, which are coupled by a radio frequency field. In the first scenario, we have simulated the dynamic behaviour of a simple rf output coupler for a BEC in a harmonic trap. The two different condensate states are assumed to possess different trapping properties, with one hyperfine state being trapped by a harmonic magnetic potential, while the second feels no trapping potential and escapes from the trapping region, accelerating under the influence of gravity. We have studied output coupling for different detunings and different coupling strengths, and we have shown that the combined effect of the static magnetic field and the rf coupling field leads to spatially localized output coupling. We have also shown that a regime exists, in which a more or less continuous stream of matter is released from a trap, and we have looked at the effects of bandwidth. We have modelled noisy coupling fields, implemented using stochastic phase jumps of the coupling field, and we showed that this destroys the coherence properties of outcoupled beams. The introduction of gravity into the system of coupled GP equations gave rise to certain technical problems in the numerical implementation of our sim105
CHAPTER 7. CONCLUSION
106
ulations, including the limitation on length of simulations due to population reaching the grid edge, falling under the influence of gravity. To overcome this limitation, we investigated two different methods of “absorbing” and removing population density from the system at the grid edges without interfering with the simulation core. We found that both of the two approaches have advantages in different parameter regimes. The second major part of this work was the implementation of a dynamical model of the time-orbiting potential (TOP) magnetic trap. Since computational resources limited us to two spatial dimensions, we modelled this in the plane of the TOP field orbit, which is perpendicular to gravitational acceleration. In both the weak and, in particular, in the strong coupling regime, which is currently typical for the BEC experiment at Otago, surprisingly, we discovered a rich dynamic vortex behaviour, when a trapped and an untrapped atomic hyperfine state are rf-coupled in a TOP trap. We developed a model, explaining vortex formation in TOP trap two-component coupling in terms of two superposed density profiles of different density and phase. In the strong coupling regime, where pseudo-adiabatical following (as described in section 6.4.2) occurs through the entire condensate cloud, the interpretation of vortex formation was possible applying our model. Our exact model fails in the weak coupling regime, where pseudo-adiabatic following does not occur through the entire cloud. A qualitative interpretation for the observed spatial features within this regime could be given nevertheless, and the weak coupling regime can be interpreted as a transition from the distinct spatial coupling behaviour, which is characteristic for the TOP trap in the strong coupling/ pseudo-adiabatic following regime, to the weak coupling regime, where the TOP trap (for increasing TOP rotation frequency ωT OP Ω0 ) becomes more and more like a simple harmonic trap. Future work In TOP simulations within the strong coupling regime, which can be interpreted using our model of vortex formation, we found short lived vortex pairs, spread out to rings of near-zero density in situations, when the phase of the two coupled condensate components happened to be near-homogeneous, i.e.
CHAPTER 7. CONCLUSION
107
“flat”. This happens after integer numbers of TOP field rotations as outlined in section 6.4.3. It has been suggested that these rings are solitonic in nature. However, we did not have sufficient time to pursue a profound investigation of this matter. At this time, we leave the question unanswered, and we encourage future work in this direction. Future work can also go into further investigation of TOP trap dynamics, and in particular, into a more complete model of vortex formation including the weak coupling regime. A matter of interest is the extension of this work to three spatial dimensions, where vortex points of our xy plane projection become vortex lines. A two dimensional treatment of the situation appears to be justified, since the third spatial dimension may be thought of as “frozen” out by a strong confinement along the z axis of the TOP trap, but a three dimensional treatment is of great interest and will likely yield further insight into the problem. Furthermore, the immediate effects of vortex formation and dissipation on the condensate states could be investigated within a finite temperature model in order to understand “heating” effects, which have been described in [27]. Dissipation of vortices created and annihilated by the mechanisms described in this work may give an explanation for these observations of heating, since heating and vortex annihilation seem to be closely related [23].
This thesis is also available electronically from the author’s internet pages at http://hubble.physik.uni-konstanz.de/jkrueger/ or http://max.yi.org/jkrueger/ The author’s email address is
[email protected] (July 2000.)
Appendix A Adiabatic Elimination Many dynamic physical systems depend on a number of variables that evolve differently in time. In some situations, and given suitable parameter regimes, it is possible and justified to make simplifications by approximating their behaviour. For example, a system may respond very slowly to interacting fields, so that a rapidly revolving field has little or no effect compared to a slowly varying or constant interaction. In the following paragraphs, application of the so called Adiabatic Elimination will be demonstrated on the example of phase evolution in a coupled two-state system [5][6].
A.1
Two-State System in Rotating Wave Approximation
The general state of a two-state system is a superposition of its eigenstates {|0i, |1i}: |Ψ(t)i = a(t)|0i + b(t)|1i (A.1) Assume a classical magnetic field B(t) = B0 cos(ωt)
(A.2)
causing a magnetic dipole interaction ˆ = m · B0 . h ¯Ω 108
(A.3)
APPENDIX A. ADIABATIC ELIMINATION
109
In a matrix representation |0i =
"
1 0
#
|1i =
"
0 1
#
(A.4)
(A.1) satisfies Schr¨odinger equation i¯h
∂ |Ψ(t)i = H|Ψ(t)i ∂t
(A.5)
with a Hamiltonian matrix H =h ¯
"
0 Ω cos(ωt) ∗ Ω cos(ωt) ω1
#
.
(A.6)
Now this is simplified by moving into a frame, rotating with the frequency of the interacting field. Rotating |Ψ0 (t)i = Γ|Ψ(t)i
(A.7)
leads to the rotated Schr¨odinger equation ∂ 0 ∂Γ ¯ 0 |Ψ0 (t)i i¯h |Ψ (t)i = i¯h + ΓH Γ−1 |Ψ0 (t)i = H ∂t ∂t with ¯0 = h H ¯
"
0 Ω cos(ωt)e−iωt Ω∗ cos(ωt)eiωt −∆
#
,
(A.8)
(A.9)
where ∆ = ω1 − ω. The time varying terms can now be expanded Ω cos(ωt)e−iωt =
Ω 1 + e−2iωt , 2
(A.10)
yielding a constant and a rapidly oscillating term. In the Rotating Wave Approximation this term is ignored since its average value vanishes. Thus, in the rotating frame, the evolution equation simplifies to # " Ω 0 ∂ 0 2 ¯ 0 |Ψ0 (t)i. i¯h |Ψ (t)i = h |Ψ0 (t)i = H (A.11) ¯ Ω∗ ∂t −∆ 2
APPENDIX A. ADIABATIC ELIMINATION
A.2
110
Adiabatic Elimination
For high detunings |∆| |Ω| and for a situation where the entire population is initially in the ground state, a(0) = 1, b(0) = 0, we can now adiabatically eliminate terms from the coupled Sch¨odinger equations (they represent the system in the rotating frame under the rotating wave approximation) Ω ∂a(t) =h ¯ b(t) ∂t 2 ∗ ∂b(t) Ω i¯h =h ¯ a(t) − ∆b(t) . ∂t 2 Now b(t) can be formally integrated: i¯h
Ω∗ ∂b(t) = i∆b(t) − i a(t) ∂t 2 ∗ Ω ∂ (b(t)e−i∆t ) = −i a(t)e−i∆t ∂t 2 Z t ∗ Ω a(t)ei∆s ds b(t) = −i ei∆t 2 0 Equation A.16 can be integrated by parts and one obtains ! t Z t Ω i∆t ∂a(s) e−i∆s e−i∆s b(t) = i e a(s) + ds 2 i∆ 0 ds i∆ 0
(A.12) (A.13)
(A.14) (A.15) (A.16)
(A.17)
The second term in the equation above can be neglected since the rate of change in the lower state population almost vanishes under our set of assumptions. One retains the following expression: Ω Ω b(t) = [a(t) − a(0)ei∆t ] ≈ a(t) (A.18) 2∆ 2∆ The last step is possible in regimes with high detunings, since the exponential term represents another fast oscillation that averages to zero. Now, putting this result back into A.12, one obtains Ω Ω2 a(t) ∂a(t) = −i b(t) = −i ∂t 2 4∆ This represents the phase evolution of
(A.19)
Ω2 t (A.20) 4∆ for the lower state in the rotating frame. Influences of the upper state have been successfully eliminated from the evolution equation for the lower state and thus the two-state problem has been approximated by a single state problem. Φ(t) = −
APPENDIX A. ADIABATIC ELIMINATION
A.3
111
Accuracy of the approximation
In order to understand how precise the approximation made in the previous section is, the magnitude of next order correction terms is required. Equation A.11 has eigenvectors (dressed states, denoted as {|−i, |+i}) and eigenvalues as follows: H0 |±i = h ¯ ω± |±i (A.21) eigenvalues are given by ∆ 1p 2 ± ∆ + |Ω|2 (A.22) 2 2 # " # " − sin( 2θ ) cos( 2θ ) and e− = in the {|0i, |1i} = sin( 2θ ) cos( 2θ ) ω± =
with eigenvectors e+
basis. Here, tan(θ) = rotation. We have
|Ω| ∆
represents the angle of the diagonalising matrix
|Ψ(t)i = a(t)|0i + b(t)|1i = a0 (t)|−i + b0 (t)|+i
(A.23)
Vectorising the state coefficients in the original base as ~a and in the eigenvector base as ~a0 , the basetransformation ~a0 = L~a is achieved by the rotation matrix # " cos( 2θ ) sin( 2θ ) . (A.24) L= − sin( 2θ ) cos( 2θ ) To obtain an expression for the evolution of the original basis, one transforms into the diagonalised basis, where evolution is very simple, decribed by matrix # " eiω− t 0 , (A.25) U (t) = 0 e−iω+ t and after that, one transforms back into the original basis. So ~a(t) = LT U (t)L~a(0).
(A.26)
needs to be calculated within the parameters applied before (for high de|Ω| tunings). the approximations cos( 2θ ) ≈ 1, sin( 2θ ) ≈ 2∆ and ω± ≈ Using |Ω|2 ∆ ∆ ± 2 1 + 2∆2 the first order approximation obtained is 2 2 |Ω|2 2 |Ω|2 |Ω|2 |Ω|2 i(∆+ |Ω| )t |Ω| −i 4∆ t −i t i(∆+ )t 4∆ 4∆ 4∆ e + 4∆2 e = e + 4∆2 e ~a(t) = Ω 2 2 |Ω|2 |Ω|2 |Ω| |Ω| i(∆+ 4∆ )t −i 4∆ t Ω −i 4∆ t Ω i(∆+ 4∆ )t e − e − e e 2∆ 2∆ 2∆ (A.27)
APPENDIX A. ADIABATIC ELIMINATION
112
Ignoring the fast rotating parts that average to zero, one gets the same results as before and it can be seen that the correction term for the ground state Ω , truely small. coefficient a(t) is second order in ∆
Appendix B Simulation and Experiment Since we have done the numerical work in this thesis mostly in dimensionless computational units, it is not always easy to relate our results to the experimental “reality”. In this appendix, we are going to outline the most important specifications of the TOP trap in the University of Otago BEC experiment. Radial and axial trap frequencies are ωtrap,r = 2π · 90.164Hz ωtrap,z = 2π · 255.02Hz
(B.1)
The TOP frequency is ωT OP = 2π · 7kHz, so that the ratio ωT OP /ωtrap is 77.6. Thus, our simulations with ratios of 40 and 60 are in the same order of magnitude. The orbiting TOP bias field has a magnitude of BT OP = 4G and a gradient of Bq0 = 200G/cm, leading to a radius of the circle of death of r0 =
BT OP = 0.2mm Bq0
(B.2)
With trapped condensate clouds of a radial diameter of 20µm, the ratio between radius of the circle of death and radius of the condensate cloud is r = 0.05. r0
(B.3)
This is quite different from the factor of 0.35 in our simulations, which was chosen as a compromise in a trade-off between simulation speed and similarity with the experiment. Another major difference between theory and experiment is the C value, the magnitude of the nonlinearity. In the experiment, this value 113
APPENDIX B. SIMULATION AND EXPERIMENT
114
is orders of magnitude larger than in our simulations, where we chose C = 200 (and C = 0 for analytical purposes in some simulations). Higher C values, due to large atom numbers and three spatial dimensions, lead to larger condensate clouds and a magnetic field gradient across the condensate, which is much larger compared with the cloud diameter than in the simulation. This causes spatially more selective coupling. Note that the differences between experiment and simulation described above do not change the fundamental processes in the TOP trap, so that our results should be applicable to the experimental situation. Observation of the vortices created by the coupling process may be very difficult, however, for two reasons. First, vortices occur on time scales short compared with the TOP field rotation frequency, so that an excellent control of the fields is necessary in order to observe “static” snapshot pictures of the dynamic process in time-of-flight imaging methods. This also applies, to a lesser extent, to experiments in the weak coupling regime. Second, vortices occur in the weakly confined xy plane (presumeably in shape of vortex lines in three spatial dimensions), which is the plane with the smallest expansion rate during time-of-flight. These unfavourable circumstances will make experimental observations of the effects described in this work difficult.
Bibliography [1] L. Allen and J. H. Eberly. Optical resonance and two-level atoms. John Wiley & Sons, New York, 1975. [2] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell. Observation of Bose-Einstein condensation in a dilute atomic vapor. Science, 269:198, 1995. [3] R. J. Ballagh, K. Burnett, and T. F. Scott. Theory of an Output Coupler for Bose-Einstein Condensed Atoms. Phys. Rev. Lett., 78:1607, 1997. [4] Peter Bance. Evaporative Cooling of Cs in a TOP trap: Prospects for BEC. PhD thesis, University of Oxford, 1998. [5] P. B. Blakie and R. J. Ballagh. Meanfield treatment of Bragg scattering from a Bose-Einstein condensate. cond-mat/9912422, 1999. [6] P. Blair Blakie. Adiabatic elimination in a 2-state system. Seminar 11/99, unpublished. [7] P. Blair Blakie. Eigenstates and Excitations of Bose-Einstein Condensates. BSc (Hons) thesis, University of Otago, 1997. [8] Immanuel Bloch, Theodor W. H¨ansch, and Tillman Esslinger. Atom Laser with cw Output Coupler. Phys. Rev. Lett., 82:3008, 1999. [9] B. H. Bransden and C. J. Joachain. Introduction to Quantum Mechanics. Longman Scientific & Technical/ John Wiley, New York, 1989. [10] B. M. Caradoc-Davies, R. J. Ballagh, and K. Burnett. Coherent Dynamics of Vortex Formation in Trapped Bose-Einstein Condensates. Phys. Rev. Lett., 83:895, 1999. 115
BIBLIOGRAPHY
116
[11] C. Cohen-Tannoudji, B. Diu, and F. Laloe. Quantum Mechanics. John Wiley & Sons, Paris, 1977. [12] C. Cohen-Tannoudji, J. Dupont-Roc, and Grynberg G. Atom-Photon Interactions. John Wiley & Sons, New York, 1988. [13] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari. Theory of Bose-Einstein condensation in trapped gases. Rev. Mod. Phys., 71:463, 1999. [14] A. Eschmann, R. J. Ballagh, and B. M. Caradoc-Davies. Formation of Ramsey fringes in double Bose-Einstein condensates. J. Opt. B: Quantum Semiclass. Opt., 1:383–386, 1999. [15] R. P. Feynman, F. L. Vernon, and R. W. Hellwarth. J. Appl. Phys., 28:49, 1957. [16] Fox, R. J. Ballagh, and A. Wall. In preparation. [17] C. W. Gardiner. Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences. Springer-Verlag, Berlin, 1983. [18] C. W. Gardiner. Particle-number-conserving Bogoliubov method which demonstrates the validity of the time-dependent Gross-Pitaevskii equation for a highly condensed Bose gas. Phys. Rev. A, 56:1414, 1997. [19] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell. Dynamics of Component Separation in a Binary Mixture of BoseEinstein Condensates. Phys. Rev. Lett., 81:1539, 1998. [20] D. S. Hall, M. R. Matthews, C. E. Wieman, and E. A. Cornell. Measurements of Relative Phase in Two-Component Bose-Einstein condensates. Phys. Rev. Lett., 81:1543, 1998. [21] K. Helmerson, W. D. Phillips, K. Burnett, and D. Hutchinson. Atom lasers. Physics World, 12(8), 1999. [22] B. Jackson, J. F. McCann, and C. S. Adams. Output coupling and flow of a dilute Bose-Einstein condensate. J. Phys. B, 31:4489, 1998.
BIBLIOGRAPHY
117
[23] B. Jackson, J. F. McCann, and C. S. Adams. Dissipation and Vortex Creation in Bose-Einstein Condensed Gases. cond-mat/9912365v2, 2000. [24] W. Ketterle and N. J. van Druten. Evaporative cooling of trapped atoms. Advances in Atomic, Molecular and Optical Physics, 37:181, 1996. [25] P. L. Knight and L. Allen. Concepts of Quantum Optics. Pergamon Press, Oxford, 1983. [26] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard. Vortex formation in a stirred Bose-Einstein condensate. cond-mat/9912015, 1999. [27] J. L. Martin, C. R. McKenzie, N. R. Thomas, D. M. Warrington, and A. C. Wilson. Production of two simultaneously trapped Bose-Einstein condensates by RF coupling. cond-mat/9912045, 1999. [28] Jocelyn L. Martin. Magnetic Trapping and Cooling in Cs. PhD thesis, University of Oxford, 1997. [29] Karl-Peter Marzlin and Weiping Zhang. Quantized circular motion of a trapped Bose-Einstein condensate: Coherent rotation and vortices. Phys. Rev. A, 57(6):4761–4769, 1996. [30] M.-O. Mewes, M. R. Andrews, D. M. Kurn, D. S. Furfee, C. G. Townsend, and W. Ketterle. Output Coupler for Bose-Einstein Condensed Atoms. Phys. Rev. Lett., 78(4):582, 1997. [31] S. A. Morgan. A Gapless Theory of Bose-Einstein Condensation in Dilute Gases at Finite Temperatures. PhD thesis, University of Oxford, 1999. article of same title on e-print server at cond-mat/9911278. [32] S. A. Morgan, R. J. Ballagh, and K. Burnett. Solitary wave solutions to nonlinear Schr¨odinger equations. Phys. Rev. A, 55:4338–4345, 1997. [33] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C. Cambridge University Press, Cambridge, 1988. [34] N.-P. Proukakis. Microscopic Mean Field Theory of Trapped BoseEinstein Condensates. PhD thesis, University of Oxford, 1997.
BIBLIOGRAPHY
118
[35] H. Steck, M. Naraschewski, and H. Wallis. Output of a pulsed atom laser. Phys. Rev. Lett., 80(1):1, January 1998. [36] D. R. Tilley and J. Tilley. Superfluidity and Superconductivity. Institute of Physics Publishing, Bristol and Philadelphia, 1990. [37] E. B. Treacy. Adiabatic Inversion with Light Pulses. 27A(7):421, 1968.
Phys. Lett.,
[38] Joseph T. Verdeyen. Laser Electronics. Prentice Hall, New Jersey, 1995. [39] H. M. Wiseman. Defining the (atom) laser. Phys. Rev. A, 56:2068, 1997. [40] W. Zhang and D. F. Walls. Gravitational and collective effects in an output coupler for a Bose-Einstein condensate in an atomic trap. Phys. Rev. A, 58:4248, 1998. [41] Weiping Zhang, Karl-Peter Marzlin, Leon Tribe, and Barry Sanders. Collisional and Collapse Dynamics of a Twin Bose-Einstein Condensate with Negative Scattering Length. in preparation, 1999.
