Relaxation, chaos, and thermalization in a three-mode model of a Bose–Einstein condensate

The emergence of new quantum simulators have allowed for a better understanding, description, and control of quantum many-body systems out of equilibrium. Progressively, approaches are found to explain and to take advantage of a variety of different factors that affect the dynamics of these complex systems, including range and strength of the interactions [1, 2], choice of initial state [3], presence of disorder [4], onset of quantum chaos [5, 6], and proximity to critical points [7]. Different behaviors have been identified at different time scales [8, 9], protocols to reach quantum speed limits have been engineered [10, 11], and conditions for isolated quantum systems to relax to equilibrium and to thermalize have been established [12–16].

Relaxation in an isolated quantum system implies the approach of a set observables to a stationary value, deviations of which are very rare and negligible for long-time averages. The conditions required for relaxation to occur have been discussed in [17–27]. Thermalization is associated with the fact that the state of the system reached at long times cannot seemingly be distinguished from an ad hoc defined thermal distribution [13–15, 28] (ad hoc as it is defined for the particular closed system).

Experiments with cold atoms have a prominent place in studies of relaxation and thermalization due to their high level of isolation [12, 16, 29–36]. Their access to precise coherence manipulation and to the preparation of desired initial states are essential for the investigation of quantum many-body dynamics [37]. In this context, a good example is a set of experiments with a quasi-one dimensional (quasi-1D) Bose–Einstein condensate (BEC) on an atom chip [11, 38, 39], which successfully performed coherent transfer between motional states of the transverse trapping potential. This allowed for the preparation of the condensate in a coherent superposition of the two lowest motional states. This superposition then evolved in the trapping potential.

The details of the dynamics of the above mentioned quasi-1D BEC are not yet entirely understood. At short times, the evolution of the initial superposition presents oscillations of the atomic density distribution [39], which agrees with simulations based on a quasi-1D Gross–Pitaevskii equation (GPE). At longer time, the density distribution relaxes to a steady state [40]. We note here that there are theoretical studies on three coherent modes, which display many interesting phenomena [41–43]. Nevertheless, the damping of the oscillations is not captured by the mean-field (GPE) approximation with the physical parameters of the system. That is the GPE approach predicts perpetual oscillations of the atomic density distribution. In this article, we investigate how a simple quantum many-body model can provide an explanation for the relaxation of this isolated system. The model shows relaxation and thermalization, and hence it is a testbed for the analysis of theoretical bounds on relaxation times.

Before discussing the quantum many-body model, we consider first a semiclassical model based on three modes. Similarly to the GPE, it accounts for the initial oscillations, but cannot explain their decay. A semiclassical two-mode model is even worse, being incapable of qualitatively describing the initial oscillations. The system investigated in [39, 40] consisted of a degenerate gas of several hundred ⁸⁷Rb atoms, which justified the use of the GPE. Mean-field approximations effectively describe various phenomena in BEC and in many cases it is preferred over many-body approaches, which are computationally more involved and often intractable. However, mean-field approximations are by construction blind to the microscopic properties of individual atoms and do not account for collisions or quantum fluctuations.

Our three-mode quantum many-body model initially prepared in the two lowest modes accounts for both the oscillations and their damping. As we show, the decay of the oscillations occurs as a pure quantum phenomenon and provides a neat example of relaxation in an isolated quantum system. By extrapolating the number of atoms reachable by our numerical tests to the number of atoms used in the experiment, we extract a value for the damping time. It is larger than the damping time in the experiment, but reproduces the decay of oscillations qualitatively.

Note that when referring to the decay of the oscillations, we employ the words damping, decay, and relaxation on an equal footing. However, strictly speaking, the isolated three-mode model cannot account for damping processes as in conventional open quantum system approaches [44, 45], where information is irreversibly lost to the environment. Our quantum model experiences dephasing, similar to what is termed collapse and revival of the wave function in quantum optics [46]. The collapse occurs since the initial state is a superposition of exact Hamiltonian many-body eigenstates with eigenenergies that are anharmonic. This anharmonicity of many-body eigenstates is the source of the 'collapse'-dephasing. Since the system has strictly speaking a discrete spectrum, in addition to dephasing, an 'approximate revival' of the intial state is expected to occur at the time scales of the order of the inverse of the smallest/typical gap between neighboring levels in the energy spectrum (for the discussion in the context of thermalization see [47]). Still, the parallel of 'dephasing' and thermalization in open systems can be drawn. Due to its complexity, our system plays the role of its own environment. Specifically, the second excited mode can be understood as a minimal environment, while the system corresponds to the initially populated ground and first excited modes. This is an intuitive interpretation drawn from the fact that the second excited mode is not initially macroscopically occupied and its consideration represents the first step towards a more general approximation: one can consider that all Bogoliubov excitations are the environment for the two lowest modes in the same spirit as done when interpreting the Bose polaron motion as that of a quantum Brownian particle [48].

We relate the decay of the oscillations with the fragmentation (loss of coherence) of the condensate and with the approach of the quantum three-mode model to the chaotic regime. The connection between damping and fragmentation is in accordance with numerical studies done for a BEC in a quasi-1D bosonic Josephson junction using the multi-configurational time-dependent Hartree method for bosons (MCTDHB) [49]. We also show that while the two-mode model can account for the loss of coherence and fragmentation, it does not show a quantum chaotic regime.

Quantum chaos is associated with level repulsion and highly delocalized eigenstates, both of which guarantee the fast relaxation and thermalization of systems perturbed far from equilibrium [13–15]. In the scenario of isolated quantum systems, fragmentation, relaxation, and thermalization are caused by the interparticle interaction, rather than by couplings with an external thermal bath. They therefore reinforce the fact that the mean-field approximation is not valid for long times.

The paper is organized as follows. Section 2 introduces the three-mode many-body model considered. Sections 3 and 4 analyze how its properties change as the interaction strength increases from zero. In section 3, the dynamics under the two- and three-mode semiclassical models are compared with the quantum three-mode model. Two kinds of oscillations, their decay, and the phenomenon of fragmentation are discussed. Section 4 addresses the onset of chaos and its connection with the decay of the oscillations. Conclusions are given in section 5. We also present the study of the quantum dynamics for the two-mode model in appendix A and discuss some previous results about the conditions for relaxation in isolated quantum systems in appendix B.

The second-quantized Hamiltonian for N ultracold bosons in an external potential V(x) is given by

$$\begin{eqnarray}&&H=\int {\rm{d}}{x}{{\rm{\Psi }}}^{\dagger }({x})\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+V({x})\right){\rm{\Psi }}({x})+\displaystyle \frac{{g}_{3{\rm{D}}}}{2}\int {\rm{d}}{x}{{\rm{\Psi }}}^{\dagger }({x}){{\rm{\Psi }}}^{\dagger }({x}){\rm{\Psi }}({x}){\rm{\Psi }}({x}),\end{eqnarray} \tag{ 1 }$$

where Ψ(x) represents the field operator, ${\bf{x}}\in {{\mathbb{R}}}^{3}$, ${\hslash }$ is the Planck constant (which will be set to 1), m is the mass of the particles, and g_3D is the coupling constant in three dimensions governing the contact interactions. The hats on the operators are omitted to simplify the notation. Let us first consider that the trapping potential V(x) is parabolic in the x and z directions, characterized by the trapping frequencies ω_x,z. We assume the trapping in the y direction is slightly anharmonic. In such case all single-particle eigenenergies are discrete and show accidental degeneracies given by ${E}_{{j}_{x},j,{j}_{z}}=(1/2+{j}_{x}){\omega }_{x}+(1/2+{j}_{z}){\omega }_{z}+{E}_{j}$, with j_x,z integers and E_j the energy of jth level of the separable single-particle Hamiltonian in the y direction.

For clarity of the explanation, we first introduce the one dimensional model. We assume that ω_x,z ≫ E₀₁ = E₁−E₀, so that we can consider dynamics only in the transverse (y) direction. We aim for a model able to describe the dynamics of an initial state where a BEC is prepared in a coherent superposition of the lowest motional states along y. We expand the field operator in terms of eigenfunctions of the single-particle Hamiltonian ψ_j(η), η = x, y, z i.e. ${\rm{\Psi }}={\sum }_{{j}_{x},j,{j}_{z}}{a}_{{j}_{x},j,{j}_{z}}{\psi }_{{j}_{x}}(x){\psi }_{j}(y){\psi }_{{j}_{z}}(z)$, with ${a}_{{j}_{x},j,{j}_{z}}$ the annihilation operator. As the atoms are tightly confined in x and z the access to excited states in these directions is in practice forbidden, so the dynamics is frozen in those directions. Thus, as we only use in practice the operators a_0j0, from here on we only keep the j index. We note that there is a symmetry in the y direction, as the Hamiltonian expressed in this basis is invariant by the reflection y → −y. According to our premises, the atoms initially occupy macroscopically the two lowest modes in the y direction. We truncate the expansion of the field operator in the third mode

$$\begin{eqnarray}&&{\rm{\Psi }}={a}_{0}{\psi }_{0}+{a}_{1}{\psi }_{1}+{a}_{2}{\psi }_{2}\end{eqnarray} \tag{ 2 }$$

which produces the Hamiltonian

$$\begin{eqnarray}\begin{array}{rcl}{H}_{3{\rm{m}}} & = & \displaystyle \sum _{i}{n}_{i}{E}_{i}+\displaystyle \sum _{i\ne j}{U}_{{ij}}[{({a}_{i}^{\dagger })}^{2}{a}_{j}^{2}+{({a}_{j}^{\dagger })}^{2}{a}_{i}^{2}]+\displaystyle \sum _{i}{U}_{i}{n}_{i}({n}_{i}-1)\\ & & +\ 4\displaystyle \sum _{i\ne j}{U}_{{ij}}{n}_{i}{n}_{j}+2{U}_{1120}[{({a}_{1}^{\dagger })}^{2}{a}_{0}{a}_{2}+{a}_{1}^{2}{a}_{0}^{\dagger }{a}_{2}^{\dagger }]+4{U}_{1120}[{n}_{1}{a}_{0}^{\dagger }{a}_{2}+{n}_{1}{a}_{2}^{\dagger }{a}_{0}]\\ & & +\ 2{U}_{0222}[{a}_{0}^{\dagger }{n}_{2}{a}_{2}+{a}_{2}^{\dagger }{a}_{0}{n}_{2}]+2{U}_{0002}[{a}_{2}^{\dagger }{n}_{0}{a}_{0}+{a}_{2}{a}_{0}^{\dagger }{n}_{0}],\end{array}\end{eqnarray} \tag{ 3 }$$

where ${U}_{{ijkl}}=(g/2)\int {\rm{d}}{y}\,{\psi }_{i}{\psi }_{j}{\psi }_{k}{\psi }_{l}$, ${U}_{{ij}}=(g/2)\int {\rm{d}}{y}\,{\psi }_{i}^{2}{\psi }_{j}^{2}$, ${U}_{i}=(g/2)\int {\rm{d}}{y}\,{\psi }_{i}^{4}$ (all wave functions are defined as real), and E_i is the energy of level i. Here, g is the effective quasi-1D coupling constant along the y direction. Figure 1 illustrates the parameters of the Hamiltonian (3). The caption explains what each parameter denotes and the processes that the arrows represent. We use the expansion in three modes, because this is the minimal model that contains all important virtual process. In appendix A, we introduce the two-mode model, which we compare with the results of the three-mode model obtained below. In the conventional approach to a BEC one assumes that fluctuations to modes other than that in which condensation occurs are negligible. Thus one neglects all quadratic terms with j > 0. In the three-mode model, however, a term like ${a}_{0}^{\dagger }{a}_{2}^{\dagger }{({a}_{1})}^{2}$ cannot be neglected. This is because it involves operators associated to modes macroscopically occupied, as initially $\langle {a}_{j}\rangle$ and $\langle {a}_{j}^{\dagger }\rangle$ are $\sim \sqrt{{N}_{j}}$, j = 0,1. Thus, one has to neglect all terms involving bilinear terms with j > 1. But we keep all bilinear terms on j = 2 that include two operators associated to the zeroth or first mode. We mention here that the aforementioned reflection symmetry y → −y translates into $[{H}_{3{\rm{m}}},P]$ with $P={(-1)}^{{n}_{1}}$, with ${n}_{1}\,=\,{a}_{1}^{\dagger }{a}_{1}$. As we discuss later, identification of this symmetry is important in the study of quantum chaos in section 4.

Zoom In Zoom Out Reset image size

To be able to model the experiments in [11, 38–40] we have to relax the one-dimensional assumption. In those experiments, a cigar-shaped BEC is produced in an elongated potential, i.e. the potential is weakly confining along x and more confining along y and very tightly confined along z. The assumption that the dynamics in the z direction is frozen is still valid, as the excited states in that direction have very large energies. We assume that the wave function in the longitudinal direction x is well described with the Thomas–Fermi profile, TF(x). In this paper we consider that initially part of the population is transferred to the first excited state. Particularly, in all examples we assume that half of the population is initially excited. In such a case, as we discuss later, the period of the density oscillations that occur even in the non-interacting case is given by E₀₁ (see section 3.1). On the other hand the excitations along x are much lower in energy, since ω_x < E₀₁. These excitations can easily occur in the system, but they are much slower. For this reason we assume that along the evolution the system remains in TF (x) along x and study only the dynamics in y. With this, a Hamiltonian formally identical to equation (3) can be obtained, with a different expression of the coefficients. They have to be calculated taking into account that N₀ and N₁ correspond to the total population when integrating the corresponding excited mode in y and TF(x), and thus all U_ijkl include the integration in x. This procedure gives rise to the interaction parameters gathered in table 1 for a system with N = 700. We term these values ${U}_{{ijkl}}^{\exp }$ as they are close to typical experimental values [11, 38–40]. We remark that, as the trapping potential along y is slightly anharmonic, the energy differences ${E}_{01}={E}_{1}-{E}_{0}$ and ${E}_{12}={E}_{2}-{E}_{1}$ are not equal.

In the numerical examples below, the interaction strengths are varied, so we take U_ijkl to be a constant g multiplied by the value from the table, that is ${U}_{{ijkl}}=g\times {U}_{{ijkl}}^{\exp }$. Since the geometry of the trapping potential does not change, the orbitals ψ_i do not change either. This means that the coupling constant g^num that we consider in the numerical examples is proportional to the experimental coupling constant, g^num = g × g^exp. Then g = 1 implies that g^num = g^exp.

With the Hamiltonian in equation (3), we derive the equations of motion. In the Heisenberg picture, they are

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{{{\rm{d}}{a}}_{j}}{{\rm{d}}{t}}=[{a}_{j},{H}_{3m}],\hspace{1cm}j=0,1,2,\end{eqnarray} \tag{ 4 }$$

which leads to

$$\begin{eqnarray}\begin{array}{rcl}{\rm{i}}\displaystyle \frac{{{\rm{d}}{a}}_{0}}{{\rm{d}}{t}} & = & 2{U}_{01}{a}_{0}^{\dagger }{a}_{1}^{2}+2{U}_{02}{a}_{0}^{\dagger }{a}_{2}^{2}+4{U}_{01}{n}_{1}{a}_{0}+4{U}_{02}{n}_{2}{a}_{0}+2{U}_{0}{n}_{0}{a}_{0}+{E}_{0}{a}_{0}\\ & & +\ 2{U}_{0112}({a}_{2}^{\dagger }{a}_{1}^{2}+2{n}_{1}{a}_{2})+2{U}_{0222}{n}_{2}{a}_{2}+2{U}_{0002}({a}_{0}^{2}{a}_{2}^{\dagger }+2{n}_{0}{a}_{2}),\end{array}\end{eqnarray} \tag{ 5 }$$

an equation similar to equation (5) for a₂, and

$$\begin{eqnarray}\begin{array}{rcl}{\rm{i}}\displaystyle \frac{{{\rm{d}}{a}}_{1}}{{\rm{d}}{t}} & = & 2{U}_{01}{a}_{1}^{\dagger }{a}_{0}^{2}+2{U}_{12}{a}_{1}^{\dagger }{a}_{2}^{2}+4{U}_{01}{n}_{0}{a}_{1}+4{U}_{12}{n}_{2}{a}_{1}+2{U}_{1}{n}_{1}{a}_{1}+{E}_{1}{a}_{1}\\ & & +\ 4{U}_{0112}({a}_{0}{a}_{2}{a}_{1}^{\dagger }+{a}_{0}^{\dagger }{a}_{2}{a}_{1}+{a}_{2}^{\dagger }{a}_{0}{a}_{1}).\end{array}\end{eqnarray} \tag{ 6 }$$

3.1. Semiclassical dynamics

To get physical insight into the role of the different parameters of equation (3) and the processes represented in figure 1, we consider a semiclassical version of the above equations of motion. For this, we neglect quantum fluctuations and treat the operators as c-numbers, ${\alpha }_{{i}}=\sqrt{{N}_{{i}}}\exp ({i}{\phi }_{{i}})$, where N_i is the amplitude and ϕ_i is the phase of the fields α_i (for details on this procedure see [50]). In this way, we obtain

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{d}}{N}}_{0}}{{\rm{d}}{t}} & = & -4{U}_{01}{N}_{0}{N}_{1}{\rm{\sin }}2({\phi }_{0}-{\phi }_{1})-4{U}_{02}{N}_{0}{N}_{2}{\rm{\sin }}2({\phi }_{0}-{\phi }_{2})\\ & & +\ 4{U}_{0112}{N}_{1}\sqrt{{N}_{0}{N}_{2}}[{\rm{\sin }}(2{\phi }_{1}-{\phi }_{2}-{\phi }_{0})+2{\rm{\sin }}({\phi }_{2}-{\phi }_{0})]\\ & & +\ 4{U}_{0222}{N}_{2}\sqrt{{N}_{0}{N}_{2}}{\rm{\sin }}({\phi }_{2}-{\phi }_{0})+4{U}_{0002}{N}_{0}\sqrt{{N}_{0}{N}_{2}}{\rm{\sin }}({\phi }_{2}-{\phi }_{0}),\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}{\phi }_{0}}{{\rm{d}}{t}} & = & -2{U}_{01}{N}_{1}{\rm{\cos }}2({\phi }_{0}-{\phi }_{1})-2{U}_{02}{N}_{2}{\rm{\cos }}2({\phi }_{0}-{\phi }_{2})\\ & & -\ 4{U}_{01}{N}_{1}-4{U}_{02}{N}_{2}-2{U}_{0}{N}_{0}-{E}_{0}\\ & & -\ 2{U}_{0112}{N}_{1}\sqrt{{N}_{2}/{N}_{0}}[{\rm{\cos }}(2{\phi }_{1}-{\phi }_{2}-{\phi }_{0})+2{\rm{\cos }}({\phi }_{2}-{\phi }_{0})]\\ & & -\ 2{U}_{0222}{N}_{2}\sqrt{{N}_{2}/{N}_{0}}{\rm{\cos }}({\phi }_{2}-{\phi }_{0})-6{U}_{0002}{N}_{0}\sqrt{{N}_{2}/{N}_{0}}{\rm{\cos }}({\phi }_{2}-{\phi }_{0}),\end{array}\end{eqnarray} \tag{ 8 }$$

with similar equations for N₂ and ϕ₂, and

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{d}}{N}}_{1}}{{\rm{d}}{t}} & = & -4{U}_{01}{N}_{0}{N}_{1}{\rm{\sin }}2({\phi }_{1}-{\phi }_{0})-4{U}_{12}{N}_{1}{N}_{2}{\rm{\sin }}2({\phi }_{1}-{\phi }_{2})\\ & & -\ 8{U}_{0112}{N}_{1}\sqrt{{N}_{0}{N}_{2}}{\rm{\sin }}(2{\phi }_{1}-{\phi }_{2}-{\phi }_{0}),\end{array}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}{t}} & = & -2{U}_{01}{N}_{0}{\rm{\cos }}2({\phi }_{1}-{\phi }_{0})-2{U}_{12}{N}_{2}{\rm{\cos }}2({\phi }_{1}-{\phi }_{2})\\ & & -\ 4{U}_{01}{N}_{0}-4{U}_{12}{N}_{2}-2{U}_{1}{N}_{1}-{E}_{1}\\ & & -\ 4{U}_{0112}\sqrt{{N}_{0}{N}_{2}}[{\rm{\cos }}(2{\phi }_{1}-{\phi }_{2}-{\phi }_{0})+2{\rm{\cos }}({\phi }_{0}-{\phi }_{2})].\end{array}\end{eqnarray} \tag{ 10 }$$

In figures 2(a) and (b) we show two exemplary dynamical evolutions of the amplitudes for the three modes. We assume that the initial condition is such that half of the atoms occupy the ground mode and the other half occupy the first excited mode. In particular, the initial state has (N₀, N₁, N₂) = (100, 100, 0) atoms and all relative phases equal to zero. Two observations are made from the two panels. First, although the third mode is not populated initially, its gets significantly populated during the evolution, as we expected. Second, the dynamics presents two types of oscillations with different timescales. A fast oscillation with a period ${T}_{\mathrm{fast}}\approx 0.5\,\mathrm{ms}$ and a slow oscillation with a longer period ${T}_{\mathrm{slow}}\approx 5\,\mathrm{ms}$.

Zoom In Zoom Out Reset image size

From inspection of equations (7)–(10) one deduces that for g = 0 (non-interacting limit, all U coefficient vanish) the amplitudes N_i are constant and the phases ϕ_i grow with a rate E_i. For the initial condition considered here one observes density oscillations in the numerical simulations (not shown here) which are only due to these running phases. Thus, this is the main origin of the density oscillations observed in the density plots at finite g (see figure 3). Because of that we adimensionalize the time in all figures with ΔE₀₁. For small g, inspection of equation (7) shows that N₀ will oscillate slightly around the initial value with amplitude proportional to U₀₁ and period proportional to E₀₁ (see first lines in equations (7) and (9)). All other terms are small because N₂ is much smaller. In our numerical simulations with small g (not shown) we observe that this is the only oscillation present in the populations. As g is made larger, there is a part of population that occupies the second mode, so that N₂ is no longer negligible. In such case we observe a second type of oscillation, which is slower and has an amplitude proportional to U₀₁₁₂.

Zoom In Zoom Out Reset image size

For comparison, we show in figures 2(c) and (d) the amplitudes of the ground and first excited modes for a semiclassical two-mode model starting with the same initial state as in figures 2(a) and (b). The evolution gets much simpler, as only the fast oscillation remains. The semiclassical equations for the two-mode model are trivially obtained by neglecting in equations (7)–(10) all variables and parameters associated with mode 2. They are equivalent to the equations of an oscillator such as the bosonic Josephson junction [51].

In figure 3(a) we depict the evolution of the corresponding density $| \psi (t,z){| }^{2}$ for the example shown in figure 2(a). This evolution resembles qualitatively the initial density oscillations observed in the experiment and also reproduced with a quasi-1D GPE description of the dynamics along y [39, 40]. On the other hand, the density evolution shown in figure 3(b), which corresponds to the case in figure 2(c), has no similarity with the experiment. These results confirm that the two-mode model is insufficient to describe even qualitatively the dynamics of the system and that the inclusion of the third mode is necessary to comprehend the complexity of the dynamics. As discussed in section 2, the two-mode model offers an oversimplified picture of the system, as it neglects important processes that populate significantly the second mode. We checked that involving more than three modes in the semiclassical description does not bring significant changes to the evolution.

The semiclassical three-mode model captures the initial oscillations as present also in the mean-field description. However, just as the GPE simulation, it does not describe any decay of these oscillations, as seen in figures 2(a) and (b) where the oscillations continue over time. In contrast, the many-body model discussed next accounts for a decay of the oscillations. We do not attempt here for a full study of the nonlinear dynamics in the semiclassical equations, as our goals relies in the quantum dynamics discussed in what follows.

3.2. Quantum many-body dynamics

To simulate the many-body dynamics, we perform the exact diagonalization of the many-body Hamiltonian in equation (3). The dimension ${ \mathcal D }$ of the Hamiltonian matrix is ${ \mathcal D }=(N+m-1)!/[N!(m-1)!]$, where m = 3 is the number of modes. For large N, ${ \mathcal D }\approx {N}^{m-1}$, as shown by using Stirling formula.

The analysis of the system's time evolution via exact diagonalization is very general and can be adapted to different models, e.g. Lipkin–Meshkov–Glick [52–54] or Bose–Hubbard Hamiltonian [55]. This approach has the advantage to be relatively simple. By comparison, other methods to describe many-body systems, such as the MCTDHB [56, 57], can include more features, but the physical phenomena at the origin of the observed features can be difficult to identify. Exact diagonalization is, however, limited by the exponential growth of the dimension ${ \mathcal D }$ with each additional mode or particle. For the three-mode model, we simulate the evolution with a total number of atoms up to N = 200.

The system is initialized in a coherent superposition of ground and first excited states, which corresponds to the experimental initial state in [39, 40],

$$\begin{eqnarray}&&| {\psi }_{\mathrm{ini}}\rangle =\displaystyle \frac{1}{\sqrt{N!}}{({c}_{0}{a}_{0}^{\dagger }+{c}_{1}{a}_{1}^{\dagger })}^{N}| \mathrm{vac}\rangle ,\end{eqnarray} \tag{ 11 }$$

where $| \mathrm{vac}\rangle$ is the vaccum and $| {c}_{0}{| }^{2}+| {c}_{1}{| }^{2}\,=\,1$. The initial average population of the two first modes is $\langle {n}_{0}(0)\rangle =| {c}_{0}{| }^{2}$ and $\langle {n}_{1}(0)\rangle =| {c}_{1}{| }^{2}$, with $\langle {n}_{i}(t)\rangle =\langle \psi (t)| {a}_{i}^{\dagger }{a}_{i}| \psi (t)\rangle$. We present results for an initial coherent state with c₀ = c₁, but no qualitative differences are observed when ${c}_{0}\ne {c}_{1}$. As expected, this initial state permits to reproduce faithfully the semiclassical results for small interacting systems. Even for largely interacting system, it reproduces the semiclassical results for the first few oscillations (see figure 2). This is not the case if one takes a Fock initial state, i.e. $| {\psi }_{\mathrm{ini}}\rangle =1/{(\sqrt{{N}_{0}!}\sqrt{{N}_{1}!})({a}_{0}^{\dagger })}^{{N}_{0}}{({a}_{1}^{\dagger })}^{{N}_{1}}| \mathrm{vac}\rangle$ (generally, we denote a Fock vector as $| {\varphi }_{k} \rangle =1/{(\sqrt{{N}_{0}!}\sqrt{{N}_{1}!}\sqrt{{N}_{2}!})({a}_{0}^{\dagger })}^{{N}_{0}}{({a}_{1}^{\dagger })}^{{N}_{1}}{({a}_{2}^{\dagger })}^{{N}_{2}}| {\rm{v}}{\rm{a}}{\rm{c}} \rangle =| {N}_{0},{N}_{1},{N}_{2} \rangle$). The reason for not using an initial Fock state is that it has a definite number of particles occupying each mode with all uncertainty translated into the relative phase, which will be undetermined. This initial state does not correspond to the initial state prepared in the experiment, where the initial relative phase is on average zero (note that we take c₀ = c₁).

To circumvent our limitation to relatively small system sizes, we increase the effective interaction constant to reach values of the product gN that are close to those found experimentally [11, 39, 40], where an example is given in table 1, and hence correspond to g = 1 and N = 700. One needs to keep in mind, however, that keeping gN constant, but varying the interaction parameter and atom number does not necessarily guarantee the same physical scenarios [49]. For a fixed and small value of gN, a small interaction constant g with a big value of N ensures that the semiclassical description is valid, while a small number of atoms N with large interactions leads to the strongly correlated quantum regime. To extend results obtained for low atom numbers (and high g) to the experimental case of high N (and low g), a mapping to a known problem (e.g. a Bose–Hubbard model in a lattice) could provide some insight, but such mapping is not always trivial.

In figures 2(e) and (f), we present two examples of the evolution of the occupation of the ground and first excited modes for N = 200 for two values of g. In figures 3(c) and (d), we show the corresponding density evolution. The figures make it evident that the many-body evolution differs from the semiclassical one. For the quantum model, the system damps to an equilibrium state after a few oscillations. The damping occurs earlier as gN is made larger. We have not found any revival even for the longest numerical simulations performed (typically 15 oscillations in terms of ΔE₀₁/2π, with the longest simulations up to t_max = 30 ΔE₀₁/2π).

3.3. Damping of the oscillations

The rest of the paper is devoted to investigating theoretically the origin of the damping of the oscillations. We evaluate numerically the damping time τ as a function of gN. To this end, we simulate the system for a number of particles N ranging between 40 and 200 and we also vary g. But before proceeding with this evaluation, it is beneficial to elaborate on some related points.

First, we note that the decay of the oscillations of all modes occupations, $\langle {a}_{i}^{\dagger }{a}_{i}\rangle$ for i = 0, 1, 2, are accompanied by a decay to zero of the coherence terms $\langle {a}_{i}^{\dagger }{a}_{j}\rangle$ with $i\ne j$. This parallel is verified by comparing the decay of the oscillations in figures 2(e), (f), and 3(d) with figures 4(a) and (b). The vanishing of the off-diagonal terms is associated with the fragmentation of the BEC, as discussed next.

Zoom In Zoom Out Reset image size

The phenomenon of fragmentation can be understood as follows. For a BEC in a trap, when there is only one large eigenvalue of the one-body density matrix (OBDM), $\rho =| {\rm{\Psi }}\rangle \langle {\rm{\Psi }}|$, the semiclassical Gross–Pitaevskii approach is appropriate. In this case, the depletion cloud of atoms which are not occupying the condensate is small. In contrast, the presence of more than one large eigenvalue indicates that the depletion cloud is large and that many atoms are not Bose–Einstein condensed [58]. For instance, in the context of double-well potentials (and generally with two-mode models), when there are two, and only two, large eigenvalues, the system is said to be fragmented [59, 60]. This means that the atoms occupying the two distinct modes of the system cease to be coherent, while coherence may still exist among the atoms occupying each individual mode. In general, fragmentation corresponds to the separation of an initially fully condensed state into two or more independent condensed parts. This scenario can be pictured as a double-well potential with an infinitely large barrier between the wells, so that the system is effectively cut in two halves, with no coherence among them. Of course, a single-shot experiment can show fringes and interference between the two condensates (see discussion in pg 343 of [61], where interference experiments in the Fock regime in a double well is discussed). We clarify here that it is in this sense that we talk about loss of coherence.

In figure 4(c), we show the time evolution of the two largest eigenvalues of the OBDM for different values of gN, taking N = 200. Initially, there is only one large eigenvalue. It decreases in time, while the second largest eigenvalue increases. After the damping occurs, one finds three eigenvalues that are significantly different from zero (only two are shown in figure 4(c), because the third one is much smaller, although it is also non-negligible, as the sum of the first two does not amount to 1). This indicates some degree of fragmentation in the two lowest modes. Since the eigenvalue corresponding to the third mode is much smaller than the other two, it should not become of order 1 as N increases, but should instead converge to zero. Therefore, its occupation should not be macroscopic in the thermodynamic limit (large N). This sustains our intuition that the third mode acts as an environment for the other two.

The loss of coherence, just as damping, occurs earlier in time as gN is made larger (compare figures 2(e), (f), and 4(c)). After relaxation, the system is found in a Fock state, that is, a state with a determined value of atoms occupying each mode. The link between damping and fragmentation has also been pointed out in two-mode models [49, 62].

We estimate the damping time τ from the evolution of the eigenvalues of the OBDM. In figure 5, we plot our numerical estimates for τ as a function of gN. Two numerical criteria are used to determine the damping time. In figure 5(a), τ is the time at which the largest eigenvalue of the OBDM gets smaller than 0.98. In figure 5(b), τ is the time when the largest eigenvalue becomes smaller than 0.85.

Zoom In Zoom Out Reset image size

For each N, we observe that the dependence of τ on gN can be fitted to a function of the form $\tau =a(N)\exp [b(N)x]$, with $x={\mathrm{log}}_{10}({gN})$. For every value of N, we fit the parameters a(N) and b(N) from the results of the numerical simulations performed with a large number of values of gN. We present the results for $A(N)=\mathrm{ln}[a(N)]$ and b(N) in figures 5(c) and (d). The corresponding fitted curves are the solid lines in figures 5(a) and (b). For every curve corresponding to a different N, we define g_damping as that in which the damping time is reduced to two oscillations in terms of ΔE₀₁/2π (see figures 5(a) and (b)). With the time adimensionalization we used, this corresponds to τ = 2. We use two oscillations as a criteria to define g_damping because we observe that, in this way, for g < g_damping the damping time increases significantly as one decreases g. This allows to distinguish from the region with g > g_damping, where the damping time is reduced strongly. For the large values of g and the atom numbers considered in figure 5, all the curves for different N show very similar behaviors.

Our numerical studies are limited to numbers of atoms up to N = 200, but one can extrapolate the damping times to the case with g = 1 and N = 700, which is the lowest number of atoms performed in the experiment (see our convention as described in the end of section 2). To this end, we fit a straight line to the coefficients a(N) and b(N) (see fitting lines in figures 5(c) and (d)). We use these fitted behaviors to extrapolate the behavior to larger number of atoms to estimate the coefficients a(700) and b(700) which correspond to the curve expected for N = 700. The result is the thick red curve in figures 5(a) and (b), where the estimated coefficient A(700) ≈ 35 and 25, respectively, while b(700) ≈ −9 and −6, respectively. According to these fittings, the damping occurs around ${\mathrm{log}}_{10}({gN})=3.6$ [3.8] with N = 700 for the criterion used in figures 5(a) [(b)]. This corresponds to g_damping ≈ 5 (g_damping ≈ 9), while in the experiment it happens at g_damping = 1 (following our convention). The damping time predicted for g = 1 and N = 700 from the curves of figures 5(a) and (b) is of the order of thousands of oscillations, much larger than the one observed in the experiment (τ_expΔE₀₁ ≤ 15) [39, 40]. Thus, even though the three-mode quantum model describes a damping of the density oscillations similarly to the experimental observations [40], the damping timescale it predicts differs from the experimental one. Some other effects may cause the shorter damping time seen in the experiment, such as dephasing dynamics in the longitudinal direction, perpendicular to the one-dimensional plane that we consider here.

In the next section, we explore the relationship between the onset of quantum chaos, the decay of the oscillations, and the fragmentation of the condensate. We numerically link the finite values of τ with the approach of the quantum model to the chaotic regime.

Isolated many-body quantum systems perturbed far from equilibrium relax quickly to a new equilibrium despite the absence of external couplings. The driving mechanism for equilibration is the internal coupling between particles [14]. While both integrable and chaotic systems undergo a similar process, relaxation to thermal equilibrium is expected only for chaotic systems.

The many-body three-mode model undergoes a transition from integrability to chaos as the interaction strength g increases from zero. This is in contrast with the many-body two-mode model, which is integrable for any value of the interaction, as discussed in appendix A.

In the many-body three-mode model, as the number of atoms N increases, smaller values of g are needed to move the system away from the integrable limit. This behavior mirrors our findings for the damping time, which also decreases as gN gets larger.

4.1. Quantum chaos

Quantum chaos refers to signatures observed at the quantum level that indicate that the classical counterpart of the system is chaotic. The concept has been extended to any quantum system that exhibits those properties even if it does not have a classical limit. A main signature of quantum chaos is level repulsion and the consequent rigidity of the spectrum.

4.1.1. Level spacing distribution

There are different ways to detect level repulsion and therefore the crossover from integrability to quantum chaos [63]. The most commonly used quantity is the distribution P of the spacings s between neighboring unfolded levels. In integrable models, the levels can cross and the distribution is usually Poisson,

$$\begin{eqnarray*}&&{P}_{P}(s)=\exp (-s),\end{eqnarray*}$$

although variations may be found. This may occur, for example, in systems with an excessive number of degeneracies or in 'picket-fence' spectra where the eigenvalues are nearly equally spaced. A typical example for the latter is the case of uncoupled harmonic oscillators [64, 65]. In chaotic systems, on the other hand, crossings are avoided and P(s) follows the Wigner–Dyson distribution, as predicted by random matrix theory. In systems described by Hamiltonian matrices that are real and symmetric, the shape of this distribution is given by

$$\begin{eqnarray*}&&{P}_{{\rm{WD}}}(s)=\displaystyle \frac{\pi s}{2}\exp \left(-\displaystyle \frac{\pi {s}^{2}}{4}\right).\end{eqnarray*}$$

In figure 6, we compare the level spacing distribution of the three-mode quantum model for different values of g and N. To get a meaningful distribution, the levels need to be separated by symmetry sector [66]. Following the description of equation (3), we separate the eigenvalues by the parity of the eigenstates. The two top rows of figure 6 are obtained for N = 140. When the interaction strength is small, g < 10, and the model is close to integrability, the level spacing distribution is not even Poisson. The distributions for g = 0.03, 0.1 suggest a 'picket-fence' spectrum [67]. For large interaction, g > 40, the transition to the Wigner–Dyson distribution is clear.

Zoom In Zoom Out Reset image size

The second and third rows of figure 6 compare P(s) for two different choices of N. As the number of atoms increases, the transition to chaos occurs for smaller values of the interaction [68]. This is evident by contrasting the panels with strong interactions (g = 40 and g = 80) for N = 140 with those for N = 220. This indicates that when N is very large, infinitesimal interactions may suffice for the onset of quantum chaos.

The two bottom panels of figure 6 show results for chaos indicators β and η. These are measures of the proximity of P(s) to Poisson or to Wigner–Dyson distributions. The indicator β is obtained by fitting P(s) with the Brody distribution [69],

$$\begin{eqnarray*}&&{P}_{B}(s)=(\beta +1){{bs}}^{\beta }\exp (-{{bs}}^{\beta +1}),b={\left[{\rm{\Gamma }}\left(\displaystyle \frac{\beta +2}{\beta +1}\right)\right]}^{\beta +1},\end{eqnarray*}$$

where Γ is Euler's gamma function. When β = 0, the distribution is Poisson and for β = 1, P(s) has the Wigner–Dyson shape. The indicator η was introduced in [70] and is defined as

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\int }_{0}^{{s}_{0}}[P(s)-{P}_{{\rm{WD}}}(s)]{\rm{d}}{s}}{{\int }_{0}^{{s}_{0}}[{P}_{P}(s)-{P}_{{\rm{WD}}}(s)]{\rm{d}}{s}},\end{eqnarray} \tag{ 12 }$$

where s₀ is the first intersection point of P_P(s) and P_WD(s). For a Poisson distribution, $\eta \to 1$, and for the Wigner–Dyson distribution, $\eta \to 0$.

In figure 6, the results for β and η for g < 10 need to be taken with care. For the numbers of atoms accessible to us, this range of interaction strengths leads to shapes other than P_P, P_WD, or any intermediate distribution between the two, as seen in the first row of figure 6. The fact that for g < 10, the indicator β (η) increases (decreases) as g and N get smaller simply indicates that we move away from the Poisson distribution, but this is not accompanied by an approach to P_WD. We instead approach the integrable point of three uncoupled oscillators, $H\sim {\sum }_{i}{n}_{i}{E}_{i}$.

For the numbers of atoms considered here, the transition from Poisson to Wigner–Dyson is well captured by the chaos indicators when g > 10. The plots for β and η in figure 6 reinforce our statement above that the transition to chaos happens for smaller values of g as N increases.

4.1.2. Structure of the eigenstates

The emergence of random matrix statistics is tightly connected with the appearance of chaotic eigenstates, that is states that are highly delocalized and fill the energy shell [71–73]. To measure the level of delocalization of the eigenstates $| {\psi }_{\nu }\rangle$, one can use quantities such as the participation ratio,

$$\begin{eqnarray}&&{\mathrm{PR}}^{(\nu )}=\displaystyle \frac{1}{{\sum }_{j}| {C}_{j}^{(\nu )}{| }^{4}},\end{eqnarray} \tag{ 13 }$$

where ${C}_{j}^{(\nu )}=\langle {\varphi }_{j}| {\psi }_{\nu }\rangle$ is the overlap between the eigenstate $| {\psi }_{\nu }\rangle$ and the basis vector $| {\varphi }_{j}\rangle$. PR is large when the eigenstate is delocalized in the chosen basis. The choice of basis for the analysis of the structure of the eigenstates is physically motivated. For the three-mode model, we select the Fock basis, $| {\varphi }_{j}\rangle =| {N}_{0},{N}_{1},{N}_{2}\rangle$. In the absence of interaction, when the eigenstates coincide with the basis vectors, PR = 1.

Each panel of the two top rows of figure 7 show the values of PR for all eigenstates. Different values of g are considered. The level of delocalization increases significantly with the interaction strength. One also notices that the highest values of PR occur close to the middle of the spectrum. This reflects the shape of the density of states ρ, shown in the bottom row of figure 7 for comparison⁵ .The density of states peaks close to the middle of the spectrum, where the largest concentration of eigenstates is found. This is the region where we expect the eigenstates to be more delocalized states, while at the borders, PR is smaller.

Zoom In Zoom Out Reset image size

The middle row of figure 7 illustrates the consequence of the transition to chaos. For g < 40 and thus away from the chaotic regime, there are large fluctuations in the values of PR. This implies that eigenstates very close in energy can have very different levels of delocalization. In contrast, in the chaotic region (g = 80), the structures of the eigenstates become very similar, especially close to the middle of the spectrum, where PR becomes a smoother function of energy. At this point, the states approach random vectors. The similarity between eigenstates very close in energy is what guarantees the validity of the eigenstate thermalization hypothesis (ETH) and the viability of thermalization [68, 74, 75], as discussed next.

4.2. Thermalization

The analysis of the onset of thermalization involves two steps. First one needs to ensure that the system equilibrates. Next, we verify whether the equilibrium is thermal or not.

4.2.1. Equilibration

How the isolated system reaches equilibrium is the subject of the broad field of nonequilibrium quantum dynamics to which the previous section and several other works have been devoted to, including studies about pre-thermalization [34, 76]. A brief discussion about the subject is presented in appendix B. In this section, we are concerned with the equilibrium point itself.

One can say that isolated quantum many-body systems without too many degeneracies equilibrate, because revivals become rare and take exceedingly long times to happen as the system size increases. For all practical purposes, the coherences are irreversibly lost. The systems equilibrate in a probabilistic sense. To better explain what we mean by this, consider a general observable O evolving in time according to the equation

$$\begin{eqnarray}\begin{array}{rcl}O(t) & = & \langle \psi (0)| {{\rm{e}}}^{{\rm{i}}{H}{t}}{{O}{\rm{e}}}^{-{\rm{i}}{H}{t}}| \psi (0)\rangle \\ & = & \ \displaystyle \sum _{\nu \ne \mu =1}^{{ \mathcal D }}{C}_{\mathrm{ini}}^{\mu * }{C}_{\mathrm{ini}}^{\nu }{O}_{\mu \nu }{{\rm{e}}}^{i({E}_{\mu }-{E}_{\nu })t}+\displaystyle \sum _{\nu =1}^{{ \mathcal D }}| {C}_{\mathrm{ini}}^{\nu }{| }^{2}{O}_{\nu \nu },\end{array}\end{eqnarray} \tag{ 14 }$$

where $| {\psi }_{\nu }\rangle$ and E_ν are the eigenstates and eigenvalues of H, 'ini' indicates the initial state, ${O}_{\mu \nu }=\langle {\psi }_{\mu }| O| {\psi }_{\nu }\rangle$, and O_νν is the eigenstate expectation value (EEV) of O. After a transient time, the system is said to have reached a new equilibrium if O(t) simply fluctuates around the infinite-time average,

$$\begin{eqnarray}&&\overline{O}={O}_{\mathrm{DE}}=\displaystyle \sum _{\nu =1}^{{ \mathcal D }}| {C}_{\mathrm{ini}}^{\nu }{| }^{2}{O}_{\nu \nu },\end{eqnarray} \tag{ 15 }$$

and remains very close to this value for most times. Since the infinite-time average only involves the diagonal matrix elements O_νν, this average is often referred to as 'diagonal ensemble' (DE) average.

To talk about equilibration, it is therefore essential that the fluctuations around O_DE be small and decrease with system size [17–27]. Equilibration does not require chaos in the sense of level repulsion, but it needs highly delocalized eigenstates, delocalized initial states, and not too many degeneracies.

As shown in figures 2(e) and (f), the three-mode model relaxes to a Fock state with a fixed number of atoms occupying each mode, that is $\langle \psi (t)| {n}_{\mathrm{0,1,2}}| \psi (t)\rangle$ decays to ${\overline{n}}_{\mathrm{0,1,2}}$. The fluctuations after equilibration are small and decrease with g, as one sees by comparing figures 2(e) and (f).

4.2.2. Thermal equilibrium

After equilibration, the observable will have reached thermal equilibrium if its infinite-time average coincides with a thermodynamic average, that is if

$$\begin{eqnarray}&&{O}_{\mathrm{DE}}={O}_{\mathrm{ME}},\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{O}_{\mathrm{ME}}\equiv \displaystyle \frac{1}{{{ \mathcal N }}_{{E}_{\mathrm{ini}},\delta E}}\quad \displaystyle \sum _{\begin{array}{c}\nu \\ | {E}_{\mathrm{ini}}-{E}_{\nu }| \lt \delta E\end{array}}{O}_{\nu \nu }\end{eqnarray} \tag{ 17 }$$

is the average over a microcanonical ensemble. In the equation above, the sum in ν is performed over all energy eigenbases within the window δE around the energy E_ini of the initial state, that is $| {E}_{\mathrm{ini}}-{E}_{\nu }| \lt \delta E$, and ${{ \mathcal N }}_{{E}_{\mathrm{ini}},\delta E}$ is the total number of these eigenstates. Equation (16) holds when O_νν for eigenstates close in energy coincide with the microcanonical average, an idea that is at the heart of statistical mechanics and has become known as ETH.

When studying thermalization in finite systems, we investigate how close the left and right sides of equation (16) are and whether they approach each other as the system size increases. This is guaranteed to happen when the eigenstates are nearly random vectors. All random vectors are equivalent, since their components are simply random numbers. Thus, O_νν computed with one random vector is very similar to the result for any other random vector, apart from small fluctuations that decrease with system size.

In full random matrices, all eigenstates are random vectors, in which case thermalization is trivial. In realistic systems, the eigenstates away from the borders of the spectrum approach random vectors as the system moves toward the chaotic regime (see the discussion about the results for PR when g = 80 in figure 7). This is paralleled by the behavior of the EEV for the occupations of the three modes depicted in figure 8. As g increases, the fluctuations decrease and the EEVs show a smoother behavior with energy, especially close to the middle of the spectrum.

Zoom In Zoom Out Reset image size

To quantify the proximity of the EEV to the microcanonical average, we compute [74]

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{O}^{\mathrm{ME}}=\displaystyle \frac{{\sum }_{\nu }| {O}_{\nu \nu }-{O}_{\mathrm{ME}}| }{{\sum }_{\nu }{O}_{\nu \nu }},\end{eqnarray} \tag{ 18 }$$

where for the three-mode model, $O=\langle {n}_{0}\rangle ,\langle {n}_{1}\rangle ,\langle {n}_{2}\rangle$. The sum includes only the eigenstates within the microcanonical window [E − δE, E + δE]. In figure 9, we choose E very close to the middle of the spectrum and δE = 0.5, so that the microcanonical window contains approximately 10² levels. Provided there is a reasonable number of levels inside the window, the precise value of δE does not affect the results. Similarly to what we find for the chaos indicators in figure 6, ${{\rm{\Delta }}}_{O}^{\mathrm{ME}}$ in figures 9(a) and (b) decreases with g and also with N, suggesting that the fluctuations vanish in the thermodynamic limit.

Zoom In Zoom Out Reset image size

A more stringent demonstration of the vanishing of the fluctuations for strong interactions and large numbers of particles is made with the normalized extremal fluctuation of O, defined as [75],

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{O}^{\max -\min }=\left|\displaystyle \frac{\max O-\min O}{{O}_{\mathrm{ME}}}\right|.\end{eqnarray} \tag{ 19 }$$

The maximum ($\max O$) and minimum ($\min O$) values of the EEV are obtained for the eigenstates within the microcanonical window. The results are shown in figures 9(c) and (d) and mirror those from figures 9(a) and (b).

In figure 9, our choice of the window of energy in the middle of the spectrum implies infinite temperature. Studies of the dependence of the size of the fluctuations on temperature can also be done [75]. The fluctuations are expected to decrease as the temperature increases.

The small fluctuations of the EEV, which happens for chaotic eigenstates, are strong indications that equation (16) should hold. But for this to be indeed the case, the initial state needs to probe those chaotic states. We can then single out conditions that guarantee the onset of thermalization: the initial state is highly delocalized, so that equilibration can take place; the initial state has significant overlaps with chaotic eigenstates, that is E_ini falls within the chaotic region of the spectrum; and the width of the energy distribution of the initial state is smaller than or equal to the microcanonical window δE [14].

In figure 10, we finally compare the infinite-time average for the initial states chosen according to equation (11) with the microcanonical average. For this, we compute the relative difference,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{O}^{\mathrm{DE}-\mathrm{ME}}=\left|\displaystyle \frac{{O}_{\mathrm{DE}}-{O}_{\mathrm{ME}}}{{O}_{\mathrm{DE}}}\right|.\end{eqnarray} \tag{ 20 }$$

The two averages get indeed closer as g and N increase, confirming our expectations that thermalization should take place.

Zoom In Zoom Out Reset image size

In addition to strong interactions and large numbers of atoms, the energy E_ini of the initial state also plays a role in pushing the system toward thermal equilibrium. The vertical lines in figures 7 and 8 mark the position of E_ini. One sees that it moves closer to the middle of the spectrum as g increases. This further contributes to the viability of thermalization. Theoretically, we could also study the dependence of ${{\rm{\Delta }}}_{O}^{\mathrm{DE}-\mathrm{ME}}$ on the energy of the initial state for fixed g's and N's [77]. Experimentally, we are restricted to the initial states that can be actually prepared.

We chose not to show the results for ${{\rm{\Delta }}}_{\langle {n}_{2}\rangle }^{\mathrm{DE}-\mathrm{ME}}$ in figure 10. For the selected initial state only modes 0 and 1 are initially populated, so when g is small, the discrepancy between $\langle {n}_{2}{\rangle }_{\mathrm{DE}}$ and $\langle {n}_{2}{\rangle }_{\mathrm{ME}}$ is very large. However, the difference decreases rapidly as the interaction increases and shows results similar to those for $\langle {n}_{0}\rangle$ and $\langle {n}_{1}\rangle$ when g > 20.

We present in appendix A the study of the quantum dynamics for the two-mode model. While this model is insufficient to describe the system, it is interesting to emphasize differences and similarities with the three-model model. We mention two points. (i) Similarly to the three-mode model, with two modes one also finds damping of the oscillations. (ii) Interestingly, with two-modes there is absence of a transition to the quantum chaos regime. In contrast, two-mode model exhibits an excited state quantum phase transition (ESQPT), as expected from the similarity with the Lipkin–Meshkov–Glick Hamiltonian.

4.3. Quantum chaos and damping: extrapolation to large $N$

We now have the tools to compare the emergence of the damping of the oscillations with the onset of quantum chaos. For this, we choose thresholds for the damping time and chaos indicator β. For each N, we find the values of g at which the damping is so strong that the damping time τ in figure 5 is smaller than 2. We use this convention to get a value for g_damping, which we defined in section 3.3. Using this convention, we get a set of values of g_damping as a function of N for the criterion used in figure 5(a) and another one for the criterion used in figure 5(b). The values of g_damping versus N are plotted in figure 11. For each N, we also obtain the value of g for which β in figure 6 is larger than 0.3. We call this g_chaos, as it indicates that the system has already moved away from the integrable point and is approaching the chaotic regime. The behavior of the curves for g versus N extracted from τ and from β is very similar: the larger the number of atoms is, the smaller the interaction needs to be for damping and chaos.

Zoom In Zoom Out Reset image size

We note, however, that damping does not require the onset of chaos, as characterized by a Wigner–Dyson distribution. Damping can take place provided we do not encounter an excessive number of degeneracies or commensurate phases. Quantum chaos is a stronger condition to guarantee that not only the system relaxes, but it also reaches an equilibrium described by the Gibbs ensemble. This is why we chose as threshold for the chaos indicator β > 0.3 instead of a value closer to 1.

Similarly to what we did in figure 5, by fitting a curve to each set of data in figure 11, we extrapolate our results to N = 700, which is the typical number of atoms in the experiments. This leads to a value of g ∼ 10. Both analysis performed here, based on the damping time and on the approach to chaos, show that the strong damping described by the quantum model takes place at larger g than the damping observed experimentally [40].

We have shown that the three-mode quantum many-body model is a minimal model to qualitatively describe both the atomic density distribution oscillations and their damping. This behavior is qualitatively similar to the one observed experimentally with a quasi-1D BEC prepared in a coherent superposition of its two lowest motional states [11, 38, 39]. This system is isolated, so it does not include a mechanism for damping through an environment. Yet, one can make a system-environment analogy by viewing the second excited mode, which is essential for the decay of the oscillations, as a minimal environment, and the ground and first excited modes as constituting the system.

To characterize the observed decay of the oscillations, we employed the exact diagonalization of the many-body Hamiltonian for a number N of atoms ranging from 40 to 220 and a range of the interaction strength g. We showed that the damping time decreases as gN increases. The model also undergoes a transition to the quantum chaos regime when g becomes sufficiently strong. This value decreases as N increases. A key finding of this paper is the link established between the decay of the oscillations, the loss of coherence (fragmentation), and the approach to chaos.

The extrapolation of our results to the smallest number of atoms considered in the experiments (N = 700) reveals that, despite qualitatively reproducing the decay of the oscillations, the many-body three-mode model predicts damping times that are larger than those observed experimentally. We conjecture that this may be due to the fact that the experimental system is not a true quasi-1D system, but a cigar shaped condensate. For large interactions, phenomena occurring in the elongated direction may be the cause of an extra damping mechanism, which makes the damping time shorter. Whether this mechanism is the twin-atom generation processs described in [78] is out of the scope of this paper and a question to be investigated as an outlook.

The three-mode model offers a good example for studies of relaxation and thermalization in isolated quantum many-body systems. We have numerically shown that thermalization can indeed take place as g increases. The viability of thermalization is tightly connected with the onset of chaos. We expect that similar results can be found in other three-mode many-body models, as e.g. three bosonic species with coherent couplings in a trap or ultracold atoms in three-wells as in [79]. The role of the interaction energies that lead to the transfers between modes in our system would be played by the coherent coupling between species in the first model and by the tunneling energies between wells in the second one. The initial condition in these cases would be a coherent superposition of two of the species for the first model and two of the wells for the second one.

As a final remark, we mention a new study [80] about the conditions required to prepare an initial state (in general a Hamiltonian protocol) that does not equilibrate, thus introducing the concept of resilience against equilibration. This suggests a link between the area of nonequilibrium quantum dynamics and that of quantum resource theory. In the system studied here, an interesting outlook would be to study the resilience of possible initial coherent states.

We thank encouraging and fruitful discussions with B Julia-Diaz, T Berrada, C Gogolin, P Grzybowski and J F Schaft. We acknowledge the Spanish Ministry MINECO (National Plan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No. SEV-2015-0522, FPI), European Social Fund, Fundació Cellex, Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341 and CERCA/Program), ERC AdG OSYRIS, ERC advanced Grant QuantumRelax, EU FETPRO QUIC, and the National Science Centre, Poland-Symfonia Grant No. 2016/20/W/ST4/00314. MB was supported by the EU through the Marie Sklodowska Curie grant ETAB (ga no.656530) and by the FWF through the Lise Meitner grant CoPaNeq (M2088-M27). LFS was supported by the NSF Grant No. DMR-1603418.

If only two modes are considered, we approximate the field operator $\hat{{\rm{\Psi }}}$ describing the condensate by

$$\begin{eqnarray}&&\hat{{\rm{\Psi }}}\simeq {a}_{0}{\psi }_{0}+{a}_{1}{\psi }_{1},\end{eqnarray} \tag{ A.1 }$$

where the ψ_i are the two lower-lying eigenstates of the non-interacting part of the Hamiltonian (taken to be real and normalized to $\int | {\psi }_{i}{| }^{2}{\rm{d}}y=1$) and the ${\hat{a}}_{i}$ are annihilation operators associated with the modes, fulfilling the commutation relation $[{\hat{a}}_{i},{\hat{a}}_{j}^{\dagger }]={\delta }_{{ij}}$. Following the approach of [53], we obtain the effective two-mode Hamiltonian

$$\begin{eqnarray}&&{\hat{H}}_{2{\rm{m}}}=\displaystyle \frac{{\rm{\Delta }}E}{2}\,({a}_{1}^{\dagger }{a}_{1}-{a}_{0}^{\dagger }{a}_{0})+\displaystyle \frac{U}{4}{({a}_{1}^{\dagger }{a}_{1}-{a}_{0}^{\dagger }{a}_{0})}^{2}+{U}_{01}{({a}_{0}{a}_{1}^{\dagger }+{a}_{0}^{\dagger }{a}_{1})}^{2},\end{eqnarray} \tag{ A.2 }$$

with

$$\begin{eqnarray}&&{\rm{\Delta }}E={E}_{01}-(N-1)({U}_{00}-{U}_{11}),\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&U={U}_{00}+{U}_{11}-2{U}_{01},\quad \mathrm{and}\quad {U}_{{ij}}=\displaystyle \frac{g}{2}\int {\rm{d}}{y}| {\psi }_{i}{| }^{2}| {\psi }_{j}{| }^{2}.\end{eqnarray} \tag{ A.4 }$$

To connect with conventional approaches, let us introduce the operators ${J}_{x}=({a}_{0}{a}_{1}^{\dagger }+{a}_{0}^{\dagger }{a}_{1})/2$, ${J}_{y}=({a}_{0}{a}_{1}^{\dagger }-{a}_{0}^{\dagger }{a}_{1})/2{\rm{i}}$ and ${J}_{z}=({a}_{1}^{\dagger }{a}_{1}-{a}_{0}^{\dagger }{a}_{0})/2$, which satisfy angular momentum commutation relations. We can write equation (A.2) in the spin representation

$$\begin{eqnarray}&&{\hat{H}}_{2{\rm{m}}}={\rm{\Delta }}E\,{J}_{z}+{{UJ}}_{z}^{2}+4{U}_{01}{J}_{x}^{2}.\end{eqnarray} \tag{ A.5 }$$

In this way we show that this Hamiltonian resembles the bosonic Josephson Hamiltonian, with an additional energy offset between the two modes. The many-body dynamics and damping of the oscillations described with the two-mode model is very different from that described with three-mode model (see figures 2–4). In figure A1 we show the evolution of the mode amplitudes for the two modes and the off-diagonal correlations $\langle {a}_{1}^{\dagger }{a}_{2}\rangle$ for N = 1000 atoms. The initial condition is an equally weighted coherent superposition of the atoms in the ground and first excited states, (N₀, N₁) = (500, 500) and all relative phases equal to zero. The first observation is that the fast oscillation observed in the three-mode mode is the only one present in the two-mode model. But more importantly, as g is increased, the two-mode model also shows damping of the oscillations. This damping is qualitatively different from that observed in the three-mode case and in the experiment, as observed from figure A2, where the corresponding densities are depicted. First, the final state is different. It is also a fragmented state, but only over the two modes considered. For N = 1000 we observe that the very quick damping (damping time smaller than two oscillations) occurs also around g = 6.5, which is of the same order of the one observed for the three-mode model for N = 700. We note that, for g > 6.5, the off-diagonal elements $\langle {a}_{1}^{\dagger }{a}_{2}\rangle$ do not tend to zero anymore.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To better understand the two-mode model, we discuss its Hamiltonian, eigenvalues, and eigenstates. When U = 0, equation (A.2) represents the Lipkin–Meshkov–Glick (LMG) model [81], with the case of ${\text{}}U\ne 0$ being a generalization. This model is integrable and therefore presents no level repulsion. It is also known to exhibit an excited ESQPT.

In systems with a quantum phase transition, the gap between the ground state and the first excited state closes in the thermodynamic limit. In systems with an ESQPT [82, 83], this crossing occurs together with the clustering of the levels near the ground state and this divergence (peak) of the density of states moves to higher energies as the control parameter increases above the ground-state critical point. Concomitantly, the eigenstates that are very close to the energy of the ESQPT are highly localized leading to the slow evolution of initial states with similar energy [84].

The features of ESQPT for the Hamiltonian (A.2) are evident in figure A3. The top panels show results for the density of states, where two peaks are seen. They must be related with two different phase transitions caused by the three competing terms in equation (A.5). They emerge for g > 3 and are initially at the borders of the spectrum. We verified numerically that g > 3 is also the minimum value for which the two-mode model with N = 1000 shows damping within the longest simulations we performed (that is a time shorter than ∼30 oscillations in terms of ΔE₀₁/2π). As g increases, the two peaks approach each other (compare g = 10 and g = 20), merge together, and then separate again (compare g = 40 and g = 80). The peaks merge when only two main competing terms remain in equation (A.5).

Zoom In Zoom Out Reset image size

The bottom panels of figure A3 depict the results for the PR for all eigenstates as a function of energy. Dips in the PR occur at the same energies of the divergences of the density of states (see top and bottom panels of the figure). The dips indicate that the eigenstates around the energies of the ESQPT are very localized.

In summary, the two-mode model is significantly different from the three-mode model. Besides not being chaotic, it exhibits an ESQPT, which should affect the relaxation process.

Here, we discuss briefly the main ingredients of the body of theory which studies relaxation in isolated quantum systems (see e.g. [13, 17–27]) to highlight the connections with our discussion on quantum chaos. To this end, let us denote the evolving state through its density matrix ρ(t), with unitary dynamics dictated by a generic Hermitian Hamiltonian H, which is determined by its collection of eigenstates $\{| {\psi }_{\nu }\rangle \}$ and corresponding eigenenergies $\{{E}_{\nu }\}$. The Hilbert space is of finite dimension. Let us introduce also the dephased state as

$$\begin{eqnarray}&&\omega ({\rho }_{\mathrm{ini}})=\displaystyle \sum _{\nu }{p}_{\mathrm{ini}}^{(\nu )}| {\psi }_{\nu }\rangle \langle {\psi }_{\nu }| ,\end{eqnarray} \tag{ B.1 }$$

where ${p}_{\mathrm{ini}}^{(\nu )}=\langle {\psi }_{\nu }| {\rho }_{\mathrm{ini}}| {\psi }_{\nu }\rangle$ and ρ_ini is the initial state. When the latter is a pure state, ${p}_{\mathrm{ini}}^{(\nu )}=| {C}_{\mathrm{ini}}^{(\nu )}{| }^{2}$, as used in equations (14) and (15). According to equation (13), the PR_ini for the chosen initial state projected in the energy eigenbasis is given by

$$\begin{eqnarray}&&{\mathrm{PR}}_{\mathrm{ini}}=\displaystyle \frac{1}{{\sum }_{\nu }| {p}_{\mathrm{ini}}^{(\nu )}{| }^{2}}=\mathrm{Tr}(\omega {({\rho }_{\mathrm{ini}})}^{2}).\end{eqnarray} \tag{ B.2 }$$

If the system relaxes to equilibrium, its long-time average agrees with equation (B.1), that is

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{T\to \infty }\displaystyle \frac{1}{T}{\int }_{0}^{T}{\rm{d}}{t}\rho (t)=\omega ({\rho }_{\mathrm{ini}}).\end{eqnarray} \tag{ B.3 }$$

Equivalently, the expectation value of an arbitrary observable O tends to

$$\begin{eqnarray}&&\overline{O}=\mathrm{Tr}(\omega ({\rho }_{\mathrm{ini}})O),\end{eqnarray} \tag{ B.4 }$$

which is the expectation value in the dephased state (see also equation (15)). As explained in the main text below equation (15), equilibration requires that the temporal fluctuations around $\overline{O}$ be small and decrease with system size. Under the condition of lack of degeneracies, more precisely absence (or a negligible number) of degenerate level spacings [19, 22, 23], it has been shown that the variance of the temporal fluctuations is bounded by

$$\begin{eqnarray}&&\mathrm{Var}(O\rho ):= \overline{\mathrm{Tr}[O(\rho -\omega ({\rho }_{\mathrm{ini}}))]}\leqslant \displaystyle \frac{\parallel O{\parallel }^{2}}{{\mathrm{PR}}_{\mathrm{ini}}}.\end{eqnarray} \tag{ B.5 }$$

This means that for a given observable and under the condition mentioned above, relaxation occurs for highly delocalized initial states. In the context of quench dynamics, highly delocalized initial states emerge in systems perturbed far from equilibrium and where most eigenstates are strongly delocalized. These conditions are fulfilled by both chaotic and also interacting integrable models, as shown numerically in [26]. This justifies the sentence from the main text: 'equilibration does not require chaos in the sense of level repulsion, but it needs highly delocalized eigenstates, delocalized initial states, and not too many degeneracies.'

In [22] and others that followed, PR_ini has been named effective dimension, d_eff(ρ_ini), as one understands that is the actual dimension used by the initial state to relax to equilibrium, in contrast with the real dimension of the Hilbert state. If the effective dimension is proportional to the dimension of the Hilbert space, ${ \mathcal D }$, as it is often the case in chaotic systems, one expects that the initial state will thermalize after evolution.

In connection with the discussion presented here, we note that, recently, the phenomena of equilibration and the time scales required to equilibrate have been related to the quantum phenomena of dephasing in [47]. In this reference, the authors also estimate the equilibration time scale as roughly the inverse of the dispersion of the relevant energy gaps. As an outlook we find that such ideas can be investigated with the three-mode model.