Chaos-induced breakdown of Bose-Hubbard modeling

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.

Figure 1

Classical Bose-Hubbard dynamics. Panel (a) is for an isolated two-mode system. It shows representative Rabi-Josephson (magenta star, black arrow), near-separatrix (green circle), and self-trapped (red square) trajectories for $u=3$ and $\kappa =0$. Panels (b–f) are for a two-mode system coupled to a third detuned mode ($\kappa =0.5$), with $\mathrm{\Omega }=0.5,2.0,4.5,6.0,7.0$, respectively. They show ${\varphi }_0=0$ Poincaré sections for the same initial conditions as in (a), while the third orbital is initially empty ($n_0=0$).

Figure 2

The deviation $d$ of Eq. (4), plotted as a function of the trajectory's energy $E$ and the detuning $\mathrm{\Omega }$: (a) classical simulations and (b) quantum simulations with $N=100$, launched at the corresponding coherent state. The parameters $u=3$ and $\kappa =0.5$ are the same as in Fig. 1. Dashed line marks the separatrix energy. The detuning values of Figs. 1–1 are marked on the horizontal axis, and the energies of the plotted trajectories are marked on the vertical axis with the same marker convention. Here and throughout the manuscript, we use dimensionless energy and timescales, e.g., $\mathrm{\Omega }\to \mathrm{\Omega }/K$, $E\to E/K$, $t\to Kt$.

Figure 3

(a) The eigenspectrum $\left\{E_{\nu }\right\}$ obtained from diagonalization of the Hamiltonian Eq. (3) for $N=100$, $u=3$, $\kappa =0.5$, $\mathrm{\Omega }=2$. Each point is positioned horizontally according to $X_{\nu }=\langle {\widehat{n}}_0\rangle$ and color coded by participation number $M_{\nu }$. Strong mixing is witnessed at the separatrix energy (dashed line) and extends to the region above it. (b) Time-averaged occupation of the third mode, $\overline{X}\equiv \overline{n_0\left(t\right)}$, obtained from three-mode classical dynamics (solid blue) as a function of $E$ compared with the estimate $\overline{X}\sim {\kappa }^2/{\omega }^2\left(E\right)$ (red dashed) that assumes quasi-integrable orbits. The separatrix energy is indicated by a vertical dashed line, while the energies of the trajectories in Fig. 1 are indicated by symbols. The inset shows the stroboscopic map for the two-mode BHM in the presence of driving with frequency $\mathrm{\Omega }=2$ and intensity $A=\sqrt{\overline{X}}$. It corresponds to Fig. 1 and has the same axes.

Figure 4

(a) Time evolution of the entanglement entropy $S$ and the single-particle purity $\mathcal{P}$ for a chaotic (dashed) and for a quasi-integrable (solid) three-mode dynamics, starting from eigenstates of the unperturbed two-mode system. (b) The time-averaged entanglement entropy $\overline{S\left(t\right)}$ as a function of $E$ and $\mathrm{\Omega }$. Markers indicate the parameter values for the curves in (a). The other parameters are as in Fig. 2.

Figure 5

Comparison with Ref. [35]. (a) Quantum dynamics of $n_2$ according to the BHM (solid black), the MCTDHB (blue circles), and our simple model (dashed red). MCTDHB data is extracted from Fig. 1(d) in Ref. [35]. (b) The corresponding semiclassical dynamics with the same color code. (c) Classical phase-space structure of the isolated two-mode BHM ($\kappa =0$). (d) Corresponding Poincaré sections at ${\varphi }_0=0$ in the presence of the excited mode with detuning $\mathrm{\Omega }=5$ and coupling $\kappa =0.75$.

Figure 6

Breakdown of the BHM despite weak interaction. Comparison of population dynamics obtained from quantum dynamics using the two-mode BHM (solid black line), the MCTDHB method of Ref. [35] (blue circles), and our $2+1$ mode model Eq. (3) (dashed red line). The parameters and the initial preparation with $n_2=N$ are identical to Fig. 1 of Ref. [35]: (a) $u=1.4$, $N=20$; (b) $u=1.35$, $N=100$; (c) $u=2.26$, $N=20$; and (d) $u=2.17$, $N=100$. The parameters used in our model are $\mathrm{\Omega }=5$ in all panels, $\kappa =0.65$ in (a, b) and 0.75 in (c, d).

Figure 7

Semiclassical dynamics. (a) Classical Gross-Pitaevskii simulations corresponding to Fig. 6 with the same color code. Representative trajectories launched in the vicinity of the classical initial conditions are shown as thin dotted lines. (b) Same for Fig. 6, with no additional trajectories. (c) Semiclassical dynamics, obtained from the classical propagation of an initially Gaussian cloud with width $\sqrt{2/N}$, reproduce the quantum results of Fig. 6. (d) Similarly, semiclassical propagation reproduces Fig. 6.

Figure 8

Chaotic ergodization at strong interaction. (a) Quantum population dynamics obtained by time propagation with the Hamiltonian of Eq. (3), with $\mathrm{\Omega }=30$ and $\kappa =40$, for the same parameters as in Fig. 2 of Ref. [35]: $u=43.4$, $N=100$. (b) Classical population dynamics for the same parameters using either the two-mode BHM (black line, depicting a nearly stationary point at $n_2=N$) or the Hamiltonian of Eq. (3) (red line, exploring the entire allowed population range). (c) The two-mode phase space. (d) Poincaré sections in the $2+1$-mode phase space, as in Fig. 1. The self-trapped dimer at $n_2=N$ becomes chaotic due to the coupling to the auxiliary mode and explores the entire stochastic band, resulting in the “thermalization” of the population distribution.