Quantum thermalization via percolation

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 $\hbar$-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 $\hbar$ dependence that reflects percolation-like dynamics in energy space.

Over a century since the formulation of the Loschmidt paradox, the emergence of thermodynamic irreversibility from microscopic reversible quantum mechanics is still not fully resolved. A connection between thermalization and chaotic ergodicity is well-established for classical systems [1]. However, strict dynamical chaos is absent in isolated many-body quantum systems, which are linear and quasi-periodic. In contrast to nearly-integrable quantized systems in which thermalization is slow and intricate [2][3][4], it is generally accepted that if the underlying classical dynamics is chaotic, there will be good quantum-classical correspondence (QCC) in the energy distribution between the constituent subsystems of an isolated bi-partite quantum system [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Nevertheless, for finite values of the (dimensionless) Planck constant , anomalies may occur in which QCC is violated. One such anomaly is due to a many-body localization effect, as discussed by Anderson and followers [19][20][21]. In this Letter we consider a generic minimal model of a bi-partite Nboson system, where = 1/N plays the role of the Planck constant. We highlight a different type of anomaly which does not originate from the lack of quantum ergodicity, but from the -dependent sparsity of the quantum network of transitions. The implication is that the relaxation of a low-dimensional bosonic system towards equilibrium cannot be calculated using the traditional formulation, even when the underlying dynamics is highly chaotic. Instead, it acquires an anomalous dependence that reflects percolation-like dynamics.
Thermalization in bi-partite mesoscopic systems, assuming generic conditions, is attained via diffusive energy spreading in each of the constituent subsystems, in response to its coupling to the other. This diffusive process, described by a Fokker-Planck-Equation (FPE) [14][15][16][17], eventually leads to ergodization of the composite system over all accessible states within the initial microcanonical energy shell. The traditional Kubo estimate [22,23] for the diffusion coefficient D of the pertinent FPE is based on a Fermi-golden-rule (FGR) picture in which the rates of transitions between the energy eigenstates of either subsystem are given by first-order-perturbation-like matrix elements, but over long timescales that involve many perturbative orders. The Kubo estimate automatically implies 'restricted QCC' [24] because for short times the variance (unlike higher moments) features a robust QCC, while for long times the central limit theorem makes all the higher moments irrelevant. Consequently the FGR approximation in a chaotic system becomes accurate far beyond the naive perturbative expectation.
The 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. Below we demonstrate, using a concrete generic example, that such dense networks do not always exist. The quantum network of transitions is generally sparse, resulting in a percolation-like process of energy spreading, that is dominated by bottlenecks and preferred pathways. This leads to the failure of the Kubo estimate and the breakdown of QCC in a novel way, not related to the loss of quantum ergodicity. We show that while the thermalization process is still described by the FGR picture, resulting in a FPE, it involves an anomalous -dependent diffusioncoefficient 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 can be much slower than its classical counterpart.
While the existence of dynamical quantum anomalies have been previously postulated for driven systems [25], it was neither demonstrated for a physically appealing system, nor considered in the context of thermalization of isolated, non-driven systems. Below we demonstrate the mechanism of "quantum thermalization via percolation" using a four-mode Bose-Hubbard model.
Model system.-Consider an isolated system of N bosons in four second quantized modes. The operatorsâ j ,â † j andn j =â † jâ j annihilate, create and count particles in site j. The dynamics is generated by the Bose-Hubbard Hamiltonian  where U is the on-site interaction, and Ω couples a chain of three sites j = 1, 2, 3. The perturbation H P generates transitions to an additional j = 0 site, namely, The Hamiltonian H thus describes a bi-partite system constituting a Bose-Hubbard trimer coupled to a monomer (see schematic illustration in Fig. 1). Weak coupling between the two subsystems is assumed, i.e. ω Ω, N U . This partition is chosen because a trimer, having two classical degrees of freedom (e.g. two population differences and two relative phases, given that conservation of N renders the absolute phase insignificant), is the minimal Bose-Hubbard model which allows chaos. The interaction within the trimer is quantified by the dimensionless interaction parameter u = N U/Ω.
Quantum network of transitions.-The trimer populationx ≡n 1 +n 2 +n 3 commutes with the unperturbed (ω = 0) Hamiltonian H 0 , and therefore constitutes a good quantum number in the absence of coupling. The unperturbed spectrum as defined by the eigenstate equation H 0 |m = E m |m is plotted in Fig.1. Each unperturbed eigenstate is associated with a 'position' x m on the trimer occupation grid. Thus, Fig.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 [26, a], verified by classical Poincare sections (not shown). Eigenstates supported by chaotic phasespace regions are marked in black in Fig.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 n|H P |m . The upper inset of Fig.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 [25].
Energy spreading.-Due to the conservation of the total particle number N and the total energy H, and because the coupling is weak, the trimer's occupation x determines its energy. We thus focus our attention on the time evolution of the probability distribution P t (x), starting with a 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 Fig. 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 Fig.2 with the growth of variance Var(x) depicted in the lower panel. Similarly to the results of Refs. [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 Fig.2 corresponds to the traditional FPE description of the dynamics, with a diffusion coefficient D cl (x) that 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 dot-dashed blue line also depicts an FPE description, but with a percolation-theory resistor network estimate D qm (x) 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 percolationlike spreading process which does not correspond to the classical dynamics.
Lineshape.-Several snapshots of P t (x) during the thermalization process are plotted in Fig.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 Fig.3c. The saturation profile P ∞ (x) of the FPE is proportional, as expected, to the density of statesg(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. 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 [a].
Stochastic FGR rate equations.-The transition rates between two chaotic sub-systems are non-zero provided |E n − E m | < 1/τ , where the bandwidth 1/τ is determined by the width of the power-spectrum of the per- turbation [16]. The FGR estimate for the non-zero rates is accordingly, With these rates, the master equation for the occupation probabilities is 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 Fig.1, contribute to the thermalization process. States outside it do not participate in the dynamics. The red dashed lines in Fig.2 and Fig.3 correspond to the propagation of Eq.(4) [a]. 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 Hereg(x) is the density of states within the allowed energy shell. Unlike the textbook version of the diffusion equation, which assumes uniformg(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 [22,23], is based on a second moment calculation: where the brackets correspond to averaging over all the in-band states m in the vicinity of x. The result of the D cl (x) calculation is illustrated in Fig.4. We have verified that the obtained values of D cl (x) are robust, i.e. are not sensitive to the exact value of the micro-canonical width 1/τ . Resistor-network calculation.-As mentioned above, the FPE simulation with the standard diffusion coefficient D cl (x) fails to reproduce the true dynamics as illustrated in Fig.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 [25]. Such evaluation gives the proper weight to low-resistance, wellconnected links, as opposed to the over-estimated democratic weighing of Eq.(7). Thus, in steady state Eq.(4) is formally the same as Kirchhoff's equation where the conductances G mn , and the voltages V n , are analogous to Γ nm and p n respectively. In order to calculate the conductance of a small x segment [x 1 , x 2 ], we set I n = 0 for all internal nodes, and I n = ±I source at the endpoints. The detailed numerical procedure is provided in [a]. Solving for the voltage we deduce that the conductance of the x segment is G(x) = I source /(V 2 − V 1 ), and hence the conductivity is D qm (x) = (x 2 −x 1 )G(x). As shown in Fig.4, the resistor-network calculated diffusion coefficient D qm (x) is substantially smaller than the Kubo result D cl (x). As previously stated, the FPE simulation [a] with D qm (x), presented in Fig.2, agrees well with the quantum simulation. The agreement persists as long as the spreading is within the chaotic region of the energy shell, confirming our expectations.
Discussion.-Does an isolated many-body quantum system prepared far from equilibrium thermalize? The short answer remains that it does, as long as the underlying classical dynamics is chaotic. As shown above, thermalization and its relation to chaos are clearly demonstrated even for an embarrassingly simple four-mode Bose-Hubbard model.
However, as shown for the same model, thermalization of a quantum system with finite is quite different from the thermalization of the corresponding ' = 0' classical system. Whereas classical thermalization is captured well by linear response theory, leading to a 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, the rate of quantum thermalization, properly estimated by a resistor-network calculation, can be strikingly different from that of the classical process.
It is important to distinguish between quantum anomalies such as the Anderson-localization effect and the quantum percolation-effect, and integrability effects such as prethermalization [2][3][4] and semiclassical-localization. The former are directly related to quantization and are important for a dynamical view of Quantum Thermodynamics [27], whereas the latter are related to incomplete chaoticity and residual quasi-integrability regions in the classical mixed phase-space.

SYMMETRY SUBSPACES IN THE TETRAMER
The full dimension of the Hilbert space in a tetramer with population N is N = (N + 1)(N + 2)(N + 3)/6. 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 Fig.1). Denoting the population basis by |n 0 |n 1 , n 2 , n 3 , the totally symmetric and the totally anti-symmetric sub-spaces are spanned by the following symmetrized and antisymmetrized superpositions: 1 √ 2 |n 0 |n 1 , n 2 , n 3 ± |n 0 |n 1 , n 3 , n 2 (9) |n 0 |n 1 , n, n .
The former is for n 2 = n 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.

IDENTIFICATION OF CHAOS BY LEVEL STATISTICS AND THE BRODY PARAMETER.-
Given the parameters N ,Ω,ω,U , we find the eigenenergies of the Hamiltonian Eq.(1) (e.g., Fig.1). Dividing the spectrum to small energy intervals, we calculate the mean level spacing and the distribution P (S) of levelspacings in each of them. We then fit it to the Brody distribution [26] P q (S) = αS q exp(−βS 1+q ) (11) with α = (1 + q)β, and β = Γ 1+q [(2 + q)/1 + q)]. 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 Fig.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 Fig.1. Some lie well within the chaotic sea, while others reside in an integrable island. The saturation profiles for these states are shown in Fig.5, showing a clear connection between integrability and localization.

THE RESISTOR-NETWORK CALCULATION
In order to find the diffusion coefficient for a sparse resistor network we rewrite Kirchhoff's law Eq.(8) in a matrix form, where G is the discrete Laplacian matrix of the network, whose diagonal elements are defined as follows: In order to find the conductance of a segment [x 1 , x 2 ] we shortcut the bonds to the left of the segments, hence defining a left lead. Likewise we define a right lead. Then we place a source I 1 = 1 and a sink I 2 = −1 at two nodes on the left and right leads, and solve Kirchhoff's equation using a psaudo-inverse routine.

FGR AND FPE SIMULATIONS
The master equation Eq.(4) can be written in a matrix form as (d/dt) p = W p and has the solution In order to perform a simulation with the FPE Eq. (6) we have to discritize the continuous x variable. There are two possible strategies. One possibility is to define formally a variable n, such that dn/dx =g(x). In this variable the FPE becomes an unbiased diffusion equation: where The discrete version of Eq.(15) is a master equation with near-neighbor hopping. The rates D n are the same in both directions, and the solution is straightforward. The second strategy to solve the FPE, which looks more natural in the present context, is to stay with the x variable. One should realize that in this variable the ergodic state is not uniform. At steady state the current across each x bond is zero, satisfying where w x,x are transition rates between nodes. Selection rules forbid transitions between non-neighboring nodes, thus the W matrix contains only two diagonals at x = x ± 1. But unlike that master equation Eq.(4), here W is a non-symmetric matrix. At steady state the probability distribution is identical to the normalized density of states, hence we deduce the relation Accordingly the forward and backward transition rates that we are using in the FPE simulation are Using the above rates we can solve the FPE using Eq. (14).

SATURATION PROFILE
The local density of states (LDOS) for a given initial state (m) is defined by its overlap with the exact eigenstates (ν), plotted as a function of the exact energy E ν . Evolving the initial state m in time we define the eigenstate distribution, The previously discussed x-distribution is related to this kernel by binning together the probabilities of all the unperturbed eigenstates with the same trimer occupation, namely P t (x) = (x) n P t (n|m) 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, P t (n|m) is the time-dependent probability distribution between the eigenstates of an uncoupled trimermonomer subsystem.
The long time saturation profile of the evolving P t (n|m) distribution, can be obtained directly from the P (ν|m) of the LDOS, via the convolution formula This relation is obtained by expanding the states |n and |m of Eq. (22) in the |ν basis, noting that only diagonal terms survive after long time averaging. This gives very good agreement with the long time limit of the exact simulation, as seen in Fig.3c. It thus becomes clear that the deviation from ergodicity is related to the localized LDOS P (ν|m) of certain  Fig. 1, compared to the micro-canonical (∝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 Fig.2-Fig.4. 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).
unperturbed-eigenstate preparations |m . Several preparations with the same energy but lying in different phasespace regions are marked in Fig.1, while their associated saturation profiles are shown in Fig.5. 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 quasi-integrable regions. Preparations supported by such integrable islands have narrow LDOS 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 Fig.4). Consequently, the Anderson localization length ξ = 2πgD 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. Fig.3c) indicates that this is an Anderson-type interference effect.