Fragile Many Body Ergodicity

Weakly nonintegrable many-body systems can restore ergodicity in distinctive ways depending on the range of the interaction network in action space. Action resonances seed chaotic dynamics into the networks. Long range networks provide well connected resonances with ergodization controlled by the individual resonance chaos time scales. Short range networks instead yield a dramatic slowing down of ergodization in action space, and lead to rare resonance diffusion. We use Josephson junction chains as a paradigmatic study case. We exploit finite time average distributions to characterize the thermalizing dynamics of actions. We identify a novel action resonance diffusion regime responsible for the slowing down. We extract the diffusion coefficient of that slow process and measure its dependence on the proximity to the integrable limit. Independent measures of correlation functions confirm our findings. The observed fragile diffusion is relying on weakly chaotic dynamics in spatially isolated action resonances. It can be suppressed, and ergodization delayed, by adding weak action noise, as a proof of concept.

The conventional perception of evolving dynamical systems with a macroscopic number of degrees of freedom (DoF) is them being in a state of thermal equilibrium, i.e. ergodic. This assumes all allowed microstates having the same probability. It goes along with trajectories visiting the vicinity of all points of the available phase space (i.e. the phase space subject to constraints due to integrals of motion such as e.g. the energy), and infinite time averages equaling available phase space averages [1]. Statistical physics approaches were paved by Gibbs and Boltzmann and provide a straight connection between microcanonical dynamics and the emergence of canonical distributions [1]. The more interesting it is to study cases when this connection is not evident, eroding, or even missing. This can happen (i) for dynamics in the proximity of an integrable limit, (ii) for dynamics in the proximity of nonergodic sets of measure zero (such as periodic orbits), and (iii) for dynamics driven out of ergodicity due to additional constraints (e.g. condensation). Fermi-Pasta-Ulam-Tsingou problems [2][3][4][5][6] can be associated with (ii), and non-Gibbs states for interacting Bose lattice gases [7][8][9][10][11] with (iii). As for (i), the celebrated Kolmogorov-Arnold-Moser (KAM) theorem is available [12], but applies to systems with finite numbers of DoF and dictates weakly non-integrable dynamics to be non-ergodic on a finite measure set of invariant tori, while being ergodic on the complementary one (Arnold diffusion [13]). The KAM borders are assumed to quickly diminish with increasing DoF numbers [12]. What lies beyond those borders for macroscopic systems? The expectation that Gibbs and Boltzmann take over, was shattered by recent results on Many-Body Localization (MBL) [14,15] which show that certain quan-tum many-body systems can resist thermalization at finite distance from integrable limits. With most analytical results being non-rigorous, and computations notoriously heavy due to exploding Hilbert space dimensions, the weakly touched field of ergodization and thermalization of corresponding classical many body systems is in the focus of this work.
Networks of weakly coupled superconducting grains are one of the few paradigmatic examples of systems where the above scenaria have been considered [16][17][18]. Also, related networks of interacting anharmonic oscillators were used to argue for and show the existence of two different classes of nonintegrable perturbations of an integrable Hamiltonian H 0 ({J k }), with a countable set of actions J k (k being an integer) [19]. Nonintegrable perturbations H 1 ({J k , Θ k }) typically span long range or short range networks (LRN or SRN) in the action-angle space (here Θ k are the canonically conjugated angles). A reference action J k in that network is coupled to R k ×L ktuples of other action-angle pairs. L k is typically a single digit integer. R k however can be intensive (SRN) or extensive (LRN) [19].
Consider a typical LRNs, and translationally invariant two-body interactions L k = 3, R k ∼ N 2 with N being the volume (system size) [19]. Chaotic dynamics can develop in a given L-tuple on a time scale T Λ = 1/Λ where Λ is the typical (largest) Lyapunov exponent in the system. Chaotic dynamics develops due to nonlinear resonances which take place when ratios of network matrix elements to certain frequency differences are large [20]. In the proximity to an integrable limit the network matrix elements scale with 1/N [19]. Therefore the probability for an L-tuple to be resonant will be π r /N 1 where π r 1 is an intensive measure of the distance to the integrable limit. However, there are R k ∼ N 2 L-tuples in which one given action is involved. Therefore the probability Π k that the reference action is in resonance with at least one of its L-tuples (i.e. satisfying the resonance condition) turns exponentially close to unity: Π k ≈ 1 − (1 − π r /N ) R k ≈ 1 − e −R k πr/N for such a macroscopic LRN. It follows that LRNs thermalize homogeneously in action space, see e.g. [19].
On the contrary, SRNs show anomalously slow ergodization and thermalization dynamics in proximity to an integrable limit [18,19]. Since R k is now intensive and N -independent, the resonance probability Π k ∼ π r is small, resonances are rare, and thermalization is delayed until resonances were able to migrate through the whole system. Thermalization is expected to be a highly inhomogeneous process in action space.
In this work we quantitatively describe the dynamics of thermalization by making use of finite time average (FTA) distributions. We (a) observe a novel regime of action diffusion, and extract the diffusion coefficient as a function of the proximity to the integrable limit, (b) show the connection between the dynamics of FTA distributions and auto-correlation functions and predict and observe algebraic decay of correlations in time, and (c) finally predict the diffusion delay through action noise destroying resonances and provide computational evidence of the delay.
We consider the Hamiltonian describing the dynamics of a chain of N superconducting islands with nearest neighbor Josephson coupling in its classical limit. This model is equivalent to a 1D XY chain or simply a coupled rotor chain with rotor momenta p n and angles q n . E J controls the strength of Josephson coupling and will be compared to the energy density h = H/N . The equations of motion of Eq. (1) reaḋ We apply periodic boundary conditions p 1 = p N +1 and q 1 = q N +1 . The system has two conserved quantities: the total energy H and the total angular momentum P = N n=1 p n . We will choose P = 0 without loss of generality. A SRN limit is obtained for The actions J n ≡ p n , and the angles Θ l ≡ q n . Note that the opposite limit E J /h → ∞ (not further studied in this work) yields a LRN network as discussed in the introduction.
The microcanonical dynamics of (1) explores the available phase space Γ. The phase space average of an ob- Here X is a point in Γ. The ergodicity property is tested quantitatively by showing that the infinite time average of any observable f ( X) will be equal to its phase space average f . Lacking infinite times, we rather compute finite time averages, which depend on both the averaging time T , and the initial condition X 0 : For an ergodic system it follows for any choice of X 0 except for a subset of measure zero. Dense scanning of all initial points X 0 over Γ yields the finite time average distribution ρ(f ; T ) of the finite time averages f T ( X 0 ). It is a function of f , parametrically depends on T , and is characterized by its moments It follows that the first moment µ 1 ≡ f is invariant under variation of the averging time T . All higher moments will in general depend on T . For an ergodic system it follows In our studies ρ is close to a Gaussian distribution, which allows us to focus on µ 2 . We use the actions (momenta) p n as the relevant slow observables in the SRN proximity to the integrable limit E J → 0: f ≡ p n . With P = 0 it follows µ 1 = 0. The second moment µ 2 (T ) is then simply the variance of ρ, and further related to the momentum-momentum autocorrelation function R(t) = lim τ →∞ 1 τ τ 0 p l (τ )p l (t+τ )dτ as µ 2 (T ) = 1 T T 0 R(t)dt. Under the usual assumption that the correlation function will have an exponential decay at large enough times (with ways to even weaken the requirement) the conclusion is that µ 2 (T → ∞) ∼ 1/T .
The details of the integration methods are outlined in Ref. [18]. We numerically integrate R trajectories using symplectic integrators [21], where each initial point was chosen by setting q n = 0, drawing p n from a Maxwell distribution, constraining P = 0, rescaling all momenta such that the desired energy density h is obtained, and giving each trajectory a prethermalization run of t prethermal = 10 6 . Since all actions p n are statistically equivalent, we measure them all and add all data into one pool which is used to compute ρ. Fig. 1 shows µ 2 (T ) for h = 1, N = 1024 and 0.25 ≤ E J ≤ 1. The variance µ 2 (T ) resists decay up to some characteristic ergodization time scale T E , after which it turns decaying as expected, signalling restoration of ergodicity. Note that this time scale T E was assessed in Ref. [18] and is an intensive time scale.
A close inspection of the size dependence of µ 2 (T ) is shown in Fig. 2 for h = 1, E J = 0.7 and a variety of system sizes. We find that µ 2 (T E ≤ T ≤ T D ) loosely E J =0.7 follows a 1/ √ T diffusive decay, which is followed by the anticipated 1/T decay for T D ≤ T . The new time scale T D (N ) is evidently system-size dependent. To support that finding, we plot δ = d(log 10 µ 2 )/d(log 10 T ) versus T in Fig. 3. The curves clearly show intermediate saturation on a plateau with δ ≈ −0.5, and a subsequent decay at T D (N ) down to δ = −1. Since T D is increasing with system size, we conjecture that T D (N → ∞) → ∞, extending the 1/ √ T decay in µ 2 (T ) to infinite times for infinite size. In turn this implies that the correlation function R(t) ∼ 1/ √ t without any exponential cutoff in the same limit.
Let us discuss possible mechanisms leading to the observed behavior of µ 2 (T ) for different system sizes. We consider the occurence of a chaotic resonance to take place when first order perturbation theory for the evolution of a given rotor at site n breaks down. A simple calculation yields ∆ + n < E J and ∆ − n < E J with ∆ ± n = |p n (p n − p n±1 )|. The presence of such resonantly coupled triplets of grains along the network generate chaotic dynamics and results in a Lyapunov exponent whose inverse yields a time scale T Λ ≤ 10 on the studied interval 0.25 ≤ E J ≤ 1 [18]. It follows that T Λ T E , T D . The resonance probability can be easily computed as π r ∼ (E J /h) 2 . At variance to the LRN cases, the SRN resonances are rare and inhomogeneously distributed over the system at any time. The typical distance between consecutive chaotic triplets l r ∼ (h/E J ) 2 grows with reducing E J turning the resonances more sparse and rare [22]. The assumption of partial thermalization of actions involved in the rare resonances will still not lead to any substantial observation of the onset of ergodization, simply because resonances are rare and separated by non-chaotic (regular) regions. To onset ergodicity instead, chaotic resonances have to diffusively migrate throughout the entire system [18]. This presum-ably happens due to an incoherent detuning of the momenta of rotors in a neighbourhood of a given resonance. Once such a neighbouring rotor is sufficiently detuned, it could become resonant with its own neighbourhood forming a new resonance.
To confirm that we observe action diffusion, we rescale µ 2 → µ 2 N and T → T /N 2 -as shown in Fig.4. For a given value of E J we observe very good collapse of all curves onto one master curve for T ≥ T E . The master   Fig.4 show the turnover from diffusive 1/ √ T to asymptotic 1/T decay at the time T D . The diffusion process assumes that a diffusion coefficient D ∼ N 2 /T D can be read off a fit of T D . We found the inverse D −1 from the intersection of the fit of the diffusive 1/ √ T and the asymptotic 1/T trends (marked with green dots in Fig 4). The measured values of the diffusion coefficient D are then reported as a function of E J in Fig.5, which appears to be reasonably close to a power-law over the analyzed interval. To check whether the asymptotic behavior of D may turn into an exponential behavior rather than power law [23][24][25][26][27] we plot the data in Linear-Log scale in the inset of Fig.5. From the available data we conclude that a power law is more close to the obtained data.

curves in
Our analysis shows that the dynamics of resonances starts with a diffusion process between chaotic triplets, i.e. on a length scale l r ∼ √ DT E . After that the diffusion continues until all fluctuations stored in N/l r nonresonant patches each of the size l r , were exchanged and reached a given location in the system. This happens for T D ∼ N 2 /D [28].
If the above scenario is correct, we can expect to delay the diffusion, relaxation, and ergodization process, if we manage to efficiently destroy resonances, before they had time to diffuse. To check this prediction, we take a system with N = 512 and E j = 0.7 at h = 1. Every time interval T Λ ≈ 10 we randomly pick a site n, and increment or decrement its momentum p n by a given value ∆ p with equal probability [29]. On average a given site n 0 is reached on a time T Λ N ≈ 5000. If ∆ p ≈ E J /h we expect to efficiently detune and destroy a nonlinear resonance. The effect should become visible for T ≈ 5000. The results in Fig.6 are excellently reproducing the prediction. The kicks will also generate new resonances at another location, such that the average number of resonances will not change. Our results confirm that resonance diffusion is at the origin of the ergodization process. It is exactly this diffusive process which is efficiently harmed, destroyed and delayed by the above kicking procedure.
Combining our results with previous studies shows that weakly nonintegrable many-body systems can restore ergodicity in distinctive ways depending on the range of the interaction network in action space. It all starts with action resonances seeding chaotic dynamics into the networks. While long range networks provide well connected resonances with ergodization controlled by the characteristic individual resonance chaos time scales, short range networks instead yield a dramatic slowing down of ergodization in action space, and lead to rare resonance diffusion. We used Josephson junction chains as a paradigmatic study case and exploited finite time average dis- tributions to characterize the thermalizing dynamics of actions. The slowing down of the thermalization dynamics upon approaching the integrable limit results in a decreasing of an effective diffusion constant which is related to heat conductivity. This slowing down appears to follow a power law in the distance from the integrable limit, rather than an exponential one. We identify a novel action resonance diffusion regime responsible for the slowing down. The observed fragile diffusion is relying on weakly chaotic dynamics in spatially isolated action resonances. We were able successfully delay and suppress it by adding weak action noise, as a proof of concept. Among a number of intriguing open questions, we mention the search for further distinct classes of nonintegrable action networks (neither short nor long ranged), and the impact of quantization on the fragile short range network dynamics in the vicintiy of an integrable limit.