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

Andrea Plati, Andrea Puglisi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Andrea Plati, Andrea Puglisi. Long range correlations and slow time scales in a boundary driven granular model. Scientific Reports, 2021, 11 (1), pp.14206. 10.1038/s41598-021-93091-1. hal-04859859

HAL Id: hal-04859859

https://hal.science/hal-04859859v1

Submitted on 31 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

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

Andrea Plati<fnref n="1" />✉ & Andrea Puglisi²,³

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.

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<fnref n="1" />,². Examples, recently observed in experiments and numerical simulations, are motility induced phase transitions in bacteria³⁻⁶, collective migration in epithelial cells⁷, persistent collective rotations in granular systems⁸⁻¹⁰. An important class of CB instances includes flocking and swarming in animals, systematically studied by physicists in the last 25 years¹¹⁻¹³. 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 ξ then we can divide the system in almost independent subsystems of size ~ ξ. 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¹⁴, 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)¹⁵,¹⁶ while in other cases fully non-equilibrium tools have to be developed¹⁷⁻¹⁹.

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 will 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: (1) analytical calculations are often possible, (2) 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²⁰⁻²², non-equilibrium

<footnote n="1">¹Dipartimento di Fisica, Università di Roma Sapienza, P.le Aldo Moro 2, 00185 Rome, Italy. ²Istituto dei Sistemi Complessi-CNR and Dipartimento di Fisica, Università di Roma Sapienza, P.le Aldo Moro 2, 00185 Rome, Italy. ³INFN, University of Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy. ⁴email: andrea.plati@uniroma1.it</footnote>
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-31}$. In the following we will 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 ith 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:

$$ \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) \quad (1a) $$

$$ \dot{v}_1 = -(\gamma + \gamma_b)v_1 + \gamma v_2 + \sqrt{2\gamma_b T_1} \eta_1(t) \quad (1b) $$

$$ \dot{v}_L = -(\gamma + \gamma_b)v_L + \gamma v_{L-1} + \sqrt{2\gamma_b T_L} \eta_L(t) \quad (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') \rangle = \delta_{ij} \delta(t - t')$.

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 Lth 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.

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 Eq. (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 Eq. (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’
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a–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') = \lim_{t \to \infty} \langle v_i(t)v_i(t + t') \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 Eq. (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$.</figcaption></figure>

velocities would be given by Eq. (1). Neglecting the dynamics of positions (they do not appear at all in Eq. 1) is surely the most relevant approximation of our approach: in the SM (S5) we briefly discuss how to go beyond it.

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–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 $V = (v_1, \dots, v_L)$, $\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 Eq. (1) as a multivariate Ornstein-Uhlenbeck process obtaining the following stochastic differential equation (SDE):

$$ \dot{V} = -\hat{A}V + \hat{B}\eta(t) \qquad (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:
<figure><img src="image_5.png" /><figcaption>**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 ith 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.</figcaption></figure>

$$ \hat{A} = \gamma \begin{pmatrix} 1 + \beta & -1 & & & & \mathbf{0} \\ -1 & 2 + \alpha & -1 & & & \\ & \ddots & \ddots & \ddots & & \\ & & -1 & 2 + \alpha & -1 & \\ & & & -1 & 1 + \beta & \end{pmatrix} \quad (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 = \langle [v_j(t) - \langle v_j(t) \rangle] [v_m(s) - \langle v_m(s) \rangle] \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} = \langle v_j v_m \rangle. \quad (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. \quad (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 exists 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 jth component of the ith 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 \ge s$) and non-stationary case:

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

where:

$$ \hat{C}(t,s) = \exp(-\hat{\lambda}t)\hat{S}^+ \langle V(0), V^T(0) \rangle \hat{S} \exp(-\hat{\lambda}s) \quad (7a) $$

$$ G_{jm}(t,s) = \frac{\left(e^{-\lambda_j(t-s)} - e^{-(\lambda_j+\lambda_m)s}\right) \sum_n S_{jm}^+ S_{nm} b_n^2}{\lambda_j + \lambda_m} \quad (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 have 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') = \lim_{t \to \infty} \hat{\sigma}(t+t', t)$ and through Eqs. (6) and (7b) it is also possible to evaluate the single particle autocorrelation function $\Gamma_j(t') \equiv \sigma_{jj}(t')/\sigma_{jj}(0)$:

$$ \Gamma_j(t') = \frac{1}{\sigma_{jj}} \sum_k q_{jk} S_{kj}^+ e^{-\lambda_k t'}, \quad q_{jk} = \sum_{ls} \frac{S_{jl} S_{ks}^+ S_{sl} b_s^2}{\lambda_l + \lambda_k} \quad (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 will 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', t') \rangle = \delta(x-x')\delta(t-t')$. 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$^{22,25}$, 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(jm\Pi) \quad (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(jl\Pi) \sin(mk\Pi) [\sum_n b_n^2 \sin(ln\Pi) \sin(kn\Pi)]}{\Delta(\alpha) - \cos(k\Pi) - \cos(l\Pi)}, \quad (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, & n = 1 \\ 2\gamma_a T_a, & 1 < n < L \\ 2(\gamma + \gamma_a)T_L, & n = L. \end{cases} \quad (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=L}^{\infty} f(k\Pi) \to \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} dz ds \frac{\sin(jz) \sin(ms) \sin(z) \sin(s)}{\Delta(\alpha) - \cos(z) - \cos(s)} \quad (12) $$

we obtain the covariance matrix for the NHHP:

$$ \sigma_{jm}^{\text{NHHP}} = \frac{4T_1}{\pi^2} \Sigma_{jm}(0). \quad (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$:

$$ \sigma_{mm}^{\text{NHHP}} \sim \frac{1}{m^2} \quad (14a) $$

$$ \sigma_{1m}^{\text{NHHP}} \sim \frac{8T_1}{\pi m^3} \quad (14b) $$

$$ \zeta_{1m}^{\text{NHHP}} \sim \frac{1}{m^2} \quad (14c) $$

As explained in the SM (S2), these results are obtained by expressing $\sigma_{jm}^{\text{NHHP}}$ 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 Eqs. (10) to (13).

Equation (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 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 \quad (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):

$$ \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) \quad (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_{jm}^{\text{HHP}}(\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} \quad (17) $$

whose inverse Fourier transform for $m > j$ is proportional to an exponential with characteristic length $\sqrt{\alpha}$, so we have that $\sigma_{jm}^{\text{HHP}}(\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:
<figure><img src="image_6.png" /><figcaption>**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_{1m}^{\text{HHP}}$ 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$.</figcaption></figure>

$$ \sigma_{mm}^{\text{HHP}}(\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 \quad (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_{mm}^{\text{HHP}}(\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_{jm}^{\text{HHP}} \sim e^{-m/\xi} \quad m \gg 1, \quad \xi = \alpha^{-1/2}. \quad (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 behaves as in a critical regime where spatial correlations exhibit a power law decay. Nevertheless, we make clear that this is just an analogy and we do not 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 Eq. (1) cannot 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.

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_j^{\text{HHP}} \sim 1/\gamma_a = \tau_a. \quad (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
<figure><img src="image_7.png" /><figcaption>**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$.</figcaption></figure>

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}'$ where:

$$ \hat{A}' = \begin{pmatrix} x & 1 & 0 \\ 1 & 0 & 1 \\ \vdots & \vdots & \vdots \\ 1 & 0 & 1 \\ 0 & 1 & x \end{pmatrix} \quad (21) $$

and $x = 1 - \beta + \alpha$ so that for $\beta = 1 + \alpha$ we recover the Toeplitz case. Once defined $\lambda'_j(S'_{ij})$ as the eigenvalues (eigenvectors) of $\hat{A}'$, then $\lambda_j = \gamma(2 + \alpha) - \gamma\lambda'_j$ and $S_{jm} = S'_{jm}$. If the eigenvalues are parametrized as $\lambda'_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) \quad (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}$.

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'_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_{1,2}^{\text{out}}$ emerge converging to a common value given by $\lambda_{1,2}^{\text{out}} = \gamma(2 + \alpha - x - x^{-1})$ that, in our case, is strictly positive
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Spatial correlation function with respect to the site $j = 1$, $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 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}$.</figcaption></figure>

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 Eq. (22) and verifying that $k_j$ tends respectively to 0 or $\pi$. Consequently the $\lambda_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 Figs. 3a and 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 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_{L/2(L/2+m-1)}^{\text{NHHP}}$ remains essentially constant up to a sharp cutoff that increases by raising $L$. Regarding $\zeta_{1m}^{\text{NHHP}}$ 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 separate when plotted as a function of $m/(L/2)$ because their decay is strictly defined by $\alpha$ regardless of $L$.

**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.

Received: 23 January 2021; Accepted: 17 June 2021

Published online: 09 July 2021

## References

1. Narayan, V., Ramaswamy, S. & Menon, N. Long-lived giant number fluctuations in a swarming granular nematic. *Science* **317**, 105–108. https://doi.org/10.1126/science.1140414 (2007).

2. Kumar, N., Soni, H., Ramaswamy, S. & Sood, A. Flocking at a distance in active granular matter. *Nat. Commun.* **5**, 1–9. https://doi.org/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. https://doi.org/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. https://doi.org/10.1103/PhysRevLett.110.055701 (2013).

5. Cates, M. E. & Tailleur, J. Motility-induced phase separation. *Annu. Rev. Condens. Matter Phys.* **6**, 219–244. https://doi.org/10.1146/annurev-conmathphys-031214-014710 (2015).

6. Caprini, L., Marconi, U. M., & Puglisi, A., Spontaneous velocity alignment in motility-induced phase separation. *Phys. Rev. Lett.*. https://doi.org/10.1103/PhysRevLett.124.078001 (2020).

7. Alert, R. & Trepat, X. Physical models of collective cell migration. *Annu. Rev. Condensed Matter Phys.* **11**, 77–101. https://doi.org/10.1146/annurev-conmathphys-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. https://doi.org/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. https://doi.org/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. https://doi.org/10.1103/PhysRevE.102.012908 (2020).

11. Vicsek, T., Czíró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. https://doi.org/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. https://doi.org/10.1103/PhysRevE.58.4828 (1998).

13. Cavagna, A. *et al*. Scale-free correlations in starling flocks. *Proc. Natl. Acad. Sci.* **107**, 11865–11870. https://doi.org/10.1073/pnas.1005766107. https://www.pnas.org/content/107/26/11865.full.pdf (2010).

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. https://doi.org/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. https://doi.org/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. https://doi.org/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. https://doi.org/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. https://doi.org/10.1103/RevModPhys.87.593 (2015).

20. Rieder, Z., Lebowitz, J. L. & Lieb, E. Properties of a harmonic crystal in a stationary nonequilibrium state. *J. Math. Phys.* **8**, 1073–1078. https://doi.org/10.1063/1.1705319 (1967).

21. Lepri, S., Livi, R. & Politi, A. Thermal conduction in classical low-dimensional lattices. *Phys. Rep.* **377**, 1–80. https://doi.org/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. https://doi.org/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. https://doi.org/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. https://doi.org/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. https://doi.org/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. Res.* **2**, 033518. https://doi.org/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. https://doi.org/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 Relat. Models* **12**, 791–827. https://doi.org/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*. https://doi.org/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. https://doi.org/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. https://doi.org/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. https://doi.org/10.1126/science.aac8167. https://science.sciencemag.org/content/352/6285/604.full.pdf (2016).

33. Mura, F., Gradziuk, G. & Broedersz, C. P. Nonequilibrium scaling behavior in driven soft biological assemblies. *Phys. Rev. Lett.* **121**, 038002. https://doi.org/10.1103/PhysRevLett.121.038002 (2018).

34. Herrmann, H., Hovi, J.-P. & Luding, S. *Physics of Dry Granular Media* (Springer, 1998).

35. Brilliantov, N. V., Spahn, F., Hertzsch, J. M. & Pöschel, T. Model for collisions in granular gases. *Phys. Rev. E* **53**, 5382–5392. https://doi.org/10.1103/PhysRevE.53.5382 (1996).
36. Actually, in many contact models, the dissipative tangential force can switch 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. https://doi.org/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. https://doi.org/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. https://doi.org/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, 1990).

42. Banchi, L. & Vaia, R. Spectral problem for quasi-uniform nearest-neighbor chains. *J. Math. Phys.* **54**, 043501. https://doi.org/10.1063/1.4797477 (2013).

43. Schroeder, K. Diffusion in crystals with traps: A simple phenomenological model. *Z. Phys. B Condens. Matter* **25**, 91–95. https://doi.org/10.1007/BF01343313 (1976).

44. Gradenigo, G., Sarracino, A., Villamaina, D. & Puglisi, A. Fluctuating hydrodynamics and correlation lengths in a driven granular fluid. *J. Stat. Mech. Theory Exp.* https://doi.org/10.1088/1742-5468/2011/08/P08017 (2011).

## Acknowledgements

The authors are indebted to Marco Baldovin and Lorenzo Caprini for fruitful scientific discussions, Alessandro Manacorda for the careful reading of the manuscript and Sergio Ciuchi for suggesting some intriguing perspectives of this work. 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

A.P. conceived the theory. A.P. and A.P. contributed to the writing of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-93091-1.

**Correspondence** and requests for materials should be addressed to A.P.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021