# Large-deviations approach to thermalization: the case of harmonic chains
with conservative noise

Stefano Lepri 1,2 1 Consiglio Nazionale delle Ricerche, Istituto dei Sistemi
Complessi, via Madonna del Piano 10, I-50019 Sesto Fiorentino, Italy
2Istituto Nazionale di Fisica Nucleare, Sezione di Firenze, via G. Sansone 1
I-50019, Sesto Fiorentino, Italy<EMAIL_ADDRESS>

###### Abstract

We investigate the possibility of characterizing the different thermalization
pathways through a large-deviation approach. Specifically, we consider clean,
disordered and quasi-periodic harmonic chains under energy and momentum-
conserving noise. For their associated master equations, describing the
dynamics of normal modes energies, we compute the fluctuations of activity and
dynamical entropy in the corresponding biased ensembles. First-order dynamical
phase transition are found that originates from different activity regions in
action space. At the transitions, the steady-state in the biased ensembles
changes from extended to localized, yielding a kind of condensation in normal-
modes space. For the disordered and quasi-periodic models, we argue that the
phase-diagram has a critical point at a finite value of the disorder or
potential strength.

Keywords: Large deviations in non-equilibrium systems; Thermalization

##  1 Introduction

Thermalization, the process by which a system reaches a state of thermal
equilibrium with its surroundings, stands as one of the fundamental phenomena
governing the behavior of diverse physical systems. From quantum particles in
condensed matter systems to the dynamics of stars in astrophysical
environments, the concept of thermalization underpins our understanding of
equilibrium states and relaxation processes across various scales.
Particularly fascinating is the key role of relaxation timescales on the very
possibility of observing signature of the past events and to store and
retrieve memory of them [1].

On general grounds, the problem is usually approached from two main sides. In
the context of nonlinear dynamics, emphasis is on how energy spreads from a
particular point or region in phase space towards the largest portion of it
where energy equipartition among all degrees of freedom occurs. This is the
familiar approach pioneered by Fermi, Pasta, Ulam and Tsingou, which has a
long-standing history [2] but is still actively investigated. Closeness to
integrability usually affects strongly the equilibration in a non-trivial way
[3, 4, 5, 6, 7, 8, 9]. The situation is even more intriguing for systems with
quenched disorder, where and interplay between Anderson localization and
chaotic diffusion is at work [10, 11, 12, 13, 14].

The second approach is more statistical, in the sense that it rests on some
form of stochastic dynamics and on how relatively unlikely, low-entropy states
evolve towards maximal entropy one. This á la Boltzmann approach involves
equations for distributions derived under suitable assumptions. For instance,
in the case of nonlinear lattices this leads to Boltzmann-Peierls [15] and
wave turbulence equations [16] . For one-dimensional chains this turned out to
be useful to understand anomalous transport [17, 18, 19, 20] and the effect of
nonlinear resonances on (pre)thermalization [21, 22].

To calculate the probability of these trajectories and how atypical ones arise
dynamically one can rely on large deviation theory [23, 24, 25]. In the last
decades, such approach has rapidly grown as an effective tool to study
fluctuations statistical mechanics and stochastic processes. For dynamics, it
provides a mean to compute the probability of untypical trajectories and
observables that depart significantly from their average, typical ones.
Dynamical phase transitions may manifest as non-analiticities of the large-
deviation functions.

Our aim here is to explore the possibility to use this approach to explore the
thermalization dynamics of many-body classical systems. The idea would be to
quantify and characterize regions in phase space that yield possibly slow, or
at least sizeably different from average, relaxation pathways. In principle,
this should provide insight on the different phase-space trajectories followed
during equilibration.

A similar idea has been discussed in [26] (see also the account in [27]. There
it has been shown that the large deviation of an observable may have
singularities even in the context of the Gaussian model, a paradigmatic non-
interacting system, before and after an instantaneous temperature quench. The
associated phenomenon of condensation of fluctuations is thus a characteristic
feature of the thermalization process. Remarkably, condensation originates
from the duality between large deviation events in the given system and
typical events in a ”biased” one. We also mention a conceptually related work
on large deviations for wave-turbulence equations [28]. We mostly have in mind
the relaxation of isolated systems: as a matter of facts, large-deviation
approaches have been reported for nonequilibrium steady states close to what
we will deal here, like harmonic networks in contact with thermal baths (see
e.g [29, 30, 31] and references therein).

Since dealing with genuine non-linear forces is a formidable task, it is
convenient to adopt an intermediate ”mesoscopic” approach whereby
nonintegrable interactions are replaced by a random noise. This idea has
proved to be useful, for instance, to study the Lyapunov time-scale [32, 33]
and the onset of diffusive dynamics in almost-integrable systems [34]. Also,
such approach reproduces many features of nonlinear lattices in steady
nonequilibrium states [35, 36]. In this spirit, we focus on a specific model:
an harmonic network with momentum and energy conserving noise that we studied
in [37]. Such model has a the remarkable advantage that, under suitable
approximations, it allows to write explicitly a linear master equation
describing energy spreading in action space. We will exploit this feature to
compute the large deviation properties exactly, through standard methods.

The paper is organized as follows. In Section 2 we present our model system:
an harmonic network subject to momentum and energy conserving noise, which is
the main object of study. We first briefly recall the kinetic approach
presented in [37] that allows to write a master equation for diffusion in
action space which is the starting point of the subsequent analysis. In
Section 3 we briefly recall the methods used to compute the cumulant
generating function of the main observables, the dynamical activity and the
dynamical entropy. The first results concerns the simplest case of a chain
with translational invariance discussed in Section 4. In the second part we
examine the thermalization and large deviation properties of two paradigmatic
models, the disordered (Section 5) and quasi-periodic chain (Section 6) that
correspond, respectively, to the famous Anderson and Aubry-André -Harper
models in the context of electronic transport.

##  2 Continuous-time master equation

Let us first consider a general master equation for the probabilities
$P_{\nu}$ of a state, $\nu=1,\ldots,N$ with transition rates
$R_{\nu,{\nu^{\prime}}}$

$\dot{P}_{\nu}=\sum_{\nu^{\prime}}\left[R_{\nu,{\nu^{\prime}}}P_{\nu^{\prime}}-R_{{\nu^{\prime}},\nu}P_{\nu}\right].$
(1)

Eq. (1) can be rewritten as

$\dot{P}_{\nu}=\sum_{\nu^{\prime}\neq\nu}W_{\nu,{\nu^{\prime}}}P_{\nu^{\prime}}\qquad
W_{\nu,{\nu^{\prime}}}=R_{{\nu^{\prime}},\nu}-r_{\nu}\delta_{\nu,{\nu^{\prime}}}$
(2)

where we introduced the escape rate from state $\nu$ as

$r_{\nu}=\sum_{\nu^{\prime}\neq\nu}R_{{\nu^{\prime}},\nu}$ (3)

and $1/r_{\nu}$ the mean residence time in state $\nu$. The calculation of the
large-deviation properties for (1) is well established and will be briefly
reviewed in the next Section.

Before this, we now specify the physical context that we will deal with. Let
us consider the following quadratic Hamiltonian,

$H=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m_{i}}+\frac{1}{2}\sum_{i,j=1}^{N}q_{i}\Phi_{ij}q_{j},$
(4)

where $(q_{i},p_{i})$, $i=1,\ldots N$ are the canonically-conjugated
coordinates ad momenta (for simplicity, we will deal henceforth with the
equal-mass case, and set $m_{i}=1$). The coupling matrix $\Phi$ is semi-
positive definite and symmetric so that its eigenvalues are real and non-
negative. The equations of motion of the isolated network are

$\ddot{q}_{i}=-\sum_{j}\Phi_{ij}q_{j}\quad.$ (5)

In order to have momentum conservation it must be $\sum_{i}\Phi_{ij}=0$. In
this general formulation, the model encompasses several different setups. It
includes the standard case of periodic lattices, where $\Phi$ is the familiar
nearest-neighbor Laplacian matrix, but also the case of disordered structures
with random and or sparse connections as on the case of elastic networks [38,
39]. It can be regarded as an idealized model of, say, a large atomic cluster
or a protein in its native state.

The normal modes’ coordinates $(\mathcal{Q}_{\nu},\mathcal{P}_{\nu})$ of the
network have eigenfrequencies $\omega_{\nu}$ and are defined by the
transformation

$\mathcal{Q}_{\nu}=\sum_{l=1}^{N}q_{l}\,\chi_{l}^{\nu},\quad\mathcal{P}_{\nu}=\sum_{l=1}^{N}p_{l}\,\chi_{l}^{\nu}.$
(6)

The $\chi_{l}^{\nu}$, $\nu=1,\ldots,N$, are orthonormal and can be taken to be
real for the time being. The main observables are the mode energies
$\mathcal{E}_{\nu}=(\mathcal{P}_{\nu}^{2}+\omega_{\nu}^{2}\mathcal{Q}_{\nu}^{2})/2$

The classic problem amounts to ask how a (possibly weak) nonlinear interaction
potential added to to $H$ brings the harmonic modes to equipartition. For
conservative dynamics, upon expressing the Hamiltonian in the familiar action-
angle variables, a nonintegrable perturbation defines a network of couplings
among the unperturbed actions [40, 41]. The connectivity of such network will
determine relaxation and ergodic properties. In general, one can distinguish
such networks depending on whether the number of groups of actions linked by
the perturbation depends intensively or extensively on the number of degrees
of freedom $N$. These are termed short or long-range networks, respectively,
depending on whether the coupling range if roughly constant or increases
proportionally to $N$ [40, 41].

In the nonlinear case the coupling may involve three or more actions,
resulting in complicated nonlinear interaction on a hyper-graph structure. In
view of the difficulty of the problem, it is of interest to formulate a hybrid
dynamics for Hamiltonian (4) defined adding on top of the linear deterministic
dynamics a conservative stochastic process [42] defined as follows. Consider
the network state at time $t$. At a later time $t+\tau$ there occur a
”collision” event whereby a couple of particles $(m,n)$, $m\neq n$ is randomly
selected according to the joint probability $w_{n,m}$ and their momenta are
exchanged, $(p_{n},p_{m})\to(p_{m},p_{n})$. This move conserves energy and
momentum. Physically, it can be interpreted as a perfectly elastic rebound as
it would occur for a infinite square-well pair potential. The $w_{n,m}$ is
usually assigned to describe local interactions (e.g. nearest-neighbors) but
may include long-range terms interaction across the lattice as well. The
intervals between subsequent collision times $\tau$ are also taken as random
variables with some preassigned distribution with given, finite, average
$\langle\tau\rangle$.

In the kinetic approximation [37] the (angle-averaged) energy
$E_{\nu}=\overline{\mathcal{E}_{\nu}}$ of each normal mode obeys the equation

$\dot{E}_{\nu}=\sum_{\nu^{\prime}\neq\nu}\left[R_{\nu^{\prime},{\nu}}E_{\nu^{\prime}}-R_{{\nu^{\prime}},\nu}E_{\nu}\right]$
(7)

This can be regarded as an averaged collision operator coarse-grained over a
time scale which is longer than the typical collision time. We work with
$\gamma=\frac{1}{N\langle\tau\rangle}$ is finite, corresponds to choosing a
finite collision probability per site (as usually employed in the simulations
[36, 43, 44, 45]). The form, of equation (7) guarantees thar
$\sum_{\nu}E_{\nu}$ remains constant as it should. One major advantage is
that, the simple form of the collision rule allows to compute explicitely the
transition rates in terms of the eigenvectors [37]. First, define $V$ as the a
unit vector of components

$V_{\nu}=V_{\nu}^{(n,m)}\equiv\frac{\chi_{n}^{\nu}-\chi_{m}^{\nu}}{\sqrt{2}},$
(8)

Note that $V=V^{(n,m)}$ is a random vector as the indexes $n,m$ are chosen at
random with the prescribed rule. Then its found that

$R_{\nu,{\nu^{\prime}}}=\frac{2}{\langle\tau\rangle}\sum_{(n,m)}(V_{\nu}^{(n,m)})^{2}\,w_{n,m}(V_{\nu^{\prime}}^{(n,m)})^{2}\equiv\frac{2}{\langle\tau\rangle}\overline{V_{\nu}^{2}V_{\nu^{\prime}}^{2}}$
(9)

where the overline is a shorthand notation for the average over $w_{n,m}$ of
the random collisions. Note that
$R_{\nu,{\nu^{\prime}}}=R_{\nu^{\prime},{\nu}}$ with the escape rate from
state $\nu$ as

$r_{\nu}=\sum_{\nu^{\prime}\neq\nu}R_{{\nu^{\prime}},\nu}=\frac{2}{\langle\tau\rangle}\overline{(1-V_{\nu}^{2})V_{\nu}^{2}}.$
(10)

For a generic coupling matrix $\Phi$, such that there is no decoupling among
different subsets they are all nonvanishing. They are of the form of the well-
known Fermi golden rule, involving the squared amplitudes of the eigenmodes.

Formally, the quantities $P_{\nu}=E_{\nu}/\sum E_{\nu}$ (comprised between 0
and 1) satisfy (7). A physical interpretation would be to imagine an large
ensemble of quasiparticles with a given distribution among modes, each
$P_{\nu}$ being the istantaneous fraction of them in mode $\nu$. Altogether,
the dynamics is a random walk of quasiparticles in action space. Transitions
are determined by the eigenstate structure of the coupling matrix $\Phi$, and
not assigned a-priori.

Following the usual arguments of Markov processes, the system is ergodic and
approaches microcanonical equilibrium where $P_{\nu}=1/N$ (equipartition)
according to the properties of the ”collision operator” defined by (1) and
(16). So the thermalization problem reduces to computing its $N$ eigenvalues
$\mu_{\nu}$, $\nu=1,\ldots,N$. whose absolute values give the spectrum of
relaxation rates. The first eigenvalue is $\mu_{1}=0$ and corresponds, as
said, to the steady state of energy equipartition. On physical grounds,
$\mu_{1}$ must be non-degenerate since we expect the dynamics to be ergodic
under the above hypoteses. All the others $\mu_{\nu}$ must be strictly
negative [46]. In particular, the long-time relaxation rate is controlled by
spectral gap $|\mu_{2}-\mu_{1}|=|\mu_{2}|$, $\mu_{2}$ being referred to as
Fiedler eigenvalue in the context of diffusion on graphs. For a finite
network, the scaling with $N$ of $\mu_{2}$ can be expressed in terms of the
scaling of the spectral density. More precisely, denoting by $\rho(\mu)$ the
integrated (cumulative) density of eigenvalues (i.e. the fraction of
eigenvalues less than $\mu$), if

$\rho(\mu)\sim|\mu|^{d_{\mu}/2},\qquad\mu\to 0$ (11)

then $d_{\mu}$ is the spectral dimension of the action network. We emphasize
that this is a distinct quantity with the standard spectral dimension of the
Laplacian matrix $\Phi$ entering in the spectral problem for the Hamiltonian
(4).

The solution of master equation with initial condition $P_{\nu}(0)$, is fairly
straightforward. Upon expanding in the basis of the eigenvectors
$\psi^{(\nu)}$ of the matrix $W$ (the stochastic matrix is symmetric so there
is no need to distinguish left and right eigenvectors), one can express the
solution as

$P_{\mu}(t)=\sum_{\nu\nu^{\prime}}e^{\mu_{\nu}t}\psi^{(\nu)}_{\mu}\psi^{(\nu)}_{\nu^{\prime}}P_{\nu^{\prime}}(0).$
(12)

This can be easily solved numerically to compute for the observable of
interest.

##  3 Large-deviations approach

In this section we briefly review the methods employed in the following. The
material is standard in the related literature, and we mostly follow reference
[47].

### Activity

One quantity studied in the literature is the activity, which counts the
number of jumps in a trajectory i.e. the number $K$ of configuration changes
during a time $t$. As the name implies, it is a measure of how active a given
trajectory is: an active trajectory has many changes of configuration, an
inactive one just few or none.

The relevant operator to study is the tilted generator which for the activity
is defined by the $N\times N$ matrix with elements [48, 49]

$\left(W_{K}(s)\right)_{{\nu^{\prime}},\nu}\equiv
e^{-s}R_{{\nu^{\prime}},\nu}-r_{\nu}\,\delta_{\nu,{\nu^{\prime}}}.$ (13)

Note that his operator has $-r_{\nu}$ on the diagonal and the first terms
appears only in off-diagonal entries (see for example eqs. (44) and (71) in
Ref. [47]).

Following the general strategy, the large-deviation statistics is encoded in
the leading eigenvalue and (right and left) eigenvectors of the tilted
operator (also called scaled cumulant generating function). More precisely, we
denote by $\lambda_{K}(s,N)$ the largest eigenvalue of $W_{K}(s)$ along with
its the right eigenvector defined by $W_{K}(s)\phi_{K}=\lambda_{K}\phi_{K}$.
Upon normalization, $\phi_{K}$ can be considered as a steady probability
distribution for the ”biased” dynamics (sometimes termed $s$-ensemble). As it
is known [47] for $s=0$ the eigenvector yields the density in the steady
state, which in our case is the equipartition which is homogeneous. For $s<0$,
it probed the regime in which the mean activity $K/t$ of histories is
typically larger than in the steady state. They correspond to explored
configurations where the density is larger than the steady state distribution.
On the other hand, at $s>0$ histories with smaller $K/t$ are favored.

As an order parameter we may thus consider (see eq. (177) in [47]).

$\rho_{K}(s)=\frac{1}{N}\sum_{\nu=1}^{N}\nu{\phi_{K,\nu}}(s)$ (14)

that, in our context, measure the occupancy of the steady state and thus of
the distribution of modes that are involved in the biased dynamics.

### Entropy production

Another observable of interest is associated to dynamical complexity,
$Q_{+}=\ln Prob(history)$. In the dynamical systems literature the
corresponding large-deviation function is known as the topological pressure.
It can be computed in a similar way through the operator $W_{+}(s)$ (see e.g.
Eq. (36) in Ref. [47])

$\left(W_{+}(s)\right)_{{\nu^{\prime}},\nu}\equiv
R_{{\nu^{\prime}},\nu}^{1-s}\,r_{\nu^{\prime}}^{s}-r_{\nu}\,\delta_{\nu,{\nu^{\prime}}}$
(15)

and its the maximal eigenvalue is the topological pressure, large deviation of
the entropy production. Following the same prescription as above, we seek for
the leading eigenvalue $\lambda_{+}(s,N)$ with right eigenvector $\phi_{+}$,
$W_{+}\phi_{+}=\lambda_{+}\phi_{+}$. An order parameter $\rho_{+}(s)$ can be
defined as in (14) with $\phi_{+,\nu}$ replacing $\phi_{K,\nu}$. The average
$Q_{+}/t$ is the Kolmogorov-Sinai entropy $h_{KS}=d\lambda_{+}/ds(0)$; another
important quantity is the topological entropy $h_{top}=\lambda_{+}(1)$ which
is the growth rate of the number of possible histories with time.

##  4 Translationally-invariant chains

Let us start from one-dimensional chains with translational invariance i.e.
the case where $\Phi$ is a circulant matrix. The simplest example would be the
standard nearest-neighbor coupling, but the results apply to any circulant
$\Phi$, including the case of harmonic long-range interactions [50]. In fact,
in all those instances, the eigenvectors are the familiar lattice Fourier
modes and the vector $V$ can be given explicitly

$\chi_{l}^{\nu}=\frac{1}{\sqrt{N}}e^{ik_{\nu}l};\quad
V_{\nu}^{(n,m)}\equiv\frac{e^{ik_{\nu}n}-e^{ik_{\nu}m}}{\sqrt{2N}},$

where $k_{\nu}=\frac{2\pi\nu}{N}$ (we use the labeling $\nu=-N/2+1\ldots N/2$
throughout this Section). Thus the transition rates in formula (9) are known
exactly once $w_{n,m}$ is given. Furthermore, for any choice of $\Phi$, the
rates only depend on $w_{n,m}$. For definiteness, we will focus on the case of
nearest-neighbor random exchange, namely $m=n+1$. This is the mostly studied
choice in the literature both for short (nearest-neighbor) [42, 51] than for
long-range interactions [52]. In this case, the transition rates can be
computed explicitely since [37]

$|V_{\nu}^{(n,m)}|^{2}=\frac{1}{2N}|1-e^{ik_{\nu}}|^{2}=\frac{2}{N}\sin^{2}\frac{k_{\nu}}{2},$

so that

$R_{\nu,\nu^{\prime}}=\frac{8\gamma}{N}\sin^{2}\frac{k_{\nu}}{2}\sin^{2}\frac{k_{\nu^{\prime}}}{2}$
(16)
$r_{\nu}=\sum_{\nu^{\prime}\neq\nu}R_{{\nu^{\prime}},\nu}=4\gamma\sin^{2}\frac{k_{\nu}}{2}[1-\frac{2}{N}\sin^{2}\frac{k_{\nu}}{2}]$
(17)

Note that the escape rate $r_{\nu}$ is different for each state $\nu$. Another
interesting property is that the off-diagonal terms are of order $1/N$, and
thus the eigenvalues of the matrix $W$ are well approximated by the diagonal
elements

$\mu_{\nu}\approx-4\gamma\sin^{2}\left(\frac{\pi(\nu-1)}{N}\right),$ (18)

and the $\nu$th eigenvector is localized on $\pm\nu$ with all other components
being small of order $1/N$. The relaxation rate of a generic non-equilibrium
initial condition occurs, at long enough times, with a rate
$|\mu_{2}|\approx\gamma(\pi/N)^{2}$. Thus, the matrices $W_{K}$ and $W_{+}$
are also known analytically: for instance the off-diagonal elements of $W_{+}$
are

$(W_{+})_{\nu^{\prime}\neq\nu}=R_{{\nu^{\prime}},\nu}^{1-s}\,r_{\nu^{\prime}}^{s}=\gamma\frac{2^{2s+3}}{N^{1-s}}\sin^{2}\frac{k_{\nu}}{2}\left(\sin^{2}\frac{k_{\nu^{\prime}}}{2}\right)^{1-s}.$
(19)

Figure 1: Large deviations of activity for the translationally invariant chain
with conservative noise, $\gamma=1$. (a) the largest eigenvalue of $W_{K}(s)$
as a function of $s$, for different chain lengths $N$; the inset in an
enlargement around $s=0$ showing the eigenvalues multiplied by $N^{2}$ ; (b)
the order parameter $\rho_{K}(s)$; (c) the eigenvectors $\phi_{K}$ for $N=256$
as a function of the wavenumber $k_{\nu}$: left panel
$s-0.005,-0.01,-0.0500.1$, middle $s=0$, right $s=0.1$. The data show evidence
of a dynamical phase transition at $s=0$ at which the activity displays a
first-order jump and $\phi_{K}$ changes from extended to localized.

Let us first discuss the case of the activity, In Figure 1a we report the
scaled large deviation function $\lambda_{K}(s)$ for different sizes, obtained
computing the leading eigenvector numerically. It is shown that that large
deviation function has a discontinuity in the derivative at $s=0$
corresponding to a dynamical first-order phase transition. Indeed, increasing
$s$ leads to a sudden jump in the typical activity $d\lambda_{K}/ds$, which
corresponds to a dramatic change in the number of configurations explored by
histories with activity rate $K/t$. As far as the system size scaling is
concerned, for $s<0$ $\lambda_{K}$ is $N-$ independent, while for positive $s$
it vanishes as $1/N^{2}$ (see inset in Fig. 1a). This finite-size analysis
provide a convincing evidence that the discontinuity persists in the
thermodynamics limit.

Figure 2: Large deviation function for the dynamical entropy for the
translationally invariant chain with conservative noise, $\gamma=1$. (a) The
largest eigenvalue of $W_{+}(s)$ as a function of $s$ for different chain
lengths $N$; the inset in an enlargement around $s=0$ showing the eigenvalues
divided by $\ln N$ and the dashed line has a slope $2\gamma$; (b) the
parameter $\rho_{+}(s)$; (c) the eigenvectors $\phi_{+}$ for $N=256$ as a
function of the wavenumber $k_{\nu}$: left panel $s=-0.1$, middle $s=0$, right
$s=0.005,0.01,0.050,0.1$ .

Taking into account the fact that plane waves come in pairs with opposite
wavenumber, it is apporpriate to adapt the definition of the order parameter
(14) as

$\rho_{K}(s)=\frac{1}{N}\sum_{\nu}|k_{\nu}|\phi_{K,\nu}(s)$ (20)

(and the same for $\rho_{+}(s)$). This gives indication about the typical
wavenumbers which are mostly active along a biased trajectory. Figure 1b shows
that $\rho_{K}$ has a discontinuos jump at the transition point $s=0$,
indicating that different subset of degrees of freedom are involved in the
biased dynamics. This is confirmed by inspecting the structure of the leading
eigenvector $\phi_{K}$, which is markedly different it the two regimes, see
Figure 1c:

  * •

at $s=0$ the eigenvector $\phi_{K}$ is uniform, as expected since it must
correspond to the steady state of energy equipartition;

  * •

for $s>0$, $\phi_{K}$ has two large components at $\pm k_{1}$ plus a uniform
background of size $O(1/N)$; moreover, this shape it is independent of $s$;

  * •

for $s<0$, $\phi_{K}$ is very small at the band center $k_{\nu}=0$ were it
attains a minimum and increasingly concentrates toward zone-boundary Fourier
modes $k_{\nu}\approx\pm\pi$ upon decreasing $s$.

Thus, trajectories in the inactive phase $s>0$ consists mostly the longest-
wavelength mode which remains excited, undergoing very little exchange with
others. On the other hand, the active phase $s<0$ is mostly peaked on short-
wavelengths and involve rapid energy exchange among a finite fraction of them.
This gives an overall picture of action space, which is consistent with the
numerical simulation of the microscopic dynamics [37]. The remarkable fact
that the transition occurs at $s=0$ means that the model dynamics displays
”phase coexistence” between fast and slow relaxing modes.

To complement the analysis, we also performed a calculation of the large
deviation for the dynamical entropy. Figure 2 reports the same type of
analysis for the topological entropy. There is still a similar transition at
$s=0$, separating two regimes: one ”regular” with zero KS entropy and one
”chaotic” with finite (perhaps $N$ dependent) entropy. So the natural dynamics
present a coexistence of these two. Since the standard KS entropy $k_{KS}$ is
given by the the slope at the origin, from the data in the inset of Figure 2a
it is seen that for $s\to 0^{+}$, $k_{KS}\approx 2\gamma\log N$, meaning that
it subextensive in the system size. It is worth noticing that such scaling
occurs in other instances, like for instance the infinite-range Ising model
(see Section 6 in [47]). This may be related to the fact that the coupling
among modes is of mean-field type, see again eqs. (16). This feature is
different from genuinely chaotic dynamical systems, where typically the KS
entropy is extensive, $h_{KS}\sim N$.

At this point, we remark that there is indeed a similarity with kinetically
constrained models, as Fredrickson-Andersen dynamics and the like, which have
been thoroughly studied as toy models of glasses [48, 25]. Also there, the
large-deviation function is singular at $s=0$ and the order parameter has a
first-order jump. The interpretation is that there has two phases, an active
one for $s<0$ and an inactive one for $s>0$. Physical dynamics take place at
$s=0$, where the two dynamic phases coexist. It our case the same
interpretation applies in action rather than in physical space.

##  5 Disordered chain

Let us now turn to the simplest non-homogeneous lattice: the harmonic chain
with disorder in the pinning potential

$H=\sum_{i=1}^{N}\left[{p_{i}\over
2}^{2}+\frac{1}{2}(1+\sigma_{i})q_{i}^{2}+\frac{1}{2}(q_{i+1}-q_{i})^{2}\right]$
(21)

with $\sigma_{i}$ being i.i.d. variables with uniform distributions in
$[0,w]$, gauged by the disorder strength parameter $w$. Periodic boundary
conditions are assumed. As it is well known, the eigenstates are the
exponentially-localized Anderson modes, whose localization length decreases
with $w$ [53, 54]. Several variants in this type of model with conservative
noise have been considered earlier [55, 56, 57]. We consider the standard case
of nearest-neighbor collisions and individual realizations of the quenched
disorder $\\{\sigma_{i}\\}$.

Figure 3: Thermalization of the disordered chain with conservative noise
$\gamma=1$: (a) the spectral gap as a function of the disorder strength $w$,
the inset shows that $\mu_{2}$ is proportional to $1/N^{2}$ for $w$ large;
(b,c) time evolution of the spectral entropy $S(t)$ for two values of $w$
corresponding to weak and strong disorder respectively. Data from numerical
solution of the master equation starting from the ground-state mode
$\nu=1$.The dashed line is the law $S(t)=\frac{1}{2}\log t+const.$.

Before discussing large deviations, we illustrate some features that can be
inferred from the kinetic equations (see also [37] for a comparison with
numerics). The localization level affect drastically the connectivity in
action space as measured by the transition rates. Indeed, the
$R_{\nu,\nu^{\prime}}$ are all non-vanishing, but can be exponentially small
due to the small overlap of the corresponding eigenvectors, see again Eq. (16
). In Figure 3a we plot the spectral gap as a function of disorder. There is a
size-dependent disorder value, separating fast relaxation regime from a slower
one. Qualitatively, this correspond to a change in the structure of the matrix
$W$ from fully connected (as in the ordered case) to almost-banded. Above such
value, the spectral gap vanishes as $1/N^{2}$, see the inset in Figure 3a.

A key indicator that has been often employed to characterize thermalization is
the non-equilibrium (Shannon) spectral entropy [58]

$S(t)=-\sum_{\nu}P_{\nu}\log P_{\nu}.$ (22)

For large $t$, $S$ approaches its equipartition value $\log N$. In general $S$
will depend on the initial conditions. In Figure 3b,c we report the evolution
of $S(t)$ for two values of the disorder strength and different sizes obtained
by numerical solution of the master equation, starting from the ground-state
excitation $P_{\nu}(0)=\delta_{\nu,\nu_{0}}$. ($\nu_{0}=1$ in the
calculations). There is a crossover from an initial exponential to an
intermediate logarithmic growth. The initial regime can be accounted for by a
”mean-field” approximation [59]: assuming that energy of the initially excited
mode is evenly redistributed among the other $N-1$ ones, namely
$E_{\nu}\approx(1-E_{\nu_{0}})/(N-1)$ for $\nu\neq{\nu_{0}}$ yielding the
approximation

$S(t)\approx S_{MF}\equiv-E_{\nu_{0}}\log
E_{\nu_{0}}-(1-E_{\nu_{0}})\log\left(\frac{1-E_{\nu_{0}}}{N-1}\right)$ (23)

where $E_{\nu_{0}}(t)\approx(1-1/N)\exp(-\mu_{\nu_{0}}t)+1/N$. This
approximation works well for long-range action networks i.e. for weak
disorder. The second regime is readily understood as a one-dimensional
diffusive process, where the number of excited modes grows as a $\sqrt{t}$.
This is in agreement with the known result that for diffusion on graphs the
Shannon entropy grows as $\frac{d_{s}}{2}\log t$ if the spectrum of the
corresponding Laplacian operator shows scaling with a finite spectral
dimension $d_{s}$ [59]. It also fits with the observed wavepacket spreading
law in nonlinear disordered chains in the strongly chaotic regime [12, 60]
where indeed microscopic chaos can be effectively approximated by a stochastic
process.

Let us now turn to the calculation of the large deviation properties. In
Figure 4a we report data for the activity. There is again a signature of a
first-order dynamical transition as in the translationally-invariant case and
thus of an active and an inactive phase. However, the critical point is now
shifted to a nonzero value $s=s_{c}$. As for the finite-size analysis, note
that the curve collapse upon dividing by $N$, signaling extensivity, at
variance with the ordered case. Moreover, the data in the inset of Figure 4a
show that for $N>10^{3}$ the curves collapse very accurately also close to
$s_{c}$. The order parameter $\rho_{K}(s)$ displays the same qualitative
behavior, with a step discontinuity at $s=s_{c}$. Moreover, we observed that
$s_{c}$ increase upon increasing the disorder strength $w$. The shift of the
critical point implies that the unbiased dynamics at $s=0$ has a finite
activity rate. The structure of the eigenvector $\phi_{K}$ displays a
qualitative change at $s=s_{c}$ from a rather homogeneous profile to a
localized one (Figure 4c).

Figure 4:  Large deviations of activity for the disordered chain with
conservative noise $w=2$, $\gamma=1$: (a) the rescaled largest eigenvalue of
$W_{K}(s)$ as a function of $s$ for different chain lengths $N$; the inset in
an enlargement around $s=0$; (b) the order parameter $\rho_{K}(s)$; (c) the
eigenvectors $\phi_{K}$ for $N=2048$. The data show evidence of a dynamical
phase transition at $s=s_{c}\approx 0.0376$ at which the activity displays a
first-order jump and the leading eigenvector of $W_{K}$ changes from extended
to localized.

In Figure 5, we report the calculation for the entropy. The data confirms that
the dynamical transition is shifted to a finite $s$ where the order parameter
$\rho_{+}(s)$ undergoes a finite jump for large enough $N$. A noteworthy
difference with the respect to the ordered case is that the slope at the
origin is proportional to $N$ (inset of Figure 5), meaning that in this case
there is a finite Kolmogorov-Sinai entropy density $h_{KS}$. This feature
suggest a form of more homogeneous ”chaos” in presence of heterogeneity.
Finally, the shape of the eigenvector $\phi_{+}$ displays a qualitative change
at $s=s_{c}$ from a rather homogeneous profile to a localized one (Figure 4c),
similar to what observed for the activity.

Figure 5: Large deviations of entropy production for the disordered chain with
conservative noise $w=2$ $\gamma=1$: (a) the largest eigenvalue of $W_{+}(s)$
as a function of $s$ for different chain lengths $N$; the inset in an
enlargement around $s=0$ showing the the eigenvalues divided by $N$ the dashed
line has a slope $2\gamma$. (b) the parameter $\rho_{+}(s)$; (c) the
eigenvectors for $N=1024$ as a function of the lattice site number left
$s=-0.1$, middle $s=0$, right $s=0.1$ evidence of transition at
$s\approx-0.059$ at which the eigenvector changes from extended to localized.

What is the dynamical phase diagram of the model? We have seen that the phase
transition persist upon increasing $w$. On the other hand, for large disorder
the analysis in [37] and the data in Figure 3c suggest that, for fixed $N$,
the thermalization basically amounts to a diffusive process on a network with
finite spectral dimension (equal to one). In this limit, the problem should
thus amount to a continuous-time random walk and we do not expect any
dynamical phase transition. Indeed, for a such a process $\lambda_{+}(s)$
should be roughly Poissonian with no singularity [47].

Figure 6: Disordered chain with conservative noise: comparison of the scaled
large deviation functions for two values of the disorder strength $w=2,10$
(left and middle panels respectively). (a,b) scaled largest eigenvalue
$\lambda_{K}$ and (c,d) $\lambda_{+}$. The dashed line in panel (d) is a fit
with a Poissonian law yielding $0.38(e^{5.0s}-1)$; (e) schematic phase diagram
for the activity in the $(s,w)$ plane where the continuous line is the
transition line from active to inactive state. The red dot marks the critical
point where the transition disappear.

Based on the above limiting cases, we may argue that the transition line
$s_{c}(w)$ in the $(s,w)$ plane should terminate at a critical point. In
Figure 6a-d we show $\lambda_{K}(s)$ and $\lambda_{+}(s)$ divided by $N$ for
two different values of the disorder strength. The data confirm the
expectation: for $w=10$ the curves have smooth derivatives and there is no
sign of criticality. In Figure 6e we draw a schematic phase diagram for the
activity based on the above consideration. A similar diagram should be valid
for the entropy.

##  6 Quasi-periodic chain

To further test the approach and the above results, let us consider the chain
with a quasi-periodic on-site potential

$H=\sum_{i=1}^{N}\left[{p_{i}\over
2}^{2}+\frac{1}{2}\left[1+\Delta\cos(2\pi\beta
i+\varphi)\right]q_{i}^{2}+\frac{1}{2}(q_{i+1}-q_{i})^{2}\right].$ (24)

Here, $\Delta>0$ is strength of the local force and the parameter $\beta$
measures the degree of incommensurability of the potential with the lattice
period and $\varphi$ in an arbitrary phase. The eigenstate problem can be
recast in the form

$\epsilon\chi_{l}=\chi_{l+1}+\chi_{l-1}+\Delta\cos(2\pi\beta
l+\varphi)\,\chi_{l}$ (25)

(where $\epsilon=\omega^{2}+1$) which is readily recognized as the celebrated
Aubry-André-Harper model [61, 62, 63]. A well-known feature of the model is
that it displays a localization transition at $\Delta=2$, where the eigenmodes
change from extended to exponentially localized. For a finite, periodic chain
of length $N$ a typical choice is to set $\beta$ as a rational approximant,
for instance the ratio of consecutive Fibonacci numbers, such that for large
$N$ $\beta$ approaches the golden mean.

To compare with the disordered chain, in Figure 7a we report the spectral gap
as a function of $\Delta$. For $\Delta<2$ the eigenvectors are extended and
the spectral gap remains finite and independent of the size. Instead, as soon
as localization sets in for $\Delta>2$, $|\mu_{1}|$ vanishes as $1/N^{2}$ (see
the inset of Figure 7a) and thermalization slows down considerably. The
typical growth of the spectral entropy follow the same two basic routes as
above: $S(t)$ grows exponentially or logarithmically in time for the two cases
respectively (see Figures 7b-c) and saturates at the equipartition value.

Thus, as far as average properties are concerned, the quasiperiodic lattice
behaves very much as the disordered chain, the main difference being that the
’closing’ of the spectral gap (and the change from exponential to logarithmic
growth of $S$) occurs sharply at $\Delta=2$, independently of the lattice
size.

Figure 7: Thermalization of the the quasi-periodic chain with with
conservative noise starting from the ground-state mode $\nu=1$, $\gamma=1$;
(a) the spectral gap as a function of the potential strength $\Delta$, inset
shows that $\mu_{2}$ is proportional to $1/N^{2}$ for $\Delta>2$. (b,c) time
evolution of the spectral entropy $S(t)$ for two values $\Delta<2$ a and
$\Delta>2$ respectively. The dashed line is $\frac{1}{2}\log t$. In all
calculation we set $\beta$ as the ratio of consecutive Fibonacci numbers and
$\varphi=\pi/5$.

As for the large deviation properties, we checked that upon increasing the
strength of the potential $\Delta$ the behavior is definitely similar to the
disordered chain. For instance, the qualitative behavior of $\lambda_{+}(s)$
for $\Delta=0.5$ is the same as that in Figures 5 (data not reported). So we
argue that the overall phase diagram should be similar in the two models, and
schematically of the same form as in Figure (6e) (with $w$ replaced by
$\Delta$).

It is also of interest to study the behavior of the the large deviation
functions across the localization transition $\Delta=2$. In particular the
numerical solution shows that the dynamical phase transition does not occur
around this value. We just discuss the case of the entropy, as the activity
displays similar features. For instance, in Figure 8a,b we report the scaled
$\lambda_{+}$ for two values of $\Delta$ below and above the transition. There
is no sign of discontinuity neither in $\lambda_{+}$ nor in the order
parameter $\rho_{+}$ , as shown in Figures 8c,d which only shows a smooth
monotonic growth in the extended phase $\Delta<2$ (Figure 8c) and a maximum
around $s=0$ in the localized phase $\Delta>2$ (Figure 8d). Another difference
is in the eigenvectors profiles that do not display a sharp change from
extended to localized as in the disordered chain (compare Figure 8e-h with the
lower panels of Figure 5). Actually, for $\Delta<2$, $\phi_{+}$ changes
gradually from a peaked distribution for negative $s$ to a uniform one for
$s>0$, albeit with some structure (Figure 8e-f). Instead, for $\Delta>2$,
$\phi_{+}$ remains extended and fluctuates wildly (Figure 8g-h), displaying a
fluctuating pattern which may reflect the underlying structure of the
eigenfunctions.

Figure 8: Large deviations of dynamical entropy for the quasiperiodic chain
with conservative noise $\gamma=1$ for different sizes and two value of the
potential strength, $\Delta=1.5$ (left columns) and $\Delta=2.5$ (right
columns) ; (a,b) the scaled largest eigenvalue of $W_{+}(s)$ as a function of
$s$ for different chain lengths $N$; (c,d) the order parameter $\rho_{+}(s)$;
(e-h) the eigenvectors $\phi_{+}$ for $N=1597$ as a function of the lattice
site number (e,g) $s=-0.1$ (f,h) $s=0.1$.

##  7 Conclusions

This work is an attempt to employ large-deviations methods to get insights on
the thermalization of classical many-body systems. Exploiting the mathematical
simplicity of conservative noise dynamics, we are able to compute the large
deviation functions for two relevant quantities, the activity and dynamical
entropy. For both observables, there is evidence of a dynamical phase
transitions that correspond to trajectories of different qualitative nature.
We thus distinguish different region in action space corresponding to
different thermalization pathways.

An appealing way to interpret the transition from extended to localized in the
leading eigenvectors, is to describe it as a form of condensation [26, 27]
occurring in action/mode space. In other words, most trajectories (a
macroscopic fraction) in the localized phase concentrate on a subset of action
space. The structure of the eigenvector gives information on the actual modes
mostly involved in the dynamics.

When formulated in terms of the kinetic equations, the thermalization of such
system is basically a random walk on a large network in action space. Every
move conserves the total energy and only changes the entropy and the graph is
undirected. From this point of view, there is a close connection with the
general problem of large deviations of observables of random walks evolving on
random graphs. This issue has been investigated recently [64, 65]. Indeed,
sudden changes in degree fluctuations, similar to dynamical phase transitions,
are related to localization transition are found that may be connected with
what discussed here.

## Acknowledgements

We thank Stefano Iubini for useful discussions and for reading the manuscript.

## References

## References

  * [1] Rovelli C 2022 Entropy 24 1022
  * [2] Gallavotti G 2007 The Fermi-Pasta-Ulam problem: a status report vol 728 (Springer)
  * [3] Benettin G and Ponno A 2011 Journal of Statistical Physics 144 793 ISSN 1572-9613 URL https://doi.org/10.1007/s10955-011-0277-9
  * [4] Benettin G, Christodoulidi H and Ponno A 2013 Journal of Statistical Physics 152 195–212
  * [5] De Roeck W and Huveneers F 2015 Commun. Pure Appl. Math. 68 1532–1568
  * [6] Huveneers F 2017 Annalen der Physik 529 1600384 (Preprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.201600384) URL https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.201600384
  * [7] Fu W, Zhang Y and Zhao H 2019 Phys. Rev. E 100 052102
  * [8] Goldfriend T and Kurchan J 2019 Physical Review E 99 022146
  * [9] Baldovin M, Vulpiani A and Gradenigo G 2021 Journal of Statistical Physics 183 1–16
  * [10] Pikovsky A S and Shepelyansky D L 2008 Physical Review Letters 100 094101 (pages 4) URL http://link.aps.org/abstract/PRL/v100/e094101
  * [11] Kopidakis G, Komineas S, Flach S and Aubry S 2008 Physical Review Letters 100 084103 (pages 4) URL http://link.aps.org/abstract/PRL/v100/e084103
  * [12] Skokos C, Krimer D O, Komineas S and Flach S 2009 Physical Review E - (Statistical, Nonlinear, and Soft Matter Physics) 79 056211 (pages 12) URL http://link.aps.org/abstract/PRE/v79/e056211
  * [13] Lepri S, Schilling R and Aubry S 2010 Physical Review E 82 056602
  * [14] Basko D 2011 Annals of Physics 326 1577–1655
  * [15] Spohn H 2006 Journal of statistical physics 124 1041–1104
  * [16] Onorato M, Lvov Y V, Dematteis G and Chibbaro S 2023 Physics Reports 1040 1–36
  * [17] Pereverzev A 2003 Phys. Rev. E 68 056124
  * [18] Nickel B 2007 J. Phys. A-Math. Gen. 40 1219–1238 ISSN 1751-8113
  * [19] Lukkarinen J and Spohn H 2008 Communications on Pure and Applied Mathematics 61 1753–1786 ISSN 1097-0312
  * [20] Lukkarinen J 2016 Kinetic theory of phonons in weakly anharmonic particle chains Thermal transport in low dimensions (Springer) pp 159–214
  * [21] Onorato M, Vozella L, Proment D and Lvov Y V 2015 Proceedings of the National Academy of Sciences 112 4208–4213 ISSN 0027-8424 (Preprint https://www.pnas.org/content/112/14/4208.full.pdf) URL https://www.pnas.org/content/112/14/4208
  * [22] Huveneers F and Lukkarinen J 2020 Physical review research 2 022034
  * [23] Touchette H 2009 Physics Reports 478 1–69
  * [24] Touchette H 2018 Physica A: Statistical Mechanics and its Applications 504 5–19
  * [25] Jack R L 2020 The European Physical Journal B 93 1–22
  * [26] Zannetti M, Corberi F and Gonnella G 2014 Physical Review E 90 012143
  * [27] Corberi F and Sarracino A 2019 Entropy 21 312
  * [28] Guioth J, Bouchet F and Eyink G L 2022 Journal of Statistical Physics 189 20
  * [29] Saito K and Dhar A 2011 Physical Review E 83 041121
  * [30] Fogedby H C and Imparato A 2012 Journal of Statistical Mechanics: Theory and Experiment 2012 P04005
  * [31] Jakšić V, Pillet C A and Shirikyan A 2017 Journal of Statistical Physics 166 926–1015
  * [32] Lam K D N T and Kurchan J 2014 Journal of Statistical Physics 156 619–646
  * [33] Goldfriend T 2023 Journal of Statistical Physics 190 70
  * [34] Lopez-Piqueres J and Vasseur R 2023 Physical Review Letters 130 247101
  * [35] Basile G, Bernardin C, Jara M, Komorowski T and Olla S 2016 Thermal conductivity in harmonic lattices with random collisions Thermal transport in low dimensions (Springer) pp 215–237
  * [36] Lepri S, Mejia-Monasterio C and Politi A 2010 J. Phys. A: Math. Theor. 43 065002
  * [37] Lepri S 2023 Journal of Statistical Physics 190 16
  * [38] Bouchaud J P and Georges A 1990 Physics reports 195 127–293
  * [39] Hastings M 2003 Physical review letters 90 148702
  * [40] Mithun T, Kati Y, Danieli C and Flach S 2018 Physical review letters 120 184101
  * [41] Danieli C, Mithun T, Kati Y, Campbell D K and Flach S 2019 Phys. Rev. E 100(3) 032217 URL https://link.aps.org/doi/10.1103/PhysRevE.100.032217
  * [42] Basile G, Bernardin C and Olla S 2006 Phys. Rev. Lett. 96 204303
  * [43] Delfini L, Lepri S, Livi R, Mejia-Monasterio C and Politi A 2010 J. Phys. A: Math. Theor. 43 145001
  * [44] Iacobucci A, Legoll F, Olla S and Stoltz G 2010 J. Stat. Phys. 140 336–348 ISSN 1572-9613 URL https://doi.org/10.1007/s10955-010-9996-6
  * [45] Lepri S, Livi R and Politi A 2020 Physical Review Letters 125 040604
  * [46] Schnakenberg J 1976 Reviews of Modern physics 48 571
  * [47] Lecomte V, Appert-Rolland C and Van Wijland F 2007 Journal of statistical physics 127 51–106
  * [48] Garrahan J P, Jack R L, Lecomte V, Pitard E, van Duijvendijk K and van Wijland F 2007 Physical review letters 98 195702
  * [49] Garrahan J P, Jack R L, Lecomte V, Pitard E, van Duijvendijk K and van Wijland F 2009 Journal of Physics A: Mathematical and Theoretical 42 075007
  * [50] Andreucci F, Lepri S, Ruffo S and Trombettoni A 2023 Phys. Rev. E 108(2) 024115 URL https://link.aps.org/doi/10.1103/PhysRevE.108.024115
  * [51] Basile G, Olla S and Spohn H 2010 Archive for rational mechanics and analysis 195 171–203
  * [52] Tamaki S and Saito K 2020 Physical Review E 101 042118
  * [53] Matsuda H and Ishii K 1970 Prog. Theor. Phys. Suppl. 45 76
  * [54] Visscher W 1971 Prog. Theor. Phys. 46 729–& ISSN 0033-068X
  * [55] Bernardin C 2008 Journal of Statistical Physics 133 417–433
  * [56] Dhar A, Venkateshan K and Lebowitz J 2011 Physical Review E 83 021108
  * [57] Bernardin C, Huveneers F and Olla S 2019 Communications in Mathematical Physics 365 215–237
  * [58] Livi R, Pettini M, Ruffo S, Sparpaglione M and Vulpiani A 1985 Physical Review A 31 1039
  * [59] Muelken O, Heinzelmann S and Dolgushev M 2017 Journal of Statistical Physics 167 1233–1243
  * [60] Flach S 2010 Chemical Physics 375 548–556
  * [61] Aubry S and André G 1980 Ann. Israel Phys. Soc 3 18
  * [62] Harper P G 1955 Proceedings of the Physical Society. Section A 68 874 URL https://dx.doi.org/10.1088/0370-1298/68/10/304
  * [63] Domínguez-Castro G and Paredes R 2019 European Journal of Physics 40 045403
  * [64] De Bacco C, Guggiola A, Kühn R and Paga P 2016 Journal of Physics A: Mathematical and Theoretical 49 184003
  * [65] Coghi F, Morand J and Touchette H 2019 Physical Review E 99 022137