# Long range correlations and slow time scales in a boundary driven granular
model

Andrea Plati  Dipartimento di Fisica, Università di Roma Sapienza, P.le Aldo
Moro 2, 00185 Roma, Italy  Andrea Puglisi  Istituto dei Sistemi Complessi -
CNR and Dipartimento di Fisica, Università di Roma Sapienza, P.le Aldo Moro 2,
00185, Rome, Italy

###### Abstract

We consider a velocity field with linear viscous interactions defined on a one
dimensional lattice. Brownian baths with different parameters can be coupled
to the boundary sites and to the bulk sites, determining different kinds of
non-equilibrium steady states or free-cooling dynamics. Analytical results for
spatial and temporal correlations are provided by analytical diagonalisation
of the system’s equations in the infinite size limit. We demonstrate that
spatial correlations are scale-free and time-scales become exceedingly long
when the system is driven only at the boundaries. On the contrary, in the case
a bath is coupled to the bulk sites too, an exponential correlation decay is
found with a finite characteristic length. This is also true in the free
cooling regime, but in this case the correlation length grows diffusively in
time. We discuss the crucial role of boundary driving for long-range
correlations and slow time-scales, proposing an analogy between this
simplified dynamical model and dense vibro-fluidized granular materials.
Several generalizations and connections with the statistical physics of active
matter are also suggested.

## Introduction

The emergence of long-range order, or collective behavior (CB), in non-
equilibrium systems such as granular materials and living organisms is a
matter of great interest for fundamental physics and applications [1, 2].
Examples, recently observed in experiments and numerical simulations, are
motility induced phase transitions in bacteria [3, 4, 5, 6], collective
migration in epithelial cells [7], persistent collective rotations in granular
systems [8, 9, 10]. An important class of CB instances includes flocking and
swarming in animals, systematically studied by physicists in the last 25 years
[11, 12, 13]. The great variety of systems in which CB has been observed makes
the formulation of a rigorous and unifying definition for them a difficult
task. Generally speaking we can say that CB occurs when a many-body system
_acts as a whole_. Indeed, a common property of the previous examples is the
interplay between different length scales: the interactions act on microscopic
distances while correlations extend to macroscopic scales, comparable with the
system size. In the study of CB it is common, in fact, to look at spatial
correlation functions of the relevant fields: if this function has a typical
decay length $\xi$ then we can divide the system in almost independent
subsystems of size $\sim\xi$. If the correlation function decays without a
typical length it is said to be _scale-free_ : in this case the dynamics of
every particle is correlated with the whole system. We underline that _scale-
free_ spatial correlations appear naturally in critical phase transitions at
equilibrium [14], but a general and well established theoretical framework to
understand the appearance of long-range ordering in non-equilibrium systems is
still lacking: sometimes equilibrium-like approaches are successful (effective
Hamiltonian/temperatures) [15, 16] while in other cases fully non-equilibrium
tools have to be developed [17, 18, 19].

In this paper, we provide analytical results about the occurrence of _scale-
free_ (more precisely power law decaying) correlations in a velocity field
defined on a one dimensional lattice with interactions mediated by viscous
friction. We’ll show that this behavior is observed in the non-equilibrium
stationary state (NESS) obtained by coupling only the boundaries of the system
with a thermal bath. We call this phase Non-Homogeneously Heated Phase (NHHP).
If the particles in the bulk are also put in contact with a bath a different
regime is found, the Homogeneously Heated Phase (HHP), where the spatial
correlation is exponential with a characteristic length scale that goes to
infinity when the contact between the bulk and the bath vanishes. The NHHP is
also characterized by slow relaxation times that scale with the square of the
system size.

Lattices (particularly in 1d) bring two main advantages: (i) analytical
calculations are often possible, (ii) they help to isolate minimal ingredients
for the occurrence of the phenomenon under study. Considering just the non-
equilibrium context, 1D models have been used to study thermal conduction [20,
21, 22], non-equilibrium fluctuations [23, 24], correlations and response with
non-symmetric couplings [25], velocity alignment in active matter [26],
systems with Vicsek-like interactions [27, 28], velocity fields in granular
materials [29, 30, 31]. In the following we’ll just consider linear
interactions between variables and this allows to work in the framework of
multivariate linear stochastic processes. Despite their simplicity, this class
of models continues to be a powerful tool when dealing with dynamics driven
out of equilibrium as in biological systems [32, 33].

As discussed in the next section, our model can be thought as an extreme
simplification of a vibrated granular system at strong compression. Looking
for the emergence of a collective motion in it is then motivated also by the
recent experimental/numerical evidence of slow collective behavior in vibro-
fluidized granular materials [8, 9]. This phenomenon is not yet fully
understood and our study tackles this problem, revealing that non-homogeneous
heating and frictional interactions (i.e standard features of vibrated
granular matter) are minimal ingredients to develop a slow collective
dynamics.

The manuscript is organized as follows: In section "Model" we present our
model discussing its phenomenology and its relation with real granular systems
and previously studied non-equilibrium 1D models. Section "Results" contains
the key-steps for the calculation of the spatial correlation function in the
NHHP and in the HHP shedding light on the limit for which diverging
correlation lengths and times are obtained. We also show the validity of our
results beyond the assumptions used to perform analytical calculations.
Finally, in "Discussion" we draw conclusions and sketch some perspectives. In
the Supplemental Material (SM) details of the calculations are provided in
addition to some insights about the cooling state and the active equivalent of
our model.

## Model

### Definition and phenomenology

We consider a velocity field on a one dimensional lattice of size $L$. The
$i$th particle interacts with their nearest neighbors $j$ through a viscous
force with coefficient $\gamma$: $F_{i}=-\sum_{j}\gamma(v_{i}-v_{j})$. The
boundary (bulk) sites are coupled with an external bath defined by a drag
coefficient $\gamma_{b}$ ($\gamma_{a}$) and relative temperatures which can be
different if at the boundaries or in the bulk. Considering particles with
unitary mass the equations for the model are:

$\displaystyle\dot{v}_{i}=-(2\gamma+\gamma_{a})v_{i}+\gamma(v_{i+1}+v_{i-1})+\sqrt{2\gamma_{a}T_{a}}\eta_{i}(t)$
(1a) $\displaystyle\dot{v}_{1}=-(\gamma+\gamma_{b})v_{1}+\gamma
v_{2}+\sqrt{2\gamma_{b}T_{1}}\eta_{1}(t)$ (1b)
$\displaystyle\dot{v}_{L}=-(\gamma+\gamma_{b})v_{L}+\gamma
v_{L-1}+\sqrt{2\gamma_{b}T_{L}}\eta_{L}(t)$ (1c)

Where the first equation holds for $1<i<L$ and the $\eta_{i}(t)$s are Gaussian
white noises with unitary variance:
$\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$.

In this model, the way in which energy is supplied to the system is consistent
with the fluctuation-dissipation theorem. Indeed, for each viscous force
($\gamma_{a(b)}$) there is a stochastic counterpart at finite temperature
($T_{a(b)}$). This is actually not true for the interaction force defined by
$\gamma$ because it is related to the viscosity of the material that forms the
grains. Thus, the associated temperature (typical of the thermal agitation at
the molecular scale) can be reasonably neglected in a granular context. We
refer to NHHP when $\gamma_{a}=0$ so that just the first and the $L$th sites
are heated, while in the HHP we consider a general $\gamma_{a}\neq 0$. We note
that the HHP is not strictly spatially homogeneous because viscous
coefficients and temperatures depend on the position: we refer to it as
homogeneously heated meaning that in this phase _all_ the particles are
coupled with a bath.

Figure 1: a,b,c) Snapshots of the velocity field in the stationary state of
the two phases. We exclude the first five (really hot) sites near the
boundaries to have a more clear view of the field. Each panel shows the
vectors in linear scale and the moduli in log scale in order to better
appreciate the phenomenology of the system. Orange and blue bars discriminate
the two directions. We note that a great cluster of particles with same
direction and similar modulus is found in the NHHP only, signaling that in
terms of correlations the key parameter is $\gamma_{a}$ rather than $T_{a}$.
d) Autocorrelation times for each site defined as the time $\tau_{i}$ for
which $\Gamma_{i}(\tau_{i})=0.4$. The autocorrelation function is defined as
$\Gamma(t^{\prime})=\lim_{t\to\infty}\langle
v_{i}(t)v_{i}(t+t^{\prime})\rangle/\langle v_{i}^{2}(t)\rangle$ where the
brackets refer to a time average on the stationary state. We note that in the
NHHP the dynamics is far slower than in the HHP also when $T_{a}=0$. The
snapshots are obtained by numerical integration of Eqs. (1) with $L=50$,
$\gamma=5$, $\gamma_{b}=10$, $\gamma_{a}=\\{3,0\\}$, $T_{1}=T_{L}=0.002$,
$T_{a}=\\{0.002,0\\}$ after a time $t_{M}=10^{8}/\gamma$ and with a temporal
step $dt=0.05/\gamma$.

As we discuss in the next paragraphs, this is a linear model and a full
solution can be found in the context of multivariate stochastic processes.
Nevertheless, a numerical integration of Eqs. (1) can be useful to have a
physical insight on the phenomenology in play. In Fig. 1 we show some
instantaneous snapshots of the system in the stationary state for three
different conditions: HHP with $T_{a}\neq 0$, HHP with $T_{a}=0$ and NHHP. We
note that in the NHHP (panel c) almost all the velocities are aligned with
similar moduli while in the HHP we have smaller aligned domains with moduli
that decay sharply moving away from the boundaries when $T_{a}=0$ (panel b)
and a random configuration when $T_{a}\neq 0$ (panel a). This comparison makes
clear that - in terms of correlations - the key parameter is $\gamma_{a}$
rather than $T_{a}$: indeed a situation where the sites experience a
collective behavior (in the intuitive sense that they _act as a whole_) is
only found in the NHHP. In Fig. 1d the typical correlation time for each site
is shown and we can see that in the NHHP the dynamics is extremely slower with
respect to the other two conditions. It is worth noting that this model does
not present any directional asymmetry so the true mean value of the velocity
field (i.e. obtained by an average over long times or equivalently over all
the realizations of the noises) is zero also in the NHHP, even if the single
time configurations clearly show an explicit global alignment. The
phenomenology of the NHHP can then be described as the occurrence of slow and
collective fluctuations around the expected mean value.

### Relation with real granular systems and other models

We note that the kind of interaction used in Eqs. 1 is typical of contact
models for granular materials [34, 35]. In these models, the grains (that are
disks or spheres depending on the geometry) interact when a distance smaller
than the sum of their radius is reached. In this condition, the particles
penetrate each other and the dynamics is ruled by contact forces that are
split into a normal and tangential component with respect to the vector
connecting the centers of the grains. Both of this contributions contain a
(linear or non-linear) elastic term that depends on the normal/tangential
displacement and a dissipative one that depends on the normal/tangential
relative velocity. The latter has, in many cases, exactly the form of the
viscous interaction we use in our model [36]. In view of this we can say that
if we fix the centers of $L$ grains on the lattice sites so that they are
partially overlapped, then the dynamics of the particles’ velocities would be
given by Eqs. 1. Neglecting the dynamics of positions (they don’t appear at
all in Eqs. 1) is surely the most relevant approximation of our approach: in
the SM (S5) we briefly discuss how to go beyond it.

Figure 2: Sketch of the model and relation with higher dimensional systems. On
the left we suggest a hypothetical 2D dense granular system where particles
are roughly located on the vertices of a regular lattice. A possible mapping
from the 2D to the 1D system involves replacing the mean horizontal velocity
on the $i$th layer of the 2D system and replacing it with the $v_{i}$ of the
1D system. The dynamics in the vertical direction is neglected, an
approximation which is justified by the presence of the vertical confinement,
while the periodic boundary conditions (indicated by the dotted lines) are
representative of a ’free’ direction in which the grains can flow without
obstacles. This can be realized experimentally, for instance, in a 3D
cylindrical geometry, where the velocity of grains in the tangential direction
(with respect to the central axis of the cylinder) constitute the horizontal
velocities in the putative 2D system sketched here, see for instance [8, 9].
Red grains are in direct contact with the external source of energy coming
from the boundaries ($\gamma_{b}$,$T_{1(L)}$) while the green ones are in
contact with the bulk bath, which is switched off in the NHHP.

Nevertheless, the physics described by our model can realistically represent
the condition of permanent contacts in which dense granular matter is found in
vertically-vibrated setups. Such kind of systems are widely studied
experimentally; they consist on assemblies of grains confined in a box
vibrated with a noisy or sinusoidal signal on the $z$ direction. For low
driving energies, the particles are always arranged in a dense packing where
they vibrate in permanent contact with each other experiencing very rare and
slow rearrangements. This implies, if the geometry is narrow enough, that just
the external layers of the system are in direct contact with the vibrating
walls while the others never touch them. This last fact tell us that, in
addition to the specific form of the viscous forces and the permanent
interactions, also the way in which the external energy injection is modeled
in the NHHP resembles the conditions of a vibrated granular system in a dense
state. Moreover, if layers of particles are mapped into lattice sites, a 1D
chain can also be representative of a higher dimensional systems (see Fig. 2).
On the other hand, the HHP can be referred to a setup where all the particles
interact with the vibrating walls, as it happens for instance in vibrated
monolayers [37].

The idea of considering velocity fields defined on lattices, i.e. neglecting
the evolution of the positions and density fluctuations in the dynamics, has
been widely exploited in granular literature [29, 30, 31] especially for
dilute systems. In these previous works, however, there is no continuous
interaction, but only instantaneous collisions occurring between pairs of
neighboring grains picked up at random, at every time step. Many results have
been obtained by solving (analytically or numerically) the corresponding
master equation or performing its hydrodynamic limit, revealing that these
models are a powerful tool to investigate complex phenomena observed in
experiments and simulations of realistic granular systems such as shock waves,
anomalous transport and current fluctuations [38, 39].

To summarize motivations and background, our model reflects three main
characteristics of dense granular materials in vertically-vibrated setups i.e.
viscous forces, permanent contacts and energy injection localized at the
boundaries. It can be then considered as the high density variant of a well
established family of models previously investigated.

It is important to note that also the dilute models can exhibit long-range
correlations [38, 39]. Nevertheless, those are finite-size effects found in
the homogeneous cooling state [40] i.e. without external driving and with
conserved total momentum. As we briefly discuss in the next paragraph and more
clearly in the SM (S4), our model makes clear that there is a sharp difference
between the correlations of the cooling state and the NESS ones.

### Compact SDE formulation of the model

Defining the vectors $\bm{V}=(v_{1},\dots,v_{L})$ ,
$\bm{\eta}(t)=(\eta_{1}(t),\dots,\eta_{L}(t))$ and the adimensional parameters
$\beta=\gamma_{b}/\gamma$, $\alpha=\gamma_{a}/\gamma$ then we can rewrite Eqs.
1 as a multivariate Ornstein-Uhlenbeck process obtaining the following
stochastic differential equation (SDE):

$\dot{\bm{V}}=-\hat{A}\bm{V}+\hat{B}\bm{\eta}(t)$ (2)

where
$\hat{B}=\text{diag}(\sqrt{2\gamma_{b}T_{1}},\sqrt{2\gamma_{a}T_{a}},\dots,\sqrt{2\gamma_{a}T_{a}},\sqrt{2\gamma_{b}T_{L}})$
and:

$\hat{A}=\gamma\begin{pmatrix}1+\beta&-1&&&\bm{0}\\\ -1&2+\alpha&-1&&\\\
&\ddots&\ddots&\ddots&\\\ &&-1&2+\alpha&-1\\\
\bm{0}&&&-1&1+\beta\end{pmatrix}$ (3)

is a $L\times L$ tridiagonal symmetric matrix.

The information about space-time correlations of the system are encoded in the
two times correlation matrix $\hat{\sigma}(t,s)$ whose entries are defined as
$\sigma_{jm}(t,s)=\langle
v_{j}(t)v_{m}(s)\rangle\equiv\langle\left[v_{j}(t)-\langle
v_{j}(t)\rangle\right]\left[v_{m}(s)-\langle v_{m}(s)\rangle\right]\rangle$.
We now define the quantity of principal interest in this paper i.e. the static
spatial correlation function of the velocity field:

$\zeta_{jm}=\frac{\sigma_{jm}}{\sqrt{\sigma_{jj}\sigma_{mm}}}\quad\text{where}\quad\sigma_{jm}=\left\langle
v_{j}v_{m}\right\rangle.$ (4)

With this definition we have $\zeta_{jm}=1$ if $j=m$ or $v_{j}=v_{m}$ and
$\zeta_{jm}=0$ if $\langle v_{j}v_{m}\rangle=0$. It is then clear that our
goal is to solve Eq. 2 and find the stationary correlation matrix
$\hat{\sigma}=\lim_{t\to\infty}\hat{\sigma}(t,t)$ that exists if $\hat{A}$ is
positive semi-definite. In this conditions, regardless the symmetry of
$\hat{A}$, the correlation matrix can be found by inverting the relation [41]:

$\hat{A}\hat{\sigma}+\hat{\sigma}\hat{A}^{T}=\hat{B}\hat{B}^{T}.$ (5)

Nevertheless, a more direct way to obtain an analytic expression of
$\hat{\sigma}$ can be followed exploiting the fact that $\hat{A}$ is
symmetric. In this case there exist a unitary matrix $\hat{S}$ such that
$\hat{S}\hat{S}^{+}=\hat{I}$ and
$\hat{S}^{+}\hat{A}\hat{S}$=$\hat{S}^{+}\hat{A}^{T}\hat{S}=\hat{\lambda}=\text{diag}(\lambda_{1},\lambda_{2},\dots,\lambda_{L})$
where $\hat{I}$ is the identity matrix, the $\lambda_{j}$s are the eigenvalues
of $\hat{A}$ while $S_{ji}$ is the $j$th component of the $i$th eigenvector of
it. With these hypotheses and in the case of
$\hat{B}=\text{diag}(b_{1},\dots,b_{L})$ we can write the covariance matrix in
the two-times (with $t\geq s$) and non-stationary case:

$\hat{\sigma}(t,s)=\hat{S}\left(\hat{C}(t,s)+\hat{G}(t,s)\right)\hat{S}^{+}$
(6)

where:

$\displaystyle\hat{C}(t,s)=\exp(-\hat{\lambda}t)\hat{S}^{+}\langle\bm{V}(0),\bm{V}^{T}(0)\rangle\hat{S}\exp(-\hat{\lambda}s)$
(7a) $\displaystyle
G_{jm}(t,s)=\frac{\left(e^{-\lambda_{j}(t-s)}-e^{-(\lambda_{j}+\lambda_{m})s}\right)\sum_{n}S_{jn}^{+}S_{nm}b_{n}^{2}}{\lambda_{j}+\lambda_{m}}.$
(7b)

The first matrix represents the transient and the brackets refer to the
average over initial conditions while the NESS is described by
$\lim_{s\to\infty}G(t,s)$. Without noises, Eq. (7a) would be the solution of
Eq. 2 representing the correlations in the cooling state. We note that the two
correlation matrices has a different mathematical structure. The consequences
of that together with some properties of the cooling state are discussed in
the SM (S4) while in the next paragraphs we will neglect $\hat{C}$
concentrating on the NESS. Defining
$\hat{\sigma}(t^{\prime})=\lim_{t\to\infty}\hat{\sigma}(t+t^{\prime},t)$ and
through Eqs. (6) and (7b) it is also possible to evaluate the single particle
autocorrelation function
$\Gamma_{j}(t^{\prime})\equiv\sigma_{jj}(t^{\prime})/\sigma_{jj}(0)$:

$\Gamma_{j}(t^{\prime})=\frac{1}{\sigma_{jj}}\sum_{k}q_{jk}S^{+}_{kj}e^{-\lambda_{k}t^{\prime}},\quad
q_{jk}=\sum_{ls}\frac{S_{jl}S^{+}_{ks}S_{sl}b^{2}_{s}}{\lambda_{l}+\lambda_{k}}$
(8)

from which is clear that, as expected for a linear system, the autocorrelation
function is a sum of exponential terms with different characteristic times
that are given by the inverse of the eigenvalues $\tau_{k}=1/\lambda_{k}$.

We will derive $\sigma_{jm}$ in a specific case where the diagonalisation of
$\hat{A}$ can be done analytically and then follow a numerical technique of
diagonalisation [42] to show the robustness of our main results i.e. power law
decay of spatial correlations. Before doing that, we briefly review what
techniques have been used to solve similar problems highlighting the
differences with the present case.

These kinds of lattice models, and also more complex ones (with higher
dimension and second order dynamics), when translational invariance holds, can
be mapped in a system of independent equations for the modes in the Bravais
lattice allowing a full solution [6]. However, our model (both NHHP and HHP)
has not periodic boundary conditions and the bath parameters depend on the
particular site position. Assuming translational invariance would mean giving
up some crucial aspects of our investigation. To keep a reasonable connection
with dense granular matter it is important to have a source of energy that
acts differently at the boundary and in the bulk of the system. Nevertheless,
in the next section we’ll discuss some common aspects between the HHP and
translational invariant systems.

We also point out that the continuous limit of Eq. 1a leads to the following
equation for the velocity field:
$\partial_{t}v(x,t)=-\gamma_{a}v(x,t)+\partial_{xx}v(x,t)+\sqrt{2T_{a}\gamma_{a}}\xi(x,t)$
with
$\langle\xi(x,t)\xi(x^{\prime},t^{\prime})\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime})$.
Equations of this form applied on a density field describe a diffusion process
with traps and noise. The variation of the field at the point $x$ is indeed
given by a noise, a diffusive term and a loss term ($-\gamma_{a}v(x,t)$) that
represents the possibility for the particles to be permanently trapped. These
processes can be used to describe the dynamics of mobile defects in crystals
where translational invariance is assumed and the problem can be easily solved
in Fourier space [43]. Our case where external thermostats are necessary to
keep stationary the system and break translational invariance is different. In
the general case with space-dependent parameters, correlations can be studied
diagonalising the matrix $\hat{A}$ or by exploiting Eq. 5 combined with
physical constraint on $\hat{\sigma}$. The former strategy, used by us and
recently applied in [25, 22], when possible, is more convenient because it
gives access also to time-dependent properties. The latter has been used to
study temperature profiles in non-equilibrium harmonic chains [20] . It is
important to stress that a crucial difference between the present work and the
aforementioned ones is that we deal with interactions acting on relative
velocities and not (only) on displacements. Indeed, we have a direct
competition between baths $\gamma_{a(b)}$ and interaction $\gamma$ in
$\hat{A}$, while in heated harmonic chains only the coupling constants appear
in the interaction matrix.

#### Toeplitz condition

In order to obtain an explicit form of Eq. 6 we consider the case of
$\gamma_{b}=\gamma+\gamma_{a}$ so that $\beta=1+\alpha$ making $\hat{A}$ a
Toeplitz matrix whose eigenvalues and eigenvectors are respectively:

$\lambda_{j}=\gamma(2+\alpha-2\cos(j\Pi)),\quad
S_{jm}=\sqrt{\frac{2\Pi}{\pi}}\sin\left(jm\Pi\right)$ (9)

where $\Pi=\pi/(L+1)$. Replacing these in Eq. (7b) and taking $t=s\to\infty$,
Eq. 6 becomes:

$\sigma_{jm}(\alpha)=\frac{2\Pi^{2}}{\gamma\pi^{2}}\sum_{lk}\frac{\sin\left(jl\Pi\right)\sin\left(mk\Pi\right)\left[\sum_{n}b^{2}_{n}\sin(ln\Pi)\sin(kn\Pi)\right]}{\Delta(\alpha)-\cos\left(k\Pi\right)-\cos\left(l\Pi\right)},$
(10)

where $\Delta(\alpha)=2+\alpha$. The sums run from 1 to $L$ and:

$b_{n}^{2}=\begin{cases}2(\gamma+\gamma_{a})T_{1},\quad n=1\\\
2\gamma_{a}T_{a},\quad 1<n<L\\\ 2(\gamma+\gamma_{a})T_{L},\quad n=L.\\\
\end{cases}$ (11)

We point out that Eq. (10) is symmetric with respect the center of the lattice
(i.e. $\sigma_{1m}=\sigma_{L(L+1-m)}$) if the coefficients $b_{n}$ are too.

## Results

### Power-Law correlations and slow time scales in the NHHP

We first study the NHHP so we put $\gamma_{a}=0$ and use the Toeplitz
condition that now reads $\gamma_{b}=\gamma$ so $\beta=1$. Exploiting the
limit for large systems ($L\gg 1$), we can exchange sums with integrals as
$\Pi\sum_{k=1}^{k=L}f(k\Pi)\rightarrow\int_{0}^{\pi}dzf(z)$. We note that in
Eq. (10), when $\gamma_{a}=0$, the sum over $n$ is actually made of two terms.
The one multiplied by $\gamma_{b}T_{L}$ has a sign that depends on the parity
of $l$ and $k$ and this brings to a subleading contribution if one considers
$L\gg 1$ and $j,m\ll L$ (see S1 in the SM). Neglecting it and defining

$\Sigma_{jm}(\alpha)=\int_{0}^{\pi}dzds\frac{\sin(jz)\sin(ms)\sin(z)\sin(s)}{\Delta(\alpha)-\cos\left(z\right)-\cos\left(s\right)}$
(12)

we obtain the covariance matrix for the NHHP:

$\sigma^{\text{NHHP}}_{jm}=\frac{4T_{1}}{\pi^{2}}\Sigma_{jm}(0).$ (13)

The integral contained in $\Sigma_{jm}(0)$ is difficult to be explicitly
evaluated but the following asymptotic behaviors can be derived in the limit
$L\gg m\gg 1$:

$\displaystyle\sigma^{\text{NHHP}}_{mm}\sim\frac{1}{m^{2}}$ (14a)
$\displaystyle\sigma^{\text{NHHP}}_{1m}\sim\frac{8T_{1}}{\pi m^{3}}$ (14b)
$\displaystyle\zeta^{\text{NHHP}}_{1m}\sim\frac{1}{m^{2}}$ (14c)

As explained in the SM (S2), these results are obtained by expressing
$\sigma^{\text{NHHP}}_{jm}$ as a power series of $(jm)^{-1}$ by multiple
integration by parts and estimating opportune upper bounds. The limit $L\gg
m\gg 1$ is important because we want to study the asymptotic behavior of the
correlations in the range for which they are not affected by the opposite
boundary of the system. This is the reason why we predict just a decay for the
variance $\sigma_{mm}$ even if it must grow approaching the $L$th site if
$T_{L}\neq 0$. This growth for large $m$ is given by the term proportional to
$\gamma_{b}T_{L}$ that we have neglected going from Eq. (10) to Eq. (13).

Eq. (14c) clearly states that the bulk sites are correlated with the first
(heated) one by a power law decay with exponent 2. Regarding the correlations
between particles in the bulk, they show a decay even slower than a power law.
We discuss them in the last paragraph of this section. Regarding time scales,
looking at Eq. (8) and at the specific form of the eigenvalues of $\hat{A}$ in
Eq. (9) for $\alpha=0$, we see that, when $j/L\ll 1$, the slowest time scales
in the single particle autocorrelation function behave as:

$\tau_{j}^{\text{NHHP}}=1/\lambda_{j}\sim\tau L^{2}$ (15)

where $\tau=1/\gamma$. We note that the emergence of characteristic times that
scale with the system size together with _scale free_ correlations is fully
consistent. Thus, the information that influences the dynamics of every
particle comes from all across the system and so the time to receive it must
increase with the system size.

### Finite Correlation Length and Times in the HHP

The emergence of _scale free_ correlations is often considered a remarkable
fact in physical systems. Nevertheless, we are now dealing with a model so it
is important to understand if this result is found just by an algebraic
coincidence or if it is consistent with the usual framework in which _scale
free_ correlations are understood i.e. a particular limit for which a finite
correlation length diverges. The study of the HHP comes into play to provide
an evidence of this last scenario. We point out that by studying the HHP with
periodic boundary conditions, and therefore assuming translational invariance
(i.e. extending Eq. 1a to all the particles in the system), it is quite easy
to derive an exponential decay for the stationary spatial correlation
function. This can be done by expressing Eq. 1a in the Bravais lattice or by
studying the continuous limit of $\dot{\sigma}_{jm}=\langle
v_{j}\dot{v}_{m}+v_{m}\dot{v}_{j}\rangle=0$. Nevertheless, we want to study
the passage from the HHP to the NHHP when $\gamma_{a}\to 0$ so we proceed with
space dependent parameters from Eq. (10). This expression, in the HHP,
contains all the contributions given by Eq. (11). Performing the large system
limit and taking into account just the leading terms we arrive at the
following expression for the covariance matrix in the HHP (see S3 in the SM
for details):

$\begin{split}\sigma_{jm}^{\text{HHP}}(\alpha)=\frac{2\alpha
T_{a}}{\pi}\int_{0}^{\pi}dz\frac{\sin(jz)\sin(mz)}{\Delta(\alpha)-2\cos(z)}+\frac{4T_{1}}{\pi^{2}}\left[1+\alpha\left(1-\frac{T_{a}}{T_{1}}\right)\right]\Sigma_{jm}(\alpha)\end{split}$
(16)

where we see that for $\alpha=0$ Eq. (13) is recovered. It is important to
note that trying to express the above equation as a power series of $(m)^{-1}$
one finds that all the coefficients are zero signaling a decay faster than
every power law. In order to go straight to the result we consider homogeneous
amplitude of noises i.e. $T_{1}=T_{L}=T_{a}\gamma_{a}/(\gamma+\gamma_{a})$ so
that the second term of Eq. (16) vanishes. In this condition the matrix
$\hat{B}$ is proportional to the identity so the system can reach
thermodynamic equilibrium. We then take the Fourier transform
$\tilde{\sigma}_{j\omega}(\alpha)=\int dm\exp(i\omega
m)\sigma^{\text{HHP}}_{jm}(\alpha)$ and study the limit $\omega\ll 1$ ($m\gg
1$):

$\tilde{\sigma}_{j\omega}(\alpha)\propto\int_{0}^{\pi}dz\frac{\delta(\omega-z)\sin(jz)}{\Delta(\alpha)-2\cos(z)}\sim\frac{\sin(j\omega)}{\alpha+\omega^{2}}$
(17)

whose inverse Fourier transform for $m>j$ is proportional to an exponential
with characteristic length $\sqrt{\alpha}$, so we have that
$\sigma^{\text{HHP}}_{jm}(\alpha)\sim\exp(-\sqrt{\alpha}m)$. This last result
is valid for a generic $j\ll L$ so it holds also for particles in the bulk. We
note that $\alpha\to 0$ is a singular limit because the pole of the last term
of the above equation tends to the real axis. Regarding variances that we need
to calculate $\zeta_{jm}$ we can write :

$\sigma^{\text{HHP}}_{mm}(\alpha)=\frac{2\alpha
T_{a}}{\pi}\int_{0}^{\pi}dz\frac{\sin^{2}(mz)}{\Delta(\alpha)-2\cos(z)}=T_{a}\sqrt{\frac{\alpha}{4+\alpha}}+o(m^{-1}),\quad
m\gg 1$ (18)

as we expect, in the HHP the asymptotic temperature is a constant that we
explicitly calculate in the SM (S3). We point out that this variance has two
reasonable limiting cases: for $\alpha=0$ it is $o(m^{-1})$ consistently with
the NHHP while $\lim_{\alpha\to\infty}\sigma^{\text{HHP}}_{mm}(\alpha)=T_{a}$
representing the condition for which the external bath overcomes the
interaction so that the variables are in equilibrium with the thermostats.

From this and by the definition of Eq. 4 we can conclude that spatial
correlations in the HHP follow an exponential decay with a finite
characteristic length scale $\xi$:

$\zeta^{\text{HHP}}_{jm}\sim e^{-m/\xi}\quad m\gg 1,\quad\xi=\alpha^{-1/2}.$
(19)

In the SM (S3) we show that this trend holds also without equal noise
amplitudes so it is not strictly related to the equilibrium condition. We note
that looking at this result in the framework of critical phenomena we would
have a critical point at $\alpha_{c}=0$ and a correlation length that diverges
as $\xi\sim(\alpha-\alpha_{c})^{-\nu}$ with a critical exponent $\nu=1/2$.
This critical point would then coincide with the NHHP. Indeed, in this phase,
the system behave as in a critical regime where spatial correlations exhibit a
power low decay. Nevertheless, we make clear that this is just an analogy and
we don’t interpret our results as a phase transition. Moreover, it is
important to remind that an equivalent equilibrium phase transition governed
by temperature could not occur because we are considering a 1D system. In
equilibrium cases there is actually a transition at zero temperature but it
coincides with a physical state with no dynamics. In other words, the model
described by Eqs. 1 can’t be mapped into an Ising or Heisenberg-like
Hamiltonian system maintaining the same properties. We also note that the same
scaling relation between correlation length and characteristic time of the
bath has also been found in dilute granular systems with an hydrodynamic
approach [44] and in dense active systems [26]. Nevertheless in these two
translational invariant systems the equivalent limit for $\alpha=0$ is
meaningless because in the first case it removes the driving while in the
second one it implies a deterministic constant self propulsion. In Fig. 3 we
show that the exponential to power law crossover and the scaling for $\xi$
derived in the large system limit are clearly visible also for finite size
lattices.

Figure 3: a) Spatial correlation function calculated via Eq. (4). The entries
of $\hat{\sigma}$ are obtained from Eq. (6) with $t=s\gg 1$ and diagonalising
$\hat{A}$. The parameters of the system are: $L=500$, $\gamma=5$,
$\beta=1+\alpha$ (i.e. Toeplitz condition) and $\alpha\in[0.002,5]$. We
observe an exponential decay with a growing correlation length that turns into
a power law when $\alpha=0$. b) Scaling of the correlation length obtained
from an exponential fit of $\zeta^{\text{HHP}}_{1m}$ for different
combinations of parameters. We can see that the relation $\xi=\alpha^{-1/2}$
does not depend on the microscopic details of the system. Quasi-Toeplitz cases
are discussed in the next paragraph. In both panels we used
$T_{1}=T_{a}=0.001$ and $T_{L}=0$.

In order to discuss also the characteristic time scales in the HHP, we note
from Eq. (9) that $\lambda_{j}>\gamma_{a}$ $\forall$ $j$ and so for finite
$\alpha$ and $j/L\ll 1$ we have that:

$\tau^{\text{HHP}}_{j}\sim 1/\gamma_{a}=\tau_{a}.$ (20)

This result is consistent with the fact that being correlated with a finite
fraction of the system implies a finite time to receive the information that
effectively determines the dynamics.

To conclude the comparison between HHP and NHHP, we stress that the difference
between the two phases is originated in the structure of the eigenvalues of
$\hat{A}$. In particular, for both space and time correlations, the crucial
ingredient is that the spectrum of $\hat{A}$ accumulates in $\gamma_{a}$ for
$L\gg 1$ (Eq. (9)). Consequently it accumulates to a finite value in the HHP
and to zero in the NHHP. The crossover between the two phases is then governed
by the limit $\alpha\to 0$ that brings to diverging correlation lengths and
times.

### Beyond the Toeplitz case

Up to now we have considered the special case $\beta=1+\alpha$ for which
$\hat{A}$ is a uniform Toeplitz matrix. Now we want to study the system with a
general viscous constant $\gamma_{b}\neq\gamma+\gamma_{a}$ at the boundaries.
Are the results obtained in the previous paragraphs still valid also in this
more general case? In order to answer this question, we follow a procedure,
systematically explained in [42], to diagonalise quasi-uniform Toeplitz
matrices i.e. matrices that deviates from the Toeplitz form just for few
external borders. It does not give an analytical expression of the eigenvalues
and eigenvectors but assures some constraints on their form and allows to find
their values by numerically solving a set of transcendental equations. In
order to uniform our notation with [42] we note that
$\hat{A}=\gamma(2+\alpha)\hat{I}-\gamma\hat{A}^{\prime}$ where:

$\hat{A}^{\prime}=\begin{pmatrix}x&1&&&\bm{0}\\\ 1&0&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&0&1\\\ \bm{0}&&&1&x\end{pmatrix}$ (21)

and $x=1-\beta+\alpha$ so that for $\beta=1+\alpha$ we recover the Toeplitz
case. Once defined $\lambda^{\prime}_{j}$ ($S^{\prime}_{ij}$) as the
eigenvalues (eigenvectors) of $\hat{A}^{\prime}$, then
$\lambda_{j}=\gamma(2+\alpha)-\gamma\lambda^{\prime}_{j}$ and
$S_{jm}=S^{\prime}_{jm}$. If the eigenvalues are parametrized as
$\lambda^{\prime}_{j}=2\cos(k_{j})$ then we can find them by solving:

$k_{j}=\frac{\pi
j+2\phi(k_{j})}{L+1},\quad\phi(k)=k-\tan^{-1}\left(\frac{\sin(k)}{\cos(k)-x}\right)$
(22)

that determine the allowed values of $k_{j}$. The entries of the eigenvector
matrix $\hat{S}$ can then be directly obtained starting from the numerical
solution of Eq. (22) [42].

Figure 4: a) Spatial correlation function for different quasi-Toeplitz cases
in both HHP and NHHP. We can see that the two phases are stable also for large
values of negative $x$. The entries of $\hat{\sigma}$ are obtained from Eq.
(6) for $t=s\gg 1$ and diagonalising $\hat{A}$. b) Spectra of $\hat{A}$ for
different values of $x$ and $\alpha=2.1$. The spectra always accumulate at the
boundary of the band $[\gamma_{a},4\gamma+\gamma_{a}]$ and out-of-band
eigenvalues can occur only for $|x|>1$. We also note that in the range of
interest for the NHHP ($x\in[-\infty,1]$) the spectra are always positive
assuring the stability of the system. In both panel we used $L=500$,
$\gamma=5$, $T_{1}=T_{a}=0.001$ and $T_{L}=0$.

Once calculated all the $\lambda_{j}$ and the $S_{jm}$ we can use Eq. (7b) in
the stationary case to obtain the covariance matrix and consequently the
correlation functions. In Fig. 4 we show the correlation function for some
quasi-Toeplitz cases for both the HHP and the NHHP finding the same asymptotic
behavior obtained for the Toeplitz one in Fig. 3a. Also the scaling for $\xi$
in the HHP does not change (see Fig. 3b). We note that the difference in terms
of parameters between Toeplitz and quasi-Toeplitz cases is that in the former
we have just one adimensional ratio between viscous constants i.e.
$\alpha=\gamma_{a}/\gamma$ while in the latter we can independently fix
$\beta=\gamma_{b}/\gamma$ and $\alpha$.

Given the form with which eigenvalues are parametrized they can take values
only in the band $\lambda^{\prime}_{j}\in[-2,2]$ and equivalently
$\lambda_{j}\in[\gamma_{a},4\gamma+\gamma_{a}]$. Nevertheless, for absolute
values of $x$ large enough, out-of-band eigenvalues can occur [42]. This fact
would compromise the existence of a stationary state in the NHHP because
$\hat{A}$ would cease to be positive semi-definite. A more refined inspection
of the spectral properties is then needed. Being $\beta>0$ by definition we
are sure that $x\in[-\infty,1)$ in the NHHP. For $L\gg 1$ and $|x|>1$ two out-
of-band eigenvalues $\lambda^{\text{out}}_{1,2}$ emerge converging to a common
value given by $\lambda^{\text{out}}_{1,2}=\gamma(2+\alpha-x-x^{-1})$ that, in
our case, is strictly positive preventing any problem of stability (see Fig.
4b). Moreover, as shown in the same panel, we can see that the spectrum of
$\hat{A}$ always accumulates at the boundary of the band independently from
the value of $x$. This is also clear by taking $j/L\ll 1$ or $\sim 1$ in Eqs.
(22) and verifying that $k_{j}$ tends respectively to $0$ or $\pi$.
Consequently the $\lambda^{\prime}_{j}$s always accumulate in $2$ and the
$\lambda_{j}$s in $\gamma_{a}$. This generalizes our result about the power
law decay in the NHHP (i.e. with $\gamma_{a}=0$) for any $\gamma_{b}>0$
because, as explained in the previous paragraphs, its origin relies in the
accumulation of the $\lambda_{j}$ spectrum in zero (see also Fig. 4).

### Correlations in the bulk and finite size effects

In previous paragraphs we focused on the correlation function with respect to
the first site $\zeta_{1m}$ in the limit $L\gg m\gg 1$. These conditions,
particularly in the NHHP, were crucial ingredients for calculations. Moreover,
in Fig. 3a and Fig. 4a we have always shown the correlation function in the
case of $T_{L}=0$ in order to treat cases more compatible with our
calculations where the terms proportional to $T_{L}\sim\mathcal{O}(1/L)$ are
neglected. In this condition the only source of stochasticity is the bath on
the first site so the finite size effects do not substantially affect the
shape of $\zeta_{1m}$. Thus, the power law regime in the NHHP spans almost all
the system size.

Here we want to discuss the behavior of spatial correlations between particles
in the bulk (i.e. $\zeta_{jm}$ with $1\ll j,m\ll L$) and the finite size
effects for $T_{L}\neq 0$. In Fig. 5 we show $\zeta_{j(j+m-1)}$ with $j=1,L/2$
for different values of $L$ and $\alpha$. In all the cases we have
$T_{1}=T_{a}=T_{L}\neq 0$. The correlation function with respect $L/2$ is
representative for the bulk and we can see from Fig. 5 that in the HHP it
presents an exponential decay with a correlation length independent from $L$
while in the NHHP it decays slower than a power law:
$\zeta^{\text{NHHP}}_{L/2(L/2+m-1)}$ remains essentially constant up to a
sharp cutoff that increases by raising $L$. Regarding
$\zeta^{\text{NHHP}}_{1m}$ for $T_{L}\neq 0$, we can still observe the power
law decay $\sim m^{-2}$ predicted in the previous paragraphs but with a sharp
cutoff that occurs when $m$ is large enough and depending on $L$. In Fig. 5b
we show the same curves as a function of $m/(L/2)$ and we note that the
cutoffs of the correlation functions in the NHHP collapse signaling that their
size scales linearly with $L$. In other words this confirms that, also when
the boundary effects affect the shape of $\zeta_{jm}$, the NHHP presents
_scale-free_ correlations. Indeed the only typical correlation length that one
can define grows with system size. As we expect, the correlation functions in
the HHP separates when plotted as a function of $m/(L/2)$ because their decay
is strictly defined by $\alpha$ regardless of $L$.

Figure 5: a) Spatial correlation function with respect to the site
$j=1,\frac{L}{2}$ for $\beta=2$, $T_{1}=T_{a}=T_{L}=0.001$ and different
values of $\alpha$ and $L$. The entries of $\hat{\sigma}$ are obtained from
from Eq. (6) for $t=s\gg 1$ and diagonalising $\hat{A}$. b) Same curves shown
in the left panel but as a function of the rescaled distance $m/(L/2)$. The
collapse of the cutoffs is a signature of _scale-free_ correlations[13]

## Discussion

We studied spatial and temporal correlations in the NESS reached by a velocity
field with viscous interactions defined on the lattice and coupled with
Brownian baths. The model reproduces three main characteristics of vibrated
granular matter at high density i.e. dissipative forces, permanent contacts
and non-homogeneous energy injection. The typical correlation lengths and
times have a finite characteristic scale when the bulk particles are coupled
to an external bath (HHP regime); however such a scale diverges with the
system size, as in a _scale-free_ scenario, when the thermal bath is removed
from the bulk particles and kept acting on the boundary sites only (NHHP
regime). Solving this model as a diagonalisable multivariate Ornstein-
Uhlenbeck process, we unveiled the role of non-homogeneous heating in the
development of slow and collective dynamics. We conclude that keeping the bath
only at the boundaries allows to have a driven NESS in which the internal
(deterministic) dynamics - and the corresponding propagation of information
and fluctuations - is not hindered by external disturbances. From a
mathematical point of view this is reflected in the spectral properties of the
interaction matrix that accumulates in zero also in the presence of noises at
the boundaries of the lattice. Our findings provide an example of a mechanism
for which power law decays of correlations can occur out of equilibrium,
shedding light on the emergence of collective behavior in dense granular
matter. Further investigations of this model, considering both harmonic and
viscous interactions, are promising steps towards the understanding of more
general non-equilibrium systems such as active matter and biological
assemblies.

## Supplemental Materials: Details of calculations

### S1: Subleading terms in the large system limit

Here we show how performing the large system limit ($L\gg 1$) subleading terms
$\sim 1/L$ occurs. Starting from Eq.(10) we consider the contribution
proportional to $b_{L}^{2}$:

$b_{L}^{2}\Pi^{2}\sum_{lk}\frac{\sin\left(jl\Pi\right)\sin\left(mk\Pi\right)\sin(lL\Pi)\sin(kL\Pi)}{\Delta(\alpha)-\cos\left(k\Pi\right)-\cos\left(l\Pi\right)}$
(23)

where $\Pi=\pi/(L+1)$ and we note that:
$\sin(lL\Pi)\sin(kL\Pi)=\left(-1\right)^{k+l+2}\sin(l\Pi)\sin(k\Pi)$.
Considering a generic function $f$ we can write

$\Pi^{2}\sum_{lk}(-1)^{k+l+2}f(jl\Pi,mk\Pi)=\Pi^{2}\sum_{nh}\big{[}f(2jn\Pi,2mh\Pi)-f(2jn\Pi+j\Pi,2mh\Pi)+f(2jn\Pi+j\Pi,2mh\Pi+m\Pi)\\\
-f(2jn\Pi,2mh\Pi+m\Pi)\big{]}$ (24)

that taking the large system limit $L\gg 1$ and replacing sums with integrals
as $\Pi\sum_{m=0}^{m=L/2}f(2m\Pi)\rightarrow\frac{1}{4}\int_{0}^{\pi}dxf(x)$
becomes:

$\frac{1}{4}\int_{0}^{\pi}dzds\left[f(jz,ms)-f(jz+j\Pi,ms)+f(jz+j\Pi,ms+m\Pi)-f(jz,ms+m\Pi)\right]\sim\mathcal{O}(1/L),\quad
L\gg 1,\quad m\vee j\ll L$ (25)

because all the terms at the zeroth order vanish in the integrand. This
explains why it is possible to neglect the term proportional to $b_{L}^{2}$ in
Eq. (10) once the large system limit is taken and for $j\vee m$ small enough.
This is consistent with the idea that the effect of the bath acting on the Lth
site can be neglected only if $\sigma_{jm}$ is calculated for sites that are
far away from $L$.

### S2: Covariance matrix in the NHHP

Here we give some details about the calculations necessary to derive the
asymptotic predictions of Eqs. (14) from Eq. (13). To do so we start from the
latter equation in a form more suitable for next calculations:

$\sigma^{\text{NHHP}}_{jm}=\lim_{L\to\infty}\frac{4T_{1}}{\pi^{2}}\int_{\frac{\pi}{L+1}}^{\frac{\pi
L}{L+1}}dz\int_{\frac{\pi}{L+1}}^{\frac{\pi
L}{L+1}}ds\sin(jz)\sin(ms)g(z,s)\quad\text{where}\quad
g(z,s)=\frac{\sin(z)\sin(s)}{2-\cos(z)-\cos(s)}.$ (26)

In this expression we have shown the explicit form of the large $L$ limit
because the integrand of the function $g$ is a function of both $z$ and $s$
that is singular in the point $(0,0)$. Indeed, its right value in the origin
comes from the limit for large $L$ of the integration domain
$[\frac{\pi}{L+1},\frac{\pi L}{L+1}]\times[\frac{\pi}{L+1},\frac{\pi L}{L+1}]$
in the $zs$ plane. More specifically we have that $0\leq g(z,s)\leq 1$
$\forall z,s\in[0,\pi]$ and that $\lim_{z\to 0}g(z^{a},z^{b})\sim z^{a-b}$ if
$a\geq b$. In the remainder, we consider the integration intervals as
$[\frac{\pi}{L+1},\pi]$ because the singularity is just in the origin.
Integrating two times by parts and noting that $g(\pi,s)=g(z,\pi)=0$ $\forall$
$z,s$ we have:

$\sigma^{\text{NHHP}}_{jm}=\lim_{L\to\infty}\frac{4T_{1}}{\pi^{2}jm}\Biggl{[}\cos\left(\frac{j\pi}{L+1}\right)\cos\left(\frac{m\pi}{L+1}\right)g\left(\frac{\pi}{L+1},\frac{\pi}{L+1}\right)+\cos\left(\frac{m\pi}{L+1}\right)\int_{\frac{\pi}{L+1}}^{\pi}dz\cos\left(jz\right)\partial_{z}g\left(z,\frac{\pi}{L+1}\right)\\\
+\cos\left(\frac{j\pi}{L+1}\right)\int_{\frac{\pi}{L+1}}^{\pi}ds\cos\left(ms\right)\partial_{s}g\left(\frac{\pi}{L+1},s\right)+\int_{\frac{\pi}{L+1}}^{\pi}dsdz\cos\left(jz\right)\cos\left(ms\right)\partial_{zs}g(z,s)\Biggr{]}.$
(27)

We want to show that $\sigma^{\text{NHHP}}_{jm}\sim(jm)^{-1}$ so we have to
demonstrate that the sum of the terms in the square brackets is
$\mathcal{O}(1)$ for $m,j\gg 1$ in the large $L$ limit. The first term clearly
tends to $1$ when $L\to\infty$ regardless the value of $j$ and $m$ (remember
that $j,m\ll L$). Reintroducing $\Pi=\pi/(L+1)$ we can express Eq. (27) as:

$\sigma^{\text{NHHP}}_{jm}\sim\frac{4T_{1}}{\pi^{2}jm}\left[1+C_{jm}\right]\quad\text{where}\quad
C_{jm}=\lim_{L\to\infty}\left[\cos(m\Pi)I_{j}+\cos(j\Pi)I_{m}+I_{jm}\right]$
(28)

and where $I_{j}$, $I_{m}$ and $I_{jm}$ are respectively the integrals of the
second, third and fourth term in the square brackets of Eq. (27). The estimate
of the asymptotic behavior of such integrals is not trivial because of the
presence of the derivatives of $g(z,s)$ that diverge in the origin. We then
proceed by estimating upper bounds. It is important to note that, in order to
demonstrate $\sigma^{\text{NHHP}}_{jm}\sim(jm)^{-1}$, requiring
$C_{jm}\sim\mathcal{O}(1)$ or $|C_{jm}|\leq 1$ is not enough because it would
bring contributions as $-1\pm o(1/j)$ that imply the emergence of a faster
decay. The right thing to do is instead to show that $|C_{jm}|\leq c$ with
$c<1$. In this way, we could be sure that $C_{jm}$ cannot cancel 1 in Eq.
(28). Starting by $I_{j}$, we define $u(z)=\partial_{z}g(z,\frac{\pi}{L+1})$
and rewrite it as:

$I_{j}=\int_{\frac{\pi}{L+1}}^{\pi+\frac{\pi}{L+1}}dz\cos\left(jz\right)u(z)+\mathcal{O}(1/L)$
(29)

Now we note that the interval of integration is much larger than the period
$T_{j}=\frac{2\pi}{j}$ of the cosine so we can split it in a sum of
contributions over consecutive periods. Without loss of generality we can
assume $j$ even and exploit the periodicity of the cosine obtaining:

$I_{j}=\sum_{k=1}^{k=j/2}\int_{(k-1)T_{j}+\Pi}^{kT_{j}+\Pi}dz\cos(jz)u(z)=\frac{1}{j}\int_{\Pi}^{2\pi+\Pi}dx\cos(x)\sum_{k=1}^{k=j/2}u\left(\frac{x}{j}+(k-1)T_{j}\right)$
(30)

where we have have changed variable as $x=jz+2\pi(k-1)$ and reintroduced the
symbol $\Pi=\frac{\pi}{L+1}$. Now we use the fact that $T_{j}\ll 1$ to
exchange the sum over $k$ with an integral as
$\sum_{k}f\left((k-1)T_{j}\right)\to T_{j}^{-1}\int d\phi_{j}f(\phi_{j})$ and
return to an expression with $g$:

$I_{j}=\frac{1}{2\pi}\int_{\Pi}^{2\pi+\Pi}dx\cos(x)\int_{0}^{\pi-\frac{2\pi}{j}}d\phi_{j}u\left(\frac{x}{j}+\phi_{j}\right)=\frac{1}{2\pi}\int_{\Pi}^{2\pi+\Pi}dx\cos(x)\left[g\left(\frac{x}{j}+\pi-\frac{2\pi}{j},\Pi\right)-g\left(\frac{x}{j},\Pi\right)\right].$
(31)

The function $g$ can be regularly expanded in series around the point
$(\pi,0)$. Doing this, it’s easy to verify that the integral of the first term
in the brackets gives $\mathcal{O}(1/j)$ contributions. We can’t perform such
an estimate for $g(x/j,\Pi)$ because the derivatives near the origin are not
well defined. Nevertheless, we know that $g(x/j,\Pi)\in[0,1]$ $\forall$
$x\in[\Pi,2\pi/+\Pi]$ if $j$ is sufficiently large so we can estimate an upper
bound for $I_{j}$ (and $I_{m}$) as: $\lim_{L\to\infty}|I_{j(m)}|\leq 1/\pi$
for $j\gg 1$. This happens because, given $T$ a $2\pi$-large interval with
$T_{+(-)}$ the sub-interval where the cosine is positive(negative) and
$g(x)\in[0,1]$ if $x\in T$, we can write:

$\Bigg{|}\int_{T}\cos(x)g(x)\Bigg{|}=\Bigg{|}\bigg{|}\int_{T_{+}}\cos(x)g(x)\bigg{|}-\bigg{|}\int_{T_{-}}\cos(x)g(x)\bigg{|}\Bigg{|}\leq\frac{1}{2}\int_{T}\big{|}\cos(x)\big{|}=2$
(32)

With the same kind of calculations leading to Eq. (31) we obtain:

$I_{jm}=\frac{1}{4\pi^{2}}\int_{\Pi}^{2\pi+\Pi}dxdy\cos(x)\cos(y)g\left(\frac{x}{j},\frac{y}{m}\right)+\mathcal{O}((mj)^{-1}).$
(33)

Using inequalities similar to the ones of Eq. (32) but for 2D integrals we
estimate the upper bound of Eq. (33) as $\lim_{L\to\infty}|I_{jm}|\leq
2/\pi^{2}$ for $j,m\gg 1$. Putting together these results in the definition of
$C_{jm}$ of Eq. (28) we are sure that in the large $L$ limit:

$|C_{jm}|\leq\lim_{L\to\infty}\left[|I_{j}|+|I_{m}|+|I_{jm}|\right]=\frac{2}{\pi}\left(1+\frac{1}{\pi}\right)\simeq
0.83926<1\quad\text{for}\quad j,m\gg 1$ (34)

We conclude that $\sigma^{\text{NHHP}}_{jm}\sim(jm)^{-1}$ from which Eq. (14a)
is straightforward.

It is important to note that, in order to obtain Eqs. (30) and (31), we need
both $j$ and $m\gg 1$. So we have to use another way to estimate the
asymptotic behavior of $\sigma^{\text{NHHP}}_{1m}$. It can be rewritten as

$\sigma^{\text{NHHP}}_{1m}=\frac{4T_{1}}{\pi^{2}}\int_{0}^{\pi}dzds\sin(ms)g_{1}(s,z)\quad\text{where}\quad
g_{1}(s,z)=\frac{\sin^{2}(z)\sin(s)}{2-\cos(z)-\cos(s)}$ (35)

and $g_{1}$ is regular in the origin because $\lim_{z\to
0}g_{1}(z^{a},z^{b})=0$ $\forall$ $a,b>0$. We can perform the integral over
$z$ obtaining
$\int_{0}^{\pi}dzg_{1}(z,s)=\pi\left[-2+\cos(s)+\sqrt{6-2\cos(s)}\sin(s/2)\right]\sin(s)$
where the first two terms in the brackets vanish when also the integral over
$s$ is performed ($m$ is an integer). We have now that
$\sigma^{\text{NHHP}}_{1m}=\frac{4T_{1}}{\pi^{2}}\int_{0}^{\pi}ds\sin(ms)f(s)$
where $f(s)=\sin(s)\left[\sqrt{6-2\cos(s)}\sin(s/2)\right]$. Integrating four
times by parts and noting that $f(0)=f(\pi)=f^{\prime\prime}(\pi)=0$ while
$f^{\prime\prime}(0)=2$ we obtain:

$\sigma^{\text{NHHP}}_{1m}=\frac{8T_{1}}{\pi m^{3}}+R_{m}\sim\frac{8T_{1}}{\pi
m^{3}}+\mathcal{O}(m^{-5})\quad m\gg 1$ (36)

where $R_{m}=(m)^{-4}(\pi)^{-1}\int_{0}^{\pi}ds\sin(ms)f^{(4)}(s)$ so
$|R_{m}|\leq(m)^{-4}(\pi)^{-1}|\text{max}(f^{(4)}(s))|\int_{0}^{\pi}ds|\sin(ms)|=2(m)^{-5}(\pi)^{-1}|\text{max}(f^{(4)}(s))|\simeq
19(m)^{-5}(\pi)^{-1}$ . The last quantity needed for the Eqs. (14) is
$\sigma^{\text{NHHP}}_{11}=\int_{0}^{\pi}dzds\sin(z)\sin(s)g(z,s)=\pi^{2}-8\pi/3$
that is finite and does not depends on $m$ so the asymptotic behavior for
$\zeta_{1m}$ directly follows from the ones derived for Eqs. (14a) and (14b).

### S3: Covariance matrix in the HHP

In order to derive Eq. (16) from Eq. (10) we have to discuss the contributions
coming from the sum $\sum_{n}b_{n}^{2}\sin(ln\Pi)\sin(kn\Pi)$ that compares in
the latter. As explained in the first appendix, the term proportional to
$b^{2}_{L}$ gives a subleading term $\mathcal{O}(1/L)$ in the large system
limit while the one proportional to $b_{1}^{2}$ gives
$4T_{1}(1+\alpha)(\pi)^{-2}\Sigma_{jm}(\alpha)$. Regarding the other
contributions, we exploit orthogonality to express the remaining sum as:

$\sum_{n=2}^{n=L-1}\sin(ln\Pi)\sin(kn\Pi)=\frac{L+1}{2}\delta_{kl}-\sin(l\Pi)\sin(k\Pi)-\sin(lL\Pi)\sin(kL\Pi)$
(37)

where again the last term gives $\mathcal{O}(1/L)$ for $L\gg 1$. Thus, using
this equation and neglecting subleading terms, Eq. (10) becomes:

$\sigma_{jm}(\alpha)=\Pi^{2}\sum_{lk}\frac{\sin\left(jl\Pi\right)\sin\left(mk\Pi\right)}{\Delta(\alpha)-\cos\left(k\Pi\right)-\cos\left(l\Pi\right)}\left[\frac{2\alpha
T_{a}(L+1)}{\pi^{2}}\delta_{kl}+\frac{4T_{1}}{\pi^{2}}\left(1+\alpha\left(1-\frac{T_{a}}{T_{1}}\right)\right)\sin(l\Pi)\sin(k\Pi)\right]$
(38)

that in the large system limit gives Eq. (16).

In the main text we proceed from Eq. (16) by considering constant amplitude of
noise i.e. $T_{1}=T_{a}\gamma_{a}/(\gamma+\gamma_{a})$. In this way the term
proportional to $\Sigma(\alpha)$ vanishes and one can shorten calculations
concentrating just on the integral over $z$. To verify that the asymptotic
behavior of Eq. (19) holds also without constant amplitude of noises we have
to show that $\Sigma_{jm}(\alpha)$ does not decay slower than
$\exp({-\sqrt{\alpha}m})$. We then consider the fourier transform
$\tilde{\Sigma}_{j\omega}(\alpha)=\int dm\exp(i\omega m)\Sigma_{jm}(\alpha)$
for small $\omega$:

$\tilde{\Sigma}_{j\omega}\sim\int_{0}^{\pi}dz\frac{\sin(jz)\sin(z)\omega}{1+\alpha-\cos(z)+\frac{\omega^{2}}{2}}\quad\text{so}\quad\Sigma_{jm}\sim\int_{0}^{\pi}dz\frac{\sin(jz)\sin(z)}{1+\alpha-\cos(z)}\exp(-m\sqrt{2(1+\alpha-\cos(z))})$
(39)

and for this last expression is simple to show that
$|\Sigma_{jm}|\leq\frac{\pi}{\alpha}\exp(-\sqrt{2\alpha}m)$. Then we are sure
that its behavior for large $m$ will be subleading with respect to
$\exp(-\sqrt{\alpha}m)$.

To complete the discussion about the exponential decay in the HHP we need to
evaluate the result of Eq. (18). We then write such integral after one
integration by parts obtaining:

$\frac{2\alpha
T_{a}}{\pi}\int_{0}^{\pi}dz\frac{\sin^{2}(mz)}{\Delta(\alpha)-2\cos(z)}=\frac{2\alpha
T_{a}}{\pi}\left[\frac{\pi}{2(4+\alpha)}-\int_{0}^{\pi}dz\frac{z\sin(z)}{(\Delta(\alpha)-2\cos(z))^{2}}-\int_{0}^{\pi}dz\frac{\sin(mz)\sin(z)}{2m(\Delta(\alpha)-2\cos(z))^{2}}\right]$
(40)

from which we have that
$\sigma^{\text{HHP}}_{mm}(\alpha)=T_{a}\sqrt{\frac{\alpha}{4+\alpha}}+o(m^{-1})$.

### S4: Spatial correlation in the cooling state

An important question that often arise in granular systems regards the
relation between the properties of the cooling dynamics and the one of the
NESS obtained with the injection of energy. In our case we obtain the cooling
state by switching off all the temperatures in the lattice (matrix $\hat{B}$
with all zero entries). In this situation the covariance matrix is simply
given by Eq. (7a). Where the brackets $\langle\rangle$ refer to a mean on the
initial condition. Exploiting the symmetricity of $\hat{A}$ we can rewrite it
as:

$\sigma_{jm}(t,s)=\sum_{nhkl}S_{hn}e^{-\lambda_{n}t}S^{+}_{nk}\langle
v_{k}(0)v_{l}(0)\rangle S_{lh}e^{-\lambda_{h}s}S^{+}_{hj}$ (41)

Keeping initial conditions identically and independently distributed around
$0$ with the variance 1 so that $\langle v_{k}(0)v_{l}(0)\rangle=\delta_{kl}$
and exploiting orthogonality of the eigenvectors we have:

$\sigma_{jm}(t,s)=\sum_{n}S_{jn}e^{-\lambda_{n}(t+s)}S^{+}_{nm}$ (42)

That in the Toeplitz case for $t=s$ becomes:

$\sigma_{jm}(t)=\frac{\exp(-2(2\gamma+\gamma_{a})t)\Pi}{\pi}\sum_{n}\sin\left(jn\Pi\right)\sin\left(nm\Pi\right)\exp\left(4\gamma
t\cos\left(n\Pi\right)\right).$ (43)

where we note that for $t=0$ $\sigma_{jm}(0)=\delta_{jm}$ as imposed by the
initial state. The same uncorrelated condition, expected for non-iteracting
systems, is also obtained with $\gamma=0$. Another important properties of the
$\sigma_{ij}(t)$ is that the dependence on $\gamma_{a}$ is factored out from
the sum so, when calculating
$\zeta_{jm}=\sigma_{jm}/\sqrt{\sigma_{jj}\sigma_{mm}}$, it simplifies.
Moreover, also the dependence from $\gamma$ can be removed just by using the
adimensional time $\tilde{t}=\gamma t$. To conclude, during the cooling the
behavior of spatial correlations is crucially different from the one observed
in the two heated phases studied in the main text. In particular, the
parameter $\alpha$ does not play a crucial role as in the NESS. This is an
intriguing result because we found that an external source of energy makes
something more than just keeping alive the dynamics that characterizes the
system when it cools down.

Figure 6: Spatial correlation function in the cooling state after different
times $\tilde{t}$. We observe a collapse by rescaling the horizontal axis by
$\sqrt{\tilde{t}}$.

In Fig. 6 we show $\zeta_{1x}(\tilde{t})$ for different times $\tilde{t}$ and
we clearly observe that it presents a finite cutoff that grows with the delay
time $\tilde{t}$. We can understand it by thinking that the information is
propagating through the system in time. In Fig. 6b we show how rescaling the
space with $\sqrt{\tilde{t}}$ all the curves collapse. So the information
propagates as $\xi(t)\propto\sqrt{\gamma t}$. This result is fully consistent
with diffusion-like coarsening dynamics of vortices, found in other models for
granular velocity fields [45, 29, 46]. In those models however the cooling
state is closer to "dilute" situations where interactions are sequences of
separate binary collisions.

### S5: Reintroduction of space and connection with active matter

Although it is reasonably justified from empirical observations, neglecting
the positional dynamics remains the main approximation of our model. A way to
reintroduce it in our description is to consider a harmonic potential between
nearest neighbors in the lattice. The equation of motion for each particle
would then be of this form

$\displaystyle\dot{x}_{i}=v_{i}$ (44a)
$\displaystyle\dot{v}_{i}=-(\gamma_{a(b)}+2\gamma)v_{i}-2kx_{i}+k(x_{i+1}+x_{i-1})+\gamma(v_{i+1}+v_{i-1})+\sqrt{2T_{a(i)}\gamma_{a(b)}}\xi_{i}(t)$
(44b)

where we consider again a bath on the boundaries characterized by
($\gamma_{b}$, $T_{1(L)}$) and a bath on the bulk ($\gamma_{a}$, $T_{a}$).

It is interesting to note that we can obtain equations of the same form when
considering a 1D chain of (overdamped) active particles with harmonic
interactions, where self-propulsion is modeled using a colored noise $\eta$
(AOUP):

$\displaystyle\dot{x}_{i}=-k(x_{i}-x_{i+1})-k(x_{i}-x_{i-1})+\eta_{i}(t)$
(45a)
$\displaystyle\dot{\eta}_{i}=-\gamma_{a}\eta_{i}+\sqrt{2T_{a}\gamma_{a}}\xi_{i}(t)$
(45b)

where $\xi_{i}$ are Gaussian white noises with unitary variance. Time-deriving
the first of these equations and following standard manipulations, we get
[47]:

$\displaystyle\dot{x}_{i}=v_{i}$ (46a)
$\displaystyle\dot{v}_{i}=-2k\gamma_{a}x_{i}-(\gamma_{a}+2k)v_{i}+k\gamma_{a}(x_{i+1}+x_{i-1})+k(v_{i+1}+v_{i-1})+\sqrt{2T_{a}\gamma_{a}}\xi_{i}(t)$
(46b)

which are formally equivalent to Eqs. (44). If we consider the particles fixed
on the lattice and neglect the positional dynamics we find the analogous of
the granular case studied in the main with a transition in $\gamma_{a}=0$.
While in the granular chain removing the bath on the bulk has a specific and
realistic physical condition (granular materials are often driven only through
boundaries) in the active case it seems meaningless. A self-propelled harmonic
chain modeled by Eqs. (46) has been studied taking account the positional
dynamics and assuming spatially homogeneous self-propulsion [26]. The authors
perform calculations based on translational invariance (they solve the system
in the Bravais reciprocal lattice). This assumption is crucial and it is also
the main difference with our approach in which we are interested in the effect
of non-homogeneous heating. The interesting connection with our investigation
is that they found a correlation length that scales as
$\xi\sim\sqrt{1/\gamma_{a}}$ as in our case [26].

The study of correlations in this kind of 1D systems with both positional
dynamics and non-homogeneous heating is, up to our knowledge, still lacking.
We are currently working in this direction.

## References

  * [1] Narayan, V., Ramaswamy, S. & Menon, N.  Long-lived giant number fluctuations in a swarming granular nematic.  _Science_ 317, 105–108, DOI: 10.1126/science.1140414 (2007).
  * [2] Kumar, N., Soni, H., Ramaswamy, S. & Sood, A.  Flocking at a distance in active granular matter.  _Nat. Comm._ 5, 1–9, DOI: 10.1038/ncomms5688 (2014).
  * [3] Fily, Y. & Marchetti, M. C.  Athermal phase separation of self-propelled particles with no alignment.  _Phys. Rev. Lett._ 108, 235702, DOI: 10.1103/PhysRevLett.108.235702 (2012).
  * [4] Redner, G. S., Hagan, M. F. & Baskaran, A.  Structure and dynamics of a phase-separating active colloidal fluid.  _Phys. Rev. Lett._ 110, 055701, DOI: 10.1103/PhysRevLett.110.055701 (2013).
  * [5] Cates, M. E. & Tailleur, J.  Motility-induced phase separation.  _Annu. Rev. Condens. Matter Phys._ 6, 219–244, DOI: 10.1146/annurev-conmatphys-031214-014710 (2015).
  * [6] Caprini, L., Marini Bettolo Marconi, U. & Puglisi, A.  Spontaneous velocity alignment in motility-induced phase separation.  _Phys. Rev. Lett._ 124, 078001, DOI: 10.1103/PhysRevLett.124.078001 (2020).
  * [7] Alert, R. & Trepat, X.  Physical models of collective cell migration.  _Annual Review of Condensed Matter Physics_ 11, 77–101, DOI: 10.1146/annurev-conmatphys-031218-013516 (2020).
  * [8] Scalliet, C., Gnoli, A., Puglisi, A. & Vulpiani, A.  Cages and anomalous diffusion in vibrated dense granular media.  _Phys. Rev. Lett._ 114, 198001, DOI: 10.1103/PhysRevLett.114.198001 (2015).
  * [9] Plati, A., Baldassarri, A., Gnoli, A., Gradenigo, G. & Puglisi, A.  Dynamical collective memory in fluidized granular materials.  _Phys. Rev. Lett._ 123, 038002, DOI: 10.1103/PhysRevLett.123.038002 (2019).
  * [10] Plati, A. & Puglisi, A.  Slow time scales in a dense vibrofluidized granular material.  _Phys. Rev. E_ 102, 012908, DOI: 10.1103/PhysRevE.102.012908 (2020).
  * [11] Vicsek, T., Czirók, A., Ben-Jacob, E., Cohen, I. & Shochet, O.  Novel type of phase transition in a system of self-driven particles.  _Phys. Rev. Lett._ 75, 1226, DOI: 10.1103/PhysRevLett.75.1226 (1995).
  * [12] Toner, J. & Tu, Y.  Flocks, herds, and schools: A quantitative theory of flocking.  _Phys. Rev. E_ 58, 4828, DOI: 10.1103/PhysRevE.58.4828 (1998).
  * [13] Cavagna, A. _et al._ Scale-free correlations in starling flocks.  _Proceedings of the National Academy of Sciences_ 107, 11865–11870, DOI: 10.1073/pnas.1005766107 (2010).  https://www.pnas.org/content/107/26/11865.full.pdf.
  * [14] Ma, S.-K.  _Modern Theory of Critical Phenomena_ (Routledge, 2018).
  * [15] Cavagna, A. _et al._ Dynamical renormalization group approach to the collective behavior of swarms.  _Phys. Rev. Lett._ 123, 268001, DOI: 10.1103/PhysRevLett.123.268001 (2019).
  * [16] Gradenigo, G., Ferrero, E. E., Bertin, E. & Barrat, J.-L.  Edwards thermodynamics for a driven athermal system with dry friction.  _Phys. Rev. Lett._ 115, 140601, DOI: 10.1103/PhysRevLett.115.140601 (2015).
  * [17] Garrido, P. L., Lebowitz, J. L., Maes, C. & Spohn, H.  Long-range correlations for conservative dynamics.  _Phys. Rev. A_ 42, 1954, DOI: 10.1103/PhysRevA.42.1954 (1990).
  * [18] Grinstein, G., Lee, D.-H. & Sachdev, S.  Conservation laws, anisotropy, and “self-organized criticality”in noisy nonequilibrium systems.  _Phys. Rev. Lett._ 64, 1927, DOI: 10.1103/PhysRevLett.64.1927 (1990).
  * [19] Bertini, L., De Sole, A., Gabrielli, D., Jona-Lasinio, G. & Landim, C.  Macroscopic fluctuation theory.  _Rev. Mod. Phys._ 87, 593, DOI: 10.1103/RevModPhys.87.593 (2015).
  * [20] Rieder, Z., Lebowitz, J. L. & Lieb, E.  Properties of a harmonic crystal in a stationary nonequilibrium state.  _Journal of Mathematical Physics_ 8, 1073–1078, DOI: 10.1063/1.1705319 (1967).  https://doi.org/10.1063/1.1705319.
  * [21] Lepri, S., Livi, R. & Politi, A.  Thermal conduction in classical low-dimensional lattices.  _Physics reports_ 377, 1–80, DOI: 10.1016/S0370-1573(02)00558-6 (2003).
  * [22] Falasco, G., Baiesi, M., Molinaro, L., Conti, L. & Baldovin, F.  Energy repartition for a harmonic chain with local reservoirs.  _Phys. Rev. E_ 92, 022129, DOI: 10.1103/PhysRevE.92.022129 (2015).
  * [23] Derrida, B.  An exactly soluble non-equilibrium system: the asymmetric simple exclusion process.  _Phys. Rep._ 301, 65–83, DOI: 10.1016/S0370-1573(98)00006-4 (1998).
  * [24] Prados, A., Lasanta, A. & Hurtado, P. I.  Large fluctuations in driven dissipative media.  _Phys. Rev. Lett._ 107, 140601, DOI: 10.1103/PhysRevLett.107.140601 (2011).
  * [25] Ishiwata, R., Yaguchi, R. & Sugiyama, Y.  Correlations and responses for a system of $n$ coupled linear oscillators with asymmetric interactions.  _Phys. Rev. E_ 102, 012150, DOI: 10.1103/PhysRevE.102.012150 (2020).
  * [26] Caprini, L. & Marconi, U. M. B.  Time-dependent properties of interacting active matter: Dynamical behavior of one-dimensional systems of self-propelled particles.  _Phys. Rev. Research_ 2, 033518, DOI: 10.1103/PhysRevResearch.2.033518 (2020).
  * [27] Manacorda, A. & Puglisi, A.  Lattice model to derive the fluctuating hydrodynamics of active particles with inertia.  _Phys. Rev. Lett._ 119, 208003, DOI: 10.1103/PhysRevLett.119.208003 (2017).
  * [28] Buttà, P., Flandoli, F., Ottobre, M. & Zegarlinski, B.  A non-linear kinetic model of self-propelled particles with multiple equilibria.  _Kinetic and Related Models_ 12, 791–827, DOI: 10.3934/krm.2019031 (2019).
  * [29] Baldassarri, A., Marini Bettolo Marconi, U. & Puglisi, A.  Cooling of a lattice granular fluid as an ordering process.  _Phys. Rev. E_ 65, 051301, DOI: 10.1103/PhysRevE.65.051301 (2002).
  * [30] Lasanta, A., Manacorda, A., Prados, A. & Puglisi, A.  Fluctuating hydrodynamics and mesoscopic effects of spatial correlations in dissipative systems with conserved momentum.  _New J. Phys._ 17, 083039, DOI: 10.1088/1367-2630/17/8/083039 (2015).
  * [31] Baldassarri, A., Puglisi, A. & Prados, A.  Hydrodynamics of granular particles on a line.  _Phys. Rev. E_ 97, 062905, DOI: 10.1103/PhysRevE.97.062905 (2018).
  * [32] Battle, C. _et al._ Broken detailed balance at mesoscopic scales in active biological systems.  _Science_ 352, 604–607, DOI: 10.1126/science.aac8167 (2016).  https://science.sciencemag.org/content/352/6285/604.full.pdf.
  * [33] Mura, F., Gradziuk, G. & Broedersz, C. P.  Nonequilibrium scaling behavior in driven soft biological assemblies.  _Phys. Rev. Lett._ 121, 038002, DOI: 10.1103/PhysRevLett.121.038002 (2018).
  * [34] Herrmann, H., Hovi, J.-P. & Luding, S.  _Physics of Dry Granular Media_ (Springer, Dordrecht, 1998).
  * [35] Brilliantov, N. V., Spahn, F., Hertzsch, J. M. & Pöschel, T.  Model for collisions in granular gases.  _Physical Review E_ 53, 5382–5392, DOI: 10.1103/PhysRevE.53.5382 (1996).
  * [36] Actually, in many contact models, the dissipative tangential force can swicth from a viscous form to a (non-linear) coulomb one if the normal force between the two grains is small enough. here we are considering cases where dense granular matter is confined by a container and an external field (such as gravity) and we assume that the particles are sufficiently compressed to consider just the linear viscous term.
  * [37] Puglisi, A., Gnoli, A., Gradenigo, G., Sarracino, A. & Villamaina, D.  Structure factors in granular experiments with homogeneous fluidization.  _J. Chem. Phys._ 136, 014704, DOI: 10.1063/1.3673876 (2012).
  * [38] Plata, C., Manacorda, A., Lasanta, A., Puglisi, A. & Prados, A.  Lattice models for granular-like velocity fields: finite-size effects.  _J. Stat. Mech._ 2016, 093203, DOI: 10.1088/1742-5468/2016/09/093203 (2016).
  * [39] Manacorda, A., Plata, C. A., Lasanta, A., Puglisi, A. & Prados, A.  Lattice models for granular-like velocity fields: hydrodynamic description.  _J. Stat. Phys._ 164, 810–841, DOI: 10.1007/s10955-016-1575-z (2016).
  * [40] Puglisi, A.  _Transport and Fluctuations in Granular Fluids: From Boltzmann Equation to Hydrodynamics, Diffusion and Motor Effects_ (Springer, 2014).
  * [41] Gardiner, C.  _Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences_ (Springer-Verlag, Berlin, 1990).
  * [42] Banchi, L. & Vaia, R.  Spectral problem for quasi-uniform nearest-neighbor chains.  _Journal of Mathematical Physics_ 54, 043501, DOI: 10.1063/1.4797477 (2013).
  * [43] Schroeder, K.  Diffusion in crystals with traps: A simple phenomenological model.  _Zeitschrift fur Physik B Condensed Matter_ 25, 91–95, DOI: 10.1007/BF01343313 (1976).
  * [44] Gradenigo, G., Sarracino, A., Villamaina, D. & Puglisi, A.  Fluctuating hydrodynamics and correlation lengths in a driven granular fluid.  _Journal of Statistical Mechanics: Theory and Experiment_ 2011, DOI: 10.1088/1742-5468/2011/08/P08017 (2011).
  * [45] Van Noije, T., Ernst, M., Brito, R. & Orza, J.  Mesoscopic theory of granular fluids.  _Phys. Rev. Lett._ 79, 411 (1997).
  * [46] Baldassarri, A., Puglisi, A. & Sarracino, A.  Coarsening in granular systems.  _C. R. Physique_ 16, 291–302 (2015).
  * [47] Maggi, C., Marconi, U. M. B., Gnan, N. & Di Leonardo, R.  Multidimensional stationary probability distribution for interacting active particles.  _Scientific reports_ 5, 10742 (2015).

## Acknowledgements

The authors are indebted to Marco Baldovin and Lorenzo Caprini for fruitful
scientific discussions and to Alessandro Manacorda for the careful reading of
the manuscript. The authors acknowledge the financial support of Regione Lazio
through the Grant “Progetti Gruppi di Ricerca” N. 85-2017-15257 and from the
MIUR PRIN2017 project 201798CZLJ

## Author contributions statement

A. Pl. conceived the theory. A. Pl. and A. Pu. contributed to the writing of
the manuscript.

## Competing interests

The authors declare no competing financial interests.