Solvable stationary non equilibrium states

We consider the one dimensional boundary driven harmonic model and its continuous version, both introduced in \cite{FGK}. By combining duality and integrability the authors of \cite{FG} obtained the invariant measures in a combinatorial representation. Here we give an integral representation of the invariant measures which turns out to be a convex combination of inhomogeneous product of geometric distributions for the discrete model and a convex combination of inhomogeneous product of exponential distributions for the continuous one. The mean values of the geometric and of the exponential variables are distributed according to the order statistics of i.i.d. uniform random variables on a suitable interval fixed by the boundary sources. The result is obtained solving exactly the stationary condition written in terms of the joint generating function. The method has an interest in itself and can be generalized to study other models. We briefly discuss some applications.


Introduction
Stationary non equilibrium states (SNS) have a rich and complex structure.A natural way to generate a SNS using stochastic interacting particle systems is to put the system, that is evolving on a lattice, in contact with external sources.This is a toy model for a thermodynamic system with external reservoirs.The Markov process obtained with this procedure is typically non-reversible when the reservoirs have different parameters and its invariant measure is the SNS.Due to the non reversibility, such measure is typically difficult to be computed and has long range correlations [3,8].
From a macroscopic point of view, for a few one dimensional solvable models it is possible to get a description of the fluctuations of the SNS by an exact computation of the density large deviations rate functional.This is obtained either by using combinatorial representation of the invariant measure [8] or by the variational dynamic approach of Macroscopic Fluctuation Theory (MFT) [3].Among the solvable models there are the symmetric exclusion process (SEP) and the Kipnis-Marchioro-Presutti (KMP) model [4,20] and more generally all the models having a constant diffusion and a quadratic mobility in the hydrodynamic scaling limit.Due to the presence of long range correlations, the rate functionals are non-local and can be written in terms of the maximization (for SEP) or minimization (for KMP) of an auxiliary function.A problem of interest is the interpretation of the auxiliary function.In the case of the KMP model it has been conjectured in [2] that the auxiliary function can be interpreted as a hidden temperature and the minimization as a contraction principle.This conjecture is solved in [7] where a joint energy-temperature dynamics has been constructed; as a consequence the invariant measure of the boundary driven case is written as a convex combination of inhomogeneous product of exponential distributions whose mean values are distributed according to the invariant measure of an auxiliary opinion model.
From a microscopic point of view, before the most recent developments, to our knowledge, there were essentially a few models with long-range correlations for which the description of the stationary measure was explicit.This is the class of open exclusion type processes, for which it is available a matrix product ansatz (see [10,22]).It was exactly this explicit knowledge that made possible to obtain the density large deviation function by a microscopic computation [11] and then to verify the agreement with the variational structure of MFT [1].Since MFT is believed to have a large degree of universality, as the theory describing fluctuations in diffusive systems, it is therefore important to have additional models of which the SNS is known.Furthermore, for stationary non-equilibrium states a general structure does not exist as it is the case for equilibrium, where one has instead the Boltzmann-Gibbs distribution.
In a series of recent works [15][16][17], two new integrable models have been introduced.These are the family of harmonic models, a class of interacting particle systems, and a suitable continuous version, that can be interpreted as a model for heat conduction.The latter is obtained as a scaling limit of the discrete one.Both families of models are parametrized by a (spin) value s > 0. The integrability of the systems relies on an algebraic description of the generator and the link with integrable systems in quantum spin chains, as is the case also for the class of exclusion processes (see e.g.[22]).Besides sharing the same algebraic description, these two models are also in a duality relation via a moment duality function [15].Both models are of zero range type, i.e., the rate at which particles or energy is transferred from one site to another depends just on the number of particles or amount of energy present on the departure site.However, differently from the classic zero range models, here there are transitions of multiple particles and the boundary driven SNS are not of product type.In Corollary 2.9 of [16] a closed formula 1  of combinatorial nature for the stationary measure has been obtained for the class of harmonic models.The derivation relies on techniques inspired by integrable systems and is based on a direct mapping between non-equilibrium and equilibrium [18,21].A similar study has been done in [15] for the family of integrable heat conduction models, for which moments of the stationary measure have been found via stochastic duality.For both classes of models an explicit description of the longrange correlations has been shown.As shown in [5] all these models have constant diffusion and quadratic mobility and are therefore good candidates for having a mixture of product distribution as invariant measure.
In this paper we provide the probabilistic description of the SNS of these two models from a microscopic perspective.We consider the special case s = 1/2 for the (spin) value.We prove that for this pair of models the invariant measure can be written, like for the KMP model, as a mixture of products of inhomogeneous distributions.Furthermore, for the models considered here, the mixing measure can be explicitly characterized in terms of the order statistics of i.i.d.uniform random 1 When constructing the mixed measures of this paper, we used this formula to check that a product of geometric with mixing measure given by the ordered statistics of i.i.d.uniforms was indeed reproducing the correlation functions in [16].
variables.This probabilistic interpretation sheds light on how the structure of longrange correlations of the SNS is rooted in the correlated structure of the mixing measure.
In the harmonic model of parameter s = 1/2 considered here, at each site of a graph there is a non-negative integer number of particles.When on a vertex x there are η x particles, then k ≤ η x particles can jump across each edge exiting from x with rate 1/k.We consider a one-dimensional lattice with left and right extrema coupled to reservoirs having densities 0 < ρ A ≤ ρ B < +∞.When ρ A = ρ B the model is reversible and its invariant measure is of product type with each marginal being geometric with mean equal to the density of the external reservoirs.When ρ A < ρ B we prove that the invariant measure is a mixture of inhomogeneous product of geometric distributions.The law of the mean values of the inhomogeneous geometric distributions is the order statistics of independent uniform random variables in the interval [ρ A , ρ B ].This is a natural representation, since the computation of the integral over the hidden parameters does not give a transparent expression, being written in terms of hypergeometric functions.For the continuous model we have a similar representation, the heat baths attached at the end points of the bulk have temperatures 0 < T A ≤ T B < +∞ and the geometric distributions have to be substituted by the exponential ones.
Concerning the methodology, our result is proved writing the stationarity condition of the master equation in terms of the joint generating function.This allows a direct verification of the mixed structure of the SNS via a telescopic property.In this paper we apply the method just to two models in order to give a direct and clear presentation.We plan to give a systematic study in the future.We believe the mixed structure with random temperatures/chemical potentials of the stationary measure to be common to several open models of interacting particles, like for example the exclusion process.This is related also to the fact that the corresponding large deviations rate functionals can be written equivalently in terms of infimum (see [2,9,12]).See [14] for results in this direction for the symmetric exclusion process.
Note added: After this article was submitted reference [6] appeared on the arXiv.It contains the mixed measure for the harmonic model with general spin s, which is obtained by a constructive approach that allows to identify the ordered Dirichlet process as the mixing measure.It further contains a direct proof that the measure for s=1/2 derived here coincides with the one derived in [16] (see Appendix A of [6] for the comparison).
2. The discrete Harmonic model with parameter s = 1/2 2.1.The model.We consider a one-dimensional lattice consisting of N sites (the bulk) Λ N := {1, . . ., N } and two ghost lattice sites (the boundaries) ∂Λ N := {0, N + 1} to which we associate two parameters 0 < β A < β B < 1, respectively.On each lattice site we can have an arbitrarily large number of particles and we denote by η x ∈ N 0 the number (possibly zero) of particles at x ∈ Λ N .We consider a continuous-time Markov chain {η(t), t ≥ 0} whose state space is the set Ω N = N ΛN 0 of configurations η = (η 1 , . . ., η N ), with η x ∈ N 0 being the number of particles at site x ∈ Λ N .The stochastic dynamics has a bulk and a boundary part which are described in terms of the generator L N defined below.For any x ∈ Λ N , we denote by δ x ∈ Ω N the configuration defined by δ x (y) = 0 when y = x and δ x (x) = 1.We have (2.1) The bulk generator applied to bounded functions reads: Furthermore, the boundary part which encodes the interaction with the reservoirs is given by: 2.2.Invariant measure.For a generic measure µ on Ω N the stationarity condition µL N = 0 reads as follows: . Our result is the following: Theorem 2.1.The invariant measure of the process with generator (2.1) is given by (2.3) In the above statement we make explicit the dependence of the invariant measure on the parameters ρ A , ρ B , N , while in the rest of the paper we omit such dependence.For simplicity of notation we used the symbol η for a configuration of particles but in order to be compatible with our notation for vectors we remark that in (2.3) η ≡ η should be interpreted as a vector.
In order to better illustrate the result we first give the proof for the case of only one site (Section 3.1) and then generalize it for the case of N sites in Section 3.2.The basic telescoping mechanism is active already in the N = 1 case.

Proof of Theorem 2.1
We introduce the moment generating function of the geometric distribution G m : Like before, given m and λ, we define F m (λ) := N x=1 F mx (λ x ).3.1.The case N = 1.In the case that our lattice is composed by one single node which is in contact with two external reservoirs, the state space of the process Ω 1 is the set of natural numbers.We denote by η 1 ∈ N 0 a generic element of the state space and the generator L 1 (from (2.1) for N = 1) is given by where 0 < β A < β B < 1 are the parameters associated to the two external reservoirs.The stationarity condition for the invariant measure µ is which must be satisfied for all η 1 ∈ N 0 .Theorem 2.1 says that for N = 1 the invariant measure is a mixture of geometric distributions, i.e., Note that in the limit ρ A → ρ B we recover the special equilibrium case, where the invariant measure is just a geometric distribution of mean ρ B .Instead of checking the validity of (3.1) for each η 1 ∈ N 0 , it will be convenient to multiply both sides of (3.1) by λ η 1 and sum over η 1 .In this way, we get an equality between generating functions for each value of λ which is equivalent to the whole set of conditions (3.1).In the sequel we will use the following elementary formulas: We write separately each one of the terms that are obtained by inserting (3.2) into (3.1) and computing the generating function.The first term gives where we used (3.3).Similarly, exchanging the order of summation, the other terms give All in all, by adding the terms, we get that the stationarity condition (3.1) is equivalent to By a direct computation we have the following simple relation for m → F m (λ), the antiderivative of m → F m (λ) : Then, in terms of the antiderivative, (3.4) is rewritten as where F ′ m (λ) denotes the derivative with respect to the parameter m.Performing the integral, apart the common (λ − 1) factor, we get ), which is clearly zero.This concludes the proof of Theorem 2.1 for N = 1.

The general case.
In this section we give the proof of Theorem 2.1 for general N .We now consider the full stationarity condition (2.2) which also contains the bulk terms.With computations similar to the ones done in the previous section we obtain that the stationarity condition (2.2) is equivalent to where we have defined and Using (3.5) the above condition can be also written as where, as usual in this paper, we denote One can check that the integrals are vanishing for each x ∈ {1, 2, . . ., N }.To verify this, let us call O ρA,ρB N −1,x the collection of N − 1 ordered variables m y , with y = x, i.e., where the variable m x is missing; we call m x a generic element of O ρA,ρB N −1,x .Then, by applying Fubini theorem, we get where again λ x is obtained from the vector λ by removing the component λ x .The integral over the variable m x on the right hand side of the above equation can now be performed and we are left with ) which is clearly zero since the term inside the squared parenthesis is identically zero.This concludes the proof of Theorem 2.1.

4.
Integrable heat conduction model with parameter s = 1/2 4.1.The model.In this section we show that the same approach based on a direct computation of the joint generating function holds for a related model.The model was introduced in [17] as a scaling limit of the harmonic model and further generalized in [15].The setting is as in the previous section, namely a one dimensional lattice Λ N with two extra ghost sites representing the reservoirs.Here we denote by z x ∈ R + the arbitrary quantity of energy at site x ∈ Λ N , and by z = (z 1 , . . ., z N ) a generic configuration in Ω N = R ΛN + , i.e., the state space.The generator of the stochastic dynamics is given as the superposition of a bulk part and a boundary part, described below: ) whose action on functions f : Ω N → R that are bounded and Lipschitz is and where we recall that, as in the discrete case, for x ∈ Λ N , δ x is the configuration with δ x (y) = 0 for y = x and δ x (x) = 1.Above T A (respectively, T B ) is the temperature associated to the left (respectively, right) reservoir whose purpose is to destroy the conservation of energy by imposing heat conduction from one side of the chain to the other.When T A = T B = T there is no transport of energy, the model is reversible and its invariant measure is of product type with each marginal being exponential with mean equal to the temperature T of the external reservoirs.Note that since 1/α is not integrable at zero, this is a jump process with a dense set of jumps.We do not address here the delicate issues related to the definition of the process.
4.2.Invariant measure.The stationarity condition imposes that the density µ of the invariant measure satistisfies Our result is the following: Theorem 4.1.The invariant measure of the process with generator (4.1) is given by Here, again, for simplicity of notation we call z a configuration of energies but in order to be compatible with our vector-notation we remark that in (4.3) z ≡ z should be interpreted as a vector.The strategy of the proof is similar to the previous one, namely we consider N = 1 first and then we show the result for a general finite chain of N sites.Below, in order to alleviate the notation for the invariant measure, we drop the dependence on the parameters T A , T B and N .The inner integrals can be computed using "Feynman's trick" which, for a, b > 0, leads to ∞ 0 e −ax − e −bx x dx = log b a , so that we have The key observation regarding F m , the antiderivative of m → F m (t), is the following This allows to write (5.3) as where F ′ m (t) denotes the derivative with respect to the parameter m.As before, by inspection the left hand side of the previous equation is zero and the proof of Theorem 4.1 for N = 1 is concluded.Computing the inner integrals, we get which can be written in terms of the antiderivative F m using equation (5.4) We show that each term of the above sum is zero.To this aim we apply Fubini theorem to the x th term to separate the integral in m x , i.e., which follows from the association of ordered statistics of i.i.d.random variables, as shown in Section 5 of [13].Notice that also for an exponential distribution of mean m, it holds that if m ≤ m ′ then E m E m ′ and so the statement is also true for the continuous model.
We remark that positive correlation inequalities have also been obtained [19] for other models having convex quadratic mobility, such as the symmetric inclusion process and the Brownian energy process.It could be interesting to further investigate whether association is true as well.
Another consequence of the representation of the invariant measure as a mixture of independent random variables is the proof of the large deviation principle for the density profile which can be deduced by a combination of two large deviation principles: one for the order statistics and another for independent inhomogeneous random variables with an additional application of the contraction principle.The heuristic argument is outlined in [2], Section 3.2.A rigorous proof is given in [6], where it is computed the pressure, the density large deviation functional and their additivity principle.This provided a rigorous proof for the expression of the density large deviation function for the whole class of harmonic models obtaining a rate function in accordance with the result of the MFT [3] for systems with convex quadratic mobility and constant diffusion.

1 |O
x F mx−1 (t x ) − 2F mx (t x ) + F mx+1 (t x ) F ′ mx (t x ) (t y ),where O TA ,TB N −1,x has the same meaning as before, namely the collection of N − 1 ordered variables m y with y = x.Computing the inner integral on the right hand side we obtain zero and the proof of Theorem 4.1 is concluded.Since m ≤ m ′ implies G m G m ′ (where the symbol indicates stochastic domination), it follows that g and h are non-decreasing functions.As a consequence, in order to prove (6.1) it is sufficient to show that O ρ A ,ρ B N g(m) • h(m)dm ≥