Many-body dynamical localization and thermalization

TL;DR

This paper demonstrates quantum dynamical localization in a system of two weakly-coupled chaotic Bose-Hubbard units, revealing localization effects in thermalization processes beyond classical chaos descriptions.

Contribution

The study introduces an intrinsic localization measure based on initial condition memory, applicable to high-dimensional quantum systems where classical correspondence breaks down.

Findings

Quantum localization observed in thermalization of coupled chaotic subsystems

Localization persists near the mobility edge, suppressing ergodization

Weak localization effects are identified despite high system dimensionality

Abstract

We show that a quantum dynamical localization effect can be observed in a generic thermalization process of two weakly-coupled chaotic subsystems. Specifically, our model consists of the minimal experimentally relevant subsystems that exhibit chaos, which are 3-site Bose-Hubbard units. Due to the high dimensionality of the composite 6-site system, the quantum localization effect is weak and can not be resolved merely by the breakdown of quantum-to-classical correspondence. Instead, we adopt an intrinsic definition of localization as the memory of initial conditions, that is not related to the underlying classical dynamics. We discuss the dynamics in the chaotic sea, and in the vicinity of the mobility edge, beyond which ergodization is suppressed.

Equations

$${\mathcal{H}}_{\alpha }=\frac{U}{2}\sum _{j=1}^3{\hat{n}}_{\alpha j}^2-\frac{K}{2}\sum _{j=1,3}\left({\hat{a}}_{\alpha j}^{†}{\hat{a}}_{\alpha 2}+\text{H.c.}\right).$$

$$\mathcal{H}$$

$${\mathcal{H}}_c$$

$$u$$

$$v$$

$${\mathcal{H}}^{cl}={\mathcal{H}}_L^{cl}+{\mathcal{H}}_R^{cl}+\frac{v}{4}\sum _{j=1}^3{I}_{Lj}{I}_{Rj},$$

$${\mathcal{H}}_{\alpha }^{cl}=\frac{u}{4}\sum _{j=1}^3{I}_{\alpha j}^2-\sum _{j=1,3}\sqrt{I_{\alpha j}{I}_{\alpha 2}}\cos ({\phi }_{\alpha j}-{\phi }_{\alpha 2}).$$

$$E$$

$$x$$

$$P_{\text{erg}}(x)\propto g(E_0+x/2)g(E_0-x/2).$$

$$P_t(E,x|E_0,x_0)=|⟨E,x\left|{\mathrm{e}}^{-i\mathcal{H}t}\right|E_0,x_0⟩{|}^2.$$

$$P_t(x)=\sum _{x<x^{\prime }<x+\delta x}{P}_t(E,x^{\prime }|E_0,x_0).$$

$$P_{\infty }(x)=\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{P}_t(x)dt.$$

$$P_t(X)=\sum _{n,m}⟨X|{\mathcal{E}}_n⟩⟨{\mathcal{E}}_n|X_0⟩⟨X_0|{\mathcal{E}}_m⟩⟨{\mathcal{E}}_m|X⟩{\mathrm{e}}^{i({\mathcal{E}}_m-{\mathcal{E}}_n)t}.$$

$$P_{\text{sat}}(X)=\sum _n|⟨X|{\mathcal{E}}_n⟩{|}^2|⟨{\mathcal{E}}_n|X_0⟩{|}^2.$$

$${\Delta }_t(x)=P_t(x)-P_{\text{erg}}(x)$$

$$f_t=\sum _{\Delta >0}{\Delta }_t(x)\in [0,1]$$

$${⟨|x|⟩}_{\infty }\equiv \sum _x|x|P_{\infty }(x)$$

$$D_E(x)=\frac{v^2}{8}\sum _{j,k=1}^3{\int }_{-\infty }^{\infty }{\omega }^2{S}_{jk}(\omega ,{\epsilon }_L)S_{kj}(\omega ,{\epsilon }_R)\frac{d\omega }{2\pi },$$

$$\displaystyle M\equiv\Big{[}\sum_{x}P(x)^{2}\Big{]}^{-1}~{}.$$

$$\begin{array}{rl}& |\text{sym}⟩=\frac{1}{2}(|1⟩+\sqrt{2}|2⟩+|3⟩),\\ & |\text{anti}⟩=\frac{1}{2}(|1⟩-\sqrt{2}|2⟩+|3⟩),\\ & |\text{dark}⟩=\frac{1}{\sqrt{2}}(|1⟩-|3⟩).\end{array}$$

$${\omega }_{\text{sym}}=-{\omega }_0,{\omega }_{\text{anti}}=+{\omega }_0,{\omega }_{\text{dark}}=0,$$

$$\begin{array}{rl}& {\hat{a}}_{\text{sym}}=\frac{1}{2}({\hat{a}}_1+\sqrt{2}{\hat{a}}_2+{\hat{a}}_3),\\ & {\hat{a}}_{\text{anti}}=\frac{1}{2}({\hat{a}}_1-\sqrt{2}{\hat{a}}_2+{\hat{a}}_3),\\ & {\hat{a}}_{\text{dark}}=\frac{1}{\sqrt{2}}({\hat{a}}_1-{\hat{a}}_3),\end{array}$$

$$H_{\text{trimer}}={\omega }_0({\hat{n}}_{\text{anti}}-{\hat{n}}_{\text{sym}}).$$

$${\rho }_{jk}^{(1)}=\frac{1}{N}⟨{\hat{a}}_j^{†}{\hat{a}}_k⟩.$$

$$\gamma =\text{purity}=\text{trace}\left\{[{\rho }^{(1)}{]}^2\right\}\in [1/3,1].$$

$$r_n=\frac{\text{min}\left\{s_n,s_{n+1}\right\}}{\text{max}\left\{s_n,s_{n+1}\right\}}\in [0,1]$$

$$Y_j(\omega )\equiv \underset{T\to \infty }{\lim }\frac{1}{\sqrt{T}}{\int }_0^T{I}_j(t)e^{-i\omega t}dt.$$

$$S_{jj}(\omega ;\epsilon )=|Y_j(\omega ){|}^2.$$

$$C_{jk}(\tau )\equiv ⟨I_j(t+\tau )I_k(t)⟩.$$

$$S_{jk}(\omega ;\epsilon )=Y_j^{\ast }(\omega )Y_k(\omega ),$$

$$\dot{x}=-\frac{v}{4}\sum _{j=1}^3[{\dot{I}}_{Rj}(t)I_{Lj}(t)-{\dot{I}}_{Lj}(t)I_{Rj}(t)].$$

Full text

I Introduction

The study of dynamical localization in a low-dimensional chaotic system was pioneered in the publication “Chaos, Quantum Recurrences, and Anderson Localization” by Fishman, Grempel and Prange Fishman , which had been motivated by the puzzling numerical observation of quantum-suppressed chaotic diffusion, by Casati, Chirikov, Izrailev, and Ford Casati . For an extended period, most research efforts focused on the quantum kicked-rotor system, regarded as the quantized version of the Chirikov standard map (see standard_map and references within). This system has one degree of freedom, and the driving (kicking) effectively adds an extra half freedom.

Dynamical localization.– In an Anderson tight-binding disordered chain all the eigenstates are exponentially localized, and therefore the spreading of a wave packet in space ($x$) is always suppressed after some break time ($t^{*}$). The break time is related to the localization length of the eigenstates. Similarly, in the kicked-rotor problem, the variable “$x$” is the angular momentum. While there is no disorder in this one+half degree of freedom dynamical system, the chaotic classical dynamics of “$x$” is stochastic-like, provided the kicks are strong enough. Thus, in analogy to Anderson localization, the diffusive spreading is suppressed in quantum simulations after some break time $t^{*}$, indicating that the Floquet eigenstates are localized. The relation between the break time and the localization length can be formalized Chirikov ; Shep ; Dittrich ; Cohen , and both are related to the disorder strength (in the Anderson model) or to the kicking strength (in the kicked-rotor system).

Localization in higher dimensions.– Strong Anderson localization is not necessarily present in higher dimensions. This is explained by scaling theory and renormalization group methods. In general, there is a mobility edge that separates regions where strong localization prevails from those where spreading is unaffected by quantization. Nevertheless, in the latter case it may be possible to obtain a weak-localization effect, meaning that the probability of return to the initial region is affected. Weak localization is thus a memory effect, expressed as a dependence of the outcome on initial conditions.

Localization in finite systems.– Most condensed matter physics literature concerns disordered infinite-volume systems, and the weak-localization effect is associated with time-reversal symmetry. The paradigmatic kicked rotor is also formally identical to an infinite-volume chain, because the range of angular momentum is unbounded, and the energy is not a constant of motion. Still, the notion of quantum localization has been extended also to the realm of undriven finite systems that have a finite-volume energy surface. In this context we can define localization as the absence of ergodicity. If the system is ergodic, an initial cloud of evolving points (“classical wave packet”) is expected to eventually be smeared uniformly over the entire energy surface, approaching a microcanonical-like distribution. Quantum localization can therefore be defined as the lack of ergodicity in quantum systems for which the corresponding classical dynamics is ergodic. This means that even after a very long time the quantum probability distribution is not microcanonically uniform. Heller Heller has realized that (weak) dynamical localization due to interference is generically implied by short-time classical recurrences.

Classical localization.– Idealized hard chaos is an exception; a generic finite undriven system features a structured mixed phase space that supports both chaotic and quasiregular motion. Ergodization over the entire energy surface is not guaranteed, or it might be extremely slow, involving exponentially long timescales. Consequently, the practical definition of localization becomes more subtle. If we define localization as a lack of ergodicity, then it should be clear that it may arise also in a purely classical context, due to fragmentation of the energy surface. This fragmentation can happen either due to energetic barriers between different regions (e.g., self-trapping selftrap ), or due to dynamical obstacles of the Kolmogorov-Arnold-Moser (KAM) type KAM .

Quantum localization.– Coming back to the definition of quantum localization, Heller Heller has suggested defining a measure for it by comparing the quantum and the classical exploration volumes. Quantum ergodicity, in this perspective, means that the ratio of volumes is identical to the value predicted by random matrix theory (RMT): somewhat less than unity due to quantum fluctuations. Quantum localization then means that this ratio is significantly smaller than the RMT expectation.

Intrinsic localization measure.– In the realm of high-dimensional finite systems, the standard localization measures, which are based on contrasting quantum versus classical behavior, typically provide inconclusive results. Due to the mixed, complicated structure of the phase space, quantum-to-classical correspondence is a priori poor, and is not in accordance with the predictions of RMT. Hence, simply comparing quantum to classical spreading volumes would lead to the conclusion that quantum localization always shows up. Correspondence-based measures thus trivialize the concept of quantum dynamical localization and furthermore, can not be extended to systems such as spin chains, where the classical limit is meaningless. It is therefore highly desirable to find an intrinsic measure for localization that is not prejudiced by the existence of a classical limit. Below, we adopt such a correspondence-free definition of dynamical localization.

Additional subtleties.– Additional subtleties are involved in the calculation of the standard measure of localization. It is impossible in practice to verify that a classical system is not ergodic. Very long simulations can be carried out, but their results cannot be trusted: on the one hand, the exploration process might involve exponentially long timescales, e.g., due to Arnold diffusion basko ; on the other hand, long simulations might be infected by numerical inaccuracies. These complications do not arise in the quantum context, as the $t=\infty$ saturation profile of the evolving wave packet can be anticipated via diagonalization of the Hamiltonian (see Sec. III below). Consequently, a time-dependent quantum simulation is not even required in order to determine whether the system is quantum-ergodic. The entire analysis of quantum localization can be formulated based on the dependence of known asymptotic quantum distributions on the initial preparation, without any need for dynamical simulations.

Thermalization.– The typical thermalization setup consists of two weakly coupled subsystems. In introductory statistical mechanics textbooks, one of the subsystems is assumed to be a large chaotic environment. It is then argued that the ergodization of the combined system implies equilibration of the smaller subsystem (which does not have to be chaotic). In point of fact, since chaotic ergodization necessitates at least two degrees of freedom, either one or both subsystems can be chaotic, but neither have to be large. It is then possible to define a reaction coordinate $x$ that quantifies the exchange of energy or particles between two such moderately sized subsystems, and to show that the associated dynamics obeys a Fokker-Planck equation trmPRL . The question arises whether upon quantization the implied thermalization process is arrested due to dynamical localization.

Minimal setup.– In a recent paper Chris we looked for a dynamical localization effect in the many-body thermalization of a minimal Bose-Hubbard model. Namely, we have considered an isolated 4-site system with interacting bosons, where three sites are strongly coupled (the trimer) and an additional site (monomer) is weakly coupled to them. The combined model has 3 degrees of freedom (corresponding to occupation differences and relative phases between the sites). The $x$ coordinate was the occupation of the trimer. If the dynamics is treated classically, the $x$ variable thermalizes within a region that is determined by the boundaries of the chaotic sea; i.e., the $x$ distribution obeys, as expected, a Fokker-Plank equation. By contrast, the quantum dynamics exhibits dynamical localization if $x$ starts at the peripheral regions (i.e., if the simulation is initiated with a large occupation imbalance). We have developed a theory for finding the mobility edge that separates the $x$ region of quantum ergodicity from the regions where localization manifests itself. For that purpose we have combined the break-time phenomenology of Chirikov and followers Chirikov ; Shep ; Dittrich ; Cohen with the phase space exploration phenomenology of Heller Heller .

Generic setup.– In this work we explore dynamical localization in the more generic case of trimer-trimer thermalization. The main conceptual difference from the trimer-monomer model is that in the present case both subsystems are chaotic. The motivation for this study is related to the possible implication of such an observation: namely, manifestation of dynamical localization in the trimer-trimer system implies that it is relevant for the analysis of thermalization in larger disordered arrays, since such arrays can be viewed as chains of weakly interacting subsystems, some (or all) of which are chaotic. It should be re-emphasized that, unlike the disorder-induced many-body localization MBL ; MBLb , the cause of dynamical localization in this work is interaction-induced chaos on the microscopic scale, and the thermodynamic limit is not an issue.

Effective $\hbar$.– The Bose-Hubbard model is formally equivalent to that of coupled oscillators, and therefore one may use the traditional distinction between classical and quantum descriptions of the same system. In this context the number of particles per subsystem plays the role of inverse $\hbar$; namely, proper scaling of the dynamical variables gives an effective dimensionless Planck constant ${\hbar_{\text{eff}}=1/N_{\text{subsystem}}}$.

The localization question.– By now the manifestation of dynamical localization in low-dimensional chaos, such as the kicked-rotor, is hardly surprising. But as we go up in dimensionality the observation of a quantum localization effect becomes questionable. The break-time phenomenology suggests quantum localization whenever the classical exploration becomes slow in some sense (see Chris ), but the argument is rather speculative and has never been tested for high-dimensional chaos, except for the trimer-monomer system. It has to be realized that the trimer-trimer system is much more challenging for dynamical localization because the peripheral regions of the chaotic sea are vast, unlike those of the trimer-monomer system. Indeed, as shown below, dynamical localization in a trimer-trimer system is a rather weak effect that requires a refined numerical procedure for its detection. A major challenge is to distinguish a “novel type” of dynamical localization that should not be confused with energetic metastability (self-trapping), with perturbative localization, or with classical localization due to KAM structures. This novel type of dynamical localization appears within the chaotic sea and cannot be explained in a simple way.

Outline.– In Sec. II we introduce the model Hamiltonian, characterize the trimer as a subsystem, and define the required representation for the discussion of the thermalization process. Numerical demonstration of dynamical localization is provided in Sec. III. The ergodicity measure whose objective is to identify the mobility edges of the chaotic sea is defined in Sec. IV; it is used for the analysis of classical localization in Sec. V and of quantum localization in Sec. VI. Finally, in Sec. VII we emphasize the nontriviality of the observed dynamical localization, as opposed to perturbative localization on the one hand, and semiclassical nonergodicity on the other.

II The model

Our building blocks are a pair of Bose-Hubbard trimers, labeled by $\alpha=L,R$. Each trimer is described by the Hamiltonian

[TABLE]

The operators $\hat{a}_{\alpha j}^{\dagger}$ and $\hat{a}_{\alpha j}$, respectively, create or annihilate a particle on site $\alpha j$, while $\hat{n}_{\alpha j}=\hat{a}_{\alpha j}^{\dagger}\hat{a}_{\alpha j}$ counts particles in this site. The parameter $K$ is the tunneling strength between neighboring sites, while $U$ describes interaction between two particles. In the absence of intertrimer coupling (as well as for our particular choice of coupling; see below), the populations $N_{\alpha}=\sum_{j}\hat{n}_{\alpha j}$ are separately conserved. The total number of particles is $N=N_{L}+N_{R}$. In our simulations we assume for simplicity that ${N_{L}=N_{R}=N/2}$.

In order to simulate a thermalization process, the two trimer subsystems are weakly coupled. For this purpose we employ a nearest-neighbor interaction of particles that occupy the same location,

[TABLE]

with $V$ denoting the nearest-neighbor interaction strength. This extended Bose-Hubbard model Dutta15 ; Dey17 can be realized using the long-range dipole-dipole interaction between lattice-BECs of particles with electric or magnetic dipole moments Lahaye09 ; Trefzger11 .

The dimensionless parameters of the model are

[TABLE]

Without loss of generality we set the units of time such that the hopping frequency is ${K=1}$. This holds for all the numerical results that will be presented.

In what follows, we refer to the exact many-body evolution under the Hamiltonian of Eq.​ (2) as the quantum dynamics. Its large-$N$ classical limit is obtained meanfield by replacing the field operators $\hat{a}_{\alpha j}$ by $c$-numbers ${a_{\alpha j}=\sqrt{N_{\alpha}I_{\alpha j}}\mathrm{e}^{i\varphi_{\alpha j}}}$, where the normalized occupations ${I_{\alpha j}\in[0,1]}$ and the phases ${\varphi_{\alpha j}\in[0,2\pi)}$ are canonical action-angle variables. Since the bosonic Hamiltonian is formally equivalent to that of a coupled-oscillators system, one may use the traditional distinction between classical and quantum descriptions of the same system meanfield . The commutation relations imply that the dimensionless Planck constant is ${\hbar_{\text{eff}}=1/N_{\alpha}=2/N}$. The classical form of the Hamiltonian is obtained formally in a straightforward manner from the second-quantized Hamiltonian using the above-defined substitutions, and the prescription ${\mathcal{H}^{cl}=\hbar_{\text{eff}}\mathcal{H}}$. This leads to

[TABLE]

with

[TABLE]

Thus we have only two dimensionless parameters $u$ and $v$ that determine the (scaled) classical dynamics, while in the quantum version we have the additional dimensionless parameter $\hbar_{\text{eff}}$.

II.1 The trimer

The chaotic dynamics of at least one subsystem is a necessary ingredient in thermalization. The characterization of chaos for a single trimer is therefore essential before considering the coupling between two such subsystems. The detailed description of the trimer and the methodology of chaos determination are deferred to Appendices A and B, respectively. Our main objective is to set the range of energies in which the trimer is classically chaotic. The straightforward way to identify chaos is to run classical simulations throughout phase space for each value of the energy. This procedure is inefficient and time-consuming. However, it is well established that classical chaos can be determined from the level spacing statistics of the quantum spectrum HaakeTextbook . Thus, via one diagonalization of the quantum Hamiltonian we obtain the complete information on the chaotic energy range for any value of the characteristic parameters. From the same diagonalization we can extract additional information on the underlying chaos by considering the quantum eigenstates. According to the eigenstates thermalization hypothesis (ETH), a somewhat misleading term in the present context, the eigenstates that dwell within the classically chaotic phase space regions are expected to be ergodically spread throughout them Olshanii . This ergodicity is reflected in the expectation value of observables, which attain their thermal values.

Moderate-strength interactions in the trimer give rise to chaos. However, stronger interaction may also cause self-trapping, which is localization within an isolated part of the energy surface that is surrounded by a (higher or lower energy) forbidden region. The splitting of the energy surface into separate, disconnected regions can be reflected in a splitting in the distribution of the site occupation.

When all the particles are condensed within a single one-particle state (either a site state or some orbital that constitutes a superposition of site states), the resulting many-body state is known as a “coherent state”; it is represented by a Gaussian-like distribution in phase space.The spread of the eigenstate in phase space can be measured by its purity (see definition in Appendix A). The inverse purity $\gamma^{-1}$ estimates the number of one-particle states that participate in the formation of the many-body state. High purity of order unity implies that only a single state is occupied. Low purity implies that the occupation is fragmented; for a trimer eigenstate the lowest purity is $1/3$ (as there are three independent one-particle states).

The spectrum of an isolated trimer (say, ${\alpha=R}$) is illustrated in Fig.​ 1(a) for ${u=6}$. Each point represents an eigenstate $\left|\varepsilon_{R}\right\rangle$ of the Hamiltonian $\mathcal{H}_{R}$, positioned horizontally according to its eigenenergy ${\varepsilon\equiv\left\langle\mathcal{H}_{R}\right\rangle/N_{R}}$, and vertically according to the normalized expectation value of the middle-site occupation, ${I_{2}\equiv\left\langle\hat{n}_{2}\right\rangle/N_{R}}$. The split of the distribution of occupations in Fig.​ 1(a) provides a clear numerical identification of the self-trapping border, indicated by a vertical dotted line.

The points in Fig.​ 1(a) are color-coded by the inverse one-particle purity of the eigenstates $\gamma^{-1}$. Detailed analysis of its implications is provided in Appendix A. In the chaotic energy region one expects low purity, because eigenstates are spread over large chaotic regions of the energy surface. But the precise determination of the chaotic region requires further effort. Fig.​ 1(b) displays the $r$ measure of level-spacing statistics (see Appendix B for definition and a detailed analysis). The value ${r\approx 0.530}$ indicates RMT statistics, the Gaussian orthogonal ensemble (GOE) in particular, that is expected for a chaotic repulsion of levels in time-reversal-invariant systems. This is opposed to the value ${r\approx 0.386}$ that indicates Poissonian statistics, innate to quasi-integrable regions. Consequently, we identify the chaotic range ${0.26<\varepsilon<1.23}$, based on the halfway mark ${\left\langle r\right\rangle>(0.386+0.530)/2}$.

The chaotic range has been further verified by direct classical phase space analysis as in trmPRL , inspecting the Poincaré sections generated by trajectories at different energies. Here we provide another piece of proof that will become useful in later analysis. In the classical limit, the Hamiltonian Eq.​ (7) describes the dynamics of coupled nonlinear oscillators. From it we can derive the classical equations of motion for the occupations, and then propagate them. For motion on an energy surface ${\mathcal{H}^{cl}_{\alpha}(I_{1,2,3},\varphi_{1,2,3})=\varepsilon}$ we can obtain classical power spectra of the site occupations (see Appendix C). This is shown in Fig.​ 1(c) for the central trimer site $I_{2}$; the picture for $I_{1,3}$ is qualitatively similar. Within the central energy region the classical motion has a wide frequency content, which is the hallmark of chaos (see Appendix B).

II.2 Two uncoupled trimers

Considering the two trimers together, yet without interaction, we obtain the energy landscape illustrated in Fig.​ 2. Each point represents an eigenstate $\left|E,x\right\rangle$ of ${\mathcal{H}_{R}+\mathcal{H}_{L}}$, positioned according to

[TABLE]

Note again that scaled units are used for these variables; i.e., the hopping frequency is set to ${K=1}$, and the trimer eigenenergies are normalized per particle. Note also that in our previous work on the thermalization of a trimer-monomer system Chris , particle exchange was allowed and $x$ was the occupation difference. The physics of dynamical localization does not depend on the nature of $x$.

Let us label the single trimer density of states by $g(\varepsilon)$. Considering two uncoupled trimers, their joint density of states $\tilde{g}(E_{0};x)$ with respect to $x$, within an energy shell $(E_{0},E_{0}+dE)$, is given by the product ${g(E_{0}+x/2)\,g(E_{0}-x/2)dE}$, indicated by color in Fig.​ 2. An ergodic distribution is defined as the microcanonical uniform occupation of all states. Therefore, up to a normalization constant, this distribution is

[TABLE]

III Spreading and Ergodization

Below we describe the dynamics in terms of a probability distribution $P_{t}(x)$. For an idealized thermalization process of coupled chaotic subsystems one expects this function to obey a Fokker-Plank equation trmPRL , leading to thermalization such that $P_{\infty}(x)\sim P_{\text{erg}}(x)$. Here we are dealing with a complicated nonhomogeneous mixed chaotic dynamics; therefore we cannot expect a simple diffusive behavior. Still we do expect that the spreading will look quasistochastic, meaning that the spreading has the same rate if we plot it against the scaled time $\tilde{t}=v^{2}t$. This should be contrasted with ballistic transport, for which ${\text{var}(x)\sim vt}$. The quasistochastic time dependence is valid classically whenever the “diffusion” along $x$ is slow compared to the ergodization within the subsystems (this is the weak-coupling assumption). In the quantum case it can be regarded as a Fermi-golden-rule picture (FGR). We shall explain below that for very small $v$ the FGR breaks down, leading to trivial perturbative localization, which is not of much interest for us. Our focus is on couplings $v$ that can be regraded as classically weak, but quantum mechanically large. This is the regime where semiclassical perspective is most appropriate, and this is also the regime where dynamical localization effects are not trivial.

III.1 Quantum dynamics of coupled trimers

Preparing the double-trimer system in one of the coupling-free eigenstates $\left|E_{0},x_{0}\right\rangle$ described above, its exact quantum dynamics is analyzed via the time-dependent probability distribution

[TABLE]

The $x$ range is divided into bins of arbitrary, but small width $\delta x$. For a given initial state $\left|E_{0},x_{0}\right\rangle$ we define the coarse-grained distribution

[TABLE]

Finally, we can average over a sufficiently long time to obtain a relatively smooth saturation profile. Formally we write

[TABLE]

In fact, the exact saturation profile can be obtained directly from the eigenstates of the Hamiltonian. Using the short notation $\left|X\right\rangle\equiv\left|E,x\right\rangle$, and expanding each state in the basis $\left|\mathcal{E}_{n}\right\rangle$ of the full Hamiltonian $\mathcal{H}$, the time-dependent distribution takes the form

[TABLE]

For a nondegenerate spectrum, upon time averaging, the oscillating terms cancel out and one obtains the result ${P_{\infty}(X)=P_{\text{sat}}(X)}$, where

[TABLE]

In the following, the term quantum saturation profile shall always refer to the exact result of Eq.​ (14).

III.2 Semiclassical dynamics of coupled trimers

In the semiclassical picture, a quantum wave packet is represented by a “cloud” of phase space points initially residing on the energy surface $(E_{0},x)$ (see Appendix C). This approach is similar to the truncated Wigner approximation. Coupling the two trimers causes the cloud to spread within the accessible energy shell. For weak intertrimer coupling, if the dynamics is strictly ergodic within the energy shell, the distribution $P_{\text{erg}}(x)$ of Eq.​ (10) is obtained after a sufficiently long time, regardless of the details of the initial preparation.

The semiclassical picture requires a large number of time-dependent simulations. We assume an initial preparation that is represented by a microcanonical cloud of 1000 points, spread uniformly over $I_{\alpha j}$ and $\varphi_{\alpha j}$ values corresponding to ${(E_{0},x_{0})}$ as determined by the unperturbed Hamiltonians $\mathcal{H}^{cl}_{\alpha}$ of uncoupled trimers. The full numerical procedure is described in Appendix C. Individual points of the cloud are propagated under the canonical equations of motion ${db_{\alpha j}/dt=-i\partial\mathcal{H}^{cl}/\partial b_{\alpha j}^{*}}$, where ${b_{\alpha j}\equiv a_{\alpha j}/\sqrt{N_{\alpha}}=\sqrt{I_{\alpha j}}\mathrm{e}^{i\varphi_{\alpha j}}}$. The final distribution of points $P_{t}(x)$ is obtained by counting the number of points falling within each bin $(x,x+\delta x)$ at a time $t$.

III.3 Quantum versus semiclassical simulations

Figure 3 displays the results of representative quantum and semiclassical simulations. The initial energy, $E_{0}$, is the same for all simulations, while the initial trimer energy difference, $x_{0}$, is different for each row of panels (these initial conditions are marked in Fig.​ 2). The value of $v$ is different for each column. The probability distribution $P_{\infty}(x)$ in the long-time limit [based on time evolution for semiclassical simulations, or on the saturation profile of Eq.​ (14) for quantum simulations] can be compared with the ergodic distribution $P_{\text{erg}}(x)$ that is indicated by the black line. We see that semiclassical simulations come close to the ergodic distribution for ${x_{0}=0.01,0.54,0.91}$ but not for ${x_{0}=1.19}$, irrespective of $v$. By contrast, quantum simulations are sensitive to $v$, and either approach the ergodic distribution for a more limited range of $x$ than their semiclassical counterparts, or are completely nonergodic. In the next section our objective is to quantify this observation.

IV Ergodicity measures

IV.1 Microcanonical ergodicity measure

Attempting to quantify ergodicity, the first inclination is to define a measure that would compare $P_{\infty}(x)$ to $P_{\text{erg}}(x)$. For that purpose we choose the overlap function

[TABLE]

A limiting value of $f_{\infty}\rightarrow 0$ indicates an ergodic simulation, whereas $f_{\infty}\rightarrow 1$ is attained when the two distributions have almost no overlap, implying localization. Though very simple, the above ergodicity measure gives a reasonably good estimate that reflects the observed localization in Fig.​ 3. For a smooth distribution, as is the case for semiclassical simulations with large clouds, or quantum simulations with a large $N$, this measure becomes independent of the bin size $\delta x$.

Results for semiclassical simulations are presented in Fig.​ 4. In panel (a) we plot $f_{t}$ against the scaled time ${\tilde{t}=v^{2}t}$ for $E_{0}=0.74$ and for different initial conditions $x_{0}$. The dynamics display an initial transient followed by a rapid decay. A steady state is reached by some of the simulations already for $\tilde{t}\sim 100$. At $\tilde{t}\sim 450$ all simulations with ${x_{0}<1}$ have reached the same value $f_{t}\approx 0.1$. The transient time becomes exponentially large for ${x_{0}>1}$. In panel (b) we show $f_{t}$ at selected time instances for the full range of $x_{0}$ values. This plot demonstrates clearly that there is some kind of “mobility edge” at ${x_{0}\sim 1}$, beyond which the system does not ergodize within a physically reasonable timescale. The different symbols establish the quasistochastic scaling with respect to the strength of the coupling.

Evaluation of the ergodicity measure $f_{\infty}$ for semiclassical simulations with different initial energies is presented in the top row of Fig.​ 5. These results are also incorporated into Fig.​ 6. We see that almost all simulations, other than those lying initially within (or close to) the self-trapped region, reach $f_{\infty}\approx 0.1$, indicating a near-complete ergodization.

IV.2 An intrinsic ergodicity measure

It should be kept in mind that in a mixed phase space system, containing both chaotic and quasi-integrable regions, the energy surface as a whole is not ergodic. It is only the connected “chaotic sea” region that can be regarded as ergodic. Furthermore, in the quantum case some peripheral regions of the sea may not be accessible, if the wave packet has to percolate via sub-Planck corridors. The naive comparison of $P_{\infty}(x)$ and $P_{\rm erg}(x)$ is hence insufficient. Instead, we would like to adopt an intrinsic definition of ergodicity that does not depend on our prejudice regarding the phase space structure, and does not involve $P_{\text{erg}}(x)$. Accordingly, below we use the term ergodic region for the region in which there is no memory of the initial conditions. It is the region where $P_{\infty}(x)$ becomes independent of $x_{0}$. By contrast, in non-ergodic regions the obtained $P_{\infty}(x)$ depends on $x_{0}$, with the particular extreme case of localization, where $P_{\infty}(x)$ is peaked around the initial value $x=x_{0}$.

In principle, a full comparison of the distributions $P_{\infty}(x)$ for different $x_{0}$ is required for our memory-based intrinsic ergodicity measure. In practice, for the purpose of quantitative analysis it is more effective to consider just a single moment of the long-time distribution, and examine its sensitivity to $x_{0}$. Due to the system’s mirror symmetry, and since we set ${N_{R}=N_{L}}$, we choose for this purpose the function

[TABLE]

For plotting purposes this measure is normalized by $x_{\text{max}}$, which is the highest classically allowed value of $x$ for the given energy $E_{0}$. The results for semiclassical simulations are presented in the bottom row of Fig.​ 5, and we see that they are well correlated with the microcanonical $f$ measure. Before considering the quantum results, which are also displayed in Fig.​ 5, we discuss the issue of classical localization.

V Classical localization

The ergodic phase space region in Fig.​ 6 consists of the $(E,x)$ domains where $f_{\infty}$ is low. Outside of this region we distinguish between several classical localization mechanisms (see Appendix A).

The first possibility is quasi-integrability due to quasilinearity. A single trimer is not chaotic in the linear, low energy region of its spectrum. When two such trimers are coupled (the lower corner of Fig.​ 6), if $v$ is small enough the dynamics is likely to remain quasi-integrable, and therefore thermalization would be avoided. In fact, we have to restrict all such statements, and point out that due to Arnold diffusion we can always have thermalization after an exponentially long time.

The second possibility is self-trapping which is localization within an isolated part of the energy surface that is surrounded by a (higher or lower energy) forbidden region. A single trimer becomes self-trapped when its energy is above a well-defined border and is therefore dominated by interaction. When two self-trapped trimers are coupled (the upper corner of Fig.​ 6), or when a self-trapped trimer is coupled to a quasilinear one (the left and right corners of Fig.​ 6), energy exchange is suppressed due to the regular motion in each trimer.

By contrast, if a self-trapped trimer is coupled to a chaotic one, energy exchange should occur. This is formally like coupling an integrable oscillator to a chaotic environment for which thermalization is expected (the situation is similar to that of a harmonic oscillator under stochastic driving, or, in reverse, a system of chaotic oscillators under harmonic driving). However, in our coupled trimer system it looks like we still have classical localization (the upper right strip of Fig.​ 6). We claim that in this region the ergodic time is finite, yet has become exponentially large; an indication for that was already given in Fig.​ 4. Below we give an expanded explanation.

A quantitative way to estimate the rate of thermalization is to use the the Fokker-Plank picture of trmPRL . The diffusion coefficient $D$ for spreading in $x$ (within a narrow energy shell centered at $E$) is given by the classical version of the Fermi-golden-rule (see Appendix C):

[TABLE]

where $\varepsilon_{L,R}=E\pm x/2$. This expression assumes that the two subsystems are independent and weakly coupled, and at least one of them has to be chaotic. The spectral density $S_{jk}(\omega;\varepsilon)$ is the Fourier transform of the cross-correlation function $C_{jk}(\tau)=\left\langle I_{j}(t+\tau)I_{k}(t)\right\rangle$. This linear-response expression implies that the dynamics for different $v$ values should be quasistochastic, meaning that the spreading is the same irrespective of $v$ when plotted against the scaled time $\tilde{t}=v^{2}t$.

The results for $D(x)$ are provided in Fig.​ 7. A probability distribution of $D$ values (generated by using different phase space trajectories with $v=0$; see Appendix C) is presented as a function of $\varepsilon_{R}$, while $\varepsilon_{L}$ belongs to the chaotic range. Whenever the ”R” trimer is self-trapped ($\varepsilon_{R}>1.31$), the range of frequencies associated with the motion becomes narrow, as in Fig.​ 1(c). This leads to an exponentially small result. Also for small $\varepsilon_{R}$ the diffusion coefficient becomes small, but not exponentially small.

Summarizing, thermalization in the classical context is generally expected whenever either one or both subsystems are chaotic. In our previous work on the trimer-monomer system, ergodization is obtained whenever the trimer is chaotic despite the integrability of the monomer. However, here we have encountered a situation where the expected thermalization process is suppressed due to a vanishingly small value of $D$. We note that having quasistochastic spreading is a weaker claim than having diffusion: in regions of mixed phase space the quasistochastic spreading is likely to be anomalous. Therefore we prefer to regard $D$ in the present context as a measure of the coupling between the subsystems and not as a transport coefficient.

VI Quantum localization

We apply the same analysis as in the classical case for the quantum saturation distributions $P_{\text{sat}}(x)$, which are obtained from Eq.​ (14) with different initial values $x_{0}$. The results are presented in Fig.​ 5, and incorporated into Fig.​ 6. In the later figure, each pixel is color-coded according to the value of its ergodicity measure $f_{\infty}$. On the basis of these results we make several observations: (i) Unlike in the classical case, there is no $x$ region where the quantum spreading is strictly ergodic. The $f$ measure is everywhere significantly larger than zero. (ii) In the quantum case there is a strong $v$ dependence. The classical limit is approached as $v$ is increased, but even for a large $v$ ergodicity is not attained. (iii) Even in the central, low-$x_{0}$ region we observe some memory of initial conditions. For example, focusing on $E_{0}=0.74$, only for the largest $v$ do we observe a tiny region (${|x_{0}|<0.25}$) where the saturation is $x_{0}$-independent. (iv) Contrary to the classical case, ergodization is strongly suppressed for large $|x_{0}|$ even in the low-energy regime (small $E_{0}$).

There is a subtle point that should be emphasized. In the classical case, there is a vast $x_{0}$ region where the numerical simulations yield ergodic distributions, as indicated by the low $f$ measure and by the $x_{0}$-independence of $\left\langle|x|\right\rangle_{\infty}$. However, we cannot rule out the possibility of a slow classical ergodization over exponentially long timescales. Thus, it is possible that the observed picture of classical localization is only a transient. The same reservation does not apply in the quantum case, since the saturation profile after an infinitely long time can be obtained via direct diagonalization using Eq.​ (14).

VI.1 Perturbative localization

When the hopping between two unperturbed states of the decoupled system is small compared with their energy difference, elementary perturbation theory stipulates that the true eigenstates of the coupled system remain localized (no mixing). For time-dependent dynamics, this means that the system cannot make a transition between the unperturbed eigenstates. This type of perturbative localization is trivial, and we would like to verify that the localization we observe is more profound. For this purpose, we calculate and display in Fig.​ 8(a) the $x$-basis participation number:

[TABLE]

Whenever the localization is due to perturbative localization we expect to observe $M\sim 2$, with $M=2$ being the minimal value for $x_{0}\neq 0$ due to the mirror symmetry. This perturbative localization takes place whenever the $v$-related couplings are small compared to the mean level spacing. In our system the density of states becomes smaller at peripheral (large $|x|$) regions; hence localization there might be of this trivial type.

VI.2 Localization due to break time

The saturation values of $M$, displayed in Fig.​ 8, are calculated for the same simulations as in the central, $E_{0}=0.741$ column of Fig.​ 5. In regions where $M_{\text{sat}}$ is significantly larger than unity, simple semiclassical reasoning suggests ergodization. Thus, naively, localization is not expected there. A more sophisticated semiclassical reasoning does suggest the possibility of a dynamical localization effect due to a break time in the quantum-to-classical correspondence. We have refined and tested this approach for the trimer-monomer system Chris , where we have demonstrated that the quantum break time (and hence the regions where localization shows up) can be deduced from careful semiclassical analysis. In the present study, we do not aim to repeat the same type of analysis for the double-trimer configuration, but rather to determine whether nontrivial dynamical localization shows up.

A general argument that supports our expectation to observe nontrivial dynamical localization suggests the existence of a mobility edge based on a semiclassical picture. According to the break-time phenomenology, the quantum dynamics fails in following the classical spreading if the classical exploration of phase space becomes sub-ballistic. This is always the case for diffusion in one dimension, and also (marginally) always the case in two dimensions, but not in higher dimensions. Still, in higher dimensions we can expect the formation of a mobility edge if the ballistic exploration becomes slow enough. (Note that the term “exploration” is not a synonym for “spreading”. Standard diffusive spreading of a cloud in high dimensions involves a ballistic exploration of phase space cells.) We argue that a slow rate of exploration is typical for systems with mixed phase space. At the peripheral regions of the chaotic sea the timescales for migration between different regions become very slow. This is indeed verified in the present example, where we witness, as $|x_{0}|$ approaches the phase space boundary, a very small $D$ or even exponentially small $D$. Respectively, at low energies the small $D$ arises due to quasi-integrability, while for large energies the exponentially small $D$ is due to self-trapping. Indeed in the quantum simulations, see Fig.​ 6, we witness a mobility edge in both cases.

Looking again at Fig.​ 8(b) we clearly see that in the quantum case there is a residual memory of $x_{0}$, which is a signature of a nontrivial dynamical localization effect. This lack of quantum ergodicity persists also in regions where $M$ is significant larger than 2. The region where we can identify a small quantum-ergodic domain that features $x_{0}$-independent saturation appears only for relatively large coupling ($v=0.2$), and it is rather small (${|x_{0}|<0.25}$) as opposed to the classically observed ergodicity.

VI.3 The approach to saturation

Remember that the whole discussion of quantum localization does not require any simulations, since the formal saturation profile $P_{\infty}(x)$ can be obtained via Eq.​ (14). Optionally, we could have used Eq.​ (11) with some large $t$. However, as explained below, the two procedures are not practically equivalent.

Besides the quasistochastic timescale of the semiclassical dynamics, we have another (much longer) timescale that has to do with the quasidegeneracies of the spectrum (see Appendix D). Due to the mirror symmetry of the double-trimer system, it is clear that all eigenstates have a definite parity and therefore the $x$-projected saturation profile satisfies $P_{\text{sat}}(-x)~{}=~{}P_{\text{sat}}(x)$. Contrary to this expectation, some long-time-averaged simulations give ${P_{\infty}}(-x)~{}\neq~{}{P_{\infty}}(x)$. This happens due to quasidegeneracies in the perturbed spectrum that are present for small $v$ and can be barely resolved in a finite-time simulation. These quasidegeneracies are remnants of the exact $v=0$ degeneracy of the $\left|E,\pm x\right\rangle$ states.

The time evolution of $M_{\infty}/M_{\text{sat}}$ is provided in Fig.​ 8(b). This ratio can be regarded as a measure for the relative contribution of quasidegeneracies. Whenever the the initial state significantly overlaps quasidegenerate pairs, we get for long time ${P_{t}(x)=2P_{\text{sat}}(x)}$ and therefore $M_{t}/M_{\text{sat}}=1/2$. A much larger timescale is required in order to recover ${P_{\infty}(x)=P_{\text{sat}}(x)}$. This discussion clarifies our choice of $\left\langle|x|\right\rangle_{\infty}$ rather than $\left\langle x\right\rangle_{\infty}$ as a measure for the characterization of the spreading profile. Strictly, the latter will be always get to zero in the infinite time limit irrespective of localization, merely due to the mirror symmetry of the system.

VII Summary

Localization is commonly viewed as an interference phenomenon that leads to a breakdown of quantum-classical correspondence. Such definition is not intrinsic: it requires an a priori definition of some reference space and a measure that compares the actual quantum spreading to a different (classical) dynamics. This point of view is quite frustrating when dealing with small quantized systems, because in practice the correspondence is quite bad to begin with. We therefore prefer to view localization as a lack of ergodicity, and to provide a measure for intrinsic ergodicity. We have defined ergodicity as the case of having some space within which the saturation profile becomes independent of the initial conditions. Accordingly, we could analyze separately the ergodicity in the classical case and in the quantum case, and then independently compare the two.

One should remember that several mechanisms can lead to localization. The most trivial one, as in the case of a single trimer with large $u$, is self-trapping. This is classical localization due to energetic stability within an isolated region of the energy surface. A more complex type of classical localization arises due to dynamical stability of quasilinear motion. The latter case is endangered by Arnold diffusion, which, in practice, can be ignored due to unrealistic time scales.

In realistic models of the type we study, where we have mixed phase space, the microcanonical ETH assumption Olshanii is typically violated, because vast regions of the energy shell are not accessible. Even if the classical dynamics is ergodic on very long time scales, the quantum will fail to penetrate peripheral regions due to the dynamical or perturbative localization effect, and therefore the expectation value of the observable will not approach the microcanonical one as expected from ETH. Our definition of quantum ergodicity does not make any assumption about the ergodic distribution (microcanonical or not).

Our aim in this work was to highlight the emergence of the quantum dynamical localization effect in the thermalization process of weakly coupled subsystems. The term weak coupling is used in a classical sense; quantum mechanically, the coupling might lead to mixing of many levels. In fact, the extreme case of weak coupling in the quantum sense is not very interesting at all, because it leads to a perturbative localization, indicated by ${M\sim 2}$ participation numbers (i.e., no spreading happens). However, we do observe clear signs of quantum localization even when the coupling is not trivially small, so that $M$ is large. We were thus able to demonstrate a quantum mechanical loss of intrinsic ergodicity despite having nonperturbative quantum mixing. This is the nontrivial type of quantum dynamical localization that we were looking for.

Acknowledgment.– The present line of study has been inspired by discussions with Shmuel Fishman who passed away recently (2-April-2019). This research was supported by the Israel Science Foundation (Grant No. 283/18).

Appendix A The eigenstates of a single trimer

The Hamiltonian of a one-particle linear trimer (${N_{\alpha}=1}$ and ${u=0}$) is diagonal in the basis

[TABLE]

with the corresponding energies

[TABLE]

where ${\omega_{0}=K/\sqrt{2}}$. We can define operators which annihilate a particle in one of these orbitals,

[TABLE]

Using the corresponding occupation operators ${n=a^{\dagger}a}$ the many-body Hamiltonian can be rewritten as

[TABLE]

The classical values ${n_{j}/N_{\alpha}}$, or the quantum expectation values ${\left\langle n_{j}\right\rangle/N_{\alpha}}$, are denoted by $I_{j}$.

If the spectrum of the trimer Hamiltonian was presented on a 2D grid ${(\varepsilon,I_{\text{dark}})}$ it would form a triangular lattice of points representing all possible occupations, with corners at ${(-\omega_{0},0)}$, ${(0,1)}$, and ${(\omega_{0},0)}$. These corners correspond to coherent states, meaning that all particles occupy a single orbital. More generally, each eigenstate can be characterized by its one-particle purity. For this purpose we define the one-particle probability matrix

[TABLE]

The purity is then defined as

[TABLE]

The inverse purity tells us what is the number of one-particle states (whether orbitals or site states) that “participate” in the formation of the many-body state, ranging from 1 (e.g., for coherent states) to 3 (maximally mixed states).

Once the interaction is turned on (${u>0}$), the situation changes. The energy landscape is illustrated in Fig.​ 9. At the low energy range one has remnants of (almost) unperturbed eigenstates (purity $\approx 1$). In the upper energy range one should distinguish between two borders whose determination is discussed in Appendix B. The first one is the edge of the chaotic sea: up to this energy each eigenstate resembles a random wave that dwells on the classically chaotic part of the energy surface (purity $\approx 1/3$). The second border is the threshold for self-trapping: above this value the energy surface is composed of 3 disjoint regions around $I_{j}\sim 1$ with ${j=1,2,3}$, with classical clouds being unable to migrate between these regions. Quantum mechanically, we see the formation of states belonging to two distinct types: states that predominantly occupy the middle site (purity $\approx 1$), or cat superpositions of states that are self-trapped in either site 1 or site 3 (purity $\approx 1/2$).

Appendix B Chaotic range determination

Chaotic classical motion can be identified by directly observing the dynamics of individual trajectories on various energy surfaces, for example, via plotting Poincaré sections. Another possibility is to inspect the power spectrum for a variable of interest: a defining property of a chaotic motion is that the Fourier transform of a fluctuating variable exhibits a continuous frequency spectrum. By contrast, the Fourier spectrum of a regular trajectory is made of zero-width delta peaks.

Alternatively, the underlying classical chaos is indicted in the spectral statistics of the many-body system. Let $\varepsilon_{n}$ be an ordered set of energy levels, and ${s_{n}\equiv\varepsilon_{n+1}-\varepsilon_{n}}$ the nearest-neighbor level spacings. Following Ref.  r_measure we consider the ratios

[TABLE]

For classically integrable systems with two or more degrees of freedom, the quantum levels typically tend to cluster. If a parameter in the Hamiltonian is varied, level crossing may be observed. The levels appear to arise from uncorrelated events in a Poissonian random process, resulting in the probability distribution ${P(r)=2/(1+r^{2})}$ r_measure . Its average value ${\left\langle r\right\rangle_{\text{Poisson}}\approx 0.386}$ is therefore considered representative for integrable classical motion.

In the opposite case, for classically nonintegrable systems whose phase space is dominated by chaos, the levels are correlated and a strong repulsion is observed. For such systems the spectral statistics are typically modeled by random matrix theory methods HaakeTextbook . In particular, the Gaussian orthogonal ensemble (GOE) is employed for time-reversal-invariant systems. A numerical estimation of the distribution of $r$ yields ${\left\langle r\right\rangle_{\text{GOE}}\approx 0.530}$ r_measure ; this value is indicative of fully chaotic classical motion.

A practical determination of the chaotic energy range is achieved by locating the energies at which the distribution of $r$ values agrees to an equal degree with both limiting cases. In Fig.​ 1(b) we plot the local average $\left\langle r\right\rangle$, taken over many levels within a small energy window, based on the spectrum of the trimer Hamiltonian, Eq.​ (1). Chaos is predominant for ${0.26<\varepsilon<1.23}$, since in this region of the spectrum we have ${\left\langle r\right\rangle>(\left\langle r\right\rangle_{\text{Poisson}}+\left\langle r\right\rangle_{\text{GOE}})/2}$.

Remarks.– One should be careful when considering energies lying in the parts of the spectrum where the density of states is low. The number of levels should be sufficient for meaningful statistics, otherwise the local average $\left\langle r\right\rangle$ becomes strongly dependent on the size of the averaging window. In Fig.​ 1(b) those bad regions were removed.

Additionally, when evaluating the distribution of $r$ values one must take into account only levels with the same values for all conserved observables except the Hamiltonian. For the trimer this implies using only states of same parity (see Appendix D).

The trimer spectrum.– A detailed examination of Fig.​ 1(b) shows that fully chaotic classical dynamics is expected in the range ${0.45<\varepsilon<1}$, where ${\left\langle r\right\rangle\approx\left\langle r\right\rangle_{\text{GOE}}}$. Indeed, as verified by Poincaré sections, a typical phase space trajectory explores the entire energy surface. For other energies, where $\left\langle r\right\rangle$ is lower, the phase space is mixed; here regular “islands” exist within a chaotic “sea”. In both the low- and high-energy edges of the spectrum we have ${\left\langle r\right\rangle\approx\left\langle r\right\rangle_{\text{Poisson}}}$, indicating quasi-integrability, as explained below.

The trimer landscape.– The trimer energy landscape is presented in Fig.​ 9. Close to the ground state, ${I_{1}\approx I_{2}\approx I_{3}}$ everywhere on the energy surface. The classical Hamiltonian, Eq.​ (7), is predominated by its linear term; without a sufficiently strong nonlinearity chaos cannot develop. At high energies the situation is opposite: here the nonlinearity is too strong, causing the phase space to fragment into isolated regions above ${\varepsilon=1.31}$ (self-trapping). A phase space trajectory remains forever trapped in its initial region, and hence can never explore the entire energy surface. This reduced exploration means that chaos cannot develop, and hence quasi-integrability is restored. [In fact, a weaker form of self-trapping appears already for ${\varepsilon>1}$, when some trajectories become bounded to one side of the line ${I_{1}=I_{3}}$ for very long times. This is the cause of the sudden increase in the dispersion of $D(x)$ in Fig.​ 7.]

Appendix C Classical simulations

Constructing a semiclassical cloud.– The initial condition for any semiclassical simulation is a cloud of 1000 points in the phase space of the double trimer, uniformly distributed over the entire energy surface $(\varepsilon_{L0},\varepsilon_{R0})$, or equivalently, $(E_{0},x_{0})$. Below we outline the procedure used to generate this distribution.

Consider a single trimer and let the set ${P=\left\{I_{1,2,3},\varphi_{1,2,3}\right\}}$ represent a phase space point belonging to an energy surface $\varepsilon_{0}$. First, a large number of points is generated by making random draws from a uniform distribution, taking ${\varphi_{1,2,3}\in[0,2\pi)}$ and ${I_{1,2}\in[0,1]}$. The third occupation is calculated from ${I_{3}=1-I_{1}-I_{2}}$; all points for which ${I_{3}<0}$ are nonphysical and thus discarded. Next, we use Eq.​ (7) to evaluate the energies $\varepsilon_{P}$ associated with each point. Whenever the result is too far from the required value (in practice, when ${|\varepsilon_{P}-\varepsilon_{0}|>10^{-4}}$), that point is discarded. Finally, after collecting 1000 good points for each trimer, we pair them at random (taking one from each trimer) and evaluate the complex amplitudes ${b_{\alpha j}=\sqrt{I_{\alpha j}}\mathrm{e}^{i\varphi_{\alpha j}}}$. Dynamics of the entire semiclassical cloud are generated by independently propagating each point under the equations of motion ${db_{\alpha j}/dt=-i\partial\mathcal{H}^{cl}/\partial b_{\alpha j}^{*}}$, where $\mathcal{H}^{cl}$ is given by Eq.​ (6).

Trimer power spectrum.– Consider a trimer phase space trajectory moving on an energy surface $\varepsilon$. The trajectory generates a fluctuating variable $I_{j}(t)$, for which we define a scaled Fourier transform

[TABLE]

The fluctuations of $I_{j}(t)$ are quantified by its power spectrum,

[TABLE]

When the phase space is fully chaotic, each trajectory eventually explores the entire energy surface, and therefore has the same frequency content. However, a numerical estimate of $Y_{j}(\omega)$ is always limited by the total simulation time $T$ and the resolution $dt$; to reduce inaccuracies in $S_{jj}(\omega;\varepsilon)$, one can average over the spectra generated by several trajectories. Furthermore, in the case of a mixed phase space each trajectory may have a different frequency content, depending on the explored phase space regions. A microcanonical averaging (i.e., over a cloud of trajectories initially uniformly spread over the energy surface) gives a smoothed power spectrum representative of the typical motion on the energy surface. The power spectrum of $I_{2}(t)$, generated by such a cloud containing 1000 trajectories for each value of $\varepsilon$, is presented in Fig.​ 1(c).

More generally we define also cross-correlation functions,

[TABLE]

Their Fourier transforms are

[TABLE]

Diffusion coefficient.– The classical equations of motion for the double-trimer system yield

[TABLE]

The spreading in $x$ can be estimated from

[TABLE]

where $\langle\cdots\rangle$ implies a microcanonical average. Assuming that the motion is chaotic, it makes sense to regard the $I(t)$ functions as stationary noisy functions. For small enough $v$ we can write the right-hand side in terms of cross-correlation functions of the ”L” and ”R” uncoupled trimers:

[TABLE]

Assuming a short correlation time, and disregarding drift (which can be always neglected for short times), we deduce that

[TABLE]

From the definition ${\text{var}(x)=2Dt}$ we arrive at Eq.​ (17). In practice, $D$ is evaluated from finite trajectories. Since the trimer phase space is mixed, some of those trajectories are quasi-integrable and some are chaotic. Thus, from an ensemble of trajectories we get some distribution of $D$ values, as shown in Fig.​ 7.

The theoretical estimate of Eq.​ (17), where formally $v\rightarrow 0$, can be tested by using semiclassical simulations with a small $v$. We evaluate ${\partial\,\text{var}(x)/\partial t=2D}$ for the time range ${0.01<v^{2}t<0.1}$, where to a good approximation the variance grows linearly with time. The results are shown in Fig.​ 10 (symbols) as a function of $\varepsilon_{R}$, with $\varepsilon_{L}$ being kept constant. The theoretical $D$ value is given by the mean of the $D$ distribution for the same values of $\varepsilon_{L,R}$ (line). The agreement is excellent across a wide range of energies, even when the “R” trimer is within the self-trapping region (${\varepsilon_{R}>1.31}$).

Appendix D Symmetry properties of eigenstates

A single trimer, shown in the inset of Fig.​ 1(a), has one nontrivial symmetry operation: the exchange of its edge sites. Hence, the eigenstates of the trimer Hamiltonian $\mathcal{H}_{\alpha}$ of Eq.​ (1) are either symmetric or antisymmetric under this operation. The Hilbert space in which they reside is split into two subspaces, each containing only states of the same parity, with the symmetry imposing $I_{1}=I_{3}$ for all eigenstates.

The double trimer, shown in the inset of Fig.​ 2, contains three nontrivial symmetry operations: $P_{1}$ that exchanges edge sites in both trimers, ${(\alpha 1\leftrightarrow\alpha 3,\ \alpha=L,R)}$; $P_{2}$ that exchanges all corresponding sites between the two trimers, ${(Lj\leftrightarrow Rj,\ j=1,2,3)}$; and their product $P_{1}P_{2}$. We label the states symmetric or antisymmetric under $P_{1}$ ($P_{2}$) by $S$ or $A$ ($s$ or $a$), respectively. The resulting four irreducible representations (correspondingly labeled by $Ss$, $Sa$, $As$, and $Aa$) are associated with four Hilbert subspaces. A new basis composed of superpositions of Fock states, each one symmetric or antisymmetric under $P1$ and $P2$, takes the Hamiltonian of Eq.​ (2) into a block-diagonal form, with approximately equal-sized blocks. This transformation significantly reduces the computational difficulty of exact diagonalization, making it possible to reach higher particle numbers.

When discussing energy exchange between weakly coupled trimers, the most intuitive energy basis is the set of direct products of the single-trimer eigenstates, $\left|E,x\right\rangle=\left|\varepsilon_{n}\right\rangle\left|\varepsilon_{m}\right\rangle$, with the mean energy $E=\varepsilon_{n}/2+\varepsilon_{m}/2$ and the imbalance $x=\varepsilon_{n}-\varepsilon_{m}$. A generic direct product (for $x\neq 0$) is a superposition of two same-energy states $\left|E\right\rangle$ belonging to one of the combined subspaces $S=Ss+Sa$ or $A=As+As$,

[TABLE]

The special case of $x=0$ consists of the states

[TABLE]

The spectral properties of all four subspaces are similar, and hence the dynamics generated by them will also be similar. We consider only eigenstates belonging to the combined $A=As+As$ subspace, thus avoiding the complications induced by the $\left|E,0\right\rangle$ states.

For $v=0$ both the direct products $\left|E_{0},\pm x_{0}\right\rangle$ and their superpositions $\left|E_{0},x_{0}\right\rangle\pm\left|E_{0},-x_{0}\right\rangle$ are energy eigenstates. For a finite, but very small $v$ the spectrum consists of pairs of near-degenerate states that evolve from these two superpositions. In reverse, for a sufficiently small $v$ the states $\left|E_{0},\pm x_{0}\right\rangle$ are mainly combinations of two near-degenerate eigenstates, and hence the flow of probability from $x_{0}$ to $-x_{0}$ is extremely slow. Even for very long simulation times we may observe the probability distribution spreading over a wide range of $x$ yet not crossing over to the other side of the $x$ axis. This creates two distinct time scales: one for spreading (by diffusion), and the other for saturation (by tunneling).