Splitting a Bose Condensate by Mohammadsadegh Mashayekhi A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Physics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August, 2009 c© Mohammadsadegh Mashayekhi 2009 Abstract The primary goal of this thesis is to find the distribution function of number of particles imbalance for ultra cold bosons trapped in a double well poten- tial elongated in one direction for different strength of interaction between particles. This distribution function has been found to be Gaussian distribu- tion function in two different limits. The first limit is weak interaction limit, where we only consider one energy level per well that is called ”two mode approximation” regime. The second limit is when the interaction energy is in the same order as gap between energy levels where in this case, we consider finite number of levels per well. The standard deviation of number of parti- cles distribution also has been found in both limits to be proportional to √ N where N is the total number of particles. In addition, some numerical work is done to find these distribution functions and the results are in agreement with analytical results. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . v Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 0.1 Experimental View . . . . . . . . . . . . . . . . . . . . . . . . 3 0.2 Schematic Picture of the Problem . . . . . . . . . . . . . . . 7 0.3 Javaneinen-Ivanov Model . . . . . . . . . . . . . . . . . . . . 9 1 ”Two-Mode Approximation” Regime . . . . . . . . . . . . . 15 1.1 Analytical Approach . . . . . . . . . . . . . . . . . . . . . . . 15 1.2 J-I Model Vs Josephson Model . . . . . . . . . . . . . . . . . 20 1.3 Numerical Approach . . . . . . . . . . . . . . . . . . . . . . . 23 2 Finite Number of Energy Levels . . . . . . . . . . . . . . . . . 28 2.1 Two Energy Levels per Well . . . . . . . . . . . . . . . . . . . 28 2.2 Finite Number of Energy Levels per Well . . . . . . . . . . . 40 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 iii List of Figures 1 Double Well Potential . . . . . . . . . . . . . . . . . . . . . . 3 2 The Experiment Process . . . . . . . . . . . . . . . . . . . . . 6 1.1 Distribution Function for δ = 0.3 . . . . . . . . . . . . . . . . 24 1.2 Distribution Function for δ = 1 . . . . . . . . . . . . . . . . . 24 1.3 Distribution Function for δ = 3 . . . . . . . . . . . . . . . . . 25 1.4 The Difference Between Gaussian Function and Distribution Function for δ = 0.3 . . . . . . . . . . . . . . . . . . . . . . . . 25 1.5 The Difference Between Gaussian Function and Distribution Function for δ = 1 . . . . . . . . . . . . . . . . . . . . . . . . 27 1.6 The Difference Between Gaussian Function and Distribution Function for δ = 3 . . . . . . . . . . . . . . . . . . . . . . . . 27 iv Acknowledgements Firstly, I would like to thank my research supervisor, Dr. Ian Affleck, for his guidance and support during my Master program and for inspiring me with his work ethic and enthusiasm for his science. I would also like to thank Dr. Fei Zhou, who generously agreed to be my second reader, for which I am very grateful. Many thanks also to Dr. Joseph Thywissen for many helpful discussions. I would also like to thank my parents and my sisters, for their endless emotional support throughout my life. Many thanks as well to my close friend, Hamed Karimi and to all my other good friends, who have been like a family to me while I was so far away from home. v Dedication To my parents and my sisters, To my Father whose support and encouregment have been always with me. To my Mother who brings calm and hope to my life and I would not have been here without her. To my Sisters for their endless emotional support throughout my life. vi Introduction Interferometers are used to measure different physical quantities such as fun- damental constants or characterize atomic and molecular properties. We can have matter wave interferometers as well as optical interferometers. Using this kind of interferometers will get more accurate measurement results be- cause of the sensitivity of atoms to the environment. To make such a matter wave interferometer, we should split a wave packet of atoms into two or more parts coherently and then look at the interference pattern of these different parts [1]. People have found that they can have better interferometers if they use Bose-Einstein condensates as wave packet of matter because condensates has a uniform phase. On the other hand, all the particles in these condensates have the same phase. The relation between thermal matter wave interfer- ometers and Bose-Einstein condensates is the same as the relation between optical interferometers and lasers. So the result would be more accurate. Recently, physicists became interested in this kind of interferometers and they have done many experiments by using these devices. They first confine bosons by using magnetic and optical traps and cool down the particles to very low temperature (in order of 100 nK) to reach Bose-Einstein conden- sation regime. Then, they split this ultra-cold condensate coherently to two condensates and after that at some point, they release the condensates and look at the phase interference of them at different times. The splitting pro- cess can be done with different methods [2]. The first method is called ”optical double well potential” when in this case lasers are used to split the condensate. The second method is called ”purely magnetic double well potential” and in this method magnetic atom chips are used to change the potential. Atom chips are micro-fabricated de- vices that integrate optical elements such as mirrors, waveguides and so on 1 Introduction for matter waves and allow precise alignment of these elements. There are also some other methods to do splitting process in literature [1]. It should be pointed out that in both cases, after coherent splitting, at the end we would have two separated condensates with well-defined relative phase. Many papers has been published on the dynamics of splitting process [3, 4]. They consider a double well potential with some bosonic atoms in each well. Atoms can tunnel from one side to the other side with some probability that is proportional to the shape of the barrier during splitting. The inter- action between particles in each well and particles in different wells are also taken into account. By considering these two type of energy(tunneling energy and interaction energy) we can see beautiful phenomena such as Josephson effect. In the next section, I will talk a little about experimental aspect of this problem and I will show the procedure of this experiment. In section(0.2), I have shown different regimes for this problem according to diffrenet en- ergy scales of different materials. In section(0.3), I follow the calculations by Javaneinen-Ivanov and review their reults that are published in their paper [4]. Chapter 1 mostly contains our calculations in the same regime as J-I model, ”two mode approxiamtion”. So that in Section(1.1), I rederive the results by J-I using a little different calculations and in section(1.2), I com- pare J-I model and Josephson model and in writing this section, again I have used the discussions in J-I paper [4]. The last section of this chapter includes our numerical work in the same regime. In the last chapter, we will find the distribution function and its standard deviation for finite number of energy levels per well. In section(2.1), we start with a simpler case where each well has only two energy levels and in sec- tion(2.2), we generalize the result we found in section(2.1) to finite number of levels in each well. 2 0.1. Experimental View -10 -5 5 10 70 80 90 100 110 Figure 1: The shape of External Potential during splitting process 0.1 Experimental View In this phenomena we have two condensates connected to each other through tunneling. The Josephson Hamiltonian would be: H = −EJcosφ+ Ec 2 M2 (1) where EJ is the tunneling energy and Ec is the interaction energy. Here φ is the phase difference between two condensates (φ = φL − φR) and M is the difference of number of particles in two wells (M = NL−NR). According to the Josephson theory, these two quantities are conjugate of each other and they have following commutation relation: [φ,M ] = i (2) So, we can write uncertainty relation for them as: ∆φ∆M ≥ 1 2 (3) where ∆φ is the standard deviation of φ (∆φ = √ < φ2 >). It means that if we know the relative phase of the two condensates exactly, there would be a large number fluctuation through tunneling between two condensates and if we fix the number of the particles in each well, we would have two completely incoherent condensates. Let us consider again the splitting process. At the beginning, about 103− 105 ultra-cold bosons are trapped in a quasi one-dimensional harmonic 3 0.1. Experimental View trap with following potential (Fig.1): Vext = 1 2 mω2xx 2 + 1 2 mω2yy 2 + 1 2 mω2zz 2 (4) where in the above equation, I assumed the frequency is the same for both transverse directions and ωx, ωy ωz. So, we have an elongated trap in z-direction. The typical numbers for these frequencies are ωz = 15Hz and ωperp = 1.38kHz [5]. Then we start to raise a barrier at the middle of the trap by changing the potential in one of the transverse directions say x. If we split the well coherently, at the end of the process, we would have two condensates with a well-defined phase. So, in this case φ is fixed and ∆φ ' 0. As a result, according to uncertainty relation, we have a huge ∆M that means a lots of tunneling trough barrier and dramatic change in the number of the particles in each well. By continuing raising the barrier, since the tunneling energy decreases by increasing height of the barrier, at some point we can ignore tunneling. At this stage, the number of particles in each well is fixed and we have ∆M = 0. Again according to uncertainty relation, ∆φ would be large or on the other hand, the condensates are completely incoherent. This picture makes sense because at this point, two condensates are completely separated. Let us go back to Josephson Hamiltonian and consider small phase dif- ference between two condensates in Josephson Hamiltonian. In this limit, we can expand cosφ term and keep terms up to second order in φ. The constant term is irrelevant and is just a shift in energy. The next term is a quadratic term in φ. The new Hamiltonian could be written as: H = EJ 2 φ2 + Ec 2 M2 (5) The above Hamiltonian has the shape of Harmonic oscillator Hamiltonian and by comparing with the SHO Hamiltonian, the Josephson frequency of this oscillation would be: h̄ωJ = √ EcEJ (6) 4 0.1. Experimental View The full version of this frequency also could be found by considering the full Josephson Hamiltonian in Eq.(1) would be: h̄ωJ = √ EJ(Ec + 4EJ/N2) (7) The typical numbers for these parameters are EJ/h = 100kHz, Ec/h = 0.06Hz and ωJ = 2pi × 80Hz for N ∼ 9.6 × 103 [5]. We will come back to this frequency later when we discuss Javaneinen-Ivanov model. By looking at Eq.(6), it is obvious that by raising the barrier, EJ decreases and ωJ decreases. Now, the question is ”when do we have adiabatic change in this system?”. Adiabaticity in this system means that during the splitting process, system remains in its ground state. Suppose we change the potential in time scale τ . The other time scale is the time that the system needs to adjust its ground state to the new potential trough oscillation between the wells. This time scale is proportional to inverse of Josephson frequency, 1/ωJ . Thus, the adiabaticity condition can be written as: [4, 5]: 1 ωJ τ (8) Since by raising the barrier, 1/ωJ becomes larger and larger, usually peo- ple consider exponentially decreasing behavior for ωJ in time as [5]: ωJ = ω 0 Je − t τ (9) where ω0J is the Josephson frequency at t = 0. Using above equation, we can rewrite the adiabaticity condition as: |ω̇J | ω2J (10) Now, let me explain the experiment by knowing the adiabaticity criterion. We start with a harmonic trap elongated in z-direction and split the trap in x-direction. At the beginning, we have large tunneling amplitude EJ and large ωJ and the adiabaticity condition is satisfied. By continuing splitting, ωJ becomes smaller and smaller and as a result, 1/ωJ becomes larger and larger. At some time that we call it t∗, we reach the point where 1 ωJ(t∗) = τ (11) 5 0.1. Experimental View Figure 2: Different steps of splitting a Bose condensate, The process starts with splitting the trap up to the point where adiabaticity condition fails and then raising the barrier very fast to have a good resolution picture and at the end counting the particles in each side [5] At this time, we miss the adiabaticity, because after that, the system can not adjust itself to remain in ground state. Here is the point that we stop splitting and do measurement. The typical value of τ is a few miliseconds [5]. As it was explained earlier, φ and M are conjugate to each other. So, we can find the distribution function of one of them and find some information about the other one as well. Sometimes people are interested to look at phase interference. In this case, they turn off the external trap and let the con- densates evolve for some time interval and then find the interference pattern of these condensates on the screen and they can observe phase evolution [1, 2]. Similarly, we can look at the number of particles distribution function. In this case, after reaching the point where we loose adiabaticity, we raise the barrier very fast to some large height to ensure that the condensates are completely separated or EJ ' 0. After that, we count the particles in each well and find the number difference , M (see Fig.2). By repeating this exper- iment many times, we can find the distribution function for M and find the number of particles fluctuation ∆M that according to uncertainty relation is related to ∆φ [5]. It should be pointed out that this kind of measurement is equivalent to collapse of the wave function of the system to one of the 6 0.2. Schematic Picture of the Problem eigenstates of the number difference operator. I will talk more about this when I introduced J-I model. In this work we are interested in the second type of experiments and we focus on the number distribution function and its standard deviation. 0.2 Schematic Picture of the Problem Before going through theoretical calculation, I would like to discuss a little about different regimes that we can have in this problem. There are different energy scales in this system that we should take into account to understand the physics of the problem better. To refresh your mind, we have a cigar- shaped harmonic trap, elongated in z-direction and confined in two other directions. I call the z-direction as longitudinal direction and x,y-direction as transverse directions. So we have two important energy scales related to external potential shape. Since the length of the trap is proportional to 1/ω, we should have following relation between ωz and ωperp for this cigar-shaped trap: ωz ωperp (12) where ωperp = √ ω2x + ω 2 y (13) We can solve harmonic potential problem and find the levels of the en- ergy both in longitudinal and transverse direction. The gap between levels in z-direction would be ωz and the gap between the transverse levels would be ωperp. Note that since particles are free to oscillate only in z-direction, we can say roughly that the kinetic energy is in order of ωz. Another important energy scale is interaction energy, Eint. This energy varies for different substances and we have variety of gases with a wide range of interaction energies [6]. It is good to define the ratio of interaction energy to kinetic energy: γ = Eint Ekin (14) 7 0.2. Schematic Picture of the Problem This ratio could be very small, γ 1 or in order of one, γ ≈ 1 like rubidium atoms or very large, γ ∼ 103 − 104 for sodium atoms [7]. Because of the confinement of the trap in transverse direction, most of the time, ωperp is so much greater than all other energy scales (someone can increase ωperp as much as he wants by making trap more narrow in transverse direction). So, we can always assume that: ωperp Eint (15) It is important to notice since the temperature is very low in practice (about 50-100 nK), we ignore temperature in our theoretical calculations. So, always we consider bosons as a condensate. For this reason, bosons can go to excited states only because of Eint. Eq.(15) tells us that bosons always remains in ground state of transverse degree of freedom. On the other hand, they can not see the higher transverse levels. For the rest of this thesis, we assume bosons are in the ground state of transverse direction. But for longitudinal degree of freedom, we can have different regimes by comparing interaction energy and ωz. If Eint ωz, we are in ”two mode approximation” regime where we only consider ground state in two wells. Javaneinen-Ivanov model is based on this approximation. The distribution function of number of particles imbalance in this limit is found to be Gaus- sian as J-I and we showed. Also, the standard deviation for this distribution is proportional to √ N where N is the total number of bosons. I will show the calculation in detail in next chapter. The second regime is when Eint is in the same order as ωz. In this case, bosons can see finite number of levels in longitudinal direction. By general- izing the J-I model, we showed that even in this case, we still have Gaussian distribution function and again the standard deviation is proportional to√ N . The proof will come up in chapter 2 The third and the last regime is when Eint ωz. Here we should con- sider infinite number of levels. We call this regime as continuum limit regime. Unfortunately we can not use the same generalization to find distribution in 8 0.3. Javaneinen-Ivanov Model this regime, because the equations becomes unsolvable and we can not use the same assumptions we already use for two other regimes to simplify the equations. In this case, the suggestion is to use some continuum limit theory like ”Luttinger model”. We tried to use this model to find the standard devi- ation in continuum limit, but unfortunatly we ended up with some divergent integrals that we could not regularize them. I think the problem is that in this model, we only consider low-lying excitations. But, in the procedure to find the syandard deviation, we have some integrals over all the frequencies and this integrals will give us divergent results. Because I was in shortage of time, unfortunately we could not find a substitude for this model to treat this problem in continuum limit. 0.3 Javaneinen-Ivanov Model In this section I will show part of J-I model and derive their results by following their paper in 1999 [4]. As I mentioned in previous section, in this model, they consider a symmetric double well potential at the time when the barrier is high enough to have a small tunneling rate. The single-particle ground state for this system, ψg, is even and belongs to both wells. The first excited state, ψe, is an odd function and again belongs to both wells. By considering the linear combination of these states we can have states ψl and ψr localized only in one well: ψl,r = 1√ 2 (ψg ± ψe) (16) They also consider the delta-like interaction potential as: U(~r1, ~r2) = 4pih̄2a3D m δ(~r1, ~r2) (17) where a3D is the s-wave scattering length and m is the atomic mass. I would like to emaphasize this point here that as we know, delta-function potential is ultra-violet divergent in 3D. So, it is not suitable to use above interaction for such a 3-dimensional system, although J-I have mentioned above interac- tion in their paper. I think they were a little careless in using this interaction. There are some options to get rid of this divergency. One way is to use the interaction that is smeared out a bit (on a scale small compared to the average 9 0.3. Javaneinen-Ivanov Model inter-particle separation) such as Gaussian function or Lorentzian function with full half at maximum smaller than average inter-particle distance. In this case, again we will end up with somthing similar to Hamiltonian(18), even though, the definition of κ would be a little different from what we have now in Eq.(19). On the other hand, we can consider following definition as a perturbative definition. The other option is to use the interaction that Olshani has used in his famous paper (see Eq(1) in his paper) [8]. This interaction regularizes the ultra-violet divergency due to delta-function interaction, In this case, we would have definitions for κ a bit different from J-I definition. But, we will end up with the same Hamiltonian(18). In this model, J-I consider Hilbert space spanned by only two one-particle states and they call this regime as ”two mode approximation”. The many- particle Hamiltonian could be written as: H = 1 2 (e + g)(a † eae + a † gag) + 1 2 (e − g)(a†eae − a†gag) + κ(a†ea † eaeae + a † ga † gagag + a † ea † eagag + a † ga † gaeae + 4a † ea † gaeag) (18) where ag and ae are boson operators for ground and excited state wave func- tions and the other constants are defined for example as: e = ∫ d3rψe(~r) [ − 1 2m ∇2 + Vext(~r) ] ψe(~r) κ = 2pia m ∫ d3r|ψe/g(~r)|4 (19) where Vext is the external harmonic potential in Eq.(4). Note that in writing interaction terms, they have assumed that tunneling probability is small, so, they can ignore the overlap of ψl and ψr in following equation: |ψg,e|2 = 1 2 (|ψl|2 + |ψr|2 ± 2ψlψr) (20) Then for all the interaction terms in Eq.(18), we consider the same in- teraction strength κ. To write everything in terms of left and right side 10 0.3. Javaneinen-Ivanov Model operators, they define annihilation operators for particles in left and right wells as: al,r = 1√ 2 (ag ± ae) (21) The other point is that, the energy difference between ground and excited states, e − g that exist in Hamiltonian (18), could be interpreted as the single-particle tunneling rate and is called δ. κN is also the interaction energy per particle. By plugging Eq.(21) into Eq.(18), we will have a simpler Hamiltonian: H = −δ 2 (a†lar + a † ral)− 4κa†lala†rar (22) Notice that in deriving above equation, we have used the conservation of the total number of the particles in this system so that we have ignored the terms that are only function of N where: N = a†lal + a † rar (23) Consider taking a large number of particles (thermodynamic limit). In that limit, the Hamiltonian and ground state energy scale proportional to N . By looking at Eq.(22) we can roughly find how our parameters scale with N . δ scales like 1 and κ scales like 1/N . We also define κ̃ as: κ̃ = Nκ (24) where κ̃ is held fixed as N → ∞. To check that if this assumption is rea- sonable, it is good to look at experimental parameters. I will discuss later in section (1.2) in detail the comparison between J-I model and Josephson model and we will see that these models are equivalent in some limits if we have below identifications between parameters of these two models: EJ ↔ Nδ 2 , Ec ↔ 2κ (25) By using above equations, we can also find the identification for Nκ/δ that is: Nκ δ ↔ N 2Ec 4EJ (26) 11 0.3. Javaneinen-Ivanov Model Now, we can use experimental data that I mentioned earlier in section (0.1) for Ec and EJ and find a rough estimate for Nκ/δ. For EJ/h = 100kHz, Ec/h = 0.06Hz and N ∼ 104 we will find [5]: Nκ δ ∼ 15 (27) that shows our assumption that Nκ/δ scales proportional to 1 is a reasonable assumption. Another way to see the scaling behavior of κ with N is to look at its definition in Eq.(19). As we explained before, ψ(~r) is the wavefunction in either of the traps. If the size of the unsplit trap is V , |ψ(~r)|2 is of order 1/V . If we formally took N and V to infinity, holding N/V fixed (the usual thermodynamic limit), we would conclude that κ scales proportional to 1/N . Since the number of the particles is conserved and we consider only two single-particle states for above Hamiltonian, we can take into account N + 1 states: |k >≡ |N 2 + k >l |N 2 − k >r k = −N 2 , ..., N 2 (28) as a complete basis and specially we can expand the wavefunction of the system in terms of these states as: |ψ >= N 2∑ k=−N 2 Ck|k > (29) where in the above equation k is NL−NR 2 compared to M that is equal to NL−NR. Note that in Eq.(28), the first and the second brackets show num- ber of particles in left and right wells respectively. Now, I go back to the comment I mentioned in section(0.1) about the col- lapse of the wavefunction when you measure the number of the particles. We know that if we expand a wavefunction in terms of eigenstates of a physical quantity and try to measure that quantity of the system, the wavefunction collapses to one of the eigenstates of that physical operator. Here, we have the same situation when we measure the number of particles in each side. By looking at Eq.(29), it is clear that after each measurement, the wavefunc- tion collapses to one the states |k > and the probability to collapse to this state is |Ck|2. Thus, we conclude that the distribution function of number of particles imbalance is composed by these |Ck|2’s. Finally, it is enough to 12 0.3. Javaneinen-Ivanov Model find these coefficients and then we have our desired distribution function. We can easily find the action of H on the system wavefunction: Hψ = H ∑ k Ck|k > = ∑ k [−δ 2 (√ (N/2 + k) (N/2− k + 1)Ck−1 + √ (N/2 + k + 1) (N/2− k)Ck+1 ) − 4κ (N/2 + k) (N/2− k)Ck ] |k > (30) Note that H-matrix is tridiagonal in this basis. At this point, J-I try to diagonalize this matrix by using some approximations that simplify Eq.(30). When N is large, we expect that Ck will be negligible except for a range of k obeying |k| N . So, we can expand the square roots in Eq.(30). As we shall see this range is of order √ N . Ck is non-negligible for many values of k (still √ N 1). So, we can consider C’s as continuous function of k and expand C(k ± 1) in terms of C(k) as: C(k ± 1) = C(k)± C ′(k) + 1 2 C ′′(k) + · · · (31) After substituting above expansion in Eq.(30) and keeping the leading order in 1/N , we would find the following time-independent Schroedinger equation: [ −Nδ 4 d2 dk2 + δk N d dk + ( δ N + 4κ ) k2 ] C(k) = C(k) (32) By considering the scaling of the parameters respect to N and assuming that k scales like √ N , the first and the third terms in left hand side of the above equation are in order of 1 but the second term scales like 1/N . So, for large N limit, we can ignore the second term. Finally, we have:[ −Nδ 4 d2 dk2 + ( δ N + 4κ ) k2 ] C(k) = C(k) (33) This equation is again similar to Schroedinger equation for the simple harmonic oscillator. The angular frequency that J-I found would be: ω = √ δ(δ + 4κN) (34) 13 0.3. Javaneinen-Ivanov Model In addition, as we know, the ground state of the harmonic oscillator problem is Gaussian in x and p spaces. So, for our problem, we also have a Gaussian ground state in k and d/dk that is proportional to the phase difference φ. The standard deviation of this Gaussian ground state is: σg = √ N 2 ( δ δ + 4κN )1/4 (35) What we are looking for is the probability distribution function and its standard devaition. This probability function would be the square of the ground state function. So, the standard deviation of the probability function would be: ∆k = √ N 2 ( δ δ + 4κN )1/4 (36) Note that, if you take teh square of a Gaussian function, its standard deviation scales with 1/ √ 2 factor. Thus, the distribution function of the number of particles imbalance would be Gaussian and its standard deviation scales proportional to √ N . We derived exactly the same results with a little different calculations in section(1.1). I will show our calculations in detail in next chapter. Also, we found this distribution function numerically and compared with a Gaussian function. The plots of this numerical work will come up at the end of next chapter. 14 Chapter 1 ”Two-Mode Approximation” Regime In section(0.3), I talked about J-I model and derived their results for distribu- tion function of number of particles imbalance in ”two mode approximation” regime. In section(1.1), I will derive the similar results by a little different calculation and show that our results are exactly the same as J-I results. In section(1.2), I will talk more about relation between tight-binding Hamilto- nian (22) and Josephson Hamiltonian and discuss about the frequency and number fluctuations in these two models. In the last section, I will find the distribution function by exact diagonalization of the tight-binding Hamilto- nian and show that it is a Gaussian function. 1.1 Analytical Approach Let’s return to Eq.(18) and this time, first consider the effect of H on only one of the eigenstates of the number imbalance operator, |k >. Similar to Eq.(30) we have: Hψ = −δ 2 (√ (N/2 + k) (N/2− k + 1)|k − 1 > + √ (N/2 + k + 1) (N/2− k)|k + 1 > ) −4κ̃ N (N/2 + k) (N/2− k) |k > (1.1) As I also mentioned in section(0.3), for large N limit, Ck’s are negligible except for a range of k where |k| N . But still Ck is non-negligible for a very large interval that I will show later that this range of k is in order of √ N that is still a very large number. So, first I ignore the constants 1 respect to k, N in square roots of Eq.(1.1) and then, expand the square roots. Let me 15 1.1. Analytical Approach first focus on δ-dependent part, H0. The coefficient of |k − 1 > state can be expanded as: −δ 2 √ (N/2 + k)(N/2− k + 1)|k − 1 >≈ −δ 2 (N/2− k2/N)|k − 1 > = −δN 4 |k − 1 > +δk 2 2N |k − 1 > (1.2) where the first term in the last equation is proportional to N and the second term is proportional to 1, because k scales like √ N . In addition we can expand states |k ± 1 > in terms of |k >, that will give us: |k ± 1 >= |k > ± d dk |k > + · · · (1.3) Then, in Eq.(1.2), I approximate the second term with: δk2 2N |k > (1.4) and ignore the higher order terms because the next term in this expansion is proportional to 1/k or 1/ √ N and for large N limit, it tends to zero. By similar justification for coefficient of |k + 1 > state, finally we can write: H0|k >≈ −δ 2 [ N 2 (|k + 1 > +|k − 1 >)− 2k 2 N |k > ] (1.5) Considering all of the above discussions, we can write time-independent Schroedinger equation for Ck’s as: −δN 4 [Ck+1 + Ck−1] + [δ + 4κ̃] k2 N Ck − κ̃NCk = ECk (1.6) Note that since the total number of particles is conserved, the last term in the left hand side of the above equation is constant and could be absorbed in E by changing the origin of the energy. Let me now add the following term to the both sides of the Eq.(1.6): Nδ 2 Ck (1.7) 16 1.1. Analytical Approach Thus, we can rewrite the Schroedinger equation as: −δN 4 [Ck+1 + Ck−1 − 2Ck] + [δ + 4κ̃]k 2 N Ck = [E + κ̃N + Nδ 2 ]Ck (1.8) Next, we approximate Ck by a continuous function, C(k) and use: Ck+1 + Ck−1 − 2Ck ≈ C(k + 1) + C(k − 1)− 2C(k) ≈ d 2C dk2 (1.9) where C(k) was Taylor expanded to second order, assuming that it varies smoothly where k changes by ±1. This assumption is reasonable at large N because, as we shall see, C(k) only varies significantly when k changes by amount of order √ N . Thus, the discrete Schroedinger equation reduces to a continuous one, the Harmonic Oscillator equation: −Nδ 2 d2C dk2 + (δ + 4κ̃ N ) k2C(k) = C(k) (1.10) where = E + κ̃N + Nδ 2 (1.11) We may define a rescaled variable: x ≡ k√ N (1.12) So that the equation becomes −δ 4 d2C dx2 + (δ + 4κ̃)x2C(x) = C(x) (1.13) Then, the frequency of this oscillator would be: ω = √ δ(δ + 4κ̃) (1.14) and the ground state wavefunction is a Gaussian: C(x) = (4(δ + 4κ̃) pi2δ )1/8 exp [ − √ δ + 4κ̃ δ x2 ] (1.15) 17 1.1. Analytical Approach Then square of the wavefunction determines the probability distribution function for x (i.e. k): |C(x)|2 = (4(δ + 4κ̃) pi2δ )1/4 exp [ −2 √ δ + 4κ̃ δ x2 ] (1.16) Thus, the distribution is Gaussian with standard deviation: √ < x2 > = 1 2 ( δ δ + 4κ̃ )1/4 (1.17) Note that in above equation, x scales proportional to N0 = 1 that is consistent with the assumption we made before. Then, we can write the standard deviation for k as: √ < k2 > = √ < (NL −N/2)2 > = √ N 2 ( δ δ + 4κ̃ )1/4 (1.18) and equivalently: √ < (NL −NR)2 > = √ N ( δ δ + 4κ̃ )1/4 (1.19) where by comparing Eq.(1.18) with Eq.(36), it is clear that our standard deviation is exactly the same as J-I one. It is good to talk a little about extreme limits for Hamiltonian (22). Let’s start by zero tunneling. In this limit, the Hamiltonian would be: H = −4κa†lala†rar (1.20) then the states that we defined in Eq.(28), are eigenstates of the above Hamil- tonian and the eigenenergies could be found easily to be: ωk = κ(4k 2 −N2) (1.21) For ground state, from above equation, k should be zero. It means that the number of particles in both wells is the same and is equal to N/2. In ad- dition, obviously there is no fluctuations in number of particles in the absence of tunneling term. This picture is consistent with our number fluctuation in 18 1.1. Analytical Approach Eq.(1.19) when we put δ = 0 in that equation. The other limit that is more interesting, is the non-interacting case. In this limit, the Hamiltonian reduces to: H = δ 2 (a†eae − a†gag) (1.22) where I have written again H in terms of creation operators for one-particle ground and excited states. The above Hamiltonian shows that, in this limit, the ground state of the system, is where all the particles are in ψg. It can be shown as: 1√ N ! (a†g) N |vac > = 1√ N ! ( a†l + a † r√ 2 )N |vac > = 1 2N/2 N/2∑ k=−N/2 √( N N/2 + k ) |k > (1.23) Note that the last term in the above equation, shows binomial distribu- tion function and it could be shown easily that this distribution function tends to Gaussian distribution function in large N limit. It should be pointed out that again in this limit, the ground state of the system, has an equal number of particles in each well. The difference of this limit and previous one is that in the first case (zero-tunneling), the num- ber fluctuation is zero but in the second limit (no interaction), the number fluctuation exists and according to Eq.(1.19), is equal to √ N when we put κ̃ = 0 in that equation. We call this distribution function as Poissonian, because ∆N = √ N . By turning on the interaction, the coefficient of √ N in Eq.(1.19) becomes smaller and smaller up to the point where interaction term dominates and we reach the first limit where the number fluctuations would be zero. Thus, for all the values of κ and δ, other than non-interacting one, the distribution function would be ”Sub-Poissonian” for the number of particles imbalance, M = NL −NR. There are various corrections to the results we found for finite N . For example, we could calculate the leading correction in 1/N by keeping terms to 19 1.2. J-I Model Vs Josephson Model next order in 1/N in Eq.(1.1) and going to next order in the Taylor expansion in Eq.(1.9) which would give a d4/dk4 term in the Hamiltonian. 1.2 J-I Model Vs Josephson Model In section(0.1), I talked about Josephson Hamiltonian: HJ = −EJcosφ+ 2Eck2 (1.24) where I changed M in Eq.(1) to k in above equation. I mentioned that in this model, EJ is called the Josephson coupling energy and Ec the capacitive energy. I showed also J-I model in this chapter with Hamiltonian: HJI = −δ 2 (a†lar + a † ral)− 4κa†lala†rar (1.25) Now, it would be interesting to look at the similarities and differences between these two models. I should point out that in this section I mainly follow the discussion in J-I paper [4]. First of all, by using the eigenstates of number imbalance operator, k, in Eq.(28), we can write the eigenstates of the phase difference operator as: |φp >= 1√ N + 1 N/2∑ k=−N/2 eikφp |k > (1.26) and we choose φp to be: φp = 2pip N + 1 , p = −N/2, · · · , N/2 (1.27) These eigenstates also form a complete basis and we can span the same Hilbert space as |k > basis. The phase operator in this notation, could be written as: φ̂ = ∑ p |φp > φp < φp| (1.28) and number of particles imbalance operator could be written in |k > basis as: 20 1.2. J-I Model Vs Josephson Model k̂ = ∑ k k|k >< k| (1.29) By using Eqs.(1.28) and (1.29), the spectral representation of the two- mode version of the Josephson Hamiltonian(1.24) would be: HJ = −1 2 EJ ∑ k (|N/2 + k + 1 >< N/2 + k|+H.c.) + 2Ec ∑ k k2|N/2 + k >< N/2 + k| (1.30) We can also write J-I Hamiltonian(1.25), in these basis as: HJI = −δ 2 ∑ k (√ (N/2 + k + 1)(N/2− k)|N/2 + k + 1 >< N/2 + k|+H.c. ) + 4κ ∑ k k2|N/2 + k >< N/2 + k| (1.31) By comparing Eqs.(1.30) and (1.31), it is obvious that the difference between these two models is small. These two Hamiltonian will give ap- proximately the same results, if we consider below identifications between parameters of these two models: EJ ↔ Nδ 2 , Ec ↔ 2κ (1.32) Another way to see the similarity between these two models, is to write Josephson Hamiltonian in Harmonic Oscillator form. To achieve this goal, we first consider: k ↔ −i d dφ (1.33) where k and φ are conjugate operators. In addition, we write the wavefunc- tion of the system as: ψ(φ) = ∑ k Cke ikφ (1.34) 21 1.2. J-I Model Vs Josephson Model Note that we choose above expression for wavefunction in terms of φ, because it should be periodic in φ. Now, let’s derive the time-independent Schroedinger equation: HJψ = Eψ (1.35) By using following equation: cosφeikφ = 1 2 [ ei(k−1)φ + ei(k+1)φ ] (1.36) and changing the dummy variable in sum over k, we can rewrite the Schroedinger equation as: −EJ 2 [Ck+1 + Ck−1] + 2Eck2Ck = ECk (1.37) If we follow the same steps as we did for Eq.(1.6), we will end up with the following Schroedinger equation for Josephson Hamiltonian:(−Nδ 4 d2 dk2 + 4κ̃k2 N ) C(k) = C(k) (1.38) Above equation is very similar to Eq.(1.10) for J-I model. So, it tells us that again the distribution function is a Gaussian function and in the presence of both tunneling and interaction, the standard deviation scales proportional to √ N . The only difference is δk2/N term that does not appear in above equa- tion. The problem that we encounter is that for Josephson Hamiltonian the frequency from Eq.(1.38) would be: ω = √ 4δκ̃ (1.39) and the number fluctuation is given by: √ < (NL −NR)2 > = √ N ( δ 4κ̃ )1/4 (1.40) For non-interacting case (κ̃ = 0), the number fluctuation diverges that is not the behaviour we expect. So, it seems that Josephson model does not work well for small interactions. But, in the other limit, we see the expected behavior that is zero number fluctuation. 22 1.3. Numerical Approach 1.3 Numerical Approach In this section, I will show our numerical results for number of particles imbalance distribution function. Let’s return to Eq.(1.1) that shows the action of J-I Hamiltonian on state |k >. We can simply find matrix elements of H in this basis as: Hk,k = −4κ (N2 4 − k2 ) (1.41) Hk,k+1 = −δ 2 √ (N/2− k)(N/2 + k + 1) = Hk+1,k (1.42) where k goes from −N/2 to N/2. Note that H is tridiagonal and symmetric in this basis. So, we can simply find the eigenvalues and eigenvectors of this system by using mathematical programs. We used Matlab program to diag- onalize this Hamiltonian. As I explained in detail in section(0.3), what we need are the coefficients of the expansion of the ground state of the system in terms of |k > states. Since we wrote Hamiltonian in |k > basis, it is enough to find the ground state wavefunction in this basis that is a vector with N + 1 elements. Let me call these elements, Ck where k goes from −N/2 to N/2. Then, the plot of |Ck|2 versus k would be our desired distribution function. We found this distribution function for different values of tunneling am- plitude, δ. In all the graphs, we keep Nκ = κ̃ = 1. In Figs.(1.1),(1.2),(1.3), you see the distribution function for δ = 0.3, 1, 3 and for different total num- ber of particles, N = 5000, 6000, · · · , 10000. Note that we chose δ and κ̃ approximately in same order to see the effect of tunneling and interaction together. We also depicted the real Gaussian function to compare our results with this function. Another point is that we have scaled all the graphs by 1/ √ N . As you can see, it is very hard to distinguish between different graphs in this figures and all the plots match the Gaussian function. So, it shows that the distribution function for different values of δ and for different number of particles is Gaussian. In inset plot, the enlarged view of graphs is shown 23 1.3. Numerical Approach Figure 1.1: The distribution functions for δ = 0.3 and κ̃ = 1 for different number of particles form 5000 to 10000 Figure 1.2: The distribution functions for δ = 0.3 and κ̃ = 1 for different number of particles form 5000 to 10000 24 1.3. Numerical Approach Figure 1.3: The distribution functions for δ = 3 and κ̃ = 1 for different number of particles form 5000 to 10000 Figure 1.4: The diffrence between Gaussian function and the distribution func- tions for δ = 3 and κ̃ = 1 for different number of particles form 5000 to 10000 25 1.3. Numerical Approach where different colors show different number of particles. As we discussed before, when N goes to infinity, according to ”Cen- teral limit theorem”, any random distribution function tends to Gaussian distribution function. To see the effect of finite number of particles more clearly, we also plot the diffrence between distribution function for finite number of particles and infinite number, that is a Gaussian function, for N = 5000, · · · , 10000 that they are shown in Figs.(1.4), (1.5),(1.6) for differ- ent values of δ. Note that in all the graphs the vertical axes is in order of 10−5 that indicates the small corrections due to finite size effects. 26 1.3. Numerical Approach Figure 1.5: The diffrence between Gaussian function and the distribution func- tions for δ = 3 and κ̃ = 1 for different number of particles form 5000 to 10000 Figure 1.6: The diffrence between Gaussian function and the distribution func- tions for δ = 3 and κ̃ = 1 for different number of particles form 5000 to 10000 27 Chapter 2 Finite Number of Energy Levels In this chapter, I will talk about distribution function of number of particles imbalance when Eint is in the same order as longitudinal frequency of trap. In contrast with ”Two mode approximation” regime, here bosons can see finite number of energy levels in each well. I will show in this chapter that even in this case, the distribution function is Gaussian and the standard deviation of number of particles imbalance scales with √ N . As the first step, In section(2.1), I will consider two levels of energy in each well that means totally we have 4 energy levels. In section(2.2), I will generalize the calculation in section(2.1) for finite number of energy levels. 2.1 Two Energy Levels per Well In this section, we consider two energy levels per well. I show the parame- ters and operators related to each level with labels 1, 2. The most general Hamiltonian for such a system, could be written as: H = [1n̂1L + 2n̂2L + (L→ R)] + [δ11a † 1La1R + δ22a † 2La2R + δ12(a † 1La2R + a † 2La1R) + h.c.] + [κ11n̂1Ln̂1L + κ22n̂2Ln̂2L + κ12n̂1Ln̂2L + (∆1122a † 1La † 1La2La2L + h.c.) + (∆1112n̂1La † 1La2L + h.c.) + (∆2221n̂2La † 1La2L + h.c.) + (L→ R)]. (2.1) where aiL/R annihilates a boson in state i in the left/right well and n̂iL/R = a†iL/RaiL/R. Note that the ∆ijkl results from collisions between bosons, in different states k and l (in the same well) which sends them into states i and j. Note that, we have a symmetrical trap. On the other hand, we have 28 2.1. Two Energy Levels per Well assumed parity symmetry for this problem (L and R are equivalent). Additionally, in the above Hamiltonian only the energy difference 1− 2 matters at T = 0. We can always choose the δij to be real by redefining the relative phases of left and right annihilation operators. The ∆ijkl can in general be complex. This gives 13 parameters for this problem. Pay atten- tion that we have ignored the interaction between bosons in different wells, because these interactions are magnetic dipole-dipole interactions and are very much smaller than other energy scales of the problem for those matters that we are interested in [9]. In the large N limit (thermodynamic limit), the Hamiltonian and ground state energy scale proportional to the total number of particles. So, by looking at Eq.(2.1), we can simply find how our parameters scale with N . 1, 2, δ11, δ22, δ12 scale like 1 and κ11, κ22, κ12,∆1122,∆1112,∆2221 scale pro- portional to 1/N . It should be pointed out that this assumed scaling of parameters with N , may not always be true, but it leads to nice theoretical behaviour. Note that, one of the differences between this regime and two mode regime is that in the latter, the total number of particles for each energy level is fixed and equal to N . But in the former, the total number of particles in each en- ergy level, N1 = N1L + N1R and N2 = N2L + N2R, are not fixed, because we have tunneling between two levels. Thus, we can not use the same trick as two mode regime to ignore terms only proportional to N as a change in energy origin. As a result, we should have more complicated equations in this regime. The number of bosons in the left and right trap in each level could be shown as: Ni,L,R = Ni 2 ± 2ki (2.2) We also need to define another quantum number for this problem. We define it as: m = M− < M > (2.3) where 29 2.1. Two Energy Levels per Well M = N1 −N2 2 (2.4) and < M > is the average value of M . On the other hand, M is the number of particles difference between different levels. This number could be changed through tunneling between different levels. To find a rough idea about how does < M > scale with N , we can only consider energy terms and κij terms in Hamiltonian(2.1) and minimize the energy of the system. The value of M at which the energy of the system is minimum, could be easily found to be: M0 = 4δ11 − 4δ22 + 41 − 42 + κ11N − κ22N −2κ11 + 2κ12 − 2κ22 (2.5) Note that, all the terms in numerator of above equation scale proportional to 1 and all the terms in denominator scale like 1/N . Thus, M roughly scales proportional to N . Now, number of particles in each energy level could be written in terms of < M > and m as: N1 = N + 2 < M > 2 +m ≡ α +m (2.6) N2 = N − 2 < M > 2 −m ≡ β −m (2.7) where I have defined two new parameters, α, β for simplicity. Note that they are constants and both scale proportional to N . In the next step, we should consider a complete basis that could span the Hilbert space and then write the Hamiltonian of the system in this basis and solve the problem. As I discussed before, in this problem, we have four kinds of creation operators, a†i,L,R, that they act on different spaces. On the other hand, they commute with each other. In addition, since the total number of particles, N , is fixed, we have a constraint on these operators. Thus, totally we need three quantum numbers to describe this system. So, we consider the following complete set of states: |k1, k2,M〉 (2.8) We can expand the eigenfunctions of our system in terms of these states: 30 2.1. Two Energy Levels per Well |Ψ〉 = Ni/4∑ ki=−Ni/4 N/2∑ M=−N/2 Ck1,k2,M |k1, k2,M〉 (2.9) where Ck1,k2,M ’s are coefficients of this expansion. By using the same justification as previous chapter, in large N limit, we expect that Ck1,k2,M be negligible except for the range of ki and m in order of √ N . This assumption is consistent with what I will show later. Note that again this assumption may not always be true but it helps us to have a nice theory. The action of Hamiltonian(2.1) on each of the states |k1, k2,M〉 would be: H |k1, k2,M〉 = [1(α+m) + 2(β −m)] |k1, k2,M〉 + δ11( √ α 2 + m 2 − 2k1 √ α 2 + m 2 + 2k1 + 1) |k1 + 2, k2,M〉 + δ11( √ α 2 + m 2 − 2k1 + 1 √ α 2 + m 2 + 2k1) |k1 − 2, k2,M〉 + δ22( √ β 2 − m 2 − 2k2 √ β 2 − m 2 + 2k2 + 1) |k1, k2 + 2,M〉 + δ22( √ β 2 − m 2 − 2k2 + 1 √ β 2 − m 2 + 2k2) |k1, k2 − 2,M〉 + δ12( √ β 2 − m 2 − 2k2 √ α 2 + m 2 + 2k1 + 1) |k1 + 1, k2 + 1,M + 1〉 + δ12( √ α 2 + m 2 + 2k1 √ β 2 − m 2 − 2k2 + 1) |k1 − 1, k2 − 1,M − 1〉 + δ12( √ α 2 + m 2 − 2k1 √ β 2 − m 2 + 2k2 + 1) |k1 + 1, k2 + 1,M − 1〉 + δ12( √ β 2 − m 2 + 2k2 √ α 2 + m 2 − 2k1 + 1) |k1 − 1, k2 − 1,M + 1〉 + κ11[( α 2 + m 2 + 2k1)2 + ( α 2 + m 2 − 2k1)2] |k1, k2,M〉 + κ22[( β 2 − m 2 + 2k2)2 + ( β 2 − m 2 − 2k2)2] |k1, k2,M〉 + κ12[( α 2 + m 2 + 2k1)( β 2 − m 2 + 2k2) + ( α 2 + m 2 − 2k1)(β2 − m 2 − 2k2)] |k1, k2,M〉 + ∆1122 √ β 2 − m 2 + 2k2 √ β 2 − m 2 + 2k2 − 1 √ α 2 + m 2 + 2k1 + 1× × √ α 2 + m 2 + 2k1 + 2 |k1 + 2, k2 − 2,M + 2〉 31 2.1. Two Energy Levels per Well + ∆1122 √ β 2 − m 2 + 2k2 + 1 √ β 2 − m 2 + 2k2 + 2 √ α 2 + m 2 + 2k1× × √ α 2 + m 2 + 2k1 − 1 |k1 − 2, k2 + 2,M − 2〉 + ∆1122 √ β 2 − m 2 − 2k2 √ β 2 − m 2 − 2k2 − 1 √ α 2 + m 2 − 2k1 + 1× × √ α 2 + m 2 − 2k1 + 2 |k1 − 2, k2 + 2,M + 2〉 + ∆1122 √ β 2 − m 2 − 2k2 + 1 √ β 2 − m 2 − 2k2 + 2 √ α 2 + m 2 − 2k1× × √ α 2 + m 2 − 2k1 − 1 |k1 + 2, k2 − 2,M − 2〉 + ∆1112 √ β 2 − m 2 + 2k2 √ α 2 + m 2 + 2k1 + 1( α 2 + m 2 + 2k1 + 1) |k1 + 1, k2 − 1,M + 1〉 + ∆1112 √ β 2 − m 2 + 2k2 + 1 √ α 2 + m 2 + 2k1( α 2 + m 2 + 2k1) |k1 − 1, k2 + 1,M − 1〉 + ∆1112 √ β 2 − m 2 − 2k2 √ α 2 + m 2 − 2k1 + 1(α2 + m 2 − 2k1 + 1) |k1 − 1, k2 + 1,M + 1〉 + ∆1112 √ β 2 − m 2 − 2k2 + 1 √ α 2 + m 2 − 2k1(α2 + m 2 − 2k1) |k1 + 1, k2 − 1,M − 1〉 + ∆2221 √ β 2 − m 2 + 2k2 √ α 2 + m 2 + 2k1 + 1( β 2 − m 2 + 2k2 − 1) |k1 + 1, k2 − 1,M + 1〉 + ∆2221 √ β 2 − m 2 + 2k2 + 1 √ α 2 + m 2 + 2k1( β 2 − m 2 + 2k2) |k1 − 1, k2 + 1,M − 1〉 + ∆2221 √ β 2 − m 2 − 2k2 √ α 2 + m 2 − 2k1 + 1(β2 − m 2 − 2k2 − 1) |k1 − 1, k2 + 1,M + 1〉 + ∆2221 √ β 2 − m 2 − 2k2 + 1 √ α 2 + m 2 − 2k1(β2 − m 2 − 2k2) |k1 + 1, k2 − 1,M − 1〉 (2.10) Note that, from the action of creation and annihilation operators on Hamiltonian(2.1), we would have 17 different states in above terms. Now, we try to expand the square roots in coefficients of above equation, by using the fact that m and ki scale proportional to √ N . Again, we neglect the constants 1, 2 in some of the square roots respect to α, β that scale with N . For example, we can write this expansion for one of the coefficients as: 32 2.1. Two Energy Levels per Well √ α 2 + m 2 − 2k1 √ α 2 + m 2 + 2k1 + 1 ≈ α 2 (1 + 1 2 m α − 2k1 α − 1 8 m2 α2 − 2k 2 1 α2 )(1 + 1 2 m α + 2k1 α − 1 8 m2 α2 − 2k 2 1 α2 ) ≈ α 2 (1 + m α − 8k 2 1 α2 ) (2.11) where we only keep terms up to order 1 in the expansion. Additionally, we can Taylor expand 17 states around |k1, k2,M〉 state. For example we have: |k1 + 1, k2 + 1,M + 1〉 = |k1, k2 + 1,M + 1〉+ ∂k1 |k1, k2 + 1,M + 1〉+ 1 2 ∂2k1 |k1, k2 + 1,M + 1〉 = |k1, k2,M + 1〉+ ∂k2 |k1, k2,M + 1〉+ 1 2 ∂2k2 |k1, k2,M + 1〉 + ∂k1 |k1, k2,M + 1〉+ ∂k1∂k2 |k1, k2,M + 1〉+ 1 2 ∂2k1 |k1, k2,M + 1〉 = |k1, k2,M〉+ (1 2 ∂2m + 1 2 ∂2k1 + 1 2 ∂2k2) |k1, k2,M〉+ (∂m + ∂k1 + ∂k2) |k1, k2,M〉 (2.12) where again we only keep terms up to order 1. Pay attention that ∂m, ∂ki are in order of 1/ √ N . After we expand all the coefficients and states, we can write the time- independent Schroedinger equation for such a system by considering the ac- tion of Hamiltonian(2.1) on eigenfunction |Ψ〉 in Eq.(2.9). Note that after Taylor expanding the states, we only have one state in our equation, |k1, k2,M〉. Thus, after ignoring sum over ki and M , the Schroedinger equation would be:( m[1 − 2 + ακ11 − βκ22 − 12ακ12 + 1 2 βκ12 + δ11 − δ22 + δ12 √ αβ( 1 α − 1 β ) + (β − α)∆1122 + 14 √ αβ(∆1112 −∆2221) + √ αβ/4( 1 α − 1 β )(α∆1112 + β∆2221)] +m2[ 1 2 κ11 + 1 2 κ22 − 12κ12 − δ12 4 ( 1 α + 1 β )2 √ αβ −∆1122 −∆1112 √ αβ α 8 ( 1 α + 1 β )2 33 2.1. Two Energy Levels per Well + √ αβ/4( 1 α − 1 β )(∆1112 −∆2221)−∆2221β8 √ αβ( 1 α + 1 β )2] + k21[8κ11 − 8δ11 α − 4δ12 α2 √ αβ + 6 α √ αβ∆1112 − 2β α2 √ αβ∆2221] + k22[8κ22 − 8δ22 β − 4δ12 β2 √ αβ + 6 β √ αβ∆2221 − 2α β2 √ αβ∆1112] + k1k2[8κ12 − 8 αβ δ12 √ αβ + 16∆1122 + 12 β √ αβ∆1112 + 12 α √ αβ∆2221] ) Ck + ( ∂2m[δ12 √ αβ + 2αβ∆1122 + √ αβ/4α∆1112 + √ αβ/4β∆2221] + ∂2k1 [2αδ11 + δ12 √ αβ + 2αβ∆1122 + √ αβ/4α∆1112 + √ αβ/4β∆2221] + ∂2k2 [2βδ22 + δ12 √ αβ + 2αβ∆1122 + √ αβ/4α∆1112 + √ αβ/4β∆2221] + ∂k1∂k2 [2δ12 √ αβ − 4αβ∆1122 − α √ αβ∆1112 − β √ αβ∆2221] ) Ck = ÉCk (2.13) where I have used a short notation for coefficients as: Ck1,k2,M ≡ Ck (2.14) In above Schroedinger equation, all the terms are quadratic in ki,m, ∂ki and ∂m except the first term that is linear term in m. To get rid of this term, we can use complete square method and omit the linear term in m and consider a shift in origin of m. On the other hand, if we show the terms proportional to m and m2 as γ1m+ γ2m 2, we can write: γ1m+ γ2m 2 = γ2(m 2 + γ1 γ2 m) = γ2(m+ γ1 2γ2 )2 − γ 2 1 4γ2 (2.15) where in the last equation, the second term is constant and could be consid- ered as a shift in energy by redefining energy of the system as below: E = É + γ21 4γ2 (2.16) After changing m to m − γ1 2γ2 , that is equal to shift of the probability profile in m without any change in the scale and shape of the profile, all the terms in Schroedinger equation would be quadratic. In the rest of this section, I will try to show that this Schroedinger equa- tion is similar to Harmonic Oscillator Schroedinger equation and as a result, 34 2.1. Two Energy Levels per Well the ground state of the system would be Gaussian. To achieve this goal, I rewrite Eq.(2.13) in following matrix form: ∂ ∂X T B ∂ ∂X Ck +X TAXCk = ECk (2.17) where A, B are 3×3 symmetric matrices (look at Eq.(2.13)) and A elements scale like 1/N and B elements sclae proportional to N. In addition, X vector, has three components m, k1, k2. To simplify above differential equation, we use several transformations. First of all, we should diagonalize the derivative term. So, we use the following orthonormal transformation: B = UDU−1 (2.18) where U is the matrix that is formed by eigenvectors of B, and D is the diagonal matrix that is formed by eigenvalues of B. Notice that again D ele- ments scale proportional to N, because the eigenvalues of a matrix scale like the elements of that matrix. In addition, U elements scale like 1, because it is orthonormal transformation matrix and all the eigenvectors are normalized to one. By plugging above transformation into Eq.(2.17), we would have: ∂ ∂X T UDU−1 ∂ ∂X Ck +X TAXCk = ECk (2.19) Now, we define a new variable, Y = UX. So we can also write U−1 ∂ ∂X = ∂ ∂Y and X = U−1Y. After changing the variables , the new equation would be: ∂ ∂Y T D ∂ ∂Y Ck +Y TUAU−1YCk = ECk (2.20) We again redefine Y to absorb the D matrix as: Z ≡ 1√ D Y (2.21) so that: √ D ∂ ∂Y = ∂ ∂Z and Y = √ DZ (2.22) Note that, since D is a diagonal matrix, the square root of this matrix is another diagonal matrix with elements that are the square root of original matrix elements. Finally we have: 35 2.1. Two Energy Levels per Well ∂ ∂Z T ∂ ∂Z Ck + Z T √ DUAU−1 √ DZCk = ECk (2.23) We call √ DUAU−1 √ D matrix as S, where its elements scale like 1. At this stage, we should diagonalize S. We again use another orthonormal transformation: S = WD̃W−1 (2.24) where D̃ is the diagonal form of S and W is the matrix, formed by eigenvec- tors of S. Similarly, D̃ and W elements scale like 1. After substituting above transformation, the differential equation has the following form: ∂ ∂T T ∂ ∂T Ck +T T D̃TCk = ECk (2.25) where T = W−1Z. It sould be noticed that the whole procedure to reach T from X contains following change of variables: T = W−1Z = W−1 1√ D Y = W−1 1√ D UX ≡ ΓX (2.26) where Γ elements scale proportional to 1/ √ N . Note that Eq.(2.25) looks like the Harmonic Oscillator equation. Let’s compare this equation with SHO Schroedinger equation: p2 2m ψ + 1 2 mω2x2ψ = Eψ (2.27) where its ground state is a Gaussian function: ψ0(x) = ( mω pih̄ ) 1 4 e− x2 2h̄/mω ∝ e− x 2 2σ2 (2.28) where σ is the standard deviation of the probability distribution function. By comparing the kinetic term in Eq.(2.27) and derivative term in Eq.(2.25), we find that the effective mass for our system scales proportional to 1. Let me call the diagonal elements of D̃ matrix as λ1, λ2, λ3, where they are all in order of 1. Now, by comparing the potential term in Eq.(2.27) and the 36 2.1. Two Energy Levels per Well second term in Eq.(2.25), we find that ω2 should be proportional to λi. So, ωi is in order of 1. Notice that, the ground state of our problem is also a Gaussian func- tion in terms of elements of T matrix. Since, the standard deviation from Eq.(2.28) is proportional to √ 1/mω, the corresponding standard deviation for T elements, ∆ti, scales proportional to 1. Let’s consider T elements as t1, t2, t3. The relation between ti and ki and m could be written by using Eq.(2.26): t1 = α1m t2 = α2k1 + α3k2 t3 = α4k1 + α5k2 (2.29) where αi are elements of Γ matrix and are all in order of 1/ √ N . It should be pointed out that this special shape of transformation comes from the fact that there is no term proportional to mki or ∂m∂ki in Eq.(2.13). On the other hand, there is no term that mixes ki with m. Since standard deviation of T elements, ∆ti are in order of 1, as a result, ∆m ∝ ∆t1 1α1 ∝ √ N . Similarly, ki could be written in terms of t2 and t3 and αi. So, ∆ki also scale proportional to √ N . Note that the scaling is consistent with our early assumption for how ki and m scale with N . The behavior of NL −NR with N could be found by having this point in your mind that NL −NR = 4(k1 + k2). So, we could write: ∆(NL −NR)2 ∝ ∆(k1 + k2)2 = ∆k21 + ∆k22 + 2 < k1k2 > −2 < k1 >< k2 > (2.30) Thus, since ∆ki and ki scale like √ N , ∆(NL − NR)2 scales like N . So, σk1+k2 scales like √ N . This is the quantity that could be measured in exper- iments. As we mentioned in this section and previous chapter, the important point of our calculations is considering a large number of particles. But, we know that by increasing the interaction strength and the number of particles, more 37 2.1. Two Energy Levels per Well layers of energy would be occupied. Now, the question is that for what max- imum number of particles we can still apply ”two mode approximation”. On the other hand, what should be the range of the number of particles to have only two modes occupied and still have enough particles to have Gaussian distribution. If you refer to the calculations we have done, for minimum number of particles, we could say that as far as we have N 1, our assumptions and approximations are valid. It is clear that for smaller number of particles, we would have deviation from Gaussian distribution. So, the larger number of particles, the less deviation from Gaussian distribution. I guess some num- ber of particles like few thousands of bosons is enough to satisfy the N 1 condition. To find the maximum number of particles, we could look at the change in interaction energy of the system when we add a particle to the system and compare it with the energy gap between energy levels. On the other hand, if particles have enough interaction energy compared to the energy gap, they could go to higher energy levels and in this case, the ”two mode approximation” is no longer valid. So, the maximum number of particles, could be found by looking at the criterion where the change in interaction energy is equal to energy gap between ground state and first excited state in each well. Let’s consider we have N particles in the ground state of each well and we want to add another particle. If we add this particle to ground state, the change in the interaction energy could be found from Hamiltonian (2.1) to be: κ11((N + 1) 2 −N2) ∼ 2κ11N (2.31) But, if we add the particle to first excited state, the change in the inter- action energy would be κ12N . So, our desired criterion would be: (2κ11 − κ12)N = 2 − 1 (2.32) Let’s do some simple calculation to find the maximum number of particles for different values of system parameters. κ11 and κ12 could be defined similar to Eq.(19) as: 38 2.1. Two Energy Levels per Well κ11 = 2piah̄2 m ∫ d3r|ψ1(~r)|4 κ12 = 2piah̄2 m ∫ d3r|ψ1(~r)|2|ψ2(~r)|2 (2.33) For simplicity, we consider the wave functions of simple harmonic motion for ground and excited states of each well. So, we have: ψ1(~r) = α⊥√ pi αz√ pi 1 2 e− α2⊥(x 2+y2) 2 e− α2zz 2 2 ψ2(~r) = α⊥√ pi ( 2αz√ pi ) 1 2 e− α2⊥(x 2+y2) 2 e− α2zz 2 2 αzz (2.34) where: αz/⊥ = √ mωz/⊥ h̄ (2.35) Note that in ψ2, the particle is in the ground state of transverse directions and in the first excited state of the longitudinal direction. By substituting Eqs.(2.34) in Eqs.(2.33) and taking the integral we would have: κ11 = 8h̄2a mpi α2perpαz × 0.156664 κ12 = 8h̄2a mpi α2perpαz × 0.078332 (2.36) By knowing that energy gap is equal to h̄ωz and plugging above equations into Eq.(2.32), the maximum number of particles for ”two mode approxima- tion” to be valid would be: Nmax = h̄ωz 8h̄2a mpi α2perpαz × 0.234996 (2.37) or after simplification: Nmax = pi 8a× 0.234996 ωz ωperp √ h̄ mωz (2.38) 39 2.2. Finite Number of Energy Levels per Well Thus, by knowing the frequencies of the trap and scattering length and mass of the substance we use, we can estimate the maximum number of par- ticles and if it was also large enough to satisfy N 1 condition, then ”two mode approximation” could be applied for that problem. It is good to find this number for the experiment done by Joseph Thy- wissen in Toronto University [5]. In that experiment, ωz = 15Hz and ωperp = 1.38kHz. They have used Rb atoms with m ' 1.42 × 10−25kg. The typical scattering length for this kind of experiment is a = 100a0 where a0 is the Bohr atom radius that is a0 = 5.2917×10−11m. For these numbers, Nmax would be around 25 that is very small number. So, it means that we can not apply ”two mode approximation” for this experiment. This con- clusion makes sense because, the chemical potential for Thywissen group’s experiment is µ/h = 600Hz. Thus, by considering ωz = 15Hz, roughly we could say that about 40 energy levels are occupied in each well and is so beyond ”two mode approximation” limit. 2.2 Finite Number of Energy Levels per Well In this section, I will generalize what I did in section(2.1). Here, we consider n energy levels per well, where n is a finite number and n N . So, totally we have 2n energy levels. Our Hamiltonian in this case has the following form: H = n∑ i=1 i(niL + niR) + n∑ i,j=1 [δij(a † iLajR) + h.c.] n∑ i,j,k,l ∆ijkl [ (a†iLa † jLakLalL + h.c.) + (a † iRa † jRakRalR + h.c.) ] (2.39) where again in large N limit, i and δij scale like 1 and ∆ijkl scales propor- tional to 1/N . In this case, we have 2n kinds of creation operators for n levels in left and right wells, where all of them commute with each other. We have still 40 2.2. Finite Number of Energy Levels per Well the conservation of total number of particles, N , that put a constraint on creation operators. So, in this case, we need 2n − 1 quantum numbers to describe the system. Similar to Eq.(2.2), we define ki as the number of particles imbalance for one energy level in two wells: Ni,L,R = N n ± ki (2.40) where i goes from 1 to n in above equation. Note that, we have assumed that the most probable case is when each level has N/n bosons. The other n − 1 quantum numbers, could be defined similar to Eq.(2.4) as: Mi = Ni −Ni+1 2 (2.41) where i goes from 1 to n−1. In fact, Mi is the number of particles imbalance for two consecutive levels. By knowing these quantum numbers, the complete set of states for span- ning the Hilbert space would be: |k1, · · · , kn,M1, · · · ,Mn−1 > (2.42) For each Mi, similar to Eq.(2.3), we consider an average value, < Mi > and the deviation from this number as: mi = Mi− < Mi > (2.43) The most important part of this generalization is writing the number of particles in each level, Ni, in terms of Mi and N . From Eq.(2.41) we have: n−1∑ i=1 Mi = N1 −Nn 2 (2.44) so, N1 = Nn + 2 n−1∑ i=1 Mi (2.45) 41 2.2. Finite Number of Energy Levels per Well Similarly, for N2 we have: N2 = Nn + 2 n−1∑ i=2 Mi (2.46) Thus, in general we can write: Nj = Nn + 2 n−1∑ i=j Mi (2.47) where j goes from 1 to n− 1. We can also show the conservation of total number of particles as: n∑ j=1 = N (2.48) Plugging Eq.(2.47) into above equation, will give us: N = nNn + 2 n−1∑ k=1 n−1∑ i=k Mi (2.49) So, we can find Nn from above equation: Nn = 1 n (N − 2 n−1∑ k=1 n−1∑ i=k Mi) (2.50) After substituting above equation in Eq.(2.47), we have our desired rela- tion: Nj = 1 n (N − 2 n−1∑ k=1 n−1∑ i=k Mi) + 2 n−1∑ i=j Mi (2.51) where again i goes from 1 to n− 1. Above equation means that we can write Nj as a linear combination of Mi with some coefficients that they are at most in order of n 2 that is still so much smaller than N . By substituting Eq.(2.43) in above equation, it could be written as: 42 2.2. Finite Number of Energy Levels per Well Nj = αj + n−1∑ i=1 γijMi (2.52) where αj’s are some constants that include N and < Mi >. With similar justification as section(2.1), < Mi > scales proportional to N . γij’s are co- efficients of linear combination and are finite numbers and |γij| N . In analogy with section(2.1), we assume ki and mi scale proportional to √ N . Note that, we can recover Eqs.(2.6) and (2.7) by putting n = 2 in above equations. The action of Hamiltonian(2.39) on state (2.42) would be an equation similar to Eq.(2.10) where in square roots we have terms in αj and Mi. Thus, we can expand all the square roots and keep only terms in order of 1. In addition, we can also Taylor expand the obtained states around state (2.42) similar to Eq.(2.12). Since all the parameters and variables scale similar to n = 2 state, the time-independent Schroedinger equation for finite n should be something similar to Eq.(2.17). But this time, A and B are (2n−1)× (2n−1) matrices and X is a (2n − 1)-components vector. We can use the same orthonormal transformations and the same change of variables to achieve a Schroedinger equation like Eq.(2.25) for 2n− 1 Harmonic Oscillators. Finally, with the same justification, the standard deviation of ki and Mi scale proportional to √ N . It should be noticed that for all of steps in above calculations, n N . Note that for n ∼ N or continuum limit, we can not expand square roots in coefficients, because there is no guaranty that for all the γij’s, |γij| N . So, this method could not be used for continuum limit and we should use another techniques to find the distribution function of particles in that regime. I also would like to emphasize this point that even here in the presence of several variables for our distribution function, we can still use the central limit theorem. On the other hand, according to this theory, in the large N limit, the distribution function of different variables would be Gaussian that is consistent with the results we have found. Another point is that since we have assumed some special scaling of variables withN in all the parts of the problem, we ended up with Gaussian distribution function for number of particles and as I have mentioned earlier, this assumptions may 43 2.2. Finite Number of Energy Levels per Well not always be true and we chose them to have a nice theoretical behaviour. 44 Bibliography [1] Y. Shin,Experiments with Bose-Einstein Condensates in a Double Well Potential, Phd thesis,2000 cua.mit.edu/ketterle_group/Theses/PhD_Thesis_Y_Shin.pdf [2] Y. Shin, M. Saba,W. Ketterle .et al, Atom Interferometry with Bose- Einstein Condensates in a Double-Well Potential, Phys Rev. Lett Vol. 92,2004. [3] For example you can look at :Decoherence Dynamics in Low- Dimensional Cold Atom Interferometers, A. Burkov, M. Lukin and E. Demler, Phys. Rev. Lett. 98, 200404, 2007; Intrinsic Dephasing in One Dimensional Ultracold Atom Interferometers, R. Bistritzer, E. Altman, cond-mat/0609047v1, 2006;Dynamic Splitting of a Bose-Einstein Con- densate, C. Menotti, J. R . Anglin, J. I. Cirac and P. Zoller, Phys. Rev. A, Vol. 63, 023601, 2001. [4] J. Javanainen and M. Ivanov, Splitting a Trap Containing a Bose- Einstein Condensate- Atom umber Fluctuation,Phys. Rev. A, Vol. 60, 1999. [5] J. Thywissen, Number squeezing a Bose Josephson Junction, Unpub- lished, 2007. [6] B. Paredes, I. Bloch .et al. Tonk-Girardeau Gas of Ultracold Atoms in an Optical Lattice,Nature, 2004. [7] F. Dalfovo, S. Giorgin, L. Pitaevski, S. Stringari, Theory of Bose- Einstein Condensation in Trapped gases, Rev. of Mod. Phys. 71,No. 3,1999. [8] M. Olshanii, Atomic Scattering in the Presence of an External Con- finement and a Gas of Impenetrable Bosons, Phys, Rev, Lett. Vol. 81, 1998. 45 Bibliography [9] S. Giovanazzi, A. Gorlitz and T. Pfau, Tunning the Dipolar Interaction in Quantum Gases, Phys. Rev. Lett, Vol. 89, 2002. 46
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Splitting a Bose condensate
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Splitting a Bose condensate Mashayekhi, Mohammadsadegh 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Splitting a Bose condensate |
Creator |
Mashayekhi, Mohammadsadegh |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | The primary goal of this thesis is to find the distribution function of number of particles imbalance for ultra cold bosons trapped in a double well potential elongated in one direction for different strength of interaction between particles. This distribution function has been found to be Gaussian distribution function in two different limits. The first limit is weak interaction limit, where we only consider one energy level per well that is called "two mode approximation" regime. The second limit is when the interaction energy is in the same order as gap between energy levels where in this case, we consider finite number of levels per well. The standard deviation of number of particles distribution also has been found in both limits to be proportional to the square root of the total number of particles. In addition, some numerical work is done to find these distribution functions and the results are in agreement with analytical results. |
Extent | 712848 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-09-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0067676 |
URI | http://hdl.handle.net/2429/12647 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_fall_mashayekhi_mohammadsadegh.pdf [ 696.14kB ]
- Metadata
- JSON: 24-1.0067676.json
- JSON-LD: 24-1.0067676-ld.json
- RDF/XML (Pretty): 24-1.0067676-rdf.xml
- RDF/JSON: 24-1.0067676-rdf.json
- Turtle: 24-1.0067676-turtle.txt
- N-Triples: 24-1.0067676-rdf-ntriples.txt
- Original Record: 24-1.0067676-source.json
- Full Text
- 24-1.0067676-fulltext.txt
- Citation
- 24-1.0067676.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0067676/manifest