Edwards thermodynamics for a driven athermal system with dry friction

We obtain, using semi-analytical transfer operator techniques, the Edwards thermodynamics of a one-dimensional model of blocks connected by harmonic springs and subjected to dry friction. The theory is able to reproduce the linear divergence of the correlation length as a function of energy density observed in direct numerical simulations of the model under tapping dynamics. We further characterize analytically this divergence using a Gaussian approximation for the distribution of mechanically stable configurations, and show that it is related to the existence of a peculiar infinite temperature critical point.

We obtain, using semi-analytical transfer operator techniques, the Edwards thermodynamics of a one-dimensional model of blocks connected by harmonic springs and subjected to dry friction. The theory is able to reproduce the linear divergence of the correlation length as a function of energy density observed in direct numerical simulations of the model under tapping dynamics. We further characterize analytically this divergence using a Gaussian approximation for the distribution of mechanically stable configurations, and show that it is related to the existence of a peculiar infinite temperature critical point.
Although systems governed by dissipative interactions do not obey equilibrium statistical mechanics, there have been several attempts to describe such systems with effective equilibrium-like theories. A paradigmatic example is the problem of amorphous packings of frictional grains, for which an effective thermodynamic was proposed by Edwards and coworkers [1][2][3][4][5]. This approach relies on the basic assumption that all mechanically stable packings of grains occupying the same volume have the same probability. This is expected if the system is repeatedly perturbed with "extensive operations" [3], like a shaking of the grains followed by a fast relaxation to a blocked (mechanically stable) configuration. One can then build an effective thermodynamics by determining all mechanically stable configurations (MSCs) of the grains, and computing the mean values of physical observables from flat averages over accessible blocked configurations. The predicted mean values can then be compared to dynamical averages obtained from a given "tapping" protocol which samples blocked configurations.
For athermal systems in which an energy is defined, Edwards' prescription can be formulated as follows [6]. One postulates the existence of an effective temperature T Ed = β −1 Ed such that the probability of a blocked configuration C of energy E(C) takes the form where Z is a generalized partition function; F(C) = 1 if C is a MSC and F(C) = 0 otherwise (only blocked configurations have a nonzero probability). This constraint is non-Hamiltonian, in the sense that it gives a zero probability to (mechanically unstable) configurations having a finite energy, whereas they would have a finite probability in canonical equilibrium. At first sight, Eq. (1) looks like a harmless generalization of equilibrium statistical mechanics, by simply restricting the set of accessible configurations. For instance, introducing an upper bound |x i | < X max on harmonic oscillators x i does not deeply affect their statistical properties. However, the nontrivial point is that the constraints arising from MSC are often much more complex than simply introducing a bound on individual variables, in particular the constraint of mechanical stability may itself introduce strong correlations in the system. Notice that, beyond volume and energy, one may also take into account other quantities when building an Edwards-type thermodynamics [7][8][9][10][11] (e.g., the stress tensor). Several attempts have been done to test the Edwards scenario, not only in packings of grains, experimentally [12][13][14][15] and numerically [16][17][18][19][20], but also, in abstract models like spin and lattice gas [21][22][23][24][25][26], and in glass and spin-glass models [6,[27][28][29][30]. Typically, one uses a specific tapping protocol to sample blocked states, and compares the dynamical average of the observables to the thermodynamic averages obtained from Eq. (1). Although it has been shown explicitly in some cases that Edwards approach is not exact [25,26], it is generally believed to be a reasonably good description in many cases [5]. The main difficulty with the Edwards measure Eq. (1) is that the partition function Z, from which all thermodynamic quantities can be derived, is very complicated to compute due to the complexity of the function F(C) characterizing blocked states [31][32][33][34][35]. Standard approaches are then either to consider abstract models [21][22][23][24][25][26], which are far from any realistic system but simple enough to allow for an explicit solution, or to resort to mean-field [36] or more involved [31] approximations, which still capture part of the interesting phenomenology, but (at least partly) miss relevant information about spatial correlations in the system.
In this Letter, we introduce a realistic model in which Edwards thermodynamics can be computed exactly. We investigate a one-dimensional model of frictional blocks connected by harmonic springs, subjected to a tapping dynamics. Due to the one-dimensional geometry, statistical properties can be computed semi-analytically in the thermodynamic limit using a transfer operator method. Our numerical simulation and theoretical results lead both to an infinite temperature critical point, with a correlation length diverging linearly with the stored energy density -a directly measurable quantities in numerical simulations. We analytically confirm these results using a Gaussian approximation for the joint probability dis-tribution of spring elongations, and further characterize this critical point in terms of the divergence of energy and length fluctuations.
Simulations-Our model is represented by a onedimensional chain of massive blocks connected by N harmonic springs sliding on a horizontal plane [37][38][39][40][41]. Each particle is subjected to dry (Coulomb) friction. The position of the i th -mass is denoted as x i . The effect of dry friction is twofold: when a block is sliding it is subjected to a dissipative force proportional to the dynamic friction coefficient, f i,diss = −µ d mg sgn(ẋ i ); when at rest, in order to make it move one has to apply a force larger than the static friction |f i | > mgµ s . We denote ξ i = x i − x i−1 − l 0 the elongation of the spring connecting block i and block i − 1; l 0 is the rest length of the spring. Taking into account an external force f ext i the total force exerted on a block is . We checked that, over the range of forces explored, no finite size effects appear changing the size from 64 to 2048 blocks. Dropping the tildes, the dimensionless equation of motion reads while the condition to start motion becomes simply |ξ i+1 − ξ i + f ext i | > µ s . We identify the "blocked" configurations as those that in absence of external force are mechanically stable: ∀ i,ẋ i = 0 and |ξ i+1 − ξ i | < µ s . We then define the following tapping dynamics: the external forces f ext i are switched on in Eq.(2) and act during a given period of time τ , after which they are switched off and the system relaxes to a MSC. This procedure, that we call driving cycle, is repeated a large number of times to sample MSCs. At each cycle, the forces f ext i are drawn (randomly for each site i) from a distribution A driving protocol is determined by fixing the parameters ρ and σ. For a given protocol, one can then vary the intensity F and duration τ of the driving. Each MSC is characterized by the typical value of the energy stored by the springs e = (1/2N ) N i=1 ξ 2 i . For each tapping protocol, the average energy e(F, τ ) of the MCSs is found to depend only on the product F τ (see SM).
To characterize the MSC we focus on the correlation function C(|i−j|) = ξ i ξ j between the elongations of the springs at position i and j in the chain. Since this function is trivial (C(|i − j| = C 0 δ ij ) in a thermal harmonic chain at all temperatures, any appearance of correlations is a signature of the unusual statistics associated with the non-Hamiltonian constraints. Correlation functions C(|i − j|) measured for different tapping protocols are shown in Fig. 5. For a given tapping protocol we find that the extent of correlations increases when the average energy of the MSC increases. We extract the correlation length (e) for each case as the distance (i.e., number of springs) |i − j| at which the measured correlation function decays below a conventional threshold C * = 0.2. The insets of Fig. 5a), 5b) and 5c) show the collapse of the correlation function when the x axis is rescaled with (e). Fig. 5d) shows, for all protocols studied, the correlation length growing linearly with the energy (e) ∼ e; the higher the energy, the more the system is correlated. This result may look puzzling at first sight, since commonly the larger is the energy the smaller is the extent of correlations. The key point to understand the physics of our system is the non-Hamiltonian nature of the constraints defining the MSC. If for a certain i we have ξ i 1 (a situation that is typical of a high energy MCS), then, to fulfill the frictional constraint |ξ i+1 − ξ i | < µ s , ξ i+1 must be close to ξ i . The same argument relates ξ i+2 to ξ i+1 , and so on. Therefore, correlations between spring extensions build up in the MSC. We also compare the correlation functions characterizing the blocked states and those at the end of the driving phase (when the force is switched off) and show that the spring-spring correlations are not due to the external driving, but come entirely from the

constraint imposed by static friction (see SM).
It is interesting to notice that the distribution of spring lengths is Gaussian at all energies. This is shown in Fig. 2, both in the case (e) < 1, namely when there is no correlation between springs elongations, and in the case (e) > 1.
Effective theory: transfer operators-Given that MSC are defined in our model by the constraint |ξ i+1 −ξ i | < µ s ∀i, from Eq. (1) the probability of a configuration ξ = (ξ 1 , . . . , ξ N ) reads All the properties of the system can be obtained from the partition sum Z = ∞ −∞ dξ 1 . . . dξ N P (ξ). Using the change of variables ξ i = µ s ξ i , the partition function depends only (up to an irrelevant prefactor) on the product β Ed µ 2 s ; hence, all thermodynamic quantities are functions of T Ed /µ 2 s . We consider periodic boundary conditions for the chain, without imposing any constraint on its total length. For convenience, we fix the rest length to l 0 = ∞, allowing us to take as domain of integration ξ i ∈ (−∞, ∞), while avoiding crossings of masses.
Using Eq.
The operator T has a maximum positive eigenvalue λ max (β Ed , µ s ), which can be computed numerically discretizing the domain of ξ, and using a complete orthonormal basis in L 2 . All relevant thermodynamic observables are computed in the same way (see SM). The free-energy is obtained as f = β −1 Ed ln(λ max (β Ed , µ s )) while the energy reads e = ∂(β Ed f )/∂β Ed = − λ max |∂T /∂β Ed |λ max /λ max . In the following, we compare results from theory and simulations by tuning the temperature T Ed such that the energy e takes the same value as in the numerics.
The behavior of energy as a function of T Ed from the transfer operator approach is shown in Fig. 3. We find two regimes separated by a crossover that depends on µ s : for T Ed µ 2 s there is an "equilibrium-like" regime where e ∼ T Ed while for T Ed µ 2 s one finds e ∼ √ T Ed . The transfer operator approach allows us to compute also the probability distribution p(ξ) of the elongation of a single spring (see SM). The theoretical result for p(ξ) is compared in Fig. 2 with the one estimated numerically from the MCSs, showing good agreement. We find that p(ξ) is Gaussian in all regimes, even when correlations are present.
We also compute the correlation C(|i − j|) = ξ i ξ j , which is very close to an exponential form for all values of T Ed , C(|i − j|) ∝ e −|i−j|/ (β Ed ,µ) (see SM). When T Ed µ 2 s both the correlation length and the energy e grow as √ T Ed (see inset of Fig. 3), implying ∼ e. We thus recover from Edwards thermodynamics the scaling behavior of correlation length with energy observed in the simulated tapping dynamics. This is a remarkable success of the Edwards approach for this system. Conversely, there is almost no correlation between neighboring springs ( < 1), in the "equilibrium-like" low energy regime (e ∼ T Ed ).
We further show in the second inset of Fig. 3 that a di-rect measure of T Ed (within a protocol dependent factor) is obtained from the dissipated energy per tapping cycle and particle, e diss = µ d τ 0 sgn(ẋ i (t))ẋ i (t)dt . Indeed, in the simulations e diss is found to have the same scaling with e as the temperature obtained within the transfer operator approach (i.e., e ∼ e diss if e µ 2 s , e ∼ √ e diss if e µ 2 s ). The dissipated energy can therefore be interpreted as the analog of the thermal energy that allows the system to sample the configuration space [42].
Effective theory: Gaussian ansatz-To obtain approximate analytical expressions for the thermodynamic quantities, we replace the Heaviside function in Eq. (4) by a Gaussian function, yielding Z ∝ Dξ e −S(ξ) with an effective Hamiltonian This effective Hamiltonian corresponds to a positive definite quadratic form S(ξ) = ξ T Aξ/2, where A is a symmetric real Toeplitz matrix (see SM). The matrix A can be exactly diagonalized, yielding analytical expressions for energy, entropy, correlation function and correlation length. The mean energy per particle reads e(T Ed , µ s ) = 1 2 from which we recover that the crossover point between the behaviors e ∼ T Ed and e ∼ √ T Ed is T * Ed ≈ µ 2 s . In Fig. 3, Eq. (9) is compared with the result from the transfer operator Eq. (5), showing a semi-quantitative agreement. We also find that the correlation function is ξ i ξ j ∼ e −|i−j|/ (T Ed ,µs) with (T Ed , µ s ) such that ∼ √ T /µ s in the limit T Ed µ 2 s . This result can be recovered from a field-theoretic viewpoint, by taking a continuous limit in Eq. (7), yielding with a "mass" term m 2 = 2µ 2 s β Ed . The correlation function of such a Gaussian field theory is known [43] to be ξ(x)ξ(y) ∼ e −m|x−y| , so that we recover a correlation length ∼ β −1/2 Ed . This field-theoretic formulation confirms the presence of an infinite temperature critical point, since the mass term goes to zero at infinite temperature. To inspect the critical exponents associated with this critical point, we study the fluctuations of the total energy of the chain, δE = 1 2 N i=1 (ξ 2 i − ξ 2 i ) and the fluctuations of its total length δL = i ξ i . We find (see SM) that the variance of both energy and length diverge linearly with temperature (or, equivalently, as 2 ), Finally, we compute the entropy density s = −∂f /∂T Ed . We find that it saturates at high temperature to a finite value, lim T Ed →∞ s(T Ed , µ s ) = 1 2 ln 2 + ln µ s . This saturation results from the presence of long-range correlation at infinite temperature. This can be confirmed by contrast, computing the "mean-field" entropy density s mf = − dξ p(ξ) ln p(ξ), with p(ξ) the distribution of a single spring elongation ξ. We find that s mf , which discards correlations, diverges like ln T Ed at infinite temperature (see SM), at odds with the saturation of the entropy s.
Conclusions-The present study provides a clearcut example of how an effective thermodynamic theory can successfully describe an athermal dissipative system. The most remarkable difference between standard equilibrium thermodynamics and the effective theory we have presented is the presence of an infinite temperature critical point, with an associated divergence of the correlation length as ∼ T 1/2 Ed (or ∼ e). As seen in the field-theoretic formulation Eq. (9), this infinite temperature critical point results from the long-range correlation generated by static friction in the blocked states. The difference with standard equilibrium systems is that the gradient term in the effective Hamiltonian does not come from an energetic interaction, but from a non-Hamiltonian constraint. Its coefficient is strictly temperature independent, while the coefficient of the energetic term scales inversely with temperature. While temperature-independent terms could also be present at equilibrium (e.g., entropic constraints such as excluded volume), they are usually purely local and do not involve gradient terms. Hence, in spite of its simplicity, our model exhibits a phenomenology clearly distinct from that of equilibrium systems, and the field-theoretic formulation suggests that the results should be quite robust to changes in the details of the model. Future work should investigate this issue in more details.
We acknowledge Financial support from ERC Grant No. ADG20110209. GG thanks A. Cavagna, A. Puglisi and A. Vulpiani for useful discussions and comments.

Details on the simulations
The tapping protocol is defined in Eq. (3) in the main text. The motivation for introducing (annealed) disorder in the tapping protocol is that, when taking ρ = 1 and σ = 0 (i.e., without any disorder in the driving), we found undamped waves which travel across the chain causing an undesired and artificial collapse of the "ensemble" of visited MSC to a very small and specific subset. As mentioned in the main text, for the numerical simulations of the tapping dynamics we consider open boundary conditions. This choice is done in order to avoid constraints on springs which may induce trivial correlations. For instance if one fixes N i=1 ξ i = N l 0 , then also has N 2 l 2 0 = L 2 = ij ξ i ξ j , which is a constraint on the two point correlation function. With open boundary conditions the first and the last blocks are connected to a single spring and their equations of motions are respectivelyẍ Note that in order to have N springs with open boundary conditions we need N + 1 blocks.

Dependence of the energy on driving parameters
The average energy of MSC depends in principle on the four parameters which charaterize the driving cycle: the duration τ of forcing, the intensity of the mean applied force F , its standard deviation σ as well as the fraction ρ of pulled particles. The parameters ρ and σ define the tapping protocol, while F and τ characterize the intensity of the driving. We have found that for a fixed tapping protocol, the energy e of the dynamically reached MSC does not depends separately on F and τ , but is uniquely determined by the product F τ . This is clearly shown by the set of data reported in Fig. 4.

One-cycle experiment: test for the origin of friction induced correlations
In the manuscript we presented data on the "springspring" correlation function ξ i ξ j obtained in the mechanically stable configurations (MSCs). Due to the nonrandom type of the external force F used to sample different MSCs the reader is allowed to wonder if the correlations we measure in the MSCs are really induced by friction at 100%, or a part of them is induced by the same F . The numerical evidence that it is not the case comes from the "one-cycle driving experiment" we are going to discuss. The one-cycle driving experiment is analogous to the tapping dynamics discussed in the main text apart from one thing: at the beginning of each driving cycle the positions of blocks are randomized. More precisely all the x i are reset to x i = i l 0 + dx, where dx is normally distributed and l 0 is the rest length of springs. During the elementary driving cycle of the dynamics a fraction ρ of the particles is pulled for a durantion τ with a constant force F , the force is then switched off and the system relaxes to a MSC. We average over many iterations two correlation functions: a) the forced correlation ξ i ξ j F , which is measured just before switching-off the force at τ ; b) The correlation in the mechanically stable configurations reached after the quench, ξ i ξ j MSC . Our aim in studying ξ i ξ j F is to see how much of the correlation is provided by the external force. On the other hand we must also be sure that when measuring ξ i ξ j F there is no memory of the correlation in ξ i ξ j MSC at the previous step. It is for this reason that at the beginning of each driving cycle the positions of block are randomized. In this way we make ourselves sure that if some correlation is measured in ξ i ξ j F , it comes solely from the external force. The result is illustrated in Fig. 5. The continuous line is ξ i ξ j MSC while the (red) circles represent ξ i ξ j F : being the latter identically zero we can conclude that no correlation is induced by the force and that the correlations we observe in the system come solely from the constraint imposed by static friction on the MSCs.

Transfer operators
The core of the exact solution of the effective thermodynamics of our model is represented by the spectral properties of the linear operator T , acting in the space of square integrable functions L 2 , T : L 2 → L 2 . The operator T is defined as follows: with T (x, y) = e − β Ed As can be easily checked dxdy|T (x, y)| 2 < ∞, so that T is a Hilbert-Schmidt operator, and hence compact and bounded, which guarantees the existence of a maximum positive eigenvalue. The kernel is symmetric and real, T (x, y) = T * (y, x), so that the operator is also selfadjoint, which guarantees the spectral theorem to hold: a complete set of eigenvalues and eigenvectors exist for our operator [1].

Correlation function
By means of the transfer operator approach we can compute several quantities of interest: the one that is most relevant for our study is the correlation between the elongation of far apart springs. Let us consider for instance the elongations ξ i and ξ j , with j > i. In terms of our effective theory the correlation between these two elongations can be written as: which in terms of operators becomes where D is the diagonal operator defined as: Exploiting the fact that the operator T has a complete spectrum of orthonormal eigenvectors we can use the completeness relation 1 = a |a a|, where |a is the eigenvector of T relative to the eigenvalue a, to simplify Eq. (5). One can write the trace in Eq. (5) as Tr(T ) = a a|T |a and then insert between each pair of neighboring operators the completeness 1 = a |a a|. The result is which in the thermodynamic limit becomes lim N →∞ where f b (x) is the eigenfunction related to the eigenvalue b. The correlation functions obtained at different values of the Edwards temperature T Ed by means of a numerical determination of eigenvalues and eigenvectors of T can be found in Fig. 6: at all temperatures, we numerically observe an exponential decay, meaning that the second largest eigenvalue dominates the sum in Eq. (8).

Energy
The average energy of the harmonic chain with N particles is defined as: (9) where in the last line of Eq. (9) we passed the derivative sign ∂/∂β Ed under the trace sign. The "matrix elements" of the operator ∂T /∂β Ed are simply defined as the derivative with respect to β Ed of the matrix elements of T : Exploiting again the definition of the trace Tr[T N ] = a a|T N |a and by inserting the completeness relation 1 = a |a a|, we can write which in the thermodynamic limit yields where f λmax (y) is the eigenfunction related to the maximum eigenvalue.

Marginal distributions
By means of the transfer matrix technique is easy to obtain marginal distributions. For instance we are interested in the probability distribution of the springs elongation, which reads where, again, f λmax (x) is the eigenfunction related to the maximum eigenvalue. This distribution, which has a Gaussian shape, is plotted for different values of T Ed in Fig.2 in the Letter.

Gaussian approximation
An explicit computation of how the internal energy depends on the temperature is possible if one assumes a Gaussian smoothening of the Heaviside function enforcing the static friction constraint: With this approximation the whole partition function reads where ξ represents the vector ξ = (ξ 1 , . . . , ξ N ) and A is a X × N symmetric real Toeplitz matrix of the kind For this kind of matrix the eigenvalues can be explicilty written as Free energy, energy and entropy The free-energy reads then as In the continuous limit the free-energy can be exactly computed and reads:  (20) Note that this method to compute the free-energy is different from the transfer operator method. Here we do not write the partition function as the trace of the N th power of a transfer operator, but we diagonalize exactly the quadratic form appearing in Eq. (15). The partition function is then exactly computed for all N using the Gaussian integral formula, and all eigenvalues of A contribute to the free-energy (while only the largest eigenvalue of the transfer operator contributes to the free energy in the thermodynamic limit).
The energy per degree of freedom can also be computed: Combining Eqs. (20) and (21), the entropy per degree of freedom can be simply obtained as s(β Ed ) = β Ed e(β Ed )− β Ed f (β Ed ). One can check that s(β Ed ) → 1 2 ln 2 + ln µ s when T Ed → ∞. Also, energy fluctuations are given by and at high temperature they diverge as [δE] 2 /N ∼ T Ed , consistently with Eq. (10) in the main text.

Correlation function
When the probability of a configuration of the system is a multivariate gaussian of the kind we know that the value of the average correlation between ξ i and ξ j is For the problem we are studying A is a tridiagonal Toeplitz matrix, for which the explicit formulae for the elements of the inverse is [1]: where d = A 0 /(2|A 1 |), and U n (x) is a Chebyshev polynomial of the second kind of degree n in the variable x.
For such polynomials, a useful identity is known: U n (x) = sin θ(n + 1) sin θ , where θ is the complex number such that x = cos(θ) and 0 ≤ (θ) < π. In our case we have where the last inequality on the right is true for every finite value of the temperature and of the friction coefficient. Therefore the number θ must be a purely imaginary one, i.e., θ = iα, with α > 0. For our problem we have therefore that the Chebyshev polynomials are defined according to the two following identities: U n (1 + µ 2 s β Ed ) = sinh α(n + 1) sinh α 1 + µ 2 s β Ed = cosh(α) Because our problem is in the continuum we are interested in the limit n → ∞ of the expression in Eq. (25). Let us notice the following: where the correlation length in the last row is defined as α = −1 . In the limit where µ 2 s β Ed 1, namely when the temperature is large compared to the friction coefficient, the correlation length must be large and we can estimate its asymptotic scaling: cosh( −1 ) = µ 2 s β Ed + 1 µ 2 s β Ed 1 =⇒ µ 2 s β Ed + 1 ∼ 1 + By recalling that A 1 = −1/(2µ 2 s ), we can finally write the two-point correlation function as: where the correlation length is (µ s , T Ed ) = 1 arcosh(1 + µ 2 s /T Ed ) .
which scales for T Ed → ∞ as (µ s , T Ed ) ∼ √ T Ed /µ s . The variance of the fluctuations of the total length L = N i=1 read with C(r) defined by C(|i − j|) = ξ i ξ j . Using Eq. (31), we get [δL] 2 ∼ N e ∞ 0 dr exp(−r/ ) ∼ N e ∼ N T Ed (34) thus recovering the second equality in Eq. (10) of the main text.

Spring length distribution
The probability P (ξ) is a multivariate Gaussian, so that the marginal distribution p(ξ) = dξ 2 . . . dξ N P (ξ) is also a Gaussian: with The "mean-field" entropy s mf introduced in the text can be computed as which yields the asymptotic result The fact that the entropy of Eq. 37 in the infinite temperature limit diverges, and is then wrong compared to Eq.10 in the Letter, is consistent with the fact that formula in Eq. 37 is approximation which does not take into account the correlations.