University of Southern California Dissertations and Theses

Chaos and thermalization in the one-dimensional Bose-Hubbard model in the classical-field approximation

(USC Thesis Other)

Chaos and thermalization in the one-dimensional Bose-Hubbard model in the classical-field approximation

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content CHAOS AND THERMALIZATION IN THE ONE-DIMENSIONAL BOSE-HUBBARD MODEL IN THE CLASSICAL-FIELD APPROXIMATION by Amy C. Cassidy 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) May 2009 Copyright 2009 Amy C. Cassidy Acknowlegements There are a many people without whose support over the last six years I would not have completed this project. First of all, I would like to thank my advisor, Maxim for all that he has taught me, for his commitment to my development and for space to grow. Others whom I would like to acknowledge are: Vanja, for asking questions; Tameem, Richard and Yi-Chun, my classmates who got me through the first year; Ilya, for guidance at a critical time; Noah and Marcos, for physics conversations and tennis matches; Katie, for countless hours working together to encourage more women to pursue physics; Sam, for perspective; Courtney, for a lifelong friendship; Karla, for reminding me to stay bal- anced; Stephan, for getting me started and continued support; Werner and the depart- ment at USC, for flexibility; Bala and the department at UMass Boston, for welcoming me these last two years. Finally I would like to thank my family: my father who encour- aged my intellectual curiosity; my mother, who taught me to be organized and persevere; Sara and Andrew for many words of wisdom; and Linda for encouragement and support. ii Table of Contents Acknowlegements ii List of Tables vi List of Figures vii Abstract x Chapter 1 : Introduction 1 1.1 Introduction to Dynamical Systems . . . . . . . . . . . . . . . . . . 1 1.2 Dynamical Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Integrable Systems . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 Chaos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Statistical Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.1 Statistical Ensembles . . . . . . . . . . . . . . . . . . . . . 7 1.4 FPU Model and Results . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4.1 KAM Theorem . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4.2 Solitons and the Korteweg-de-Vries Equation . . . . . . . . 14 1.4.3 Chirikov Criterion . . . . . . . . . . . . . . . . . . . . . . . 15 1.5 Relationship between Chaos and Thermalization . . . . . . . . . . . 19 1.5.1 Thermalization in Classical Field Models . . . . . . . . . . 21 1.6 Quantum Degenerate Gases - Ultracold Atoms . . . . . . . . . . . . 23 1.6.1 Bose-Einstein Condensation . . . . . . . . . . . . . . . . . 24 1.6.2 Bosons in Optical Lattices . . . . . . . . . . . . . . . . . . 25 1.6.3 Atom Chips . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.6.4 Classical Field Model of Bose Gas . . . . . . . . . . . . . . 28 1.6.5 Chaos and Integrability in Quantum Systems . . . . . . . . 30 1.6.6 Constrained equilibrium . . . . . . . . . . . . . . . . . . . 32 1.7 Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Chapter 2 : Thermalization and Chaos in the 1D Bose-Hubbard Model 35 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 iii 2.2 BHM: Hamiltonian and Equations of Motion . . . . . . . . . . . . . 35 2.2.1 Real-Space Hamiltonian . . . . . . . . . . . . . . . . . . . 36 2.2.2 Momentum-Space Representation . . . . . . . . . . . . . . 36 2.2.3 Action-Angle Representation . . . . . . . . . . . . . . . . . 38 2.2.4 Validity of the Classical-Field Theory . . . . . . . . . . . . 39 2.2.5 Nearby Integrable Models . . . . . . . . . . . . . . . . . . 40 2.3 Time Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4 Chaos: Calculating Lyapunov Exponents . . . . . . . . . . . . . . . 43 2.5 Thermalization: Calculating Spectral Entropy . . . . . . . . . . . . 44 2.5.1 Conserved Quantities . . . . . . . . . . . . . . . . . . . . . 44 2.5.2 Hartree-Fock . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.5.3 Minimization of the Grand Potential . . . . . . . . . . . . . 46 2.5.4 Spectral Entropy . . . . . . . . . . . . . . . . . . . . . . . 49 2.6 Results: N=21 Sites, Three-Mode Initial Conditions . . . . . . . . . 50 2.6.1 Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.7 Chaos Threshold for Different Lattice Sizes . . . . . . . . . . . . . 55 2.8 Two Parametric Theory of the Chaos Threshold . . . . . . . . . . . 57 2.8.1 Thermalization Times and Slow Relaxation . . . . . . . . . 62 Chapter 3 : Resonance Model and Failure of Chirikov’s Criterion 66 3.1 Perturbation Theory Study of BHM . . . . . . . . . . . . . . . . . . 66 3.1.1 Dynamics of an Initially Unpopulated Mode . . . . . . . . 67 3.1.2 Nonlinear Frequencies in the Real-Space Dynamics . . . . . 71 3.2 Chirikov’s Criterion . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.2.1 Classification of Resonances . . . . . . . . . . . . . . . . . 74 3.2.2 First-Order Resonances in BOA . . . . . . . . . . . . . . . 76 3.2.3 1-R Resonances in BOA.Ω p =Ω q . . . . . . . . . . . . . . 84 3.2.4 Second Order Resonances in BOA . . . . . . . . . . . . . . 85 3.3 Failure of Chirikov’s Criterion . . . . . . . . . . . . . . . . . . . . 92 3.3.1 Thermodynamic Limit . . . . . . . . . . . . . . . . . . . . 96 3.3.2 Continuous Limit . . . . . . . . . . . . . . . . . . . . . . . 96 Chapter 4 : Conserved Quantities of the Ablowitz-Ladik Lattice 98 4.1 Integrable Discrete Nonlinear Schr¨ odinger (IDNLS) Equation . . . 99 4.2 Conserved Quantities of the AL equation . . . . . . . . . . . . . . . 101 4.2.1 AL equation . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.2.2 Direct Scattering Problem . . . . . . . . . . . . . . . . . . 107 4.2.3 Conservation Laws . . . . . . . . . . . . . . . . . . . . . . 111 4.2.4 Conserved Quantities . . . . . . . . . . . . . . . . . . . . . 116 4.3 AL and BHM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 iv Chapter 5 : Outlook and Conlusion 120 5.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.2 Open Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Bibliography 123 Appendix A : Thermodynamic Distribution within Hartree-Fock 128 A.1 Hartree Fock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 A.2 Minimization of the Grand Potential . . . . . . . . . . . . . . . . . 131 v List of Tables 2.1 List of Parameters and Variables . . . . . . . . . . . . . . . . . . . 38 2.2 Relationship between Parameters . . . . . . . . . . . . . . . . . . . 40 2.3 Thermodynamic Limit and Scaling . . . . . . . . . . . . . . . . . . 54 3.1 Resonance Parameters . . . . . . . . . . . . . . . . . . . . . . . . 75 3.2 Bare Detuning: Generic Values . . . . . . . . . . . . . . . . . . . . 76 vi List of Figures 1.1 Time dependence of normal modes for cubic forces with β = 8. (Fermi et al., 1974) . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 Space-time diagram of soliton trajectories. u 0 = cos(πx). δ = 0.022. T R is the recurrence time (Zabusky and Kuskal, 1965). . . . . 15 1.3 Schematic of Chirikov’s criterion of the onset of chaos for the kicked rotor. J,ψ and J 0 ,ψ 0 are coordinates generated by two different canonical transformations. ΔJ and ΔJ 0 are width of the separa- tricies. J−J 0 represents the distance between the resonances . . . . 17 1.4 Velocity distribution of rubidium atoms. Left: prior to condensa- tion. Center: just after condensation. Right: after further cooling (Anderson et al., 1995) . . . . . . . . . . . . . . . . . . . . . . . . 24 1.5 f(p ex ), momentum distribution for (a)γ 0 =4,τ=34 ms. t blue = 15τ, t red = 30τ (b)γ 0 =1,τ=13 ms.t blue = 15τ,t red = 40τ (c)γ 0 =0.62, τ=13 ms. t blue = 15τ, t red = 40τ. Green is the initial momen- tum distribution averaged over the first period. (Kinoshita et al., 2006) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.1 Time evolution of the wavefunction at the center of the box for a typical run for κ = 0.009,0.09,0.36,0.9 with identical initial conditions. (a) Time dynamics of the real part of the wavefunc- tion, Re ψ x=0 (t). (b) Frequencies of the real-space wavefunction, |FT[ψ x=0 (t)]| 2 . For κ = 0.009,0.09 the descendants of the unper- turbed frequencies generated by the first, “integrable”term of the Hamiltonian (2.5) are labeled. . . . . . . . . . . . . . . . . . . . . 42 2.2 Initial, final and Hartree-Fock thermal momentum distributions for κ = 0.09,0.45,1.8, starting from the same initial state. N s = 21. The initial state is a representative state and the final state is time- averaged. ² T is the total energy per particle. . . . . . . . . . . . . . 49 vii 2.3 Ensemble-averaged finite-time maximal Lyapunov exponent,λ, and normalized spectral entropy, η, as a function of the nonlinearity, κ. N s =21. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.4 Normalized spectral entropy of the final time-averaged state versus finite-time maximal Lyapunov exponent for each of the 100 initial condition used to compute the averaged value for κ = 0.36, 0.54, 0.72,0.9. N s =21. . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.5 Fluctuations in kinetic energy. Normalized standard deviation/mean for N s =21, 42, 63, 84 lattice sites and N 0 = 21. Data points are for five sample runs with equivalent initial conditions for different lattice sizes andκ=2.69. . . . . . . . . . . . . . . . . . . . . . . . 55 2.6 Averaged Finite Time Lyapunov exponent, λ[J], for three different system sizes, N s = 11,21,41. For each κ, the same energy-per- particle was used for each lattice size. The error bars represent one standard deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.7 (a) Contour lines of the averaged FTMLE versus the nonlinearity, κ, and energy-per-particle, ² T = (H − H 0 )/N a , where H is the Hamiltonian (2.5), and H 0 = −2J + (1/2)μ 0 is the ground state value of H. The solid contour line corresponds to λ c = 0.01. The diagonal solid line represents the set of energies and nonlinearities used in Fig. 2.6. (b) Contour lines of the averaged normalized spec- tral entropy versus the nonlinearity,κ and energy-per-particle. Solid contour lines correspond to η = 0.68 andη = 0.36. For reference, the threshold line from the FTMLE in (a) is plotted (dashed line). N s =11. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 2.8 Data points used for interpolation for contour plots in Fig. 2.7 super- imposed on data for Lyapunov exponent. . . . . . . . . . . . . . . 59 2.9 Contour lines of the standared deviation of the NSE versus the non- linearity, κ and energy-per-particle. Solid contour line corresponds to atη =0.36 andη =0.68 from Fig. 2.7(b). N s =11. . . . . . . . 61 2.10 Averaged normalized spectral entropy, η, for three different system sizes,N s =11,21,41. For eachκ, the same energy-per-particle was used for each lattice size. . . . . . . . . . . . . . . . . . . . . . . . 62 2.11 Sample time dependence of normalized spectral entropy, η. κ = 0.09,² T = 0.081J,N s = 21. Inset: Initial (dashed red), final (solid blue) and thermal (dotted black) momentum distributions. . . . . . 63 viii 2.12 Sample time dependence of normalized spectral entropy,η, for two different initial states. κ = 0.54,² T = 0.19J for both. N s = 21. Inset: Initial (dashed red), final (solid blue) and thermal (dotted black) momentum distributions. . . . . . . . . . . . . . . . . . . . 64 3.1 Frequencies of the real-space wavefunction at the center of the box, |FT[ψ x=0 (t)] 2 | forξ =0.1,1,2,3. . . . . . . . . . . . . . . . . . . 72 3.2 Real and imaginary parts of the fixed points of first-order Hamilto- nian (3.28). Δ 01 =−4, ¯ I 1 =1/3. . . . . . . . . . . . . . . . . . . 80 3.3 Action-angle phase-space plots for variables ˜ I n , ˜ θ n for the first-order resonant Hamiltonian (3.28). (n+2,n+1)→ (n,n+3) (Δ 01 = −4). ¯ I 1 = 1/3. (a) ξ = 1 (b) ξ = 3.5 (c) ξ = ξ c1 = 4.3 (d) ξ = 7 > ξ c2 . The black contour line corresponds to ˜ h n = 0 and the white contour line to ˜ h n ( ˜ I n , ˜ θ n )= ˜ h n ( ˜ I ∗ n3 ,−π). . . . . . . . . . . . 81 3.4 Real and imaginary parts of the solutions of ˜ h n = 0 for the first- order Hamiltonian (3.28). Δ 01 = −4 (n + 1,n + 2) → (n,n + 3), ¯ I 1 = 1/3. For ξ ≤ ξ c1 , ˜ I sep2 gives the separatrix for the the resonance at( ˜ I ∗ n2 , ˜ θ n =±π). Forξ ≥ ξ c1 , ˜ I sep1 gives the separatrix for the the resonance at( ˜ I ∗ n1 , ˜ θ n =0). . . . . . . . . . . . . . . . . 83 3.5 Phase-space plots for the second order Hamiltonian. (a)ξ =0.5(κ= 0.16 forN s = 11). There is no resonance in the range 0≤ ˜ I n ≤ 1 for ξ < ξ c = 1. (b) ξ = 1.5 (κ = 0.49 forN s = 11). ¯ I 2 = 1/3, Δ 02 = 2. Forξ > ξ c the resonant value is given by (3.72) and the separatrix with is given by (3.73). . . . . . . . . . . . . . . . . . 91 3.6 Resonant values and separatrix width for second-order resonance for Δ 02 =2[(n+1,n−1)→(n,n)] and forΔ 02 =8[(n+2,n−2)→ (n,n)]. ¯ I = 0.33. Vertical lines: Δ 20 = 2: ξ c1 = 1/3 ¯ I 2 = 1,Δ 02 = 8: ξ c2 = 4/3 ¯ I 2 = 4. The upper x-axis gives the corresponding values ofκ forN s =11. . . . . . . . . . . . . . . . . . . . . . . . 92 ix Abstract One of the fundamental assertions of statistical mechanics is that the time average of a physical observable is equivalent to the average over phase space, with microcanonical measure. A system for which this is true is said to be ergodic and dynamical properties can be calculated from static phase-space averages. Dynamics of a system which is fully integrable, that is has as many conserved quantities as degrees of freedom, is constrained to a reduced phase space and thus not ergodic, although it may relax to a modified equilibrium. In this thesis, we present a comprehensive study of chaos and thermalization of the one-dimensional Bose-Hubbard Model (BHM) within the classical field approximation. This model describes the dynamics of quantum degenerate gases in a lattice for sufficient occupation of every momentum mode and weak two-body scattering, and is of interest because of experimental advances of cooling and trapping alkali atoms in the quantum degenerate regime. We study the chaos and its relation to thermalization. Two quantitative measures are compared: the ensemble-averaged Finite-time Maximal Lyapunov exponent, a mea- sures of chaos and the normalized spectral entropy, a measure of the distance between the numerical time-averaged momentum distribution and the one predicted by thermo- dynamics. A threshold for chaos is found, which depends on two parameters, the nonlin- earity and the total energy-per-particle. Below the threshold, the dynamics are regular, x while far above the threshold, complete thermalization is observed, as measured by the normalized spectral entropy. We study individual resonances in the Bose-Hubbard model to determine the cri- terion for chaos. The criterion based on Chirikov’s method of overlapping resonances diverges in the thermodynamic limit, in contrast to the criterion parameters inferred from numerical calculations, signifying the failure of the standard Chirikov’s approach. The Ablowitz-Ladik lattice is one of several integrable models that are close to the BHM. We outline the method of Inverse Scattering Transform and generate the integrals of motion of the Ablowitz-Ladik lattice. Furthermore, we discuss the possible role of these quantities in the relaxation dynamics of the BHM. xi Chapter 1 Introduction 1.1 Introduction to Dynamical Systems The study of thermalization in nonlinear systems dates back to the early numerical stud- ies of coupled anharmonic oscillators by Fermi, Pasta, and Ulam (FPU). At the time it was expected that any amount of nonlinearity, no matter how small, in a system with many degrees of freedom would lead to thermal behavior. Quite unexpectedly, thermalization was not observed. Later work found a threshold in the nonlinear cou- pling strength, with energy equipartition occurring above the threshold. The absence of thermalization for small nonlinearities has been explained in terms of the presence of a threshold for the onset of widespread chaos determined by Chirikov’s criterion of overlapping resonances and also in terms of the closeness to a fully integrable model, the Korteweg-de-Vries equation. Additional studies on thermalization and approach to equilibrium have been carried out in many classical field theories. Fermi, Pasta, Ulam (FPU) and Tsingou performed one of the first numerical exper- iments in the 1950’s (Fermi et al., 1974). Their discovery led to significant work and important discoveries in the field of nonlinear dynamics. The motivation for the experi- ment was to study the approach to equilibrium of a nonlinear system. The contemporary expectation was that an system with a large number of degrees of freedom would exhibit thermal behavior in the presence of any non-linearity, no matter how small. In particular 1 they were looking to observe the Fourier heat law for conduction. They numerical inte- grated a system of one-dimensional anharmonically coupled oscillators with quadratic and cubic forces. The results found by FPU were completely suprising. In order to understand the prevalent view at the time and the significance of their results, let us take some time to review important concepts in Hamiltonian systems and statistical mechanics. 1.2 Dynamical Systems Given a system whose state is completely known at some time t 0 , what can be said about the state at some future timet>t 0 ? What about its state at some time in the past, t<t 0 ? Classical mechanics addresses these questions. We begin with a short review of Hamiltonian dynamics, referring to Tabor (1989). Hamiltonian Systems Consider a generic Hamiltonian of an N-particle system, H({q i ,p i }) where {q i } are the coordinates of the N-particles and {p i } are the corre- sponding momentum. The equations of motion are given by ˙ q i =− ∂H ∂p i , ˙ p i = ∂H ∂q i . (1.1) These equations, called the Hamilton’s equations of motion, uniquely determine the time evolution of the state, from the initial state, which is specified by a complete set of values{q i ,p i }. Note that canonical coordinates and momentum satisfy Liouville’s theorem, so that a volume element in phase space moves like an incompressible fluid under the Hamiltonian flow. 2 Poisson Brackets For any dynamical quantitiesf,g in a system with canonical coor- dinates{q i ,p i }, the Poisson brackets are defined as {f,g}≡ X i ∂f ∂p i ∂g ∂q i − ∂f ∂q i ∂g ∂p i . From this one can write the time dependence of a functionf(p,q,t) as df dt = X i µ ∂f ∂q i ∂q i ∂t + ∂f ∂p i ∂p i ∂t ¶ + ∂f ∂t = X i µ ∂f ∂q i ∂H ∂p i − ∂f ∂p i ∂H ∂q i ¶ + ∂f ∂t ={H,f}+ ∂f ∂t . Thus any quantity that does not explicity depend on time is a constant of motion if the Poisson bracket of that quantity vanishes. The Poisson brackets are antisymmetric and satisfy the identity {f,gh}=g{f,h}+{f,g}h. Canonical Transformations Oftentimes it is more convenient to work in one coor- dinate system than another. A canonical transformation is one in which the canonical form of Hamilton’s equations is preserved. The phase volume is preserved under a canonical transformation. Thus the Jacobian, must be unity. Consider a transforma- tion from one set of phase-space variables(p i ,q i ) to a new set of variables(P i ,Q i ) the preservation of phase volume is expressed as Z dp 1 dq 1 Z dp 2 dq 2 ... Z dp N dq N = Z dP 1 dQ 1 Z dP 2 dQ 2 ... Z dP N dQ N 3 and the Jacobian must satisfy the condition ∂(p 1 ,...,p N ,q 1 ,...,q N ) ∂(P 1 ,...,P N ,Q 1 ,...,Q N ) = ∂(P 1 ,...,P N ,Q 1 ,...,Q N ) ∂(p 1 ,...,p N ,q 1 ,...,q N ) =1 (Tabor, 1989). One of the representations of the canonical transformation, is through the type 2 gen- erating functionΦ(θ, ˜ I), which transformsI,θ→ ˜ I, ˜ θ. The Hamiltonian is transformed according to e H( ˜ I, ˜ θ)=H(I,θ)+ ∂Φ ∂t I = ∂Φ ∂θ ˜ θ = ∂Φ ∂ ˜ I . (1.2) 1.2.1 Integrable Systems An integrable system has as many independent constants of the motion as degrees of freedom. Stated another way, a systems with n degrees of freedom is completely inte- grable is there existn integrals of motion which are in involution. {F i ,H}=0, i=1···n, {F i } are constants of the motion {F i ,F j }=0, i=1···n,j =1···n {F i } are in involution Ford makes the distinction that these must be in fact well-behaved integrals of motion (Ford, 1992). All one dimensional systems are integrable (Tabor, 1989). Exam- ples of integrable systems with many degrees of freedom include the Toda lattice (infinite-dimensional countable), the Korteweg-de-Vries equation, sine-Gordon model and continuous nonlinear Schr¨ odinger equation (infinite-dimensional continuous). 4 1.2.2 Chaos The term chaos is commonly used to describe systems where the motion appears ran- dom, erratic or unpredictable. It is sometimes called ”deterministic chaos” to emphasize that it refers to dynamical systems in which the motion is governed by deterministic equations of motion. A key feature of chaotic systems is extreme sensitivity to initial conditions. Imagine two identical systems that are governed by the same equations of motion and with initial conditions that are only slightly different. In a chaotic system, the initial different will grow rapidly and after some time the states of the two systems will be entirely different from one another. In contrast, in a regular system, the differ- ence will grow slowly, so that two remain highly correlated. The difference between chaotic and regular motion is not only quantitative, but qualitative. The very functional form of the divergence is different: chaotic trajectories diverge exponentially while reg- ular trajectories diverge linearly. Lyapunov Exponents A chaotic system is characterized by local instability, specifi- cally by exponential divergence of trajectories neighboring in phase space. This rate of divergence of is measured by the Lyapunov exponent. The maximal Lyapunov Exponent (MLE) of a system is defined as λ(x 0 )= lim t→∞ lim e x 0 →x 0 ln ke x(t)−x(t)k ke x 0 −x 0 k where x(t) ande x(t) are two phase space trajectories (Tabor, 1989). Chaotic motion is characterized by positive MLE, λ > 0, in which case trajectories that are initially close in phase space will diverge exponentially. On the other hand, a zero MLE,λ = 0, indicates regular motion and linear divergence. The Lyapunov exponent is a function of an initial state and is thus a measure of local chaos. A system may have a mixed 5 phase space, consisting of both chaotic and regular regions (Zaslavsky, 1999). Thus initial states in the regular regions will have zero Lyapunov exponent while initial states in chaotic regions have a positive Lyapunov exponent. As a system transitions from regular to globally chaotic behavior, the volume of phase space corresponding to chaotic regions grows while the regular regions shrink. When the chaotic regions dominate the phase space, there can still exist regular regions, called islands of stability, which are surrounded by the chaotic sea. 1.3 Statistical Mechanics Thermodynamics is a phenomenological theory that has very successfully described the equilibrium properties of isolated systems with many degrees of freedom. In this sec- tion we give an overview of basic concepts of statistical mechanics. The field arose in the study of macroscopic systems of atoms with N 10 23 particles and correspondingly large volumes. One of the great success of thermodynamics is the Maxwell-Boltzmann distribution is for ideal gases. It can be derived in two ways: (1) through the Boltzmann transport equation and (2) through the most probable distribution. Systems are described not by position and momentum of ever particle in the system, but by probability distri- bution functions from which macroscopic properties such as pressure and temperature and volume can be determined. The aim of statistical mechanics is to derive the laws of thermodynamics from molecular dynamics. For a system with N particles with canonical coordinates q 1 ,q 2 ,...,q 3N and conju- gate momentaq 1 ,q 2 ,...,q 3N , the6N dimensional space spanned by these coordinates is the phase space or Γ space. A representative point is a point in this Γ space that spec- ifies the position and momentum of each of the N particles. The representative point will also be called a microstate. For each microstate, there is a corresponding macrostate 6 that specifies properties of the system such as density, temperature, and pressure. Many microstates and in fact, an infinite number, correspond to the same macrostate. The dynamics of the individualN particles is governed by the Hamiltonian,H(p,q) where(p,q)=(p 1 ,p 2 ,...,p 3N ,q 1 ,q 2 ,...,q 3N ). The equations of motion are given by the Hamilton’s equations of motion, ˙ q i = ∂H ∂p i ˙ p i =− ∂H ∂q i . These equations become quickly intractable for systems with large numbers of degrees of freedom. Instead systems will be studies by various ensembles which give the distri- bution of representative points in phase space. 1.3.1 Statistical Ensembles A statistical ensemble is a collection of microstates which correspond to the same macrostate. Geometrically the ensemble can be described by the distribution of rep- resentative points in theΓ space with the density distribution functionρ(p,q,t) which is defined so that ρ(p,q,t)d 3N pd 3N q is the number of representative points in phase space volumed 3N pd 3N q at time t. Liouville Theorem dρ dt + 3N X i=1 µ dρ dp i ˙ p i + dρ dq i ˙ q i ¶ =0 The interpretation of the Liouville theorem is that the distribution of representative points in phase space move like an incompressible fluid. 7 Ensemble averages The ensemble average of an observableO is give by hOi= R d 3N pd 3N qO(p,q)ρ(p,q,t) R d 3N pd 3N qρ(p,q,t) The time dependence ofO comes from the time-dependence ofρ. Postulate of Equal a Priori probability For a system in thermodynamic equilibrium, it is equally likely to be in any microstate (of the same phase volume) that satisfies the macroscopic conditions of the system. For an isolated system, the distribution is described by the microcanonical ensem- ble, which corresponds to all microstates with the same energy, volume and number of particles. The density of representative points in phase space is given by ρ(p,q)= 8 < : const. E <H(p,q)<E +Δ 0 otherwise The question arises, what ensemble describes a system that is not in isolation, but instead is in equilibrium with another, larger system? A system in contact with a heat reservoir such that the temperature, volume and number of particles are constant is described by the canonical ensemble. Is is also possible to have a system where par- ticles can be exchanged with a larger system. Such a system, that is in contact with a heat reservoir and a particle reservoir is described by the grand canonical ensemble, in which temperature, volume and chemical potential are kept constant. 8 The most probable value of an observable is given by the value of O that the most members of the ensemble have. The most probable value and ensemble average are close if the mean square fluctuation hO 2 i−hOi 2 hOi 2 ¿1 is small. Ergodic theorem Under certain conditions, a representative point in phase space will pass arbitrarily close to any other point in the accessible phase space, if one waits a sufficiently long time. Following from the postulate of equal a priori probability, the time average of some physical observable is equal to the ensemble average over phase- space, that is lim T→∞ 1 T Z T 0 O(x,t)dt=hO(x,t)i Statistical mechanics does not specify whether a system is ergodic or not and there is no generic test for ergodicity. For an ergodic system, a single trajectory will uni- formly cover the phase space. An integrable system will not be ergodic in the full phase space due to the additional conserved quantities that act as constraints. Ergodicity says nothing about the time-scale involved in covering the entire phase space. For typical thermodynamic systems, the size of the phase space is immense. A stronger condition than ergodicity is mixing (Tabor, 1989). In the long time limit, the values of a macro- scopic observable in a system that is mixing will equal the ensemble average, without a need to time average over the trajectory as in an ergodic system. 9 1.4 FPU Model and Results Let us now return to the numerical experiments of FPU. Recall that the expectation was that a nonlinear system with many degrees of freedom would exhibit thermal behavior for any non-linearity, no matter how small. The stystem studies was a chain of one- dimensional anharmonically coupled oscillators with quadratic and cubic forces. The Hamiltonians of these two models can be written as a sum of H =H 0 +H 1 , (1.3) where the integrable Hamiltonian H 0 = N−1 X n=1 p 2 n + N−1 X n=1 (x n+1 −x n ) 2 , (1.4) is weakly perturbed by the non-integrable HamiltonianH 1 , which for theα−model is H 1 = α 3 N−1 X n=1 (x n+1 −x n ) 3 (1.5) and for theβ−model is H 1 = β 4 N−1 X n=1 (x n+1 −x n ) 4 . (1.6) The displacement of particlen from equilibrium isx n andα andβ are nonlinear inter- action parameters. The resulting equations of motions for theα-model are ¨ x n =(x n+1 −2x n +x n−1 )+α[(x n+1 −x n ) 2 −(x n −x n−1 ) 2 ] (1.7) 10 and for theβ-model ¨ x n =(x n+1 −2x n +x n−1 )+β[(x n+1 −x n ) 3 −(x n −x n−1 ) 3 ]. (1.8) The problem can be analyzed in terms of normal modes, A k (t), which are the Fourier modes of the displacement,x n (t), A k (t)= r 2 N N X n=1 x n (t)sin µ πkn N ¶ (1.9) In terms of the normal modes,H 0 is a sum of harmonic oscillators andH 1 is a perturbing term that couples the oscillators. H 0 = 1 2 X k ³ ˙ A 2 k +ω 2 k A 2 k ´ ω k =2 ¯ ¯ ¯ ¯ sin µ πk N ¶¯ ¯ ¯ ¯ For the linear system, that is when α = β = 0, there is no interaction between normal modes, so that modes that are initially populated remain populated and modes that are initially unpopulated will remain unpopulated. What happens hen nonlinearity is added? Would the coupling between the modes cause the energy to spread from a single mode to all of the other modes in the system? For the systems studied, thermalization would be marked by equipartition of energy among all of the modes, A k E k = A 0 k E 0 k for all k,k 0 . The expectation was that thermal behavior would be observed. In Fig. 1.1 the time dependence of the normal modes is plotted forN =32 oscillators with cubic forces and β = 8. The initial condition is given by a sine wave and the velocity is zero. The ends are fixed. The model preserves symmetry so that the effective number of particles is 16 and even modes have zero energy. As can be seen from the plot only a few modes are 11 Figure 1.1: Time dependence of normal modes for cubic forces with β = 8. (Fermi et al., 1974) active in the dynamics and there are recurrences. Additionally the period of recurrence was found to decrease with increasing nonlinearity. Later work found a super period, where almost all of energy (99 % ) returns to the initial modes after about 80 000 T 1 . (Tuck and Menzel, 1972). The failure of FPU to thermalize was very puzzling when first discovered and it was some time before explanations came out to explain the observed phenomena. The first explanation has to do with closeness to integrable systems and the second has to do with a stochasticity threshold with respect to the linear system. 12 1.4.1 KAM Theorem A theorem outlined by Komologorov and subsequently proved by Arnold and Moser provided one resolution to the apparent paradox of FPU (Kolmogorov, 1954; Moser, 1962; Arnold, 1963). Consider an integrable Hamiltonian that is weakly perturbed: H = H 0 (I) + ²H 1 (I,θ) where H 0 is integrable with constants of motion I and fre- quencies ω i = ∂H 0 ∂I and H 1 is periodic in the original angle variables, H 1 (I,θ +2π) = H 1 (I,θ). The motion of the unperturbed Hamiltonian corresponds to motion on an n−dimensional torus. It is assumed that the Hamiltonian is analytic on the complex domain and that the unperturbed Hamiltonian is non-degenerate, det ¯ ¯ ¯ ¯ ∂ 2 H 0 ∂I i ∂I j ¯ ¯ ¯ ¯ 6=0. Consider a frequency vector of the unperturbed Hamiltonian, ω ∗ that is incommen- surate (ω ∗ ·k 6= 0 for all integer k i ). The corresponding motion of the unperturbed system is on the torusT 0 (ω ∗ ). One statement of the KAM theorem is Theorem 1.4.1 KAM Theorem IfH 1 is small enough, then for almost allω ∗ , there exists an invariant torus T(ω ∗ ) of the perturbed system such that T(ω ∗ ) is close to T 0 (ω ∗ ). (Arnold and Avez, 1968). Stated informally, KAM showed that if the perturbation is sufficiently small, then almost all of the tori of the unperturbed motion are preserved and the resulting motion is quasi-periodic. 13 1.4.2 Solitons and the Korteweg-de-Vries Equation Soon after the proof of the KAM theorem, Zabusky and Kruskal discovered solitary wave solutions to the Korteweg-de-Vries (KdV) equation, which they termed ”soli- tons”(Zabusky and Kuskal, 1965). Solitons are solitary wave that preserve their shape both under free propagation and after collisions. Zabusky and Kruskal also show that that the continuum limit of theβ− FPU model is close to the KdV equation. In KdV the speed of the solitons depends only on their amplitude. The KdV equation is given by u t +uu x +δ 2 u xxx =0 (1.10) and was first found as a description of the motion of shallow water waves. Starting with a cosine pulse, the negative slope regions of u steepen, and then oscillations develop on the steep front and grow in amplitude, and finally each solitary wave or ”soliton” moves at a constant speed, which is proportional to the amplitude. The trajectories of interacting solitons are shown in Fig. 1.2. The solid (dashed) lines represent the odd- (even-)numbered solitons. From the trajectories, it is clear that when solitons interact, they emerge with the same speed and thus shape. The dotted lines represent interactions, during which the joint amplitude is less than the sum of the individual amplitudes due to the non-linear interaction. At the recurrence timeT R , all of the solitons arrive at one point in space and almost reconstruct the initial state. KdV was later shown to be fully integrable by the method of Inverse Scattering Transform (Gardner et al., 1967). The discovery of the solitons of KdV provides one route for explaining the absence of thermal behavior in FPU. The KAM theorem pro- vides another explanation. The nonlinearities of the FPU studies were insufficient to break the KAM tori and thus the motion remained quasi-periodic. Indeed later work showed that for larger nonlinearities, FPU did exhibit thermal behavior 14 Figure 1.2: Space-time diagram of soliton trajectories. u 0 = cos(πx). δ = 0.022. T R is the recurrence time (Zabusky and Kuskal, 1965). 1.4.3 Chirikov Criterion A method for predicting the criterion for the onset of chaos in the FPU experiments was the theory of overlapping resonances developed by Chirikov in the context of plasma dynamics, and later applied to FPU (Izrailev and Chirikov, 1966). 15 Consider a one-dimensional non-linear oscillator perturbed by an external periodic force. Write the unperturbed system in action-angle variables (I,θ) and the external field as a Fourier series. H =H 0 (I)+² X m,n V mn (I)e i(mθ+nφ) whereH 0 is the unperturbed Hamiltonian with frequenciesω(I)= ∂H 0 ∂I ,φ=Ωt+φ 0 is the phase of the external force andθ =ωt+φ 0 . Resonances occur for the set ofI r such thatω(I r )l =Ωk. The main ingredients are then to: 1. Assume that when the system is near a resonance, the resonance dominates the motion. Thus resonance are studied in isolation. 2. Make a canonical transformation from(I,θ)→ (J,ψ) into the rotating reference frame of the resonance. (ψ measures deviations from resonance). The old and new coordinates are related by J = 1 l (I−I r ), ψ =lθ−kφ. 3. Integrate out the fast phaseφ motion. 4. Expand the Hamiltonian aboutI =I r and keep terms up to second order inH 0 (I) and zeroth order inV nm (I). The resulting Hamiltonian of this process is, H r = J 2 2M +²V l,−k cos(ψ) which is the Hamiltonian of a simple pendulum, with “mass” M ≡ l 2 ³ ∂ 2 H 0 ∂I 2 ´ I=I r . In 16 − 2 0 2 4 6 Angle, ψ − 1.0 − 0.5 0.0 0.5 1.0 Action, J J − 6 − 4 − 2 0 2 4 6 ψ ' − 1.0 − 0.5 0.0 0.5 1.0 J' J' J- J' Figure 1.3: Schematic of Chirikov’s criterion of the onset of chaos for the kicked rotor. J,ψ and J 0 ,ψ 0 are coordinates generated by two different canonical transformations. ΔJ andΔJ 0 are width of the separatricies. J−J 0 represents the distance between the resonances Fig. 1.3 a schematic is plotted forJ,ψ andJ 0 ,ψ 0 coordinates generated by two different canonical transformations corresponding to two different resonances. The separatricies, which separate bounded and unbounded motion, are given by the black contour lines. An initial state inside the separatrix is captured by the resonance and both the action variable and the angle are bounded in phase-space. Outside of the, the angle is unbounded. When the resonances are well-separated in phase-space, that is J 0 − J À ΔJ, the motion near an individual resonance is dominated by that resonance. As the strength of the external drive increases, the width of the separatrix,ΔJ, grows. Eventually the width of the two separatricies becomes comparable to the distance between the resonances. For two neighboring resonances, with resonant values I A and I B , the “distance” between 17 the resonances is given by I B − I A . For separatrix width ΔI, the ratio of these two quantities, K≡ separatrix width distance between resonances = ΔI I B −I A , gives the condition for onset of chaos in the system K&1. Crossing the threshold corresponds to the overlap of the separatrices of the two reso- nances in phase space. so that the action variable is free to travel between resonances and thus explore all of phase space. Complete chaos emerges when all of the neighbor- ing resonances become coupled. Chirikov applied this criterion to the FPU system. For a chain of N oscillators, and an excitation of momentum modek, the conditions for chaos are given by: Low modes k¿N : 3β µ ∂x ∂z ¶ 2 m ∼ 3 k (1.11) High modes (k≈N) N−k¿N : 3β µ ∂x ∂z ¶ 2 m ∼ 3π 2 N 2 µ k N ¶ 2 (1.12) This criterion predicts that for initial excitation of low modes, there is only have stochas- ticity for large perturbations, while for the high mode have stochasticity even for small non-linearities whenN is large. Data from FPU experiments show a threshold for ther- malization that is close to the prediction of Chirikov. 18 1.5 Relationship between Chaos and Thermalization Statistical mechanics dictates the of equilibrium state of a system, but it does not tell us whether a system will thermalize or not. How is it that statistical behavior can arise from dynamics? Consider a system in equilibrium, with some constraint. Lift the constraint and let it evolve. What will the equilibrium state be, if in fact it reaches equilibrium? According to the second law of thermodynamics, when a constraint is lifted, the system moves to a state of greater entropy. To explain this further, consider a simple example a box of volume V with N particles, initially confined to 1/2 of the box. When Boltzmann first derived statistical mechanics for an ideal gas through kinetic equations, there were several objections that were raised (Zaslavsky, 1999): Objection 1: Zermelo’s Paradox (Recurrence) Poincar´ e’s recurrence theorem states that after sufficiently long times, any trajectory will pass arbitrarily close to any phase space point, including initial state. This would contradict the expected increase in entropy. Objection 2: Loschmidt’s Paradox (Reversibility) The equations of motion are invariant under time-reserval. If one simply reverses the velocities of particles, one goes back to the original state, a process that decreases entropy. Microscopic Origins of Macroscopic Irreversibility The second objection to the can be restated, how is it that the microscopic equations of motion are reversible, but the macroscopic behavior is irreversible? Lebowitz argues that this question was satisfac- torily settled years ago by Thomson, Maxwell and Boltzmann (Lebowitz, 1999). The essential ingredients of understanding this, according to Lebowitz are 19 1. the vast difference in scales between microstates and macrostates 2. initial conditions are special 3. significance of probabilities Consider an isolated classical system of N particles. Let X be the microstate that completely specifies the system, X = (r 1 ,p 1 ,r 2 ,p 2 ,...,r N ,p N ) in phase space Γ. Let M represent the macrostate of a system and let Γ M be the the region of the phase spaceΓ that correspond to macrostateM. Note that there are many microstatesX that correspond to the same macrostateM. The number of microstates that correspond to is so large as to make certain macrostates extremely unlikely. Consider the example of the gas of particles initially confined to half of a box. The initial state is indeed a special state. Once the constraint is lifted, the volume in phase space corresponding to all of the particles in one-half of the box is vastly smaller than the volume of phase space that corresponds to a roughly equal distribution of the particles between the two halves of the box. Thus, the probability that it will return to the initial state is effectively zero. Boltzmann made a rough estimate of the Poincar´ e recurrence times and found them to be much larger than the lifetime of the universe and thus irrelevant. This interpretation of the second law of thermodynamics explains why FPU expected to see thermal behavior in their system. However, it fails to account for the thermalization threshold found in FPU. Additionally, thermal behavior has been observed in systems with at few degrees of freedom, with this reasoning does not apply. Zavslasky argues that chaotic dynamics introduce mixing properties in a system (Zaslavsky, 1999) that well resolves these paradoxes. Furthermore, he calculates the distribution of Poincar´ e recurrence times and concludes that they are irrelevant. While there is not a current consensus on the origins of statistical laws, this discussion high- lights some relevant questions, namely, 20 • Is chaos necessary for thermalization? • What is the role of chaos in thermalization in systems with many degrees of free- dom? 1.5.1 Thermalization in Classical Field Models Since the FPU studies on anharmonic oscillators, further studies on thermalization and approach to equilibrium have been carried out in several classical field theories, including recent studies on the classicalφ 4 model (Boyanovsky et al., 2004), nonlinear Klein-Gordon equation (NLKG) (Gerhardt et al., 2002), nonlinear Schr¨ odinger equation (NLSE) (Villain and Lewenstein, 2000; Herbst and Ablowitz, 1989), discrete nonlin- ear Schr¨ odinger equation (DNLS) (Herbst and Ablowitz, 1989; Ablowitz et al., 1993) equivalent to the Bose-Hubbard model, and Integrable Discrete Non-Linear Schr¨ odinger equation (IDNLS)(Herbst and Ablowitz, 1989). No conventional thermalization is expected in the NLSE and IDNLS, which are both integrable. In NLKG, like in FPU, the ability of the system to reach thermal equilibrium in the course of time evolution emerges only when the degree of nonlinearity exceeds a certain critical value (see (Izrailev and Chirikov, 1966; Livi et al., 1985) for the thermal- ization threshold in FPU). On the contrary, theφ 4 model eventually reaches equilibrium regardless of how small the nonlinearity is. There are several other studies on thermalization and chaos in system with a large number of degrees of freedom that are highly relevant to the work presented here. Livi et al. (1985) investigated the equipartition threshold in the FPUβ model in the thermo- dynamic limit. ForN oscillators,64≤N ≤512, the thermodynamic limit is simulated by initially exciting a block of modes, Δn, such that Δn N remains constant. The thresh- old for equipartition of energy is found to be independent of the number of degrees of 21 freedom with respect to the relevant control parameter, the energy density, with a critical value of² c '0.35. They also calculate the Asymptotic Reynolds number, R, given by ¿ O(NL) O(L) À space = β N N−1 X i=1 (φ i+1 −φ i ) 2 →R t→∞ which is a measure of ratio of strength of nonlinear to linear terms. Again, there is universal behavior withR c '0.03, consistent with findings for energy density. There is evidence that the threshold energy is independent of the mode excited, when a narrow range of energies is initially excited. Note that very long equipartition times (t → ∞) are not ruled out and they conclude that the results are relevant for long, but finite times. It is significant to note that the result that the threshold for FPU remains in the thermodynamics limit contradict the predictions of Chirikov’s criterion of overlapping resonances. Another study by the same group focuses on the relationship between chaotic dynamics and statistical mechanics in two nonlinear Hamiltonian systems, the FPU model of nonlinearly coupled oscillators and coupled rotators (Livi et al., 1987). For both systems thermodynamic quantities are computed analytically using ensemble the- ory and compared with dynamical results from numerical simulations. For the FPU model, there is qualitative agreement between ensemble-averages and time-averages, independent of the stochasticity. That is the system is ergodic in both the chaotic and regular region. For the rotator model, there is good agreement between the ensemble- averages and time-averages at low temperatures but not at high temperatures where the system is strongly chaotic. This result is explained in terms of localized chaos. In con- clusion it is possible that a system (a) is not chaotic, but is ergodic for some physically “relevant” quantities and also (b) is chaotic, but some observable are not ergodic. This 22 study highlights the open questions in the relationship between stochasticity and ther- malization, particularly in systems with many degrees of freedom. 1.6 Quantum Degenerate Gases - Ultracold Atoms Advances in the cooling and trapping of alkali atoms into the quantum degenerate regime has led to an explosion of experimental and theoretical studies of ultracold atoms. In recent years Nobel prizes have been awarded for advances in laser cooling techniques (Phillips and Metcalf, 1982; Chu et al., 1985; Aspect et al., 1988) and the subsequent observation of Bose-Einstein condensation (BEC) (Davis et al., 1995; Anderson et al., 1995). The manipulation of atoms by electric and magentic fields offers unprecedented control over parameters in the system and the ability to address fundamental questions in physics. Numerous proposals have been put forth for quantum simulators and appli- cations have arisen in precision measurements. Studies in ultracold atoms have lead to fruitful collaborations across fields, such as condensed matter and quantum informa- tion. One application of these advances is using matter-wave interferometers based on ultracold atomic systems for high precision sensing of accelerations and gravitational fields (Gustavson et al., 2000; Durfee et al., 2006; Weiss et al., 1994; Wicht et al., 2002). Fundamental questions in physics related to out-of-equilibrium dynamics and thermalization in classical and quantum integrable systems have also been studied in one-dimensional ultracold atoms. Quasi-one-dimensional systems have been realized in optical lattices (Paredes et al., 2004; Kinoshita et al., 2006) and on atom chips, where BEC’s have been created and manipulated (Esteve et al., 2006; Schumm et al., 2005; Wang et al., 2005). 23 Figure 1.4: Velocity distribution of rubidium atoms. Left: prior to condensation. Center: just after condensation. Right: after further cooling (Anderson et al., 1995) 1.6.1 Bose-Einstein Condensation Advances in laser cooling and trapping led to the realization of Bose-Einstein conden- sation (BEC) in alkali atoms (Davis et al., 1995; Anderson et al., 1995; Bradley et al., 1995). BEC is a phase of matter, first proposed by Bose (1924) and Einstein (1925), in which there is macroscopic occupation of a single quantum state. Seventy year later BEC was created in Rb 87 gas (Anderson et al., 1995), sodium (Davis et al., 1995) and Li 7 (Bradley et al., 1995). Since the initial experiments, BEC has been observed in twelve species of alkali atoms as well as in Bose molecules (Yukalov, 2009). BEC was created by confining and cooling atoms to microkelvin temperatures with a magneto- optical trap (MOT), followed by evaporative cooling to nanokelvin temperatures. In Fig. 1.4 the velocity distribution of rubidium atoms is show prior to and after conden- sation. The velocity distribution of rubidium atoms is measured by turning off the con- fining trap, allowing the atoms to expand and performing a time-of-flight measurement. The leftmost plot shows the velocity distribution just before condensation. The center 24 plot is the velocity distribution just after condensation, where the sharp peak in veloc- ity distribution clearly indicates the presence of the condensate. In rightmost plot, the system has been cooled further such that most of the atoms are in the condensate. The presence of the condensate was confirmed by the anisotropic velocity distribution due to the magnetic trap in contrast with the isotropic, thermal velocity distribution. BEC’s demonstrate long-range phase coherence, confirmed experimentally by the observation of interference between two independent condensates (Andrews et al., 1997). Historically BEC was defined for a uniform, ideal gas as the macroscopic occupation of a single quantum state in the thermodynamic limit, lim N→∞, V→∞, N/V→const N 0 V >0 whereN is the total number of particles,V is the volume andN 0 is the occupation num- ber of a single quantum state (Yukalov, 2009). The question arises as to to define BEC in a non-uniform system such as in presence of trap. Penrose and Onsager propose that a condensate is present when the largest eigenvalue of the single-particle density matrix is extensive (Penrose and Onsager, 1956). The Penrose-Onsager scheme is more general than idea of off-diagonal long-range order and applicable to both uniform systems and trapped systems. There is no true condensate in finite or 1D systems, which will be the focus of this work, however there can be a quasi-condensate in 1D, when the coherence length is much larger than the de Broglie wavelength (Castin, 2004). 1.6.2 Bosons in Optical Lattices The versatility offered by optical lattices allows one to control parameters such as the interaction strength, lattice spacing and the dimensionality of the system. In particular, 25 one-dimensional systems in cold atoms have been realized in optical lattices by tight confinement in two dimensions. The dynamics of ultracold bosons in optical lattices can be described the Bose- Hubbard model (BHM) . The Bose-Hubbard Hamiltonian is H =−J X hi,ji ˆ b † i ˆ b j + X i ² i ˆ n i + 1 2 U X i ˆ n i (ˆ n i −1). To derive this Hamiltonian, one begins with the Hamiltonian for bosonic atoms in an external potential H = Z d 3 xψ † (x) µ − ~ 2 2m ∇ 2 +V 0 (x)+V T (x) ¶ ψ(x) + 1 2 4πa s ~ 2 m Z d 3 xψ † (x)ψ † (x)ψ(x)ψ(x) whereψ(x) is a boson field operator,V 0 (x) is the potential of the optical lattice,V T (x) is the trapping potential anda s is the s-wave scattering length. For single atoms, energy eigenfunctions are Bloch wave functions. If the energies involved are much less than the excitation energy to the second band then a single band model is justified. Wave functions localized at an individual lattice site, w(x), which are called Wannier wave function are introduced and the energy eigenfunctions are expanded in the Wannier basis ψ(x)= X i b i w(x−x i ). The Bose-Hubbard Hamiltonian follows from this expansion, where the hopping energy between matrix elements is given by J = Z d 3 xw ∗ (x−x i ) · − ~ 2 2m ∇ 2 +V 0 (x) ¸ w(x−x j ), 26 the on-site repulsion is U = 4πa s ~ 2 m Z d 3 x|w(x)| 4 , and the energy offset due to the lattice is ² i = Z d 3 xV T (x)|w(x−x i )| 2 ≈V T (x i ). In the last step the trapping potential is assumed to be approximately constant over the spatial variation of a single Wannier function. One expects a zero-temperature quantum phase transition from the superfluid (SF) state to the Mott insulator (MI) state and as the depth of the lattice is increased for integer fillings. The SF state supports long-range phase coherence while in the MI state, the atoms are localized and there is no phase coherence. This transition was observed in ultracold atoms by Greiner et al. (2002). 1.6.3 Atom Chips Quasi-one-dimensional systems have also been realized on atom chips. The miniatur- ization and integration of matter-wave optics has led to the development of atom chips (Folman et al., 2002; Fortagh and Zimmermann, 2007). It is now possible to confine, manipulate and measure atoms on a single device using electric, magnetic and optical fields. Bose-Einstein condenstation has been created in magnetic microtraps (Ott et al., 2001; Hansel et al., 2001). The traps are highly elongated and the one-dimensional regime is realized when the transvere confining potential,~ω ⊥ is much greater than the relevant energy scales of the system, the thermal energy, k B T and chemical potential, μ. Esteve et al. (2006) realized both the ideal Bose Gas as well as the quasicondenstate in a quasi-one-dimensional trap on an atom chip. Theoretical work has investigated the transition from the 1D Bose gas to the quasicondensate (Bouchoule et al., 2007) as well 27 as the growth of the quasicondensate (Proukakis et al., 2006). Other experiments in one-dimensional traps include the demonstration of the first phase-preserving matter- wave beam-splitter on an atom chip (Schumm et al., 2005) and of an atom Michelson interferometer on an atom chip (Wang et al., 2005). 1.6.4 Classical Field Model of Bose Gas In this thesis we numerically study the Bose-Hubbard model, presented earlier, within the classical-field approximation. The classical field approximation is equivalent to the first-order mean-field approximation. In this section we outline the validity of the clas- sical field approach for studying the dynamics of interacting Bose gases. The dynamics of a BEC can be well-described by the Gross-Piteavskii equation (GPE) (Gross, 1961; Pitaevskii, 1961) i~ ∂Ψ(r,t) ∂t = µ − ~ 2 2m ∇+V(r)+g|Ψ(r,t)| 2 ¶ Ψ(r,t), where the coupling constant is given by g = 4π~ 2 a/m (a =scattering length, m = mass). The Gross-Pitaevskii equation is a mean-field approximation and is equivalent to the continuous nonlinear Schr¨ odinger equation which is integrable. The GPE has been used extensively to describe the dynamics of the condensate in three-dimensional systems. Several studies have looked at the applicability of the mean-field or classical field description beyond the dynamics of the condensate. Kagan and Svistunov studied the evolution of an interacting Bose gas from a strongly non-equilibrium state towards con- denstation. They demonstrated that the classical-field description accurately describes a weakly interacting Bose gas in the absence of a condensate provided that the occu- pation numbers of the initially occupied state are much greater than unity (Kagan and 28 Svistunov, 1997). Given this condition, the time evolution of a state can be accurately described by the diagonal elements of the statistical matrix in the coherent state repre- sentation. Castin studies the classical field model for one-dimensional weakly interacting Bose gases (Castin, 2004). The classical field model is generated by replacing the quantum mechanical operator ˆ ψ(z) with a complex field ψ(z). For the interacting Bose gas, the state of the classical field is governed by a single parameter, χ= ~ 2 ρ 2 mk B T ρg k B T where ρ is the mean density, T is the temperature, m is the mass, and g is the interac- tion parameter. Castin calculates the correlation functions g 1 (z) = D ˆ ψ † (z) ˆ ψ(0) E and g 2 (z) = D ˆ ψ † (z) ˆ ψ † (0) ˆ ψ(0) ˆ ψ(z) E . The contrast, C ≡ g 2 (0)/g 2 1 (0) drops off quickly as χ increases and then slowly approaches unity for χ À 1, indicating that density fluctuations are suppressed for largeχ. The conditions for the validity of the classical field model forχÀ1 are summarized as : 1. large occupation numbers,k B T Àμ 2. gas is degenerate,ρλÀ1 3. weakly interacting regime,ρξÀ1 where ξ = (~/mμ) 1/2 is the healing length. Note that the condition that χ À 1 auto- matically satisfies conditions (2) and (3). In summary, the classical field model is a good approximation for weakly-interacting particles of a degenerate gas for large occupation number, in which case fluctuations are suppressed. 29 Mishmash and Carr study the correspondence between the mean-field and the fully quantum BHM in the dynamics of atoms in 1D optical lattices (Mishmash and Carr, 2008). The mean-field BHM is equivalent to the discrete nonlinear Schr¨ odinger equation (DNLS). They numerically investigate the analogs of dark soliton of DLNS in BHM and use the time-evolving block decimation algorithm (TEBD) developed by Vidal (Vidal, 2004) to carry out the full quantum calculations. 1.6.5 Chaos and Integrability in Quantum Systems Access to one-dimensional systems of ultracold atoms in optical lattices has led to real- ization of some known integrable models and observed effects of integrability in the dynamics of these models. We focus on the effects of integrability in bosonic systems in optical lattices. The Lieb-Liniger model is a completely integrable quantum description of one-dimensional bosons with two-bodyδ-interactions (Lieb, 1963; Lieb and Liniger, 1963). The Hamiltonian is ˆ H =− N X i=1 ∂ 2 ∂x 2 i +2c X hi,ji δ(x i −x j ), (1.13) where the interaction, governed by the parameterc is repulsive. The Lieb-Liniger model has been solved via Bethe Ansatz. In the limit of infinitely strongδ− repulsions,c→∞, the hard-core bosons, also known as a Tonks-Girardeau gas, map to non-interacting fermions (Girardeau, 1960). While these models were proposed a half-century ago, the Tonks-Girardeau gas was only recently realized experimentally in Rb 87 atoms that were strongly confined in two directions in an optical lattice to create one-dimensional tubes (Paredes et al., 2004). By applying a shallow lattice in the longitudinal direction, the effective mass and thus interaction strength were increased in order to reach the Tonks- Giradeau regime. 30 Figure 1.5: f(p ex ), momentum distribution for (a) γ 0 =4, τ=34 ms. t blue = 15τ, t red = 30τ (b) γ 0 =1, τ=13 ms.t blue = 15τ, t red = 40τ (c) γ 0 =0.62, τ=13 ms. t blue = 15τ, t red = 40τ. Green is the initial momentum distribution averaged over the first period. (Kinoshita et al., 2006) Later experiments observed the effects of integrability on thermalization in a one- dimensional Bose gas. Kinoshita et al. (2006) demonstrate the first experimental evi- dence for the lack of thermalization in a many-body system with a large number of degrees of freedom for bosons in optical lattices. A gas of interacting bosons was pre- pared out-of-equilibrium by applying a laser pulse to a one-dimensional Bose-Einstein condensate in an optical lattice. For both strongly- and weakly-interacting bosons, the expanded momentum distribution retains the initial double peak structure. Even with the background harmonic potential, the system is integrable in the limit of infinite-strength repulsion. It was expected that a system with finite interactions, which is believed to be non-integrable in the presence of a harmonic trap, would reach thermal equilibrium. However, the absence of thermalization occurred even for finite interactions. Figure 1.5 shows the expanded momentum distribution for three different coupling strengths. From 31 the three peak structure in Fig. 1.5(a)-(b) it is clear that the gas has not thermalized for the Tonks-Girardeau limit (γ 0 = 4) and the intermediate regime (γ 0 = 1) even after thousands of collisions have occurred between atoms. 1.6.6 Constrained equilibrium One question that arises from these experiments on hard-core bosons is: do integrable systems, which don’t relax to the usual thermodynamic equilibrium distribution attain some other steady state? Numerical studies on one-dimensional hard-core bosons in a lattice addressed the relaxation dynamics of a fully-integrable quantum system (Rigol et al., 2007). The method to derive the steady-state distribution is to maximize the entropy subject to the constraints of the system, which include all of the conserved quantities (Jaynes, 1957a,b). In this approach the many-body density matrix is given by ˆ ρ=Z −1 exp " − X m λ m ˆ I m # where{ ˆ I m },{λ m } are the Lagrange multipliers which are determined by the initial con- ditions. This distribution is called the generalized Gibbs ensemble or fully-constrained thermodynamic ensemble. One-dimensional hard-core bosons on a lattice can be mapped to free-fermions via a Jordan-Wigner transformation. The conserved quantities are moments of the fermionic momentum distribution. Rigol et al. solve analytically for the density matrix with the constraints from the fermionic momentum distribution. The results of numerical simu- lations confirmed that when the system is prepared in the ground state of a small box and then allowed to expand in a larger box, it reaches a steady state, which is in agreement 32 with the analytic results for the fully constrained system, rather than the grand canon- ical thermodynamic distribution. Additionally for an initial state with two momentum peaks, the two peaks structure remains after many oscillations, in agreement with the experiment performed by Kinoshita et al. (2006). 1.7 Outline of Thesis In this work we present a comprehensive study of chaos and thermalization in the 1D Bose-Hubbard model within the classical-field approximation. We study the threshold for chaos and its relation to thermalization. Two quantitative measures of thermaliz- ability are compared: the Finite-time Maximal Lyapunov exponents (FTMLE) and the normalized spectral entropy (NSE). The FTMLE, averaged over phase space, converges to the maximal Lyapunov Exponent, the standard measure of chaos. A positive MLE indicates that points that are initially close in phase-space diverge exponentially, rather than linearly. The spectral entropy measures the distance between the time-averaged momentum distribution of the numerical results and the momentum distribution pre- dicted by thermodynamics, within the independent mode approximation. We investigate the dependence of the averaged FTMLE and normalized spectral entropy on a dimen- sionless nonlinearity parameter and the energy-per-particle, both of which are finite in the thermodynamic limit. The BHM is found to have a threshold for chaos which depends on the nonlinearity and the energy-per-particle. We study the size scaling of the Lyapunov exponent and normalized spectral entropy. Furthermore we study resonances in the Bose-Hubbard model to find the Chrikov criterion for chaos. The criterion predicted by the Chrikov criterion is different from the one inferred from numerical calculations, signifying the failure of the standard Chirkov’s approach. 33 There are at least three near-by integrable models: the Ablowitz-Ladik lattice, the continuous nonlinear Schr¨ odinger equation and the noninteracting model. We outline the method of Inverse Scattering Transform and generate all of the integrals of motion of the closely related, fully integrable model of Ablowitz-Ladik. Furthermore, we dis- cuss the possible role of these conserved quantities in relaxation in the BHM. We con- jecture that the presense of quasi-conserved quantities may alter the scaling of the chaos criterion. 34 Chapter 2 Thermalization and Chaos in the 1D Bose-Hubbard Model 2.1 Introduction One of the fundamental assertions of statistical mechanics is that the time average of a physical observable is equivalent to the average over phase-space, with microcanonical measure. A system for which this is true is said to be ergodic and one can calculate dynamical properties of the system from static phase-space averages. While this is believed to be true, because of the success of statistical mechanics in accurately pre- dicting experimental results, many open questions remain. Is ergodicity sufficient to ensure the accuracy of statstical mechanical predictions for times that are relevant for observations? 2.2 BHM: Hamiltonian and Equations of Motion We study the dynamics of an interacting one-dimensional Bose gas on a lattice (1D Bose-Hubbard model (BHM)) (Jaksch et al., 1998) with periodic boundary conditions in the classical field approximation. The Hamiltonian of the system of interest can be studied in many different forms through canonical transformations. Each equivalent rep- resentation has a different Hamiltonian, canonical coordinates and equations of motion. A well-choosen canonical transformation can pose the problem in a way which is more 35 intuitive. Here we present three different representations: The real-space representa- tion, the momentum-space representation in terms of classical fields and the action-angle representation. 2.2.1 Real-Space Hamiltonian In real space, the Hamiltonian is H =−J X s ˜ ψ ∗ s ³ ˜ ψ s+1 + ˜ ψ s−1 ´ + μ 0 N s 2 X s | ˜ ψ s | 4 . (2.1) The equations of motion are given by ∂ ∂t ˜ ψ s =− i ~ ∂H ∂ ˜ ψ ∗ s =−i h −J ³ ˜ ψ s+1 + ˜ ψ s−1 −2 ˜ ψ s ´ +μ 0 N s | ˜ ψ s | 2 ˜ ψ s i (2.2) ∂ ∂t ˜ ψ ∗ s = i ~ ∂H ∂ ˜ ψ s = i h −J ³ ˜ ψ ∗ s+1 + ˜ ψ ∗ s−1 −2 ˜ ψ ∗ s ´ +μ 0 N s | ˜ ψ s | 2 ˜ ψ ∗ s i , (2.3) and the canonical pairs are Q s = ψ s ,P s = i~ψ ∗ s . These equations of motion are equivalent to the discrete nonlinear Schr¨ odinger equation (DNLS). The time-evolution of the fields is carried out in real-space in all of the calculations. 2.2.2 Momentum-Space Representation Another set of canonical coordinates, the momentum-space fields, ψ n = ψ(k n ), are related to the real-space field ˜ ψ s = ˜ ψ(x s ), by ψ n = 1 √ N s Ns X s=1 ˜ ψ s e −ikn·xs ˜ ψ s = 1 √ N s N s X n=1 ψ n e ik n ·x s . (2.4) 36 wherek n =2πn/L andx s =sa. In the momentum-space representation, the Hamiltonian is H = X m ³ ~ω m |ψ m | 2 − μ 0 2 |ψ m | 4 ´ + μ 0 2 X l,i,j ψ ∗ l ψ ∗ i ψ j ψ l+i−j (2.5) where the sum carries the restrictions: j 6= l,i. Several canonical transformations have been performed in order to write the Hamiltonian in this form. The canonical pairs are Q n = ψ n ,P n = i~ψ ∗ n . and the equations of motion are given by ∂ ∂t ψ n =− i ~ ∂H ∂ψ ∗ n =−i ³ ω n − μ 0 ~ |ψ n | 2 ´ ψ n −i μ 0 ~ X i,j ψ ∗ i ψ j ψ n+i−j (2.6) ∂ ∂t ψ n =− i ~ ∂H ∂ψ ∗ n =−i ³ ω n − μ 0 ~ |ψ n | 2 ´ ψ n −i μ 0 ~ X i,j ψ ∗ i ψ j ψ n+i−j (2.7) with the restrictions j 6= n,i on the sums and the indices span the range n, i, j = 0,±1,±2, ...,± Ns−1 2 (N s is supposed to be odd). The bare frequency of each momentum mode is given by ~ω n = ~ 2 ma 2 µ 1−cos µ 2πn N s ¶¶ =2J µ 1−cos µ 2πn N s ¶¶ . (2.8) The coupling constant is μ 0 = UN a /N s . Here J and U are the nearest-neighbor site- hopping and on-site repulsion constants of the standard Bose-Hubbard model respec- tively, andN a is the number of atoms. Time propagation is performed in real space, while the output and analysis of the numerical calculations focus on the momentum fields. 37 Table 2.1: List of Parameters and Variables a = lattice spacing J = B-H nearest-neighbor kinetic energy κ = nonlinearity parameter L = length of lattice μ 0 = coupling constant N s = number of lattice sites = number of momentum modes N a = number of atoms τ tal = Talbot time U = B-H on-site repulsion energy ~˜ ω 1 = ground state energy of non-interacting model with quadratic dispersion ξ = size-dependent nonlinearity parameter 2.2.3 Action-Angle Representation Equivalently, the momentum-space Hamiltonian can be written in terms of action-angle variables, by performing serveral canonical transformation on the momentum-space rep- resentation. It is this action-angle representation that will be the starting point of the resonant approximations and studies of individual resonances. The Hamiltonian is H = X m ³ ω m I m − μ 0 2~ 2 I 2 m ´ + μ 0 2~ 2 X m,l,i,j (I m I l I i I j ) 1/2 e i(θm+θ l −θ i −θ j ) (2.9) where the sum carries the restrictions: m + l = i + j; m 6= i,j; l 6= i,j. The momentum wavefunction canonical variables, {(ψ k ,i~ψ ∗ k )} are related to the action- angle canonical variables {(I k ,θ k )} by ψ k = q I k ~ e −iθ k , ψ ∗ k = q I k ~ e iθ k . In this 38 form, the Hamiltonian can be seen as the sum of an integrable term and a perturbation, H({I k ,θ k })=H 0 ({I k })+V({I k ,θ k }), where the integrable Hamiltonian is H 0 = X m ³ ω m I m − μ 0 2~ 2 I 2 m ´ = X m ³ ~ω m |ψ m | 2 − μ 0 2 |ψ m | 4 ´ . (2.10) Throughout the wavefunctionψ n is normalized to unity: P n |ψ n | 2 =1. Dimensionless nonlinearity parameter We define the dimensionless nonlinearity parameter, κ≡ μ 0 J ≡ U(N at /N s ) 2 J(N a /N s ) (2.11) whose physical meaning is the ratio between the typical interaction energy per site U(N a /N s ) 2 and the kinetic energy per site JN a /N s . Note that this parameter governs both the strength of the nonlinearity and the strength of the perturbation from the inte- grable Hamiltonian (2.10). In Tables 2.1 and 2.2, variables are listed and the relationship between relevant parameters are summarized. 2.2.4 Validity of the Classical-Field Theory Based on the studies of the validity of the classical-field theory for Bose gases (Castin, 2004; Kagan and Svistunov, 1997) discussed earlier, the classical-field approximation will apply for the lattice site occupations satisfying N a N s À max(κ, 1) max ·µ N s Δn ¶ , 1 ¸ , 39 Table 2.2: Relationship between Parameters ~˜ ω 1 = ~ 2 2m µ 2π L ¶ 2 =J µ 2π N s ¶ 2 τ tal = 2π ˜ ω 1 μ 0 = gN a L J = ~ 2 2ma 2 U = g a =μ 0 N s N a U/J = 2mag ~ 2 = 2mL 2 ~ 2 N s N a μ 0 κ= μ 0 J ≡ U(N at /N s ) 2 J(N a /N s ) ξ≡ μ 0 ~˜ ω 1 =κ µ N s 2π ¶ 2 ~ω n = ~ 2 ma 2 · 1−cos µ 2πn N s ¶¸ =2J · 1−cos µ 2πn N s ¶¸ where Δn is the typical width of the momentum distribution. We note that the Mott regime,N a = integer×N s ,Δn=N s ,U/J ≥2.2N a /N s (Hamer and Kogut, 1979), lies well outside of the above criteria. 2.2.5 Nearby Integrable Models There are several known intergrable models that are limiting cases of the BHM. These include 40 1. Continuous Nonlinear Schr¨ odinger Equation: In the contiuum limit, the 1D BHM becomes H = Z L 0 dzΨ ∗ ½ − ~ 2 2m Δ ¾ Ψ+ Z L 0 dz 1 2 g|Ψ| 4 . 2. Linear Model: In the non-interacting limit,κ→0, the 1D BHM becomes a sum of harmonic oscillators H = X m ~ω m |ψ m | 2 . 3. Independent Mode: If the interating term of vanishes, the 1D BHM becomes a sum of decoupled nonlinear oscillators, H = X m ~ω m |ψ m | 2 − μ 0 2 |ψ m | 4 with nonlinear frequencies given byΩ m =~ω m − μ 0 2 |ψ m | 2 . 4. Ablowitz-Ladik Lattice: An alternate discretization of the NLS yields the Ablowitz-Ladik lattice, with Hamiltonian, H =− X n ¡ q n q ∗ n+1 +q ∗ n q n+1 ¢ − 4 σ X n ln ³ 1− σ 2 |q n | 2 ´ , which will be discussed further in later chapters. 2.3 Time Dynamics First we study the time dynamics of the 1D BHM on a lattice withN s =21 modes. The system is prepared in a state that is narrowly distributed in momentum space and evolves according to the classical equations of motion. Initially the lowest three momentum 41 modes are occupied, the minimum number of modes required by selection rules for non-trivial processes leading to population of initially unoccupied modes. In Fig. 2.1, -0.5 -0.25 0 0.25 0.5 (a) κ=0.009 -0.5 -0.25 0 0.25 0.5 κ=0.09 -0.5 -0.25 0 0.25 0.5 Re ψ x=0 κ=0.36 -0.5 -0.25 0 0.25 0.5 0 2 4 6 8 10 κ=0.9 time [τ talbot ] 0 0.08 0.16 (b) Ω 0 Ω 1+,1- κ=0.009 0 0.08 Ω 0 Ω 1+,1- κ=0.09 0 0.08 Power Spectrum of ψ x=0 (t) κ=0.36 0 0.08 -15 -10 -5 0 κ=0.9 Frequency, ω [ω ~ 1 ] Figure 2.1: Time evolution of the wavefunction at the center of the box for a typical run for κ = 0.009,0.09,0.36,0.9 with identical initial conditions. (a) Time dynamics of the real part of the wavefunction, Re ψ x=0 (t). (b) Frequencies of the real-space wavefunction,|FT[ψ x=0 (t)]| 2 . For κ = 0.009,0.09 the descendants of the unperturbed frequencies generated by the first, “integrable”term of the Hamiltonian (2.5) are labeled. the time dynamics and power spectrum of the wavefunction at the center of the box, ψ x=0 (t), are plotted for various interaction strengths for a typical initial state. As seen in Fig. 2.1(a), the time evolutions of the zero momentum mode is quasi-periodic for weak interactions, with a few easily identifiable frequencies entering the dynamics, which is confirmed by the power spectrum in Fig. 2.1(b). As the nonlinearity increases, more frequencies determine the dynamics and for sufficiently large nonlinearity the motion loses its quasi-periodic character and appears to be chaotic. The clear distinction between quasi-periodic and seemingly chaotic behavior of the time dynamics leads to the following questions: 1. Is the motion really chaotic? 2. If it is, where is the chaos threshold as one increases the nonlinearityκ? 3. When chaotic, does the system reach thermal equilibrium? 42 In order to answers these questions, it is necessary to define appropriate measures of chaos and thermalization. 2.4 Chaos: Calculating Lyapunov Exponents The standard signature of the chaotic nature of a region in phase-space is that the sep- aration between trajectories that are initially close grows exponentially with time, for typical trajectories, as captured by a positive maximal Lyapunov exponent (MLE). In regular regions the separation grows linearly (Chirikov, 1979), resulting in zero MLE. As we increaseκ in our system, we expect the phase space to change from being dom- inated by regular regions for small κ to being dominated by chaotic regions for large κ. In the present section, we use the MLEs to quantify this transition to chaos, which, as we will see in the subsequent section, coincides with a relatively broad change from unthermalizability to complete thermalizability. Consider two trajectoriesx(t) ande x(t) with initial pointsx 0 ande x 0 , respectively. The separationδx(t) =e x(t)−x(t) initially satisfies a linear differential equation, and the duration of this linear regime grows without bound as the initial separatione x 0 −x 0 goes to zero. The finite-time maximal Lyapunov exponent (FTMLE) corresponding to the phase-space pointx 0 (Eckhardt and Yao, 1993; V oglis and Contopoulos, 1994) is given by λ t fin (x 0 )= lim e x 0 →x 0 1 t fin ln ke x(t fin )−x(t fin )k ke x 0 −x 0 k . (2.12) The limitt fin →∞ gives the MLE,λ ∞ (x 0 ), but the FTMLE are themselves of intrinsic interest (Eckhardt and Yao, 1993; V oglis and Contopoulos, 1994; Contopoulos et al., 1978; Contopoulos and V oglis, 1997). We chose a convenient quantum mechanical 43 metric,ke x−xk 2 = P n | e ψ n −ψ n | 2 . This metric becomes Euclidian under the canonical transformation Q n =(2~) 1/2 Re(ψ n ) P n =(2~) 1/2 Im(ψ n ) ke x−xk 2 = 1 2~ X n £ (Q 0 n −Q n ) 2 +(P 0 n −P n ) 2 ¤ where the sum runs fromn=−(N s −1)/2 ton=(N s −1)/2. 2.5 Thermalization: Calculating Spectral Entropy In order to measure thermalization, it is necessary to make thermodynamic predictions and compare the dynamics from the propagation of the equations of motion, with the thermodynamic state. In this section a method for calculating the thermodynamic state, within the Hartree-Fock approximation, is laid out. A full account is given in Appendix A. Additionally the spectral entropy is defined, which is a quantitiative measure of the difference between the time-averaged dynamical state and the expected thermodynamic state. 2.5.1 Conserved Quantities A treatment of the thermodynamic state must take into account all conserved quantities. The known conserved quantities of the 1D BHM are the energy and norm. Conservation of norm, P N n=1 |ψ n | 2 , is associated with U(1) symmetry in the real/imaginary plane represented by the tranformation{ψ n ,ψ ∗ n } → {ψ n e iθ ,ψ ∗ n e −iθ }. Unlike the continuous NLSE, the total momentum is not conserved. 44 2.5.2 Hartree-Fock Within Hartree-Fock the form of the density distribution function is taken to be Gaussian and the thermal expectation value of the Grand Potental, hFi=hHi−ThSi−μhN a i, (2.13) is minimized, where N a is the norm. The density distribution function with two-body interactions has the form, σ HF = 1 Z exp Ã X n,n 0 −α n,n 0ψ n ψ ∗ n 0 ! = 1 Z exp Ã − X n α n I n ! . (2.14) We use the independent mode approximation so that in the second step, the off-diagonal elements are taken to be zero. Theα n coefficients are unknown and are determined by the condition of minimizing the grand potential. The density distribution function is normalized, so that the integration ofσ HF over all of phase space is 1, by Z = 1 (2π~) N Z d N θ Z d N Ie − P n α n I n = N Y i=1 1 ~α i . (2.15) The expectation value of each term in the Grand Potential is calculated, using the Hartree-Fock density distribution function σ HF . The expectation value of a generic observable is given by hOi= 1 (2π~) N Z d N θ Z d N IO({I n ,θ n })σ HF (2.16) = 1 (2π) N N Y i=1 α i Z d N θ Z d N IOe − P n α n I n (2.17) The expectation values of the relevant observables are listed below. 45 Norm hN a i= * 1 ~ X n I n + = N X n=1 1 ~α n (2.18) Hamiltonian - Kinetic Term hH 0 i= * X n I n ω n + = N X n=1 ω n α n (2.19) Hamiltonian - Interaction Term hH I i= * μ 0 2~ 2 X m,p,q,r (I m I p I q I r ) 1/2 δ m+p,q+r e −i(θm+θp−θq−θr) + = μ 0 ~ 2 X m,p 1 α m 1 α p (2.20) Entropy S =− 1 (2π~) N Z d N θ Z d N σ HF logσ HF =N− X j log(~α j ) (2.21) 2.5.3 Minimization of the Grand Potential The thermal expectation value of the Grand Potential within Hartree-Fock is given by hFi= X m ω m α m + μ 0 ~ 2 X m,p 1 α m 1 α p −T Ã N− X m log(~α m ) ! −μ X m 1 ~α m (2.22) Taking the variation with respect toα n , and setting it equal to zero gives δhFi δα n =− ω n α 2 n −2 μ 0 ~ 2 1 α 2 n X m 1 α m + T α n + μ ~α 2 n =0. (2.23) 46 Using P m α −1 m =~N a and solving forα n , α n = 1 ~T [~ω n +2μ 0 N a −μ] (2.24) The thermal expectation values of the occupation of momentum moden become hI n i= 1 α n = ~T ~ω n +2μ 0 N a −μ . (2.25) In general, the coefficientsμ andT are unknown and are determined by imposing con- straints on the norm and energy, which come from the dynamical code. The constraints are N a =hN a i= 1 ~ X n hI n i E T =hHi= X n ω n hI n i+ μ 0 ~ 2 X m,n hI m ihI n i= X n ω n hI n i+μ 0 N 2 a (2.26) Beginning with the expression forhI n i, we can solve forT in terms ofμ,N a and energy, T =ω n hI n i+2 μ 0 ~ N a hI n i− μ ~ hI n i. (2.27) Summing overn, T = 1 N £ E k +2μ 0 N 2 a −μN a ¤ (2.28) where E k ≡ P n ω n hI n i. This expression for T can be substituted back into the con- traints to reduce the system to two equations with two unknowns. Using the expression for temperature and normalization condition, a single constraint remains to be solved, 1 N X n [E k +2μ 0 N 2 a −μN a ] [~ω n +2μ 0 N a −μ] −N a =0. (2.29) 47 The Hartree-Fock approximation is known to overestimate the interaction energy in the regime of strong interactions. For sufficiently largeμ 0 , the Hartree-Fock interaction energy,μ 0 N 2 a becomes greater than the total energy resulting in negative kinetic energy, where the kinetic energy isE k =E T −μ 0 N 2 a . For this reason, we determine the temper- atureT and the chemical potentialμ using the time-averaged numerical kinetic energy (along with the norm) instead of the total energy. The quantityE k in the thermal distri- bution is fixed to the time-averaged kinetic energy of the final state from the dynamical code. We fix the norm to its numerical value, and subsequently solve for the norm to find all parameters. In this way, the total energy is never used in the constraints. Additionally the solutions must satisfy the physical constraint that I n ≥ 0 for all n, which leads to bounds on μ. For T > 0, the condition such that the denominator is greater than zero for all n, is ω n −2μ 0 N a −μ > 0, which leads to an upper bound for μ, μ < 2μ 0 N a . There is a critical kinetic energy that corresponds to infinite tem- perature, which leads to equal population of all the modes,hI n i = N −1 s . The critical kinetic energy, which separates the positive and negative temperature regime can be cal- culated asE k−cr = P n 1 Ns ω n = 1 Ns P n N 2 (1−cos ¡ 2πn N ¢ = N 2 . ForE k > E k−cr the temerature is negative, and the lower bound on μ is ω n +2μ 0 N a −μ < 0 for all n or μ > 2N 2 s +2μ 0 N a . Close to the critical kinetic energy, both the temperature and the chemical potential diverge. By expanding the norm in powers ofω n /(μ−2μ 0 ), an esti- mate for the chemical potential whenE k =E k−cr ±ε isμ≈2μ 0 +E 2 k−cr ±E 2 k−cr /(2ε). The temperature and the chemical potential were computed individually for each initial condition used. In Fig. 2.2, the initial and time-averaged momentum distributions of a representa- tive state are plotted for κ = 0.09,0.36 and 0.9, along with the thermal Hartree-Fock predictions,h|ψ n | 2 i=(T/N a )/(~ω n +2μ 0 N a −μ). 48 0 0.2 0.4 0.6 0.8 1 κ=0.09 ε T =0.15 J initial final thermal 0 0.2 0.4 0.6 0.8 Momentum distribution, |ψ(n)| 2 κ=0.45 ε T =0.41 J initial final thermal 0 0.2 0.4 0.6 0.8 -8 -6 -4 -2 0 2 4 6 8 Momentum index, n κ=1.8 ε T =1.4 J initial final thermal Figure 2.2: Initial, final and Hartree-Fock thermal momentum distributions for κ = 0.09,0.45,1.8, starting from the same initial state. N s = 21. The initial state is a representative state and the final state is time-averaged. ² T is the total energy per particle. 2.5.4 Spectral Entropy For coupled anharmonic oscillators, as in the FPU study, energy equipartition among the normal momentum modes signified thermalization. In the BHM, the additional conservation of the norm modifies the quantity that is equipartitioned. To determine the best measure for the equipartition we use the variational Hartree-Fock Hamiltonian (Castin, 2001; Ohberg and Stenholm, 1997), H HF = P n ~ω HF n |ψ n | 2 , where the set of Hartree-Fock energies {~ω HF n } was regarded as the variational field. This procedure gives~ω HF n =~ω n +2μ 0 N a −μ, whereμ is the chemical potential. The new quantity to be equipartitioned is the distribution of the Hartree-Fock energy, q n (t)= |ψ n (t)| 2 ~ω HF n P n 0 |ψ n 0(t)| 2 ~ω HF n 0 (2.30) 49 A quantitative measure of the distance from thermodynamic equilibrium is the spec- tral entropy S(t)=− X n q n (t)lnq n (t). (2.31) In thermal equillibrium,q n is equipartioned, maximizing the spectral entropy at a value S max =log(N s ). A more convenient quantity to study is the normalized spectral entropy (Livi et al., 1985), η(t)= S max −S(t) S max −S(0) , (2.32) which is unity att=0 and vanishes as the system approaches thermal equilibrium. 2.6 Results: N=21 Sites, Three-Mode Initial Conditions Initially, we study the FTMLE for a class of initial conditions where only thek =0,±1 modes are occupied. In this subspace we sample uniformly from the intersection of the microcanonical shells in energy and norm; the energy is chosen to be the infinite- temperature energy of the subsystem, and the norm is 1. For each value ofκ, we sample 100 points, which we set as the initial pointsx 0 . To each initial point we add a small random vector, as little as machine precision allows, to obtain the correspondinge x 0 ’s. Each pair we propagate for a timet fin , short enough to ensure linearity of the evolution of δx(t) but long enough to be able to clearly distinguish chaotic trajectories from regular ones on a plot oflnδx(t) versust: the former are straight lines of positive slope, while the latter are logarithm-like (Contopoulos et al., 1978). We also verify that the average of the FTMLE’s over the ensemble of initial conditions does not depend ont fin as long as both criteria above are satisfied. In Fig. 2.3 the averaged FTMLEs are plotted as a function of the interaction strength. There is a distinct regime with zero Lyapunov exponent for small κ . 0.5 and a strongly chaotic regime for κ & 1 where all initial 50 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 0 0.1 0.2 0.3 0.4 spectral entropy, η Lyapunov exponent, λ [J] nonlinearity parameter, κ η λ Figure 2.3: Ensemble-averaged finite-time maximal Lyapunov exponent,λ, and normal- ized spectral entropy,η, as a function of the nonlinearity,κ. N s =21. conditions have positive exponent. Next we consider the relation between chaos and thermalization in the system. In Fig. 2.3 the spectral entropy of the final time-averaged state, also averaged over 100 initial states (drawn from the same ensemble that was used for the Lyapunov expo- nent calculation) is plotted for each value of κ. For large nonlinearities, κ & 1, the normalized spectral entropy goes to zero, indicating remarkable agreement between the final state and the thermal predictions. Note that this corresponds to chaos threshold observed previously. Forκ. .5 the normalized spectral entropy is above.5 signifying that during the time evolution the state of the system remains close to the initial state. In the Fig. 2.4, the normalized spectral entropy is plotted versus the FTMLE for each of the 100 individual runs forκ = 0.36,0.54,0.72,0.9. As seen in the plot, an individual initial state with larger FTMLE tends to have lower spectral entropy, i.e. to relax to a state which is closer to the thermal one. Beginning at κ ≈ 0.5, where the averaged FTMLE is substantially non-zero, some of the initial states thermalize completely. 51 0 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 0.08 0.1 0.12 spectral entropy, η Lyapunov exponent, λ [J] κ=0.18 κ=0.36 κ=0.54 κ=0.72 κ=0.9 Figure 2.4: Normalized spectral entropy of the final time-averaged state versus finite- time maximal Lyapunov exponent for each of the 100 initial condition used to compute the averaged value forκ=0.36,0.54,0.72,0.9. N s =21. Method of Calculating Spectral Entropy For an individual realization, we calculate μ 0 and T given the norm and the final kinetic energy. The spectral entropy is cal- culated for (1) the initial momentum distribution and (2) the final mometum distribu- tion,h|ψ n (t)| 2 i, which is the time-averaged population of momentum moden from the dynamical code. 2.6.1 Fluctuations In order to confirm that the system is thermal whenη→0, we investigate the scaling of fluctuations of the kinetic energy for largeκ. For a system withN particles and volumeV in thermal equillibrium, the fluctuations of extensive observables scale with the size of the system as V −1/2 , in the thermody- namic limit,N →∞,V →∞,N/V = const. Consider an extensive quantity,O. The temporal fluctuations ofO are given by σ 2 O = O 2 ® −hOi 2 . 52 wherehOi = 1 T R t+ΔT/2 t−ΔT/2 dt O(t), and the relative fluctuations are σ O /hOi. Now con- sider two idential systems,A andB, with corresponding extensive observablesO A and O B , such thatO =O A +O B . Then σ 2 O = O 2 ® −hOi 2 . = O 2 A +O 2 B +2O A O B ® −hO A +O B i 2 = O 2 A ® + O 2 B ® +2hO A O B i−hO A i 2 −hO B i 2 −2hO A ihO B i =σ 2 O A +σ 2 O B +2hO A O B i−2hO A ihO B i For a system in thermal equilibrium, the two parts of the system are decorrelated, so thathO A O B i=hO A ihO B i and the last two terms cancel. For a general rescaling of the systemN 0 = αN and the extensive observablehOi 0 = αhOi, the fluctuations scale as σ 0 O = √ ασ O and the relative fluctuations scale asσ 0 O /hOi=α −1/2 σ O /hOi. In contrast, consider the case where two identical systems in identical initial states are concatenated in a regime where the behavior is regular. In this case, there will be strong correlations between the two parts of the system. In the extreme case of O A (t) = O B (t), then σ O =2σ O A =2σ O B and the relative fluctuationsσ 0 O /hOi will be constant, independent of the size of the system. For the system under consideration, the thermodynamic limit is taken by scaling the number of atoms and length asN 0 at = αN at ,N 0 s = αN s , while the interaction parameter μ 0 and the lattice spacing, a, remain constant (U, J = constant as well). In order to simulate the thermodynamic limit, the initial conditions are generated by concatenating α duplicates of the real-space wavefunction of the reference lattice. This is equivalent to generating an initial state with momentum modesψ 0 k , given byψ 0 0 =ψ 0 ,ψ 0 ±α =ψ ±1 , from the initial state in the N s0 lattice, ψ k . Generating the initial conditions in this way preserves the average energy per particle, in units of J = ~ 2 /2ma 2 . A small 53 Table 2.3: Thermodynamic Limit and Scaling ~˜ ω 0 1 = ~ 2 2m µ 2π L 0 ¶ 2 = ~˜ ω 1 α 2 τ 0 talbot = 2π ˜ ω 1 =α 2 τ talbot μ 0 0 = gN at L =μ 0 J 0 = ~ 2 2ma 2 =J U 0 = g a =U U 0 /J 0 =U/J ξ 0 = μ 0 ~˜ ω 1 = N at N s 4π 2 U J =α 2 ξ ~ω 0 n =2J 0 ½ 1−cos µ 2πn N 0 s ¶¾ =2J ½ 1−cos µ 2πn αN s ¶¾ perturbation is added to the initial wavefunction to break the symmetry associated with the translational invariance. For the thermodynamic limit, we want to take the case where N at → ∞,L → ∞, with N at /L = const, with the lattice spacing, a, and interaction strength, g, remaining constant. Consider the case where L 0 = αL, N 0 at = αN at . The scaling of relevant parameters is given in Table 2.3. We study the standard deviation of the fluctuations for systems with N s = 21 sites andα=2,3,4. To compare fluctuations for different lengths, we calculate: ¯ σ(N s ,N s0 )≡ σ E k (N s )/E k (N s ) σ E k (N s0 )/E k (N s0 ) (2.33) 54 0.1 1 1 σ(N s , N s0 ) N s /N s0 κ = 2.69 (N s /N s0 ) -1/2 Figure 2.5: Fluctuations in kinetic energy. Normalized standard deviation/mean for N s =21, 42, 63, 84 lattice sites and N 0 = 21. Data points are for five sample runs with equivalent initial conditions for different lattice sizes andκ=2.69. where σ E (N s ) and E k (N s ) are the standard deviation and time average of the kinetic energy for a chain with N s lattice sites. The reference lattice size is N s0 = 21 and the calculation begins after the system has already thermalized. In Fig. 2.5 we plot ¯ σ(N s ,N s0 ) as a function of α = N s /N s0 for various lattice sizes on a log-log scale for κ = 2.69. As can be seen clearly from the plot, the fluctuations scale as N −1/2 s , indicating that the fluctuations are indeed thermal. 2.7 Chaos Threshold for Different Lattice Sizes Let us start from the notion that the parameterκ introduced in (2.11) is the only dimen- sionless combination of the parameters of the problem that remains finite in the ther- modynamic limit, N s → ∞, N a /N s = const,J = const,U = const. Curiously, the chaos threshold for N s = 21 is at κ ≈ .5, i.e. κ ∼ 1. Another observation comes from a related work (Villain and Lewenstein, 2000) on chaos threshold in NLSE with 55 hard-wall boundary conditions. The authors find that the boundary between regular and chaotic motions of momentum mode,n, is given by(μ 0 |ψ n | 2 )/(~ω 1 n)∼ 1, where~ω 1 is the lowest excitation energy, e.g. the energy of the first excited mode in the case of the Hamiltonian (2.5). Assuming that the shape of the momentum distribution|ψ n | 2 as a function of n/N s should be fixed in the thermodynamic limit, the left-hand side of the above relationship also remains finite. These observations lead to a conjecture that the chaos criterion involves only intensive parameters and observables, i.e. those that are finite in the thermodynamic limit. Our test for the above conjecture is based on the 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Nonlinearity parameter, κ 0.0 0.1 0.2 0.3 0.4 Lyapunov exponent, λ [J] N s =11 N s =21 N s =41 Figure 2.6: Averaged Finite Time Lyapunov exponent, λ[J], for three different system sizes,N s =11,21,41. For eachκ, the same energy-per-particle was used for each lattice size. The error bars represent one standard deviation. fact that for a chaotic motion the majority of the trajectories cover the whole available phase space, and as a result the LE becomes (for a given set of parameters) a function of the energy only. This implies that for the same energy-per-particle and the same non- linearity parameter,κ, Lyapunov exponents for different lattice sizes should be similar. In Fig. 2.6 the time-averaged finite-time maximal Lypunov exponent is plotted for three different lattices, N s = 11, 21, and 41. For each κ, the same energy-per-particle (in units of J) is used for all three lattices. From the plot it is indeed evident that the LE is 56 universal with respect to the size of the lattice and that the values of the LE forN s =11 already give a very good estimate of both the value of the LE and the threshold. 2.8 Two Parametric Theory of the Chaos Threshold The universality observed above suggests the most relevant pair of variables for map- ping the chaos threshold, namelyκ and the total energy-per-particle, ² T /J. In order to test these parameters, we independently vary the nonlinearity, κ and the energy-per- particle, ² T . For each κ and ² T , we generate ten initial states with microcanonical weight in the reduced phase space of three (n = 0,±1 ) or five (n = 0,±1,±2 ) momentum modes. It is necessary to generate initial states with five momentum modes n=0,±1,±2 because there is an upper limit on the energy of an initial state with only three modes occupied. The finite-time maximal Lyapunov exponent and normalized spectral entropy are calculated for each of the ten realizations and averaged over this ensemble. For the rest of this work we will use the term Lyapunov exponent (LE) to denote the ensemble-averaged finite time maximal Lyapunov exponent and normalized spectral entropy (NSE) to denote the ensemble-averaged normalized spectral entropy, unless otherwise specified. The total data represents over 3200 runs. In Fig. 2.7(a) contour lines of the LE forN s =11 are plotted versus the nonlinearity parameter and total energy-per-particle. One can observe an initial plateau in the LE for λ. λ c = 0.02, given by the solid line. The threshold depends on both the nonlinearity and the total energy-per-particle. Based on parameter regime investigated, it appears that the threshold persists for smallκ no matter how much energy is present. For small energies it is unclear whether the threshold will persist or vanish forκÀ1. After crossing the critical line the LE increases with uniform slope. The critical line resembles a hyperbola with the point of the closest approach to the origin at (κ, ² T )∼ 57 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Nonlinearity parameter, κ 0.0 0.2 0.4 0.6 0.8 1.0 Total energy per particle, ǫ T [J] λ =0.01 0.00 0.04 0.10 0.16 0.22 0.28 0.34 0.40 0.46 Lyapunov exponent, λ [J] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Nonlinearity parameter, κ 0.0 0.2 0.4 0.6 0.8 1.0 Total energy per particle, ǫ T [J] λ =0.01 η =0.36 η =0.68 0.00 0.12 0.24 0.36 0.48 0.60 0.72 0.84 0.96 Normalized Spectral Entropy, η Figure 2.7: (a) Contour lines of the averaged FTMLE versus the nonlinearity, κ, and energy-per-particle, ² T = (H−H 0 )/N a , whereH is the Hamiltonian (2.5), and H 0 = −2J +(1/2)μ 0 is the ground state value of H. The solid contour line corresponds to λ c =0.01. The diagonal solid line represents the set of energies and nonlinearities used in Fig. 2.6. (b) Contour lines of the averaged normalized spectral entropy versus the nonlinearity,κ and energy-per-particle. Solid contour lines correspond toη = 0.68 and η = 0.36. For reference, the threshold line from the FTMLE in (a) is plotted (dashed line). N s =11. (0.5, 0.2J), so that the hopping parameterJ appears to be a relevant energy scale. This is probably not an accident: for² T ÀJ the dispersion lawω n begins to deviate from the (quadratic) dispersion law of the continuous NLSE with periodic boundary conditions, which is integrable. 58 The normalized spectral entropy was calculated for the same set of data runs and is plotted in Fig. 2.7(b). There are two solid contour lines, atη = 0.68 andη = 0.36. The second contour line at η = 0.36 follows closely the dotted line, which is the threshold from the Lyapunov exponent. It is apparent that the two plots have the same general features, and that there is a strong correspondence between the presense of chaos and thermalization in the BHM. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Nonlinearity parameter, κ 0.0 0.2 0.4 0.6 0.8 1.0 Total energy per particle, ǫ T [J] Figure 2.8: Data points used for interpolation for contour plots in Fig. 2.7 superimposed on data for Lyapunov exponent. A few features deserve discussion. In the region of η & 0.68,λ = 0, which is enclosed by the first contour line in the NSE and the x- and y-axes, the system relaxes to a state that is closer to the thermal state even though the Lyapunov exponent is zero. The time dependence of the spectral entropy reveals that this relaxation takes place very quickly, after which the spectral entropy remains flat for many fundamental cycles, indicating that no further relaxation will take place (See Fig. 2.11). This raises a few questions: To what state does the system relax in this region? It is possible to describe it by a constrained ensemble, where 59 the constrained quantities are the conserved quantities of near-by intergrable systems? We will return to these questions later. For λ À λ c , the region of strong chaos, the majority of initial states relax to the thermal state and all final states are close to the thermal. It is likely in this region, that full relaxation is not seen due to slow relaxation times and the spectral entropy would vanish for longer propagation times. We consider the two limits, κ → 0,² T ∼ J and ² T /J → 0,κ ≈ 3. In the limit of small κ, the chaos threshold and NSE contour line η = 0.36 overlap and for κ → 0, ² T & 0.6J converge to a value that is independent of the total energy-per-particle. For the parameter region explored there is no indication that the threshold will vanish, every for very large energies. This suggests a dependence on the ratio of the nonlinear to linear terms, similar to the critical Reynolds number found in FPU (Livi et al., 1985). In the opposite limit of ² T /J → 0,κ & 1.5 the behavior is quite different. While the Lypunov exponent is zero, there is significant relaxation in the momentum distribution. It is important to note that in the limit that² T → 0 the initial state approaches the state where only the n = 0 mode is populated, which is also the thermal state and thus the normalized spectral entropy is not well-defined in that limit. For this reason, data is plotted for ² T > 0.02J, the lowest energies simulated. However even for ² T ∼ 0.02J, relaxation is visible in the momentum distributions. In Fig. 2.9 the standard deviation of the normalized spectral entropy is plotted along with the contour lines at η = 0.36 and η = 0.68 from Fig. 2.7(b). The contour line at η = 0.36 follows closely the chaos threshold from Fig. 2.7(a). Far above the threshold, where theη → 0, the standard deviation is also small indicating that most of the states themalize, as expected. Belowη =0.68, (in the region bounded by the axes), relaxation is minimal and the standard deviation is small, indicating that most initial states will not thermalize. In the vacinity of the threshold (λ = λ c , which is close to η = 0.36), the 60 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Nonlinearity parameter, κ 0.0 0.2 0.4 0.6 0.8 1.0 Total energy per particle, ǫ T [J] η =0.36 η =0.68 0.00 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 Standard Deviation of η Figure 2.9: Contour lines of the standared deviation of the NSE versus the nonlinearity, κ and energy-per-particle. Solid contour line corresponds to at η = 0.36 andη = 0.68 from Fig. 2.7(b). N s =11. standard deviation is larger. We conjecture that this is because there is a large spread in the amount of relaxation expected for different initial states with the same parameters and/or that some states have not fully relaxed due to insufficient propagation times as a result of multiple relaxation time scales. In Fig. 2.10 the normalized spectral entropy is plotted for three different lattice lengths, N s = 11,21,41 for the same energy-per-particle at each κ, using the same energy values as in Fig. 2.6. While the main features are similar for the different lat- tice sizes, the size scaling of the NSE is not as universal as the scaling of Lyapuonv exponent. For small κ’s with the same energy-per-particle, more relaxation is seen in larger lattices. This suggests that the number of modes in involved in the dynamics may plays a role in relaxation. In addition the standard deviation of the NSE is larger than for the Lyapunov exponents. In contrast to the size scaling of the Lyapunov exponents, the variance of the NSE increases with smaller chains. We conjecture that the large vari- ance can be attributed to multiple relaxation scales. For example, forN s =11,κ=1.5, 61 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Nonlinearity parameter, κ 0.0 0.2 0.4 0.6 0.8 Normalized Spectral Entropy, η N s =11 N s =21 N s =41 Figure 2.10: Averaged normalized spectral entropy, η, for three different system sizes, N s =11,21,41. For eachκ, the same energy-per-particle was used for each lattice size. individual states reveal that while most relax fully to the thermal state there is a single state that remains very far from thermal, which is the cause of the large variance. The large variance of the spectral entropy is the reason that there are more features in the contour plot of the normalized spectral entropy compared with the contour plot of the Lyapunov exponent. Repeating the simulations for longer times would likely smooth some of the features of the NSE contour plot, and decrease the variance in regions where it is currently large. 2.8.1 Thermalization Times and Slow Relaxation The time dependence of the normalized spectral entropy is plotted in Fig. 2.11 forκ = 0.09 and² T =0.081J, which is in the non-chaotic region. The initial, thermal and final momentum distributions are plotted in the inset. The time-dependent spectral entropy is calculated from a running average over the momentum distribution and plotted in units of τ tal , the talbot time, which is the period associated with the lowest frequency of the non-interacting system with quadratic dispersion. After an initial relaxation during the 62 first talbot time, the spectral entropy saturates and remains flat for close to1000τ tal . The momentum distributions confirm that there is some relaxation, but that the state remains far from the thermal state. 0 0.2 0.4 0.6 0.8 1 0.01 0.1 1 10 100 1000 Norm. Spectral Entropy, η(t) t [τ tal ] κ=0.09 ε T =0.081 0 0.2 0.4 0.6 0.8 1 -4 -3 -2 -1 0 1 2 3 4 |ψ k | 2 k Figure 2.11: Sample time dependence of normalized spectral entropy,η. κ=0.09,² T = 0.081J,N s =21. Inset: Initial (dashed red), final (solid blue) and thermal (dotted black) momentum distributions. In Fig. 2.12 the time dependence of the normalized spectral entropy is plotted for two different initial states that both have κ = 0.54 and total energy ² T = 0.19J. In Fig. 2.12(a) the normalized spectral entropy vanishes indicating that the state relaxes to the thermal state, which is also seen in the final momentum distribution. The normalized spectral entropy drops in several stages suggesting that there are multiple relaxations time scales. In Fig. 2.12(b), the state does not fully relax during the observed propa- gation time. After the initial relaxation, which is very similar to the previous case, the normalized spectral entropy slowly relaxes further but does not vanish in the observed time. Both states have a positive finite-time maximal Lyapunov exponent and thus are in the chaotic regime. These observations brings up several questions. Given states in the chaotic sea with the same total energy and strength of nonlinearity, why do some fully relax while other do not? Will these states fully thermalize for longer propagation times? What are the relevant time scales? What governs the slow relaxation times? 63 0 0.2 0.4 0.6 0.8 1 0.001 0.01 0.1 1 10 100 Norm. Spectral Entropy, η(t) t [τ tal ] κ=0.54 ε T =0.19 0 0.2 0.4 0.6 0.8 1 -4 -3 -2 -1 0 1 2 3 4 |ψ k | 2 k 0 0.2 0.4 0.6 0.8 1 0.001 0.01 0.1 1 10 100 Norm. Spectral Entropy, η(t) t [τ tal ] κ=0.54 ε T =0.19 0 0.2 0.4 0.6 0.8 1 -4 -3 -2 -1 0 1 2 3 4 |ψ k | 2 k Figure 2.12: Sample time dependence of normalized spectral entropy,η, for two differ- ent initial states. κ = 0.54,² T = 0.19J for both. N s = 21. Inset: Initial (dashed red), final (solid blue) and thermal (dotted black) momentum distributions. Comparing the momentum distributions for both plots, the initial momentum distri- bution is almost symmetric in Fig. 2.12(a) so that the total quasi-momentum of the initial state is close to zero. Total quasi-momentum is not a conserved quantity of the BHM, although it is a conserved quantity of the noninteracting model, the continous model (in which case is becomes the true momentum) and the Ablowitz-Ladik discretization of the NLS. The total quasi-momentum is zero in the thermal state. For smallκ’s, there is very little redistribution among the momentum modes and thus the total quasi-mementum is well conserved. We call the total quasi-momentum a “quasi-conserved quantity”because it is not actually conserved in the BHM, but it conserved in the nearby integrable models 64 and thus is expected to be conserved in the BHM when is “close”to one of the integrable limits. Proposed future work includes the investigation of the role of the conserved quantities of the nearby integrable models in the dynamics of the BHM. While these plots are sample runs, the pattern just described is observed in other indi- vidual runs for different values ofκ and² T in the chaotic region. Relaxation occurs on multiple time scales and the propagation times used in the simulations are long enough for the fast relaxation, but are not always long enough for the slow relaxation. For a given set of parameters, (κ,² T ) there are different slow relaxation times for different initial states. Insufficient propagation times are one possible reason for large variation of the individual NSE’s observed in Fig. 2.9. In the strongly chaotic regime, it is expected that the the normalized spectral entropy will converge to zero for longer propagation times. However it is also possible that forλ≈λ c , (η≈0.36) some initial states will not fully relax, even for very long times. Furthermore, for0.68& η& 0.36 it is likely that the variance will remain large. It is clear from Fig. 2.11 that some states do not relax, even for very long times. In summary, we have observed a threshold for chaos in the BHM, which depends on two parameters, the strength of the nonlinearity, κ and the total energy-per-particle, ² T /J. Far above the threshold, the state relaxes to the one predicted by statistical mechanics. Below the chaos threshold, we observe relaxation to a non-thermal steady- state. For small nonlinearities, κ’s the chaos threshold and absense of thermalization persist even for large energy-per-particle, ² T ∼ J. For regions just above the thresh- old, there are multiple relaxation times, with different intitial states relaxing on different time scales. These observations bring up several questions: What is the origin of the chaos threshold? What governs the long relaxation times? Is the nonthermal steady- state affected by the conserved quantities of the nearby integrable systems? 65 Chapter 3 Resonance Model and Failure of Chirikov’s Criterion In this chapter we study a Hamiltonian coupling a single set of modes for initial states with few modes occupied. First we use a perturbation theory expansion that is valid for small nonlinearities. Second we study nonlinear resonances of the one-dimensional mean-field Bose-Hubbard model to predict the chaos criterion by Chirikov’s method of overlapping resonances. 3.1 Perturbation Theory Study of BHM To study a perturbation theory expansion we introduce the a prefactor ² in the driving term, which will be set to unity in the end. The full Hamiltonian in the momentum-space wavefunction representation is: H = X m ³ ~ω m |ψ m | 2 − μ 0 2 |ψ m | 4 ´ +² μ 0 2 X m 0 ,l 0 ,i 0 ,j 0 ψ ∗ m 0ψ ∗ l 0ψ i 0ψ j 0 (3.1) where the sum carries the restrictions: m 0 +l 0 =i 0 +j 0 ;m 0 6=i 0 ,j 0 ;l 0 6=i 0 ,j 0 . We define the unperturbed Hamiltonian as H 0 = X m ³ ~ω m |ψ m | 2 − μ 0 2 |ψ m | 4 ´ , (3.2) 66 which consists of decoupled nonlinear oscillators. The frequencies of the unperturbed Hamiltonian are Ω n ≡ω n −(μ 0 /~)|ψ n | 2 = ˜ ω 1 n 2 −(μ 0 /~)|ψ n | 2 . The equations of motion of the full Hamiltonian are given by: ∂ ∂t ψ n =− i ~ ∂H ∂ψ ∗ n =−i ³ ω n − μ 0 ~ |ψ n | 2 ´ ψ n −i² μ 0 ~ X l,i,j;i6=l,n ψ ∗ l ψ i ψ n+l−i (3.3) =−iΩ n ψ n −i² μ 0 ~ X l,i;i6=n,l ψ ∗ l ψ i ψ n+l−i (3.4) Now we make the following assumption: (1) All modes except n are treated as inde- pendent oscillators with equations of motion given by the unperturbed Hamiltonian, ˙ ψ l =−iΩψ l which has solutions ofψ l (t) = ¯ ψ l e −iΩ l t and (2) the nonlinear frequency of mode n associated with the unperturbed Hamiltonian is fixed, Ω n = const. With these restrictions, the equations of motion of moden become ∂ ∂t ψ n =−iΩ n ψ n −i² μ 0 ~ X l,j;j6=n,l ¯ ψ l ¯ ψ j ¯ ψ n+l−j e −i(Ω j +Ω n+l−j −Ω l )t . (3.5) 3.1.1 Dynamics of an Initially Unpopulated Mode First we study the maximum value of initially unpopulated mode in the perturba- tion theory expansion. We consider the initial state where a block of modes from [−N max ,N max ] are equally populated and study the time dynamics of mode Q = N max +1. In order to do this, we seek a solution of d dt ψ n +iΩ n ψ n =f(t) (3.6) 67 We want to find an integrating factorh(t), satisfying ˙ h(t)=iΩ n h(t) such that d dt [ψ n (t)h(t)]=h(t) d dt ψ n +iΩ n h(t)ψ n =f(t)h(t). (3.7) The solution ish(t)=exp(iΩ n t). Letf(t)= P α g(α)exp(iΔ α t) where X α = X l,j;j6=n,l g(α)=− μ 0 ~ ¯ ψ l ¯ ψ j ¯ ψ n+l−j Δ α =Ω l −Ω j −Ω n+l−j . (3.8) Next we integrate over time Z τ 0 dt d dt ¡ ψ n (t)e iΩnt ¢ = Z τ 0 dt X α g(α)e i(Ωn+Δα)t ψ n (τ)e iΩ n τ −ψ n (0)= X α g(α) i(Ω n +Δ α ) ¡ e i(Ω n +Δ α )τ −1 ¢ ψ n (t)= X α g(α) i(Ω n +Δ α ) ¡ e iΔ α t −e −iΩ n t ¢ (3.9) The time-averaged value ofψ n (t) is given by |ψ n | 2 = lim T→∞ 1 T Z T 0 dt|ψ n (t)| 2 = lim T→∞ 1 T Z T 0 dt X α,β g(α)g(β) ¡ e i(Δα−Δ β )t −e i(Δα−Ωn)t −e −i(Δ β −Ωn)t −1 ¢ (Ω n +Δ α )(Ω n +Δ β ) = X α,β g(α)g(β) (Ω n +Δ α )(Ω n +Δ β ) [δ(Δ α −Δ β )−δ(Δ α −Ω n )−δ(Δ β −Ω n )−1] = X α,β g(α)g(β) (Ω n +Δ α )(Ω n +Δ β ) [δ(Δ α −Δ β )−1], (3.10) 68 We introduce the dimesionless variables, ξ = μ 0 ~˜ ω 1 , Ã ξ =κ µ N s 2π ¶ 2 ! which is directly proportional to the nonlinearity parameter, κ, but scales with the size of the system in the thermodynamic limit. We also introduce a new time scale,τ = ˜ ω 1 t. Occupation of mode N max +1 for Quadratic Dispersion We are interested in the population of mode Q ≡ N max +1 when initially modes (−N max ,N max ) are equally occupied with population ¯ ψ. When only low momentum modes are excited, the true cosine dispersion laws can be approximated by a quadratic dispersion relation,ω n = ˜ ω 1 n 2 and the unperturbed fre- quencies can be written asΩ l =ω l (l 2 −ξ| ¯ ψ| 2 ). Dropping termsO(ξ) in the denominator, the equation of motion forψ Q can be written as ψ Q (t)=−²ξ˜ ω 1 e −iΩ Q t ¯ ψ 3 X |j|,|l|<Q;j6=l e i(Ω Q +Ω l −Ω j −Ω Q+l−j )t −1 Ω Q +Ω l −Ω j −Ω Q+l−j ψ Q (τ)=−²ξ ¯ ψ 3 X |j|,|l|<Q j6=l e i(l 2 −j 2 −(Q+l−j) 2 )τ −e −iQ 2 τ Q 2 +l 2 −j 2 −(Q+l−j) 2 ψ Q (τ)=−²ξ ¯ ψ 3 X |j|,|i|<Q ¯ ψ j ¯ ψ Q+l−j e i(2(i−Q)(j−Q)−Q 2 )τ −e −iQ 2 τ 2(i−Q)(j−Q) (3.11) where we make the substitutionl = j +i−Q in the last expression. Consider the case where modes 0,±1 are initially occupied with equal occupation, ¯ ψ and Q = 2. The 69 nonzero terms in the sum correspond to modes(2,0↔1,1), (2,−1↔1,0), (2,−1↔ 0,1). The time evolution becomes: ψ 2 (t)=−e −iΩ 2 t μ 0 ¯ ψ 3 ~ · e i(Ω 2 +Ω 0 −2Ω 1 )t −1 Ω 2 +Ω 0 −2Ω 1 +2 e i(Ω 2 +Ω −1 −Ω 0 −Ω 1 )t −1 Ω 2 +Ω −1 −Ω 0 −Ω 1 ¸ ψ 2 (τ)=ξ ¯ ψ 3 e −i4τ µ e i2τ −1 2 +2 e i4τ −1 4 ¶ =ξ ¯ ψ 3 µ 1+e −i2τ −2e −i4τ 2 ¶ (3.12) So that forN max =1,the time evolution of|ψ 2 (τ)| 2 is: |ψ 2 (τ)| 2 =ξ 2 ¯ ψ 6 µ 3−cos(2τ)−2cos(4τ) 2 ¶ (3.13) Occupation of mode N max + 1 for Cosine Dispersion In the previous case, we assumed a quadratic dispersion, which corresponds to the free space. The true dis- persion of the BHM, which is a lattice model, is cosine, ω n =−˜ ω 1 N 2 s 2π 2 µ 1−N 2 s cos µ 2πn N s ¶ . ¶ (3.14) and the time evolution of mode population|ψ 2 (τ)| 2 becomes ψ 2 (t)=−ξ ¯ ψ 3 e iN 2 s cos(4π/N s )t ˜ ω 1 N 2 s Ã e itN 2 s (cos(4π/N s )+1−2cos(2π/N s )) −1 −cos(4π/N s )−1+2cos(2π/N s ) +2 e itN 2 s (cos(4π/N s )−1) −1 −cos(4π/N s )+1 ! DefiningΔ 1 ≡ N 2 s (cos(4π/N s )+1−2cos(2π/N s ))/˜ ω 1 ,Δ 2 ≡ N 2 s (cos(4π/N s )− 1)/˜ ω 1 , the |ψ 2 (τ)| 2 =ξ 2 ¯ ψ 6 µ e iτΔ 1 −1 Δ 1 +2 e iτΔ 2 −1 Δ 2 ¶µ e −iτΔ 1 −1 Δ 1 +2 e −iτΔ 2 −1 Δ 2 ¶ (3.15) 70 |ψ 2 (τ)| 2 =ξ 2 ¯ ψ 6 Ã e iτΔ 1 −1 Δ 1 +2 e iτΔ 2 −1 Δ 2 !Ã e −iτΔ 1 −1 Δ 1 +2 e −iτΔ 2 −1 Δ 2 ! =ξ 2 ¯ ψ 6 Ã 2(1−cos(τΔ 1 )) Δ 2 1 +8 1−cos(τΔ 2 ) Δ 2 2 + 4(cos(τ(Δ 1 −Δ 2 ))−cos(τΔ 1 )−cos(τΔ 2 )+1) Δ 1 Δ 2 ! (3.16) This expression for the time evolution of the modes is expected to be accurate for small values of ξ. Comparison of these predictions with numerics shows very good agreement for smallξ. 3.1.2 Nonlinear Frequencies in the Real-Space Dynamics For small values of ξ, the perturbation theory expansion accurately predicts the time- dynamics of the real-space wavefunctions. Next we consider the same model as the nonlinear coupling increases. In Fig. 3.1 we plot the modulus-squared of the Fourier transform of the real-space wavefunction at the center of the box,ψ x=0 (t), is plotted for N s = 21 and ξ = 0.1,1,2,3 for an initial state with momentum modes ψ n=0 ,ψ n=±1 occupied. Using a non-interacting model (no coupling between modes, but nonlinear terms associated with each mode), the predicted value of each of the resonances isΩ n = ω n − μ 0 ~ 2 I n . Taking I n as the time-averaged value of I n (t). The driving term is of the forme i(Ω m +Ω n −Ω l ) , so additional frequencies are expected to be linear combinations of three other frequencies. 71 0 0.05 0.1 0.15 0.2 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 |FT of ψ x=0 (t)| 2 ω/ω ~ 1 Ω 0 Ω 1+ Ω 1− Ω 0 Ω 1+ Ω 1− Ω α Ω β Ω χ Ω δ Ω ε Ω φ Ω γ Ω 0 Ω 1+ Ω 1− Ω α Ω β Ω χ Ω δ Ω ε Ω φ Ω γ Ω 0 Ω 1+ Ω 1− Ω χ Ω δ Ω ε Ω φ Ω γ ξ=0.1 ξ=1 ξ=2 ξ=3 Figure 3.1: Frequencies of the real-space wavefunction at the center of the box, |FT[ψ x=0 (t)] 2 | forξ =0.1,1,2,3. Ω α =2Ω 0 −Ω 1- Ω β =2Ω 0 −Ω 1+ Ω χ =Ω 0 +Ω 1+ −Ω 1- Ω δ =Ω 0 +Ω 1- −Ω 1+ Ω ² =2Ω 1+ −Ω 0 Ω φ =Ω 0 −Ω 1+ −Ω 1- Ω γ =2Ω 1- −Ω 0 3.2 Chirikov’s Criterion Next, we study the Bose-Hubbard Model in the Born-Oppenheimer approximation (BOA). In the BOA, the populations of all but one mode are fixed. The fixed modes rotate in phase-space with constant frequency. The motion of the free mode is coupled to the the fixed modes and motion is governed by the equations of motion. 72 In the next section, we deduce the Chirikov chaos criterion by studying single res- onances within the Born Oppenheimer approximation. Within the resonant approxima- tion it is assumed that near a resonance, the resonance dominates the motion, so that driving terms can be studies independently. Within the BOA, we assume that the action variable of the single mode under study is small compared to the other modes, p,q,r, and that the other modes can be well-described by the integrable Hamiltonian, H 0 that describes independent modes, H 0 = X m ³ ω m − μ 0 2~ 2 I m ´ I m . Solving the equations of motion for the action variable, ˙ I m =− ∂H 0 ∂θm givesI m = const≡ ¯ I m , while the equation of motion for the angle variable, ˙ θ m = ∂H 0 ∂I m gives θ m = ³ ω m − μ 0 ~ 2 ¯ I m ´ t≡Ω m ( ¯ I m )t. (3.17) Thus the action variable of modes p,q,r is fixed and the angle rotates with constant frequency. Summary of assumptions: Throughout this section, we assume that 1. Single resonance approximation: we isolate individual driving terms of the Hamil- tonian and study the resonances of these models. 2. Born-Oppenheimer approximation: The population of all modes, except the one under study are fixed. The fixed modes have frequenciesΩ m = ¡ ω m − μ 0 ~ 2 ¯ I m ¢ . 3. Quadratic Dispersion: For low-energy modes the dispersion ω n = 2˜ ω 1 (N s /2π) 2 [1−cos(2πn/N s )] can be approximated by a quadratic dispersion, ω n = ˜ ω 1 n 2 . 73 4. Equal population of fixed modes: The populations of the fixed modes are taken to be equal in order to reduce the number of parameters. 3.2.1 Classification of Resonances We classify the single resonances into three categories. We use the notation (p,q) → (n,r), where the first two modes are the “feeding” modes and the second two modes are the “filling modes”. Mode n is the mode of interest, whose dynamics are governed by the resonant Hamiltonian, while p,q and r are three different modes with fixed action variables. Additionally, the momentum indices must satisfy p+q = r +n. The res- onances are classified according to the order of the exponent of mode n in the driving term. The three classes of Hamiltonians are: 1. First-order Resonance(p,q)→(r,n) The two feeding modesp,q fill two different modes,r,n. The first-order resonant Hamiltonian is H n =ω n I n − μ 0 2~ 2 I 2 n + 4μ 0 ~ 2 (I n I p I q I r ) 1/2 cos(θ n +θ r −θ p −θ q ). (3.18) 2. 1-R Resonance(p,p)→(r,n) A single feeding modep fills two different modes,r,n. The 1-R resonant Hamil- tonian is H n =ω n I n − μ 0 2~ 2 I 2 n + 2μ 0 ~ 2 ¡ I n I r I 2 p ¢ 1/2 cos(θ n +θ r −2θ p ). (3.19) The 1-R resonance differs from the first-order resonance only in the prefactor in front of the driving term. 74 Table 3.1: Resonance Parameters Order Modes Bare Detuning Mean Occupation of Fixed Modes First p,q→r,n p+q =r+n Δ 01 ≡p 2 +q 2 −r 2 −n 2 ¯ I 1 ≡( ¯ I p ¯ I q ¯ I r ) 3/2 1R p,p→r,n 2p=r+n Δ 01r ≡2p 2 −r 2 −n 2 ¯ I 1r ≡( ¯ I p 2 ¯ I r ) 3/2 Second p,q→n,n p+q =2n Δ 02 ≡p 2 +q 2 −2n 2 ¯ I 2 ≡( ¯ I p ¯ I q ) 1/2 Order Drive Frequency Nonlinear Detuning First ν 1 ≡Ω p +Ω q −Ω r =ω p +ω q −ω r − 2μ 0 ~ ( ¯ I p + ¯ I p − ¯ I r ) Δ 1 ≡(ν 1 −ω n )/˜ ω 1 =Δ 01 −ξ ¯ I 1 1R ν 1r ≡2Ω p −Ω r =2ω p −ω r − 2μ 0 ~ (2 ¯ I p − ¯ I r ) Δ 1r ≡(ν 1r −ω n )/˜ ω 1 =Δ 01r −ξ ¯ I 1r Second ν 2 ≡Ω p +Ω q =ω p +ω q − 2μ 0 ~ ( ¯ I p + ¯ I q ) Δ 2 ≡(ν 2 −2ω n )/(2˜ ω 1 ) =Δ 02 /2−ξ ¯ I 2 3. Second-order Resonance(p,q)→(n,n) Two different feeding modes, p,q fill the same mode, n. The second order reso- nant Hamiltonian is H n =ω n I n − μ 0 2~ 2 I 2 n + 2μ 0 ~ 2 I n (I p I q ) 1/2 cos(2θ n −θ p −θ q ) (3.20) In Table 3.1 the important parameters and definitions are listed for the three classes of resonances. The parameters include: the geometric mean of the fixed modes, ¯ I, the bare detuning, Δ 0 , and the frequency of the drive ν. For the nonlinear detuning, the second equality holds when all the fixed modes have the same population. For a generic driving term of each type, the bare detuning,Δ 0 is always negative for 1-R terms, always 75 Table 3.2: Bare Detuning: Generic Values Order Modes Bare Detuning First n+m,n+l→n,n+l+m Δ 01 ≡(n+m) 2 +(n+l) 2 −(n+m+l) 2 −n 2 =−2ml 1R n+m,n+m→n,n+2m Δ 01r ≡2(n+m) 2 −(n+2m) 2 −n 2 =−2m 2 Second n+m,n−m→n,n Δ 02 ≡(n+m) 2 +(n−m) 2 −2n 2 =2m 2 positive for 2nd order terms and can be positive or negative for 1st order terms (see Table 3.2). Next we will study the resonances in each class of Hamiltonians. 3.2.2 First-Order Resonances in BOA Consider a first-order resonance wheren+r = p+q,p,q 6= n,r. Fix the populations of the feeding modes, p,q and the filling mode r, so that the only variables are the population and angle of mode n. The frequencies of modes p,q,r are fixed to their values in the unperturbed Hamiltonian (3.2). Solving ˙ θ n = ∂H 0 ∂I n for the angle gives θ n =(ω n − μ 0 ~ 2 ¯ I n )t≡Ω n ( ¯ I n )t. (3.21) 76 In this approximation, which we call the Born-Oppenheimer approximation (BOA), we keep the eight terms corresponding to (p,q) → (r,n) from the full Hamiltonian. The Hamiltonian for first-order resonances in the Born-Oppenheimer approximation is: H n =ω n I n − μ 0 2~ 2 I 2 n + 4μ 0 ~ 2 (I n I p I q I r ) 1/2 cos(θ n +θ r −θ p −θ q ) =ω n I n − μ 0 2~ 2 I 2 n + 4μ 0 ~ 2 (I n I p I q I r ) 1/2 cos(θ n −νt+φ 1 ) (3.22) whereν≡Ω p +Ω q −Ω r . Now we let ¯ I 1 ≡( ¯ I p ¯ I q ¯ I r ) 1/3 , divide by ˜ ω 1 , set~=1 and use ξ =μ 0 /(~˜ ω 1 ) to get h n ≡ H n ˜ ω 1 = ω n ˜ ω 1 I n − ξ 2 I 2 n +4ξ ¡ ¯ I 1 ¢ 1/2 cos(θ n −νt+φ 1 ). (3.23) We seek a canonical transformation to a rotating reference frame (I n ,θ n → ˜ I n , ˜ θ n ), where the new angle, ˜ θ n is slowly varying. Introducing the type 2 generating function Φ=(θ n −νt+φ 1 ) ˜ I n , (3.24) the new canonical variables are given by I n ≡ ∂Φ ∂θ n = ˜ I n (3.25) ˜ θ n ≡ ∂Φ ∂ ˜ I n =θ n −νt+φ 1 . (3.26) (3.27) The Hamiltonian, which transforms according to ˜ h n =h n + ∂Φ ∂t , becomes e h n =(ω n −ν 1 )/˜ ω 1 ˜ I n − ξ 2 ˜ I 2 n +4ξ ¯ I 3/2 1 ˜ I 1/2 n cos ˜ θ n . (3.28) 77 The corresponding equations of motion are ˙ ˜ I n =− ∂ ˜ h n ∂ ˜ θ n =4ξ e I 1/2 n ¡ ¯ I 1 ¢ 3/2 sin ˜ θ n (3.29) ˙ ˜ θ n = ∂ ˜ h n ∂ ˜ I n =−Δ 1 −ξ ˜ I n +2ξ ¯ I 3/2 1 ˜ I 1/2 n cos ˜ θ n , (3.30) whereΔ 1 ≡ν 1 −ω n . Resonance Condition Resonance occurs at the stationary points, that is ˙ ˜ I n ¯ ¯ ¯ ( ˜ I ∗ n , ˜ θ ∗ n ) =0 ⇒ 4ξ( ˜ I ∗ n ) 1/2 ¯ I 3/2 1 sin ˜ θ ∗ n =0 (3.31) ˙ ˜ θ n ¯ ¯ ¯ ( ˜ I ∗ n , ˜ θ ∗ n ) =0 ⇒ −Δ 1 −ξ ˜ I n +2ξ ¯ I 3/2 1 ( ˜ I ∗ n ) 1/2 cos ˜ θ ∗ n =0. (3.32) The first condition is satisfied by (a) ˜ I ∗ n =0 or (b)sin ˜ θ ∗ n =0. For ˜ I n =0 the phase is not well-defined and corresponds to a stationary point, even though (3.32) is not satisfied. The resonances of this model will correspond to taking sin ˜ θ ∗ n = 0 and solving (3.32) for ˜ I ∗ n , −Δ 1 ˜ I 1/2 n −ξ ˜ I 3/2 n =∓2ξ ¯ I 3/2 1 . (3.33) Squaring both sides and rearranging gives, ˜ I 3 n +2 Δ 1 ξ ˜ I 2 n + Δ 2 1 ξ 2 ˜ I n −4 ¯ I 3 1 =0. (3.34) 78 The solution to this cubic equation gives the fixed points, ˜ I ∗ n1 =− 2Δ 1 3ξ +A+B (3.35) ˜ I ∗ n2 =− 2Δ 1 3ξ − 1 2 (A+B)+i √ 3 2 (A−B) (3.36) ˜ I ∗ n3 =− 2Δ 1 3ξ − 1 2 (A+B)−i √ 3 2 (A−B), (3.37) where A= 2 4 µ Δ 1 3ξ ¶ 3 +2 ¯ I 3 1 +2 ¯ I 3/2 1 Ã µ Δ 1 3ξ ¶ 3 + ¯ I 3 1 ! 1/2 3 5 1/3 (3.38) B = 2 4 µ Δ 1 3ξ ¶ 3 +2 ¯ I 3 1 −2 ¯ I 3/2 1 Ã µ Δ 1 3ξ ¶ 3 + ¯ I 3 1 ! 1/2 3 5 1/3 . (3.39) The real and imaginary parts of the roots are plotted in Fig. 3.2 for ¯ I 1 =1/3 andΔ 01 = −4, which corresponds to(n+2,n+1)→ (n,n+3), which is the negative detuning closest to zero. For the initial condition studied by the previous perturbation expansion, where a block of consecutive modes are initially occupied, the bare detuning is always negative. Thus, we still study the case of negative detuning although it is possible to have a first-order resonance with positive or negative values. The value ofξ the separates regions with three real solutions to one real solution is given by ³ Δ 1 3ξ c ´ 3 + ¯ I 3 1 =0 or Δ 1 3ξ c = Δ 01 −ξ c ¯ I 1 3ξ c =− ¯ I 1 ⇒ ξ c2 ≡− Δ 01 2 ¯ I 1 (3.40) As will be shown below there is another relevant critical value of ξ at lower values, so this critical value is defined as ξ c2 . There are three physical fixed points ( ˜ I ∗ n ≥ 0 for ξ≤ξ c2 and one physical fixed point forξ >ξ c2 . 79 0 5 10 15 20 ξ 0 1 2 3 4 5 Re I ξ c =6 0 5 10 15 20 ξ − 0.4 0.0 0.4 Im I ξ c =6 I ∗ n1 I ∗ n2 I ∗ n3 Figure 3.2: Real and imaginary parts of the fixed points of first-order Hamiltonian (3.28). Δ 01 =−4, ¯ I 1 =1/3. Separatrix Width There are two separatricies: one is associated with ˜ h n = 0 another passes through the fixed point( ˜ I ∗ n3 , ˜ θ ∗ n ±π) and is defined by the contour ˜ h n = ˜ h n ( ˜ I ∗ n3 ,±π). In this section we calculate the maximum value of ˜ I n for the separatrix defined by ˜ h n = 0. The maximum value of the separatrix occurs at ˜ θ n =±mπ for integerm. We thus solve ˜ h n =0 forcos( ˜ θ n )=±1. −Δ 1 ˜ I n − 1 2 ξ ˜ I 2 n =∓2ξ ¯ I 3/2 1 ˜ I 1/2 n . (3.41) Squaring both sides and rearranging gives, ˜ I 3 n +4 Δ 1 ξ ˜ I 2 n +4 Δ 2 1 ξ 2 ˜ I n −64 ¯ I 3 1 =0. (3.42) 80 Figure 3.3: Action-angle phase-space plots for variables ˜ I n , ˜ θ n for the first-order reso- nant Hamiltonian (3.28). (n+2,n+1)→(n,n+3)(Δ 01 =−4). ¯ I 1 =1/3. (a)ξ =1 (b) ξ = 3.5 (c) ξ = ξ c1 = 4.3 (d) ξ = 7 > ξ c2 . The black contour line corresponds to ˜ h n =0 and the white contour line to ˜ h n ( ˜ I n , ˜ θ n )= ˜ h n ( ˜ I ∗ n3 ,−π). Solving this cubic equation gives the roots, ˜ I sep1 =− 4Δ 1 3ξ +A+B (3.43) ˜ I sep2 =− 4Δ 1 3ξ − 1 2 (A+B)+i √ 3 2 (A−B) (3.44) ˜ I sep3 =− 4Δ 1 3ξ − 1 2 (A+B)−i √ 3 2 (A−B) (3.45) 81 where A= 2 4 8 µ Δ 1 3ξ ¶ 3 +4 ¯ I 3 1 +16 √ 2 ¯ I 3/2 1 Ã µ Δ 1 3ξ ¶ 3 +2 ¯ I 3 1 ! 1/2 3 5 1/3 (3.46) B = 2 4 8 µ Δ 1 3ξ ¶ 3 +4 ¯ I 3 1 −16 √ 2 ¯ I 3/2 1 Ã µ Δ 1 3ξ ¶ 3 +2 ¯ I 3 1 ! 1/2 3 5 1/3 . (3.47) The second critical value ofξ is given by ³ Δ 1 3ξ c ´ 3 +2 ¯ I 3 1 =0 or Δ 1 3ξ c = Δ 01 −ξ c ¯ I 1 3ξ c =− 3 √ 2 ¯ I 1 ⇒ ξ c1 ≡− Δ 01 (1+3 3 √ 2) ¯ I 1 (3.48) The maximum value of the separatrix occurs at ˜ θ n =±π, (cos( ˜ θ n ) =−1) forξ ≤ ξ c1 and at ˜ θ n ) = 0±2π, (cos( ˜ θ n = +1) for ξ ≥ ξ c1 . In Fig. 3.4 the real and imaginary parts of these solutions are plotted. From the fixed points of the Hamiltonian, solutions of ˜ h n = 0 and sample contour plots, the phase-space of the first-order Hamiltonian can be characterized as follows: ξ <ξ c1 There are two resonances. One is at( ˜ I n , ˜ θ n =) = ( ˜ I ∗ n2 ,±π) with the separatrix defined by the contour at ˜ h n = 0. This resonance can be seen in the phase-space diagrams of Fig. 3.3(a)-(b). The height of the separatrix is given by I sep2 which increases from zero atξ = 0 to a maximum value ˜ I n2−max = 2 3 √ 2 ¯ I 1 atξ = ξ c2 as seen in Fig. 3.4. For ¯ I 1 = 1/3 the maximum value is ˜ I n2−max ≈ 0.84. The second resonance is given by ( ˜ I ∗ n1 , ˜ θ n = 0). The separatrix is defined by the contour that passes through the saddle point( ˜ I ∗ n3 ,±π). This resonance diverges asξ→ 0 and remains above 1 in this parameter range is thus not physically accessible. However, the lower separatrix comes into the accessible phase space as seen in Fig. 3.3(b).Due to the normalization, the occupation of ˜ I n is bounded by one. 82 0 5 10 15 20 ξ 0 2 4 6 8 10 Re I ξ c =4.3 0 5 10 15 20 ξ − 1 0 1 Im I ξ c =4.3 I sep1 I sep2 I sep3 Figure 3.4: Real and imaginary parts of the solutions of ˜ h n = 0 for the first-order Hamiltonian (3.28). Δ 01 =−4(n+1,n+2)→(n,n+3), ¯ I 1 =1/3. Forξ≤ξ c1 , ˜ I sep2 gives the separatrix for the the resonance at( ˜ I ∗ n2 , ˜ θ n =±π). Forξ≥ξ c1 , ˜ I sep1 gives the separatrix for the the resonance at( ˜ I ∗ n1 , ˜ θ n =0). Larger values of ˜ I n are plotted in Fig. 3.3 to gain a of deeper understanding of the phase space. ξ =ξ c1 At the first critical point, the separatricies of the two resonances overlap as show in Fig. 3.3(c). At this critical point, the lower bound of the separatrix defined by ˜ h n = ˜ h n ( ˜ I ∗ n3 ,−π) touches zero. Starting with an arbitrarily small occupation, the entire physical phase space becomes accessible to ˜ I n once the two resonances touch. ξ c1 <ξ <ξ c2 The two resonances at( ˜ I ∗ n1 ,0) and( ˜ I ∗ n2 ,±π) remain, with the separatrix of the first defined by ˜ h n = 0 and the separatrix of the second resonance defined by ˜ h n = ˜ h n ( ˜ I ∗ n3 ,−π). The association of the separatricies have switched from the case whereξ <ξ c1 . 83 ξ≥ξ c2 At the second critical point the fixed points( ˜ I ∗ n2 ,±π) and( ˜ I ∗ n3 ,±π) merge and annihilate and a single resonance at( ˜ I n , ˜ θ n )=( ˜ I ∗ n1 ,0) persists. Aboveξ c , there is a single resonance at( ˜ I n , ˜ θ n ) = ( ˜ I ∗ n1 ,0), as show in Fig. 3.3(d) and the separatrix is defined by the contour ˜ h n = 0. The height of the separatrix is given by ˜ I sep1 in Fig. 3.4. 3.2.3 1-R Resonances in BOA.Ω p =Ω q Consider a first-order resonance where the feeding modes are the same, n +r = 2p. In this case, we keep four terms from the full Hamiltonian, and the Hamiltonian for the first-order resonances (1-R) in H n =ω n I n − μ 0 2~ 2 I 2 n + 2μ 0 ~ 2 ¡ I n I 2 p I r ¢ 1/2 cos(θ n +θ r −2θ p ) (3.49) After making the transformation to a rotating reference frame, introducingν 1r =2Ω p − Ω r , fixing the populations of modes I p ,I r , and dividing by ω n , the Hamiltonian, ˜ h n ≡ e H n /ω n becomes ˜ h n =−Δ 1r ˜ I n − ξ 2 ˜ I 2 n +2ξ ˜ I 1/2 n ¯ I 3/2 1r cos ˜ θ n (3.50) whereΔ 1r ≡ν 1r −ω n . The corresponding equations of motion are ˙ ˜ I n =− ∂ ˜ h n ∂ ˜ θ n =2ξ ˜ I 1/2 n ¯ I 3/2 1r sin ˜ θ n (3.51) ˙ ˜ θ n = ∂ ˜ h n ∂ ˜ I n =−Δ−ξ ˜ I n +2ξ Ã ¯ I 3/2 1r ˜ I n ! 1/2 cos ˜ θ n (3.52) The Hamiltonian for the 1R resonances differs from the first-order Hamiltonian only by prefactors. The resonances of the 1-R Hamiltonian will thus be very similar in form 84 to those of the first-order resonances, and can be determined from the previous analysis by a proper rescaling of ¯ I 1 . 3.2.4 Second Order Resonances in BOA For the second order case, we consider a single filling moden and restrict the resonance to modes that satisfyp+q = 2n. The Hamiltonian for a single second order resonance in the Born Oppenheimer approximation is H n =ω n I n − μ 0 2~ 2 I 2 n + 2μ 0 ~ 2 I n (I p I q ) 1/2 cos(2θ n −θ p −θ q ) =ω n I n − μ 0 2~ 2 I 2 n + 2μ 0 ~ 2 I n ¡ ¯ I p ¯ I q ¢ 1/2 cos(2θ n −ν 2 t) (3.53) where in the second line we setI p andI q to their unperturbed values and ν 2 ≡Ω p ( ¯ I p )+Ω q ( ¯ I q )=ω p +ω q − μ 0 ~ 2 ¡ ¯ I p + ¯ I q ¢ . Next we make a canonical transformation to a rotating reference frame (I n ,θ n → ˜ I n , ˜ θ n ), through the type 2 generating function Φ= 1 2 (2θ n −ν 2 t+φ 1 ) ˜ I n . (3.54) The new canonical variables are determined by I n = ∂Φ ∂θ n = ˜ I n (3.55) ˜ θ n = ∂Φ ∂ ˜ I n = 1 2 (2θ n −ν 2 t+φ 1 ), (3.56) 85 and the Hamiltonian transforms according to ˜ H n ( ˜ I n ˜ θ n )=H n (I,θ n )+ ∂Φ ∂t becomes e H n = µ ω n − 1 2 ν 2 ¶ ˜ I n − μ 0 2~ 2 ˜ I 2 n + 2μ 0 ~ 2 ˜ I n ¡ ¯ I p ¯ I q ¢ 1/2 cos2 ˜ θ n (3.57) Next we divide by ˜ ω 1 , which is equivalent to rescaling time by a factorτ = ˜ ω 1 t, intro- duceξ≡μ 0 /~˜ ω 1 and set~=1. Furthermore, we define ¯ I 2 ≡( ¯ I p ¯ I q ) 1/2 as the geometric mean of the filling modes ¯ I p , ¯ I q andΔ 2 ≡( 1 2 ν 2 −ω n )/˜ ω 1 . e h n ≡ e H n ˜ ω 1 =−Δ 2 ˜ I n − ξ 2 ˜ I 2 n +2ξ ˜ I n ¯ I 2 cos2 ˜ θ n (3.58) The corresponding equations of motion are ˙ ˜ I n =− ∂ e h n ∂ ˜ θ n =4ξ ˜ I n ¯ I 2 sin2 ˜ θ n (3.59) ˙ ˜ θ n = ∂ e h n ∂ ˜ I n =−Δ 2 −ξ ˜ I n +2ξ ¯ I 2 cos2 ˜ θ n . (3.60) If we furthermore assume that the filling modes have equal populations, ¯ I p = ¯ I q = ¯ I 2 and that the dispersion is quadratic,ω n =n 2 ˜ ω 1 then the detuning reduces to Δ 2 = 1 2˜ ω 1 (ν−2ω n )= 1 2˜ ω 1 £ (p 2 +q 2 −2n 2 )˜ ω 1 −ξ ¡ ¯ I p + ¯ I q ¢¤ (3.61) = Δ 02 2 −ξ ¯ I 2 . (3.62) 86 Resonances Conditions Unlike the first-order resonances there is a simple expression for the fixed points of the second order resonances. The conditions for the fixed points are ˙ ˜ I n | ( ˜ I ∗ n , ˜ θ ∗ n ) = 0, ˙ ˜ θ n | ( ˜ I ∗ n , ˜ θ ∗ n ) =0. ˙ ˜ I n ¯ ¯ ¯ ( ˜ I ∗ n , ˜ θ ∗ n ) =0 ⇒ 4ξ ˜ I ∗ n ¯ I 2 sin2 ˜ θ ∗ n =0 (3.63) ˙ ˜ θ n ¯ ¯ ¯ ( ˜ I ∗ n , ˜ θ ∗ n ) =0 ⇒ −Δ 2 −ξ ˜ I ∗ n +2ξ ¯ I 2 cos2 ˜ θ ∗ n =0 (3.64) 1. Case 1: ˜ I ∗ n =0. The condition for ˙ ˜ I n = 0 is satisfied by ˜ I ∗ n = 0 and the condition on ˜ θ ∗ n such that ˙ ˜ θ n =0 is cos2 ˜ θ ∗ n = Δ 2 2ξ ¯ I 2 = Δ 20 4ξ ¯ I 2 − 1 2 . Note that ifI n = 0, the phase of the moden is not well-defined. In this case, the motion of the action variable is always stationary. Thus a second-order resonant Hamiltonian can never populate an initially unoccupied mode. However, once a small seed is present, the mode can grow and enter the dynamics. 2. Case 2: sin2 ˜ θ ∗ n =0 (cos2 ˜ θ ∗ n =±1) The condition for ˙ ˜ I n = 0 is satisfied bysin2 ˜ θ ∗ n = 0 and the condition on ˜ I ∗ n such that ˙ ˜ θ n =0 is ˜ I ∗ n =− Δ 2 ξ ±2 ¯ I 2 =− Δ 02 2ξ + ¯ I 2 ±2 ¯ I 2 (for cos2 ˜ θ ∗ n =±1). 87 Classification of Stationary Points Following Tabor (1989) we outline a method for determining the stability of the sta- tionary points by linearizing about the fixed points. Consider a general second order differential equation, written as a pair of first-order differential equations, ˙ x=f(x,y) (3.65) ˙ y =g(x,y) (3.66) that has stationary points{(x ∗ ,y ∗ )|f(x ∗ ,y ∗ )=0, g(x ∗ ,y ∗ )=0}. The stability of these points can be determined by linearizing about the stationary points. 0 @ Δ˙ x Δ˙ y 1 A = 0 @ ∂f ∂x | x ∗ ,y ∗ ∂f ∂y | x ∗ ,y ∗ ∂g ∂x | x ∗ ,y ∗ ∂g ∂y | x ∗ ,y ∗ 1 A 0 @ Δx Δy 1 A . (3.67) The solution to these equations of motion, ˙ x = Ax, isx = c 1 e λ 1 t v 1 +c 2 e λ 2 t v 2 where λ j is an eigenvalue ofA with corresponding eigenvectorv j . The stability of the fixed points is determined from the eigenvalues,λ j = a j +ib j , a j ,b j ∈R. The fixed points are classified as follows, 1. center: a 1,2 =0 2. spiral: b 1,2 6=0. Stable fora 1,2 <0. Unstable fora 1,2 >0. 3. node: b 1,2 =0. Stable fora 1,2 <0. Unstable fora 1,2 >0 4. saddle point: b 1,2 =0, a 1 <0, a 2 >0 The linearized equations of motion for the second order Hamiltonian (3.58) are A= 0 @ − ∂ 2˜ h n ∂ ˜ θ n ∂ ˜ I n | ˜ I ∗ n , ˜ θ ∗ n − ∂ 2˜ h n ∂ ˜ θ 2 n | ˜ I ∗ n , ˜ θ ∗ n ∂ 2˜ h n ∂ ˜ I 2 n | ˜ I ∗ n , ˜ θ ∗ n ∂ 2˜ h n ∂ ˜ I n ∂ ˜ θ n | ˜ I ∗ n , ˜ θ ∗ n 1 A = 0 @ 4ξ ¯ I 2 sin2 ˜ θ ∗ n 8ξ ˜ I n ¯ I 2 cos2 ˜ θ ∗ n −ξ −4ξ ¯ I 2 sin2 ˜ θ ∗ n 1 A . (3.68) 88 For each set of fixed points, we calculate the eigenvalues of the matrix, given by det|A−λI|=0 to determine the type of fixed point. 1. Fixed Point 1: ( ˜ I ∗ n , ˜ θ ∗ n )=(0, 1 2 cos −1 ¡ Δ 2 /(2ξ ¯ I 2 ) ¢ ) The eigenvalues are given by λ=±4ξ ¯ I 2 " 1− µ Δ 2 2ξ ¯ I 2 ¶ 2 # 1/2 =±2ξ " ¡ 2 ¯ I 2 ¢ 2 − µ Δ 2 ξ ¶ 2 # 1/2 . (3.69) These eigenvalues are real for|Δ 2 | < 2ξ ¯ I 2 . Given that the bare detuning,Δ 02 , is always positive for the second-order Hamiltonian (see Table 3.2) and the typical population of the fixed modes, ¯ I 2 is also positive, this condition leads to a critical value of the nonlinearity parameter,ξ, ξ c = Δ 02 6 ¯ I 2 (3.70) Forξ≥ ξ c this fixed point is a saddle point and the separatrix passes through this point. 2. Fixed Point 2: ( ˜ I ∗ n , ˜ θ ∗ n )=(−Δ 2 /ξ+2 ¯ I 2 ,mπ), m= integer λ=±2 √ 2ξ ¯ I 1/2 2 · Δ 2 ξ −2 ¯ I 2 ¸ 1/2 . (3.71) For ˜ I ∗ n = −Δ 2 /ξ + 2 ¯ I > 0 the fixed point is a center, and for ˜ I ∗ n < 0 it is a saddle point. We exclude the unphysical values ˜ I ∗ n < 0. The stationary points for ˜ I n > 0 is atcos(2 ˜ θ n ) = +1, are the physical resonances for the second order driving terms, with resonant value given by I res =− Δ 02 2ξ +3 ¯ I 2 . (3.72) 89 Furthermore the condition for the existence of this type of resonance isξ >ξ c . 3. Fixed Point 3: ( ˜ I ∗ n , ˜ θ ∗ n )=(−Δ 2 /ξ+2 ¯ I,(2m+1)π/2), m= integer Forcos(2 ˜ θ ∗ n )=−1, the value of ˜ I n is always negative and thus not relevant to the current analysis. Separatrix Height The separatrix passes through the fixed point at ˜ I ∗ n =0, ˜ θ ∗ n = 1 2 cos −1 ¡ Δ 2 /(2ξ ¯ I 2 ) ¢ when ˜ θ ∗ n is real. The contour plots of ˜ h n confirm that the separatrix passes through this fixed point. The maximum value of ˜ I n along the separatrix can be found by solving for ˜ I n when ˜ h n =0. Given ˜ h n =−Δ 2 ˜ I n − ξ 2 ˜ I 2 n +2ξ ˜ I n ¯ I 2 cos2 ˜ θ n =0, ˜ I n =− 2Δ 2 ξ +4 ¯ I 2 cos2 ˜ θ n =− Δ 02 ξ +2 ¯ I 2 (1+2cos2 ˜ θ n ). The maximum height of the separatrix occurs at ˜ θ n =0 and is given by ˜ I sep =− Δ 02 ξ +6 ¯ I 2 . (3.73) In Fig. 3.5(a)-(b) the phase plots for the second order Hamiltonian is plotted for nonlinearities above and below the critical value. For ξ ≤ ξ c , there is no resonance in the phase-space that corresponds to physical values of ˜ I n . In Fig. 3.5(b), whereξ > ξ c , there is a resonance and the separatrix is defined by the contour ˜ h n =0 which is plotted in black. The resonance values and separatrix height are labeled by the values given by (3.72) and (3.73), respectively. 90 Figure 3.5: Phase-space plots for the second order Hamiltonian. (a) ξ = 0.5 (κ = 0.16 forN s = 11). There is no resonance in the range 0≤ ˜ I n ≤ 1 forξ < ξ c = 1. (b) ξ = 1.5 (κ = 0.49 forN s = 11). ¯ I 2 = 1/3, Δ 02 = 2. Forξ > ξ c the resonant value is given by (3.72) and the separatrix with is given by (3.73). The values of the action variable at resonance and the separatrix height are plotted in Fig. 3.6 for the two second-order Hamiltonians as a function of the nonlinearity, ξ for the two lowest bare detunings. The two detunings are Δ 02 = 2, corresponding to (n+1,n−1) → (n,n) and ξ c = 2 ¯ I 2 /3, and Δ 02 = 8 for (n+2,n−2) → (n,n) 91 0 2 4 6 8 10 ξ − 0.5 0.0 0.5 1.0 1.5 I ξ c1 ξ c2 I res , 0 = 2 I sep , 0 = 2 I res , 0 = 8 I sep , 0 = 8 0.0 0.5 1.0 1.5 2.0 2.5 3.0 κ − 0.5 0.0 0.5 1.0 1.5 Figure 3.6: Resonant values and separatrix width for second-order resonance forΔ 02 = 2[(n+1,n−1)→(n,n)] and forΔ 02 =8[(n+2,n−2)→(n,n)]. ¯ I =0.33. Vertical lines: Δ 20 = 2: ξ c1 = 1/3 ¯ I 2 = 1, Δ 02 = 8: ξ c2 = 4/3 ¯ I 2 = 4. The upper x-axis gives the corresponding values ofκ forN s =11. andξ c = 4 ¯ I 2 /3. It is clear from the plot that for the second order resonant Hamiltonian there there is a critical value of the nonlinearity such that below that value there are no physical resonances present. 3.3 Failure of Chirikov’s Criterion The Chirikov criterion for the onset of global chaos is governed by the ratio of the width of the separatrix and the distance between the resonances for individual resonances. When the width of the separatricies becomes comparable to the distance between the resonances, the Chirikov parameterK satisfies the inequality, K≡ separatrix width distance between resonances = ΔI I res1 − ¯ I res2 &1, (3.74) 92 and the system is predicted to be chaotic. We should note that the study of resonances in the BHM is quite different from the case of the kicked rotor. In particular, it is not possible to independently vary the fre- quency of the drive and the strength of the drive in the BHM. In addition, in the canonical transformation to the rotating reference frame the action variable is unchanged, ˜ I n =I n , so that the resonant values of ˜ I n are not well-spaced in the original action coordinates. Furthermore in the BHM, for the first order resonances, a single driving term has reso- nances at multiple action values. Note that for a generic driving term, (n+p,n+m−p) → (n,n+m), the bare detuning is given by Δ 0 =(n+p) 2 +(n+m−p) 2 −n 2 −(n+m) 2 = =n 2 +2np+p 2 +n 2 +m 2 +p 2 −2np−2mp+2nm − £ n 2 +n 2 +2nm+m 2 ¤ =2p 2 −2mp, (3.75) which is always constant. We present several ways to deduce the criterion governing the threshold for the mean-field Bose-Hubbard model and see that the various methods are in agreement. Onset of Second-Order Resonances First we consider the case of the second-order resonance, where the analytic expression is simple. From the previous analysis for the second order resonance, the separatrix width is given by, I sep = −Δ 0 ξ +6 ¯ I. (3.76) 93 The distance between resonances is proportional toξ −1 , ΔI res =I res (Δ b )−I res (Δ a )= Δ 0a −Δ 0b ξ . (3.77) For the lowest two resonances, a naive comparison of the separatrix height and distance between the resonances gives, K≡ I sep ΔI res = −2/ξ+6 ¯ I 6/ξ (3.78) and chaos exists for K&1⇒ ¯ Iξ& 5 3 (3.79) However, but looking at Fig. 3.6 , we see that once the second resonance appears, the two resonances overlap. Thus we take the appearance of the second resonance as the criterion for the second-order case, which gives ξ >ξ c2 = 4 3 ¯ I or ¯ Iξ& 4 3 . (3.80) Overlap of resonances within first-order resonant Hamiltonian For the first-order resonance, the resonances of ˜ I ∗ n2 grow from zero and exist for any ξ > 0. However, the separatrix of the resonance is small and thus even though the resonances may over- lap, the population is still expected to be confined in a narrow region of phase-space. However, there are additional resonances at larger values of ˜ I n , which are initially inac- cessible and move down from above. For the first-order resonance we use criterion that the two separatricies of the same first-order Hamiltonian overlap - i.e. all of phase space 94 is accessible to a single mode. The condition for is being able to explore the entire phase-space for a givenI n is thus ξ≥ξ c1 ≡− Δ 01 (1+3 3 √ 2) ¯ I 1 (3.81) which is equivalent to ¯ Iξ≥− Δ 01 (1+3 3 √ 2) . (3.82) For the lowest bare detuning,Δ 01 =−4, the criterion is ¯ Iξ≥ 4 (1+3 3 √ 2) ≈.0.84 (3.83) which is the comparable to the results of for the second order resonance. Dimensional Analysis An alternate way of coming to this conclusion is to consider dimensional grounds. The starting point is to assume: 1. Typical mode occupation: ¯ I ∼ 1 Δn 2. Resonant approximation: only include resonant terms in equations of motion. 3. Quadratic dispersion instead of cosine: ω n = ˜ ω 1 n 2 These are the same ingredients for the previous analytic analysis. A quadratic dispersion leads to translational invariance. The momentum distribution can be shifted with no effect. The resonant equations are invariant under a shift in ’n’, thus the frequency scale is set by ˜ ω 1 . As a consequence of the resonant approximation, the equations of motion only couple neighboring modes. Thus,μ 0 andΔn can only enter as μ 0 Δn . Thus the only parameters areμ 0 /Δn and ˜ ω 1 . The only dimensionless combination of these quantities isμ 0 /(~˜ ω 1 Δn)=ξ/Δn, so this parameter must be what governs the Chirikov threshold. 95 Signficantly, all of these methods give the same functional form of the criterion and are identical up to numerical factors. 3.3.1 Thermodynamic Limit Next we consider how K scales in the thermodynamic limit, (N s → ∞ while U = const,J = const,N a /N s = const). Additionally, the width of the momentum dis- tribution remains fixed, Δn/N s = const. In this limit, the nonlinearity parameter κ is constant while ξ ∼ N 2 s . The typical occupation of the fixed modes is given by ¯ I ∼(Δn) −1 ∼(N s ) −1 We find that the Chirikov parameter, K≈ξ ¯ I ∼ N 2 s N s ∼N s (3.84) diverges in the thermodynamic limit, predicting that there is no threshold in that limit andK is always greater than one indicating that the system is always chaotic. However, this is not what we have observed numerically. At a minimum Chirikov predicts the threshold to scales linearly with the size of the system and our numerics indicate that the threshold depends on parameters that are independent of the system size. 3.3.2 Continuous Limit To check this threshold, we also consider the continuum limit, in which the length of the system, normalization and interaction parameter are fixed, while the distance between lattice sites goes to zero: L,N a ,μ 0 = const, a→0,N s →∞. K∼ U J µ N a Δn ¶ N s = μ 0 N s N a 2mL 2 ~ 2 N 2 s µ N a Δn ¶ N s = 2μ 0 mL 2 ~ 2 µ 1 Δn ¶ ∼ 1 N s →0 (3.85) 96 In the continuum limit, the number of modes in the initial momentum distribution,Δn, must diverge, so that N a = Δn| ¯ ψ n | 2 remains constant. Alternatively we can say that Δn/N s is fixed. From either perspective, it is clear thatK vanishes in the continuum limit, predicting regular motion, as expected due to the known integrability of the con- tinuous NLS equation. This leads to the question: Why does the Chirikov analysis fail in this case? Several possible explanations are: 1. non-quadraticity of the spectrum 2. break down of the resonant approximation 3. interferences between resonances. First there are possible corrections due to the lattice and that we have assumed a quadratic dispersion instead of a cosine dispersion. We conjecture that the Chirikov analysis is incorrect because the the resonant approximation is wrong and off-resonant terms are significant. That is the motion is not dominated by single driving terms. It is also possible that the single resonance approximation is incorrect and a multi-resonance model is necessary for recovering the correct scaling. 97 Chapter 4 Conserved Quantities of the Ablowitz-Ladik Lattice From the numerical work on the Bose-Hubbard Model, we have seen regions where although individual trajectories are chaotic the system does not relax to the expected thermal state. Additionally there are regions in the parameter space of nonlinearity κ and energy-per-particle, ² T where the system relaxes even though chaos is not present. Additionally it has been found that the slow relaxation times for individual states with the same nonlinearity and energy-per-particle vary widely. What is the origin of these phenomena? There are at least three nearby integrable systems of the BHM. These are the noninteracting case (κ = 0), the continuum limit, which is the continuous nonlin- ear Schr¨ odinger equation and the integrable discrete nonlinear Schr¨ odinger equation (IDNLS) also know as the Ablowitz-Ladik (AL) lattice, which we discuss in this sec- tion. Looking at the relaxation of individual realizations, it appears that initial states which have a higher quasi-momentum also have a higher final spectral entropy. The total quasi-momentum is not a conserved quantity of the BHM. It is however, a con- served quantity of all three nearby integrable systems. In some cases the normalized spectral entropy many not vanish because the system has a very slow relaxation time. In other cases it may be that it will never fully relax because of the nearby conserved quantities. Are the slow relaxation times governed by the conserved quantities of the nearby integrable systems? In the regions where it is clear that steady-state is not the 98 thermal state, what governs the steady-state? This evidence suggests that the integrals of motion of the nearby conserved quantities play a role in the relaxation dynamics of the BHM. In this section, we derive the conserved quantities of the AL. 4.1 Integrable Discrete Nonlinear Schr¨ odinger (IDNLS) Equation The continuous Nonlinear Schr¨ odinger (NLS) Equation is given by, i˙ q+q xx +σ|q| 2 q =0 (4.1) where ˙ q = dq/dt and q x = dq/dx, is know to be completely integrable. The corre- sponding Hamiltonian is H(q,q ∗ )=− Z L 0 · q x q ∗ x + 1 2 σ(qq ∗ ) 2 ¸ dx, (4.2) with canonical pairs q,q ∗ . Periodic boundary conditions are assumed. The mean-field Bose-Hubbard model is a discretization of the NLS equation. In real space, the Hamil- tonian can be written as: H =− X n ¡ q n q ∗ n+1 +q ∗ n q n+1 −2|q n | 2 ¢ + σ 2 X n |q n | 4 (4.3) with the Poisson brackets {q m ,q ∗ n }=iδ m,n , {q m ,q n }={q ∗ m ,q ∗ n }=0. (4.4) 99 The equations of motion, which are given by ˙ q n ={H,q n }, are i˙ q n =−(q n+1 +q n−1 −2q n )+σ|q n | 2 q n . (4.5) This equation, also know as the Diagonal Discrete Nonlinear Schr¨ odinger (DDNLS) Equation, does not preserve the integrability of the continuous case. An alternate dis- cretization, which is integrable (Ablowitz and Ladik, 1976), has equations of motion: i˙ q n =−(q n+1 +q n−1 −2q n )+σ|q n | 2 (q n+1 +q n−1 )/2 (4.6) and is suitably called the Integrable Discrete Nonlinear Schr¨ odinger (IDNLS) Equation, also know as the Ablowitz-Ladik lattice (AL). This equation can be derived from the Hamiltonian (Scharf and Bishop, 1991; Herbst et al., 1994) H =− X n ¡ q n q ∗ n+1 +q ∗ n q n+1 ¢ − 4 σ X n ln ³ 1− σ 2 |q n | 2 ´ (4.7) with the nonstandard Poisson brackets {q m ,q ∗ n }=i(1− σ 2 |q n | 2 )δ m,n , {q m ,q n }={q ∗ m ,q ∗ n }=0. (4.8) The equations of motion derived in the usual way, ˙ q m ={H,q m }. Using the following properties of Poisson brackets, {a,bc}=b{a,c}+{a,b}c (4.9) {f(a),b}= ∂f ∂a {a,b}, (4.10) 100 the AL equation can be derived from the equations of motion of Hamiltonian (4.7), ˙ q m ={H,q m } =− X n ¡ {q n q ∗ n+1 ,q m }+{q ∗ n q n+1 ,q m } ¢ − 4 σ X n {ln ³ 1− σ 2 |q n | 2 ´ ,q m } =− X n ¡ q n {q ∗ n+1 ,q m }+{q ∗ n ,q m }q n+1 ¢ − 4 σ X n · ∂ ∂q ∗ n ln ³ 1− σ 2 |q n | 2 ´ ¸ {q ∗ n ,q m } =− X n h q n ·(−i)(1− σ 2 |q n+1 | 2 )δ m,n+1 +q n+1 ·(−i)(1− σ 2 |q n | 2 )δ m,n i + 4 σ X n (σ/2)q n 1− σ 2 |q n | 2 ·(−i)(1− σ 2 |q n | 2 )δ m,n =i h q m−1 (1− σ 2 |q m | 2 )+q m+1 (1− σ 2 |q m | 2 ) i −2iq m =i(q m−1 +q m+1 −2q m )−i σ 2 |q m | 2 (q m−1 +q m+1 ). This equation is completely integrable and has an infinite number of conserved quantities (for the infinite lattice) and can be solved by the method of Inverse Scatter- ing Transform (IST) developed by Gardner, Greene, Kruskal and Mira (Gardner et al., 1967). 4.2 Conserved Quantities of the AL equation In this section, we outline the approach of inverse scattering and the work of Ablowitz and Ladik for calculating the conserved quantities of the AL equation (Ablowitz and Ladik, 1976; Ablowitz and Segur, 1981). The method of inverse scattering is analogous to the Discrete Fourier Transform. In the direct scattering problem, scattering data is derived from initial data, the potential. The time evolutions of the scattering data is simple. From the scattering data at some time ’t’, the potential can be calculated via the (non-trivial) inverse scattering transform. 101 In order to calculate the conserved quantities is it only necessary to calculate the scat- tering data. From the scattering data it is clear which quantities are time independent. The following is an outline of Ablowitz and Ladik (1976) with details filled in using (Ablowitz and Segur, 1981). For the AL equation, we taking S n (t) = T n (t) = 0 in (Ablowitz and Ladik, 1976). The canonical pairs are R n and Q n . In the end we set R n =±αQ ∗ n . To begin, consider the generalized eigenvalue problem, V 1,n+1 =zV 1,n +Q n (t)V 2,n V 2,n+1 = 1 z V 2,n +R n (t)V 1,n , (4.11) which can be equivalently expressed as 0 @ ˆ E −Q n (t) −R n (t) ˆ E 1 A 0 @ V 1 V 2 1 A = 0 @ z 0 0 1/z 1 A 0 @ V 1 V 2 1 A (4.12) where ˆ E is the shift operator: ˆ EX n =X n+1 . In this form, one can see the analogy with the Schr¨ odinger equation in quantum mechanics withQ n andR n playing the role of the potential. The potentials correspond to the real-space classical fieldsψ n ,ψ ∗ n in the AL. The time-dependence is postulated to have the form ˙ V 1,n =A n (t)V 1,n +B n (t)V 2,n ˙ V 2,n =C n (t)V 1,n +D n (t)V 2,n . (4.13) 102 The associated equations of motion for A n ,B n ,C n ,D n are generated by cross- differentiating (4.11) and (4.13) ∂ ∂t (V 1,n+1 )=z ˙ V 1,n + ˙ Q n V 2,n +Q n ˙ V 2,n =z(A n V 1,n +B n V 2,n )+ ˙ Q n V 2,n +Q n (C n V 1,n +D n V 2,n ) ∂ ∂t (V 2,n+1 )= 1 z ˙ V 2,n + ˙ R n V 1,n +R n ˙ V 1,n = 1 z (C n V 1,n +D n V 2,n )+ ˙ R n V 1,n +R n (A n V 1,n +B n V 2,n ) µ ∂V 1,n 0 ∂t ¶ n 0 =n+1 =A n+1 V 1,n+1 +B n+1 V 2,n+1 =A n+1 (zV 1,n +Q n V 2,n )+B n+1 µ 1 z V 2,n +R n V 1,n ¶ µ ∂V 2,n 0 ∂t ¶ n 0 =n+1 =C n+1 V 1,n+1 +D n+1 V 2,n+1 =C n+1 (zV 1,n +Q n V 2,n )+D n+1 µ 1 z V 2,n +R n V 1,n ¶ (4.14) where the eigenvalue, z is time-invariant and the explicit time-dependence has been dropped. Requiring ∂ ∂t V i,n+1 = ¡ ∂ ∂t V i,n 0 ¢ n 0 =n+1 , z(A n V 1,n +B n V 2,n )+ ˙ Q n V 2,n +Q n (C n V 1,n +D n V 2,n ) = A n+1 (zV 1,n +Q n V 2,n )+B n+1 µ 1 z V 2,n +R n V 1,n ¶ 1 z (C n V 1,n +D n V 2,n )+ ˙ R n V 1,n +R n (A n V 1,n +B n V 2,n ) = C n+1 (zV 1,n +Q n V 2,n )+D n+1 µ 1 z V 2,n +R n V 1,n ¶ (4.15) 103 and equating the coefficients of V i,n for each equation results in the following time- evolution equations, zΔ n A n +R n B n+1 −Q n C n =0 (1/z)B n+1 −zB n +Q n (A n+1 −D n ) = ˙ Q n zC n+1 −(1/z)C n −R n (A n −D n+1 ) = ˙ R n (1/z)Δ n D n +Q n C n+1 −R n B n =0, (4.16) whereΔ n X n ≡X n+1 −X n . The linear dispersion relation and the above equations suggest the following expan- sions for the functions of (4.13): A n =z 2 A (2) n +A (0) n B n =zB (1) n + 1 z B (-1) n C n =zC (1) n + 1 z C (-1) n D n =D (0) n + 1 z 2 D (-2) n (4.17) 104 The coefficients of these expansions are determined by subsituting back into (4.16), z 3 Δ n A (2) n +zΔ n A (0) n +zR n B (1) n+1 +(1/z)R n B (-1) n+1 −zQ n C (1) n −(1/z)Q n C (-1) n =0 B (1) n+1 +(1/z 2 )B (-1) n+1 −z 2 B (1) n −B (-1) n +Q n (z 2 A (2) n+1 +A (0) n+1 −D (0) n −(1/z 2 )D (-2) n ) = ˙ Q n z 2 C (1) n+1 +C (-1) n+1 −C (1) n −(1/z 2 )C (-1) n −R n (z 2 A (2) n +A (0) n −D (0) n+1 −(1/z 2 )D (-2) n+1 ) = ˙ R n (1/z)Δ n D (0) n +(1/z 3 )Δ n D (-2) n +zQ n C (1) n+1 −(1/z)Q n C (-1) n+1 −zR n B (1) n −(1/z)R n B (-1) n =0, (4.18) and solving in powers of ’z’: O(z 3 ): Δ n A (2) n =0 ⇒ A (2) n+1 =A (2) n ≡A (2) O ¡ 1 z 3 ¢ : Δ n D (-2) n =0 ⇒ D (-2) n+1 =D (-2) n ≡D (-2) O(z 2 ): B (1) n −Q n A (2) n =0 ⇒ B (1) n =Q n A (2) C (1) n+1 −R n A (2) n =0 ⇒ C (1) n =R n−1 A (2) O ¡ 1 z 2 ¢ : B (-1) n+1 −Q n D (-2) n = ⇒ B (-1) n =Q n−1 D (-2) C (-1) n −R n D (-2) n =0 ⇒ C (1) n =R n D (-2) O(z): Δ n A (0) n =Q n C (1) n −R n B (1) n+1 =−A (2) (Q n+1 R n −Q n R n−1 ) ⇒ A (0) n =−A (2) Q n R n−1 +A (0) O ¡ 1 z ¢ : Δ n D (0) n =R n B (-1) n −Q n C (-1) n+1 =−D (-2) (Q n R n+1 −Q n−1 R n ) ⇒ D (0) n =−D (2) Q n−1 R n +D (0) . (4.19) 105 The zeroeth order yields the time-dependence of the potentials, O ¡ z 0 ¢ : ˙ Q n =B (1) n+1 −B (-1) 1 +Q n ³ A (0) n+1 −D (0) n ´ ˙ Q n = ¡ A (2) Q n+1 −D (-2) Q n−1 ¢ (1−Q n R n )+Q n ¡ A (0) −D (0) ¢ ˙ R n =C (-1) n+1 −C (1) n −R n ³ A (0) n −D (0) n+1 ´ ˙ R n = ¡ D (2) R n+1 −A (-2) R n−1 ¢ (1−Q n R n )−R n ¡ A (0) −D (0) ¢ . (4.20) In summary, the coefficients of the time-evolution equations are given by A n =A (2) (z 2 −Q n R n−1 )+A (0) B n =A (2) Q n z−D (-2) Q n−1 (1/z) C n =A (2) R n−1 z−D (-2) R n (1/z) D n =D (-2) ((1/z 2 )−Q n−1 R n )+D (0) , (4.21) while the time evolutions ofQ n andR n are governed by ˙ Q n = ¡ A (2) Q n+1 −D (-2) Q n−1 ¢ (1−Q n R n )+Q n ¡ A (0) −D (0) ¢ (4.22a) ˙ R n = ¡ D (2) R n+1 −A (-2) R n−1 ¢ (1−Q n R n )−R n ¡ A (0) −D (0) ¢ (4.22b) 4.2.1 AL equation To obtain the AL equation, we let R n = ±αQ ∗ n , A (2) = −D (-2) = i, and A (0) = −D (0) =−i. In this case equation (4.22a) becomes ˙ Q n =i £ Q n+1 +Q n−1 −2Q n ∓α|Q n | 2 (Q n+1 +Q n−1 ) ¤ , (4.23) 106 which is the AL equation. With these substitutions, the time-dependent function A n ,...,D n obey A n =i(z 2 ∓αQ n Q ∗ n−1 −1) B n =i(zQ n −(1/z)Q n−1 ) C n =±iα(zQ ∗ n−1 −(1/z)Q ∗ n ) D n =i(1−(1/z 2 )±αQ n−1 Q ∗ n ). (4.24) 4.2.2 Direct Scattering Problem Asymptotic Solution of Generalized Eigenvalue Problem First of all we assume that Q n and R n are on compact support, that is that as|n| → ∞, Q n ,R n → 0. In the limit|n|→∞, the generalized eigenvalue problem becomes, |n|→∞: 8 < : V 1,n+1 = zV 1,n V 2,n+1 = 1 z V 2,n . (4.25) To solve this direct scattering problem, define the time-independent eigenfunctions, φ n , ¯ φ n ,ψ n , ¯ ψ n which have the asymptotic forms: n→−∞: φ n ∼ 0 @ 1 0 1 A z n , ¯ φ n ∼ 0 @ 0 −1 1 A z −n n→+∞: ψ n ∼ 0 @ 0 1 1 A z −n , ¯ ψ n ∼ 0 @ 1 0 1 A z n . (4.26) 107 These eigenfunctions satisfy the generalized eigenvalue problem at a fixed time, which will be taken to be t=0. They do not solve the time-evolution equations, which we will return to later. One can establish by induction, that z −n φ n , z n ψ n are polynomial in powers of 1 z and are analytic for|z|>1 z n ¯ φ n , z −n ¯ ψ n are polynomial in powers ofz and are analytic for|z|<1 (|z|<∞) (4.27) given thatQ n andR n are on compact support. Wronskian and Linear Independence The Wronskian of a set of functions {φ 1 ,...φ n } is defined by: W(φ 1 ,...φ n )≡ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ φ 1 φ 2 ··· φ n φ 0 1 φ 0 2 ··· φ 0 n . . . . . . . . . . . . φ (n−1) 1 φ (n−1) 2 ··· φ (n−1) n ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ (4.28) If the Wronskian is non-zero in some interval, then the functions are linearly indepen- dent in that interval. 108 For this problem, the Wronskian is defined asW n (u,v) = u 1,n v 2,n −u 2,n v 1,n . The Wronskian of the functions that obey (4.11) can be found by W n+1 (u,v)=u 1,n+1 v 2,n+1 −u 2,n+1 v 1,n+1 = (zu 1,n +Q n u 2,n )((1/z)v 2,n +R n v 1,n ) −((1/z)u 2,n +R n u 1,n )(zv 1,n +Q n v 2,n ) =u 1,n v 2,n +zR n u 1,n v 1,n +(1/z)Q n u 2,n v 2,n +Q n R n u 2,n v 1,n −(u 2,n v 2,n +(1/z)Q n u 2,n v 2,n +zR n u 1,n v 1,n +Q n R n u 2,n v 1,n ) = (1−R n Q n )(u 1,n v 2,n −u 2,n v 1,n ) = (1−R n Q n )W n (u,v) (4.29) Using the asymptotic forms ofφ n , ¯ φ n ,ψ n , ¯ ψ n , one gets lim n→−∞ W n ( ¯ φ,φ)=− ¯ φ 2,n φ 1,n =−(-z −n )(z n )=1 lim n→+∞ W n ( ¯ ψ,ψ)=− ¯ ψ 1,n ψ 2,n =(z n )(z −n )=1. (4.30) By induction, one can show that the Wronskian of the eigenfunctions obey W n ( ¯ φ,φ)= n Y i=−∞ (1−R i Q i ), W n ( ¯ ψ,ψ)= ∞ Y i=n+1 (1−R i Q i ) −1 (4.31) forz on the unit circle. ForR n =−Q ∗ n ,W n is positive-definite so that ¯ φ n ,φ n are linearly independent. Otherwise assume that initiallyR n andQ n are less than one. Likewise one can show that ¯ ψ n ,ψ n are linearly independent. Thus one can define the scattering data by φ n =a(z,t) ¯ ψ n +b(z,t)ψ n ¯ φ n = −¯ a(z,t)ψ n + ¯ b(z,t) ¯ ψ n (4.32) 109 wheret is a parameter. Time Dependence The eigenfunctions φ, ¯ φ,ψ, ¯ ψ satisfy the generalized eigenvalue problem, but not the time evolution equations (4.13). In the asymptotic limit (4.13) becomes ˙ V 1,n = A ± V 1,n ˙ V 2,n = D ± V 2,n (4.33) where A ± = lim n→±∞ (A n )=z 2 A (2) +A (0) D ± = lim n→±∞ (D n )=(1/z 2 )D (2) +D (0) . (4.34) Define a new set of eigenfunctions which will satisfy both the generalized eigenvalue problem and the time-dependence φ (t) n =φ n exp(A t) ψ (t) n =ψ n exp(D + t) ¯ φ (t) n = ¯ φ n exp(D t) ψ (t) n =ψ n exp(A + t) (4.35) Both{ ¯ φ (t) n ,φ (t) n } and{ ¯ ψ (t) n ,ψ (t) n } are linearly independent, so we may write φ (t) n =a 0 ¯ ψ (t) n +b 0 ψ (t) n ¯ φ (t) n = −¯ a 0 ψ (t) n + ¯ b 0 ¯ ψ (t) n (4.36) wherea 0 ,¯ a 0 ,b 0 , ¯ b 0 equal the scattering coefficients,a,¯ a,b, ¯ b at t=0. Substituting (4.35) into (4.36) gives φ n =a 0 e [(A + −A − )t] ¯ ψ+b 0 e [(D + −A − )t] ψ n ¯ φ n = −¯ a 0 e [(D + −D − )t] ψ n + ¯ b 0 e [(A + −D − )t] ¯ ψ n (4.37) 110 Comparing with (4.32) and noting thatA + = A − andD + = D − , the scattering coeffi- cients are a=a 0 , b=b 0 e [(D + −A − )t] ¯ a= ¯ a 0 , ¯ b= ¯ b 0 e [(D + −A − )t] . (4.38) The most significant result here, for our purposes, is that the coefficients a,¯ a are con- stant. In the next section, we use this to derive the conserved quantities of the AL equation. 4.2.3 Conservation Laws The conservation laws can be derived by considering the asymptotic form of¯ a(z). ¯ φ n ∼−¯ a(z) 0 @ 0 1 1 A z −n + ¯ b(z) 0 @ 1 0 1 A z n ⇒ ¯ a(z)∼−z n ¯ φ 2,n (4.39) Next substitute ¯ φ n into the eigenvalue problem (4.11), z n−1 ¯ φ 1,n+1 =z n ¯ φ 1,n +z n−1 Q n ¯ φ 2,n (4.40a) z n ¯ φ 2,n+1 =z n−1 ¯ φ 2,n +z n R n ¯ φ 1,n . (4.40b) Solve (4.40b) forz n ¯ φ 1,n and substitute into (4.40a), to eliminate ¯ φ 1,n z n ¯ φ 1,n = 1 R n ¡ z n ¯ φ 2,n+1 −z n−1 ¯ φ 2,n ¢ = 1 zR n Δ n ¡ z n ¯ φ 2,n ¢ z n−2 1 R n+1 Δ n ¡ z n+1 ¯ φ 2,n+1 ¢ = 1 zR n Δ n ¡ z n ¯ φ 2,n ¢ +z n−1 Q n ¯ φ 2,n 1 z 2n+3 R n+1 Δ n ¡ z n+1 ¯ φ 2,n+1 ¢ − 1 z 2n+1 R n Δ n ¡ z n ¯ φ 2,n ¢ =Q n z −2n−1 ¡ z n ¯ φ 2,n ¢ , (4.41) 111 which can be written as, Δ n µ Δ n (z n ¯ φ 2,n ) z 2n+1 R n ¶ =Q n z −2n−1 (z n ¯ φ 2,n ). (4.42) Define −z n ¯ φ 2,n = n Y k=−∞ g k (4.43) such that Δ n (−z n ¯ φ 2,n ) = Δ n ( Q n k=−∞ g k ) = (g n+1 −1) Q n k=−∞ g k and substitute into (4.42) to get the recursion relation g n+1 (g n+2 −1)−z 2 R n+1 R n (g n+1 −1)=z 2 R n+1 Q n , (4.44) which can be re-written as R n−1 g n (g n+1 −1)−z 2 R n (g n −1)=z 2 R n R n−1 Q n−1 . (4.45) Expandg n andg n+1 in powers ofz 2 , g n =g (0) n +z 2 g (1) n +z 4 g (2) n +...= X m=0 g (m) n z 2m g n+1 =g (0) n+1 +z 2 g (1) n+1 +z 4 g (2) n+1 +...= X m=0 g (m) n+1 z 2m . (4.46) 112 Substitute into the recursion relation forg n ,g n+1 , R n−1 X m=0 g (m) n z 2m ( X l=0 g (l) n+1 z 2l −1)−z 2 R n ( X m=0 g (m) n z 2m −1) =z 2 R n R n−1 Q n−1 R n−1 " X m,l=0 z 2(m+l) g (m) n g (l) n+1 − X m=0 z 2m g (m) n # −R n X m=0 z 2(m+1) g (m) n =z 2 R n (R n−1 Q n−1 −1). (4.47) Solving recursively forg (m) n in orders ofz 2 , gives O(z 0 ): g (0) n =1 O(z 2 ): g (1) n =R n−1 Q n−2 O(z 4 ): g (2) n =R n−1 Q n−3 (1−R n−2 Q n−2 ) (4.48) For orders ofz 2p ,p≥2, this recursion relation reduces to R n−1 " p X m=0 g (m) n g (p−m) n+1 −g (p) n # =R n−1 " p−1 X m=1 g (m) n g (p−m) n+1 +g (p) n+1 # =R n g (p−1) n (4.49) or g (p) n+1 = R n R n−1 g (p−1) n − p−1 X m=1 g (m) n g (p−m) n+1 . (4.50) 113 which gives expressions for higher orders of the coefficients, O(z 6 ): g (3) n = R n−1 R n−2 g (2) n−1 −g (1) n−1 g (2) n −g (2) n−1 g (1) n =R n−1 Q n−4 ¡ 1−R n−3 Q n−3 −R n−2 Q n−2 +R n−3 Q n−3 R n−2 Q n−2 ¢ −R n−1 R n−2 Q 2 n−3 (1−R n−2 Q n−2 ) O(z 8 ): g (4) n = R n−1 R n−2 g (3) n−1 −g (1) n−1 g (3) n −g (2) n−1 g (2) n −g (3) n−1 g (1) n O(z 10 ): g (5) n = R n−1 R n−2 g (4) n−1 −g (1) n−1 g (4) n −g (2) n−1 g (3) n −g (3) n−1 g (2) n −g (4) n−1 g (1) n (4.51) Next use the expansion log(1+x)= ∞ X n=1 x n n (4.52) to write out the series expansions, log[¯ a(z)]= lim n→∞ log Ã n Y k=−∞ g k ! = ∞ X n=−∞ log(g n )= ∞ X n=−∞ log " ∞ X m=0 g (m) n z 2m # = ∞ X n=−∞ log " 1+ ∞ X m=1 g (m) n z 2m # = ∞ X n=−∞ ∞ X p=1 1 p " ∞ X m=1 g (m) n z 2m # p . (4.53) 114 Since ¯ a(z) is a constant of motion, each coefficients of z 2α in the above expansion of the must also be time-independent. These coefficients are the conserved quantities. Expanding the first few terms gives, log[¯ a(z)]= ∞ X n=−∞ ∞ X p=1 1 p £ g (1) n z 2 +g (2) n z 4 +g (3) n z 6 +g (4) n z 8 +g (5) n z 10 +O(z 12 ) ¤ p = ∞ X n=−∞ g (1) n z 2 + µ g (2) n + 1 2 £ g (1) n ¤ 2 ¶ z 4 + µ g (3) n +g (1) n g (2) n + 1 3 £ g (1) n ¤ 3 ¶ z 6 + µ g (4) n + 1 2 £ g (2) n ¤ 2 +g (1) n g (3) n + £ g (1) n ¤ 2 g (2) n + 1 4 £ g (1) n ¤ 4 ¶ z 8 + µ g (5) n +g (4) n g (1) n +g (3) n g (2) n + £ g (2) n ¤ 2 g (1) n +g (3) n £ g (1) n ¤ 2 +g (2) n £ g (1) n ¤ 3 ¶ z 10 +O(z 12 ). (4.54) 115 4.2.4 Conserved Quantities The first few conserved quantities are: C 1 = X n g (1) n = X n R n Q n−1 C 2 = X n g (2) n + 1 2 £ g (1) n ¤ 2 = X n R n Q n−2 (1−R n−1 Q n−1 )+ 1 2 R 2 n Q 2 n−1 C 3 = X n g (3) n +g (1) n g (2) n + 1 3 £ g (1) n ¤ 3 = X n R n Q n−3 (1−R n−2 Q n−2 −R n−1 Q n−1 +R n−2 Q n−2 R n−1 Q n−1 ) −R n R n−1 Q 2 n−2 (1−R n−1 Q n−1 )+R 2 n Q n−1 Q n−2 (1−R n−1 Q n−1 ) + 1 3 (R n Q n−1 ) 3 C 4 = X n g (4) n + 1 2 £ g (2) n ¤ 2 +g (1) n g (3) n + £ g (1) n ¤ 2 g (2) n + 1 4 £ g (1) n ¤ 4 C 5 = X n g (5) n +g (4) n g (1) n +g (3) n g (2) n + £ g (2) n ¤ 2 g (1) n +g (3) n £ g (1) n ¤ 2 +g (2) n £ g (1) n ¤ 3 (4.55) Further conserved quantities can be determined by using the recursion relations and expansions given. From looking at these equations, we notice a few patterns. First of all, the leading order inα ofC m is P n ψ ∗ n ψ n−m , and the largest order ofα isα m . 4.3 AL and BHM To obtain the AL, 4.23 we let R n = αQ ∗ n . Furthermore, to write it in a more familiar form, we write in terms of the real-space fields, ˜ ψ n =Q n , and ˜ ψ ∗ n =Q ∗ n . ∂ ∂t ˜ ψ n =i h³ ˜ ψ n+1 + ˜ ψ n−1 −2 ˜ ψ n ´ −α| ˜ ψ n | 2 ( ˜ ψ n+1 + ˜ ψ n−1 ) i (4.56) 116 Recall that in real-space the BHM has the form ∂ ∂t ˜ ψ n =−i h −J ³ ˜ ψ n+1 + ˜ ψ n−1 −2 ˜ ψ n ´ +μ 0 N s | ˜ ψ n | 2 ˜ ψ i (4.57) By defining a new time scale,τ =Jt, ∂ ∂τ ˜ ψ n =i · ³ ˜ ψ n+1 + ˜ ψ n−1 −2 ˜ ψ n ´ − μ 0 N s J | ˜ ψ n | 2 ˜ ψ n ¸ (4.58) Note that here n is an index of real-space fields, following the notation of Ablowitz and Ladik (1976), while in earlier chapters it was used as a momentum index. From these two forms, it can be seen that these two equations are mathematically close for α =2μ 0 N s /J ≡2κN s and identical in the limitψ n+1 +ψ n−1 →2ψ n . The first few coefficients ofg n are g (1) n =α ˜ ψ ∗ n−1 ˜ ψ n−2 g (2) n =α ˜ ψ ∗ n−1 ˜ ψ n−3 ³ 1−α| ˜ ψ n−2 | 2 ´ g (3) n = α ˜ ψ ∗ n−1 α ˜ ψ ∗ n−2 g (2) n−1 −g (1) n−1 g (2) n −g (2) n−1 g (1) n =α ˜ ψ ∗ n−1 ˜ ψ n−4 ³ 1−α| ˜ ψ n−3 | 2 −α| ˜ ψ n−2 | 2 +α 2 | ˜ ψ n−3 | 2 | ˜ ψ n−2 | 2 ´ −α 2 ˜ ψ ∗ n−1 ˜ ψ ∗ n−2 ˜ ψ 2 n−3 ³ 1−α| ˜ ψ n−2 | 2 ´ . (4.59) 117 In terms of the real-space fields{ e ψ n , e ψ ∗ n }, the first few conserved quantities are C 1 = X n α ˜ ψ ∗ n ˜ ψ n−1 C 2 = X n g (2) n + 1 2 £ g (1) n ¤ 2 = X n α ˜ ψ ∗ n ˜ ψ n−2 ³ 1−α| ˜ ψ n−1 | 2 ´ + 1 2 α 2 ( ˜ ψ ∗ n ) 2 ˜ ψ 2 n−1 C 3 = X n g (3) n +g (1) n g (2) n + 1 3 £ g (1) n ¤ 3 = X n α ˜ ψ ∗ n ˜ ψ n−3 ³ 1−α| ˜ ψ n−2 | 2 −α| ˜ ψ n−1 | 2 +α 2 | ˜ ψ n−2 | 2 | ˜ ψ n−1 | 2 ´ −α 2 ˜ ψ ∗ n ˜ ψ ∗ n−1 ˜ ψ 2 n−2 ³ 1−α| ˜ ψ n−1 | 2 ´ +α[ ˜ ψ ∗ n ] 2 ˜ ψ n−1 ˜ ψ n−2 ³ 1−α| ˜ ψ n−1 | 2 ´ + 1 3 ³ α ˜ ψ ∗ n ˜ ψ n−1 ´ 3 . (4.60) In terms of the momentum-space fields{ψ p ,ψ ∗ p },C 1 andC 2 are C 1 =α X p |ψ p | 2 e −2πip/N s =α X p |ψ p | 2 cos µ 2πp N s ¶ −iα X p |ψ p | 2 sin µ 2πp N s ¶ C 2 =α X p |ψ p | 2 e −4πip/N s − α 2 N s X p,q,r ψ ∗ p ψ ∗ q ψ r ψ p+q−r e −2πi(p+r)/N s + α 2 2N s X p,q,r ψ ∗ p ψ ∗ q ψ r ψ p+q−r e −2πi(p+q)/Ns . (4.61) As seen in earlier chapters, there is a threshold for chaos in the BHM. Below this threshold there are initial states that do not thermalize, but nevertheless relax to a steady- state. For fully integrable systems it is known that this steady-state can be described by a constrained thermodynamic ensemble that accounts for all of the integrals of motion. Anecdotal evidence suggests that the initial total quasi-momentum plays a role in the 118 dynamics of the BHM for small nonlinearities, κ. Is it possible that the other integrals of motion of the AL affect the dynamics? Is the nonthermal steady-state governed by a constained ensemble that takes into account the conserved quantities of nearby inte- grable models? There are other nearby integrable models, which include the noninter- acting case and the continuum limit. The imaginary part of the first conserved quantity of AL, ImC 1 , is analogous to the total quasi-momentum and ReC 1 to the total kinetic energy. In the noninteracting the momentum distribution is conserved, so any higher moments of the momentum distribution is also conserved. If the mapping between the BHM and AL fields corresponded to simply equating the fields (ψ n,BHM = ψ n,AL ), then conservation of C 1 in BHM would not distinguish between being close to the nonin- teracting case and close to AL. In contrast, the second and third terms of C 2 are not conserved in the noninteracing model. However, the mapping between the BHM fields and the AL is nontrivial and the nonlinear corrections to the mapping are expected to allow one to distinguish between effects of quasi-conserved quantities of AL and of the noninteracting case. A derivation of the mapping between the BHM and AL fields and a study of the role of the conserved quantities of AL in the dynamics of BHM is the subject of proposed future work. 119 Chapter 5 Outlook and Conlusion 5.1 Summary of Results One of the fundamental assertions of statistical mechanics is that the time average of a physical observable is equivalent to the average over phase-space, with microcanonical measure. A system for which this is true is said to be ergodic and one can calculate dynamical properties of the system from static phase-space averages. Dynamics of a system which is fully integrable, that is has as many conserved quantities as degrees of freedom, is constrained to a reduced phase space and thus not ergodic, although it may relax to a modified equilibrium. What happens as one moves away from the fully integrable case? In this work we have studied the relationship between chaos and thermalization in the one-dimensional Bose-Hubbard model in the classical-field approximation. We have compared two quantitative measures of chaos and thermalization: (1) the finite-time maximal Lyapunov exponent, averaged over a microcanonical ensemble and (2) the normalized spectral entropy, which is a measure of equipartition of a modified energy in the independent mode approximation. There is a strong correspondence between the Lyapunov exponents and normalized spectral entropy. We find a threshold for chaos and a corresponding broad transition from incomplete to complete thermalization. The stochasticity threshold is governed by two parameters: the strength of the nonlinearityκ and the average energy-per-particle,² T . Both of these parameters are finite in the thermodynamic limit and suggest that the threshold will 120 survive in that limit. Far above the threshold, in the strongly chaotic regime, relaxation to the thermal state is complete. In this region, the fluctuations in kinetic energy scale as √ N s confirming their thermal nature. We study the size scaling of the Lyapunov exponent and find that it is universal with respect to the size of the lattice. For small nonlinearities, the stochasticity and thermalization thresholds are finite for the range of energies studied and don’t show tendencies to vanish at high energies. In the vicinity of the threshold the relationship between chaos and thermalization is complex. There is a transient regime supporing both chaotic and regular trajecto- ries, so that the Lyapunov exponent is non-zero, but full thermalization does not occur. Remarkably, in this region individual initial states with larger Lyapunov exponent tend to relax closer to the thermal state. There is also a region where although the Lyapunov exponent is zero, there is significant relaxation towards the thermal state. We conjecture that this redistribution in the phase-space is due to the quasi-regular dynamics governed by the nearby Ablowitz-Ladik lattice, which is integrable. Above ² T ' 0.6J, both the stochasticity threshold and thermalization threshold overlap very closely and appear to depend only on the nonlinearity strength. This suggests that there is a critical nonlinear- ity,κ c '0.2, such that the system is regular forκ<κ c independent of the energy of the system. In the opposite limit of small energy and large nonlinearity, there is a separation of thresholds, where almost complete relaxation to the thermal distribution is observed in the absence of chaos. An analysis of resonances of the BHM gives a Chirikov’s criterion for the chaos threshold that depends on parameters which vanish in the thermodynamic limit. This conclusion is confirmed on dimensional grounds. The criterion predicted by the Chrikov criterion is different from the one inferred from numerical calculations, signifying the failure of the standard Chirikov’s approach. 121 Anecdotal evidence suggests that the total quasi-momentum may play a role in the relaxation dynamics. The quasi-momentum is strictly conserved in the nearby fully inte- grable Ablowitz-Ladik model, as well as in the non-interacting and continuum limits. There are at least three known near-by integrable models: the Ablowitz-Ladik lat- tice, the continuous nonlinear Schr¨ odinger equation and the noninteracting model. We outline the method of Inverse Scattering Transform and generate all of the integrals of motion of the closely related, fully integrable model of Ablowitz-Ladik. We conjec- ture that the presense of quasi-conserved quantities may alter the scaling of the chaos criterion. 5.2 Open Questions These observations lead to many questions that deserve further investigation. • What is the reason for the failure of the Chirikov criterion to accurately predict the chaos threshold? Is it related to interference between resonances due to near-by integrable systems? Could the proper scaling be recovered in a multiple-resonance model? • What is the underlying theory that governs the threshold? • Forκ&1, how does the chaos threshold scale in the thermodynamic limit? Is the number of modes involved relevant? • What governs the slow relaxation times where they appear? Is it related to the conserved quantities of the near-by integrable systems? If this is the case, which near-by integrable system? 122 • For those states that show no signs of thermalization, what governs the steady- state? Can these states be described by a constrained ensemble? Which “quasi- conserved quantities”are the relevant for the constrained ensemble? 123 Bibliography Ablowitz, M., and H. Segur, 1981, Solitons and the Inverse Scattering Transform (Philadelphia, PA: SIAM). Ablowitz, M. J., and J. F. Ladik, 1976, Journal of Mathematical Physics 17(6), 1011. Ablowitz, M. J., C. Schober, and B. M. Herbst, 1993, Physical Review Letters 71(17), 2683. Anderson, M. H., J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, 1995, Science 269(5221), 198. Andrews, M. R., C. G. Townsend, H.-J. Miesner, D. S. Durfee, D. M. Kurn, and W. Ket- terle, 1997, Science 275(5300), 637. Arnold, V ., 1963, Russian Math. Survey 18, 85. Arnold, V ., and A. Avez, 1968, Ergodic Problems of Classical Mechanics (Benjamin). Aspect, A., E. Arimondo, R. Kaiser, N. Vansteenkiste, and C. Cohen-Tannoudji, 1988, Physical Review Letters 61(7), 826. Bose, S., 1924, Z. Phys. 26, 178. Bouchoule, I., K. V . Kheruntsyan, and G. V . Shlyapnikov, 2007, Physical Review A 75(3), 031606. Boyanovsky, D., C. Destri, and H. J. de Vega, 2004, Physical Review D 69(4), 045003. Bradley, C. C., C. A. Sackett, J. J. Tollett, and R. G. Hulet, 1995, Physical Review Letters 75(9), 1687. Castin, Y ., 2001, in Coherent Atomic Matter Waves, edited by R. Kaiser, C. Westbrook, and F. David (Springer-Verlag). Castin, Y ., 2004, cond-mat/0407118 . Chirikov, B. V ., 1979, Physics Reports 52(5), 263. 124 Chu, S., L. Hollberg, J. E. Bjorkholm, A. Cable, and A. Ashkin, 1985, Physical Review Letters 55(1), 48. Contopoulos, G., L. Galgani, and A. Giorgilli, 1978, Physical Review A 18(3), 1183. Contopoulos, G., and N. V oglis, 1997, Astronomy and Astrophysics 317(1), 73. Davis, K. B., M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, 1995, Physical Review Letters 75(22), 3969. Durfee, D. S., Y . K. Shaham, and M. A. Kasevich, 2006, Physical Review Letters 97(24), 240801. Eckhardt, B., and D. Yao, 1993, Physica D: Nonlinear Phenomena 65(1-2), 100. Einstein, A., 1925, Sitzber. Kgl. Preuss. Akad. Wiss. 3. Esteve, J., J.-B. Trebbia, T. Schumm, A. Aspect, C. I. Westbrook, and I. Bouchoule, 2006, Physical Review Letters 96(13), 130403. Fermi, E., J. Pasta, and S. Ulam, 1974, Appl. Math 15, 143. Folman, R., P. Kr¨ oger, J. Schmiedmayer, J. Denschlag, and C. Henkel, 2002, Advances In Atomic, Molecular, and Optical Physics 48, 263. Ford, J., 1992, Physics Reports 213(5), 271. Fortagh, J., and C. Zimmermann, 2007, Reviews of Modern Physics 79(1), 235. Gardner, C. S., J. M. Greene, M. D. Kruskal, and R. M. Miura, 1967, Physical Review Letters 19(19), 1095. Gerhardt, G. J. L., M. Frichembruder, F. B. Rizzato, and S. R. Lopes, 2002, Chaos, Solitons & Fractals 13(6), 1269. Girardeau, M., 1960, Journal of Mathematical Physics 1(6), 516. Greiner, M., O. Mandel, T. Esslinger, T. W. Hansch, and I. Bloch, 2002, Nature 415(6867), 39. Gross, E., 1961, Il Nuovo Cimento (1955-1965) 20(3), 454. Gustavson, T. L., A. Landragin, and M. A. Kasevich, 2000, Classical and Quantum Gravity 17(12), 2385. Hamer, C. J., and J. B. Kogut, 1979, Physical Review B 20(9), 3859. Hansel, W., P. Hommelhoff, T. W. Hansch, and J. Reichel, 2001, Nature 413(6855), 498. 125 Herbst, B. M., and M. J. Ablowitz, 1989, Physical Review Letters 62(18), 2065. Herbst, B. M., F. Varadi, and M. J. Ablowitz, 1994, Math. Comput. Simul. 37(4-5), 353. Izrailev, F., and B. Chirikov, 1966, in Dokl. Akad. Nauk SSSR, volume 166, p. 57. Jaksch, D., C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, 1998, Physical Review Letters 81(15), 3108. Jaynes, E. T., 1957a, Physical Review 106(4), 620. Jaynes, E. T., 1957b, Physical Review 108(2), 171. Kagan, Y ., and B. V . Svistunov, 1997, Physical Review Letters 79(18), 3331. Kinoshita, T., T. Wenger, and D. S. Weiss, 2006, Nature 440(7086), 900. Kolmogorov, A., 1954, Dokl. Akad. Nayk SSSR V ol 48, 527. Lebowitz, J. L., 1999, Physica A 263(1-4), 516. Lieb, E. H., 1963, Physical Review 130(4), 1616. Lieb, E. H., and W. Liniger, 1963, Physical Review 130(4), 1605. Livi, R., M. Pettini, S. Ruffo, M. Sparpaglione, and A. Vulpiani, 1985, Physical Review A 31(2), 1039. Livi, R., M. Pettini, S. Ruffo, and A. Vulpiani, 1987, Journal of Statistical Physics 48(3), 539. Mishmash, R. V ., and L. D. Carr, 2008, 0810.2593 . Moser, J., 1962, Nachr. Akad. Wiss. Goett. II, Math.-Phys. Kl 2, 1. Ohberg, P., and S. Stenholm, 1997, Journal of Physics B: Atomic, Molecular and Optical Physics 30(12), 2749. Ott, H., J. Fortagh, G. Schlotterbeck, A. Grossmann, and C. Zimmermann, 2001, Phys- ical Review Letters 87(23), 230401. Paredes, B., A. Widera, V . Murg, O. Mandel, S. F¨ olling, I. Cirac, G. V . Shlyapnikov, T. W. H¨ ansch, and I. Bloch, 2004, Nature 429(6989), 277. Penrose, O., and L. Onsager, 1956, Physical Review 104(3), 576. Phillips, W. D., and H. Metcalf, 1982, Physical Review Letters 48(9), 596. Pitaevskii, L., 1961, Soviet Physics JETP-USSR 13(2), 451. 126 Proukakis, N. P., J. Schmiedmayer, and H. T. C. Stoof, 2006, Physical Review A 73(5), 053603. Rigol, M., V . Dunjko, V . Yurovsky, and M. Olshanii, 2007, Physical Review Letters 98(5), 050405. Scharf, R., and A. R. Bishop, 1991, Physical Review A 43(12), 6535. Schumm, T., S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer, and P. Kruger, 2005, Nat. Phys. 1(1), 57. Tabor, M., 1989, Chaos and integrability in nonlinear dynamics: an introduction (Wiley). Tuck, J., and M. Menzel, 1972, Adv. Math 9, 399. Vidal, G., 2004, Physical Review Letters 93(4), 040502. Villain, P., and M. Lewenstein, 2000, Physical Review A 62(4), 043601. V oglis, N., and G. J. Contopoulos, 1994, Journal of Physics A: Mathematical and Gen- eral 27(14), 4899. Wang, Y .-J., D. Z. Anderson, V . M. Bright, E. A. Cornell, Q. Diot, T. Kishimoto, M. Prentiss, R. A. Saravanan, S. R. Segal, and S. Wu, 2005, Physical Review Let- ters 94(9), 090405. Weiss, D. S., B. C. Young, and S. Chu, 1994, Applied Physics B: Lasers and Optics 59(3), 217. Wicht, A., J. M. Hensley, E. Sarajlic, and S. Chu, 2002, Physica Scripta T102, 82. Yukalov, V ., 2009, Laser Physics 19(1), 1. Zabusky, N. J., and M. D. Kuskal, 1965, Physical Review Letters 15(6), 240. Zaslavsky, G. M., 1999, Physics Today 52(8), 39. 127 Appendix A Thermodynamic Distribution within Hartree-Fock A.1 Hartree Fock In this section, the following integrals are used: Z ∞ 0 dxe −αx = 1 α Z ∞ 0 dxxe −αx = 1 α 2 Z ∞ 0 dxx 2 e −αx = 2 α 3 (A.1) Note thatN is the number of degrees of freedom anda=L/N is the lattice spacing. Within Hartree-Fock the form of the density distribution function is taken to be Gaussian and the thermal expectation value of the Grand Potential, hFi=hHi−ThSi−μhN a i, (A.2) is minimized, where N a is the norm. The density distribution function with two-body interactions has the form, σ HF = 1 Z exp Ã X n,n 0 −α n,n 0ψ n ψ ∗ n 0 ! = 1 Z exp Ã − X n α n I n ! . (A.3) 128 In the second step, we assume that off-diagonal terms are zero (justify). Theα n coeffi- cients are unknown and are determined by the condition of minimizing the grand poten- tial. Calculate Z, so that the integration ofσ HF over all of phase space is 1. Z = 1 (2π~) N Z d N θ Z d N Ie − P n α n I n = 1 (2π~) N Z 2π 0 dθ 1 Z 2π 0 dθ 2 ··· Z 2π 0 dθ N Z ∞ 0 dI 1 Z ∞ 0 dI 2 ··· Z ∞ 0 dI N e − P n α n I n = 1 ~ N Z ∞ 0 dI 1 e −α 1 I 1 Z ∞ 0 dI 2 e −α 2 I 2 ··· Z ∞ 0 dI N e −α N I N = N Y i=1 1 ~α i (A.4) The expectation values of each term in the Grand Potential is calculated, using the Hartree-Fock density distribution function, σ HF . The expectation value of a generic observable is given by hOi= 1 (2π~) N Z d N θ Z d N IO({I n ,θ n })σ HF = 1 (2π) N N Y i=1 α i Z d N θ Z d N IOe − P n α n I n (A.5) The expectation values of the relevant observable are calculated below. 129 Norm hNormi= 1 (2π~) N N Y i=1 ~α i Z d N θ Z d N I " 1 ~ X n I n # e − P m αmIm = N Y i=1 α i X n 1 ~ Z ∞ 0 dI 1 e −α 1 I 1 ··· Z ∞ 0 dI n I n e −αnIn ··· Z ∞ 0 dI N e −α N I N = N X n=1 1 ~α n Hamiltonian - Kinetic Term hH 0 i= 1 (2π~) N N Y i=1 ~α i Z d N θ Z d N I " X n I n ω n # e − P m α m I m = N X n=1 ω n α n Hamiltonian - Interaction Term hH I i= 1 (2π~) N N Y i=1 ~α i Z d N θ Z d N I " μ 0 2~ 2 X m,p,q,r (I m I p I q I r ) 1/2 δ m+p,q+r e −i(θm+θp−θq−θr) # e − P n αnIn = 1 (2π) N N Y i=1 α i (2π) N Z d N I μ 0 2~ 2 Ã X m I 2 m +2 X m6=p I m I p ! e (− P n αnIn) = μ 0 2~ 2 N Y i=1 α i Ã N Y j=1 1 α j X m 2 α 2 m +2 X m6=p 1 α m 1 α p ! = μ 0 ~ 2 Ã X m 1 α 2 m + X m6=p 1 α m 1 α p ! = μ 0 ~ 2 X m,p 1 α m 1 α p (A.6) 130 Entropy S =− 1 (2π~) N Z d N θ Z d N σ HF logσ HF = 1 (2π~) N Z d N θ Z d N 1 Z e − P n α n I n Ã − X n α n I n +logZ ! = 1 ~ N N Y i=1 (~α i ) Z d N e − P n α n I n Ã − X n α n I n +log N Y j=1 1 ~α j ! = X n α n 1 α n − X j log(~α j ) =N− X j log(~α j ) (A.7) A.2 Minimization of the Grand Potential The thermal expectation value of the Grand Potential within Hartree-Fock is given by hFi= X m ω m α m + μ 0 ~ 2 X m,p 1 α m 1 α p −T Ã N− X m log(~α m ) ! −μ X m 1 ~α m (A.8) Taking the variation with respect toα n , and setting it equal to zero gives δhFi δα n =− ω n α 2 n −2 μ 0 ~ 2 1 α 2 n X m 1 α m + T α n + μ ~α 2 n =0 (A.9) Using P m α −1 m =~N a and solving forα n , α n = 1 ~T [~ω n +2μ 0 N a −μ] (A.10) The thermal expectation values of the occupation of momentum moden become hI n i= 1 α n = ~T ~ω n +2μ 0 N a −μ . (A.11) 131 In general, the coefficientsμ andT are unknown and are determined by imposing con- straints on the norm and energy, which come from the dynamical code. The constraints are N a =hN a i= 1 ~ X n hI n i E T =hHi= X n ω n hI n i+ μ 0 ~ 2 X m,n hI m ihI n i= X n ω n hI n i+μ 0 N 2 a (A.12) Beginning with the expression forhI n i, we can solve forT in terms ofμ,N a and energy. T =ω n hI n i+2 μ 0 ~ N a hI n i− μ ~ hI n i (A.13) Summing overn T = 1 N £ E k +2μ 0 N 2 a −μN a ¤ (A.14) where E k ≡ P n ω n hI n i. This expression for T can be substituted back into the con- straints to reduce the system to two equations with two unknowns. Using the expression for temperature and normalization condition, a single constraint remains to be solved, 1 N X n [E k +2μ 0 N 2 a −μN a ] [~ω n +2μ 0 N a −μ] −N a =0 (A.15) 132

Abstract (if available)

Abstract One of the fundamental assertions of statistical mechanics is that the time average of a physical observable is equivalent to the average over phase space, with microcanonical measure. A system for which this is true is said to be ergodic and dynamical properties can be calculated from static phase-space averages. Dynamics of a system which is fully integrable, that is has as many conserved quantities as degrees of freedom, is constrained to a reduced phase space and thus not ergodic, although it may relax to a modified equilibrium.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Molecular dynamics and virial expansions for equilibrium thermodynamics in the solar interior

PDF

Modeling graphene: magnetic, transport and optical properties

PDF

Collision dynamics of Bose Einstein condensates: a computational-theoretic exploration at the new frontier of cold atom physics in space

Asset Metadata

Creator Cassidy, Amy C. (author)

Core Title Chaos and thermalization in the one-dimensional Bose-Hubbard model in the classical-field approximation

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 05/08/2009

Defense Date 03/12/2009

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

Tag Ablowitz-Ladik,classical chaos,discrete nonlinear Schrodinger equation,near-integrability,OAI-PMH Harvest,one-dimensional Bose-Hubbard,thermalization

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Haas, Stephan W. (committee chair), Bozler, Hans M. (committee member), Dappen, Werner (committee member), Krylov, Anna I. (committee member), Olshanii, Maxim (committee member)

Creator Email amy_c_cassidy@yahoo.com,amycassi@usc.edu

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

Unique identifier UC194859

Identifier etd-Cassidy-2749 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-244558 (legacy record id),usctheses-m2217 (legacy record id)

Legacy Identifier etd-Cassidy-2749.pdf

Dmrecord 244558

Document Type Dissertation

Rights Cassidy, Amy C.

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