Non-homogeneous heating induces scale-free correlations and slow time scales in a granular-like velocity field

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 non-homogeneous energy injection for long-range correlations and slow time-scales, proposing an analogy between this simplified dynamical model and recent experiments with 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 ξ 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 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 while in other cases fully non-equilibrium tools have to be developed [16][17][18] .
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 [19][20][21] , non-equilibrium fluctuations 22,23 , correlations and response with non-symmetric couplings 24 , velocity alignment in active matter 25 , systems with Vicsek-like interactions 26,27 , velocity fields in granular materials [28][29][30] . 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 31,32 .
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 γ: The boundary (bulk) sites are coupled with an external bath defined by a drag coefficient γ b (γ 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: Where the first equation holds for 1 < i < L and the η i (t)s are Gaussian white noises with unitary variance: η i (t)η j (t ) = δ i j δ (t − t ). We refer to NHHP when γ a = 0 so that just the first and the Lth sites are heated, while in the HHP we consider a general γ a = 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 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 = 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 = 0 (panel a). This comparison makes clear that -in terms of correlations -the key parameter is γ 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 33,34 . 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 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 γ a rather than T a . d) Autocorrelation times for each site defined as the time τ i for which Γ i (τ i ) = 0.4. The autocorrelation function is defined as 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, γ = 5, γ b = 10, γ a = {3, 0}, T 1 = T L = 0.002, T a = {0.002, 0} after a time t M = 10 8 /γ and with a temporal step dt = 0.05/γ. use in our model 35 . 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. 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 36 .
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 [28][29][30] . 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. The physical difference is that the model considered here is closer to the situation where grains are all in contact and under continuous uninterrupted compression, so that interaction is persistent and not separable into instantaneous collisions. 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 37,38 .
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 continuous interaction variant of a well established family of models previously investigated.
It is important to note that also these dilute models exhibit long-range correlations 37,38 . Nevertheless, those are finite-size effects found in the homogeneous cooling state 39 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 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 (γ b ,T 1(L) ) while the green ones are in contact with the bulk bath, which is switched off in the NHHP.
between the correlations of the cooling state and the NESS ones.

Compact SDE formulation of the model
Defining the vectors V V V = (v 1 , . . . , v L ) , η η η(t) = (η 1 (t), . . . , η L (t)) and the adimensional parameters β = γ b /γ, α = γ a /γ then we can rewrite Eqs. 1 as a multivariate Ornstein-Uhlenbeck process obtaining the following stochastic differential equation (SDE): The information about space-time correlations of the system are encoded in the two times correlation matrixσ (t, s) whose entries are defined as σ jm (t, . We now define the quantity of principal interest in this paper i.e. the static spatial correlation function of the velocity field: With this definition we have ζ jm = 1 if j = m or v j = v m and ζ jm = 0 if v j v m = 0. It is then clear that our goal is to solve Eq. 2 and find the stationary correlation matrixσ = lim t→∞σ (t,t) that exists ifÂ is positive semi-definite. In this conditions, regardless the symmetry ofÂ, the correlation matrix can be found by inverting the relation 40 : Nevertheless, a more direct way to obtain an analytic expression ofσ can be followed exploiting the fact thatÂ is symmetric. In this case there exist a unitary matrixŜ such thatŜŜ + =Î andŜ +ÂŜ =Ŝ +ÂTŜ =λ = diag(λ 1 , λ 2 , . . . , λ L ) whereÎ is the identity matrix, the λ j s are the eigenvalues ofÂ while S jm is the jth component of the ith eigenvector of it. With these hypotheses and in the case ofB = diag(b 1 , . . . , b L ) we can write the covariance matrix in the two-times (with t ≥ s) and non-stationary case: The first matrix represents the transient and the brackets refer to the average over initial conditions while the NESS is described by lim s→∞ 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Ĉ concentrating on the NESS. Definingσ (t ) = lim t→∞σ (t,t + t ) and through Eqs. (6) and (7b) it is also possible to evaluate the single particle 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 τ k = 1/λ k . We will derive σ jm in a specific case where the diagonalisation ofÂ can be done analytically and then follow a numerical technique of diagonalisation 41 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.
In the general case with space-dependent parameters, correlations can be studied diagonalising the matrixÂ or by exploiting Eq. 5 combined with physical constraint onσ . The former strategy, used by us and recently applied in 21,24 , 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 19 . 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 γ a(b) and interaction γ inÂ, 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 γ b = γ + γ a so that β = 1 + α makingÂ a Toeplitz matrix whose eigenvalues and eigenvectors are respectively: where Π = π/(L + 1). Replacing these in Eq. (7b) and taking t = s → ∞, Eq. 6 becomes:

5/16
where ∆(α) = 2 + α. The sums run from 1 to L and: We point out that Eq. (10) is symmetric with respect the center of the lattice (i.e. σ 1m = σ L(L+1−m) ) if the coefficients b n are too.

Power-Law correlations and slow time scales in the NHHP
We first study the NHHP so we put γ a = 0 and use the Toeplitz condition that now reads γ b = γ so β = 1. Exploiting the limit for large systems (L 1), we can exchange sums with integrals as Π ∑ k=L k=1 f (kΠ) → π 0 dz f (z). We note that in Eq. (10), when γ a = 0, the sum over n is actually made of two terms. The one multiplied by γ 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 1 and j, m L (see S1 in the SM). Neglecting it and defining we obtain the covariance matrix for the NHHP: The integral contained in Σ jm (0) is difficult to be explicitly evaluated but the following asymptotic behaviors can be derived in the limit L m 1: As explained in the SM (S2), these results are obtained by expressing σ NHHP jm as a power series of ( jm) −1 by multiple integration by parts and estimating opportune upper bounds. The limit L m 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 σ mm even if it must grow approaching the Lth site if T L = 0. This growth for large m is given by the term proportional to γ 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Â in Eq. (9) for α = 0, we see that, when j/L 1, the slowest time scales in the single particle autocorrelation function behave as: where τ = 1/γ. 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 critical limit in which a finite correlation length diverges. The study of the HHP comes into play to provide an evidence of this last scenario. In the HHP we have γ a = 0 so 6/16 the sum over n of Eq. (10) 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): where we see that for α = 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 find 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 γ a /(γ + γ a ) so that the second term of Eq. (16) vanishes. We then take the Fourier transformσ jω (α) = dm exp(iωm)σ HHP jm (α) and study the limit ω 1 (m 1): whose inverse Fourier transform for m > j is proportional to an exponential with characteristic length √ α, so we have that σ HHP jm (α) ∼ exp(− √ αm). This last result is valid for a generic j L so it holds also for particles in the bulk. We note that α → 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 ζ jm we can write : 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 to two reasonable limiting cases: for α = 0 it is o(m −1 ) consistently with the NHHP while lim α→∞ σ HHP mm (α) = 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 ξ : In the SM (S3) we show that this trend holds also without equal noise amplitudes. We note that looking at this result in the framework of critical phenomena we can say that we have a critical point for α c = 0 and the correlation length diverges as ξ ∼ (α − α c ) −ν with a critical exponent ν = 1/2. Remarkably, this critical point coincides with the NHHP. Consequently, the power law decay found in this phase is explained by the fact that, when heated just on the boundaries, our system is intrinsically in a critical regime. Nevertheless, it is important to remind that there is not an equivalent equilibrium phase transition governed by temperature because we are considering a 1D system, where nucleation of phases can't occur. 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 noise has also been found in dilute granular systems with an hydrodynamic approach 42 and in dense active systems 25 . Nevertheless in these two cases the equivalent limit for α = 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 ξ 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 λ j > γ a ∀ j and so for finite α and j/L 1 we have that: 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Â. In particular, for both space and time correlations, the crucial ingredient is that the spectrum ofÂ accumulates in γ a for L 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 α → 0 that brings to diverging correlation lengths and times. We observe an exponential decay with a growing correlation length that turns into a power law when α = 0. b) Scaling of the correlation length obtained from an exponential fit of ζ HHP 1m for different combinations of parameters. We can see that the relation ξ = α −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.

Beyond the Toeplitz case
Up to now we have considered the special case β = 1 + α for whichÂ is a uniform Toeplitz matrix. Now we want to study the system with a general viscous constant γ b = γ + γ 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 41 , 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 41 we note thatÂ = γ(2 + α)Î − γÂ where: and x = 1 − β + α so that for β = 1 + α we recover the Toeplitz case. Once defined λ j (S i j ) as the eigenvalues (eigenvectors) ofÂ , then λ j = γ(2 + α) − γλ j and S jm = S jm . If the eigenvalues are parametrized as λ j = 2 cos(k j ) then we can find them by solving: that determine the allowed values of k j . The entries of the eigenvector matrixŜ can then be directly obtained starting from the numerical solution of Eq. (22) 41 .
Once calculated all the λ 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 ξ 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. α = γ a /γ while in the latter we can independently fix β = γ b /γ and α.  41 . This fact would compromise the existence of a stationary state in the NHHP becauseÂ would cease to be positive semi-definite. A more refined inspection of the spectral properties is then needed. Being β > 0 by definition we are sure that x ∈ [−∞, 1) in the NHHP. For L 1 and |x| > 1 two out-of-band eigenvalues λ out 1,2 emerge converging to a common value given by λ out 1,2 = γ(2 + α − 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Â always accumulates at the boundary of the band independently from the value of x. This is also clear by taking j/L 1 or ∼ 1 in Eqs. (22) and verifying that k j tends respectively to 0 or π. Consequently the λ j s always accumulate in 2 and the λ j s in γ a . This generalizes our result about the power law decay in the NHHP (i.e. with γ a = 0) for any γ b > 0 because, as explained in the previous paragraphs, its origin relies in the accumulation of the λ 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 ζ 1m in the limit L m 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 ∼ 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 ζ 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. ζ jm with 1 j, m L) and the finite size effects for T L = 0. In Fig. 5 we show ζ j( j+m−1) with j = 1, L/2 for different values of L and α. In all the cases we have T 1 = T a = T L = 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: ζ NHHP L/2(L/2+m−1) remains essentially constant up to a sharp cutoff that increases by raising L. Regarding ζ NHHP 1m for T L = 0, we can still observe the power law decay ∼ 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 ζ 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 α 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. When the bath on the bulk particles is removed (while keeping switched on the boundary bath), the system undergoes a transition between a dynamical state where correlations lengths and times have a characteristic scale (the HHP) to a scale-free scenario where they grow with the lattice size (the NHHP). 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 1) subleading terms ∼ 1/L occurs. Starting from Eq.(10) we consider the contribution proportional to b 2 L : where Π = π/(L + 1) and we note that: sin(lLΠ) sin(kLΠ) = (−1) k+l+2 sin(lΠ) sin(kΠ). Considering a generic function f we can write 10/16 that taking the large system limit L 1 and replacing sums with integrals as Π ∑ m=L/2 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 2 L in Eq. (10) once the large system limit is taken and for j ∨ 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 σ 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: ds sin( jz) sin(ms)g(z, s) where g(z, s) = sin(z) sin(s) 2 − cos(z) − cos(s) .
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 in the zs plane. More specifically we have that 0 ≤ g(z, s) ≤ 1 ∀z, s ∈ [0, π] and that lim z→0 g(z a , z b ) ∼ z a−b if a ≥ b. In the remainder, we consider the integration intervals as [ π L+1 , π] because the singularity is just in the origin. Integrating two times by parts and noting that g(π, s) = g(z, π) = 0 ∀ z, s we have: dsdz cos ( jz) cos (ms) ∂ zs g(z, s) . (27) We want to show that σ NHHP jm ∼ ( jm) −1 so we have to demonstrate that the sum of the terms in the square brackets is O(1) for m, j 1 in the large L limit. The first term clearly tends to 1 when L → ∞ regardless the value of j and m (remember that j, m L). Reintroducing Π = π/(L + 1) we can express Eq. (27) as: 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 σ NHHP jm ∼ ( jm) −1 , requiring C jm ∼ O(1) or |C jm | ≤ 1 is not enough because it would bring contributions as −1 ± o(1/ j) that imply the emergence of a faster decay. The right thing to do is instead to show that |C jm | ≤ 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) = ∂ z g(z, π L+1 ) and rewrite it as: Now we note that the interval of integration is much larger than the period T j = 2π 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: where we have have changed variable as x = jz + 2π(k − 1) and reintroduced the symbol Π = π L+1 . Now we use the fact that T j 1 to exchange the sum over k with an integral as ∑ k f ((k − 1)T j ) → T −1 j dφ j f (φ j ) and return to an expression with g: The function g can be regularly expanded in series around the point (π, 0). Doing this, it's easy to verify that the integral of the first term in the brackets gives O(1/ j) contributions. We can't perform such an estimate for g(x/ j, Π) because the derivatives near the origin are not well defined. Nevertheless, we know that g(x/ j, Π) ∈ [0, 1] ∀ x ∈ [Π, 2π/ + Π] if j is sufficiently large so we can estimate an upper bound for I j (and I m ) as: lim L→∞ |I j(m) | ≤ 1/π for j 1. This happens because, given T a 2π-large interval with T +(−) the sub-interval where the cosine is positive(negative) and g(x) ∈ [0, 1] if x ∈ T , we can write: With the same kind of calculations leading to Eq. (31) we obtain: Using inequalities similar to the ones of Eq. (32) but for 2D integrals we estimate the upper bound of Eq. (33) as lim L→∞ |I jm | ≤ 2/π 2 for j, m 1. Putting together these results in the definition of C jm of Eq. (28) we are sure that in the large L limit: We conclude that σ NHHP jm ∼ ( 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 1. So we have to use another way to estimate the asymptotic behavior of σ NHHP 1m . It can be rewritten as and g 1 is regular in the origin because lim z→0 g 1 (z a , z b ) = 0 ∀ a, b > 0. We can perform the integral over z obtaining π 0 dzg 1 (z, s) = π −2 + cos(s) + 6 − 2 cos(s) sin(s/2) 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 σ NHHP Integrating four times by parts and noting that f (0) = f (π) = f (π) = 0 while f (0) = 2 we obtain: where R m = (m) −4 (π) −1 π 0 ds sin(ms) f (4) (s) so |R m | ≤ (m) −4 (π) −1 |max( f (4) (s))| π 0 ds| sin(ms)| = 2(m) −5 (π) −1 |max( f (4) (s))| 19(m) −5 (π) −1 . The last quantity needed for the Eqs. (14) is σ NHHP 11 = π 0 dzds sin(z) sin(s)g(z, s) = π 2 − 8π/3 that is finite and does not depends on m so the asymptotic behavior for ζ 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 ∑ n b 2 n sin(lnΠ) sin(knΠ) that compares in the latter. As explained in the first appendix, the term proportional to b 2 L gives a subleading term O(1/L) in the large system limit while the one proportional to b 2 1 gives 4T 1 (1 + α)(π) −2 Σ jm (α). Regarding the other contributions, we exploit orthogonality to express the remaining sum as: where again the last term gives O(1/L) for L 1. Thus, using this equation and neglecting subleading terms, Eq. (10) becomes: 2αT a (L + 1) 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 γ a /(γ + γ a ). In this way the term proportional to Σ(α) vanishes and one can shorten calculations concentrating just on the integral over z. To verify that 12/16 the asymptotic behavior of Eq. (19) holds also without constant amplitude of noises we have to show that Σ jm (α) does not decay slower than exp(− √ αm). We then consider the fourier transformΣ jω (α) = dm exp(iωm)Σ jm (α) for small ω: dz sin( jz) sin(z) 1 + α − cos(z) exp(−m 2(1 + α − cos(z))) (39) and for this last expression is simple to show that |Σ jm | ≤ π α exp(− √ 2αm). Then we are sure that its behavior for large m will be subleading with respect to exp(− √ α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: 2αT a π π 0 dz sin 2 (mz) sin(mz) sin(z) 2m(∆(α) − 2 cos(z)) 2 (40) from which we have that σ HHP mm (α) = T a α 4+α + 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 (matrixB with all zero entries). In this situation the covariance matrix is simply given by Eq. (7a). Where the brackets refer to a mean on the initial condition. Exploiting the symmetricity ofÂ we can rewrite it as: Keeping initial conditions identically and independently distributed around 0 with the variance 1 so that v k (0)v l (0) = δ kl and exploiting orthogonality of the eigenvectors we have: σ jm (t, s) = ∑ n S jn e −λ n (t+s) S + nm (42) That in the Toeplitz case for t = s becomes: σ jm (t) = exp(−2(2γ + γ a )t)Π π ∑ n sin ( jnΠ) sin (nmΠ) exp (4γt cos (nΠ)) .
where we note that for t = 0 σ jm (0) = δ jm as imposed by the initial state. The same uncorrelated condition, expected for non-iteracting systems, is also obtained with γ = 0. Another important properties of the σ i j (t) is that the dependence on γ a is factored out from the sum so, when calculating ζ jm = σ jm / √ σ j j σ mm , it simplifies. Moreover, also the dependence from γ can be removed just by using the adimensional timet = γ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 α 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.
In Fig. 6 we show ζ 1x (t) for different timest and we clearly observe that it presents a finite cutoff that grows with the delay timet. 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 √t all the curves collapse. So the information propagates as ξ (t) ∝ √ γt. This result is fully consistent with diffusion-like coarsening dynamics of vortices, found in other models for granular velocity fields 28,43,44 . 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 13/16 Figure 6. Spatial correlation function in the cooling state after different timest. We observe a collapse by rescaling the horizontal axis by where we consider again a bath on the boundaries characterized by (γ b , T 1(L) ) and a bath on the bulk (γ 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 η (AOUP): where ξ i are Gaussian white noises with unitary variance. Time-deriving the first of these equations and following standard manipulations, we get 45 :ẋ 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 γ 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 25 . 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 ξ ∼ 1/γ a as in our case 25 .
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.