University of Southern California Dissertations and Theses

Properties of hard-core bosons in potential traps

(USC Thesis Other)

Properties of hard-core bosons in potential traps

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content PROPERTIES OF HARD-CORE BOSONS IN POTENTIAL TRAPS by Aditya Raghavan A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) August 2009 Copyright 2009 Aditya Raghavan Dedication I dedicate this thesis to my family and closest friends. My work would not be possible without their support. ii Acknowlegements In the years spent in Los Angeles, I have had the opportunity to interact with several individuals, which has deeply influenced me and shaped me as a person. First and foremost, my deepest appreciation goes to my thesis advisor, Stephan Haas for guiding me through these years of my PhD. He has gone out of his way to help me, on several occasions, making my journey as smooth as possible. Stephan has helped me learn the importance of every stage of research, from the small, toy models, to tech- niques of large scale simulations. His knowledge, intelligence, intuition and friendly demeanor will have a lasting impression on me. I would also like to thank Marcos Rigol, Pinaki Sengupta, Tommaso Roscilde and the late Kazumi Maki, for helping me in invaluable ways. I have worked closely with each of them during different times of my PhD and have learned various aspects of research from them. Through the years of my PhD, it has been my pleasure to have insightful discussions with fellow graduate students and post doctoral fellows. For this, I thank Kaushik Roy Choudhury, Ramakrishnan Iyer, Arnab Kundu, Amy Cassidy, Joyita Dutta, Manoj Gopalkrishnan, Noah Bray-Ali, Wen Zhang, Letian Ding and Alioscia Hamma. A special thanks to Amy for also helping me in my thesis preparation. iii Apart from my guide, I would like to thank Gene Bickers, Werner D¨ appen, Vitaly Kresin, Daniel Lidar and Chunming Wang who have been on my qualifying and disser- tation committees. A big thank you to Aiichiro Nakano for helping me with postdoctoral applications. My life would not have been the same without the love and support of my closest friends. I thank Sachin Telang for being a constant factor in my life through the years I have spent in Los Angeles. Finally, I would like to express immense gratitude and love to my father, mother and brother for their understanding and support. iv Table of Contents Dedication ii Acknowlegements iii List of Figures vii Abstract xii Chapter 1: Introduction 1 1.1 One-dimensional traps and hard-core bosons . . . . . . . . . . . . . 9 Chapter 2: Exact numerical technique to study the properties of one dimensional hard-core bosons 13 2.1 The Bose-Hubbard model and hard-core bosons . . . . . . . . . . . 14 2.2 The Jordan-Wigner transformation and the ground state of one dimen- sional hard-core bosons . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Numerical results for N=200 lattice sites . . . . . . . . . . . . . . . 20 2.4 Extension to other confinements . . . . . . . . . . . . . . . . . . . 28 Chapter 3: Time evolution of one dimensional hard-core bosons 33 3.1 Equal-time Green’s functions and time evolution . . . . . . . . . . . 34 3.2 Exact dynamics approach to study fundamental modes of a gas of one-dimensional hard-core bosons . . . . . . . . . . . . . . . . . . 35 3.3 Time-dependent traps: Splitting the Hamiltonian . . . . . . . . . . . 39 3.4 The Driven oscillator . . . . . . . . . . . . . . . . . . . . . . . . . 44 Chapter 4: Parametric excitations of hard-core bosons in a time-varying optical lattice 49 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2 Formalism to study time varying traps . . . . . . . . . . . . . . . . 50 4.3 Results forN = 50 lattice sites . . . . . . . . . . . . . . . . . . . . 52 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 v Chapter 5: The effects of disorder on one dimensional bosons in optical lattices 62 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2 Theoretical formulation of the problem . . . . . . . . . . . . . . . . 63 5.3 Quantum Monte Carlo results . . . . . . . . . . . . . . . . . . . . . 66 5.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Chapter 6: Conclusions 72 Bibliography 74 Chapter A: Eigenvalues of a Hamiltonian, H with hopping, t = 1 and trapping potential,V(x) = 0 79 vi List of Figures 2.1 Real space densities ofN b = 40(a),80(b),140(c) and corresponding momentum distribution functions forN b = 40(d),80(e),140(f) with number of lattice sites, N = 200. The densities are centered at the minimum of trap, i = 100. As the number of bosons increase, the real space density approaches pure Fock state occupation, like the one seen in (e). Here, part of the system is in a Mott insulating phae where the occupation is exactly 1 boson per lattice site over roughly 100 lattice sites. The corresponding momentum distribution functions show an increasing peak atk = 0 in the superfluid regime and a sharp decline when the Mott insulator is formed. . . . . . . . 22 2.2 Densities plotted for three system sizes,N = 100,200,500. In each case,V 0 is appropriately scaled and the densities appear to be iden- tical. However, for small lattices, one can see finite-size effects in the density profile. For larger values of N, these finite-size effects are smoothed out. . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3 Quasiparticle spectra plotted for N = 200 lattice sites. In Fig. (a), the, A = 9.0 and for large n, the spectrum is dominated by the potential energy. In Fig. (b), A = 3.0. The Mott plateau is created at a much larger value ofn and the potential energy never dominates the eigenvalue spectrum. . . . . . . . . . . . . . . . . . . . . . . . 24 2.4 The off-diagonal elements of the one-particle density matrix plotted here forN = 200 andN b = 40,120. A = 9.0 and forN b = 120 the system is in a Mott regime. The straight line fits are forx − 1 2 . As one can see, correlation between bosons follows a power-law decay. . . . 24 2.5 The difference between nearest energy levels, also called the deriva- tive spectrum is shown forA = 9.0 andN = 200. BeyondN ≈ 85, the system goes into a Mott regime which is signified by the branch- ing of the energy spectrum. In the continuum limit, this leads to a gapped system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 vii 2.6 The eigenvectors plotted for n = 1,6,120,121. For n = 1 we see a gassian-like eigenvector, much like the quantum oscillator eigen- function. For n = 6 we again see harmonic oscillator-type wave- function. For large n, the symmetry is broken by the presence of a Mott plateau. In this case, the wavefunctions are localized to the right or the left of the Mott plateau. . . . . . . . . . . . . . . . . . 27 2.7 The phase diagram for 1D HCBs. The white region is a full super- fluid phase while the black region is a fully Mott insulating phase. The Mott phase is usually in co-existence with superfluid phases on either side. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.8 Various even confining potentials for N = 200 lattice sites and strength, A = 4.0. In the limit α → ∞, we get the infinte square well potential. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.9 Eigenvalue spectra for different values of α. As α increases, the transition from linear to quadratic takes place at higher n. In the limit of large α one recovers the infinite square well where the energy follows the tight binding energy. . . . . . . . . . . . . . . . 30 2.10 Gaussian confinement plotted with a harmonic confinement, forN = 500 lattice sites with strength A = 9.0. The parameters a,b are explained in the text. . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.11 Comparison of eigenvalue spectrum of the gaussian confinement with a harmonic confinement of similar strengths of A = 9.0 with N = 500 lattice sites. In Fig. (a) we see that the spectra are same for low energies. The main difference comes at higher energies. The bounded nature of the gaussian, reduces the overall “height” of the energy spectrum as opposed to that of the harmonic confine- ment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1 Figure (a) is the oscillations of the peak of the momentum distribu- tion function,n(k) max . Figure (b) is the center of mass oscillations, x cm . In both cases N b = 20 HCBs and N = 100 lattice sites were used. In figure (a), the trap is suddenly relaxed with V 0 → 0.9V 0 , exhibiting a “breathing mode” of vibration. In figure (b), the trap is suddenly shifted to the right byδi = 5, exhibiting a dipole oscilla- tion (sloshing mode). . . . . . . . . . . . . . . . . . . . . . . . . . 37 viii 3.2 Frequency spectra plots of the center of mass oscillations of a sud- den shift in trap and the zero ofn(k) when the trap strength is sud- denly changed. The frequency spectra show dominant contributions atω S = 2 √ V 0 andω B = 4 √ V 0 . . . . . . . . . . . . . . . . . . . . . 37 3.3 The propagation of error forN b = 20 bosons on aN = 100 site lat- tice using a first order Trotter decomposition. The error is calculated as the absolute difference between expectation values of the number of bosons for an exact evolution and an approximate evolution when the trap curvature is suddenly changed fromV 0 → 0.9V 0 . It is clear from the grap that the error is linear in total timeτ. . . . . . . . . . 43 3.4 Enhanced oscillations of the center of mass motion of N b = 20 bosons onN = 100 lattice sites. The blue curve represents collec- tive oscillations when the trap is quenched or shifted to the right. The oscillations certainly do not die down in the 10 time periods in which they were measured. Here,V 0 = 9.0/(N/2) 2 = 0.0036 . . . . 46 3.5 The discrete Fourier transform of the center of mass motion ofN b = 20 bosons on N = 100 lattice sites. The frequency spectrum is dominated by two nearby frequencies close to the driving frequency, q 0 =ω 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1 Density plots of 15 hard-core bosons trapped in a N = 50-site lat- tice, plotted at timesτ = 0 and at τ = 10t/~. At τ = 0 the HCBs are in a superfluid phase with the spatial confinement given by a har- monic oscillator potential withV 0 a 2 = 4.0×10 −3 t. Forτ > 0, the curvature of the potential is sinusoidally modulated with a frequency q 0 = 1.26×10 −1 h/t and amplitudeΔV 0 = 0.1V 0 . . . . . . . . . . 50 4.2 Momentum distribution functions for at various times for a system ofN b = 15 bosons on a lattice withN = 50 sites. The momentum distribution function is characteristically peaked atk = 0. The three plots are shifted to resolve the peak heights. For small momenta, one can see the variation in the inset figure . . . . . . . . . . . . . . 51 4.3 Several plots of n(k = 0,τ) are shown here. The upper plot is for the superfluid case. The time series are divided by N b and are vertically adjusted to be well separated from each other. The scale at which oscillations take place in the Mott regime are drastically smaller. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 ix 4.4 The normalized frequency dependence of the zero-momentum peak of the momentum distribution function is shown for various total numbers of bosons N b in the trap. In (a) the system is initially in a superfluid phase, and the normalized dependence is an order of magnitude larger than in (b), where the system is initially in a Mott insulator regime. In (a), one observes a sharp peak atω 0 and satellite peaks close to 2ω 0 and 3ω 0 . These satellite peaks are shifted from 2ω 0 and 3ω 0 with increasingN b . ω 15 represents one such shift from 2ω 0 for N b = 15. In (b), the system is in a Mott insulator, and one observes a sharp resonance atω 0 with very small contributions elsewhere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.5 Contour map of the normalized frequency dependence of the zero- momentum peak of the momentum distribution functionn(k = 0,ω) for all values ofN b , the number of bosons. The contour color is in logscale and is quite noisly with several harmonics have small con- tributions. In the Mott regime, whenN b > 37, it is difficult to extract useful information from this contour plot. . . . . . . . . . . . . . . 56 4.6 (a) Contour map of the amplitude of the normalized momentum distribution function n(k = 0,ω) for all values of N b , the num- ber of bosons. (b) Differences between the effective single particle energy levels of the diagonalized Hamiltionian at τ = 0. The dif- ferent curves in (b) represent different inter-level transitions,Δε i = ε(N b +i)−ε(N b ). The contours in (a) of the excitation spectrum ofn(k = 0,ω) match the transition energies in (b). ForN b ≥ 37 a Mott plateau is present in the system. (a) shows that the system con- tinues to follow multi-level transitions in the Mott regime, however with an order of magnitude decrease in amplitude. . . . . . . . . . . 58 4.7 The dependence of the maximum amplitudes of|n(k = 0,ω)| are shown as a function of N b , using a logarithmic scale. In the super- fluid regime, the amplitudes increase linearly as a function of N b . However, in the Mott regime, one observes that the maximum ampli- tude is drastically reduced as a function of N b . The effect of this drastic decrease is seen in the scale of the contours in the Mott regime. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.1 Illustration of atoms in a potential trap. In the absence of interac- tions between the atoms, disorder leads to an immediate localization of the single-particle wave functions. . . . . . . . . . . . . . . . . 63 x 5.2 Effects of disorder on the momentum distribution n(q) and on the real-space density profilen(x) in the weak (upper panels) and strong (lower panels) coupling limit. The density profiles are vertically offset to enhance clarity of presentation. . . . . . . . . . . . . . . . 67 5.3 Dependence of the zero-momentum peak height, n(q = 0), on the disorder strength W for different values of U/t. At U/t = 5, the ground state consists of a single SF domain, whereas for the other values, there exist a finite MI region at the center of the trap. . . . . 69 5.4 Dependence of the momentum distribution function at q = 0, the stiffness and the compressibility of SF (left panels) and MI (right panels) phases on disorder strength. For the SF phase, the ground state is predominantly SF at small to moderate disorder and a Bose glass at large disorder. For the MI phase, the ground state is a Bose glass at all strengths of disorder studied. The strength of the zero-momentum peak varies non-monotonically, but the ground state always remains a Bose glass. . . . . . . . . . . . . . . . . . . 70 A.1 The convergence of the smallest and largest eigenvalues converging very quickly to ∓2. Here we have concentrated only on data for N ≤ 50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 xi Abstract Following the recent advances in controlling ultracold quantum gases that have led to the realization of boson and fermion condensation, we present computational studies of one-dimensional bosons on optical lattices. These systems have proven to be a great tool to study exciting new frontiers of condensed matter, from exotic phase transitions to nonequilibrium phenomena. In this thesis we study one-dimesional hard-core bosons, which is a limiting case of the Bose-Hubbard model when the on-site repulsion becomes infinite. Previously, this model has been well studied and is known to admit superfluid and Mott insulating phases. The time evolution of hard-core bosons in a one-dimensional trap is studied. Using an exact numerical approach, we study the ratio of the breathing mode and slosh- ing mode vibrations of the one-dimensional gas. We then study the evolution in a para- metrically driven trap. Using an exact numerical approach, the dynamics of the system is determined as the trap curvature is modulated. The response is found to be markedly different in the superfluid and Mott-insulating regimes. By measuring the frequency dependence of the zero-momentum peak of the momentum distribution function, para- metric excitations are observed, which depend on the phases present in the system. It is shown that these excitations closely match the quasi-particle spectrum, thus providing a useful tool to probe the low energy excitations of the system in its various phases. xii We then revisit the one-dimensional Bose-Hubbard model with a finite on-site repul- sion. We study the effects of disorder on interacting trapped bosons. At small to mod- erate disorder strengths it is observed that if there is a Mott plateau at the center of the trap in the clean limit, phase coherence increases as a result of turning on disorder. The localization effects due to correlation and disorder compete against each other, resulting in a partial delocalization of the particles in the Mott region. We show that in this regime delocalization can lead to an increase in phase coherence. xiii Chapter 1 Introduction Ensembles of trapped atoms are known to have a rich phenomenology, includ- ing Bose and Fermi condensation, Feschbach resonances, and superfluid-to-insulator transitions(22; 23; 38; 57; 60; 46; 47; 48). Recent experimental realizations of these systems have lead to valuable insight into correlated many-body systems. Work in this area can be dated back to the Nobel prize winning experiments that led to the first ever seen Bose-Einstein condensates (BECs) of dilute akali atoms(1; 6; 9). Subsequently, the production of artificial crystals of light(22; 23; 28; 41), consisting of several thousands of atomic microtraps with the help of counter-propagating laser beams continued to stir the excitement of experimentalists and theorists. These lattices in 1, 2 and 3 dimensions have now become a testing ground for strongly correlated quantum phases. Strong correlation effects in these dilute quantum gases can be easily studied, even though the density of these gases is several orders of magnitude smaller than their traditional condensed matter counterparts. The idea that these experiments produce clean and pure degenerate quantum systems leads one to think that optical lat- tice systems can be used as quantum simulators which could simulate the dynamics of other, complex and exotic quantum Hamiltonians. In theory, this idea was first proposed by Richard P. Feynman(15), and has since been suggested by Immanuel Bloch(4) as a future avenue of research in optical lattice systems. Apart from strongly correlated many-body systems, these experimental techniques have seen direct applications in quantum information processing(5). Spin-dependent lattice potentials allow for atoms in one internal state to experience a potential that is 1 different from those atoms in another internal state. In this manner, one could create two distinct states (two-state system) and allow for controlled contact between atoms in different states to produce a coherent superposition of these states. Describing bosonic systems from the theoretical side, the superfluid to Mott insulator transition has been studied extensively within the Bose-Hubbard model(60; 46; 47; 48). For soft-core bosons, these studies were performed using Quantum Monte Carlo(60) and density matrix renormalization group (DMRG) (29; 30) techniques. As a common element to these studies, it was found that beyond a sufficiently high density of bosons, determined by the number of bosons in the system and by the curvature of the trap poten- tial, a “Mott plateau” starts to form at the center of the trap. As the density is increased, the Mott plateau expands towards the edges of the trap, and additional plateaus appear at the trap center. Furthermore, it was shown that this finite-size rendition of a quan- tum phase transition can also be triggered by tuning the relative strength of the hopping parameter,t and the on-site repulsionU. For the one-dimensional case, the properties of hard-core bosons(46; 47; 48) in the presence of a trapping potential can be studied for large system sizes. Specifically, using the Jordan-Wigner transformation, the density profiles, the one-particle density matrix and the momentum density profiles can be calculated exactly. In this case, since U/t→∞, the superfluid to Mott insulator transition is only controlled by the strength of the trapping potential and the number of bosons. More recently, there has been a surge of interest in the non-equilibrium dynamics of one-dimensional bosons(48; 29). For hard-core bosons, non-equilibrium properties have been studied for the case when the trap is suddenly shifted(48). It was found that the resulting collective oscillations are strongly suppressed in the Mott regime, indicating that the bosons trapped in the Mott plateau do not participate in the collective modes(48). 2 The extremely low temperatures required for Bose-Einstein condensation (BEC) and ultracold atom experiments can be achieved by a combination of cooling techniques. Firstly, spin-polarized atoms are trapped in a magneto-optical trap using laser cooling. The cooling is produced by illuminating the atoms with laser beams in six orthogonal directions. The laser beams are slightly red-detuned from a natural absorption line of the atoms, so that fast moving atoms that move in the direction of one of these laser beams, resonantly absorb due to the Doppler effect and subsequently slow down. The emission of the absorbed photon occurs in any direction, and for a sufficiently large sample the change in momentum is averaged to zero. However, the recoil energy of the emitted photon limits the reduction to temperatureμK scales. In order to producenK tempera- tures, evaporative cooling is used, in which the atoms are transferred into a magnetic trap and the hottest atoms are evaporated to allow the remaining atoms to thermalize thus bringing down the overall temperature to nK. The typical temperatures and densities at which BECs are observed are usually such that an equilibrium configuration of the system would actually be in a solid phase. However, in the ultracold gases experiments, the gaseous phase is preserved by keeping the system very dilute, restricting the atoms to two-body collisions and making three-body and higher collisions much rarer events. Usually BECs produced by this method are in a spherically symmetric potential. In order to form a three-dimensional lattice, three optical standing waves are aligned orthogonal to each other, with their mutual crossing point exactly over the BEC. The standing waves are produced by focussing the laser beam to a desired waist at the posi- tion of the condensate, with a second lens and a mirror on the other side of the conden- sate. Each laser beam is polarized in a direction orthogonal to the other two and the beams typically have slightly different frequencies with respect to each other, so that any interference between orthogonal beams would average out to zero. The beam waist along with the wavelength of the laser beam, determines the number of lattice sites in 3 each orthogonal direction. For each direction it is possible to produce about 100 lattice sites, thus reaching close to a million lattice sites in 3-D. The lattice spacing isa =λ/2 whereλ is the wavelength of the laser beam. The periodic intensity maxima and minima of the standing waves produce a periodic potential in which two-level atoms are attracted or repelled from maxima depending on the detuning of the laser light with respect to the atomic resonance. For red detuned lasers, atoms are attracted towards the intensity maxima while for blue detuned, the intensity minima. The resulting three dimensional optical potential 1 seen by the atoms is proportional to the sum of the intensities of the standing waves which leads to a simple cubic geometry V 0 (x,y,z) =V 0 (sin 2 (kx)+sin 2 (ky)+sin 2 (kz)), (1.1) here k = 2π/λ is the wave-vector of the laser light and V 0 is the maximum potential depth of a single standing wave laser field. A reasonable starting point to study optical lattices theoretically is to look at the con- tinuum many-body description. We write the Hamiltonian operator in terms of bosonic fields for atoms in a given internal state H = Z ˆ Ψ † (~ r) −~ 2 2m ▽ 2 +V 0 (~ r)+V T (~ r) ˆ Ψ(~ r) ~ dr + 1 2 Z ~ dr ~ dr ′ˆ Ψ † (~ r) ˆ Ψ † ( ~ r ′ )U(~ r− ~ r ′ ) ˆ Ψ( ~ r ′ ) ˆ Ψ(~ r), (1.2) where the bosonic fields have the commutation relations h ˆ Ψ(~ r), ˆ Ψ † (~ r) i =δ(~ r− ~ r ′ ), h ˆ Ψ(~ r), ˆ Ψ( ~ r ′ ) i = 0. (1.3) 1 For details on the derivation of optical dipole traps for neutral atoms see ref (24) 4 HereV 0 (~ r) is taken to be of the form in Eq. (1.1),V T (~ r) is a trapping potential that is typically a magnetic trap. The interaction term is restricted to two body interactions due to the very dilute nature of the gas. The energies of colliding atoms is very small and typically smaller than the centrifugal barrier imposed by non-zero angular momentum collisions. Thus, only s-wave collisions are relevant. It is then possible to replaceU(~ r− ~ r ′ ) by a pseudopotential which is a delta-function potential with strength proportional to the real s-wave scattering lengtha s of the atoms U(~ r− ~ r ′ ) = 4π~ 2 a s m δ(~ r− ~ r ′ ). (1.4) The Hamiltonian is then reduced to H = Z ˆ Ψ † (~ r) −~ 2 2m ▽ 2 +V 0 (~ r)+V T (~ r) ˆ Ψ(~ r) ~ dr + 1 2 4π~ 2 a s m Z ~ dr ˆ Ψ † (~ r) ˆ Ψ † (~ r) ˆ Ψ(~ r) ˆ Ψ(~ r). (1.5) The lattice-like periodicity leads to localized wave-functions and if we assume that the s-wave scattering energy is less than the gap between the first and second band in the lattice, we can expand the bosonic fields in the basis of the Wannier functions in the lowest vibrational band, ˆ Ψ(~ r) = X i b i φ i (~ r), ˆ Ψ † (~ r) = X i b † i φ ∗ i (~ r). (1.6) The Wannier functions are complete and orthonormal, X i φ i (~ r)φ ∗ i ( ~ r ′ ) =δ(~ r− ~ r ′ ), Z ~ drφ i (~ r)φ ∗ j (~ r) =δ ij , (1.7) 5 and the creation (b † i ) and annihilation 2 (b i ) operators at site i satisfy the commutation rules h b i ,b † j i =δ ij , [b i ,b j ] = 0, (1.8) which follows from Eq. (1.3). The Hamiltonian can be written in terms of the creation and annihilation operators H =− X ij (t ij +V ij )b † i b j + X ijkl U ijkl b † i b † j b k b l , (1.9) where t ij = − Z ~ drφ ∗ i (~ r) − ~ 2 2m ▽ 2 +V 0 (~ r) ~ dr φ j (~ r), (1.10) V ij = Z ~ dr φ ∗ i (~ r)V T (~ r)φ j (~ r), U ijkl = 2π~ 2 a s 2m Z ~ dr φ ∗ i (~ r) ~ dr φ ∗ j (~ r) ~ dr φ k (~ r) ~ dr φ l (~ r). We can further assume that the Wannier functions are localized over a few lattice sites and the only significant contributions to these sums come from nearest-neighbor hop- ping, on-site repulsion and an energy offset at each site due to the external trapping potential. The Hamiltonian reduces to the Bose-Hubbard model H =−t X {ij} b † i b j +h.c. + X i n i ε i + U 2 X i n i (n i −1) (1.11) 2 We choose the widely accepted symbolsb † i ,b i to represent bosonic creation and annihilation operators without the operator “hat” to avoid redundancy. In the remainder of the text, standard bosonic creation- annihilation operators, as well as number operators (n i ) will be represented without hats. 6 wheret = t {ij} , the nearest-neighbor hopping is assumed to be constant throughout the lattice and ε i = Z ~ dr V T (~ r)|φ(~ r i )| 2 (1.12) is a site-wise energy offset that depends on the functional dependence of V T (~ r). Note that this site-wise energy offset breaks the homogeneity of the system and leads to unique physics in optical lattice systems. Since the Wannier functions are localized about ~ r i , Eqn. (1.12) can be approximated by the value ofV T (~ r) at thei th lattice site, 3 V T (~ r i ). As mentioned earlier, the trapping potentialV T is usually a standard magnetic trap. David Pritchard(43) was the first to propose a magnetic trap to be used to trap very cold atoms. Although the traps have undergone sophistication over the years, the basic idea behind this kind of trap is to have a uniform magnetic field in the ˆ z-direction, which fixes the quantization axis. Superimposed with this uniform field, gradients are added in such a way as to produce a minimum in the field atz = 0, thus localizing “weak-field seeking” particles to the bottom of the trap. The original trap conceived by Pritchard(43) involved a complicated “bottle field” written here in cylindrical polar co-ordinates. ~ B b =B 0 1+ k 2 2 z 2 − ρ 2 2 ˆ z− B 0 2 k 2 zρˆ ρ. (1.13) At the origin of this “bottle”,z = ρ = 0, the field has a minimum along the ˆ z-axis but a maximum in theρ−φ-plane. An additional quadropole field ~ B q =b q ρ[cos(2φ)ˆ ρ−sin(2φ)ˆ ρ], (1.14) 3 In other words, V T (~ r) is slowly varying in the regime where φ(~ r) is finite. Hence, R ~ drV T (~ r)|φ(~ r)| 2 ≈ V T (~ r i ) R ~ dr|φ(~ r)| 2 = V T (~ r i ). V T (~ r) need not be slowly varying globally. The main assumption here is that V T (~ r) is not a very locally varying function. (Local meaning over one neighborhood of a lattice site.) 7 is added and b q ,k are adjusted so that a minimum is created at z = 0. Near the mini- mum the field is quadratic in nature (second derivative6= 0) leading to simple harmonic oscillator energy levels for the field. The disadvantage of having magnetic traps, however is that only a small subset of the available atomic spin states, the “weak-field-seeking” states, can be trapped. This limitation has been overcome by using optical dipole traps that use the interaction of the induced dipole moment in an atom and an external electric field. Today’s optical lattice experiments rely on the oscillating electric light field from a laser, which induces an oscillating dipole moment in the atom ( ~ d) while at the same time interacts with this dipole moment in order to create a trapping potential,V dip (~ r) for the atoms V dip (~ r) =− ~ d. ~ E(~ r) =α(ω)| ~ E(~ r)| 2 , (1.15) where α is the polarizability of the atom. Typically, the intensity, I(~ r) = | ~ E(~ r)| 2 of a laser beam is a gaussian profile which, for a broad enough beam-waist, provides a shallow harmonic oscillator type potential at the peak of the intensity 4 We can then write our Hamiltonian as H =−t X {ij} b † i b j +h.c. +V T X i n i (~ r i −~ r 0 ) 2 + U 2 X i n i (n i −1), (1.16) where~ r 0 is the center of the trap andV T is the strength of the trapping potential. For a given optical potential it is possibly to numerically compute the Hubbard coef- ficientsU,t as a function of the parameters of the potential,V 0 ,λ. The Wannier functions 4 For a laser beam travelling in the z-direction, we can write the intensity as I(ρ) = I 0 e −ρ 2 /ρ 2 0 where ρ 0 is the beam-waist. Then, V dip = α(ω)I 0 e −ρ 2 /ρ 2 0 ≈ α(ω)I 0 1−ρ 2 /ρ 2 0 +1/2(ρ/ρ 0 ) 4 +... . The beam waist is usually 100 times larger than a lattice spacing, so V dip (~ r) is a shallow harmonic oscillator potential. 8 can be written as products. φ(~ r) = φ(x)φ(y)φ(z) and one can solve the 1D band- structure for a cosine potential. However, since we are interested in the ground state (lowest band), qualitative insight into the dependence of these parameters is obtained in a harmonic approximation of the cosine potential byV 0 cos(kx)≈V 0 [1− 1 2 k 2 x 2 ] With the use of the ground-state gaussian wave-functions of the oscillator potential, one can estimate reasonably well, the Hubbard-t,U parameters as t = 4 π V 0 E r 3/4 E r e −2 q V 0 Er , U = r 8 π V 0 E r 3/4 k a s E r , (1.17) wherek = 2π/λ is the wave-vector of the laser light,a s is the s-wave scattering length andE r = ~ 2 k 2 /2m is the recoil energy of atoms absorbing light of wave-vectork and is a suitable energy scale used to normalizeV 0 . 1.1 One-dimensional traps and hard-core bosons Eq. (1.17) gives a convenient form of t and U parameters to experimentally study the Bose-Hubbard model. In this exposition, we do various studies of 1D hard-core bosons which is obtained from the 1D Bose-Hubbard model, when U/t → ∞. Through the technique explained above, it is not very clear how an experimental realization of this limit would be possible. For the purpose of better understanding this, we revisit the problem for the case of a “cigar” shaped atom trap as first discussed by M. Olshanii(40). The “cigar” shaped trap is effectively an axially symmetric 2D harmonic potential with natural frequencyω ⊥ , V T = 1 2 mω 2 ⊥ x 2 +y 2 . (1.18) 9 The atoms are cooled down to below the transverse vibrational energy~ω ⊥ resulting in motion being highly restricted in the transverse direction. The atoms move freely in the z-direction, which is the “effective” dimension of the problem. The two-body interaction termU(~ r) is given by U(~ r) =g δ(~ r) ∂ ∂r (r·) , (1.19) where~ r is now the relative co-ordinate andg = 2π~ 2 a s /μ is the strength,a s is the real s-wave scattering length for the atoms, μ = m/2 is the effective mass. The (r·) term implies that the partial derivative with respect tor is taken after doing the multiplication, r · Ψ(r). This removes the 1/r divergence in the scattered wave. The reason why this pseudopotential differs from Eqn.(1.4) is that the following calculation is not for Wannier functions in the presence of an optical lattice. The wave function in the z- direction is an unbounded traveling wave. The Hamiltonian for this model can be written as(40) H = p 2 z 2μ +g δ(~ r) ∂ ∂r (r·)+H ⊥ (p x ,p y ,x,y), (1.20) where H ⊥ is the transverse Hamiltonian and contains the V T term discussed in Eqn. (1.18). The incident wave ise ikzz φ 00 (ρ), whereφ 00 is the ground-state wave-function of the 2D harmonic oscillator. The scattered wave has an asymptotic form Ψ(z,ρ)−→ |z|→∞ {e ikzz +f even e ikz|z| +f odd sign(z)e ikz|z| }φ 00 (ρ), (1.21) wheref even andf odd are the 1D scattering amplitudes. 10 In order to calculate the scattering amplitudes, Ψ(z,ρ) is expanded in the basis of eigenstates of the transverse Hamiltonian and substituted in the Schr¨ odinger equation for the total Hamiltonian with the energy taken to be E = ~ 2 k 2 z 2μ +~ω ⊥ . (1.22) Along with the asymptotic conditions and the continuity of the wave function and its derivative, one can write the wave-function in terms of a phase shift Ψ∝ sin(k z |z|+Δ(k z ))φ 00 (ρ). (1.23) The 1D scattering amplitude is defined as−∂Δ/∂k z | kz→0 + and is calculated a 1D =− a ⊥ 2a s 1−C a s a ⊥ , (1.24) with the following expressions for the scattering amplitudes f even =− 1 1+ik z a 1D −O((k z a ⊥ ) 3 ) , f odd = 0. (1.25) Here, a ⊥ is related to the width of the ground-state wavefunctions of the transverse Hamiltonian and a s is the “true” 3-D s-wave scattering length for the atoms. The con- stant C = 1.4603... is the limit of a series that arrives in this calculation(40). It can then be seen that for low velocities, the scattering amplitude is approximated by f even (k z ) =−1/(1+ik z a 1D ) for a 1D delta function potential U 1D =g 1D δ(z), (1.26) 11 with coupling strength, following Eqn.(1.24) g 1D =g|φ 00 (ρ)| 2 1−C a a ⊥ −1 . (1.27) Once we have acquired this expression for the 1D scattering potential, we see that the case of g 1D → ∞ is achievable via this technique. g 1D depends on the 3-D scattering length as well as the width of the ground-state wave-function of the transverse Hamilto- nian. Both these quantities are adjustable via a suitable choice of atoms and transverse trapping frequencies. Experiments to study reduced dimensional bosons and fermions have become very common and, most notably first done by Goerlitz et. al,(21) and followed by Schreck et al,(52). Moritz et al. have measured the ratio of the breathing mode and sloshing mode of a quasi 1D gas of bosons and shown that it follows the universal ratio for 1D. Strong evidence of fermionization of 1D strongly interacting bosons has been seen in a few experiments(57)(41). 12 Chapter 2 Exact numerical technique to study the properties of one dimensional hard-core bosons Following recent advances in experiments, it is now possible to engineer very clean systems of strongly interacting bosons and fermions in a confined trap geometry with tunable parameters(22; 23; 38; 57). In particular, lattice hard-core bosons have recently been realized experimentally in tube-like quasi-1D confining potentials(41). Lot of work has been done in the subject of 1D hard-core bosons, including non-equilibrium studies(46; 47; 48). In this chapter, we present a study of the ground state properties of 1D hard-core bosons (HCBs) on optical lattices with an external confining potential. Without loss of generality, we focus on harmonic confinements which is the most popular case in recent theoretical and experimental endeavors of trapped ultracold gases(38; 57; 60; 46; 47). Following the work done by Rigol and Muramatsu(46), we review a technique to study exact ground state properties of hard-core bosons in 1D lattices. As we will show ahead, the fermionization of strongly repelling bosons, leads to the mapping of 1D HCBs to a simple, integrable model of non-interacting fermions. In the natural basis of the fermionic operators, the Hamiltonian is then shown to be diagonal in potential energy and off-diagonal in kinetic energy. We then discuss other confining potentials and also give a specific example of the attractive Gaussian potential which is similar to 13 modern-day trapped atom experiments. The differences that arise when one chooses a harmonic confinement as an approximation to the Gaussian confinement is discussed in conclusion. 2.1 The Bose-Hubbard model and hard-core bosons A good starting point for this subject is the Bose-Hubbard model which has been dis- cussed in Chapters 1,2. Eq. (2.1) is the Bose-Hubbard model on a one dimensional lattice of sizeN. H =−t N X i=1 b † i b i+1 +b † i+1 b i + U 2 N X i=1 n i (n i −1)+V 0 N X i=1 (x i −x 0 ) 2 n i , (2.1) wheret is the hopping parameter,U is the on-site repulsion andx 0 marks the center of the trap. x i =ai, wherea is the lattice spacing andn i =b † i b i , is the number operator. For an optical lattice, the hopping parameter and on-site repulsion are consant throughout the lattice 1 . The repulsion term can be understood as the energy expenditure in having two bosons at one lattice site. As U increases, the probability of having two bosons at a particular site becomes smaller and smaller. In theU →∞ limit, the probability of two bosons sitting on one lattice site becomes zero and we get hard-core bosons (HCBs) 2 . In 1 A detailed discussion on the experimental tunability of these parameters has been done in the Intro- duction 2 There are various names to describe these bosons. Some people say impenetrable bosons while others refer to the regime where inter-particle repulsion dominates 1D bosons as the Tonks-Girardeau regime. We will follow Rigol et. al(46; 47) among others, and use the term hard-core bosons for the remainder of the text. 14 this limit, the creation (or annihilation) of two bosons at a particular lattice site is taken to be zero. b †2 i = 0, b 2 i = 0. (2.2) This leads to on-site anti-commutation relations (like fermions) but off-site, the oper- ators continue to commute (like bosons). {b † i ,b i } = 1, h b † i ,b j i = 0,i6=j. (2.3) The Hamiltonian for 1D HCBs can then be written as H =−t N X i=1 b † i b i+1 +b † i+1 b i +V 0 N X i=1 (i−i 0 ) 2 n i , (2.4) where the strength of the trap,V 0 is measured in units of lattice spacing squared,V 0 a 2 → V 0 . 2.2 The Jordan-Wigner transformation and the ground state of one dimensional hard-core bosons In order to solve the Hamiltonian in Eq. (2.4) we use the Jordan-Wigner transformation which maps the system of 1D HCBs to one of spinless fermions(26) b † i =f † i i−1 Y k=1 e −iπf † k f k , (2.5) b i = i−1 Y k=1 e iπf † k f k f i , (2.6) 15 where f † i and f i are the creation-annihilation operators for spinless fermions. They satisfy the standard anti-commutation rules n f † i ,f j o =δ ij . (2.7) One can evaluationn i =b † i b i in the representation off ′ i s. b † i b i =f † i i−1 Y k=1 e −iπf † k f k . i−1 Y l=1 e iπf † l f l f i , =f † i f i . (2.8) n i = n f i = f † i f i and the total number of quasi-particles (fermions) is just equal to the total number of bosons. The Hamiltonian then becomes H =−t N X i=1 f † i f i+1 +f † i+1 f i +V 0 N X i=1 (i−i 0 ) 2 n f i . (2.9) The Hamiltonian now reduces to a simple, diagonalizable Hamiltonian for non- interacting, spinless, fermions in a harmonic trap. Hence, the spectrum for the HCBs is exactly identical to that of non-interacting fermions. The off-diagonal terms of the density matrix, D b † i b j E have a non-trivial relation with the fermionic representation and this lead to a significant difference in the momentum distribution function of 1D HCBs as opposed to 1D fermions. The momentum distribution function of the HCBs is sharply peaked atk = 0, which implies a multiple occupation ofk = 0 momentum state which is not allowed for fermions due to the Pauli exclusion principle. 16 Following (46), the one-particle Green’s function can be written in the basis of HCBs. Then, after applying the Jordan-Wigner transformation, we get G ij =hΨ 0 |b i b † j |Ψ 0 i =hΨ f 0 | i−1 Y k=1 e iπf † k f k f i f † j j−1 Y l=1 e −iπf † l f l |Ψ f 0 i, (2.10) where|Ψ f 0 i is the non-interacting fermionic ground-state. This can be obtained by diag- onalizing equation (2.9) and filling the firstN f number of energy levels with one fermion each (47). Thus, the ground-state is given by (46) |Ψ f 0 i = N f Y m=1 N X n=1 P nm f † n |0i, (2.11) where P nm is the matrix of the first N f eigenfunctions of the Hamiltonian. The total number of fermions, N f is simply equal to the total number of bosons, and it is clear from equation (2.11) that this quantity is a fixed quantity. Here, P is aNxN f matrix P =             P 11 P 12 . . . P 1N f P 21 P 22 . . . P 2N f . . . . . . . . . . . . P N1 P N2 . . . P NN f             , (2.12) which is produced by diagonalizing Eq. (2.9) and collecting the first N f eigenvectors. Hence, P = v 1 v 2 . . . . v N f , (2.13) 17 wherev i is thei th eigenvector of H. Now, Eq. (2.10) can be written in terms of the inner product of two transformed vectors|Ψ 1 i,|Ψ 2 i, where hΨ 1 | = f † i i−1 Y k=1 e −iπf † k f k |Ψ f 0 i ! † , |Ψ 2 i =f † j i−1 Y l=1 e −iπf † l f l |Ψ f 0 i, (2.14) and, G ij =hΨ 1 |Ψ 2 i. (2.15) In order to simplify Eq. (2.15) and Eq. (2.14) we notice that f † k f k is an idempotent matrix, so f † k f k n =f † k f k . Then, one can simplify the exponentiale −iπf † k f k as e −iπf † k f k = ∞ X j=0 (−iπ) j j! (f † k f k ) j , = 1+ (−iπ) 1! f † k f k + (−iπ) 2 2! f † k f k +... (−iπ) n n! f † k f k ..., = 1+ ∞ X j=1 (−iπ) j j! f † k f k = 1+ e −iπ −1 f † k f k = 1−2f † k f k . (2.16) This results in a simplification of Eq. (2.5) i−1 Y k=1 e −iπf † k f k = i−1 Y k=1 1+ e −iπ −1 f † k f k = i−1 Y k=1 1−2f † k f k . (2.17) Although the Jordan-Wigner transformation appeared to be complicated, it is important to comment that on simplification, the exponent reduces to a simple form in terms of f † k f k . Further, f † k f k = n f k = 1,0 and so, the exponent takes the values 3 ±1. Then, b † i =∓f † i depending on how many fermions are occupied before thei th lattice point. 3 One can reach this result by simply claiming thatf † k f k = 0,1 which implies thate −iπf † k f k =±1 18 Since the fermionic ground-state has N f particles, the action of Q i−1 k=1 e −iπf † k f k changes the sign of the elements of P nm , for which n ≤ i− 1 and for all values of m. On the other hand, f i acts onhΨ f 0 | to create a particle at the i th lattice point. This is equivalent to adding a new column to the matrix P nm with every element zero, but P iN f+1 = 1. An example fori = 4 is shown here 3 Y k=1 e −iπf † k f k P =                   −P 11 −P 12 . . . −P 1N f 0 −P 21 −P 22 . . . −P 2N f 0 −P 31 −P 32 . . . −P 3N f 0 P 41 P 42 . . . P 4N f 1 . . . . . . 0 . . . . . . 0 P N1 P N2 . . . P NN f 0                   . (2.18) The same can be done for the other product term, and the creation operator for the j th position. Then, we get a simple form for the Green’s function in terms of the transformed P matrices G ij = det n P † 1 P 2 o , (2.19) where P 1 and P 2 are the matrices formed from the action of the operators mentioned above. Finally, the one-particle density matrix is given by ρ ij =hb † i b j i, =hb j b † i i(1−δ ij )+δ ij (1−b i b † i ), =G ij (1−δ ij )+δ ij (1−G ii ), =G ij +δ ij (1−2G ii ). (2.20) 19 Experiments usually probe the momentum distribution function by measuring the veloc- ity of moving particles after the trap is suddenly released. This procedure is called a time of flight measurement. The resulting data is simply the momentum distribution function. One can study this function analytically and numerically. To do so, we expand the cre- ation (annihilation) operators in the plane-wave basis, withk taking2N values in the1 st Brillouin zone from [− π N , π N ] b i = 1 √ N X k e iki ˆ b k , (2.21) Then, one can write ˆ b k = 1 √ N X i e −iki b i , (2.22) One can define, the expectation valueh ˆ b † k ˆ b k i =n k =n(k). Then, n(k) = 1 N N X i,j=1 e ik(i−j) hb † i b j i, = 1 N N X i,j=1 e ik(i−j) ρ ij . (2.23) Here we have used the ˆ b k sign to signify annihilation (creation) operators in the plane- wave basis. 2.3 Numerical results for N=200 lattice sites We extend our study to N = 200. Although it is possible to diagonalize much larger matrices, the Green’s function calculation is computationally very costly. For each ele- ment of the Green’s function matrix, one needs to calculate the determinant of aN b xN b 20 matrix produced by the product of the P † 1 P 2 . The determinant algorithm has a compu- tation cost of N 3 b which must be multiplied by N 2 , since there are N 2 elements of the Green’s function matrix. This is the costliest part of the code. The computation time for the Green’s function matrix is then,T c ∝N 3 b N 2 . The Hamiltonian depends only on one parameter, namelyV 0 /t. The overall size of the lattice does not matter and the shape of the real-space and momentum densities are invariant under scaling transformations of the lattice size, N → 2N provided V 0 /t is appropriately scaled. We have observed that V 0 /t = 4A N 2 , (2.24) ensures invariance in real-space density of the HCBs. Here,A, the strength of the con- fining potential is suitably chosen so that the major contribution of energy at the edges of the trap is the potential energy, but near the center the energy is shared between the kinetic energy and the potential energy. At the edge of the trap, the lattice index takes i = N/2 and the potential energy is V(i = N/2)/t = A. Then, ifA∼ 10, the largest energy eigenvalue (E(n =N)) is dominated byA and in effect, E N ∼A. (2.25) So, if we had to occupy every energy level with one boson, the approximate energy needed to occupy the last energy level isA. The momentum distribution function (MDF),n(k) also has the same shape, but usu- ally the peak n(k = 0) scales as √ N(47). Of course, for very small lattices, N < 20 the scale invariance is still valid, but edge effects become significant. The edges are just 2 points, however for small lattices that is 10% of the total lattice. 21 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 n(x) 0 0.5 1 1.5 2 2.5 n(k) 0 50 100 150 200 x/a 0 0.2 0.4 0.6 0.8 1 100 200 300 400 ka 0.5 1 1.5 2 2.5 (a) (b) (c) (f) (e) (d) Figure 2.1: Real space densities of N b = 40(a),80(b),140(c) and corresponding momentum distribution functions for N b = 40(d),80(e),140(f) with number of lattice sites, N = 200. The densities are centered at the minimum of trap, i = 100. As the number of bosons increase, the real space density approaches pure Fock state occupa- tion, like the one seen in (e). Here, part of the system is in a Mott insulating phae where the occupation is exactly 1 boson per lattice site over roughly 100 lattice sites. The cor- responding momentum distribution functions show an increasing peak at k = 0 in the superfluid regime and a sharp decline when the Mott insulator is formed. The momentum distribution function is calculated using the one-particle density matrix,ρ ij . The diagonal elements of this matrix is simply the expectation valuehb † i b i i, which is just the real-space density. The off-diagonal elements represent two-point cor- relations in space. As we see in Fig. (2.3), the correlations between bosons follows a power-law behavior with an exponent of 1 2 . For the case of N b = 40, the system is in a superfluid phase and the correlations are about two times larger than for N b = 120. Also, the correlations die off suddenly in both cases. In the superfluid case, the real space density spreads across only about 100 lattice sites and the remainder of the lattice has occupation identically zero. In the Mott regime, pure Fock states ofn i = 1 cut-off 22 0 0.2 0.4 0.6 0.8 1 x/N 0 0.1 0.2 0.3 0.4 0.5 n(x) N b =20, N=100 N b =40, N=200 N b =100, N=500 Figure 2.2: Densities plotted for three system sizes, N = 100,200,500. In each case, V 0 is appropriately scaled and the densities appear to be identical. However, for small lattices, one can see finite-size effects in the density profile. For larger values of N, these finite-size effects are smoothed out. the correlations across the lattice. The correlations are then reduced to the portion of the lattice that is still in a superfluid phase. Usually that portion is much smaller and in the case ofN b = 120 it is about 40 lattice sites. Experiments produce lattices of size about N = 200 lattice sites. It is enough to study systems of reasonably large size (N ∼ 10 2 ) with V 0 /t following Eq. 2.24. With known scaling, it is possible to simulate any density (N b /N) and any strength of confin- ing potential. The quasiparticle spectrum of the Hamiltonian is significant. It is simply the single particle energy levels of non-interacting fermions on a lattice with a confining potential. Since we are using the harmonic confinement, for n = 1, the groundstate eigenvector looks very similar to the gaussian wave-function of the quantum oscillator in continuum. The reason for this is that the spread in the wave-function, for n = 1 is much larger than the lattice constant, a. Clearly, the lattice effects become significant 23 0 50 100 150 200 n 0 5 10 E(n) Eigenvalue spectrum V(x) ordered E(n) ∝ n 0 50 100 150 200 n -2 -1 0 1 2 3 4 5 E(n) Eigenspectrum V(x), ordered E(n) ∝ n Figure 2.3: Quasiparticle spectra plotted for N = 200 lattice sites. In Fig. (a), the, A = 9.0 and for large n, the spectrum is dominated by the potential energy. In Fig. (b),A = 3.0. The Mott plateau is created at a much larger value ofn and the potential energy never dominates the eigenvalue spectrum. 1 10 100 |i-j| 0.01 0.1 1 10 100 ρ ij N b =40 N b = 120 Figure 2.4: The off-diagonal elements of the one-particle density matrix plotted here for N = 200 andN b = 40,120. A = 9.0 and forN b = 120 the system is in a Mott regime. The straight line fits are forx − 1 2 . As one can see, correlation between bosons follows a power-law decay. 24 0 50 100 150 200 n 0 0.05 0.1 0.15 0.2 ΔE(n) Figure 2.5: The difference between nearest energy levels, also called the derivative spectrum is shown for A = 9.0 and N = 200. Beyond N ≈ 85, the system goes into a Mott regime which is signified by the branching of the energy spectrum. In the continuum limit, this leads to a gapped system. for higher energies, when the variation in the wave-function is significant over a single lattice spacing. The spectra in Fig.( 2.3), for low energies is linear withE(n) ∝ n. This is similar to energy levels of a quantum oscillator. This is in agreement with the behavior of the eigenvector for small n. For large energies, as discuussed earlier, the energy spectrum is dominated by the potential energy. ForA = 9.0, we see the energy spectrum almost overlapping the “ordered” potential energy. Here, we have sorted the potential energy in ascending order, since our eigenvalue spectrum is also in ascending order. The point at which the transition from a linear spectrum to a quadratic spectrum is the point at which the phase transition takes place from a fully superfluid regime to one with a Mott insulator at the center of the trap. Although, the finite systems found in optical lattices have these quasi-condensates and partial transitions, it is crucial to enumerate some of the key aspects of the Mott 25 insulator. Firstly, in the Mott regime, the eigenvectors have a broken symmetry, unlike in the superfluid regime. Low energy eigenstates have even and odd parity for odd and evenn 4 . In the limit of N → ∞, the energy spacing becomes small and the spectrum becomes continuous. However, as we can see, in the Mott regime, the system will remain gapped. Another crucial aspect of this transition is the energy scale at which it happens. There are two useful quantities that drive this transition,V 0 /t andN b . Of course, exper- imentally, fine tuning of the number of bosons,N b is not yet possible. V 0 /t on the other hand is easily adjustable via experiment. For a fixed number of bosons, a shallow trap would lead to a complete superfluid regime while a deep trap would produce a Mott plateau at the center. To quantify the words “shallow” and “deep”, we notice that the onset of the quadratic dependence in the eigenvalue spectrum is the same as the onset of the Mott regime. Hence, a substantially good cross-over point is when the potential energy, V c (x) =V 0 (i c −N/2) 2 > 2t, (2.26) where 2t is the supremum of the tight binding energy E t = −2tcos(ka) or the max- imum kinetic energy when the potential energy is taken to be zero. Numerically, the diagonalization ofH withV(x) identically zero andt = 1 5 renders a highest eigenvalue that tends to the value 2 − . This limit is reached for very large lattice sizes. However, even for N = 5, the maximum eigenvalue is 1.732. For N = 200, we see from Fig. (2.2) that the cross over point in both cases (a) and (b) is close to 2. We have plotted the phase diagram of 1D HCBs in a harmonic confinement. Mott plateaus are formed in the center of the trap with the sides of the trap still superfluid. This can be seen in 4 This is based on the numbering system we have used. We have taken the groundstate to be n = 1 and so the first excited state isn = 2 and so on. 5 Details of the eigenvalue spectrum for this case is discussed in the Appendix 26 0 50 100 150 200 n 0 0.02 0.04 0.06 0.08 0.1 0.12 |ψ n | 2 n=1 n=120 n=121 Figure 2.6: The eigenvectors plotted forn = 1,6,120,121. Forn = 1 we see a gassian- like eigenvector, much like the quantum oscillator eigenfunction. For n = 6 we again see harmonic oscillator-type wave-function. For largen, the symmetry is broken by the presence of a Mott plateau. In this case, the wavefunctions are localized to the right or the left of the Mott plateau. 0 25 50 75 100 125 150 175 200 N m 0 0.0004 0.0008 0.0012 V 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 N b /N Figure 2.7: The phase diagram for 1D HCBs. The white region is a full superfluid phase while the black region is a fully Mott insulating phase. The Mott phase is usually in co-existence with superfluid phases on either side. 27 the color scale of our plot in Fig.(2.4). For a fixed value ofV 0 , one can access all parts of the phase diagram. However, for a fixed value of density,N b /N, one can invariably only access a part of the phase diagram. For example, at half filling,N b /N = 0.5, one will be able to see a Mott plateau for V 0 > 0.0009. However, it won’t be possible to access a high density Mott regime of N m > 150(blue and darker regions). N m is just the total number of lattice sites across which the Mott plateau exists, or alternatively, the total number of HCBs that are participating in the Mott phase. If we are interested in driving the phase transition by changing the density, it is desirable to choose a value of V 0 > 0.0008 or strength, A > 8.0. In this manner, it is easy to access the different phases of Mott densities (different colors in the phase diagram). 2.4 Extension to other confinements Since the potential energy is the diagonal part of the Hamiltonian, it is easy to study other confinements. We have discussed the crossover point for the phase transition, so this is a crucial quantity in determining the parameters that define the confining potential. Experimentally, it has been possible to introduce quartic confinement, where a gaus- sian confinement is superimposed with a repulsive harmonic potential, canceling the first order, quadratic term, of the gaussian confinement. We discuss briefly the general case of V(x) = V 0 x α where α is an even power. It is experimentally perceivable to further manipulate a gaussian confining potential into higher order power confinements. However, since the Taylor’s series of the exponential,e −x 2 /2 has only even powers ofx, we restrict our discussion to such even powers. We parametrize the confining potentials in terms of strengthA, as discussed for the harmonic confinement case in the previous section. A is chosen so that the edge of the confining potential has energy∼ A. Usually, A should be greater than 2t. For our 28 0 50 100 150 200 x 0 1 2 3 4 V(x) α=4 α=6 α=8 α=20 Figure 2.8: Various even confining potentials for N = 200 lattice sites and strength, A = 4.0. In the limitα→∞, we get the infinte square well potential. numerical results, we have takent = 1, without loss of generality. In such a case, typical values ofA range from 2.0 to 10.0. In general, for a powerα, we then have V(x =i) = 4A N α (i−N/2) α , (2.27) where x is discretized by the indexi. We have taken the lattice constant, a = 1. The factor 4 comes from the fact that the center of the trap is ati = N/2. The strength term can be understood if we rewrite Eq. 2.27 in terms of the scaled length V(i) =A i−N/2 N/2 α . (2.28) In the limit asα→∞ one reaches an infinite square well potential. In this limit, one gets a quasiparticle spectrum determined by a tight binding Hamiltonian for non-interacting fermions with zero on-site potential, E TB =−2tcos(ka). (2.29) 29 0 50 100 150 200 n -2 0 2 4 6 E(n) α=4 α=6 α=8 α=20 Figure 2.9: Eigenvalue spectra for different values of α. As α increases, the transition from linear to quadratic takes place at highern. In the limit of largeα one recovers the infinite square well where the energy follows the tight binding energy. Fig. (2.7) shows plots of various eigenspectra for different values ofα. As one can see, forα = 20, one retains the cos(ka) type dependence in the energy spectrum. Also, it is crucial to note that the Mott transition takes place at higher and higher values ofn with increasingα. The reason for that is asα increases, the point at which the potential energy starts becoming more than 2t moves more and more to the edge of the confining potential. In the limit of largeα and even in the case ofα = 20, the maximum energy is 2 and no Mott plateau will be formed. Asα increases, one can also see that the low lying energy eigenvalues change from the linear spectrum for harmonic confinement, to the quadratic spectrum of the tight binding Hamiltonian with zero on-site potential energy. Most confining traps in experiments are now produced by the gaussian beam pro- file of a laser beam. Theoretically, it is convenient to approximate the gaussian by a quadratic, provided the entire length of the lattice is covered by just the “tip” of the 30 0 100 200 300 400 500 n 0 2 4 6 8 V(x) V(x)=9 [1-e -a(x-250) 2 ] V(x)=b(x-250) 2 Figure 2.10: Gaussian confinement plotted with a harmonic confinement, forN = 500 lattice sites with strengthA = 9.0. The parametersa,b are explained in the text. gaussian beam. We briefly describe the case of a gaussian confining potential and the differences in the spectra of that from a harmonic confinement. We parametrize the gaussian potential, such that the small x expansion yields a quadratic confinement that is reasonable. Since we require the gaussian potential to be attractive, we choose V g (x) =A 1−e −a(x−250) 2 . (2.30) For small x, we have V g (x)≈A−A 1−a(x−250) 2 + a 2 (x−250) 4 2 +... , ≈A·a(x−250) 2 +O((x−250) 4 ). (2.31) We have chosen the gaussian such that when the gaussian goes to zero, the potential energy isA = 9.0. Also, the only other parameter in the gaussian potential, the standard 31 0 100 200 300 400 500 n -2 0 2 4 6 8 10 E(n) Gaussian confinement Harmonic confinement 0 100 200 300 400 500 n 0 0.05 0.1 ΔE(n) Gaussian confinement Harmonic confinement Figure 2.11: Comparison of eigenvalue spectrum of the gaussian confinement with a harmonic confinement of similar strengths of A = 9.0 with N = 500 lattice sites. In Fig. (a) we see that the spectra are same for low energies. The main difference comes at higher energies. The bounded nature of the gaussian, reduces the overall “height” of the energy spectrum as opposed to that of the harmonic confinement. deviation or spread is defined in such a way that the first order in the Taylor’s series expansion matches that of the quadratic potential. Then, a = 1/250 2 so that V g (x) approximatesV(x) =A·( x−250 250 ) 2 . As one can see in Fig. (2.8,2.9), the harmonic confinement is a reasonable approx- imation of the gaussian confinement. The actual point at which the phase transition takes place is shifted only slightly to a higher value of n. The difference is only about Δn = 20, which is less than 5% of the total lattice size. Also, as one can see, the shape of the eigenvalue spectrum does not change much. The main difference can be made out from Fig. (2.8) where we clearly see that the edge of the harmonic potential is about 2.5 units higher than that of the gaussian potential. This is difference is mirrored in the difference in the energy spectra. 32 Chapter 3 Time evolution of one dimensional hard-core bosons The study of non-equilibrium dynamics has been shown to be useful in quantum gases and optical lattice experiments(10). The anisotropic expansion of BECs(1; 6; 9) helped understand interparticle interaction in the BEC. Ever since, it has become increasingly more popular to study the dynamics of bosons and fermions in reduced dimensions. Theoretically, there has been a lot of work done in the area of dynamics when the system has been in an initial excited state(44). Also, sudden changes in the Hamiltonian have been simulated with interesting results(48; 45). In this chapter we present a technique first developed by M. Rigol and A. Muramatsu(47) to study the exact dynamics of 1D HCBs. Since it is possible to diag- onalize the full Hamiltonian of 1D HCBs, it is easy to extend the study to the time evolution of time independent Hamiltonians and Hamiltonians where a sudden change is introduced. We review the technique and using this technique we measure the two characteristic frequencies of 1D HCBs, namely the breathing mode frequency (ω B ) and the sloshing mode frequency (ω S ). We compare this with the theoretical value of the ratio of breathing mode and sloshing mode(38). We then go on to develop a technique to study the dynamics of 1D HCBs when there is an explicitly time dependent Hamiltonian. Experimentally, it is easy to introduce time dependence in the Hamiltonian. We use our technique to illustrate the dynamics of 1D 33 HCBs when the center of the trap is oscillated sinusoidally. This leads to enhanced oscillations. 3.1 Equal-time Green’s functions and time evolution To study the time evolution we follow the work done by M. Rigol and A. Muramatsu (47) we continue from the previous section to approach the time evolution of HCBs. The Green’s function introduced in chapter 3, equation (2.10) is a good starting point. G ij =hΨ 0 |b i b † j |Ψ 0 i. (3.1) To calculate the time evolution, we use the Schrodinger picture. The ground state wave function at any time,τ 1 |Ψ 0 (τ)i =e −iτH/~ |Ψ 0 (0)i. (3.2) Then, the equal-time Green’s function at timeτ is G ij (τ) =hΨ 0 (τ)|b i b † j |Ψ 0 (τ)i. (3.3) If we proceed to apply the Jordan-Wigner transformation, we can evaluate the non- interacting fermionic ground-state. The time-evolved state is then, |Ψ f 0 (τ)i =e −iτH f /~ |Ψ f 0 i = N f Y m=1 N X n=1 P nm (τ)f † n |0i, (3.4) 1 Throughout this exposition, we use the notation τ for time to avoid confusion with the hopping parameter in the Hamiltonian, which ist. 34 where the matrix P nm (τ) is just e −iτH f /~ P. This matrix-product can be simplified by diagonalizing the exponential matrix. If U † H f U = D, where D is some diagonal matrix, we can write,H f =UDU † and we can simplify the exponent e −iτH f /~ = ∞ X j=0 (−iτ) j H j f j!~ j = ∞ X j=0 (−iτ) j j!~ j UDU † UDU † ...UDU † =Ue −iτD/~ U † . (3.5) Then, we can write the action of the exponent on the P matrix, P(τ) =e −iτH f /~ P =Ue −iτD/~ U † P, (3.6) Following the steps in the previous section we arrive at G ij (τ) = det (P 1 (τ)) † P 2 (τ) . (3.7) This approach leads to exact dynamics of hard-core bosons in a trap. Extensive studies of exact dynamics has been done(47; 48; 45). We use this technique to study the ratio of the breathing mode and sloshing mode frequencies of a 1D HCB. 3.2 Exact dynamics approach to study fundamental modes of a gas of one-dimensional hard-core bosons We use the above technique to demonstrate experimental measurements of the breathing mode and sloshing mode of 1D quantum gases. The breathing mode is just the lowest compressional mode of the lattice. We study the compressional mode by making a 35 slight change to the potential. Experimentally, the compression is brought about by a compression of the lattice size or reducing the lattice constanta. In our system, we write our potential energy as: V(x) =V 0 a 2 N X i=1 (i−i 0 ) 2 , (3.8) wherei 0 is the center of the trap. For most of the results in the text, we have taken the lattice constant a = 1. Now, a compression of the lattice leads to more lattice points coming closer to the center of the trap keeping the trap curvature fixed. However, if we change the lattice separation, the the rate at which the curvature of the trap changes, will also change. Hence, a change in lattice separation leads to a change in curvature. If we compress the lattice, it leads to a reduction in the curvature of the trap. Then, a compression atτ = 0, froma→ka, where 0<k < 1 leads to H ii (τ < 0) =V 0 a 2 N X i=1 (i−i 0 ) 2 , H ii (τ ≥ 0) =V 0 a 2 k 2 N X i=1 (i−i 0 ) 2 . (3.9) On the other hand, the sloshing mode, also known as the dipole oscillation mode is brought about by a slight horizontal shift of the center of the trap at τ = 0. Ifi 0 is the center, a horizontal shift results ini 0 → i 0 +δi i . δi is typically a few lattice sites. Then, the Hamiltonian can be written as H ii (τ < 0) =V 0 a 2 N X i=1 (i−i 0 ) 2 , H ii (τ ≥ 0) =V 0 a 2 N X i=1 (i−i 0 −δ i ) 2 . (3.10) 36 0 100 200 300 400 500 τ 0.8 1 1.2 1.4 n(k=0,τ) (a) 0 100 200 300 400 500 τ 0 2 4 6 8 10 x cm (b) Figure 3.1: Figure (a) is the oscillations of the peak of the momentum distribution func- tion,n(k) max . Figure (b) is the center of mass oscillations,x cm . In both casesN b = 20 HCBs and N = 100 lattice sites were used. In figure (a), the trap is suddenly relaxed withV 0 → 0.9V 0 , exhibiting a “breathing mode” of vibration. In figure (b), the trap is suddenly shifted to the right byδi = 5, exhibiting a dipole oscillation (sloshing mode). 0 ω 0 2ω 0 3ω 0 4ω 0 5ω 0 6ω 0 7ω 0 ω 0 2 4 6 8 10 Frequency spectra Sloshing mode, F.T.[x cm (τ)] Breathing mode, F.T.[n(k=0,τ)] Figure 3.2: Frequency spectra plots of the center of mass oscillations of a sudden shift in trap and the zero ofn(k) when the trap strength is suddenly changed. The frequency spectra show dominant contributions atω S = 2 √ V 0 andω B = 4 √ V 0 . 37 For a 1D degenerate gas the mean-field value of (ω B /ω S ) 2 = 3. For a 3D elon- gated condensate, (ω B /ω S ) 2 = 5/2(55; 37; 7). For 1D HCBs it is anticipated to be (ω B /ω S ) 2 = 4. The sloshing mode vibration also known as dipole oscillation is seen in the real space density. The center of mass of the superfluid begins to oscillate, since the sudden change effectively moves the entire density to an excited state of the harmonic oscillator. The oscillations seen are collective oscillations of the entire gas. It should be remarked here that the 1D HCB model is an integrable model and even though the oscillations appear to die out, they do so only in a collective sense. One could, for instance attribute this to a simple interference of wavefunctions that leads to the reduction of oscillations. It has been shown that integrable systems do not “thermalize” in this sense(45; 44). The oscillations are shown to recur at long times(45), consistent with the fact that there is no dissipation in the system. The breathing mode is more difficult to observe since the center of mass does not change. We present a technique by which it would be easy to experimentally measure these oscillations. The momentum distribution function (MDF) is an easily measurable experimental quantity and we study the peak,{n(k)} max of the MDF. In principle, the peak position could change as a function of time if higher momentum modes are excited. However, we have seen for the case of the breathing mode measurement, the peak of the MDF remains atn(k = 0). In figures (3.1,3.2), we see the computational measurement of the respective modes. In figure (3.1), we see time series of (a)n(k = 0,τ) and (b)x cm (τ) forN b = 20 bosons on aN = 100 site lattice. In (a), we have changed the strength of the potentialV 0 → 0.9V 0 which amounts to a change in lattice constant from a → 0.95a, a 5% change. In (b), we have shifted the trap by δ i = 5 lattice sites, again a shift of 5%. In figure (3.2) the frequency spectra of the time series from figure (3.1) are plotted. The frequency plots 38 are scaled appropriately so that a comparison can be made. In this particular discussion, we are only interested in seeing which frequency modes are excited. The actual scale of the plot is not significant. As we can see, in figure (3.2), the breathing mode is dominated by frequency 2ω 0 , where ω 0 = 2 √ V 0 is the natural frequency of the trap. From figure (3.1)(a) it is evident that the higher mode oscillation of4ω 0 in the breathing mode is only a transient feature that dies out leaving a pure2ω 0 mode in the steady-state. To a good approximation, then we have ω 2 B ω 2 S ≈ 4. (3.11) This is in agreement with other predictions(55; 37; 7). ω 2 B ω 2 S = 4 is also the frequency ratio for a thermal gas and a gas of degenerate, non-interacting fermions(55). Apart from recent experiments, there have been several theoretical and computational studies(35). 3.3 Time-dependent traps: Splitting the Hamiltonian We proceed to discuss the case when the Hamiltonian is time-dependent. The exact approach is no longer valid because the time-dependent Hamiltonian cannot be expo- nentiated since the Hamiltonians at different times do not commute with each other. We present an approximation that helps us solve the problem. For time-dependent traps, we notice that the Hamiltonian has two parts: 1. The time-independent term which can be diagonalized atτ = 0,H 1 2. The time-dependent term is small, but diagonal,H 2 (τ) Then, the Hamiltonian for all times is written as H(τ) =H 1 +H 2 (τ)θ(τ), (3.12) 39 where theθ(τ) ensures thatH 2 switches on at timeτ = 0. We further chooseH 2 in such a way that for every value ofτ, ||H 2 || ||H 1 || < 1, (3.13) where||A|| is the 1-norm of A, i.e. the maximum column sum of A. This condition limits the technique to small driving forces. This is not a sufficient condition because the norm only looks at the maximum column sums. ForH 2 this is simply the maximum element of the matrix, since it is diagonal. However, there are perceivable situations, like at the center of the trap, where at a certain time τ, H 2 is maximum and H 1 is minimum. Generally, such points at which this happens are only a few compared to the total number of lattice points. What we mean by equation (3.13) is that the scale at whichH 2 varies is smaller than the scale ofH 1 . We split the time interval into M very small intervals of length Δτ = τ/M so that the time dependent part of the Hamiltonian, namelyH 2 , is nearly constant with respect toH 1 within the interval. In order to get an upper bound on the value ofΔτ, we further require that ||{H 2 ((k+1)Δτ)−H 2 (kΔτ)}||/||H 2 || = Δτ ||H 2 || ∂H 2 ∂τ +O Δτ 2 ||H 2 || ≪ 1. (3.14) IfH 2 (τ) has a sinusoidal dependence, such that H 2 = ΔV sin(q 0 τ), (3.15) where ΔV is the strength of the driving force. Then, equation (3.14) reduces to the condition q 0 Δτ ≪ 1. (3.16) 40 Clearly, it is then necessary to have ΔV,Δτ appropriately chosen for this technique to work. In order to evaluate the unitary propogator from thek−1 th time interval to thek th time interval, we use the Suzuki-Trotter decomposition(58; 56). To understand how this works, we start with a general formula for e A e B called the Baker-Campbell-Hausdorff formula log(e A e B ) =A+B + 1 2 [A,B]+ 1 12 [A,[A,B]]+.... (3.17) Then, we have e A+B =e A e B e − 1 2 [A,B] . (3.18) Of course, this result would be exact if the commutator[A,B] = 0. This is not the case for most driving forces. We assume that the driving force commutes with the potential, so that the norm of the commutator ||[H 1 ,H 2 ]||∝|tΔV|, (3.19) where t is the hopping parameter. Since|ΔV/t| < 1, the higher order commutation relations can be neglected 2 . Then, followiing Suzuki(56; 51), we have e − 1 2 [H 1 ,H 2 ]Δτ ∼ = 1− 1 2 O(Δτ 2 ). (3.20) We can then write the unitary propogator at thek th interval as e −iΔτ(H 1 +H 2 (kΔτ))/~ =e −iΔτH 1 /~ e −iΔτH 2 (kΔτ)/~ +O(kΔτ 2 ) =Ue −iΔτD 1 /~ U † e −iΔτH 2 (kΔτ)/~ +O(kΔτ 2 ), (3.21) 2 Some reviews assume that [H 1 ,H 2 ] commutes with bothH 1 andH 2 (51) 41 where U is a matrix that diagonalizes H 1 and D 1 is the diagonal matrix. The error at every step is of the order of Δτ 2 and afterk = τ/Δτ steps this error is proportional to Δτ. The decomposition is called First order Trotter decomposition. By splitting the Hamiltonian in this manner, we have to compute the matrixU only once at timeτ = 0. Diagonalizing the Hamiltonian in every step, would lead to a very costly code. However, in this case the computation time goes down tremendously. This result is now introduced into equation (3.6) to derive the time-evolved P-matrix at the k th time-step e −iΔτ(H 1 +H 2 (kΔτ)/~ P =Ue −iΔτD 1 /~ U † e −iΔτH 2 (kΔτ) P. (3.22) The rest of the computation follows the time evolution of the time-independent system discussed in the previous two sections. Since our computation has approximations, we need to check the errors due to these approximations. The discretization error is controlled by a suitable choice of Δτ. It may not be always clear what the choice of Δτ should be and in such a case, it would be necessary to check the condition in equation (3.14) either analytically or numerically for all values ofτ. The second error is due to the First order Trotter decomposition of the time evolution operator. Clearly the technique cannot be used for arbitrarily long times. The Trotter error is limited by the total time,τ and the discretization of time, Δτ. One can further reduce the Trotter error by doing a second-order Trotter decom- position. Like in the case of t-DMRG(51), any higher order Trotter decompositions comes with the possibility of rounding errors due to the decomposition of one matrix into several matrices at each time step. We present a second order Trotter decomposition without delving into higher orders. 42 0 50 100 150 200 τ 0 0.2% 0.4% 0.6% 0.8% 1% 1.2% 〈n e 〉−〈n a 〉 Figure 3.3: The propagation of error for N b = 20 bosons on a N = 100 site lattice using a first order Trotter decomposition. The error is calculated as the absolute differ- ence between expectation values of the number of bosons for an exact evolution and an approximate evolution when the trap curvature is suddenly changed fromV 0 → 0.9V 0 . It is clear from the grap that the error is linear in total timeτ. For the second order Trotter decomposition, the time evolution operator can be writ- ten in the symmetric form e −iΔτ(H 1 +H 2 (kΔτ))/~ =e −i Δτ 2~ H 2 (kΔτ) e −iΔτH 1 /~ e −i Δτ 2~ H 2 (kΔτ) +O(kΔτ 3 ) =e −i Δτ 2~ H 2 (kΔτ) Ue −iΔτD 1 /~ U † e −i Δτ 2~ H 2 (kΔτ) +O(kΔτ 3 ), (3.23) such that the overall error is of the orderO(Δτ 2 ). In order to calculate the error, we take a system of 20 bosons on a 1D lattice of 100 lattice sites. The system appears to be small, but experimentally, a 1D lattice far greater 43 than 100 sites is not achievable. The starting configuration is a superfluid phase. At τ = 0, we suddenly reduce the trap curvature V 0 →V 0 −ΔV 0 , (3.24) here ΔV 0 = 0.1V 0 . We then run approximated code, using a first order Trotter decom- position, to calculate the final configuration after some time τ. Since this is a time- independent evolution, we know the final configuration and can calculate the error in the real space density. The error is shown in figure (3.3). In figure (3.3) we see that the error increases linearly with time,τ up to a maximum of 1% error for τ = 200. Here τ is measuring in units of t/~ and in our computation, we have takent = ~ = 1. 3.4 The Driven oscillator We study the dynamics of the trap when it is driven. Classically, if an oscillator is driven with a time-dependent force that has no spatial dependence, we get V driving (x) =−F(t)·x, (3.25) whereF(t) is the initial force function. In general, one can add this term to a quadratic potential centered at 0, resulting in an effective potential V eff (x,τ) =V driving (τ)+ 1 2 mω 2 x 2 = 1 2 mω 2 x− F(t) mω 2 2 − F 2 (t) 2mω 2 . (3.26) 44 There are two time-dependent terms in this. Firstly, the minimum of the potential will have a time dependence and then, there will be a time-dependent shift in energy at every time step. Experimentally, it would be difficult to continuously shift the energy scale as a function of time. However, it is has been shown(4) that optical lattices can be moved across a few lattice sites with precision. Hence, it seems only plausible that this motion can be ramped to a time varying function of some kind. We discuss below the case when the center of the trap is oscillated with a sinusoidal function, ignoring the second term in the effective potential. If the trap is displaced by some distanceδ i , it has been shown in the previous sections that the center of mass of the bosons follows a damped oscillator motion (48). We now proceed to compare this motion with that of a driven oscillator with all the other parameters the same. In particular, the displacement used in the previous section was δ i = 5 lattice sites. We then drive the trap so that it oscillates betweeni 0 ±5, wherei 0 is the center of the trap. The time dependent confining potential can be written as V(x) =θ(−τ)V 0 N X i=1 n i (i−i 0 ) 2 +θ(τ)V 0 N X i=1 n i (i−i 0 −δ i sin(q 0 τ)) 2 , (3.27) withδ i = 5 andq 0 = ω 0 , the natural frequency of the trap. We have chosen the driving frequency to be the natural frequency because the natural sloshing mode vibration of 1DHCBs was shown to be dominated by the frequency,ω s =ω 0 . The second term contributes only for τ > 0 and the initial configuration is deter- mined by the first term atτ < 0. Forτ > 0 we need to be able to write the Hamiltonian in terms ofH 1 , which is time independent andH 2 which is time dependent, such that ||H 2 || ||H 1 || < 1. (3.28) 45 0 100 200 300 400 500 τ -15 -10 -5 0 5 10 15 x cm 0 100 200 300 400 500 0 Driven Quenched Figure 3.4: Enhanced oscillations of the center of mass motion of N b = 20 bosons on N = 100 lattice sites. The blue curve represents collective oscillations when the trap is quenched or shifted to the right. The oscillations certainly do not die down in the 10 time periods in which they were measured. Here,V 0 = 9.0/(N/2) 2 = 0.0036 In order to do this, we extract out a time independent term fromV(x,τ > 0) V(x,τ > 0) =V 0 N X i=1 n i (i−i 0 −δ i sin(q 0 τ)) 2 =V 0 N X i=1 n i (i−i 0 ) 2 +V 0 N X i=1 n i sin(q 0 τ) δ 2 i sin(q 0 τ)−2δ i (i−i 0 ) . (3.29) The first term, clearly matches the potential energy term atτ < 0 while the second term has all the time dependence and is also diagonal. Both terms seem to have the same strength, V 0 , but following the discussions in the previous chapter, one can carefully notice that since V 0 = 4A N 2 , (3.30) 46 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 ω 0 0.5 1 1.5 2 2.5 3 x cm (ω) Figure 3.5: The discrete Fourier transform of the center of mass motion of N b = 20 bosons onN = 100 lattice sites. The frequency spectrum is dominated by two nearby frequencies close to the driving frequency,q 0 =ω 0 . the first sum is dominated by the (i−i 0 ) 2 term. In the second sum,δ 2 i term has a 1/N 2 contribution while(i−i 0 ) term has a1/N contribution. To quantify this, for a lattice of N = 100, withA = 5.0, we have V 0 = 1 500 = 0.002. (3.31) In the first summation, at the edges of the trap (where V is maximum), we have contributions ∼ 4A = 20.0, while the second sum has maximum contributions ∼ δ 2 i 4A/N 2 = 0.05 and ∼ (−4A/N)δ i = −1.0. Then, the condition (3.28) is sat- isfied for all values of τ. We have to also be careful in choosing Δτ since there is a sin 2 (q 0 τ) term in the time dependence. Then we have to take Δτ ≪ 1 2q 0 . We seek out some resonant behavior of the bosons in the trap. In order to numeri- cally compute a resonance curve for this driven oscillator, we would have to drive the 47 oscillator for very long times, so that we can calculate the Fourier transform. How- ever, we run it for 10 time periods of the driving force and calculate a discrete Fourier transform of the time series. In figure (3.4), we see enhanced oscillations of the center of mass when the trap is driven at a frequency q 0 = ω 0 . The blue curve represents a sudden shift which results in damped oscillations as discussed in previous sections. As one can see in figure (3.5), two nearby frequencies contribute in the center of mass motion. The two frequencies are very close to the driving frequency. The enhanced oscillations show that it is possible to experimentally set up driven dipole oscillations. It isn’t, however, clear whether these oscillations will continue for long times (steady state) or whether these are transient features. In the following chapter we give a detailed discussion of parametrically driven traps and show that even in the short time scales of our numerical experiment, steady state properties of the quantum gas can be recovered. 48 Chapter 4 Parametric excitations of hard-core bosons in a time-varying optical lattice 4.1 Introduction Computational studies of time dependent confining potentials have only been done recently(30) for soft-core bosons using t-DMRG. Although experimentally, it is a fea- sible extension of known techniques, there have been few recent studies of excitation spectra of bosons(2). 1D HCBs offer an ideal testing ground to study time-dependent potentials for quantum particles and also see the effects of driving forces on both sides of the superfluid-Mott insulator transition. In this chapter we examine the dynamics of 1D lattice HCBs in the presence of a time dependent confining potential where the curvature of the potential is sinusoidally varied betweenV 0 ±ΔV 0 . We follow the numerically exact approach developed in the previous chapter which was generalized from previously used formalism for suddenly shifted potentials(48). As will be discussed below, the parametric excitations brought about by the time variation of the trap lead to a probe into the internal structure of the spectrum in both superfluid and Mott insulating phases. The experiment to study the dynamics of HCBs when the trap curvature is modulated in time is sketched Fig. (4.1). This type of sinusoidal time variation can easily be introduced into an already existing experimental setup for 1D trapped bosons. Which dynamic quantity would be easily accessible to experiments and show markedly different response in the superfluid and 49 10 20 30 40 x 0 0.1 0.2 0.3 0.4 0.5 n(x) Density at τ =10t/h - Density at τ = 0 q 0 V 0 ±δV 0 Figure 4.1: Density plots of15 hard-core bosons trapped in aN = 50-site lattice, plotted at timesτ = 0 and atτ = 10t/~. Atτ = 0 the HCBs are in a superfluid phase with the spatial confinement given by a harmonic oscillator potential withV 0 a 2 = 4.0× 10 −3 t. For τ > 0, the curvature of the potential is sinusoidally modulated with a frequency q 0 = 1.26×10 −1 h/t and amplitudeΔV 0 = 0.1V 0 Mott insulator regimes? In the following, we show that the frequency dependence of the zero-momentum peak of the momentum distribution function (MDF), n(k = 0,ω) displays such characteristic signatures for the various competing phases. 4.2 Formalism to study time varying traps The Hamiltonian for a system of a system ofN b 1D HCBs in a trap potential is given by H =−t N X i=1 ˆ b † i ˆ b i+1 +h.c. +V 0 (τ)a 2 N X i=1 ˆ n i (i−i 0 ) 2 , (4.1) where ˆ b i and ˆ b † i are boson annihilation and creation operators at lattice sites i, a is the lattice length,t is the nearest-neighbor hopping parameter, and the time dependent trap- ping potential is given by V 0 (τ) = V 0 + ΔV 0 sin(q 0 τ). Here, ΔV 0 ≪ V 0 and q 0 is the 50 -4 -2 0 2 4 ka 0 0.5 1 1.5 2 n(k) -4 -2 0 2 4 0 0.5 1 1.5 2 τ=1 τ=25 τ=51 0 π/16 π/8 ka 0 0.5 1 1.5 n(k) Figure 4.2: Momentum distribution functions for at various times for a system ofN b = 15 bosons on a lattice with N = 50 sites. The momentum distribution function is characteristically peaked atk = 0. The three plots are shifted to resolve the peak heights. For small momenta, one can see the variation in the inset figure driving frequency (Fig. 4.1). Unless mentioned otherwise, without loss of generality in the following we choose q 0 = ω 0 = 0.126h/t, where ω 0 is the characteristic fre- quency of the trap. Classically, parametric resonances are expected at integer multiples of the characteristic frequency q 0 = ω 0 (33). Due to the numerical limitation that we need to avoid Hamiltonians that vary too rapidly in time, here we focus on the case q 0 = ω 0 . Also note, that in Eq. (4.1) it is implicitly assumed that the hard-core con- straint has been applied, i.e. the operators anti-commute on-site,{ ˆ b i , ˆ b † i } = 1. Then, using the Jordan-Wigner transformation(26), one can transform the Hamiltonian into effective non-interacting fermions. To calculate dynamical observables, one needs to evaluate the equal-time Green’s function, G ij (τ). We follow the technique developed in Refs. (46; 45) to study the dynamics of 1D HCBs, which involves solving for the Green’s function in the basis of the effective non-interacting fermions. However, since in the case we are considering the Hamiltonian is explicitly time dependent, the technique developed in Refs. (46; 45) 51 needs to be generalized. First, where, U 0 is the unitary matrix which diagonalizes H 0 at timeτ = 0, and D 0 is the diagonal matrix containing the eigenvalues of H 0 . In this manner, we need to perform the diagonalization only once atτ = 0. We follow the technique outlined in the previous chapter to study time-dependent traps. Very briefly, we note that the Hamiltonian can be split into a diagonal time- dependent part, H D 1 (τ) = ΔV 0 sin(q 0 τ) P i ˆ n i (i− i 0 ) 2 and a time-independent part, H 0 , that still needs to be diagonalized. We then use a first-order Suzuki-Trotter expansion(58; 56) to write the time evolution operator atτ =kΔτ as e − i ~ (H 0 +H D 1 (τ))kΔτ =U 0 e − i ~ D 0 kΔτ U † 0 e − i ~ H D 1 (τ)kΔτ +kO(Δτ 2 ). (4.2) To perform a time-series of up toτ = 500 t ~ , we have usedM = 10,000 time-steps with Δτ = 0.05 t ~ . 4.3 Results forN = 50 lattice sites Turning now to the results of the numerical analysis, we first note that even for a rela- tively small lattice withN = 50 sites, one observes a well-defined transition from the superfluid to the Mott insulator regime. For the particular parameter set chosen here, with V 0 a 2 = 4.0× 10 −3 t, we find that a Mott plateau begins to form at the critical density ofN b = 37 bosons. Following Ref. (45), we evaluate the momentum distribution function (MDF) in terms of the one-particle density matrix at timesτ n(k,τ) = 1 N N X l,m=1 e ik(l−m) ρ lm (τ). (4.3) 52 0 100 200 300 400 500 τ 0.1 0.12 0.14 0.16 n(k=0,τ)/N b N b =1 N b =10 N b =10 0 100 200 300 400 500 τ 0.015 0.02 0.025 n(k=0,τ)/N b N b =42 N b =44 N b =46 Figure 4.3: Several plots of n(k = 0,τ) are shown here. The upper plot is for the superfluid case. The time series are divided byN b and are vertically adjusted to be well separated from each other. The scale at which oscillations take place in the Mott regime are drastically smaller. Time-of-flight measurements in experiments give a direct representation of MDFs pro- viding a link between theoretical and experimental studies. For 1D bosons, the super- fluid phase is well represented by a sharply peaked MDF at k = 0 signifying a macro- scopic occupation of the zero momentum state. The Mott insulating regime is identified by a broad peak with a spread over all momenta. The height of the peak is much lower and is present only because of the finite size of these experiments. In Fig. 4.2 we see sharply peaked MDFs for three different timesτ. Note that in this plot and the data sets, we have taken ~ = t = 1. The number of bosons here isN b = 15 which is well within the superfluid regime. In the inset of Fig. 4.2 we have plottedn(k) for small values of momenta at the same times τ. Although the MDF is sharply peaked at k = 0, we see significant variation of the zero momentum peak height as compared to the rest of the MDF. 53 Fig. 4.3 is a collection of plots of n(k = 0,τ) for various densities of bosons. The time series are normalized byN b , the number of bosons. We have added slight vertical shifts to the plots so that the various curves are not overlapping. As one can see, in the superfluid regime the magnitude of the time variation of n(k = 0) is about 5-10 times bigger. The system is in a phase of high coherence and this is reflected in an overall higher averagen(k = 0). Also, the time variation appears to be non-trivial and perhaps contains contributions from more than a few frequencies. The time variation is performed over10 time periods of the driving force. To convert this time series to a frequency series, we perform a numerical discrete Fourier transform on the data. The Fourier basis is numbered according to the discrete frequencies,ω n = nπ/τ T , whereτ T is the total time of 10 time periods. The total scale of this time variation is of the order of ΔV 0 /V 0 or 10% of theτ = 0 peak height. The time dependence is not purely sinusoidal, however a dominance of a single frequency is clear. As we will show ahead, the frequencies that dominate the oscillations of n(k = 0,τ) are close to the driving frequencyω 0 along with non-trivial frequency contributions of a smaller scale around overtones of the driving frequency. In Fig. 4.4 the frequency dependence of the MDF is shown for the superfluid (a) and for the Mott insulator (b) regimes. n(k = 0,ω), the discrete Fourier transform ofn(k = 0,τ) is a finite-time, real-valued Fourier transform of the input data,n(k = 0,τ). It has large positive and negative values since n(k = 0,τ) has cosine and sine contributions. The Fourier transform is defined for even values of discretizedω j as n(k = 0,2ω j −2) = M X i=1 n(k = 0,iΔτ)cos((w j −1)(i−1)2π/M), (4.4) 54 0 ω 0 2ω 0 3ω 0 ω 0 0.2 0.4 0.6 0.8 1 |n(k=0, ω)| /N b N b =1 N b =5 N b =10 N b =15 ω 0 2ω 0 3ω 0 ω 0 0.05 0.1 0.15 0.2 0.25 0.3 N b =40 N b =42 N b =45 ω 15 (a) (b) Figure 4.4: The normalized frequency dependence of the zero-momentum peak of the momentum distribution function is shown for various total numbers of bosonsN b in the trap. In (a) the system is initially in a superfluid phase, and the normalized dependence is an order of magnitude larger than in (b), where the system is initially in a Mott insulator regime. In (a), one observes a sharp peak atω 0 and satellite peaks close to 2ω 0 and 3ω 0 . These satellite peaks are shifted from 2ω 0 and 3ω 0 with increasing N b . ω 15 represents one such shift from 2ω 0 forN b = 15. In (b), the system is in a Mott insulator, and one observes a sharp resonance atω 0 with very small contributions elsewhere. while for odd values ofω j as n(k = 0,2ω j −1) =− M X i=1 n(k = 0,iΔτ)sin((ω j −1)(i−1)2π/M). (4.5) We have plotted the magnitude,|n(k = 0,ω)|. As can be clearly seen, the frequency dependence of the MDFs are at least an order of magnitude larger for superfluids in Fig. 4.4(a), compared to Mott insulators, Fig. 4.4(b). Further, one can clearly observe in Fig. 4.4(a) multiple frequencies dominating the spectrum of n(k = 0,ω). We see a sharp peak atω =ω 0 signifying parametric resonance at the driving frequency. This is a collective mode in which most of the superfluid participates. We also see satellite peaks close toω = 2ω 0 . In the superfluid regime, these peaks shift to lower frequencies asN b 55 10 20 30 40 N b 0.1 0.2 0.3 0.4 0.5 ω (in units of t/h) 10 -5 10 -4 10 -3 10 -2 |n(k=0,ω)|/N b Figure 4.5: Contour map of the normalized frequency dependence of the zero- momentum peak of the momentum distribution function n(k = 0,ω) for all values of N b , the number of bosons. The contour color is in logscale and is quite noisly with several harmonics have small contributions. In the Mott regime, when N b > 37, it is difficult to extract useful information from this contour plot. is increased. In the limiting case of very low density, ρ b → 0 the effects of the lattice are negligible. In this case, the single particle makes transitions from the lowest energy level to the second and third energy levels. In the effective fermion picture, the satellite peaks are transitions between single- particle energy levels. In the absence of a lattice, we know that the single-particle energy spacing of a quantum oscillator is constant and equal to ~ω 0 . However, in the presence of a lattice(46), the single particle energy spacing reduces as N b , i.e. the number of bosons in the trap, is increased. This reduction in energy spacing is mirrored in a red- shift of the satellite peaks. One also observes smaller peaks, e.g. for N b = 15 close to 3ω 0 . These are single particle transitions between multiple energy levels. In order to better understand these shifitng transitions, in Fig. 4.5 we consider a contour plot of|n(k = 0,ω)| as a function of N b and ω. Since the individual plots of |n(k = 0,ω)| are very noisy in the logscale, we can refine the contour plot by keeping the 20 highest peaks for each plot of|n(k = 0,ω)| and making the rest of the function 56 |n(k = 0,ω)| equal to 10 −12 . The number 20 was adjusted to give the clearest continu- ous contours across the scale ofN b . Larger values would result in noisier results, while smaller values will reduce the detail. In Fig. 4.6(a), we have a contour plot of|n(k = 0,ω)| with the noise removed. For each value of N b , we have kept the 20 highest peaks that contribute to the frequency dependence of the zero-momentum peak. By doing this, we can extract useful infor- mation of frequency dependence in the Mott regime. In Fig. 4.5 it was not possible to see this information because the order of magnitude of the peaks in the Mott regime is similar to the order of magnitude of the noise in the superfluid regime. Note that N b , the number of bosons, is not a dynamical quantity. Instead, we are working in the microcanonical ensemble with fixed total energyE and number of par- ticles, N b . Each plot for different N b corresponds to a different initial configuration of the trapped system. The contour plot is a collection of 50 such initial configurations for values of N b ranging from 1 to 50. In Fig. 4.6(a), the magnitude of n(k = 0,ω) is encoded into the color, using a log scale. As one can see, with increasing boson fill- ing, forN b ≤ 37, satellite peaks develop into definite contour lines, roughly parallel to the N b -axis. In Fig. 4.6(b), the effective single-particle transition energies are plotted as a function of N b . The lowest curve, Δε 1 = (ε i+1 − ε i ) is the difference between nearest energy eigenvalues of the time-independent Hamiltonian, H 0 after the Jordan- Wigner transformation. These correspond to the single-level transitions for a system of the effective spinless non-interacting fermions in a 1D lattice with a harmonic trap. Similarly,Δε 2 = (ε i+1 −ε i−1 ), etc. As it is clear from comparing Fig. 4.6(a) and Fig. 4.6(b), the difference of the effec- tive energy lavels precisely matches the contours of the magnitude ofn(k = 0,ω). Note that theω-axis in Fig. 4.6(a) is in units of t h which is the same scale in which our Hamil- tonian is numerically diagonalized. The close agreement between Fig. 4.6(a) and Fig. 57 (a) 10 20 30 40 N b 0.1 0.2 0.3 0.4 0.5 ω (in units of t/h) 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 |n(k=0,ω)|/N b (b) 0.1 0.2 0.3 0.4 0.5 10 20 30 40 Δε i N b Δε 1 Δε 2 Δε 3 Δε 4 Figure 4.6: (a) Contour map of the amplitude of the normalized momentum distribution functionn(k = 0,ω) for all values ofN b , the number of bosons. (b) Differences between the effective single particle energy levels of the diagonalized Hamiltionian atτ = 0. The different curves in (b) represent different inter-level transitions,Δε i =ε(N b +i)−ε(N b ). The contours in (a) of the excitation spectrum of n(k = 0,ω) match the transition energies in (b). ForN b ≥ 37 a Mott plateau is present in the system. (a) shows that the system continues to follow multi-level transitions in the Mott regime, however with an order of magnitude decrease in amplitude. 4.6(b) indicates that the effective single-particle excitations are due to “edge” excitations of the superfluid in the trap. Every particle added to the trap is added to the next energy level since they are similar to fermions. The edge particles have the highest energy and excitations of these particles will manifest in differences of the highest occupied energy state with the next state. These are exactly the differences shown in Fig. 4.6(b). In Fig. 4.7 the variation of the maximum peak height of the MDF,|n(k = 0,ω)| is shown. For small to moderate fillings one observes that the maximum peak height scales linearly with N b . The curve that fits the data (red line) has a slope of m = 1.013. In the Mott regime, we see that the maximum peak height of the MDF is sharply reduced, demonstrating that in the Mott regime a majority of the particles does not participate in the oscillations. Multi-level single-particle excitations are also much rarer in the Mott regime. This is also evident from Fig. 4.4(b), which shows very small peaks away from the collective mode of ω = ω 0 . However, even the very small peaks that exist 58 1 10 100 N b 0.01 0.1 1 10 n(k=0,ω) max SF MI Figure 4.7: The dependence of the maximum amplitudes of|n(k = 0,ω)| are shown as a function ofN b , using a logarithmic scale. In the superfluid regime, the amplitudes increase linearly as a function of N b . However, in the Mott regime, one observes that the maximum amplitude is drastically reduced as a function of N b . The effect of this drastic decrease is seen in the scale of the contours in the Mott regime. match closely the transition levels in the Mott regime if one compares Fig. 4.6(a) and Fig. 4.6(b). In addition, the contour plot of Fig. 4.6(a), shows horizontal lines at ω 0 and 2ω 0 . These features do not match up with Fig. 4.6(b). Instead, they are the collective modes discussed earlier, which amount to the bulk of the superfluid performing oscillations at these frequencies. The magnitude of these peaks are dominant in the plot ofn(k = 0,ω). The superfluid bulk exhibits parametric resonance much like the classical oscillator. When a classical oscillator is driven parametrically in the vicinity of frequencyω 0 , the oscillations contain frequenciesω 0 and2ω 0 , with the amplitude of oscillations at2ω 0 an order of magnitude smaller(33). This is consistent with what is observed in the contour plot. The two horizontal lines persist beyond the superfluid-to-Mott-insulator transition at N b ≥ 37, however, with smaller magnitude. In this case, part of the atoms in the 59 system participate in the Mott plateau, but the “edges” are still superfluid and therefore take part in the parametric oscillations. The bulk of the Mott insulator however does not participate in the oscillations, which is also evident from Fig 4.7. As we move deeper into the Mott regime, more and more atoms begin to join the Mott plateau, leaving a small number of particles still superfluid. This explains the dramatic drop seen in Fig 4.7 beyond the critical density ofN b = 37 bosons. 4.4 Conclusions In summary, this study suggests a new way for experimentalists to probe into the internal structure of the spectrum of 1D HCBs via parametric excitations easily brought about by a time-varying trap potential. Even for small amplitudes of the driving force, one observes two well-defined features. Firstly, collective excitations of the superfluid are seen at the driving frequency and also, to a smaller extent at twice the driving frequency. The collective mode draws analogies to collective oscillations seen in 1D HCBs(48) when a sudden shift of the trap is seen. Secondly, sharp and distinct features of the spectrum can be extracted in the vicinity of low overtones of the driving frequency,ω 0 . In the superfluid regime, these parametric excitations are seen at the driving frequency ω 0 and at 2ω 0 ,3ω 0 and 4ω 0 with successively smaller and smaller amplitudes. The extent of participation in these excitations increases linearly with increasing number of bosons, N b as long as N b ≤ 37. The multiple excitations is analogous to classical parametric resonance whereω 0 oscillations are of orderΔV 0 /V 0 and2ω 0 oscillations are of orderΔV 2 0 /V 2 0 (33). In the Mott insulator regime we see a sharp decline in parametric resonance, with most of the higher frequency excitations a few orders of magnitude smaller. Although these events are rarer, they still reflect on the internal structure of effective single-particle levels in the Mott regime. Experimentally, these might be harder 60 to observe since the total scale at which these oscillations are taking place is much smaller in the Mott regime. 61 Chapter 5 The effects of disorder on one dimensional bosons in optical lattices 5.1 Introduction Until recently, the majority of experiments with trapped ultracold atoms have been per- formed on clean systems. Indeed, the ability to create perfect (defect-free) optical lat- tices is a major advantage over condensed matter systems. On the other hand, the com- plete control over the system potentials makes it possible to introduce disorder and tune the disorder strength in a controlled fashion. The exciting possibilities of being able to study novel disorder-related phenomena such as Anderson localization and to explore novel quantum glassy phases have made the investigation of disorder in these systems an area of emerging interest. Disorder can be generated in optical lattices by exposure to speckle lasers(25; 8), adding an incommensurate lattice-forming laser(11; 49; 14), or other means(17; 19). It is thus possible to investigate different disorder-induced phenomena in a control- lable manner, in contrast e.g. to previous studies with granular superconductors or 4 He in vycor glass. The interplay between disorder and interactions in trapped Bose- Einstein condensates has recently been explored experimentally in 87 Rb, both in the continuum(36; 18; 8) and in an optical lattice(53). On the theoretical front, such sys- tems have been investigated within the frameworks of mean-field theories,(11; 53; 8; 62 Figure 5.1: Illustration of atoms in a potential trap. In the absence of interactions between the atoms, disorder leads to an immediate localization of the single-particle wave functions. 32; 16; 31), Monte Carlo simulations,(50) Bose-Fermi mapping,(13) and the transfer matrix formalism.(19) 5.2 Theoretical formulation of the problem In this chapter, we study interacting trapped bosons in one-dimensional optical lattices with a tunable random potential using a quantum Monte Carlo method. Previous inves- tigations have argued that under realistic conditions, such diagonal disorder dominates over randomness in the hopping amplitude.(11) In experiments with optical lattices, the primary source of information about the state of the system arises from the analysis of the momentum distribution function, obtained from matter-wave interference after the release of the trap and subsequent free evolution of the particles. In view of this, here we focus on the signatures of disorder on the momentum distribution. The simulations presented here predict a unique enhancement of phase coherence at small to moderate 63 strengths of disorder, when the ground state has a Mott insulating region, i.e. a Mott plateau, at the center of the trap. Since Mott insulating states are now experimentally observable, and disorder strengths are controllable by tuning the intensity of the speckle laser, these results can be straightforwardly tested using currently available experimental methods. One technical difficulty in using current setups is the large correlation length of speckle patterns (typically 8-10 times the optical lattice spacing) and the limited size of optical lattices (40-65 lattice spacings along each axis). This makes it difficult to real- ize reasonably disordered lattices. However, these difficulties are likely to be overcome in the near future with advances in realizing larger lattices and / or using more than one speckle lasers to create disordered potentials with a shorter correlation. Bosons in an optical lattice are well described by the one-band Bose-Hubbard model.(? ) In the presence of a strong periodic (lattice) potential, single-particle Wan- nier functions localized on the lattice sites form a complete basis set. Interactions are typically not strong enough to excite higher vibrational bands, justifying the single band approximation. Hence the low-energy physics of trapped bosons in a one-dimensional optical lattice can be described by the Hamiltonian H = −t L−1 X i=1 ( ˆ b † i+1 ˆ b i +h.c.)+ U 2 L X i=1 ˆ n i (ˆ n i −1)−μ 0 L X i=1 ˆ n i +V T a 2 L X i=1 (i− L 2 ) 2 ˆ n i + L X i=1 W i ˆ n i . (5.1) Here ˆ b † i ( ˆ b i ) creates (annihilates) a boson at sitei, ˆ n i = ˆ b † i ˆ b i is the number operator, and U is the strength of the repulsive on-site Hubbard interaction between bosons. The bare chemical potentialμ 0 controls the filling of the lattice,V T is the strength of the trapping potential,a denotes the lattice spacing, andW i introduces diagonal disorder in the form of a random site energy. The hopping amplitudet is set equal to unity in the simulations, 64 and all other parameters are implicitly expressed in units of t. In the experiments, W i is commonly realized by a speckle laser. Furthermore, recent measurements indicate significantly correlated disorder, with a correlation length that is typically several times longer than the optical lattice spacing, due to the diffraction-limited imaging of the speckles coming from a diffusion plate onto the trapped atoms.(36; 18) In view of this, we consider finitely correlated disorder with a random potentialW i that is distributed as W l = 0, W l W l ′ = Δδ l,l ′. (5.2) The overbar denotes averaging over disorder realizations, and Δ is the strength of the disorder correlations. The indices l are measured in units of 5a, i.e. the correlated disorder is simulated by randomly choosing clusters of 5 adjacent sites with the same on-site energy, but no correlations between the clusters. In the absence of disorder, the Bose-Hubbard model (5.1) in a trap has been exten- sively studied.(3; 60) Due to the spatially varying trapping potential, the system is never in a uniform phase. Depending on the strength of the trapping potential,V T , the on-site interaction,U, and the density of particles,n, the ground state can belong to one of two classes. At sufficiently small U/t and total density, n, it consists of a superfluid (SF) domain extending across the entire system. The site-occupation,n i varies continuously from zero at the edges to a finite maximum value at the center of the trap. The momen- tum distribution,n(q), features a sharp peak atq = 0 reflecting strong phase coherence across the entire system. The full width at half maximum of the zero-momentum peak is inversely proportional to the coherence lengthξ. In a thermodynamic SFξ diverges, whereas in the one-dimensional confined SFξ∝N andn(q = 0)∝ √ N. On the other hand, at larger U/t and sufficiently large n, the ground state contains one or more Mott insulating (MI) regions with integern i , along with domains of SF at 65 the edges of the trap and in between the MI regions (if there are more than one such regions).(3; 60) In the present study, we consider the simplest case of one MI plateau at the center of the trap and two SF domains at the edges – the local density profile consists of an extended region at the center withn i = 1 and two regions with0<n i < 1 near the trap edges. In the presence of MI region(s), the range of coherence is greatly reduced, and the momentum distribution shows a weak peak at q = 0, and a subsequent slow decrease inn(q) as a function ofq. Depending on whether there exists an MI region at the center of the trap, we find the effects of disorder are to be markedly different. 5.3 Quantum Monte Carlo results We use the stochastic series expansion quantum Monte Carlo method to simulate the disordered Bose-Hubbard model on finite-sized lattices.(? ) The density profiles,n(x), and the momentum disribution function,n(q), are measured within the grand canonical ensemble, using sufficiently large values of the inverse temperature,β, in order to obtain ground state properties. Disorder averages over 400-2000 realizations are performed, depending on the parameters. The lattice sizes are chosen to be sufficently large to ensure that the local density vanishes at the edge of the trap. We found that a lattice sizeL = 120 is necessary for large values ofU(≥ 8t) at strong disorder (W/U > 1.5), whereas L = 80 was sufficient for the other parameter sets. Results are displayed for 40 bosons in a parabolic trap withV T a 2 = 0.015t for a range of values ofU/t, ranging from the weak to strong coupling limit. Fig. 5.2 shows the consequences of a disordered site-potential on the momentum distribution functionn(q) and the real-space density profile n(x) for two values of the on-site interaction, U/t. The upper panels show the results for U/t = 2, when the ground state in the clean limit (W = 0) consists of a single SF domain extending across 66 0.001 0.01 0.1 q/2π 5 10 15 20 n(q) -40 -20 0 20 40 x/a n(x) 0.001 0.01 0.1 q/2π 0 2 4 6 8 n(q) -40 -20 0 20 40 x/a n(x) U=2t U=8t x W/U=0.0 W/U=2.5 x W/U=0.0 W/U=0.8 W/U=2.0 Figure 5.2: Effects of disorder on the momentum distributionn(q) and on the real-space density profilen(x) in the weak (upper panels) and strong (lower panels) coupling limit. The density profiles are vertically offset to enhance clarity of presentation. the trap. Adding disorder produces localization centers around which the wave func- tions get localized. This in turn reduces coherence and hence results in a smaller zero- momentum peak. At sufficiently weak disorder strengths, these effects are found to be relatively small in finite systems, and the ground state retains a predominantly SF char- acter – the deviation of the momentum distribution from its clean limit is exponentially small. However, with increasing disorder strength, the localization effects strengthen, and hence the momentum distribution starts to deviate significantly, i.e. the peak at n(q→ 0) becomes shorter and broadened. On the other hand, forU/t = 8 (lower panels), the ground state of the system, in the clean limit, has a Mott plateau at the center and the effects of disorder are considerably more interesting. The wave functions in the MI region in the absence of disorder are localized due to the strong Hubbard interactions. The onset of disorder weakens the correlation between particles, and thus leads to a partial delocalization of the strongly localized Wannier orbitals at small to moderate strengths of disorder. Phase coherence 67 hence increases as the MI region “melts”, reflected in a strongerq = 0 peak inn(q) with increasing disorder. For larger disorder potentials, the Anderson localization effects dominate, and the zero-momentum peak starts to decrease. Fig. 5.3 shows the variation of the strength of the zero-momentum peak as a function of disorder for four different values of U/t, ranging from the weak to the strong cou- pling limit. For intermediate values of the interaction strength, the peak at n(q = 0) increases at small to intermediate disorder (reflecting the melting of the central MI region), reaches a maximum, and finally decreases continuously for largerW . This is to be contrasted with the weak coupling limit, where n(q = 0) remains fairly unchanged for small W , and then decreases rapidly at larger disorder strengths. The exact posi- tion of the maximum and its value, n max (q = 0) depends on U/t, V T /t, and n. The peak position shifts to larger disorder strengths with increasing Hubbard intearction strength, indicating a competition between Mott localization and Anderson localiza- tion that depends on the ratio of W vs. U. The maximum height, n max (q = 0) and its ratio to the clean system value,n 0 (q = 0), decrease rapidly withU/t. AtU/t = 20, i.e. close to the hard-core limit, no visible increase withW is noticable inn(q = 0). What is the nature of the state that results from the onset of disorder? Does the par- tial delocalization of the MI states result in a SF or a glassy phase? To answer these questions, we analyze the effects of disorder in the absence of a trapping potential. At V T = 0, the ground state is in a pure phase, and can be conveniently characterized by measuring global observables, such as the compressibility and the stiffness. For shallow trapping potentials, the system is adiabatically connected to the periodic case, such that the domains of the different phases retain their global characteristics to a large extent. This allows us to study the effects of disorder in the pure phases and extend the results to corresponding domains in the trapped system. Let us focus on the global compressibil- ity,κ = β(hn 2 i−hni 2 ), and the stiffness,ρ s , to differentiate between various possible 68 0.0 0.5 1.0 1.5 W/U 0.8 1.2 1.6 2.0 n(q=0)/n 0 (q=0) U=5t U=8t U=12t U=20t Figure 5.3: Dependence of the zero-momentum peak height,n(q = 0), on the disorder strengthW for different values ofU/t. AtU/t = 5, the ground state consists of a single SF domain, whereas for the other values, there exist a finite MI region at the center of the trap. phases – SF, MI or Bose glass (BG). The stiffness gives the response of the system to a uniform twist in the spatial boundary as ρ s = ∂ 2 E/∂φ 2 (whereE is the ground state energy), under which the kinetic energy term in the Hamiltonian (5.1) is replaced by −t P (e −iφ ˆ b † i+1 ˆ b i +h.c.). In practice, the stiffness is simply related to the fluctuations in the winding numberw in the simulations asρ s =hw 2 i/2βL.(42) The results shown in Fig. 5.4 are for uniform SF (left panels) and MI (right panels) phases in the presence of potential disorder. The dependence of n(q = 0) on W is qualitatively similar com- pared with the case of trapped atoms. For the SF phase, n(q = 0) remains practically unchanged at smallW , and starts to decrease monotonically for larger disorder. On the other hand, in the MI phasen(q = 0) first increases and then decreases with increasing W , confirming that the disordered states are qualitatively similar to their counterparts in the presence of a trapping potential. The stiffness and compressibility data allow us to further characterize the disordered states. For the SF ground state, bothρ s andκ are 69 0.4 0.6 0.8 1.0 n(q=0)/n 0 (q=0) 1.0 1.5 2.0 2.5 n(q=0)/n(q=0) 0.0 0.5 1.0 1.5 W/U 0.0 0.5 1.0 1.5 ρ S 4∗κ 0.0 0.5 1.0 1.5 W/U 0.0 0.5 1.0 1.5 ρ S 4*κ U=5t U=9t Figure 5.4: Dependence of the momentum distribution function at q = 0, the stiffness and the compressibility of SF (left panels) and MI (right panels) phases on disorder strength. For the SF phase, the ground state is predominantly SF at small to moderate disorder and a Bose glass at large disorder. For the MI phase, the ground state is a Bose glass at all strengths of disorder studied. The strength of the zero-momentum peak varies non-monotonically, but the ground state always remains a Bose glass. finite for the range ofW over whichn(q = 0) remains close to itsW = 0 value, consis- tent with a predominantly SF character of the disordered state. At largerW , asn(q = 0) starts to deviate significantly, the compressibility remains finite, but the stiffness rapidly decreases to zero. Hence the ground state at large disorder is a compressible insulator, i.e. a Bose glass (BG). From the present data it is not clear if one needs a finite critical disorder to destroy superfluidity or simply the localization length at small disorder is larger than the system sizes considered here. For the MI phase, bothρ s andκ vanish in the clean system. With the onset of disorder, the stiffness remains zero, but the com- pressibility acquires a finite value. The disordered MI is thus a BG at all strengths of disordered considered here. 70 5.4 Concluding remarks In summary, based on Quantum Monte Carlo simulations ensembles of interacting trapped atoms are found to respond in a highly non-trivial fashion to diagonal disor- der. If the clean system consists of a single superfluid domain, the particles are simply localized by the random potential. This is reflected by progressive suppression of the q→ 0 peak in the momentum distribution function. On the other hand, if the clean sys- tem contains Mott insulating regions the response to on-site disorder is non-monotonic, and depends on the relative strength of the Hubbard interaction and the disorder poten- tial. Because of this competition there is an intermediate regime with enhanced phase coherence. An analysis of the compressibility and stiffness in the corresponding periodic systems, i.e. in the absence of a trapping potential, reveals that the resulting disordered state in either case is a Bose glass. The predicted enhancement of phase coherence by diagonal disorder suggests an interesting extension of recent experiments on optical traps in the presence of a tunable speckle laser.(18) Finally, we note that the predicted order-by-disorder phenomenon is most pro- nounced in parameter regimes which stabilize a single continuous Mott plateau with minimal superfluid tails at the trap edges. This configuration minimizes the local- ization effects within the tails and simultaneously maximizes the tunneling processes between the disorder induced superfluid lakes in the plateau region, which gives rise to the enhanced phase coherence. Accordingly, recent experiments on bichromatic opti- cal lattices at larger filling fractions do not exhibit evidence of disorder enhanced phase coherence,(14) whereas from the present study it is expected that for smaller fillings such effects should be observable. Moreover, the spatially correlated disorder provided by speckles leads to extended superfluid lakes, and hence should be more favorable for induced phase coherence than the quasi-periodic potential provided by incommensurate lasers. 71 Chapter 6 Conclusions In this thesis we have investigated various aspects of one-dimensional bosons in opti- cal traps. In particular, we have studied the properties of one-dimensional hard-core bosons. Using an exact numerical technique, the time evolution of ensembles of hard- core bosons was studied. This technique was then generalized to include time-dependent traps. The effects of parametrically driving an optical trap were studied. Finally, we examined the effects of disorder on interacting bosons with a finite on-site repulsion. Following Ref.(46), a numerical technique involving exact diagonalization after the Jordan-Wigner transformation was extended. The ground state of trapped hard-core bosons can either be a superfluid or a Mott insulator. For a fixed number of bosons, N b , we have observed that the Hamiltonian depends on one parameter, namely V 0 /t, where V 0 is the trapping potential strength and t is the hopping parameter. The energy spectrum was seen to have a linear behavior in energy index, ε(n) ∝ n, for low-lying energy states, much like the continuum quantum oscillator. For high energy states we have observed that the energy spectrum is dominated by the potential energy of the trap. This leads to a non-trivial dependence of the spectrum on n. We have shown that the cross over from a superfluid regime to a Mott insulator regime is universally seen when the potential energy of the highest occupied state is greater than2t, the maximum kinetic energy of the tight binding Hamiltonian. The possibility of other confining potentials was also studied, with particular emphasis on the gaussian confinement which has been used in recent experiments of trapped atoms. 72 In Chapter 3, we presented a study of the time evolution of one-dimensional trapped hard-core bosons. A technique to study time evolution with the help of the equal-time Green’s function was developed. The technique was then used to study the breathing mode vibration (ω B ) and sloshing mode vibration (ω S ) of one-dimensional hard-core bosons. We showed that the ratio ω B /ω S was close to the value 4. Other studies have found a similar result. This is a useful technique to determine experimentally whether a gas of one-dimensional trapped bosons was close to the hard-core constraint. We then went on to study time-dependent traps. Since the trapping potential is diag- onal, it is possible to split the Hamiltonian into a time independent part that requires diagonalization and a time-dependent part that is diagonal. A first order Trotter decom- position was used to evaluate the equal-time Green’s function. We set bounds on the strength of the time-dependent term and discussed the error due to the Trotter decompo- sition. The error was found to be linear as a function of total time. The driven oscillator was then studied using the above approximations. We observed resonant oscillations in the center of mass of the hard-core bosons, when they were driven. In Chapter 4, we continued to study time-dependent traps. In particular, paramet- rically driven traps were studied, where the time dependence was introduced in the strength term of the trapping potential, V 0 . The frequency dependence of the zero- momentum peak of the momentum distribution function was studied in detail. We found that the zero-momentum peak has frequency contributions at exactly the energy differ- ences between single particle states of the time independent Hamiltonian. This feature was seen in the superfluid and Mott insulating regimes. Apart from single particle exci- tations, collective modes were observed atω 0 and2ω 0 , whereω 0 is the natural frequency of the trap before time dependence is introduced. These collective modes, similar to classical parametric resonance, represent the bulk oscillations of the gas. The study of 73 parametrically driven traps suggests a new technique to probe into the internal struc- ture of the spectrum of one-dimensional hard-core bosons. The momentum distribution function, which is easily measurable through experiment, was shown to be a crucial tool in studying the parametric resonances of the system. The effects of disorder on interacting bosons was studied in Chapter 5. Disorder can be generated in optical lattices by exposure to speckle lasers. Ensembles of interacting trapped bosons were found to respond in a highly non-trivial fashion when a diagonal disorder was turned on. If a Mott plateau was present in the clean limit, we found that phase coherence increases as a result of turning on disorder. In this thesis we have found that systems of trapped bosons have an extremely rich phenomenology, including collective effects and exotic quantum phases. Moreover, they can serve as model systems to gain insight into the mysterious world of interacting quantum particles. 74 Bibliography [1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995). [2] for e.g., Nir Bar-Gill, Eitan E. Rowen, Gershon Kurizki, Nir Davidson, arXiv:0902.2721 [3] G. G. Batrouni, et al., Phys. Rev. Lett. 89, 117203 (2002); P. Sengupta, M. Rigol, G. G. Batrouni, P. J. H. Denteneer, and R. T. Scalettar, Phys. Rev. Lett. 95, 220402 (2005). [4] I. Bloch, Nature Physics 1, 23 - 30 (2005) [5] I. Bloch, Phys. World 17, 2529 (2004) [6] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995). [7] F. Chevy, V . Bretin, P. Rosenbusch, K.W. Madison, and J. Dalibard, Phys. Rev. Lett. 88, 250402 (2002). [8] D. Cle´ ement, A. F. Var´ on, M. Hugbart, J. A. Retter, P. Bouyer, L. Sanchez-Palencia, D. M. Gangardt, G. V . Shlyapnikov, and A. Aspect, Phys. Rev. Lett., 95, 170409 (2005). [9] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995). [10] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999). [11] B. Damski, J. Zakrzewski, L. Santos, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 91, 080403 (2003). 75 [12] B. DeMarco, and D. S. Jin, Science 285, 1703 (1999). [13] A. De Martino, M. Thorwart, R. Egger, and R. Graham, Phys. Rev. Lett., 94, 060402 (2005). [14] L. Fallani, J. E. Lye, V . Guarrera, C. Fort, and M. Inguscio, cond-mat/0603655. [15] Feynman, R. P.,Opt. News 11, 1120 (1985) [16] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989). [17] R. Folman, P. Kr¨ uger, J. Schmiedmayer, J. Denschlag, and C. Henkel, Adv. At. Mol. Opt. Phys. 48, 263 (2002). [18] C. Fort, L. Fallani, V . Guarrera, J. E. Lye, M. Modungo, D. S. Wiersma, and M. Inguscio, Phys. Rev. Lett. 95, 170410 (2005). [19] U. Gavish, and Y . Castin, Phys. Rev. Lett., 95, 020401 (2005). [20] M. D. Girardeau, E. M. Wright, and J. M. Triscari, Phys. Rev. A 63, 033601 (2001) [21] Goerlitz et al, A., Phys. Rev. Lett. 87, 130402 [22] M. Greiner, I. Bloch, O. Mandel, T. W. Hansch and T. Esslinger, Phys. Rev. Lett. 87, 160405 (2001). [23] M. Greiner, O. Mandel, T. Esslinger, T.W. HŁnsch I. Bloch, Nature 415, 3944 (2002). [24] Grimm, R., Weidem¨ uller, M. Ovchinnikov, Yu. B, Adv. At. Mol. Opt. Phys, 42, 95-170 (2000). [25] P. Horak, J. Y . Courtois, and G. Grynberg, Phys. Rev. A 58, 3953 (2000). [26] P. Jordan and E. Wigner, Z. Phys. 47, 631 (1928). [27] R. M. Kalas, and D. Blume, cond-mat/051203. [28] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature (London) 440, 900 (2006). [29] C. Kollath, A Iucci, T. Giamarchi, W. Hofstetter, U. Schollwoeck, Phys. Rev. Lett. 97, 050402 (2006). [30] C. Kollath et al., Phys. Rev. Lett 98, 180601 (2007). [31] K. V . Krutitsky, A. Pelster, and R. Graham, New J. Phys. 8, 187 (2006). 76 [32] R. C. Kuhn, C. Miniatura, D. Delande, O. Sigwarth, and C. A. M¨ uller, cond- mat/0506371. [33] L.D. Landau and E.M. Lifshitz, Mechanics 27 (1960). [34] A. Lenard, J. Math. Phys 1, 516 (1960). [35] Xia-Ji Liu, Peter D. Drummond, Hui Hu, Phys. Rev. Lett. 94, 136406 (2005). [36] J. E. Lye, L. Fallani, M. Modungo, D. S. Wiersma, C. Fort, anf M. Inguscio, Phys. Rev. Lett. 95, 070401 (2005). [37] M.-O. Mewes, M. R. Andrews, N. J. van Druten, D.M. Kurn, D. S. Durfee, C.G. Townsend, and W. Ketterle, Phys. Rev. Lett. 77, 988 (1996). [38] H. Moritz et al., Phys. Rev. Lett. 91, 250402 (2003). [39] K. M. O’Hara, et al., Science 298, 2179 (2002);De M. Greiner, et al., Nature 426, 537 (2003); S. Jochim, et al., Science 302, 2101 (2003). [40] M. Olshanii, Phys. Rev. Lett., 81, 938 (1998) [41] B. Paredes, A. Widera, V . Murg, O. Mandel, S. Folling, I. Cirac, G. V . Shlyapnikov, T. W. Hansch, and I. Bloch, Nature (London) 429, 277 (2004). [42] E.L. Pollock and D.M. Ceperly, Phys. Rev. B 36, 8343 (1987). [43] D. Pritchard, Phys. Rev. Lett., 51, 1336 (1983). [44] M. Rigol, V . Dunjko, M. Olshanii, Nature 452, 854-858 (2008). [45] M. Rigol, A. Muramatsu, Mod. Phys. Lett. B 19, 861 (2005). [46] M. Rigol, A. Muramatsu, Phys. Rev A 72, 013604 (2005). [47] M. Rigol, A. Muramatsu, Phys. Rev. Lett. 93, 230404 [48] M. Rigol, V . Rousseau, R. T. Scalletar, R. R. P. Singh, Phys. Rev. Lett 95 110402 (2005) [49] A. Sanpera, A. Kantian, L. Sanchez-Palencia, J. Zakrzewski, and M. Lewenstein, Phys. Rev. Lett., 93, 040401 (2004). [50] R. T. Scalettar, G. G. Batrouni, and G. T. Zimanyi, Phys. Rev. Lett. 66, 3144 (1991). [51] Schollwoeck and White, Effective models for low-dimensional strongly correlated systems, AIP, Melville, New York (2006) 77 [52] Schreck, F. et al., Phys. Rev. Lett. 87, 080403 (2001) [53] T. Schulte, S. Drenkelforth, J. Kruse, W. Ertmer, J. Arlt, K. Sacha, J. Zakrzewski, and M. Lewenstein, Phys. Rev. Lett. 95, 170411 (2005). [54] T. St¨ oferle, et al., Phys. Rev. Lett. 92, 130403 (2004). [55] S. Stringari, Phys. Rev. Lett. 77, 2360 (1996). [56] M. Suzuki, Prog. Theoret. Phys. 56, 1454 (1976). [57] B. L. Tolra et al., Phys. Rev. Lett. 92, 190401 (2004). [58] H. F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959). [59] see for e.g., J. H. Wilkinson, The algebraic eigenvalue problem [60] S. Wessel et al., Phys Rev. A 70, 053615 (2005). 78 Appendix A Eigenvalues of a Hamiltonian,H with hopping,t = 1 and trapping potential,V(x) = 0 Following the discussion in Chapter 2, we study the numerically acquired eigenval- ues of a tridiagonal matrix H with diagonal entries = 0 H =                0 −1 0 0 ... 0 −1 0 −1 0 ... 0 . . . . ... 0 . . . . −1 0 ... ... 0 −1 0 −1 ... ... ... 0 −1 0                . (A.1) The significance of this matrix is seen as a limiting case of V(x) = 0 or when α→∞, inV α (x) =V 0 (x−x 0 ) α . 79 We need to find the determinant, det(H−λI) in order to get the eigenvalues (λ). One can define D k (λ) = det                −λ −1 0 0 ... 0 −1 −λ −1 0 ... 0 . −1 . . ... 0 . . . . −1 0 ... ... 0 −1 −λ −1 ... ... ... 0 −1 −λ                kxk , (A.2) fork = 1,2...,N. One can further define, D 0 (λ) = 1. (A.3) Then,D k follows the recurrence relation D k (λ) =−λD k−1 (λ)−D k−2 (λ). (A.4) The characteristic polynomials are associated with the continued fractions, namely, −λ− 1 −λ− 1 −λ− 1 −λ− 1 −λ−... = 0. (A.5) 80 The smallest and largest eigenvalues approach the value∓2 and the general eigen- value spectrum follows(59) λ k =−2cos πk n+1 . (A.6) LAPACK can very easily diagonalize this tridiagonal matrix whose eigenvalues are related to the tight binding dispersion relation, E TB =−2tcos(ka), (A.7) with eigenvectors then become |ki = N X j=1 e ijka |ji, (A.8) where|ji is the original Wannier function basis in which the Hamiltonian is written. N λ 1 λ N −2cos nπ n+1 2 -1.00000000000000 1.00000000000000 1.0000 4 -1.61803398874989 1.61803398874989 1.6180 8 -1.87938524157182 1.87938524157182 1.8794 16 -1.96594619936780 1.96594619936780 1.9659 32 -1.99094384514617 1.99094384514617 1.9909 64 -1.99766445366465 1.99766445366465 1.9977 128 -1.99940693969028 1.99940693969028 1.9994 256 -1.99985057333946 1.99985057333946 1.9998 512 -1.99996249720310 1.99996249720310 1.9999 1024 -1.99999060597580 1.99999060597580 2.0000 2048 -1.99999764919967 1.99999764919967 2.0000 81 0 5 10 15 20 25 30 35 40 45 50 N -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 λ λ 1 λ N Figure A.1: The convergence of the smallest and largest eigenvalues converging very quickly to∓2. Here we have concentrated only on data forN ≤ 50. ForN > 30, we see that the smallest and largest eigenvalues approach very close to the value of∓2. In the table, we can see data for uptoN = 2048. The first column is the smallest eigenvalue, the second column is the largest eigenvalue and the third column is the evaluation of the expression−2cos nπ n+1 . 82

Abstract (if available)

Abstract Following the recent advances in controlling ultracold quantum gases that have led to the realization of boson and fermion condensation, we present computational studies of one-dimensional bosons on optical lattices. These systems have proven to be a great tool to study exciting new frontiers of condensed matter, from exotic phase transitions to nonequilibrium phenomena.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Quantum phase transitions in disordered antiferromagnets

PDF

Phase diagram of disordered quantum antiferromagnets

PDF

Modeling graphene: magnetic, transport and optical properties

PDF

Out-of-equilibrium dynamics of inhomogeneous quantum systems

$X-ray coherent diffractive Imaging of doped quantum fluids$

PDF

X-ray coherent diffractive Imaging of doped quantum fluids

Asset Metadata

Creator Raghavan, Aditya (author)

Core Title Properties of hard-core bosons in potential traps

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 08/08/2009

Defense Date 06/05/2009

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag Bose Hubbard model,bosons,condensed matter physics,hard core bosons,OAI-PMH Harvest,optical traps,Physics,quantum gases,time dependent Hamiltonian,ultracold atoms

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Haas, Stephan W. (committee chair), Bickers, Nelson Eugene, Jr. (committee member), Dappen, Werner (committee member), Kresin, Vitaly V. (committee member), Wang, Chunming (committee member)

Creator Email aditya.raghavan@gmail.com,adityara@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2549

Unique identifier UC1495631

Identifier etd-RAGHAVAN-3065 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-174843 (legacy record id),usctheses-m2549 (legacy record id)

Legacy Identifier etd-RAGHAVAN-3065.pdf

Dmrecord 174843

Document Type Dissertation

Rights Raghavan, Aditya

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu