Chaos induced breakdown of Bose-Hubbard modeling

TL;DR

This paper demonstrates that chaos can cause the breakdown of the Bose-Hubbard model's validity, even under conditions previously thought to be safe, by linking it to classical chaos theory and numerical observations.

Contribution

It introduces a simple model that explains the chaos-induced failure of the Bose-Hubbard approximation, connecting many-body effects with classical chaos analysis.

Findings

Bose-Hubbard approximation fails due to chaos even when detuned.

The model reproduces previous numerical breakdown observations.

Chaos enhancement explains the breakdown in many-body systems.

Abstract

We show that the Bose-Hubbard approximation fails due to the emergence of chaos, even when excited modes are far detuned and the standard validity condition is satisfied. This is formally identical to the Melnikov-Arnold analysis of the stochastic pump model. Previous numerical observations of Bose-Hubbard breakdown are precisely reproduced by our simple model and can be attributed to many body enhancement of chaos.

Figures

Equations

$${\hat{H}}_{\text{BH}}=\frac{U}{2}\sum _{i=1}^M{\hat{n}}_i({\hat{n}}_i-1)-\frac{K}{2}\sum _{⟨ij⟩}^M({\hat{a}}_i^{†}{\hat{a}}_j+{\hat{a}}_j^{†}{\hat{a}}_i)$$

$$u\ll \frac{\Omega }{K},u\equiv \frac{NU}{K}.$$

$$\hat{\mathcal{H}}={\hat{H}}_{\text{BH}}+\Omega {\hat{n}}_0+\frac{U}{2}{\hat{n}}_0({\hat{n}}_0-1)-\frac{\kappa }{2}\sum _{i=1}^2({\hat{a}}_i^{†}{\hat{a}}_0+\text{h.c.})$$

$$d(E,\Omega )=\frac{1}{T}{\int }_0^T|n(t)-n_{\mathrm{BH}}(t)|dt$$

$${\hat{H}}_{BH}$$

$$\hat{\psi }(x)={\Phi }_1(x){\hat{a}}_1+{\Phi }_2(x){\hat{a}}_2+{\Phi }_0(x){\hat{a}}_0$$

$$\kappa =-\int \left[\frac{{\hslash }^2}{2m}({\partial }_x{\Phi }_0{\partial }_x{\Phi }_{1,2})+{\Phi }_{0}V(x){\Phi }_{1,2}\right]dx,$$

$$\Omega =E_0-E_{1,2},$$

$$E_{1,2,0}=-\int \left[\frac{{\hslash }^2}{2m}|{\partial }_x{\Phi }_{1,2,0}{|}^2+V(x)|{\Phi }_{1,2,0}{|}^2\right]dx.$$

$$i{\dot{\hat{a}}}_1$$

$$i{\dot{\hat{a}}}_2$$

$$i{\dot{\hat{a}}}_0$$

$$i{\dot{a}}_1$$

$$i{\dot{a}}_2$$

$$i{\dot{a}}_0$$

$$|n,\varphi ,n_0,{\varphi }_0⟩=\frac{1}{\sqrt{N!}}{\left(\sum _{i=0}^2{\alpha }_i{\hat{a}}_i^{†}\right)}^N|\text{Vac}⟩,$$

Full text

I Introduction

The Bose-Hubbard model (BHM) Fisher89 ; Leggett01 is one of the most prominent tools in the study of interacting many-body systems, no less significant than its celebrated fermionic counterpart Hubbard63 ; Gutzwiller63 . Its use for treating ultracold bosons in shallow optical lattices Bloch05 ; Krutitsky16 , resulted in landmark experimental demonstrations of the superfluid to Mott insulator quantum phase transition Jaksch98 ; Orzel01 ; Greiner02a and of dynamical quantum phase revivals Greiner02b ; Will10 . These experiments and the subsequent development of quantum engineering techniques provide an unique opportunity to simulate various models and to explore the exotic phases of quantum matter using ultracold atoms Lewenstein07 ; Bloch12 .

The simplest BHM includes only two modes. Originally introduced in nuclear physics by Lipkin, Meshkov, and Glick LMG , two-mode theories apply equally well to Bose-Einstein condensates (BECs) of atoms in the same hyperfine state confined in a double-well trap, and to spinor BECs of atoms with two hyperfine states confined in a single-well trap. Two-mode BHMs are used to describe schemes for squeezing and entanglement Jo07 ; Esteve08 ; Riedel10 , Josephson oscillations and self-trapping in bosonic Josephson junctions Milburn97 ; Smerzi97 ; Albiez05 ; Gati07 ; Levy07 , the dynamical growth of quantum fluctuations Vardi01 ; Anglin01 ; Khodorkovsky08 ; Boukobza09 ; Boukobza10 ; Chuchem10 , and persistent currents and phase slips in superfluid atom circuits Wright13 ; Eckel14 ; Fialko12 ; Arwas16 ; Arwas17a ; Arwas17b .

Efforts towards realization of quantum simulators assume the validity of the BHM for describing cold atoms in optical lattices. The determination of validity criteria for the BHM is therefore of utmost importance. While Bose-Hubbard dynamics can be affected by coupling to the ever present excited Bloch bands Cederbaum ; Bergeman06 ; Trujillo-Martinez09 ; Smerzi11 ; Modugno18 ; Cataldo13 , it is generally accepted that the motion can be restricted to the lowest band of the site-chain, provided that the interaction energy is small with respect to the gap between it and the first excited band. Leggett01 ; Gati07 ; Milburn97 ; Anglin01 .

Despite this convention, the BHM was numerically shown to fail within its expected validity regime Cederbaum . Comparing multiconfigurational time-dependent Hartree for bosons (MCTDHB) simulations of the exact dynamics, with two-mode Bose-Hubbard dynamics, deviations from the BHM were observed even for tight traps where the standard validity condition is satisfied. The cause of this surprising failure remained unclear. No mechanism was offered to explain it or to attribute it to a particular dynamical feature, and no revised validity condition has been suggested.

Here, we show that even if its standard validity condition is satisfied, the BHM is still prone to failure due to the emergence of dynamical chaos. The naive assumption that high lying orbitals that do not participate in the dynamics merely renormalize the hopping elements via virtual transitions does not generally apply. Rather, chaos can be induced via the Melnikov-Arnold mechanism, as in Chirikov’s stochastic pump model Chirikov . Specifically, near-separatrix dynamics becomes stochastic due to the coupling with a high lying orbital. Beyond this zero order resonance, higher order resonances also show up and affect the dynamics. Consequently the lower band becomes entangled with higher bands, as reflected in reduced subsystem entropy measures. The previous observations of BHM breakdown Cederbaum can thus be attributed to emergence of chaos and generation of entanglement, and are fully reproduced by a simple three mode BHM.

II The BHM

II.1 Validity criteria

The experimental parameters of a linear 1D chain are the axial trap frequency $\omega_{\parallel}$, the barrier transmission coefficient $T$, and the atom number $N$. The 1D interaction strength is $\lambda_{0}=2\hbar\omega_{\perp}a_{s}$, where $a_{s}$ is the $s$-wave scattering length, and $\omega_{\perp}$ is the transverse trap frequency. The atoms mass is $m$. These parameters define three characteristic length-scales: the axial trap size $L=\sqrt{\hbar/m\omega_{\parallel}}$, the healing length $l_{c}=\sqrt{\hbar/2m\lambda_{0}n}$, and the mean distance between atoms ${d=1/n}$, where ${n=N/(2L)}$ is the average atom density. The standard BHM validity criterion assumes that the interaction energy is too small to bridge the $\Omega=\hbar\omega_{\parallel}$ gap between the lowest Bloch band and the first excited band. In terms of characteristic lengths, this means $\nu\ll 1$, where ${\nu\equiv(L/l_{c})^{2}=\lambda_{0}n/(\hbar\omega_{\parallel})}$. The $M$-site chain is then described by the tight-binding Hamiltonian,

[TABLE]

where $\hat{a}_{i}$ and $\hat{n}_{i}$ are the bosonic annihilation and number operators, $\langle ij\rangle$ denotes summation over nearest neighbors, and the effective parameters are the tunnel-splitting $K\approx\hbar\omega\sqrt{T}$, and the interaction-strength $UN\approx\lambda_{0}n$. The standard BHM validity criterion ${\nu\ll 1}$ then takes the form,

[TABLE]

where $u$ is the dimensionless interaction parameter.

II.2 Unexpected breakdown

In Ref. Cederbaum , strong deviations from the BHM were observed for $u~{}\sim 2$, even when $\nu=1/14$, leading to the conclusion: “Clearly, we have shown a failure of GP theory and the BH model within their range of expected validity”. Beyond the empirical determination that this unforeseen failure is “associated with fragmentation and correlations not captured by the standard theories”, no concrete mechanism has been suggested so far, for its origin.

Below, we establish that the underlying mechanism for the observed failure of the BHM is related to the robustness of chaos in the vicinity of separatrix regions. In order to investigate the chaos-induced breakdown of the BHM, we consider the $M=2$ (dimer) case, and emulate the effect of excited Bloch bands by adding to Eq.​​ (1) a single detuned bosonic mode:

[TABLE]

where $\Omega$ and $\kappa$ are respectively, the detuning and coupling of the auxiliary bosonic mode ‘0’. For an actual double well potential the above parameters including $\Omega$ and $\kappa$ can be estimated from the first principles as explaiend in Appendix A. In our calculation we set $\hbar=1$ and determine the units of time such that the hopping frequency is ${K=1}$.

III Dynamics

III.1 Quantum, classical, and semiclassical propagation

The motional constants of the Hamiltonian Eq.​​ (3) are the total three-mode energy $E$ and the total three-mode particle number $N$. In what follows, we distinguish between three types of propagation: (a) Quantum dynamics refers to the full many-body dynamics, i.e. the time-evolved many-body state is computed as $|\Psi(t)\rangle=e^{-i\hat{\mathcal{H}}t}|\Psi(0)\rangle$, for a given initial state $|\Psi(0)\rangle$; (b) Classical dynamics refers to the large $N$ Gross-Pitaevskii mean-field theory, in which the field operators $\hat{a}_{i}$ are replaced by classical $c$-numbers $a_{i}=\sqrt{n_{i}}\exp(i\phi_{i})$. Using $N$ conservation to eliminate the overall phase, the classical canonical variables are the dimer’s population imbalance ${n=n_{1}-n_{2}}$, the relative phase ${\phi=\phi_{1}-\phi_{2}}$, and the auxiliary mode’s population $n_{0}$ and phase $\phi_{0}$; (c) Semiclassical dynamics refers to the truncated Wigner classical propagation of a cloud of classical points that emulates the phase space distribution of the quantum state and averaging it to obtain the pertinent obsevables. For further detail on the dynamical equations, see Appendix B.

III.2 Semiclassical perspective

In Fig.​​ 1(a), we plot the pendulum-like classical phase space of the dimer Albiez05 ; Gati07 ; Chuchem10 , where $\kappa=0$. Since the isolated dimer has just one degree of freedom, its motion is necessarily integrable and in fact solvable by the algebraic Bethe ansatz Guan03 . The dimensionless interaction parameter $u$ distinguishes between three interaction regimes Leggett01 ; Gati07 ; Chuchem10 : Rabi ($u{<}1$); Josephson ($1{<}u{<}N^{2}$); and Fock ($u{>}N^{2}$). Within the Josephson interaction regime, the phase-space consists of a low energy Rabi-Josephson oscillation region Milburn97 ; Levy07 and a high-energy self-trapping region Smerzi97 ; Albiez05 ; Gati07 , that are separated by a mid-energy separatrix Chuchem10 .

In panels Fig.​​ 1(b)-(f), we set the coupling to $\kappa=0.5$ and study the effect of the third mode for various values of its detuning $\Omega$. It should be noticed that different trajectories do not have the same energy $E$, hence they do not belong to the same Poincare section. The addition of a third mode opens the way to non-integrable motion Korsch09 ; Hiller09 ; Liu07 ; Jung06 ; Penna03 ; Tikhonenkov13 ; Gallemi15 ; Dey18 ; Burkle19 resulting in stochastic regions in phase space due to non-linear resonances. Trajectories within these stochastic regions exhibit sensitive dependence on initial conditions, that can be quantified in terms of Lyapunov exponents. However, below we do not focus on the quantification of this sensitivity, but rather on the consequences of the chaotic motion.

Consider first the low $\Omega$ panels Fig.​​ 1(b)-(c). In this regime Eq.​​ (2) is not satisfied, and the trajectories are strongly affected by the coupling to the excited mode. Not only the separatrix motion is affected: in panel (b) most Rabi-Josephson trajectories become chaotic; while in panel (c) chaotic motion is obtained mainly for self-trapped trajectories.

By contrast, for the large $\Omega$ panels Fig.​​ 1(d)-(f), the validity criterion of Eq.​​ (2) is satisfied, and one expects the far-detuned mode to have a negligible effect. Indeed, this seems to be the case in most regions of phase space. However, a stochastic strip remains in the vicinity of the separatrix, even when the detuning is large. This is expected from the standard theory of non-linear resonances, and can be treated using the same procedure used for the Melnikov-Arnold analysis of the stochastic pump model Chirikov . Thus, due to chaos, the BHM breaks down for near-separatrix motion, and in the vicinity of resonances, irrespective of the standard validity condition.

III.3 Deviation measure

The deviation from the BHM is quantified by the instantaneous difference between the population imbalance $n(t)$ obtained by including the extra mode, and the same quantity $n_{\mathrm{BH}}(t)$ obtained in the two-mode approximation, averaged over $T=2\pi$,

[TABLE]

In Fig.​​ 2(a) we plot $d$ as function of the trajectory’s energy $E$ and the detuning $\Omega$. Since all simulations are launched with $n_{0}(0)=0$, the total energy $E$ equals the initial dimer energy. One observes that the deviation $d$ is large for trajectories in regions of stochastic motion. The corresponding quantum results are displayed in Fig.​​ 2(b). Initial states $|\Psi(0)\rangle$ in these quantum simulations were three-mode coherent states $|n,\phi,n_{0},\phi_{0}\rangle$, where $n,\phi,n_{0},\phi_{0}$ are the same as the initial values of the classical parameters (see Appendix B). We observe good quantum-classical agreement, with some blurring of classical features due to the finite uncertainty width of the initial coherent state.

IV Many body enhancement of chaos

The emergence of stochastic regions in Fig.​​ 1 and the deviation depicted in Fig.​​ 2 would have been obtained also if the two-mode system is driven at frequency $\Omega$. However, as shown below, the many-body aspect of incorporating an auxiliary mode in the Hamiltonian, amounts to considerable enhancement of chaos with respect to the corresponding driven system.

Consider the representative quantum spectrum in Fig.​​ 3(a). Each point represents an exact eigenstate $\nu$ of the Hamiltonian in Eq.​​ (3), positioned horizontally according to $X_{\nu}=\left\langle\nu\middle|\hat{n}_{0}\middle|\nu\right\rangle$, vertically according to its energy $E_{\nu}$, and color coded by its participation number $M_{\nu}=\left(\sum_{m}|\left\langle m\middle|\nu\right\rangle|^{4}\right)^{-1}$, where $m$ labels the unperturbed eigenstates of the uncoupled (${\kappa=0}$) system. Large $M_{\nu}$ implies that many eigenstates are mixed due to the ${\kappa\neq 0}$ coupling.

The Floquet states of any driving scenario would only mix $\left|m\right\rangle$ states with the same $n_{0}$, i.e. driving corresponds to vertical mixing in Fig.​​ 3(a). By contrast, referring to the ${(X,E)}$ diagram of the three-mode spectrum, we see that chaos in the vertical direction induces mixing also in the horizontal direction. Thus, the possibility of back action by the two mode dynamics on the auxiliary mode, results in enhanced chaos.

Considering the form of the two-mode dynamical equations ${\dot{a}_{j}=...+i(\kappa/2)a_{0}}$ (see Appendix B), the corresponding c-number driving would be obtained by substituting ${a_{0}=\sqrt{A}\exp(-i\Omega t)}$, where $A$ is a free constant parameter. To reproduce the effect of the third mode, the effective drive intensity $A$ should reflect its occupation $n_{0}$. Since from inspection of Fig.​​ 3(a) it is clear that this occupation is larger in regions of chaos due to the many body mixing, the third mode’s effect corresponds to an amplified drive intensity in these regions, i.e. to many body enhancement of chaos. This observation is somewhat reminiscenct of the dynamical enhancement of small perturbations in the nuclear physics context Flambaum .

We may attempt to evaluate the effective $A$ from the time averaged occupation in the 3rd orbital $\bar{X}=\overline{n_{0}(t)}$. If chaos was not present, we could estimate $\bar{X}$ from the equation of motion ${\dot{a}_{0}\approx i(\kappa/2)(a_{1}+a_{2})}$ leading to ${\bar{X}\sim\kappa^{2}/\omega^{2}(E)}$, where $\omega(E)$ is the frequency of the unperturbed dimer oscillations. The comparison of this estimate with the actual three-mode result in Fig.​​ 3(b), demonstrates the enhancement of $\bar{X}$ in chaotic regions. The inset in Fig.​​ 3(b) shows the dynamics of a driving scenario where $A=\bar{X}$ is the numerically extracted actual value of $\bar{X}$. Comparing with Fig.​​ 1c we realize that even this procedure still underestimates the enhanced chaos.

V Fragmentation and entanglement

V.1 One particle purity

In Fig.​​ 4a we prepare the system in representative eigenstates $|m\rangle$, with vacant 3rd orbital ($n_{0}=0$), of the unperturbed (${\kappa=0}$) Hamiltonian Eq.​​ (3), and plot the time dependence of the one-particle purity $\mathcal{P}={\mathrm{T}r}\left(\rho_{\mathrm{s}p}^{2}\right)$, where $\rho_{\mathrm{s}p}=\langle{\hat{a}}_{i}^{\dagger}{\hat{a}}_{j}\rangle$ is the reduced single particle density matrix. The one-particle purity value lies in the range $1/3<\mathcal{P}<1$ and its inverse indicates the number of modes required to capture the dynamics. Thus $\mathcal{P}\approx 1$ indicates the validity of mean-field theory, wherein all particles occupy a single orbital, $1/2<\mathcal{P}<1$ indicates two-orbital dynamics (strictly speaking, this range will also be obtained if the pertinent two orbitals project onto the auxiliary mode), while $\mathcal{P}<1/2$ clearly indicates the breakdown of the two-mode approximation. The latter is observed if the dynamics is affected by chaos.

V.2 Entanglement entropy

The addition of an excited mode also implies that entanglement could be generated between the dimer modes and the high frequency mode. The entanglement entropy is defined as $S=\text{Tr}(\rho_{\mathrm{d}}\ln\rho_{\mathrm{d}})=\text{Tr}(\rho_{0}\ln\rho_{0})$, where $\rho_{\mathrm{d}}=\text{Tr}_{0}(\rho)$ and $\rho_{0}=\text{Tr}_{\mathrm{d}}(\rho)$ are the many-body reduced density matrices of the dimer and the auxiliary mode, respectively. In Fig.​​ 4(a) we plot the time evolution of $S$ for a couple of representative simulations. Carrying out such simulations at various values of $\Omega$, for all dimer eigenstate preparations (distinguished by $E$) we plot the time-averaged $\overline{S(t)}$ in Fig.​​ 4(b). We observe that the entanglement entropy is large in region that support chaotic motion. Comparing with Fig.​​ 2 we see that chaos dominates at the vicinity of the separtrix: just below it for ${\Omega<1.5}$, and just above it for ${\Omega>1.5}$. Corresponding red regions below and above the separatrix energy in Fig.​​ 4(b) show the entanglement fingerprint of chaos.

VI Reconstruction of BHM breakdown

Returning to the results of Ref. Cederbaum , in Fig.​​ 5(a) we compare representative MCTDHB population dynamics with the quantum dynamics of our simple model (Eq.​​ (3)), for the same parameters and initial conditions. The BHM breakdown in Cederbaum is clearly captured by the $2+1$ mode model. The quantum results are reproduced by the semiclassical simulation in Fig.​​ 5(b), attained by propagating a cloud of classical trajectories from Fig.​​ 5(c)-(d). The observed failure of the BHM is therefore unequivocally related to the persistent separatrix chaos. In fact, as shown in Appendix C all the results of Ref. Cederbaum are reproduced by our model, including the thermalization of a self-trapped trajectory for large values of $u$ where the validity condition Eq.​​ (2) is violated and the two mode approximation breaks down. This thermalization too, is the result of chaotic ergodization of the type shown in Fig.​​ 1(c).

VII Summary

The description of optical lattices by the BHM should not be taken for granted. It is distinct from the truncation of a high-lying band of states in electronic systems that features a spectral gap Sowinski16 . In any $M$-site BHM the phase space is typically mixed, meaning that there are possibly vast regions of quasi-regular dynamics. Using the stochastic pump paradigm that leads to Arnold diffusion, we have shown that the effect of far-detuned modes requires to forego naive reasoning. The effect is amplified due to the many-body mixing of the eigenstates, and provides new insight for the relevance of the formal MCTDHB.

Acknowledgements

This research was supported by the Israel Science Foundation (Grant No. 283/18).

Appendix A Estimate of BHM parameters

The Hamiltonian describing the bosons in a double well trapping potential can be written as Smerzi97 ,

[TABLE]

where, $V(x)$ denotes the double well potential, $\lambda_{0}$ is the 1D interaction strength and $\hat{\psi}(x)$ is the bosonic field operator. Representing the wave functions, at the lowest quasi-degenerate orbitals (‘1’ and ‘2’) of the double well by $\Phi_{1}(x)$ and $\Phi_{2}(x)$ respectively, and at the excited orbital by $\Phi_{0}(x)$, we can expand $\hat{\psi}(x)$ as follows:

[TABLE]

where, $\hat{a}_{1,2,0}$ are the bosonic annihilation operators corresponding to the three modes ‘1’, ‘2’ and ‘0’ respectively. We define the two-mode parameters as in Smerzi97 and add the coupling parameters to the third mode :

[TABLE]

[TABLE]

where,

[TABLE]

Substituting Eq. 6 into Eq. 5 and using the definitions of the above parameters, we obtain the trimer Hamiltonian (3), where the single detuned bosonic mode ‘0’ is coupled to two modes (‘1’ and ‘2’) of the dimer.

Appendix B Dynamical equations

To study the dynamics governed by the Hamiltonian of Eq.​​ (3), we first write down the Heisenberg equations of motion $i\dot{\hat{a}}_{i}=[\hat{a}_{i},\hat{\mathcal{H}}]$, for the bosonic annihilation operators $\hat{a}_{i}$ :

[TABLE]

Two conserved quantities constrain the dynamics: One is the total energy and the other is the total number of particle resulting from $[\hat{\mathcal{H}},\hat{N}]=0$ where, $\hat{N}=\sum_{i=0}^{2}\hat{n}_{i}$.

Classical dynamics:

In the large $N$-limit the field operators can be replaced by $c$-numbers such as, $\hat{a}_{i}\rightarrow a_{i}$, $\hat{n}_{i}\rightarrow n_{i}$, where, $a_{i}=\sqrt{n_{i}}e^{i\phi_{i}}$, $n_{i}=|\alpha_{i}|^{2}$. The classical populations and phases $\{n_{i},\phi_{i}\}$ serve as conjugate dynamical variables. The equation of motion for the complex variables $a_{i}$ are thus given by,

[TABLE]

where, $u=UN/K$. Using $N$ conservation, we eliminate one degree of freedom, leaving $n=n_{1}-n_{2}$, $\phi=\phi_{1}-\phi_{2}$, $n_{0}$, and $\phi_{0}$ as the dynamical variables of our two-freedoms system. Evolving Eq.​​ (10), we obtain the classical results of Fig.​​ 1.

Quantum dynamics:

The quantum Fock states of the three-mode models are $|n,n_{0}\rangle$ where, $n_{0}\in[0,N]$ and $n\in[-(N-n_{0}),(N-n_{0})]$. The operation of the field operators is e.g. $\hat{a}_{1}^{\dagger}\hat{a}_{2}\,|n,n_{0}\rangle=\sqrt{n_{2}(n_{1}+1)}\,|n+2,n_{0}\rangle$ where, $n_{1}=(N-n_{0}+n)/2$ and $n_{2}=(N-n_{0}-n)/2$. Given some initial quantum preparation $|\Psi(0)\rangle$, we employ the Fock state representation to obtain the time evolved state $|\Psi(t)\rangle=e^{-i\hat{\mathcal{H}}t}|\Psi(0)\rangle$. In comparing quantum and classical dynamics, the initial quantum states $|\Psi(0)\rangle$ are three-mode coherent states,

[TABLE]

where $\alpha_{1}=\sqrt{(N+n)/2}$ and $\alpha_{2}=\sqrt{(N-n)/2}e^{i\phi}$ are the same as in the classical simulations and $\alpha_{0}=\sqrt{n_{0}}e^{i\phi_{0}}$ is set to zero.

Semiclassical dynamics:

The classical dynamics can only capture the quantum evolution up to the Ehrenfest time. However, averaging over many classical trajectories offers accuracy over longer timescales. In semiclassical simulations, we prepare a cloud of classical points that emulate the initial quantum phase-space distribution. Initial coherent states correspond to minimal Gaussian clouds with $2/N$ variance. Each point is then propagated classically and the values of observables are obtained by averaging over the cloud. This procedure is essentially a truncated Wigner phase-space approach.

Appendix C Comparison with previous numerical results

In Ref. Cederbaum , deviations from the two-mode BHM are observed numerically, by employing the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method, to obtain exact double-well dynamics that include the effect of excited modes. The dynamics shown in Fig. 1 of Cederbaum , that is reproduced here in Fig.​​ 6, demonstrate that the BHM breaks down even when $u$ is small with respect to $\Omega/K$ and Eq.​​ (2) is satisfied. The initial preparation with $n_{2}=N$ gives Rabi-Josephson oscillations for $u<2$ or self-trapped motion for $u>2$. For $u\approx 2$ it lies near the separatrix. Thus, Panels (a) and (b) illustrate the validity of the BHM for Rabi-Josephson oscillations, but panels (c) and (d) show its failure to depict near-separatrix motion.

The dashed red lines in Fig.​​ 6 corresponds to quantum propagation using the Hamiltonian of Eq.​​ (3). Our simple model clearly reproduces the results of the MCTDHB model (blue circles) to great accuracy. The relation between the quantum failure of the BHM and the classical persistence of chaos in the vicinity of the separatrix, is made clear by inspection of the semiclassical simulations presented in Fig.​​ 7. The dynamical instability of the classical trajectories that constitute the semiclassical cloud leads to strong deviations from the BHM in the vicinity of the separatrix, thereby reproducing the observed quantum result. We therefore deduce that the numerical observations in Ref. Cederbaum can be attributed in essence, to our mechanism of BHM failure due to chaos.

The observed agreement extends beyond the weak-interaction regime where the BHM validity condition is satisfied. In Fig. 2 of Ref. Cederbaum , the interaction strength is large enough to violate Eq.​​ (2). Thus, while the $n_{2}=N$ preparation is deep in the self trapped region in this case, so that almost no population oscillation exists in its two-mode dynamics, the excited mode can affect the entire phase-space, and in particular can transform self-trapped trajectories to chaotic ones, as in Fig.​​ 1(c). The result is the thermalization of the population distribution that was noted in Cederbaum . In Fig.​​ 8 we reproduce this result for the pertinent parameters, showing that it too is attributed to chaos.