# Thermalization of isolated harmonic networks under conservative noise

Stefano Lepri1,∗

(1Consiglio Nazionale delle Ricerche, Istituto dei Sistemi Complessi
via Madonna del piano 10, I-50019 Sesto fiorentino, Italy
∗E-mail<EMAIL_ADDRESS>
)

###### Abstract

We study a scalar harmonic network with pair interactions and a binary
collision rule, exchanging the momenta of a randomly-chosen couple of sites.
We consider the case of the isolated network where the total energy is
conserved. In the first part, we recast the dynamics as a stochastic map in
normal modes (or action-angle) coordinates and provide a geometric
interpretation of it. We formulate the problem for generic networks but, for
completeness, also reconsider the translation-invariant lattices. In the
second part, we examine the kinetic limit and its range of validity. A general
form of the linear collision operator in terms of eigenstates of the network
is given. This defines an action network, whose connectivity gives information
on the out-of-equilibrium dynamics. We present a few examples (ordered and
disordered chains and elastic networks) where the topology of connections in
action spaces can be determined in a neat way. As an application, we consider
the classic problem of relaxation to equipartition from the point of view of
the dynamics of linear actions. We compare the results based on the spectrum
of the collision operator with numerical simulation, performed with a novel
scheme based on direct solution of the equations of motion in normal modes
coordinates.

##  1 Introduction

Nonequilibrium processes in many particle systems can be cathegorized in two
wide classes: transport induced by external forces (either mechanical and
thermodynamical) and relaxation to equilibrium. In the first class, we
encounter steady-state transport, or even time-dependent states like it
happens e.g. in surface growth, pattern formation, turbulent flows etc. The
second class concerns problems like thermalization and approach to
equipartition, coarsening, relaxation after a quench and so on. Both themes
continue to receive a large attention and the goal is to categorize them in
broad classes characterized by some form of universality.

The issue of thermalization has a long-standing history, starting with the
classic Fermi-Pasta-Ulam-Tsingou problem [1] until the most recent
developments, regarding the effect of integrability and quasi-integrability on
the statistical behavior of classical systems [2, 3, 4, 5, 6, 7].

One common trait in the study of such a problem is the difficulty of treating
genuine anharmonic systems which are usually tackled by molecular dynamics
simulations. This is the case of large ensembles of coupled nonlinear
oscillators driven out of equilibrium [8, 9, 10]. Even in the weakly-nonlinear
limit, perturbative methods like the KAM theory are notoriously of little use
in thermodynamically-large systems. In the kinetic regime, approaches based on
the phonon Boltzmann or wave-kinetic equation are well developed [11]. For
one-dimensional chains this turned out to be useful to understand anomalous
transport [12, 13, 14, 15] and (pre)thermalization caused by nonlinear
resonances [16, 17]

An alternative approach relies on stochastic modeling of the mesoscopic
dynamics, where the toolbox of stochastic processes can be fruitfully used.
Considerable insights is thereby obtained on the nature of the non-equilibrium
states and macroscopic transport laws [18, 19].

An intermediate viewpoint is to consider a sort of hybrid dynamics, where a
deterministic evolution is accompanied by stochastic interactions, possibly
preserving the basic conservation laws. In the simplest case, the
deterministic dynamics is linear and solvable, while the random part can be
seen as is a microcanonical Monte-Carlo rule, ergodizing the dynamics. This is
referred to as conservative noise dynamics: in its simplest versions it
entails random exchange of momenta between particle or a random reshuffling of
a subset of particles. Such schemes, allows sometimes for exact solutions or,
at least, very efficient simulations. For fluid systems, this is the strategy
of the Multi-Particle-Collision dynamics that can be fruitfully employed in
many diverse contexts [20, 21, 22]. For the class of oscillator chains, this
class of random dynamical systems allows for an exhaustive rigorous treatment
[23, 24, 25, 26]. Indeed, large-scale hydrodynamics equations can be
demonstrated and phonon Boltzmann equation can be derived, yielding relatively
simple linear collision operators [27, 28]. Moreover, most of nonequilibrium
steady-state properties can be computed exactly [29, 30, 31, 32] and were
shown to reproduce many features of deterministic nonlinear lattices [33]. The
effect of conservative noise on nonlinear oscillator chains has also been
considered [24, 34, 35, 36].

Another ingredient is the effect of quenched disorder, leading to subtle
interplay between Anderson localization and chaotic diffusion that affects
both thermalization and transport [37, 38, 39, 40, 41]. A further motivation
to tackle this case is that classical disordered chains with nonlinear
interactions exhibit a regime analogous to that seen in quantum many-body-
localized systems [42]. In addition to intrinsic heterogeneity of individual
units, disorder can originate also from topology and structure of connections,
as it occurs for the dynamics of elastic networks (see e.g. ref.[43] and
references therein). This type of structures have an intrinsic theoretical
interest, and may serve, for instance, as toy models of macromolecules (e.g.
proteins) in their native state [44, 45]. Another possible domain of
application is the study of nanoscale heat transfer trough networked
structures [46, 47], namely as toy models of devices composed of networks of
nanowires and nanotubes.

When dealing with weakly-interacting linear systems it is natural to refer to
the unperturbed harmonic modes and their energies or actions. For regular,
homogeneous, anharmonic lattices these are the usual Fourier modes. Their
perturbed dynamics determines the evolution of relevant observables both far
[48] than close [49, 50] to equilibrium. For disordered lattices one has
instead to consider Anderson modes that are typically localized in space.

Generally, when expressed in the familiar action-angle variables, a
nonintegrable perturbation defines a network of interaction among the
unperturbed actions [51, 52]. It is reasonable to argue that the connectivity
of such network will affect 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$ [51, 52]. In the nonlinear case
the coupling may involve three or more actions, resulting in complicated
hyper-graph structure. It is thus of interest to characterize the connectivity
of such action network in the simplest setting, and this is one of the aims of
the present paper.

In the present work, we present a general class of dynamical systems with
conservative-noise: an harmonic network with general pair interaction
described by a coupling matrix and equipped with a binary collision rule,
exchanging the momenta of a couple of randomly chosen pair of particles. The
motivations are twofold. First, we would like to understand how the
conservative-noise dynamics translates in the space of collective coordinates
(i.e. the eigenstates of the underlying harmonic network). In other words, how
do local interactions in space affects the normal modes. Second, we would like
to gain some insights on the connectivity and topology of the resulting action
networks, in the simplest setting. We also focus on the problem of relaxation
to equilibrium concentrating on the dynamics of linear actions, a subject that
has not been considered so far for conservative noise dynamics.

In the first part, we recast the stochastic dynamics as a stochastic map in
normal modes (or action-angle) coordinates and provide a geometric
interpretation of it. We will formulate the problem on generic scalar networks
but, for completeness, we will reconsider also the translation-invariant
lattices. As a byproduct, we derive an efficient and novel scheme to integrate
numerically the equations of motion in the collective coordinates.

In the second part, we examine the kinetic limit and its range of validity. We
compute the linear collision operator in the normal mode basis and discuss a
few examples where the topology of connections in action spaces can be
discussed in a neat way. As an application, we study the classic problem of
relaxation to equipartition a lá Fermi-Pasta-Ulam-Tsingou, starting from an
initial condition with a few actions excited. We compare the exact dynamics
with the kinetic results. A discussion and summary is given in the last
Section.

##  2 Harmonic network with conservative noise

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}\Phi_{ij}q_{i}q_{j},$
(1)

where the matrix $\Phi$ is semi-positive definite and symmetric so that its
eigenvalues are real and non-negative. The equation of motion of the isolated
network are

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

In order to have momentum conservation it must be $\sum_{i}\Phi_{ij}=0$. For
simplicity we will deal henceforth with the equal-mass case, and set
$m_{i}=1$.

We now introduce the following stochastic process [23]. Suppose we start with
the system at time $t$. At a later time $t+\tau$ there occur a ”collision”
event defined as follows. 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})$. Clearly, this move conserves
energy and the total momentum. Physically, it can be interpreted as a
perfectly elastic collision as it would occur for a infinite square-well pair
potential. In this general formulation, the model encompasses several
different setups. It includes the standard case of regular Euclidean lattices
when $\Phi$ is the familiar nearest-neighbor Laplacian matrix. The $W_{n,m}$
can be assigned to include some form of non-local interaction across the
lattice. The intervals between subsequent collision times $\tau$ are also
taken as random variables with some preassigned distribution with given,
finite, average $\langle\tau\rangle$. A natural choice would be, for instance,
the exponential distribution
$\exp(-\tau/\langle\tau\rangle)/\langle\tau\rangle$.

The normal modes’ coordinates of the network have eigenfrequencies
$\omega_{\nu}$ and are defined by the transformation

$Q_{\nu}=\sum_{l=1}^{N}q_{l}\,\chi_{l}^{\nu},\quad
P_{\nu}=\sum_{l=1}^{N}p_{l}\,\chi_{l}^{\nu}.$ (3)

The $\chi_{l}^{\nu}$, $\nu=1,\ldots,N$, are orthornormal and can be taken to
be real for the time being. Once the Hamiltonian is expressed in these new
canonical variables, the deterministic part of the equations of motion become

$\displaystyle\dot{Q}_{\nu}\;=\;\frac{\partial H}{\partial P_{\nu}}=P_{\nu},$
(4) $\displaystyle\dot{P}_{\nu}\;=\;-\frac{\partial H}{\partial
Q_{\nu}}=-\omega_{\nu}^{2}Q_{\nu}.$ (5)

Let us now determine the dynamics associated to the random part. At each
collision there is change in normal-mode momenta $\Delta P_{\nu}$ given by

$\Delta P_{\nu}=(p_{m}-p_{n})(\chi_{n}^{\nu}-\chi_{m}^{\nu})$

while the $Q$ are unchanged. Using the inverse of (3)

$p_{n}=\sum_{\nu}P_{\nu}\chi_{n}^{\nu},\quad
p_{m}-p_{n}=\sum_{\nu}P_{\nu}(\chi_{m}^{\nu}-\chi_{n}^{\nu})$

we obtain

$\Delta
P_{\nu}=(\chi_{n}^{\nu}-\chi_{m}^{\nu})\sum_{\mu}(\chi_{m}^{\mu}-\chi_{n}^{\mu})P_{\mu}.$

In the case in which $\Phi$ conserve the total momentum, momentum and stretch
are automatically conserved since $\Delta P_{0}=0$. So one can work with a
$N-1$ variables. Altogether, in column vector notation, the collision rule can
be written as

$\displaystyle P^{\prime}=P+\Delta P=(1-2VV^{T})P,$ (6) $\displaystyle
V_{\nu}^{(n,m)}\equiv\frac{\chi_{n}^{\nu}-\chi_{m}^{\nu}}{\sqrt{2}},$ (7)

where the prime denotes the value after the collision and $V^{T}$ the
transpose. Note that $V=V^{(n,m)}$ is a random vector as the indexes $n,m$ are
chosen at random with the prescribed rule (unless needed, we do not write
explicitly the dependence on the site indexes henceforth). It can be checked
that the vector $V$ is a unit vector, $V^{T}V=1$ Thus, equation (7) has an
interesting geometrical interpretation. Indeed, $(1-2VV^{T})$ it is a
Householder matrix that describes a reflection about an hyperplane containing
the origin. The hyperplane is orthogonal to the vector $V$. So the dynamics
can be seen as a sequence of reflections of the vector $P$ around a random
hyper-plane in phase space, a transformation that conserves the vector length
$|P|^{2}=P^{T}P$. The matrix properties are well known:

  * •

it is idempotent $(1-2VV^{T})^{2}=1$;

  * •

it is self adjoint $(1-2VV^{T})^{+}=(1-2VV^{T})$;

  * •

it has has one eigenvalue equal to $-1$ and the remaining $N-1$ equal to $1$.

  * •

also $\Delta P=-2VV^{T}P=-2(P^{T}V)V$ where $R=P^{T}V=|P|\cos\varphi$ is the
projection of $P$ on $V$ so

$P^{\prime}=P-2RV$

Since $|\Delta P|/|P|$ is of order one, the change in each component is of
order $1/\sqrt{N}$.

##  3 Equations for the collective variables

Let us now consider the dynamics in the usual action variables to understand
how the above transformations affects their dynamics. Introducing

$A=i(2\Omega)^{1/2}Q+(2\Omega)^{-1/2}P,$ (8)

with $\Omega_{\nu\nu^{\prime}}=\omega_{\nu}\delta_{\nu\nu^{\prime}}$ (note
that $A_{\nu}$ are $2N$ independent complex variables). The inverse formulae
(recall that $P,Q$ are assumed to be real):

$P=\frac{1}{2}(2\Omega)^{1/2}\left(A+A^{*}\right);\quad
Q=\frac{1}{2i}(2\Omega)^{-1/2}\left(A-A^{*}\right),$

with Hamiltonian transforming to $H=A^{+}\Omega A$. Substituting into (7) we
obtain the collision map:

$A^{\prime}=A-M(A+A^{*}),$ (9)

where the matrix $M$ is defined as

$M\equiv\Omega^{-1/2}VV^{T}\Omega^{1/2};\qquad
M_{\nu,\nu^{\prime}}=\sqrt{\frac{\omega_{\nu^{\prime}}}{\omega_{\nu}}}\,V_{\nu}V_{\nu^{\prime}}$

($V$ is real). The matrix $M$ is real, not symmetric since from the definition
$M^{T}=\Omega M\Omega^{-1}$ and is idempotent $M^{2}=1$ expressing the fact
that applying twice the same transformation does not change the state.

For the collective variables, the free evolution in the time interval
$(t,t+\tau)$ between subsequent collisions is diagonal

$A(t+\tau)=e^{i\Omega\tau}A(t),\qquad A^{*}(t+\tau)=e^{-i\Omega\tau}A^{*}(t).$
(10)

Combining the two processes (9) and (10), we obtain an exact map from $t^{+}$
to $(t+\tau)^{+}$

$A\to(1-M)e^{i\Omega\tau}A-Me^{-i\Omega\tau}A^{*}.$ (11)

The evolution thus amounts to a sequence of multiplication of random matrices
originating from the randomness of both the matrix $M$ and of the collision
intervals $\tau$. This formulation is particularly useful for implementing the
numerical solution as an event-driven dynamics, not requiring any approximate
integration schemes as in the case of ordinary or stochastic differential
equations (see Appendix A).

As a further step, we seek for the equations for the action-angle variables
$I,\theta$ defined by $A_{\nu}=\sqrt{I_{\nu}}e^{i\theta_{\nu}}$. From equation
(9) and recalling that $M$ is real, we obtain the collision map

$\displaystyle
I_{\nu}^{\prime}=I_{\nu}-4\sqrt{\frac{I_{\nu}}{\omega_{\nu}}}V_{\nu}Z\cos\theta_{\nu}+\frac{4V_{\nu}^{2}}{\omega_{\nu}}Z^{2}$
$\displaystyle=I_{\nu}\left[\sin^{2}\theta_{\nu}+\left(\cos\theta_{\nu}-\frac{2V_{\nu}}{\sqrt{\omega_{\nu}I_{\nu}}}Z\right)^{2}\right]$
$\displaystyle\sin\theta_{\nu}^{\prime}=\sqrt{\frac{I_{\nu}}{I^{\prime}_{\nu}}}\,\sin\theta_{\nu}$
(12) $\displaystyle
Z\equiv\sum_{\mu}V_{\mu}\sqrt{\omega_{\mu}I_{\mu}}\cos\theta_{\mu}.$

  * •

Those equations are exact: the first one guarantees that the energy is
conserved at each collision as $\sum_{\nu}\omega_{\nu}\Delta I_{\nu}=0$. Note
however that the total action $\sum_{\nu}I_{\nu}$ is not a constant of motion.

  * •

The transformation made the collision map nonlinear, and $Z$ is global
coupling among the eigenmodes: each collision entails a change of action and
angles of all the modes.

  * •

Finally, taking into account the free evolution, we can rewrite the dynamics
as a mapping from $t^{-}$ to $(t+\tau)^{-}$ (i.e. just before two subsequent
collisions)

$\displaystyle I\to I+\Delta I(I,\theta)$
$\displaystyle\theta\to\theta+\omega\tau+\Delta\theta(I,\theta),$

where $\Delta I=I^{\prime}-I$ and $\Delta\theta$ are defined by the (3) and
$\omega$ is the vector of the eigenfrequencies $\omega_{\nu}$. It is a a kind
of random mapping or discrete Langevin equation, (but note that it is not
symplectic).

Before proceeding further, we note that the calculation can be performed also
rewriting the equations (9) as a suitable stochastic equation. The details are
in Appendix B.

##  4 Kinetic equation for the actions

Let us now consider the kinetic limit, namely the case in which very many
collisions occur on the time scale $T$, $T\gg\langle\tau\rangle$. The limit in
which $N\to\infty$, $\langle\tau\rangle\to 0$ keeping

$\gamma=\frac{1}{N\langle\tau\rangle}$ (13)

finite, corresponds to choosing a finite collision probability per node of the
network. This is the choice mostly employed in the simulations [31, 34, 36].
In this case, on the time-scale $T\sim 1/\gamma$ the number of collisions
becomes macroscopically large. Qualitatively, the phases perform a random walk
and are randomized on a faster time scale with respect to the evolution of the
actions. We thus argue that the angles $\theta_{\nu}$ are quickly randomized
on this time scale and it is legitimate to average the first of (3) over a
uniform distribution of $\theta$, which is the expected invariant measure.
Denoting by $\overline{I_{\nu}}$ such averaged actions we obtain

$\overline{I_{\nu}^{\prime}}=(1-2V_{\nu}^{2})\overline{I_{\nu}}+2\frac{V_{\nu}^{2}}{\omega_{\nu}}\sum_{\nu^{\prime}}\omega_{\nu^{\prime}}V_{\nu^{\prime}}^{2}\overline{I_{\nu^{\prime}}}.$
(14)

In other words, this last equation represents the evolution of an ensemble of
trajectories with a given distribution of initial angles, subject to the same
sequence of collisions. Equivalently, in terms of the mode energies,
$E_{\nu}=\omega_{\nu}\overline{I}_{\nu}$:

$E_{\nu}^{\prime}=(1-2V_{\nu}^{2})E_{\nu}+2V_{\nu}^{2}\sum_{\nu^{\prime}}V_{\nu^{\prime}}^{2}E_{\nu^{\prime}}\equiv
E_{\nu}+\sum_{\nu^{\prime}}K_{\nu,\nu^{\prime}}E_{\nu^{\prime}}$ (15)

which, taking into account the normalization of $V$, shows immediately that
the sum of the $E_{\nu}$ is conserved also for the averaged equation, and we
have defined the angle-averaged collision matrix

$\displaystyle
K_{\nu,\nu^{\prime}}^{(n,m)}=-2V_{\nu}^{2}\,\delta_{\nu,\nu^{\prime}}+2V_{\nu}^{2}V_{\nu^{\prime}}^{2}.$
(16)

Each of the matrices $K$ has the following properties:

  * •

It is a symmetric matrix and it is random as it depends on the couple of
colliding particles $(n,m)$. It does not seem to belong to any known ensemble
of standard random matrices.

  * •

It does not depend on the eigenfrequencies and the distribution of collision
times but only on the eigemodes’ shape and on the collision probabilities
$W_{n,m}$

  * •

It is seen immediately that $K$ admits a zero eigenvalue with uniform
eigenvector $E_{\nu}=E_{eq}$ where $E_{eq}$ is the energy value of each mode
at equipartition.

  * •

$K$ is doubly-stochastic matrix, rows and columns sums to one.

We now would like to perform an average over realization of the collision
process. The formal solution from $t$ to $t+T$ is thus given in terms of a
product of (uncorrelated) random matrices

$E(t+T)=\left(\prod_{collisions(n,m)}(1+K^{(n,m)})\right)\,E(t)$ (17)

(with $1$ denoting here the identity matrix) and the product is extended to
all the collision occurring in the time interval $t,t+T$.On the basis of the
usual arguments based on Oseledec theorem, we thus expect that the associated
Lyapunov exponents exist [53].

If each application of the collision rule yields a small change in the vector
$E$ (for instance if all elements of $K$ are small), one can approximate the
product in (17) as

$\prod_{collisions(n,m)}(1+K^{(n,m)})\approx
1+\sum_{collisions(n,m)}K^{(n,m)}.$ (18)

This approximation is akin to the well-known weak-disorder expansion, a method
used to evaluate the product of random matrices for small disorder strengths
[53]. Its accuracy here will be afterwards checked numerically for the
specific examples considered below. If we accept its validity for the time
being we can, as a further step, replace the sum in (18) with the average
$\overline{K}$ of the matrix over the random process. In other words, one
considers the averaged operator coarse-grained over a time scale which is
longer than the typical collision time:

$\overline{K}=\sum_{(n,m)}W_{n,m}K^{(n,m)},\quad
R\equiv\frac{1}{\langle\tau\rangle}\overline{K}=\gamma N\overline{K}.$ (19)

The rate $\gamma$ just sets the overall time-scale of the kinetic process. The
(constant) matrix $\overline{K}$ inherits the properties of the $K^{(n,m)}$.
Altogether, replacing the differences in (17) with a time derivative, the
relaxation processes on such time scale is given by

$\dot{E}_{\nu}=\sum_{\nu^{\prime}}\left[R_{\nu,{\nu^{\prime}}}E_{\nu^{\prime}}-R_{{\nu^{\prime}},\nu}E_{\nu}\right];\qquad
R_{\nu,{\nu^{\prime}}}\equiv\frac{2}{\langle\tau\rangle}\overline{V_{\nu}^{2}V_{\nu^{\prime}}^{2}}.$
(20)

This equation thus provides the seeked effective approximate evolution of mode
energies on a time scale where the network has undergone a large number of
collisions. Since energy is conserved, equation (20) can be seen as a Master
Equation in action space and the elements of $R$ can be interpreted as
transition probabilities from a state $E$ to $E^{\prime}$. They are of the
form of the well-known Fermi golden rule, involving the squared amplitudes of
the eigenmodes. For a generic coupling of the network, such that there is no
decoupling among different subsets they are all nonvanishing. Following the
usual arguments of Markov processes, the system is ergodic and approaches
equipartition according to the properties of the ”collision operator” defined
by (20). 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 $\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 for a quadratic Hamiltonian with such collision rules. All the others
$\mu_{\nu}$ must be strictly negative [54]. In particular, it is of interest
to look at the spectral gap $|\mu_{2}-\mu_{1}|=|\mu_{2}|$ that controls the
relaxation at long times in a finite network. The scaling with $N$ of
$\mu_{2}$ (also called the Fiedler eigenvalue in the context of diffusion on
graphs) can be recast in terms of the scaling of spectral density of the
action network. Let us denote with $\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$ (21)

one may estimate $\mu_{2}$ for a finite network from the condition
$\rho(\mu_{2})=1/N$, from which we obtain $|\mu_{2}|\sim N^{-d_{\mu}/2}$. Note
that, on general grounds, $d_{\mu}$ is a distinct quantity from the standard
spectral dimension $d_{s}$ of the original network (1), which is defined by
the same relation as (21) for the integrated density of eigenfrequencies
$\omega_{\nu}$ [55, 56].

Before proceeding further, let us comment on the validity of the weak-disorder
approximations (18). It rests on the fact that the change of the energy vector
is somehow small at each collision. This is plausibly true in the case of
extended modes, as it occur in the case of translation-invariant or weakly
disordered lattices that will be treated below: each $\chi$ is of order
$1/\sqrt{N}$ for normalization. From the very definition of the $V$, equation
(7), we see that the matrix elements of the relaxation matrix should be small
for large $N$. In other words, each local collision (which yield a change
$O(1)$ of the particle momenta) yield a $O(1/N)$ change of the mode energies.
This expectation will be made precise in the standard case of the ordered
chain (see Section 5).

On the other hand, this may not be true in general, e.g if the network admits
localized modes. To understand this issue, let us consider the most trivial
example in which $\Phi$ is diagonal with frequencies
$\omega_{\nu}=\sqrt{\Phi}$ and trivial eigenvectors
$\chi_{m}^{\nu}=\delta_{m,\nu}$ (assuming that they are labeled according to
their spatial location). The matrix $K$ can be worked out: it is all zeros
except for two diagonal elements and two elements in positions $(n,m),(m,n)$,

$K_{\mu,\nu}^{(n,m)}=\begin{cases*}-\frac{1}{2}&if $(\mu,\nu)=(m,m),(n,n)$\\\
\phantom{-}\frac{1}{2}&if $(\mu,\nu)=(m,n),(n,m)$\\\ 0&otherwise.\end{cases*}$
(22)

Thus, the matrix is very sparse with $O(1)$ entries, very different from the
case of extended plane waves. However the change involves only one mode, thus
resulting in a relatively small change of the vector $E$ even in this case.
111Actually, this extreme case is basically the well-known KMP model whereby
harmonic oscillators exchange energy stochastically [18]. But, it may be also
related to deterministic models like the linear chain of oscillators with
hard-core collisions, termed the ding-dong model [57] in its disordered
version [58]. It would correspond to the ”stochastic approximation” of the
dynamics whereby the deterministic sequence of collisions is replaced by a
random one.

For a general $\Phi$ the eigenvectors are not known analytically but can be
easily computed numerically by exact diagonalization. Then, one can compute
the matrix form of the collision operator by averaging over the chosen
distribution $W_{n,m}$. The structure of the matrix will yield information on
the structure of action network. An analytically treatable case is the one of
ordinary lattices that we will discuss below.

##  5 Translation-invariant lattice

We now consider the case of a translation-invariant lattice, focusing on a
one-dimensional chain with periodic boundary conditions. In this case of
$\Phi$ is a circulant matrix: this is a well-studied case both for short
(nearest-neighbour) [23] than for long-range interactions [59]. In both cases
the eigenvectors are the familiar lattice Fourier modes

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

where $k_{\nu}=\frac{2\pi\nu}{N}$ are the wavenumbers and

$Q_{\nu}={1\over\sqrt{N}}\sum_{l=1}^{N}q_{l}\,e^{ik_{\nu}l},\qquad
Q_{-\nu}=Q_{\nu}^{*},\qquad\nu=-{N\over 2}+1,\ldots,{N\over 2}.$ (23)

(for simplicity of notation we use matrix indexes in the same range). A
difference with the case above is that the vectors $V$ are complex. This
requires some minor modifications when expressing the equation of motion in
phononic variables. As explained in Appendix A, the collision map (11) now
reads

$A\to(1-M)e^{i\Omega\tau}A-Me^{-i\Omega\tau}\tilde{A}^{*},$ (24)

where we introduced the shorthand notation $(\tilde{A}^{*})_{\nu}=A_{-\nu}$.
This makes transparent how the collision involves scattering of phonons with
opposite propagation directions. The matrix $M$ is formally as before but is
now complex-valued

$M\equiv\Omega^{-1/2}VV^{+}\Omega^{1/2};\qquad
M_{\nu,\nu^{\prime}}=\sqrt{\frac{\omega_{\nu^{\prime}}}{\omega_{\nu}}}\,V_{\nu}V_{\nu^{\prime}}^{*}$
(25)

and it is non-Hermitean and idempotent.

The collision map in action-angle variables is:

$\displaystyle
I_{\nu}^{\prime}=I_{\nu}-4\sqrt{\frac{I_{\nu}}{\omega_{\nu}}}Re(V_{\nu}e^{-i\theta_{\nu}})Z+\frac{4|V_{\nu}|^{2}}{\omega_{\nu}}Z^{2}=$
$\displaystyle
I_{\nu}\left|1-\frac{2V_{\nu}e^{-i\theta_{\nu}}}{\sqrt{I_{\nu}\omega_{\nu}}}Z\right|^{2}$
(26)
$\displaystyle\sin\theta_{\nu}^{\prime}=\sqrt{\frac{I_{\nu}}{I^{\prime}_{\nu}}}\,\sin\theta_{\nu}-2\frac{Im(V_{\nu})}{\sqrt{\omega_{\nu}I^{\prime}_{\nu}}}Z$
(27) $\displaystyle
Z\equiv\frac{1}{2}\sum_{\mu}\sqrt{\omega_{\mu}I_{\mu}}(V_{\mu}e^{-i\theta_{\nu}}+V_{\mu}^{*}e^{i\theta_{\nu}})$

(note that $Z$ is real).

Accordingly, in the kinetic limit we can proceed with the averaging as above.
The resulting kinetic equations are equal to (14,15) and (16) with the
$V_{\nu}^{2}$ replaced by

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

With respect to the general network, for linear chains there are some
important simplifications.

  1. 1.

In the standard case of nearest-neighbor collisions $m=n+1$ so that the above
quantity is independent of $n,m$

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

and the matrix $K$ is constant. Note that this remains true even for any
fixed-distance collision rule $l=m-n$. The collision operator is determined
straightforwardly from definitions (13) and (19)

$R_{\nu,\nu^{\prime}}=\gamma
N{K}_{\nu,\nu^{\prime}}=-4\gamma\sin^{2}\frac{k_{\nu}}{2}\delta_{\nu,\nu^{\prime}}+\frac{8\gamma}{N}\sin^{2}\frac{k_{\nu}}{2}\sin^{2}\frac{k_{\nu^{\prime}}}{2}$
(29)

Since the off-diagonal terms are small for large $N$, the eigenvalues of the
matrix are well approximated by the diagonal elements giving the approximation

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

up to higher-order corrections in $1/N$ that may be computed perturbatively.
Also, the $\nu$th eigenvectors is localized on $\pm\nu$ with all other
components being small of order $1/N$. Since the spectrum has a vanishing gap
for large $N$, the relaxation rate of a generic non-equilibrium initial
condition is expected to occur on of the slowest time scale
$\mu_{2}\approx\gamma(\pi/N)^{2}$ at long times.

  2. 2.

A more general model would consist to consider a collision probability of the
form

$W_{n,m}=\,w_{n-m}$ (31)

meaning that we choose an random particle $n$ with equal probability and a
neighbor at distance $l=n-m$ with probability $w_{l}$ (with periodic boundary
condition assumed), the standard case above being $w_{l}=\delta_{l,1}$. In
this case, the matrix $K$ is no longer constant. The averaged collision
operator is expressed as

$R_{\nu,\nu^{\prime}}=\gamma
N{\overline{K}}_{\nu,\nu^{\prime}}=\gamma\sum_{l}w_{l}\left[-4\sin^{2}\frac{k_{\nu}l}{2}+\frac{8}{N}\sin^{2}\frac{k_{\nu}l}{2}\sin^{2}\frac{k_{\nu^{\prime}}l}{2}\right]$
(32)

which, as above, suggests the following approximation in terms of the diagonal
elements

$\mu_{\nu}\approx-{4\gamma}\sum_{l}w_{l}.\sin^{2}\frac{k_{\nu}l}{2}.$ (33)

We will discuss an example later on.

Before concluding this Section, we note that the choice of the $w_{l}$ is
crucial to assess the dependence of the $\mu_{\nu}$ on the wavenumber
$k_{\nu}$, as seen already at the level of the approximation (33). In
particular, introducing correlations and/or long-ranged rules may yield non-
standard relaxation and transport depending on the range and correlation
strengths. We will defer investigation of those issues to future works.

##  6 Examples

In this section we illustrate the above in some specific examples and check
for the validity of the approximations done.

###  6.1 Disordered chain

For a first test, let us consider the case of 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]$
(34)

with $\sigma_{i}$ being i.i.d. variables with uniform distributions in
$[0,w]$, and $w$ gauges the disorder strength. Periodic boundary conditions
are assumed. As it is well known, the eigenstates are the exponentially-
localized Anderson modes, whose localization length decreases with $w$ [60,
61]. We consider the standard case of nearest-neighbor collisions. Several
variants in this type of model have been considered earlier [62, 63, 64].

In figure 1 we illustrate the kinetic regimes for two cases corresponding to
relatively weak and strong disorder. We first of all compute the relaxation
spectrum as the Lyapunov exponents of the product or random matrices (17) via
the usual QR algorithm [53]. Then, we compare the result with the eigenvalues
of the averaged collision matrix as obtained from the weak-disorder expansion
(18). In the case of weak disorder, figure 1(a) , the two methods give the
same spectra, thus confirming the accuracy of the approximation. Moreover, the
eigenvalues are also very close to the diagonal elements which, in turn, are
almost indistinguishable from the one of the ordered case (see below). This
means that the collision operator is almost diagonal with very small off-
diagonal elements (see figure 1b).

The situation is different in the strong disorder case. Here, only first 25%
of Lyapunov exponents coincide with the eigenvalues coincide and they are both
very different from the diagonal elements (see figure 1c). This is presumably
due to the fact that the dynamics is comparatively less ”mixing” and stronger
deviations from the average occur on shorter times (those corresponding to
more negative Lyapunov exponents). Accordingly, the structure of the collision
matrix is comparatively less homogeneous, with relatively larger off-diagonal
elements (figure 1d).

Figure 1: Disordered chain (34), $N=50$, $\gamma=0.1$ with nearest neighbour
collisions and disorder strengths $w=0.1$ (a,b) $w=10$ (c,d). Panels (a,c)
report the Lyapunov exponents of equation (20): Circles are the numerical
values computed via the QR algorithm, crosses are from the weak-disorder
expansion (eigenvalues of the average matrix, as in equation (18)) The solid
lines are the diagonal elements of the matrix itself. The dashed in panel (a)
line are the diagonal elements in the case of no disorder $w=0$, equation
(30). Panels (b,d): color plot of the averaged collision matrix (19).

A further insight is obtained by looking at the connectivity in action space.
As said above, the transition rates are all-nonvanishing, but can be very
small due to the exponential localization of the eigenvectors. For the sake of
visualization, we draw a representation of action network as an adjacency
matrix as follows. We consider two modes $\nu,\nu^{\prime}$ to be connected if
$|R_{\nu,\nu^{\prime}}|$ is larger than some preassigned threshold and not-
connected otherwise. As seen in figure 2 there is a clear qualitative
difference between the two cases: increasing disorder the connectivity
decreases and passes from a complete-graph type to a more sparse structure.
Accordingly, the degree distributions (the histogram of the number of links of
each node) is distinctly different. For weak disorder, each action is
connected to $O(N)$ ones, while the average degree is finite in the strong
disorder case.

Figure 2: Disordered chain (34), $N=50$, $\gamma=0.1$ with nearest neighbour
collisions and disorder strengths $w=0.1$ (a,b) $w=10$ (c,d). (a,c) Graphical
representation of the action network: it corresponds to an adjacency matrix
where two nodes $\nu,\nu^{\prime}$ are consider to be connected if
$|R_{\nu,\nu^{\prime}}|>10^{-4}$. (b,d): degree distributions.

###  6.2 Mean-field chain

To test also the case of translation-invariant lattice we consider, for the
sake of an example, a standard ordered chain with nearest-neighbour with an
additional mean-field interaction:

$H=\sum_{i=1}^{N}\left[{p_{i}\over
2}^{2}+\frac{1}{2}(q_{i+1}-q_{i})^{2}\right]+\frac{k}{2N}\sum_{ij}(q_{i}-q_{j})^{2}.$
(35)

This type of chains has been discussed as a model for Bose-Einstein condensate
[65] and studied recently in a nonequilibrium setup [66, 67]. The
eigenfrequencies are

$\omega_{\nu}^{2}\;=\;2k+4\sin^{2}\left({\pi\nu\over N}\right)\quad.$ (36)

for $\nu\neq 0$ and $\omega_{0}=0$. Apart from this zero mode, the dispersion
is formally the same as the one of a chain with harmonic pinning.

As a result of the calculation in Section 5, the kinetic equations are
independent on $k$ (this will be checked later on). For the sake of testing
the validity of the kinetic approximation, we here will concentrate on case
$k=0$, which is well-studied and sometimes referred to as the acoustic chain,
since the dispersion of waves is linear for small wavenumbers. We will profit
also to test the effect of different collision rules of type (31).

We first start with the relaxation rates obtained by the kinetic
approximation. As above, we computed the spectrum numerically using the QR
decomposition method [53] and compare it with the weak-disorder expansion and
diagonal approximation. In figure 3 we report the relaxation rates computed
with the different methods in two instances of (31) , namely the standard
nearest-neighbor exchange and the one with a dichotomic choice between both
the nearest and next-to-nearest neighbor, with probabilities $s,1-s$
respectively,

$w_{l}=s\delta_{l,1}+(1-s)\delta_{l,2}.$ (37)

In the first case, the matrix $R$ is given by (29) while in the second we need
to compute the Lyapunov spectra of the product of two constant matrices,
randomly chosen with probability $s$. This correspond to the case of binary
disorder in the random matrix language [53].

As illustrated in figure 3 the three methods agree fairly well for both
collision rules. We examined also other values of $k$ (data not reported)
obtaining identical results as predicted. It should also be noticed that the
diagonal approximation does reproduce quite accurately the spectrum already
for $N=50$.

Figure 3: Lyapunov exponents of equation (20) for the case of the translation-
invariant acoustic chain (35) with $k=0$, $N=50$, $\gamma=0.1$ and two types
of for collision rule: among (a) nearest-neighbor and (b) nearest-neighbor and
next-to-nearest neighbor sites with probability $s=1/2$ each, equation (31.
Circles are the numerical values of the Lyapunov exponents computed via the QR
algorithm, crosses are from the weak-disorder expansion namely the eigenvalues
$\mu_{\nu}$ of the averaged matrix, as in equation (18). Solid lines are the
diagonal approximation given by (30) and (33). To take into account the double
degeneracy, only half of the exponents is plotted. Finally, the triangles are
the relaxation rates measured from a simulation of a chain with the same
parameters and $N=128$ as described in Section 7.

###  6.3 Elastic network

As an example of a more complicate topology, we consider a more general
elastic network [43, 68]

$H=\sum_{i=1}^{N}\left[{p_{i}\over
2}^{2}+\frac{1}{2}(q_{i+1}-q_{i})^{2}\right]+\frac{1}{2}\sum_{ij}C_{i,j}(q_{i}-q_{j})^{2},$
(38)

where $C_{i,j}=1$ if the nodes are connected and zero otherwise. As a case
study, we choose the well-known Newman-Watts-Strogatz network, generated
starting from a ring where each oscillator is coupled to its neighbor [69].
Each oscillator is then connected by to with randomly-chosen existing one with
probability $p$. Increasing the probability $p$ increases the level of
disorder and we may expect a change the normal mode structure from extended
plane-waves to localized. The spectral properties have been studied in great
detail in [70]. The main feature is the existence of a pseudo-gap with very
small eigenfrequencies below a certain value.

In figure 4 we we compare for two cases corresponding to a single realization
of $C_{i,j}$ and for relatively weak and strong disorder. We limit ourselves
to compare the eigenvalues of the averaged collision matrix from the weak-
disorder expansion (18) with the diagonal elements. In both examples, the two
are pretty close, especially the smallest ones. This means that the collision
operator is almost diagonal with very small off-diagonal elements, similar to
the case of the strongly disordered chain.

Figure 4: Watts-Newman-Strogatz elastic network (38) $N=50$, $\gamma=0.1$ with
nearest-neighbour collisions and $p=0.1$ (a,b) $p=0.8$ (c,d). Panels (a,c)
report the Lyapunov exponents of equation (20): crosses are from the weak-
disorder expansion (eigenvalues of the average matrix, as in equation (18))
The solid lines are the diagonal elements of the matrix itself. Panels (b,d):
color plot of the averaged collision matrix (19).

In figure 5 we compare the degree distribution in the real network with the
action network (as defined above) for two values of $p$. As it is well known,
the degree distribution of the is peaked around the average value [69]. On the
contrary, the action network, constructed as described above, is fully
connected, with practically all-to-all couplings.

Figure 5: Watts-Newman-Strogatz elastic network (34), $N=50$, $\gamma=0.1$
with nearest-neighbour collisions and $p=0.1$ (a,b) $p=0.8$ (c,d). Comparison
of degree distributions in the real network (a,b) and in the the action
network (c,d): it corresponds to an adjacency matrix where two nodes
$\nu,\nu^{\prime}$ are considered to be connected if
$|R_{\nu,\nu^{\prime}}|>10^{-4}$. .

##  7 Thermalization

We perform a thermalization simulation in the spirit of the Fermi-Pasta-Ulam-
Tsingou famous numerical experiments. We initialize the chain by exciting a
packet of $N_{0}$ normal modes to have a prescribed energy, namely we set

$A_{\nu}(0)=\sqrt{\frac{E_{0}}{2N_{0}\omega_{\nu}}}e^{i\alpha_{\nu}}$ (39)

for $0<\nu\leq N_{0}$ ($|\nu|\leq N_{0}/2,\nu\neq 0$ in the case of
translation-invariant chain) where $\alpha_{\nu}$ being i.i.d.random phases
uniformly distributed in $[0,2\pi]$. The remaining modes are initialized as
$A_{\nu}(0)=0$. The dynamics is evolved iterating the dynamics as prescribed,
with random Poissonian collision times whose average is given by (13).

The evolution of the mode energies $E_{\nu}$ towards equipartition is
monitored and compared with the predictions. To improve the statistical
accuracy a time-average over a fixed number of collision (typically $10^{3}$)
is performed, along with an average over different initial conditions, i.e.
different realization of $\alpha_{\nu}$ in (39). At long enough times, the
average mode energies converge to the equipartition value as they should.

###  7.1 Disordered chain

In figure 6 we report the simulation data of (34) for nearest-neighbor
collisions and two values of the disorder parameter $w$, starting with all the
energy fed into a few modes. In the upper panels we show the approach to
equipartition of some of the excited modes. After an initial transient the
decay sets to the exponential behavior with a rate that matches fairly well
the value of the eigenvalue $\mu_{2}$ (dashed lines). The mode evolution and
repartition is illustrated in the lower panels. For weaker disorder, there is
a steady transfer from the initial modes to the others that leads to a uniform
depletion of their energies. On the other hand, for stronger disorder the
transfer is in a sense, more irregular and involves modes with $\nu$ values
different from $\nu_{0}$. This is qualitative agreement with what we would
expect in view of the different structures of the collision matrices as given
in figure 1 and 2 and the different connectivity in action space.

Figure 6: Disordered chain (34), for $w=2$ (left panels) $w=10$ (right panels)
$\gamma=0.1$ with nearest-neighbor collisions. initial wavepacket lowest-
frequency $N_{0}=5$ modes starting from $\nu=0$ and $465$ respectively. Upper
panels (a,d): solid lines are $E_{\nu}(t)-E_{eq}$, dashed line corresponds to
$0.3\exp(-|\mu_{2}|t)$. Middle (b,f): snapshots of the mode distributions at
different times during the thermalization process. Lower panels (c,g):
evolution of distributions $E_{\nu}(t)$ in mode space, equipartition
corresponds to uniform color (up to fluctuations). Data are averaged over 2000
consecutive collisions and an ensemble of about 100 initial conditions, with
different choices of $\alpha_{\nu}$ in (39). Notice the difference in time
scales between left and right panels.

Another observation from figure 6 is that the time scales involved are much
smaller in the case of weaker disorder. To understand this issue, In figure 7a
we compare the integrated spectral density $\rho(\mu)$ of (34) in the case
with no disorder $w=0$ and $w=0.1$, and for different sizes. As expected, in
the ordered case, $\rho\sim\sqrt{|\mu|}$, i.e $d_{\mu}=1$. Remarkably,
introducing a weak randomness $\rho$ does approach zero for a finite $\mu$.,
i.e the spectral gap is finite, meaning that $d_{\mu}=\infty$ and that the
relaxation rate is finite and $N$-independent. This result should be however
be taken with some care in the light of 7b that reports the integrated
spectral density for stronger disorder. In the first case, the crossover to
the $\sqrt{\mu}$ occurs at very small values of $\mu$. It is thus possible
that such a crossover will occur to smaller and smaller $\mu$ values upon
decreasing the disorder, making $d_{\mu}$ converge to one in the thermodynamic
limit.

Figure 7: The integrated spectral densities $\rho(\mu)$ (i.e. the fraction of
eigenvalues less than $\mu$) of the collision operator for model (34). Plots
are for different disorder strengths and three lattice sizes
$N=512,1024,2048$; $\gamma=0.1$ case with nearest-neighbor collisions. (a)
Weak disorder: Comparison between with $w=0$ (ordered case) and $w=0.1$ , and
sizes $N=512,1024,2048$; (b) Intermediate and strong disorder, $w=4.2,40$. The
dashed line correspond to the scaling $\sqrt{|\mu|}$ expected for the ordered
lattice.

To further characterize the effect of increasing disorder, in figure 8a we
report the spectral gap for different lattice lengths as a function of $w$
along with the average inverse participation ratio (figure 8b), defined as

$IPR=\left[\frac{1}{N}\sum_{\nu}\sum_{n}\left|\chi_{n}^{\nu}\right|^{4}\right]^{-1}$
(40)

which is $O(1)$ for localized modes and $O(N)$ for extended ones. It thus
gives a rough measure of the average localization length.

From the numerical data, it is seen that there exist a value $w_{*}(N)$ such
that for $w<w_{*}$ the spectral gap is almost $N$-independent. This confirms
that the relaxation is faster for weak enough disorder. On the other hand, for
$w>w_{*}$ the spectral gap vanishes as $N^{-2}$, figure 8c. A possible
argument to estimate the crossover value $w_{*}$ is to identify it with the
point at which the localization length becomes of the order the lattice size.
Above such value some of the coefficients $V_{\nu}$ start to become
exponentially small and the collision matrix more and more sparse accordingly.
This estimate roughly correspond to the location of the maxima of $|\mu_{2}|$
in figure 8a. Indeed, using the average IPR as a measure of the localization
length, we see that $w_{*}$ should decrease with $N$, in qualitative agreement
with the leftmost shifting of the maxima in figure 8a.

The $1/N^{2}$ scaling of the spectral gap suggests that the thermalization is
ruled by some form of diffusive process in action space. This can be
understood qualitatively as follows. In the strong disorder case,
$w\to\infty$, one can to a first approximation, neglect the coupling in (34)
and consider the eigenmodes as localized on a single site. We are thus in
close to the case of independent oscillators and the matrix $K$ can be
approximated, up to exponentially small terms, by (22). So, considering the
case of nearest-neighbor collisions, we argue that, up to a re-ordering of the
eigenmodes indexes $\nu$ by their spatial positions, (20) is approximated as

$\dot{E}_{\nu}=\frac{\gamma}{2}(E_{\nu+1}+E_{\nu-1}-2E_{\nu}),$ (41)

which immediately shows that the relaxation is a Brownian diffusion process in
action space. Thus thermalization starting from a bunch of nearby modes would
entail an initial growth of the number of excited modes proportional to
$\sqrt{\gamma t}$. For a finite chain the longest time scale is given by the
smallest eigenvalue of (41) which is of order $\gamma(2\pi/N)^{2}/2$. This is
consistent with the scaling of $\mu_{2}$ as given in figure 8.

Figure 8: Disordered chain (34), $\gamma=0.1$ with nearest-neighbor
collisions. (a) The absolute value of the second eigenvalue $\mu_{2}$ of the
averaged collision matrix (19) and (b) the average inverse participation ratio
(40) as a function of disorder strength for different chain lengths; (c)
reports the dependence of the eigenvalue on the lattice size $N$ for different
$w$. The dashed line correspond to the smallest non-zero eigenvalue of (41),
$\gamma(2\pi/N)^{2}/2$.

Before passing to the next example, we mention the work [71] where a kinetic
theory in a weakly-disordered nonlinear Schrödinger chain in the regime of
homogeneous chaos. The regime there considered, is the one of weak enough
interaction for the normal modes of the linear problem to remain well resolved
but strong enough for the dynamics to be chaotic for almost all modes. In that
case, the kinetic equations are nonlinear, leading to nonlinear diffusion
[71].

###  7.2 Mean-field chain

Let us now consider again model (35). In figure 9 we report the data for the
case for the acoustic chain with $k=0$ and nearest-neighbor collisions. The
relaxation of the four lowest modes is plotted in figure 9a,b along with the
snapshots of the evolution of mode energy distributions in wavenumber space,
9c,d. The simulation data are in excellent agreement with the relaxation rates
computed within the kinetic approach.

From the plots it is seen that thermalization occurs by gradual transfer of
the energy from the excited modes towards all the others. Indeed, the energy
in the whole background increases steadily and equipartition is reached with
an overall rate $\mu_{2}\sim 1/N^{2}$ as predicted.

Figure 9: Simulations of relaxation to equipartition for the chain (35) with
nearest-neighbor collisions, $k=0$, $N=512$, $\gamma=0.1$: initial wavepacket
consisting of lowest-frequency $N_{0}=20$ modes. (a) Evolution of
$E_{\nu}-E_{eq}$ the four lowest modes Dashed lines are exponential
$95\exp(-|\mu_{\nu}|t)$ with rates $\mu_{\nu}$ given by the approximation
(30); plot of dynamics in mode space (c) and snapshots (d) of the mode
distribution during the thermalization. (b) Evolution of the first mode
$E_{1}-E_{eq}$ as a function of $\mu_{\nu}t$ for $k=0$ , $k=1$ and pure mean-
field coupling, dashed line is $95\exp(-|\mu_{2}|t)$. Data are averaged over
10000 consecutive collisions and over an ensemble of about 100 initial
conditions, with different choices of $\alpha_{\nu}$ in (39).

We also tested the case with collision probability given by (37) obtaining
similar results. In figure 3 we plot as red triangles the relaxation rates
obtained by exponential fitting of $E_{\nu}(t)$: the agreement is again very
good.

As said, the kinetic equation in this translation-invariant should be
independent of $k$, if the the collision probability $W_{n,m}$ is the same.
Thus, we simulated the model for different $k$ maintaining the collisions only
among nearest neighbors In figure 9b we confirm this expectation comparing the
relaxation of the first Fourier mode for three different cases.

###  7.3 Elastic network

For the Newman-Watts-Strogatz network (38) we limit ourselves to examine the
spectral gap of the collision operator for different disorder strengths and
different number of sites. We just give a first account, leaving a more
complete study of such class of networks to future work.

For each $N$ we consider a fixed random realization of the coupling matrix
$C_{i,j}$ and change the probability $p$. A few different realizations of the
couplings were examined with qualitatively similar results indicating that
sample-to-sample fluctuations may not be very relevant. In figure 10a we
compare the integrated spectral density $\rho(\mu)$ of (34) for different
sizes and $p$. The data are compatible with a finite spectral gap for $p>0$,
and accordingly $d_{\mu}=\infty$, yielding a finite and $N$-independent
relaxation rate.

Figure 10: The integrated spectral densities $\rho(\mu)$ of the collision
operator (i.e. the fraction of eigenvalues less than $\mu$) for the Newman-
Watts-Strogatz lattices for different $p$ and three lattice sizes
$N=512,1024,2048$; $\gamma=0.1$ case with nearest-neighbor collisions. (a)
Comparison between $p=0$ (ordered case) and $p=0.11$ , and sizes
$N=400,800,1600$; The dashed line correspond to the scaling $\sqrt{|\mu|}$
expected for the ordered lattice. (b) Intermediate and strong disorder,
$p=0.78,0.89$.

In figure (11)a,b we plot the spectral gap as a function of $p$ along with the
IPR. The situation is similar to the disordered chain in the weak-disorder
regime see again figure 8. The spectral gap is finite for any finite $p$ and
closes for $p\to 0$. We do not observe a crossover to $1/N^{2}$ or similar
scaling as in the case of the disordered chain.

Figure 11: Newman-Watts-Strogatz network (38), $\gamma=0.1$ with nearest-
neighbor collisions. (a) The absolute value of the second eigenvalue $\mu_{2}$
of the averaged collision matrix (19) and (b) the average inverse
participation ratio (40) as a function of disorder strength for different
chain lengths; (c) reports the dependence of the eigenvalue on the lattice
size $N$ for different $w$.

##  8 Spectral entropy

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

$S(t)=-\sum_{\nu}E_{\nu}\log E_{\nu};\qquad\sum_{\nu}E_{\nu}=1,$ (42)

where we have fixed the total energy to one, without loosing generality. For
large $t$, $S$ approaches its equipartition value $\log N$. In general $S$
will depend on the initial conditions.

As we have shown above, the energy transfer processes in the translation-
invariant models occurs through a gradual, global, redistribution of the
initial energy towards all the other modes. In the intermediate times, and
assuming all the energy initially in one single mode
$E_{\nu_{0}}(0)=\delta_{\nu,{\nu_{0}}}$ for simplicity, one could thus perform
a kind of ”mean-field” approximation [73].  By this we mean that the 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)$ (43)

where $E_{\nu_{0}}(t)\approx(1-1/N)\exp(-\mu_{\nu_{0}}t)+1/N$. In figure 12 we
plot the simulated time evolution of $S$ for the harmonic and disordered
chains. The above formula accounts very well for the data.

The situation is instead pretty different in the case of the strongly
disordered chain (see the full circles in figure 12b). Here, the action
network is different yielding a diffusive process. Thus the entropy grows
logarithmically in the intermediate time range. This is readily understood
since, as said above, the number of excited modes grows as a $\sqrt{t}$. This
fits with the general fact, that there exist a close relation between the
growth of $S$ and the spectrum of the master equation for diffusion on
networks. In this context, it is known that if the spectral density of the
Laplacian operator of the graph shows scaling with a finite spectral dimension
$d_{s}$, the Shannon entropy grows as $\frac{d_{s}}{2}\log t$ with time before
reaching the steady-state (uniform) value  (see e.g. [73] and the references
therein) .

Figure 12: Spectral entropy (42), $\gamma=0.1$, $N=512$ for the harmonic and
disordered chains with nearest-neighbor collisions. (a) model (35) with
$k=0,1$ and (b) (34) with two values of the disorder. Initial conditions with
only the lowest frequency mode is excited. Dashed lines correspond to the
approximation (43); dot-dashed (red and violet) horizontal lines correspond to
the equipartition value $S=\log N$.

##  9 Discussion and perspectives

The conservative-noise dynamics for binary collision is a succession of random
reflections in action space, see equations (7) and (9). Averaging over random
phases yield a kinetic equations for the actions/energies of normal modes. The
mathematical advantage of such dynamics is that the collision operator is
linear and yields a linear master equation (20). But the main merit, for our
purposes, is that it was possible to construct explicitly the action network
in term of the eigenvectors and the collision rule.

Being a stochastic model there is no doubt that equipartition will be finally
obtained and that ergodicity is insured. However, one has the possibility to
study relaxation of an arbitrary finite network by simply computing the
spectral gap and its dependence on $N$. This amounts primarily to understand
how the spectral density of the action network depends on the underlying
connections. We also note that momentum conservation does not appear to play a
major role for thermalization. For instance, both the case of strongly
disordered case and the translation-invariant lattices display the same
$1/N^{2}$ scaling of relaxation rates, despite the fact that momentum is
conserved only in the second case.

Some comments are in order to conclude.

  * •

The spectral density $\rho(\mu)$ of the collision operator for small $\mu$ is
the main ingredient determining the the equipartition time in a finite
network, see (21). A natural question is whether there exist a relation
between $d_{\mu}$ and the spectral dimension in the case in which $\Phi$ of
the harmonic network in (35) is Laplacian matrix. In general, we surmise that
the two may not be necessarily related. This is understood in the case of the
translation-invariant chain. There, $d_{\mu}$ is controlled by the collision
rule only, as seen from (33) and can, in principle, be changed independently
of the spectral dimension the Laplacian. However, we have seen that for the
disordered chain that the relation is less straightforward and this issue may
deserve a further analysis.

  * •

To what extent are the present results applicable to describe a genuine
nonlinear network? One may guess that the stochastic dynamics it is an
idealization of a regime where chaos is well developed (i.e. the maximal
Lyapunov exponent is very large) and homogeneous in space. This would occur
for very large energy densities and/or strongly nonlinear networks. If so,
correlations should decay rapidly and replacing local chaos by random
collisions would be a reasonable assumption. It is also plausible that the
overall structure of the action network may be described in the same language
as here, although with higher-order connectivities (i.e. on a hyper-graph).
This is a problem that deserves further studies.

  * •

In the translation-invariant case the kinetic equation is basically the
spatially-homogeneous version of the well-known phonon Boltzmann equation (or
Boltzmann-Peierls equation) or wave-kinetic equation for phonon distribution
[11, 16]. The main difference is that the random dynamics does not conserve
the total number (i.e $\sum_{\nu}I_{\nu}$ in our notation). As it is known, a
more general derivation yields an advection term, proportional to the phonon
group velocity [11]. For the harmonic chain with conserving noise, it has been
derived in [27] with a more rigorous approach based on the Wigner
distribution. The difference in our case is that we consider the mode energies
and, more importantly, we dealt from the very start with the homogeneous case,
whereby the energy distribution is initially spatially uniform and it remains
so a later times. The inclusion of such terms is crucial to derive the correct
hydrodynamic scaling, leading to correct (fractional) heat equation [74] on
longer time-scales (see [25] and references therein for details).

  * •

Also, in the language of phonon kinetic theory, a nonlinear potential of order
$x^{p}$ induces interaction involving $p$ phonons. To the extent in which the
random collisions here described represent an hard-core potential, namely an
”infinite-order” $p\to\infty$ potential, it is thus understandable why the
phonon modes are globally coupled. An important difference is also that, in
the standard framework, energy and momentum conservation imposes constraints
on the allowed phonon processes: (for example conditions like
$\omega_{w}=\omega_{k}\pm\omega_{q}$ for three-phonons collisions etc. ).

  * •

In this work we focused on relaxation to equipartition of the isolated
network. As it is well-known there is a close relation between relaxation
close to equilibrium and transport. The relaxation rates can be used to
predict the correlation decay, and thus the transport coefficients. The
existence of a finite spectral gap of the collision operator is an indication
of fast correlation decay which, in turn, imply normal diffusive transport. In
the model we discussed this is relatively straightforward once the relaxation
matrix is known. It remains to be investigated how its spectral properties
depend on the underlying connections of the action network.

  * •

One may wonder if the same approach can be applied to other systems with more
conserved quantities. A relevant example would be the tight-binding electron
problem, modeled as a discrete Schroedinger equation. In this case the
collision rule should preserve both energy and the total wavefunction norm.
However, it has been shown that this requires a nonlinear rule [75]. Thus the
mathematical advantages following from the linearity are lost.

The class of models considered here will allow to extend the study of non-
equilibrium properties to the case of open networks interacting with external
reservoirs [76]. We plan to continue this program in the future.

Acknowledgements: I thank Raffaella Burioni, Stefano Iubini, Francesco Piazza
and Antonio Politi for useful correspondences and discussions during
elaboration of this work.

Funding: This work did not receive any funding.

Data availability: The datasets generated during and/or analysed during the
current study are available from the corresponding author on reasonable
request.

### Appendix A: Translation-invariant model

Using (8) we can proceed with the calculation as before except for the fact
that the normal coordinates $Q,P$ are complex:

$\displaystyle
A_{\nu}=i(2\omega_{\nu})^{1/2}Q_{\nu}+(2\omega_{\nu})^{-1/2}P_{\nu}$ (44)
$\displaystyle
A_{-\nu}^{*}=-i(2\omega_{\nu})^{1/2}Q_{\nu}+(2\omega_{\nu})^{-1/2}P_{\nu},$
(45)

and the inverse formulae now read:

$P=\frac{1}{2}(2\Omega)^{1/2}\left(A+\tilde{A^{*}}\right),$

where we introduced the shorthand notation $(\tilde{A})_{\nu}=A_{-\nu}$.
Substituting into (7):

$A^{\prime}=A-\Omega^{-1/2}VV^{+}\Omega^{1/2}(A+\tilde{A}^{*})\equiv(1-M)A-M\tilde{A}^{*},$
(46)

where the matrix $M$ is given by (25). Using equation (10) for the free
evolution and noticing that $\tilde{A}^{*}=e^{-i\Omega\tau}\tilde{A}^{*}(t)$
we obtain (24).

As stated in the text, the formulation of the equation of motion in the normal
modes coordinates is very convenient for the implementation of the numerical
solution. Taking into account the form of the matrix $M$ and that
$V_{\nu}^{*}=V_{-\nu}$ one can write the collision map as

$A^{\prime}=A-U(W^{+}A+W^{T}A^{*})$ (47)

where the auxiliary vectors (of size $N$) are

$U_{\nu}=V_{\nu}\sqrt{2\over\omega_{\nu}},\qquad
W_{\nu}=V_{\nu}\sqrt{2\omega_{\nu}}$

which is convenient for memory allocation as it requires only products of one-
dimensional arrays at each collision. Free evolution (10) can be equally
implemented as a product if one-dimensional arrays since $\Omega$ is diagonal.
Therefore this scheme ensures exact energy and momentum conservation.

### Appendix B: Stochastic equations

It is useful to rewrite (11) as a stochastic equation. Writing $M\to Mdw$
where $dw(t)$ is a stochastic process, which take the value 0 or 1 at random
times, with the Poisson process satisfying $\langle
dw\rangle=dt/\langle\tau\rangle$. The infinitesimal changes in the variables
from time $t$ to $t+dt$ to leading order is

$\displaystyle dA=i\Omega Adt-M\left(A+A^{*}\right)dw,$ (48)

where $dw^{2}=dw$ has been used. For the translation-invariant case instead

$\displaystyle dA=i\Omega Adt-M\left(A+\tilde{A}^{*}\right)dw,$ (49)

(see the main text), with the same meaning of $dw$.

Equation (48) can be used to determine the equation for the action vector
$I=A\circ A^{*}$, where $\circ$ is the Hadamard product (element-wise
$A_{\nu}A^{*}_{\nu}=I_{\nu}$) with stochastic calculus. Dropping $dt^{2},dtdw$
terms and using again $dw^{2}=dw$:

$\displaystyle d(A\circ A^{*})=dA\circ A^{*}+dA^{*}\circ A+\frac{1}{2}dA\circ
dA^{*}$ (50)
$\displaystyle=[-M\left(A+A^{*}\right)\circ\left(A+A^{*}\right)+M\left(A+A^{*}\right)\circ
M\left(A+A^{*}\right)]dw,$

which contains quadratic, non-diagonal terms. For completeness, we mention
that stochastic equations have been discussed before for conservative noise,
for instance for weakly interacting anharmonic oscillators [77]. Here,
however, we deal with stochastic dynamics in action space.

In the kinetic limit, we proceed as in the main text. We first perform
averaging over uniform random phases so that only diagonal terms survive
$A_{\nu}A^{*}_{\nu}=I_{\nu}$ in the right-hand side. We get the stochastic
equation for the actions

$dI=2\left[M\circ M-\mathrm{diag}(M)\right]I\,dw,$ (51)

with $\mathrm{diag}(M)$ being the diagonal matrix with diagonal elements
$M_{\nu,\nu}$ (again the same symbol $I$ is used for the averages).
Alternatively, we can work with the the mode energies defined, in matrix
notation, as $E=\Omega I$:

$dE=KE\,dw.$ (52)

Equations (14) and (15) are of course equivalent to the above stochastic
differential equation formulation.

## References

  * [1] Giovanni Gallavotti.  The Fermi-Pasta-Ulam problem: a status report, volume 728.  Springer, 2007.
  * [2] G. Benettin, H. Christodoulidi, and A. Ponno.  The Fermi-Pasta-Ulam problem and its underlying integrable dynamics.  Journal of Statistical Physics, 152(2):195–212, Jul 2013.
  * [3] Wojciech De Roeck and Francois Huveneers.  Asymptotic localization of energy in nondisordered oscillator chains.  Commun. Pure Appl. Math., 68(9):1532–1568, 2015.
  * [4] Francois Huveneers.  Classical and quantum systems: transport due to rare events.  Annalen der Physik, 529(7):1600384, 2017.
  * [5] Weicheng Fu, Yong Zhang, and Hong Zhao.  Nonintegrability and thermalization of one-dimensional diatomic lattices.  Phys. Rev. E, 100(5):052102, 2019.
  * [6] Tomer Goldfriend and Jorge Kurchan.  Equilibration of quasi-integrable systems.  Physical Review E, 99(2):022146, 2019.
  * [7] Marco Baldovin, Angelo Vulpiani, and Giacomo Gradenigo.  Statistical mechanics of an integrable system.  Journal of Statistical Physics, 183(3):1–16, 2021.
  * [8] Abhishek Dhar.  Heat transport in low-dimensional systems.  Adv. Phys., 57:457–537, 2008.
  * [9] Stefano Lepri, editor.  Thermal transport in low dimensions: from statistical physics to nanoscale heat transfer, volume 921 of Lect. Notes Phys.  Springer-Verlag, Berlin Heidelberg, 2016.
  * [10] Giuliano Benenti, Stefano Lepri, and Roberto Livi.  Anomalous heat transport in classical many-body systems: Overview and perspectives.  Frontiers in Physics, 8:292, 2020.
  * [11] Herbert Spohn.  The phonon Boltzmann equation, properties and link to weakly anharmonic lattice dynamics.  Journal of statistical physics, 124(2):1041–1104, 2006.
  * [12] Andrey Pereverzev.  Fermi-Pasta-Ulam $\beta$ lattice: Peierls equation and anomalous heat conductivity.  Phys. Rev. E, 68(5):056124, Nov 2003.
  * [13] Bernie Nickel.  The solution to the 4-phonon Boltzmann equation for a 1d chain in a thermal gradient.  J. Phys. A-Math. Gen., 40(6):1219–1238, FEB 9 2007.
  * [14] J. Lukkarinen and H. Spohn.  Anomalous energy transport in the FPU-$\beta$ chain.  Communications on Pure and Applied Mathematics, 61(12):1753–1786, 2008.
  * [15] Jani Lukkarinen.  Kinetic theory of phonons in weakly anharmonic particle chains.  In Thermal transport in low dimensions, pages 159–214. Springer, 2016.
  * [16] Miguel Onorato, Lara Vozella, Davide Proment, and Yuri V. Lvov.  Route to thermalization in the $\alpha$-Fermi-Pasta-Ulam system.  Proceedings of the National Academy of Sciences, 112(14):4208–4213, 2015.
  * [17] Francois Huveneers and Jani Lukkarinen.  Prethermalization in a classical phonon field: Slow relaxation of the number of phonons.  Physical review research, 2(2):022034, 2020.
  * [18] C Kipnis, C Marchioro, and E Presutti.  Heat flow in an exactly solvable model.  J. Stat. Phys., 27:65, 1982.
  * [19] Bernard Derrida.  An exactly soluble non-equilibrium system: the asymmetric simple exclusion process.  Phys. Rep., 301(1):65–83, 1998.
  * [20] A. Malevanets and R. Kapral.  Mesoscopic model for solvent dynamics.  The Journal of chemical physics, 110:8605–8613, May 1999.
  * [21] R. Kapral.  Multiparticle Collision Dynamics: Simulation of Complex Systems on Mesoscales, pages 89–146.  2008\.
  * [22] Pierfrancesco Di Cintio, Roberto Livi, Stefano Lepri, and Guido Ciraolo.  Multiparticle collision simulations of two-dimensional one-component plasmas: Anomalous transport and dimensional crossovers.  Phys. Rev. E, 95:043203, Apr 2017.
  * [23] G Basile, C Bernardin, and S Olla.  Momentum conserving model with anomalous thermal conductivity in low dimensional systems.  Phys. Rev. Lett., 96:204303, 2006.
  * [24] G. Basile, L. Delfini, S. Lepri, R. Livi, S. Olla, and A. Politi.  Anomalous transport and relaxation in classical one-dimensional models.  Eur. Phys J.-Special Topics, 151:85–93, 2007.
  * [25] Giada Basile, Cedric Bernardin, Milton Jara, Tomasz Komorowski, and Stefano Olla.  Thermal conductivity in harmonic lattices with random collisions.  In Thermal transport in low dimensions, pages 215–237. Springer, 2016.
  * [26] Cedric Bernardin, Venkateshan Kannan, Joel L Lebowitz, and Jani Lukkarinen.  Harmonic systems with bulk noises.  Journal of Statistical Physics, 146(4):800–831, 2012.
  * [27] Giada Basile, Stefano Olla, and Herbert Spohn.  Energy transport in stochastically perturbed lattice dynamics.  Archive for rational mechanics and analysis, 195(1):171–203, 2010\.
  * [28] Jani Lukkarinen, Matteo Marcozzi, and Alessia Nota.  Harmonic chain with velocity flips: thermalization and kinetic theory.  Journal of Statistical Physics, 165(5):809–844, 2016.
  * [29] S. Lepri, C. Mejia-Monasterio, and A. Politi.  Stochastic model of anomalous heat transport.  J. Phys. A: Math. Theor., 42:025001, 2009.
  * [30] S. Lepri, C. Mejia-Monasterio, and A. Politi.  Dynamics of anomalous heat transport.  J. Phys. A: Math. Theor., 43:065002, 2010.
  * [31] L Delfini, S Lepri, R Livi, C Mejia-Monasterio, and A Politi.  Nonequilibrium dynamics of a stochastic model of anomalous heat transport: numerical analysis.  J. Phys. A: Math. Theor., 43(14):145001, 2010.
  * [32] Aritra Kundu, Cedric Bernardin, Keji Saito, Anupam Kundu, and Abhishek Dhar.  Fractional equation description of an open anomalous heat conduction set-up.  J. Stat. Mech: Theory Exp., 2019(1):013205, 2019.
  * [33] Stefano Lepri and Antonio Politi.  Density profiles in open superdiffusive systems.  Phys. Rev. E, 83(3):030107, 2011.
  * [34] Alessandra Iacobucci, Frederic Legoll, Stefano Olla, and Gabriel Stoltz.  Thermal conductivity of the Toda lattice with conservative noise.  J. Stat. Phys., 140(2):336–348, Jul 2010.
  * [35] C Bernardin and P Goncalves.  Anomalous fluctuations for a perturbed hamiltonian system with exponential interactions.  Communications in Mathematical Physics, 325(1):291–332, 2014.
  * [36] Stefano Lepri, Roberto Livi, and Antonio Politi.  Too close to integrable: Crossover from normal to anomalous heat diffusion.  Physical Review Letters, 125(4):040604, 2020.
  * [37] A. S. Pikovsky and D. L. Shepelyansky.  Destruction of Anderson localization by a weak nonlinearity.  Physical Review Letters, 100(9):094101, 2008.
  * [38] G. Kopidakis, S. Komineas, S. Flach, and S. Aubry.  Absence of wave packet diffusion in disordered nonlinear systems.  Physical Review Letters, 100(8):084103, 2008.
  * [39] Ch. Skokos, D. O. Krimer, S. Komineas, and S. Flach.  Delocalization of wave packets in disordered nonlinear chains.  Physical Review E - (Statistical, Nonlinear, and Soft Matter Physics), 79(5):056211, 2009.
  * [40] S Lepri, R Schilling, and S Aubry.  Asymptotic energy profile of a wave packet in disordered chains.  Physical Review E, 82(5):056602, 2010.
  * [41] D.M. Basko.  Weak chaos in the disordered nonlinear Schroedinger chain: Destruction of Anderson localization by Arnold diffusion.  Annals of Physics, 326(7):1577–1655, 2011.
  * [42] Manoj Kumar, Anupam Kundu, Manas Kulkarni, David A Huse, and Abhishek Dhar.  Transport, correlations, and chaos in a classical disordered anharmonic chain.  Physical Review E, 102(2), 2020.
  * [43] Jean-Philippe Bouchaud and Antoine Georges.  Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications.  Physics reports, 195(4-5):127–293, 1990.
  * [44] Monique M Tirion.  Large amplitude elastic motions in proteins from a single-parameter, atomic analysis.  Physical review letters, 77(9):1905, 1996.
  * [45] Brice Juanico, Y-H Sanejouand, Francesco Piazza, and Paolo De Los Rios.  Discrete breathers in nonlinear network models of proteins.  Physical review letters, 99(23):238104, 2007.
  * [46] Nahuel Freitas and Juan Pablo Paz.  Analytic solution for heat flow through a general harmonic network.  Physical Review E, 90(4):042128, 2014.
  * [47] Kezhao Xiong, Chunhua Zeng, Zonghua Liu, and Baowen Li.  Influence of the degree of a complex network on heat conduction.  Physical Review E, 98(2):022115, 2018.
  * [48] G. Benettin and A. Ponno.  Time-scales to equipartition in the Fermi–Pasta–Ulam problem: Finite-size effects and thermodynamic limit.  Journal of Statistical Physics, 144(4):793, Aug 2011.
  * [49] S Lepri.  Relaxation of classical many-body hamiltonians in one dimension.  Phys. Rev. E, 58(6):7165–7171, 1998.
  * [50] S Lepri.  Memory effects and heat transport in one-dimensional insulators.  Eur. Phys J. B, 18(3):441–446, DEC 2000.
  * [51] Thudiyangal Mithun, Yagmur Kati, Carlo Danieli, and Sergej Flach.  Weakly nonergodic dynamics in the Gross-Pitaevskii lattice.  Physical review letters, 120(18):184101, 2018.
  * [52] Carlo Danieli, Thudiyangal Mithun, Yagmur Kati, David K. Campbell, and Sergej Flach.  Dynamical glass in weakly nonintegrable Klein-Gordon chains.  Phys. Rev. E, 100:032217, Sep 2019.
  * [53] Arkady Pikovsky and Antonio Politi.  Lyapunov exponents: a tool to explore complex dynamics.  Cambridge University Press, 2016.
  * [54] Jürgen Schnakenberg.  Network theory of microscopic and macroscopic behavior of master equation systems.  Reviews of Modern physics, 48(4):571, 1976.
  * [55] Raffaella Burioni and Davide Cassi.  Universal properties of spectral dimension.  Physical review letters, 76(7):1091, 1996.
  * [56] Raffaella Burioni and Davide Cassi.  Random walks on graphs: ideas, techniques and results.  Journal of Physics A: Mathematical and General, 38(8):R45, 2005\.
  * [57] T Prosen and M Robnik.  Energy-transport and detailed verification of Fourier heat law in a chain of colliding harmonic-oscillators.  J. Phys. A-Math. Gen., 25(12):3449–3472, 1992.
  * [58] A Pikovsky.  Scaling of energy spreading in a disordered ding-dong lattice.  Journal of Statistical Mechanics: Theory and Experiment, 2020(5):053301, 2020.
  * [59] Shuji Tamaki and Keiji Saito.  Energy current correlation in solvable long-range interacting systems.  Physical Review E, 101(4):042118, 2020.
  * [60] H. Matsuda and K. Ishii.  Localization of normal modes and energy transport in the disordered harmonic chain.  Prog. Theor. Phys. Suppl., 45:76, 1970.
  * [61] WM Visscher.  Localization of normal modes and energy transport in disordered harmonic chain.  Prog. Theor. Phys., 46(3):729–&, 1971.
  * [62] Cedric Bernardin.  Thermal conductivity for a noisy disordered harmonic chain.  Journal of Statistical Physics, 133(3):417–433, 2008.
  * [63] Abhishek Dhar, K Venkateshan, and JL Lebowitz.  Heat conduction in disordered harmonic lattices with energy-conserving noise.  Physical Review E, 83(2):021108, 2011.
  * [64] Cedric Bernardin, Francois Huveneers, and Stefano Olla.  Hydrodynamic limit for a disordered harmonic chain.  Communications in Mathematical Physics, 365(1):215–237, 2019.
  * [65] Jun Yan.  Harmonic interaction model and its applications in Bose–Einstein condensation.  Journal of statistical physics, 113(3):623–634, 2003.
  * [66] Lucianno Defaveri, Carlos Olivares, and Celia Anteneodo.  Heat flux in chains of nonlocally coupled harmonic oscillators: Mean-field limit.  Phys. Rev. E, 105:054149, May 2022.
  * [67] Francesco Andreucci, Stefano Lepri, Stefano Ruffo, and Andrea Trombettoni.  Classical and quantum harmonic mean-field models coupled intensively and extensively with external baths.  SciPost Physics Core, 5(3):036, 2022.
  * [68] MB Hastings.  Random vibrational networks and the renormalization group.  Physical review letters, 90(14):148702, 2003.
  * [69] Mark EJ Newman and Duncan J Watts.  Renormalization group analysis of the small-world network model.  Physics Letters A, 263(4-6):341–346, 1999.
  * [70] Remi Monasson.  Diffusion, localization and dispersion relations on small-world lattices.  The European Physical Journal B-Condensed Matter and Complex Systems, 12(4):555–567, 1999.
  * [71] DM Basko.  Kinetic theory of nonlinear diffusion in a weakly disordered nonlinear Schroedinger chain in the regime of homogeneous chaos.  Physical Review E, 89(2):022921, 2014.
  * [72] Roberto Livi, Marco Pettini, Stefano Ruffo, Massimo Sparpaglione, and Angelo Vulpiani.  Equipartition threshold in nonlinear large hamiltonian systems: The Fermi-Pasta-Ulam model.  Physical Review A, 31(2):1039, 1985.
  * [73] Oliver Mülken, Sarah Heinzelmann, and Maxim Dolgushev.  Information dimension of stochastic processes on networks: relating entropy production to spectral properties.  Journal of Statistical Physics, 167(5):1233–1243, 2017.
  * [74] Abhishek Dhar, Anupam Kundu, and Aritra Kundu.  Anomalous heat transport in one dimensional systems: a description using non-local fractional-type diffusion equation.  Frontiers in Physics, 7:159, 2019.
  * [75] Stefano Iubini.  Coupled transport in a linear-stochastic Schroedinger equation.  Journal of Statistical Mechanics: Theory and Experiment, 2019(9):094016, 2019.
  * [76] Noé Cuneo, Jean-Pierre Eckmann, Martin Hairer, and Luc Rey-Bellet.  Non-equilibrium steady states for networks of oscillators.  Electronic Journal of Probability, 23:1–28, 2018.
  * [77] Carlangelo Liverani and Stefano Olla.  Toward the Fourier law for a weakly interacting anharmonic crystal.  Journal of the American Mathematical Society, 25(2):555–583, 2012\.