Relaxation, chaos, and thermalization in a three-mode model of a BEC

We study the complex quantum dynamics of a system of many interacting atoms in an elongated anharmonic trap. The system is initially in a Bose-Einstein condensed state, well described by Thomas-Fermi profile in the elongated direction and the ground state in the transverse directions. After a sudden quench to a coherent superposition of the ground and lowest energy transverse modes, quantum dynamics starts. We describe this process employing a three-mode many-body model. The experimental realization of this system displays decaying oscillations of the atomic density distribution. While a mean-field description predicts perpetual oscillations of the atomic density distribution, our quantum many-body model exhibits a decay of the oscillations for sufficiently strong atomic interactions. We associate this decay with the fragmentation of the condensate during the evolution. The decay and fragmentation are also linked with the approach of the many-body model to the chaotic regime. The approach to chaos lifts degeneracies and increases the complexity of the eigenstates, enabling the relaxation to equilibrium and the onset of thermalization. We verify that the damping time and quantum signatures of chaos show similar dependences on the interaction strength and on the number of atoms.


Introduction
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,13,14,15,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 longtime averages. The conditions required for relaxation to occur have been discussed in [17,18,19,20,21,22,23,24,25,26,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,14,15,28] (ad-hoc as it is defined for the particular closed system).  Figure 1.
Illustration of the trapping potential and the Hamiltonian interaction parameters. U i represents the interaction coefficient between the atoms at level i = 0, 1, 2 and E ij = E j − E i is the energy difference between the subsequent levels i and j. The arrows between levels represent interaction or transfer processes: (i) those with an interaction energy U ij are interactions between atoms at level i and level j or transfer of two atoms from level i to j or vice versa; (ii) those with an interaction energy U kl ij destroy/create two atoms at level i and j and create/destroy one atom at level i and two atoms at level k and l. Note that all transfers between the levels are led by the interactions.
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 [38,39,11], 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, which was let to evolve 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,42,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 hundreds 87 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 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,14,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 Sec. 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 Sec. V. We also present the study of the quantum dynamics for the two-mode model in App. Appendix A and discuss some previous results about the conditions for relaxation in isolated quantum systems in App. Appendix B.

Three-Mode Many-Body Model
The second-quantized Hamiltonian for N ultracold bosons in an external potential V (x) is given by where Ψ(x) represents the field operator, x ∈ R 3 , 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 jx,j,jz = (1/2 + j x )ω x + (1/2 + j z )ω z + E j , with j x,z integers and E j the energy of j-th 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 01 = E 1 − E 0 , 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. Ψ = jx,j,jz a jx,j,jz ψ jx (x)ψ j (y)ψ jz (z), with a jx,j,jz 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 which produces the Hamiltonian where U ijkl = (g/2) dy ψ i ψ j ψ k ψ l , U ij = (g/2) dy ψ 2 i ψ 2 j , U i = (g/2) dy ψ 4 i (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. For example, the term a † 0 a † 2 (a 1 ) 2 cannot be neglected when one assumes that the lowest two modes are macroscopically occupied, a j and a † j are ∼ N j , j = 0, 1. In a BEC one assumes that fluctuations to other modes are negligible, and therefore neglects all quadratic terms with j > 0. But in this case, the aforementioned term is not quadratic, thus giving a relevant contribution to the Hamiltonian. In Appendix Appendix A we introduce the two-mode model, which we will compare with the results of the three-mode model obtained below. We mention here that the aforementioned reflection symmetry y → −y translates into [H 3m , P ] with P = (−1) n 1 , with n 1 = a † 1 a 1 . As we discuss later, identification of this symmetry is important in the study of quantum chaos below.
To be able to model the experiments in [38,39,11,40] we have to relax the onedimensional 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. It is still valid the assumption that the dynamics in the z direction is frozen, 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 case, as we discuss later, the period of the density oscillations that occur even in the noninteracting case is given by E 01 (see Sec. 3.1). On the other hand the excitations along x are much lower in energy, since ω x < E 01 . 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 Eq. (3) can be obtained, with a different expression of the coefficients. They have to be calculated taking into account that N 0 and N 1 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 exp ijkl as they are close to typical experimental values [38,39,11,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 × U exp ijkl . 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 .

Semiclassical dynamics
To get physical insight into the role of the different parameters of Eq. (3) and the processes represented in Fig. 1, we consider a semiclassical version of the above equations of motion. For this, we neglect quantum fluctuations and treat the operators as c- 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 with similar equations for N 2 and φ 2 , and In Fig. 2 (a) and Fig. 2 (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 0 , N 1 , N 2 ) = (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 fast ≈ 0.5 ms and a slow oscillation with a longer period T slow ≈ 5 ms. From inspection of Eqs. (8)-(11) 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 Fig. 3). Because of that we adimensionalize the time in all figures with ∆E 01 . For g small inspection of Eq. (8) shows that N 0 will oscillate slightly around their initial value with amplitude proportional to U 01 and period proportional to E 01 [see first lines in Eqs. (8) and (9)]. All other terms are small because N 2 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 2 is no longer negligible. In such case we observe a second type of oscillation, which is slower and has an amplitude proportional to U 0112 .
For comparison, we show in Figs. 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 Figs. 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 Eq. (8-11) all variables and parameters associated with mode 2. They are equivalent to the equations of an oscillator as the bosonic Josephson junction [51].
In Fig. 3 (a) we depict the evolution of the corresponding density |ψ(t, z)| 2 for the example shown in Fig. 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 Fig. 3 (b), which corresponds to the case in Fig. 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 apprehend the complexity of the dynamics. As discussed in Sec. 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 Fig. 2 (a) and Fig. 2 (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 non-linear dynamics in the semiclassical equations, as our goals relies in the quantum dynamics discussed in what follows.

Quantum many-body dynamics
To simulate the many-body dynamics, we perform the exact diagonalization of the many-body Hamiltonian in Eq. 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,53,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 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], where |vac is the vaccum and |c 0 | 2 + |c 1 | 2 = 1. The initial average population of the two first modes is n 0 (0) = |c 0 | 2 and n 1 (0) = |c 1 | 2 , with n i (t) = ψ(t)|a † i a i |ψ(t) . We present results for an initial coherent state with c 0 = c 1 , but no qualitative differences are observed when c 0 = 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 Fig. 2). This is not the case if one takes a Fock initial state, i.e. |ψ ini = 1/( . To circumvent our limitation to relatively small system sizes, we increase the effective interaction constant to reach values of the product g N that are close to those found experimentally [39,11,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 g N 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 g N , 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 Figs. 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 Figs. 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 01 /2π, with the longest simulations up to t max = 30 ∆E 01 /2π).

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 g N . 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, a † i a i for i = 0, 1, 2, are accompanied by a decay to zero of the coherence terms a † i a j with i = j. At initial times, there is only one large eigenvalue, as expected for condensation in the coherent superposition of the ground and first excited mode. In time, the second largest eigenvalue becomes also sizeable, indicating fragmentation. [The third eigenvalue, associated with the occupation of the second excited mode, is not shown, but it also becomes nonnegligible.] This effect occurs at shorter times for larger interaction strengths. Thus, the damping in Fig. 2 (e), Fig. 2 (f), and Fig. 3 (d), the loss of coherence in panels (a) and (b), and the fragmentation in (c) occur together.
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), ρ = |Ψ Ψ|, 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 Fig. 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 Fig. 4 (c), but the third one is also non-negligible, as the sum of the first two does not amount to 1]. This indicates fragmentation in the three modes.
The loss of coherence, just as damping, occurs earlier in time as gN is made larger [compare Fig. 2 (e), Fig. 2 (f), and Fig. 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 [62,49].
We estimate the damping time τ from the evolution of the eigenvalues of the OBDM. In Fig. 5, we plot our numerical estimates for τ as a function of g N . Two numerical criteria are used to determine the damping time. In Fig. 5 (a), τ is the time at which the largest eigenvalue of the OBDM gets smaller than 0.98. In Fig. 5 (b), τ is the time when the largest eigenvalue becomes smaller than 0.85. 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 which correspond to the curve expected for N = 700. The result is the thick red curve in Figs. 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 log 10 (gN ) = 3.6 [3.8] with N = 700 for the criterion used in Fig. 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 Figs. 5 (a) and (b) is of the order of thousands of oscillations, much larger than the one observed in the experiment (τ exp ∆E 01 ≤ 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.

Onset of Quantum Chaos and Damping
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 the equilibration is the internal couplings 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 App. 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.

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.

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, P P (s) = exp(−s), 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 In Fig. 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 Eq. (3), we separate the eigenvalues by the parity of the eigenstates. The two top rows of Fig. 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.
The second and third rows of Fig. 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 Fig. 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], where Γ is the 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 where s 0 is the first intersection point of P P (s) and P W D (s). For a Poisson distribution, η → 1, and for the Wigner-Dyson distribution, η → 0. In Fig. 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 W D , or any intermediate distribution between the two, as seen in the first row of Fig. 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 W D . We instead approach the integrable point of three uncoupled oscillators, 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 Fig. 6 reinforce our statement above that the transition to chaos happens for smaller values of g as N increases.

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,72,73]. To measure the level of delocalization of the eigenstates |ψ ν , one can use quantities such as the participation ratio, where C (ν) j = ϕ j |ψ ν is the overlap between the eigenstate |ψ ν and the basis vector |ϕ j . 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, |ϕ j = |N 0 , N 1 , N 2 . In the absence of interaction, when the eigenstates coincide with the basis vectors, PR=1.
Each panel of the two top rows of Fig. 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 Fig. 7 for comparison [74]. 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.
The middle row of Fig. 7 illustrates the consequence of the transition to chaos.  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,75,76], as discussed next.

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.

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,77]. A brief discussion about the subject is presented in App. 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 where 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,18,19,20,21,22,23,24,25,26,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 Figs. 2 (e) and (f), the three-mode model relaxes to a Fock state with a fixed number of atoms occupying each mode, that is ψ(t)|n 0,1,2 |ψ(t) decays to n 0,1,2 . The fluctuations after equilibration are small and decrease with g, as one sees by comparing Fig. 2 (e) and Fig. 2 (f).

Thermal Equilibrium
After equilibration, the observable will have reached thermal equilibrium if its infinite-time average coincides with a thermodynamic average, that is if where is the average over a microcanonical ensemble and N E ini ,δE is the number of energy eigenbases in the window δE taken around the energy E ini of the initial state. 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 eigenstate thermalization hypothesis (ETH). When studying thermalization in finite systems, we investigate how close the left and right sides of Eq. (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 Fig. 7). This is paralleled by the behavior of the EEV for the occupations of the three modes depicted in Fig. 8. As g increases, the fluctuations decrease and the EEVs show a smoother behavior with energy, especially close to the middle of the spectrum. To quantify the proximity of the EEV to the microcanonical average, we compute [75] where for the three-mode model, O = n 0 , n 1 , n 2 . The sum includes only the eigenstates within the microcanonical window [E − δE, E + δE]. In Fig. 9, we choose E very close to the middle of the spectrum and δE = 0.5, so that the microcanonical window contains approximately 10 2 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 Fig. 6, ∆ ME O in Figs. 9 (a) and (b) decreases with g and also with N , suggesting that the fluctuations vanish in the thermodynamic limit.
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 [76], The  In Fig. 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 [76]. 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 Eq. (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 Fig. 10, we finally compare the infinite-time average for the initial states chosen according to Eq. (11) with the microcanonical average. For this, we compute the relative difference, The two averages get indeed closer as g and N increase, confirming our expectations that thermalization should take place. 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 Fig. 7 and Fig. 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 ∆ DE-ME O on the energy of the initial state for fixed g's and N 's [78]. Experimentally, we are restricted to the initial states that can be actually prepared.
We chose not to show the results for ∆ DE-ME n 2 in Fig. 10. For the selected initial state only modes 0 and 1 are initially populated, so when g is small, the discrepancy between n 2 DE and n 2 ME is very large. However, the difference decreases rapidly as the interaction increases and shows results similar to those for n 0 and n 1 when g > 20.
We present in App. Appendix A the study of the quantum dynamics for the twomode 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 it similarity with the Lipkin-Meshkov-Glick Hamiltonian.

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 Fig. 5 is smaller than 2. We use this convention to get a value for g damping , which we defined in Sec. 3.3. Using this convention, we get a set of values of g damping as a function of N for the criterion used in Fig. 5 (a) and another one for the criterion used in Fig. 5 (b). The values of g damping vs N are plotted in Fig. 11. For each N , we also obtain the value of g for which β in Fig. 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 vs 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.  Figure 11. For each N , values of g at which the damping time is smaller than τ = 2 (which we name as g damping ) are shown with crosses for the criterion used in Fig. 5 (a) and with circles for the criterion in Fig. 5 (b)]. We also show as a function of N , values of g at which the chaos indicator β > 0.3, which we name as g chaos and represent with squares. The solid lines correspond to fittings to the numerical results. The three curves have the same qualitative behavior. The extrapolation to N = 700 gives g ∈ [8,14].
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 Fig. 5, by fitting a curve to each set of data in Fig. 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].

Conclusion
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 [38,39,11]. 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 g N 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 [79] 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 Ref. [80]. 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 [81] 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. evolution of the mode amplitudes for the two modes and the off-diagonal correlations a † 1 a 2 for N = 1000 atoms. The initial condition is an equally weighted coherent superposition of the atoms in the ground and first excited states, (N 0 , N 1 ) = (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 Fig. 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 a † 1 a 2 do not tend to zero anymore.  Figure A1. Evolution of the mode average occupations (left column) and the offdiagonal elements of the one body density matrix (OBDM) a † 1 a 2 (right column) for the two-mode many-body model for N = 1000 atoms (same initial state as in Fig. 2). On left column, we represent n 0,(1) with a red thin (blue thick) line. From the top to the bottom panel, the interaction is increased as g = 5, 6, 6.5 and 7.
To better understand the two-mode model, we discuss its Hamiltonian, eigenvalues, and eigenstates. When U = 0, Eq. (A.2) represents the Lipkin-Meshkov-Glick (LMG) model [82], with the case of U = 0 being a generalization. This model is integrable and therefore presents no level repulsion. It is also known to exhibit an excited state quantum phase transition (ESQPT).  Figure A2. Evolution of the density profiles for N = 1000 atoms for the two-mode many-body model (same initial state as in Fig. 2). (a) to (d) correspond to g = 5, 6, 6.5 and g = 7, respectively. Damping occurs for similar gN as in the three mode model, but qualitatively the final state is different both to that reached at large times in the three-mode model and in the experiment.
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 [83,84], 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 [85]. The features of ESQPT for the Hamiltonian (A.2) are evident in Fig. 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 Eq. (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 01 /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 Eq. (A.5).
The bottom panels of Fig. 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 (cf. 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.

Appendix B. Condition for relaxation
Here, we discuss briefly the main ingredients of the body of theory which studies relaxation in isolated quantum systems (see e.g. [13,17,18,19,20,21,22,23,24,25,26,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 {|ψ ν } and corresponding eigenenergies {E ν }. The Hilbert space is of finite dimension. Let us introduce also the dephased state as where p (ν) ini = ψ ν |ρ ini |ψ ν and ρ ini is the initial state. When the latter is a pure state, p (ν) ini = |C (ν) ini | 2 , as used in Eqs. (14) and (15). According to Eq. (13), the PR ini for the chosen initial state projected in the energy eigenbasis is given by If the system relaxes to equilibrium, its long-time average agrees with Eq. (B.1), that is which is the expectation value in the dephased state [see also Eq. (15)]. As explained in the main text below Eq. (15), equilibration requires that the temporal fluctuations around 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 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 Ref. [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, 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.