Quantum thermalization: anomalous slow relaxation due to percolation-like dynamics

We highlight a dynamical anomaly in which the rate of relaxation towards thermal equilibrium in a bi-partite quantum system violates the standard linear-response (Kubo) formulation, even when the underlying dynamics is highly chaotic. This anomaly originates from an ℏ-dependent sparsity of the underlying quantum network of transitions. Using a minimal bi-partite Bose–Hubbard model as an example, we find that the relaxation rate acquires an anomalous ℏ dependence that reflects percolation-like dynamics in energy space.

LRT is related to the Fermi-golden-rule (FGR) picture in which the rates of transitions between the unperturbed eigenstates of the subsystems are given by first-order-perturbation matrix elements, but over long timescales that involve many perturbative orders. The diffusion coefficient D of the FPE is estimated from these rates by a Kubo formula [23,24]. LRT implies quantum-to-classical correspondence (QCC) in the FPE description, which is somewhat analogous to the Thomas-Reiche-Kuhn f-sum-rule, and has been termed 'restricted QCC' [25]. The argument that supports restricted QCC with regard to the FPE picture is based on the observation that for short times the variance (unlike the higher moments) features a robust QCC, while for long times the central limit theorem makes all higher moments irrelevant. Thus LRT based description becomes accurate far beyond the naive expectation. The restricted QCC assumption prevails in all current work on thermalization [5][6][7][8][9][10][11][12][13][14][15][16][17][18].
Deviations from LRT have either a classical or a quantum origin. Classical deviations result from dynamical quasi-integrability in the mixed phase space [26,27] which can make thermalization a slow and intricate process [2][3][4]. By contrast quantum anomalies are directly related to the breakdown of QCC due to the finite value of the Planck constant ℏ. One well-known example for such quantum anomaly is the loss of ergodicity due to manybody Anderson localization [20][21][22].
In this Letter we highlight a new type of quantum anomaly which does not originate from the lack of quantum ergodicity, but from the ℏ-dependent sparsity of the quantum network of transitions. The classical-Kubo-FGR picture relies critically on the existence of a dense, connected network of transitions between all the available states, so that all transitions contribute to the diffusive energy spreading process. However, such dense networks do not always exist. The quantum network of transitions is generally sparse [28], resulting in a percolation-like process of energy spreading, that is dominated by bottlenecks and preferred pathways. As a result, the Kubo formula grossly overestimates the thermalization rate and QCC is lost even when the underlying classical dynamics is highly chaotic. Any further distribution of this work must maintain attribution to the author (s) and the title of the work, journal citation and DOI.
To illustrate this point, we consider a minimal Bose-Hubbard model of a bi-partite N-boson system, where =  N 1 plays the role of the Planck constant. We show that while the thermalization process is still described by the FGR picture, resulting in an FPE, it involves an anomalous ℏ-dependent diffusion-coefficient D whose estimate requires a resistor-network calculation. Thus, while the approach to equilibrium still relies on diffusive energy flow with the same long-time stationary energy distributions, the unique mechanism of 'quantum thermalization via percolation' can be much slower than its classical counterpart. Further (technical) details regarding the resistor network calculation; the percolation-like aspect; and its ℏ dependence, are provided in the appendices.

Model system
Consider an isolated system of N bosons in four second quantized modes. The operators â j , â j † and = n a âˆĵ j j † annihilate, create and count particles in site j. The dynamics is generated by the Bose-Hubbard-Hamiltonian (BHH) where U is the on-site interaction, and Ω couples a chain of three sites j = 1, 2, 3. The perturbation  P generates transitions to an additional j = 0 site, namely .c. . Thus  describes a bi-partite system: a BHH trimer coupled to a monomer (see schematic illustration in figure 1). Weak coupling between the two subsystems is assumed (ω Ω ≪ NU , ), and the interaction within the trimer is quantified by the dimensionless interaction parameter Ω = u NU . In the classical description each site is described by conjugate action angle variables φ n ( , ) j j . The standard procedure [29] is to work with dimensionless variables. In particular the scaled occupations are n j /N, hence upon quantization the scaled Planck constant is =  N 1 . The classical limit is attained by taking the limit → ∞ N keeping NU constant. In this limit quantum fluctuations diminish and the bosonic operators can be replaced by c-numbers. The semiclassical description becomes valid if ≪  1. The above trimer plus monomer model is the minimal Bose-Hubbard configuration which allows chaos and thermalization, because the trimer subsystem is classically chaotic [30], while a dimer is not. Furthermore, this type of minimal configuration serves as the building-block for progressive thermalization of large arrays [31,32].

Quantum network of transitions
The trimer population ≡ + + x n n nˆˆ1 2 3 commutes with the unperturbed (ω = 0) Hamiltonian  0 , and therefore constitutes a good quantum number in the absence of coupling. The unperturbed spectrum as defined by the eigenstate equation is plotted in figure 1. Each unperturbed eigenstate is associated with a 'position' x m on the trimer occupation grid. Thus, figure 1 should be interpreted as specifying the unperturbed trimer spectrum for all possible trimer occupations from x = 1 to x = N. We identify the region of chaotic dynamics by a Brody parameter map [33] (see appendix B), verified by classical-Poincare sections (not shown). Eigenstates supported by chaotic phase-space regions are marked in black in figure 1.
The perturbation due to coupling with the additional mode allows transfer of particles and energy and thus generates transitions along the occupation axis x. The transition strengths are given as The upper inset of figure 1 depicts the coupling network within a narrow [x,E] window. Due to the wide distribution of transition strengths, the obtained network is glassy. This glassiness is reminiscent of the sparsity that arises in integrable systems due to selection rules [28].

Diffusive spreading
We focus our attention on the evolution of the probability distribution P t (x), starting with an initial state| 〉 m . This preparation is an eigenstate of the unperturbed Hamiltonian, but a far from equilibrium initial state for the combined system. The system's parameters are chosen such that the energy of this state (diamond blue marker in figure 1) lies within a broad chaotic phase-space window.
A representative example for the evolution of the x probability distribution in the chaotic regime is plotted in figure 2 with the growth of variance x Var( ) depicted in the lower panel. Similarly to the results of [15,16], the hallmark of chaos is stochastic-like spreading. This diffusive behavior persists until the distribution saturates the accessible energy window, thus leading to thermalization.
However, the rate in which the equilibrium distribution is approached is very far from the conventional Kubo estimate and is therefore highly non-classical. The thin solid gray line in the lower panel of corresponds to the classical result. It is evident that the standard classical prediction greatly overestimates the equilibration rate and that indeed quantum thermalization is slower due to the sparsity of the transition network. By contrast, the dotted-dashed blue line also depicts an FPE description, but with a percolation-theory resistor network estimate D x ( ) qm for the diffusion coefficient, that, as described below, takes into account the ℏ dependent transition network sparsity. We thus observe a novel anomalous process of quantum thermalization, which is stochastic and adheres to an FPE description, albeit with an underlying percolation-like spreading process which does not correspond to the classical dynamics.

Evolution of the distribution profile
Several snapshots of P t (x) during the thermalization process are plotted in figure 3, showing good agreement between the percolation-FPE and the full numerical simulation of the four-mode dynamics. By contrast, the conventional classical FPE thermalization gives far broader distributions at the same times.
An additional observation concerns the long time equilibrium distributions, plotted in figure 3(c). The saturation profile ∞ P x ( )of the FPE is proportional, as expected, to the density of states g x ( ). By contrast the exact equilibrium distribution is somewhat non-ergodic. The lack of ergodicity in the low x region of the saturation profile, is due to residual integrability within islands of the underlying mixed phase-space. It therefore disappears when the simulation is started deeper within the chaotic sea, see figure 4. In addition, there are deviations from ergodicity in the high x region due to Anderson-type localization. The former semiclassical effect and the latter quantum anomaly are both distinct from the dynamical anomaly which constitutes our main theme. For further detail on these deviations see section 8.

Stochastic FGR rate equations
The transition rates between two chaotic sub-systems are non-zero provided , where the bandwidth τ 1 is determined by the width of the power-spectrum of the perturbation [16]. The FGR estimate for the non-zero rates is accordingly With these rates, the master equation for the occupation probabilities is n m mn n m Our model is sub-minimal in the sense that the monomer is not a chaotic sub-system. Still, the dynamics is the same as for two chaotic sub-systems with τ 1 determined by the width of the energy shell. Namely Only states within this energy shell, marked by blue lines in figure 1, contribute to the thermalization process. States outside it do not participate in the dynamics. The red dashed lines in figures 2 and 3 correspond to the propagation of equation (4) (see appendix). The agreement with the full quantum simulation validates the stochastic FGR picture.

The FPE description
Coarse graining of the kinetic equations (4) results in the FPE, which is merely a diffusion equation in x space Here g x ( )is the density of states within the allowed energy shell. Unlike the textbook version of the diffusion equation, which assumes uniform g x ( )and D(x), the form of the FPE (6) reflects the simple observation that an ergodic distribution occupies uniformly all accessible eigenstates, so that the FPE ergodic saturation profile must satisfy ∝ ∞ P x g x ( )˜( ). The standard linear response estimate for the diffusion coefficient, i.e. the Kubo formula [23,24], is based on a second moment calculation:   figure 1, compared to the microcanonical (∝g x ( )) thermal distribution (square markers). The quasi-integrable region is marked in gray and an arrow in the chaotic region marks the initial state used in figures 2-5. Non-ergodicity is due to quasi-integrability at the low x region (red lines) and due to Anderson-type localization at the high x region (magenta lines). Quantum thermalization is obtained for intermediate x preparations, regardless of the precise initial conditions (blue lines).

Resistor-network calculation
As mentioned above, the FPE simulation with the standard diffusion coefficient D x ( ) cl fails to reproduce the true dynamics as illustrated in figure 2. This striking breakdown of QCC is due to the percolation-like nature of energy spreading. As appropriate for a percolation process, D(x) should be estimated from the conductivity of the 'resistor network' that is formed by the quantum transitions [28]. Such evaluation gives the proper weight to low-resistance, well-connected links, as opposed to the over-estimated democratic weighing of equation (7). Thus, in steady state equation (4) is formally the same as Kirchhoff's equation

Saturation profile
For completeness we further discuss the saturation profiles of figure 4. Given an initial state (m), we take its overlap with the exact eigenstates (ν), Evolving the initial state m in time we define the probability distribution The P t (x) distribution is related to this kernel by binning together the probabilities of all the unperturbed eigenstates with the same trimer occupation, namely where the summation is over all unperturbed states n with x n = x. Note that while ν | P m ( )is the fixed probability distribution between the exact eigenstates of the composite four-mode system, This relation is obtained by expanding the states| 〉 n and| 〉 m of equation (10) in the ν | 〉basis, assuming that the spectrum is non-degenerate; hence only diagonal terms survive after the long time averaging [35,36]. Note that whenever the Wigner surmise applies, degeneracies have measure zero due to level repulsion. We have verified that equation (11) is in very good agreement with the exact simulation, as demonstrated in figure 3(c).
It thus becomes clear that the deviation from ergodicity is related to the localization of some unperturbedeigenstate preparations| 〉 m , as reflected in the overlaps ν | P m ( ). Several preparations with the same energy but lying in different phase-space regions are marked in figure 1, while their associated saturation profiles are shown in figure 4. Preparations in the chaotic region give the micro-canonical ergodic saturation profile ∝ ∞ P x g x ( )˜( ), independently of the choice of initial state (blue lines). In the low x region of the saturation profile the localization is of semi-classical nature, due to the underlying mixed phase-space which contains remnant quasiintegrable regions. Preparations supported by such integrable islands have narrow ν | P m ( )which leads to localized saturation profiles. At the high x region, the coupling between eigenstates in different x manifolds, as quantified by the value of the diffusion coefficient D qm , becomes small (see figure 5). Consequently, the Anderson localization length ξ π =g D 2 qm is only a few sites, again resulting in localized saturation profiles (magenta lines). The deviation of the saturation profile in this region from the ergodic result of the stochastic FGR calculation (see e.g. Figure 3(c)) indicates that this is an Anderson-type interference effect.

Experimental realization
Few-mode Bose-Hubbard systems can be realized in confining potentials with toroidal shapes and tunable weak links [37][38][39][40]. Of particular relevance for the realization of bi-partite Bose-Hubbard models is the experimental generation of arbitrary and dynamical potentials in a 87 Rb Bose-Einstein condensate by means of a rapidly moving laser beam [38]. Alternatively, the interference of the rotationally-symmetric Gauss-Laguerre laser modes and optical lattices may be used to generate toroidal Bose-Hubbard systems [37] where adjustable weak links may be introduced [40] to separate the ring into two weakly-coupled subsystems. In this context, one simple configuration may be attained by tilting the lattice potential with respect to a four-node Gauss-Laguerre mode, thus generating two adjacent high barriers and two adjacent low barriers along the four-site ring, separating it into a trimer and a monomer. Equilibration can be readily detected by monitoring the populations of the two subsystems as a function of time and full relative-number distributions may be attained by multirealization measurements. As long as the constituent subsystems are weakly-connected, our observations should be independent of the details of the coupling (e.g. which sites of the two subsystems are linked) due to the generic nature of chaotic motion. The interplay between realistic dephasing and particle loss, and the chaotic dynamics will be the subject of future studies.

Discussion
All stochastic descriptions eventually fail to describe quantum coherent processes, because they inevitably lead to a a microcanonical distribution at → ∞ t , whereas the quantum evolution has an infinite memory of the initial conditions. However, the equivalence between the diagonal and the microcanonical ensembles [34][35][36] in the eigenstate thermalization hypothesis picture [5,6,9] implies that in the quantum evolution of classically chaotic systems, the memory of initial conditions is effectively lost over an ergodization period with all initial conditions leading to a microcanonical distribution. On longer timescales, quantum recurrences take place and the memory of initial conditions is regained. It is thus understood that stochastic methods should be evaluated by their ability to describe quantum dynamics within the time scale of interest, i.e. until an ergodic-like distribution for the pertinent observable is attained.
Within this ergodization time, deviations from LRT include both quantum anomalies and semiclassical integrability effects. The former are directly related to quantization and are important for a dynamical view of quantum thermodynamics [41], whereas the latter are related to incomplete chaoticity and residual quasiintegrability regions in the classical mixed phase-space.
Our main objective was to highlight a novel quantum anomaly in the thermalization process of a quantized chaotic system: a bi-partite Bose-Hubbard complex that can be regarded as the building block for thermalization of larger arrays. We have demonstrated that thermalization with finite ℏ is quite different from that of the corresponding ' =  0' classical system. Whereas classical thermalization is captured well by LRT, leading to an FPE with a Kubo estimate for the energy diffusion coefficient, this approximation fails badly upon quantization. The reason for this dynamical anomaly is the sparsity of the network of couplings between the energy eigenstates of the constituent subsystems which leads to percolation-like dynamics of the energy distribution. As a result, while an FPE description still holds (within the timescale of interest), quantum thermalization, properly described by a resistor-network calculation, can be strikingly slower than the corresponding classical process.
The full dimension of the Hilbert space in a tetramer with population N is  = The Hamiltonian of the system can be separated into blocks of smaller dimensions by considering the permutation symmetry between the external trimer sites (sites 2 and 3, in the schematic illustration inset of figure 1). Denoting the population basis by| 〉| 〉 n n n n , , The former is for ≠ n n 2 3 . We restrict the simulations to the antisymmetric subspace which includes less states and therefore allows us to use a higher number of particles. The antisymmetric subspace excluded the possibility of having zero trimer population x = 0.

Appendix B. Identification of chaos by level statistics
Given the parameters N, Ω, ω, U, we find the eigen-energies of the Hamiltonian equation (1) (e.g., figure 1). Dividing the spectrum to small energy intervals, we calculate the mean level spacing and the distribution P(S) of level-spacings in each of them. We then fit it to the Brody distribution [33]  . Here Γ denotes the Euler gamma function. A Brody parameter value of q = 0 indicates a Poissonian level-spacing distribution characteristic of the uncorrelated levels of integrable system. By contrast for q = 1 we have the Wigner level-spacing distribution, that reflects the level repulsion in the case of a quantized chaotic system. Thus, by plotting q as a function of energy we map the domain of chaotic motion, marked in black in figure 1. The result was then ascertained by inspecting classical Poincare sections in the various regions of the map.
In order to illustrate the connection between the deviation from ergodicity of the saturation profiles and the quasi-integrability islands in the mixed phase-space, we employ the initial states marked in figure 1. Some lie well within the chaotic sea, while others reside in an integrable island. The saturation profiles for these states are shown in figure 4, showing a clear connection between integrability and localization. The polynomial ν x EXP ( ) has degree ν, and equals the truncated Taylor expansion of x exp ( ). Its degree is determined by the effective dimensionality d of the network. The linear estimate D linear is what we call here D cl , and s is the sparsity parameter ( ≪ s 1 means sparse network). If the network originates from the quantization of a weakly chaotic system we expect s to be proportional to some power of  1 [43]. Accordingly the ratio ≡ g D D s qm cl reflects the sparsity of the network. In the 'sparse' limit ( ≪ s 1) the expression above resembles that of variable-range-hopping. For small ℏ the network becomes more connected (less sparse) and g s goes to unity. This crossover can be regarded as a smoothed 'percolation' transition.
Form the above discussion it should be clear that sparsity and hence the quantum anomaly diminish in the large N limit. However, it is important to realize that for thermalization of large arrays, which proceeds via progressive process that involves 'chaotic spots' [31], the relevant ℏ is determined by the number of particles per 'spot', and not by the total number of particles in the system.