Two simple systems with cold atoms: quantum chaos tests and nonequilibrium dynamics

This article is an attempt to provide a link between the quantum nonequilibrium dynamics of cold gases and fifty years of progress in the lowdimensional quantum chaos. We identify two atomic systems lying on the interface: two interacting atoms in a harmonic multimode waveguide and an interacting two-component Bose-Bose mixture in a double-well potential. In particular, we study the level spacing distribution, the wavefunction statistics, the eigenstate thermalization, and the ability to thermalize in a relaxation process as such.


Introduction
In this article, we identify two simple atomic system lying on the interface between the nonequilibrium dynamics of cold gases and low-dimensional quantum chaos. The system of two short-range-interacting atoms in a harmonic waveguide is compared to the well studiedŠeba billiard [1]. A two-component interacting Bose-Bose mixture in a double-well potential provides a realization of a quantum particle moving on a four-dimensional surface of a tensor product of two spheres under a non-integrable perturbation.
The waveguide problem allows for an exact analytic solution, in spite of the absence of a complete set of integrals of motion. The double-well problem reduces to a numerical diagonalization of a small, a-few-hundred-by-a-few-hundred matrix, a few minute calculation on a moderate laptop.
To assess the degree of quantum chaos, we utilize two standard quantumchaotic measures: level spacing statistics [2,3] and the statistics of the values of the wavefunction in the coordinate representation [4,5]. The former predicts that the spectra of quantum-chaotic systems have the statistical properties of the spectra of random matrices. The latter dictates the Gaussian statistics for the wavefunction.
From the point of view of the nonequilibrium dynamics, we study both the degree of the eigenstate thermalization and the actual thermalization in an expansion from a class of realistic initial states. The eigenstate thermalization [6,7,8,9,10] -the suppression of the eigenstate-by-eigenstate variance of quantum expectation values of simple observables -provides an ultimate upper bound for a possible deviation of the relaxed value of an observable from its thermodynamical expectation, for any initial state in principle. However, recently a new direction of research has emerged: quantum quench in many-body interacting systems [11,12,13,14,15,16,17,18,19,20,21]. In this class of problems, the initial state is inevitably decomposed into a large superposition of the eigenstates of the Hamiltonian governing the dynamics, and the discrepancy between the result of the relaxation and thermal values is expected to be diminished.
Note that a related comprehensive study of the relationship between the quantum chaos and thermalization for atoms in optical lattices has been performed in Ref. [22].

Two bosons in a circular, transversely harmonic multimode waveguide
In the case of the relative motion of two short-range-interacting atoms in a circular, transversely harmonic waveguide [23] , the full Hamiltonian can be split onto an integrable, "free motion" part and a non-integrable perturbation. The unperturbed Hamiltonian is given by the sum of the longitudinal and transverse kinetic energies and the transverse trapping energyÛ = µω 2 ⊥ ρ 2 /2, where ω ⊥ is the transverse frequency, ∆ ρ is the transverse two-dimensional Laplacian, µ = m/2 is the reduced mass, and m is the atomic mass. We will assume periodic boundary conditions along z, with a period L. In what follows, we will restrict the Hilbert space to the states of zero z-component of the angular momentum and even under the z ↔ −z reflection. Note that the zero-range interaction has no effect on the rest of the Hilbert space. The non-interacting eigenstates are products of the transverse two-dimensional zero-angular-momentum harmonic wavefunctions, labeled by the quantum number n ≥ 0, and the symmetric plane waves cos(2πlz/L), l ≥ 0. The unperturbed spectrum is therefore given by E nl = 2 ω ⊥ n + 2 (2πl/L) 2 /(2µ) The interactionV F H in (2) couples the transverse and longitudinal degrees of freedom. We assume the interaction to be of the Fermi-Huang type. It can be symbolically represented as a separable rank I perturbation [24]: with the formfactors |L = δ 3 (r) and R| = δ 3 (r)(∂/∂r)(r ·). Here, the interaction strength is V = (2π 2 a s /µ), where a s is the three-dimensional s-wave scattering length. The rigorous definition of the operator (3) reads Following the the Refs. [24,1], it can be shown that the operator (4) produces a member of a (parametrized by V ) family of self-adjoint extensions of the "free motion" Hamiltonian (2). Note that the model (1), with an unbounded along z motion, was first used to derive the effective one-dimensional atom-atom interaction in a monomode atom guide [25]. It was further extended to the multimode regime in Ref. [26], where the intermode kinetic coefficients were derived. In this article, we consider a bound-motion analogue of the multimode guide of Ref. [26], suitable for a study of the relaxation dynamics and chaotic properties.
Even though the perturbation of the form (3) destroys the complete set of integrals of motion of the unperturbed system, the problem can be solved exactly. The solution for the full system can be obtained as follows. The action of the perturbation on each wavefunction produces the left formfactor |L . As a result, any eigenstate |α of the perturbed system functionally coincides with the Green function of the unperturbed Hamiltonian taken at the energy of the eigenstate E α . The following eigenenergy equation then applies, along with an expression for the eigenstates: where | n stands for an eigenstate of the unperturbed system of energy E n . Similar expressions were obtained in Refs. [1,27] for the case of theŠeba billiard and its generalizations. Solutions of this type, for the integrable case of a spherically symmetric harmonic confinement, are presented in Ref. [28]. Substituting the above non-interacting eigenstates and eigenenergies to the Eq. (5) and summing over l, one obtains the following expression for the interacting eigenstates: where the rescaled energy ǫ α is given by n (ξ) are the Legendre polynomials, and a ⊥ = ( /µω ⊥ ) 1/2 is the size of the transverse ground state. The normalization constant C α is expressed as Extracting the regular part of the wavefunction (6) via the procedure given in Refs. [26,29], we arrive at the following transcendental equation for the eigenenergies where ζ (ν, x) is the Hurwitz ζ-function (see [29]). The summands in the sums Eqs. (6), (7), and (8) decay exponentially with n, leading to the fast converging series. Note that the imaginary parts of the two terms in the left hand side of the (8) cancel each other automatically. Similar solutions were obtained for two atoms with a zerorange interaction in a cylindrically-symmetric harmonic potential [30] (that system was analyzed numerically in Ref. [31]). Higher partial wave scatterers were analyzed in the Ref. [32].
At rational values of λ/π 2 the unperturbed energy spectrum shows degeneracies that are not fully lifted in the full [deduced from (8)] spectrum. To minimize the effect of the degeneracies, we, following Ref. [1], fix the length of the guide by (L/a ⊥ ) 2 ≡ λ = 2φ gr π 7 ≈ 9774, where φ gr = (1 + √ 5)/2 is the golden ratio.

A two-component interacting Bose gas in two coupled potential wells
Next, we are going to consider a two-component Bose-Bose mixture trapped in a two-well potential. The second-quantized Hamiltonian iŝ whereĉ † iL (ĉ † iR ) is the creation operator for the atom of the i-th type in the left(right) well, a and b denote the two species, and the factor of −2 in the definition of V is used to ensure the consistency with the conventional definition of the hopping constant in the theory of multi-well lattices with periodic boundary conditions. We fix the coupling constants U ijL(R) to U aaR = 0.3041 U, U abR = 1.0410 U, U bbR = 0.3684 U, U aaL = 0.3085 U, U abL = 1.4574 U, U bbL = 0.4524 U , where U is the overall coupling strength. The relative strengths of the interaction were chosen once and for all from a uniform random distribution between 0 and 2. The same set was used in all further numerical experiments.
The set of deviations of the number of atoms of a given type in the right well from the half of the total number of atoms of this type constitutes a convenient set of quantum numbers: The observableŝ a has been chosen to be the principal observable of interest in the eigenstate thermalization and relaxation studies. Below, we fix the numbers of atoms to N a = N b = N = 20.
Analogously to the one-specie case [33], the Hamiltonian (9) is integrable in both J ≫ U and J ≪ U limits.
In the J ≪ U limit, the eigenstates are represented by the Fock states (eigenstates ofĤ 0 , see (10)) where number of bosons of each type in each well is fixed. For our choice of the coupling constants, the energy as a function of (s a , s b ) constitutes a hyperboloid. At low energies, the spectrum ofĤ 0 is dominated by the eigenstates where the atoms of different type are localized in the opposite wells: (s a , s b ) ≈ (±N, ∓N ). At high energies, the atoms prefer to localize in the same well: (s a , s b ) ≈ (±N, ±N ).
In the opposite limit J ≫ U , the good quantum numbers are given by the numbers of particles of each type in the even and odd (with respect to the wells) one-body states. In this limit, the mean-field Hamiltonian , with the hopping as the principal component, generates a good approximation to the energy spectrum. The spectra of the full HamiltonianĤ (9), the pure-interaction HamiltonianĤ 0 (10), and the pure-hopping-mean-field HamiltonianĤ mf are shown at Fig. 9a. Note the degeneracies in the spectrum of the pure-hopping Hamiltonian H mf . Recall that in our model, the hopping constants are the same for both types of atoms; as a result, each energy level corresponding to the total of N + atoms in the even state ( where the exact, pure-interaction, and pure-hopping eigenenergy curves coincide corresponds to the maximal depolarization, ŝ a(b) ≈ ĉ † a(b)Lĉa(b)R ≈ 0, where neither of the two integrable limiting models apply. This is the point where the motion is expected to be "maximally chaotic".
The attractive feature of the two-well system is that it can be mapped to a system consisting of a single particle moving on the four-dimensional surface of a tensor product of the two two-dimensional spheres. The map is |s a , s b → |(L a = N a /2, m a = s a ), ( is the total angular momentum on the first(second) sphere, and m a (m b ) is the corresponding projection to the z-axis. The corresponding Hamiltonian reads: with the Hamiltonian (13). In the many-boson picture, these momenta are related to the number of bosons of each type as L a(b) = N a(b) /2 (both assumed to be even). For the one-specie case, similar mapping was employed in Ref. [33]. From the single-particle point of view, the Hamiltonian (13) describes an external perturbation of the free motion along the bi-sphere. For clarity, we do not include the constant "free" kinetic energy ∝ (ˆ L a ) 2 + (ˆ L b ) 2 in the Hamiltonian (13).
In the coordinate representation, a given state maps to a two-particle wavefunction (defined at the points of a tensor product of two unit spheres) given by where Y l, m (θ, φ) are spherical harmonics.

Level spacing statistics
The commonly-used measure to access the degree of quantum chaos the shape of level spacing distribution [3]. As usual, the starting point of the procedure is to determine the mean level density to be subsequently removed from the spectral data. The resulting energy spectrum-the so-called unfolded spectrum-does not contain any information specific to the system besides the dimensionless, homogeneous across the spectrum level spacing statistics. This statistics is the principal object of the investigation.
For the case of the waveguide, the mean level density was obtained from the semiclassical approximation to the unperturbed spectrum of the system. Such a simplification was possible due to the fact that for a singular perturbation, the density of states is exactly the same for the interacting and non-interacting systems. For the two-well system, the mean level density was approximated by a fit of the true energy spectrum to a seventh degree polynomial.
The integrable systems exhibit the Poisson statistics of the level spacings, with a peak at zero and an exponential decay at large spacings. One would expect to this type of statistics in a model the energy levels appear under no constraint, besides the given a priori mean level spacing. Indeed, in integrable systems, the energy levels corresponding to different sets of (additional to the energy) integrals of motion can approach each other infinitely closely. It turns out that this property along is sufficient to ensure the Poisson statistics [2,3].
On the contrary, in a quantum system with no non-trivial integrals of motion, the energy levels tend to repel. This repulsion leads to the appearance of a hole, around zero spacing point, in the level spacing distribution. It is remarkable however that generally, quantum chaotic systems with real Hamiltonians exhibit exactly the same level spacing distribution, namely the one of the Wigner-Dyson type. An ensemble of orthogonal matrices with a Gaussian distribution of the matrix elements (the socalled Gaussian Orthogonal Ensemble (GOE)) serves as a paradigm for this class of Hamiltonians [2,3].
A multi-dimensional harmonic oscillator constitutes a remarkable exception from the integrable vs. chaotic dichotomy [34]. The unfolded distribution function in this case vanishes at the zero spacing and shows a sharp decay at large spacings. Indeed, we have found that our waveguide system, with a harmonic confinement transversally and a periodic-boundary-conditions box longitudinally, is more sensitive to the effect of a non-integrable perturbation than a seemingly more elegant elongated threedimensional harmonic oscillator; thus our preference for the former.
The waveguide system shows systematic deviations from the GOE Random Matrix Theory predictions [2]. Note that the unattainability of a complete quantum chaos in singular billiards has been already addressed in the context of theŠeba billiard -flat two-dimensional rectangular billiard with a zero-range scatterer in the middle [1]. Even though theŠeba billiard shows some signatures of the quantum-chaotic behavior, both the level statistics [35,27,36]  eigenstates [37] exhibit substantial deviations from the conventional quantum-chaotic behavior. Our waveguide system belongs to the same class as theŠeba billiard: it is represented by an integrable system perturbed by a separable rank-one perturbation. For both systems, the eigenstates have the "Green's function" form (5). The results of the study of the level spacing distribution in the waveguide system are fully consistent with the analogous results for theŠeba billiard [35,27]. One can show that for a s ≫ a ⊥ (globally unitary regime), the scattering strength reaches the unitary limit for all eigenenergies of the full Hamiltonian. In this regime the distribution quickly converges to theŠeba distribution [35]. This distribution does show a gap at small level spacings but fails to reproduce the Gaussian tail predicted by the GOE. At small a s the level spacing distribution tends to the Poisson one (see Fig. 1). For the energy range E 100 ω ⊥ and the aspect ratio (L/a ⊥ ) 2 = 2φ gr π 7 theŠeba distribution is approached at a s 10 −1 a ⊥ .
The level spacing statistics for the two-well system is consistent with the GOE prediction (Fig. 2), even though the statistics is very poor.

Wavefunction statistics
Another important measure of quantum chaos is the statistics of the wavefunction. Divide the coordinate space of a quantum-chaotic system into regions larger that the  . Forty bosons in two wells. The real and imaginary parts of the fourdimensional eigenstate wavefunction living on the surface of a tensor product of two unit spheres. The "maximally chaotic" state α = 190 is used. The hopping is fixed to J/U = 2. Even though in both J ≪ U and J ≫ U limits the eigenstates are strongly polarized, for J/U = 2 no preferential direction is visible. The rest of the parameters is the same as for the Fig. 2. To enhance contrast, the shading density is made to be proportional to the forth power of the real and imaginary parts respectively.
de-Broglie wavelength but small enough so that the spatial density does not change substantially across the region. According to the Berry conjecture [4,5], the coordinate wavefunction should behave, within each region, as a sample from an ensemble of large superpositions of plane waves with random coefficients. For the waveguide system, we verified that in the globally unitary regime, a s ≫ a ⊥ , the statistics of the point-by-point variations of the eigenstate wavefunction in the spatial representation is indeed close to the Gaussian one (Fig. 3), according to Berry's prediction and a particular observation in the case of theŠeba billiard [1].
The two-well system also behaves according to the Berry conjecture. presents the distribution of the real part of the wavefunction (14) taken at random points on the bi-sphere and averaged over the eigenstates between α = 150 and α = 250. For J = 2 the distribution is indistinct from the Gaussian. Fig. 4b represents the dependence of the entropy of the distribution on the hopping constant. The entropy converges to the Berry conjecture prediction in the region between J/U ≈ 1 and J/U ≈ 5. Fig. 5 shows the wavefunction of the eigenstate α = 190 at hopping of J = 2. This state corresponds to the crossing of all three spectra of the Fig. 9a where the expectation values of the projection of the angular momentumˆ L a (along withˆ L b ) on both z and x axes is small. Accordingly, the nodal structure of the wavefunction appears completely random, with no preferred direction. For comparison, Fig. 6 shows the same state but for weak hopping (J = .2) where the motion is expected to be regular.   The spectrum of the mean-field HamiltonianĤ mf ≡V + Ĥ 0 (see (12), with the hopping as the primary Hamiltonian, is shown as well. The point α ≈ 190 where all three eigenenergy curves coincide corresponds to the maximal depolarization where the motion is "maximally chaotic". Here, the density is distributed uniformly over the bi-sphere with no preferred orientation. The labels α and α 0 correspond to the eigenstates of the full Hamiltonian (9) and the pure-interaction Hamiltonian (J/U = 0 limit) respectively. (b) Fluctuations of the observablê sa ≡ĉ † aRĉaR −Na/2. Green crosses correspond to the microcanonical fluctuations in the absence of hopping. Blue crosses correspond to the eigenstate-by-eigenstate fluctuations of the quantum expectation values, suppressed according to the eigenstate thermalization hypothesis. Red crosses give the infinite time average of the quantum expectation for a propagation from an initial state corresponding to one of the eigenstates of the HamiltonianĤ 0 . The fluctuations of the latter are further suppressed due to an additional averaging of the initial state over the eigenstates. The microcanonical average is shown as the blue dashed line.

Eigenstate thermalization phenomenon and the ability to thermalize in a relaxation from an initial state
According to the Eigenstate Thermalization Hypothesis [6,7,9,10], the ability of an isolated quantum system to thermalize follows from suppression of the eigenstateby-eigenstate fluctuations of the quantum expectation values of relevant observables α|Â|α over the eigenstates |α of the system. Indeed, if one assumes that an observableÂ reaches its thermal value in a time evolution from any initial state with an energy close to a given energy E, including the exact eigenstates of the Hamiltonian, then it immediately follows that the matrix elements α 1 |Â|α 1 and α 2 |Â|α 2 must be close to each other if the corresponding energies E 1 and E 2 are close. Consider first the system of two particles in a waveguide. As the principal observable of interest, we choose the transverse trapping energy U (see (2)). In our analysis, the potential energy U plays a role of a generic observable in an integrable system: it acts on only one of the degrees of freedom and it has both diagonal and offdiagonal non-zero matrix elements between the eigenstates of the integrable system.
Scattered grey squares at the Fig. 7a show the quantum expectation values of the transverse trapping energy, α|Û |α , in every hundreds eigenstate (6) of our system (1), relative to its microcanonical expectation value U therm. ≈ E/3 (in accordance with the energy equipartition between the two harmonic and one free-space degrees of freedom).
Even for the energies much larger than any conceivable energy scale of the system, the variation of α|Û |α is only marginally lower than the microcanonical fluctuations in the unperturbed system (2), represented by the blue circles. Figure 8 demonstrates that even though the destruction of the integrals of motion in the interacting system does narrow the distribution of the diagonal matrix elements α|Û |α and it does shift the peak towards the microcanonical prediction U therm. , the distribution remains relatively broad. We trace this effect to the predicted unattainability of the complete quantum chaos in the singular billiards [35,36,38].
Here and below, we assume the globally unitary value of the scattering length, a s = 10 6 a ⊥ .
On the contrary, the behavior of the two-well system is fully consistent with the Eigenstate Thermalization Hypothesis [6,7,9,10] (see Fig. 9b). Here, the observablê s a (the deviation number of the type-a particles in the right well from N/2) is the primary observable of interest. In the window 150 α 250, all quantum expectation values ofŝ a in the eigenstates of the full Hamiltonian (9) are substantially closer to the thermal prediction than their unperturbed counterparts given by the eigenstates |α 0 = |s a0 , s b0 of the Hamiltonian (10). Now we are going to address directly the ability of our system to thermalize from an initial state |ψ(t = 0) . In this case, the infinite time average of the quantum expectation value of an observableÂ will be given by where |α are the eigenstates of the system (see [10]).
First, let us attract our attention to the waveguide system.

Two simple systems with cold atoms: quantum chaos tests and nonequilibrium dynamics15
Consider the following initial state: Longitudinally, the state is represented by the ground state of a length Lδ hard-wall box split initially by an ideal beamsplitter with momenta ±2πl 0 /L; the box is centered at the maximal interatomic distance. The state is distributed among approximately π/δ axial modes localized about l = ±l 0 . The initial transverse state is limited to a single mode n 0 . Note that for n 0 = 0 and δ ≈ 1 the state is conceptually similar to the initial state used in the equilibration experiments [19].
The transverse trapping energy U (see (2)) remains very far from the equilibrium predictions after the relaxation. Fig. 7b shows the equilibrium value of the transverse trapping energy U for three sequences of the initial states of the type (15). For the first two, the energy is controlled via the kinetic energy given by the beamsplitter, while the transverse state is fixed to the ground state. Both show an approximately -40% deviation of the relaxed value of U from the thermal prediction. Note that the initial states of the second sequence have a much greater spread over the longitudinal modes than the first one; consistently, the energy-to-energy variation of the relaxed values of U is less than for the first sequence. In the third sequence, there is no beamsplitting, and the energy is controlled via the initial transverse energy. In this case, the deviation from the thermal prediction ranges between +50% and +5% but never reaches zero.
Another type of the initial state ρ, z|ψ(t = 0) = C ax cos πζ δ θ δ 2 − |ζ| cos (2πl 0 ζ) allows for an additional spread over the transverse modes. The transverse wavefunction vanishes at the waveguide axis. The mode occupation has a minimum at n = 0, and the occupation of the ground transverse mode tends to zero when κ 1 < κ 2 ≪ 1. The results for a sequence similar to the third sequence described above are shown at Fig. 7b as well. In spite of the significant differences in the initial distribution over both the transverse and the longitudinal modes, the deviations from the equilibrium are close to the ones for the third sequence, even though the energyto-energy variations for the former are less than for the latter. For our second system, we use the eigenstates of the unperturbed Hamiltonian (10) as the initial states of the relaxation process. The results are shown at Fig. 9b.
Here, each red cross shows the infinite time average of the quantum expectation value of the observableŝ a (the deviation of the type-a occupation of the right well from N a /2) for a relaxation process originating from the state |α 0 = |s a0 , s b0 . In the window 150 α 250 the observable thermalizes fully. Recall, that the eigenstate themalization occurs in the same window.
It is instructive to analyze our results from the point of view of the macroscopic quantum self-trapping, predicted in Refs. [33,39,40] and experimentally observed in Ref. [41], for the single-specie two-well bosonic system. From the classical field theory point of view it is a one-dimensional Hamiltonian system similar to a pendulum. Here, the population imbalance between the two wells plays a role of momentum, while the relative phase between the wells plays a role of a coordinate. Similarly to the case of the conventional pendulum, the initial states of sufficiently high momentum preserve the sign of the latter; the time evolution that starts from the states with even higher momentum shows the momentum that almost does not change in time. The states with the opposite sign of the momentum become mutually inaccessible. In terms of the original two-well system, the initial states with the population imbalance higher than a certain critical value preserve the imbalance over time.
The energy landscape for out two-specie two-well system is more complicated. Both the cosine-like dependence of the hopping energy on the relative phases and the hyperbolic shape of the interaction energy surface as a function of the two population imbalances contribute to the potential multi-connectedness of the phase space. At the same time, the initial states used to produce the Fig. 9b are the exact eigenstates of both population imbalances; they leave the relative phases completely uncertain.
We have identified a set of absolutely localized pairs of population imbalances, i.e. the pairs for which no set of imbalances of opposite signature is classically accessible for any pair of initial phases. In particular, the analysis shows that all the unperturbed eigenstates with indices α 0 < 39 or α 0 > 339 constitute the absolutely localized initial states. This result is consistent with the fully quantum calculation that shows a strong memory of the initial conditions (see Fig. 9b).

Summary of results
In this article, we analyze the behavior of two simple atomic systems within two complementary frameworks: quantum chaos and nonequilibrium dynamics. The first system is represented by two short-range-interacting bosonic atoms in a circular, transversely harmonic multimode waveguide. The second is the interacting twocomponent Bose-Bose mixture in two coupled potential wells. For both systems, we study the level spacing statistics, the statistics of the values of the wavefunction in coordinate representation, degree of the eigenstate thermalization, and the thermalizability in a relaxation from an exited initial state. For the waveguide system, the wavefunction statistics is fully consistent with the Berry conjecture [4] that predicts a Gaussian distribution for the quantum-chaotic systems. However, the rest of the tests show an incomplete chaotization. We trace this effect to the previously predicted incomplete quantum chaos in the singular billiards [35,27,36,37].
For the double-well system, we identify a quantum-chaotic region of the spectrum that spans, for our set of parameters, about a quarter of the full spectrum. In this part of the spectrum, the system exhibits both strong eigenstate thermalization and the ability to fully thermalize from an initial state.