A mass transport model with a simple non-factorized steady-state distribution

TL;DR

This paper introduces a mass transport model on a ring with a simple non-factorized steady-state distribution, analyzing its properties and measures of distance from equilibrium.

Contribution

It presents a new mass transport model with a non-factorized steady-state, explicitly derives the non-equilibrium free energy, and compares various non-equilibrium measures.

Findings

Steady-state distribution is a sum of two inhomogeneous product measures.

A symmetric redistribution yields a factorized, equilibrium measure.

Explicit expressions for non-equilibrium free energy and related parameters.

Abstract

We study a mass transport model on a ring with parallel update, where a continuous mass is randomly redistributed along distinct links of the lattice, choosing at random one of the two partitions at each time step. The redistribution process on a given link depends on the masses on both sites, at variance with the Zero Range Process and its continuous mass generalizations. We show that the steady-state distribution takes a simple non-factorized form that can be written as a sum of two inhomogeneous product measures. A factorized measure is recovered for a symmetric mass redistribution, corresponding to an equilibrium process. A non-equilibrium free energy can be explicitly defined from the partition function. We evaluate different characterizations of the `distance' to equilibrium, either dynamic or static: the mass flux, the entropy production rate, the Gibbs free-energy difference between the equilibrium and non-equilibrium stationary states, and the derivative of the non-equilibrium free energy with respect to the applied driving force. The connection between these different non-equilibrium parameters is discussed.

Equations

$$\phi (m_i^{\prime }|S_i)=\frac{v(m_i^{\prime })w(S_i-m_i^{\prime })}{v\ast w(S_i)},S_i\equiv m_i+m_{i+1}$$

$$v\ast w(S)={\int }_0^S\text{rmd}mv(m)w(S-m).$$

$$P({\mathbf{m}}^{\prime },t)=\int \text{rmd}\mathbf{m}T({\mathbf{m}}^{\prime }|\mathbf{m})P(\mathbf{m},t)$$

$$\int \text{rmd}{\mathbf{m}}^{\prime }T({\mathbf{m}}^{\prime }|\mathbf{m})=1.$$

$$T({\mathbf{m}}^{\prime }|\mathbf{m})=\frac{1}{2}{T}_1({\mathbf{m}}^{\prime }|\mathbf{m})+\frac{1}{2}{T}_2({\mathbf{m}}^{\prime }|\mathbf{m})$$

$$T_1({\mathbf{m}}^{\prime }|\mathbf{m})$$

$$T_2({\mathbf{m}}^{\prime }|\mathbf{m})$$

$$P(\mathbf{m})=\frac{1}{Z_N(M)}\left(\prod _{k=1}^{N^{\prime }}v(m_{2k})w(m_{2k+1})+\prod _{k=1}^{N^{\prime }}w(m_{2k})v(m_{2k+1})\right)\delta \left(\sum _{i=1}^N{m}_i-M\right)$$

$$P_j(\mathbf{m})=\frac{2}{Z_N(M)}{Q}_j(\mathbf{m})\delta \left(\sum _{i=1}^N{m}_i-M\right),$$

$$Q_1(\mathbf{m})=\prod _{k=1}^{N^{\prime }}v(m_{2k})w(m_{2k+1}),Q_2(\mathbf{m})=\prod _{k=1}^{N^{\prime }}w(m_{2k})v(m_{2k+1}).$$

$$Q_1({\mathbf{m}}^{\prime })+Q_2({\mathbf{m}}^{\prime })$$

$$\displaystyle\!\!=\frac{1}{2}\int\rmd\mathbf{m}\,\big{[}T_{1}(\mathbf{m}^{\prime}|\mathbf{m})+T_{2}(\mathbf{m}^{\prime}|\mathbf{m})\big{]}\,\big{[}Q_{1}(\mathbf{m})+Q_{2}(\mathbf{m})\big{]}$$

$$\int \text{rmd}\mathbf{m}{T}_k({\mathbf{m}}^{\prime }|\mathbf{m})Q_j({\mathbf{m}}^{\prime })=Q_k(\mathbf{m}).$$

$$v(m)=e^{-\beta \epsilon (m)-\beta h(m)},w(m)=e^{-\beta \epsilon (m)+\beta h(m)}$$

$$e^{-\beta \epsilon (m)}=\sqrt{v(m)w(m)},e^{-\beta h(m)}=\sqrt{\frac{v(m)}{w(m)}}.$$

$$P(\mathbf{m})=\frac{2}{Z_N(M)}{e}^{-\beta {\sum }_{i=1}^N\epsilon (m_i)}\delta \left(\sum _{i=1}^N{m}_i-M\right).$$

$$P(\mathbf{m})=\frac{2}{Z_N(M)}{e}^{-\beta {\sum }_{i=1}^N\epsilon (m_i)}\cosh \left(\sum _{i=1}^N(-1{)}^i\beta h(m_i)\right)\delta \left(\sum _{i=1}^N{m}_i-M\right).$$

$$P(\mathbf{m})=\frac{2}{Z_N(M)}{e}^{-\beta E(\mathbf{m})}\cosh [\beta H(\mathbf{m})]\delta \left(\sum _{i=1}^N{m}_i-M\right)$$

$$E(\mathbf{m})=\sum _{i=1}^N\epsilon (m_i),H(\mathbf{m})=\sum _{i=1}^N(-1{)}^{i}h(m_i).$$

$$G_{j}\binrel@{=}\binrel@@{\mathop{\kern 0.0pt=}\limits_{f\rho\,\ll\,1}}\frac{\big{(}2+(-1)^{j}\big{)}}{\big{(}\varepsilon_{0}-\mu(\rho)\big{)}^{2}}\,\rho^{2}f^{2}+\mathcal{O}\left((\rho f)^{4}\right).$$

$$Z_N(M)=\int \text{rmd}\mathbf{m}[Q_1(\mathbf{m})+Q_2(\mathbf{m})]\delta \left(\sum _{i=1}^N{m}_i-M\right)$$

$$I(\rho )=-\underset{N\to \infty }{\lim }\frac{1}{N}\ln Z(N\rho ),$$

$$Z_N(N\rho )\sim e^{-NI(\rho )}.$$

$$\delta (s)=\frac{1}{2\pi i}{\int }_{a-i\infty }^{a+i\infty }\text{rmd}\zeta e^{\zeta s}$$

$$Z_N(M)=\frac{1}{2\pi i}{\int }_{a-i\infty }^{a+i\infty }\text{rmd}\zeta e^{-\zeta M}{\left[\hat{v}(\zeta )\hat{w}(\zeta )\right]}^{N/2}$$

$$\hat{v}(\zeta )={\int }_0^{\infty }\text{rmd}m{e}^{\zeta m}v(m),\hat{w}(\zeta )={\int }_0^{\infty }\text{rmd}m{e}^{\zeta m}w(m).$$

$$Z_N(N\rho )=\frac{1}{2\pi i}{\int }_{a-i\infty }^{a+i\infty }\text{rmd}\zeta e^{N(\lambda (\zeta )-\rho \zeta )}$$

$$\lambda (\zeta )=\frac{1}{2}\ln [\hat{v}(\zeta )\hat{w}(\zeta )]$$

$$\frac{d\lambda }{d\zeta }({\zeta }^{\ast })=\rho ,$$

$$Z_N(N\rho )\sim e^{-N[\rho {\zeta }^{\ast }(\rho )-\lambda ({\zeta }^{\ast }(\rho ))]}.$$

$$I(\rho )=\rho {\zeta }^{\ast }(\rho )-\lambda ({\zeta }^{\ast }(\rho )).$$

Full text

1 Introduction

One of the goals of non-equilibrium statistical physics is to be able to describe the statistical properties of systems driven in a non-equilibrium steady state by an external non-conservative force. As no general statistical formalism is available to deal with driven systems, exactly solvable models have played an important role in the development of this field. A paradigmatic exactly solvable model is the Asymmetric Simple Exclusion Process (ASEP) [1], either with periodic [2] or open boundary conditions [3, 4, 5, 6]. Generalizations with several types of particles have also been proposed, with periodic [7, 8, 9] or open geometries [10, 11, 12, 13]. The ABC model [14], which includes three types of particles, also falls into this class. The solution of the ASEP model requires in most cases the use of matrix product states [15], often with infinite size matrices, making its analysis relatively involved. Such matrix product state solutions are required even with a periodic geometry, when the model includes several types of particles [7, 8, 9] —except if some restrictive conditions are imposed [16].

Simpler models, like the Zero-Range Process (ZRP) [1, 17, 18] and related mass transport models [19, 20, 21, 22], have also been considered, often in relation to condensation transitions [23, 24, 17]. Multispecies generalizations of these models have also been proposed [23, 25]. When the transition rates satisfy certain conditions [19, 20, 21], these models have the advantage that their steady-state distribution factorizes, making their analytical study much easier. However, in a closed geometry, they have the drawback that the distribution does not depend on the driving force, and thus remains identical to the equilibrium distribution obtained for unbiased dynamics. Note that the same property also holds for the (single-species) ASEP on a ring [2].

In this paper, we propose a class of mass transport models for which the steady-state distribution takes a simple form (a sum of two inhomogeneous product measures) and explicitly depends on the local driving force. The present model is inspired by the equilibrium model considered in [26], though it differs from the latter in several respects, notably the presence of a driving force and of a synchronous dynamics. The simple form of the steady-state probability distribution makes calculations easy, as illustrated below on several examples including the evaluation of a non-equilibrium free energy. In addition, the dependence of the distribution on the forcing allows us to compare dynamical characterizations of the ‘degree of non-equilibrium’ (mass flux and entropy production rate) with static characterizations like the difference of Gibbs free energy functional (or Kullback-Leibler divergence [27]) between the non-equilibrium distribution and the corresponding equilibrium one. We also evaluate the non-equilibrium order parameter introduced by Sasa and Tasaki [28], defined as a derivative of the non-equilibrium free energy with respect to the driving force, and discuss the relationship between these different measures of the ’distance’ to equilibrium.

2 Definition of the model

We consider a one-dimensional lattice with $N$ sites, labelled by $i=1,\dots,N$, with periodic boundary conditions ($i\pm N\equiv i$); $N$ is assumed to be even, namely $N=2N^{\prime}$ with $N^{\prime}$ integer. On each site $i$, one defines a real positive mass $m_{i}$. The model is endowed with a synchronous, discrete time dynamics111Note that although the asynchronous, continuous time dynamics is most often used in this context, synchronous dynamics has also been used in the ASEP [29, 30, 31] and in mass transport models [19, 25]. The dynamics proceeds, at each time step $t=0,1,2,\dots$, by parallel redistributions of mass between neighboring sites $i$ and $i+1$ on one of the two partitions $\mathcal{P}_{1}=\{(2k,2k+1)\}$ and $\mathcal{P}_{2}=\{(2k+1,2k+2)\}$, randomly chosen with equal probability. Once a partition $\mathcal{P}_{j}$ has been selected, all links belonging to the partition $\mathcal{P}_{j}$ are simultaneously updated. To update a link $(i,i+1)$, a new value $m_{i}^{\prime}$ of the mass on site $i$ is randomly drawn from the distribution

[TABLE]

where $v(m)$ and $w(m)$ are arbitrary positive functions, and $v*w(S)$ is the convolution product of $v$ and $w$,

[TABLE]

From mass conservation, the mass on site $i+1$ is, after redistribution, $m_{i+1}^{\prime}=S_{i}-m_{i}^{\prime}$.

3 Master equation and steady-state solution

3.1 Discrete time master equation

To describe the statistical evolution of the system under the above dynamics, we write down the corresponding master equation. A configuration of the system is given by the ordered list $\mathbf{m}=(m_{1},\dots,m_{N})$ of all the masses in the system. The probability density $P(\mathbf{m},t)$ evolves according to the discrete time master equation

[TABLE]

with $\rmd\mathbf{m}=\prod_{i=1}^{N}\rmd m_{i}$, and where $T(\mathbf{m}^{\prime}|\mathbf{m})$ is the probability (density) to jump from configuration $\mathbf{m}$ to configuration $\mathbf{m}^{\prime}$ in a single time step. This transition probability is normalized according to

[TABLE]

For the present mass transport model, the transition probability is given by

[TABLE]

where

[TABLE]

with the shorthand notations $S_{i}\equiv m_{i}+m_{i+1}$ and $S_{i}^{\prime}\equiv m_{i}^{\prime}+m_{i+1}^{\prime}$.

3.2 Steady-state distribution

In the following, we show that the distribution

[TABLE]

is a stationary solution of the master equation Eq. (3). In Eq. (8), $M$ is the (constant) total mass, and $Z_{N}(M)$ is a normalization factor. In some cases, it may be convenient to write $P(\mathbf{m})$ in the form $P(\mathbf{m})=\frac{1}{2}[P_{1}(\mathbf{m})+P_{2}(\mathbf{m})]$ with, for $j\in\{1,2\}$

[TABLE]

having defined

[TABLE]

Using Eq. (8), the master equation (3) reads, taking into account the fact that the dynamics conserves the total mass,

[TABLE]

where, to lighten notations, the Dirac delta function accounting for the total mass conservation is understood.

Expanding the r.h.s. of Eq. (11) into four terms, we evaluate these terms separately, obtaining for $j,k\in\{1,2\}$ (see A)

[TABLE]

The sum of the four contributions appearing in the r.h.s. of Eq. (11) is thus equal to $Q_{1}(\mathbf{m}^{\prime})+Q_{2}(\mathbf{m}^{\prime})$, so that Eq. (11) is satisfied. Hence the distribution $P(\mathbf{m})$ given in Eq. (8) is the stationary solution of the model.

3.3 Physical interpretation of the dynamics

Without loss of generality, one can rewrite the functions $v(m)$ and $w(m)$ as

[TABLE]

where we have defined

[TABLE]

The parameter $\beta$, to be thought of as an inverse temperature, is arbitrary here, and has only been introduced to facilitate the comparison with equilibrium. A symmetric redistribution process, obtained for $v(m)=w(m)$, corresponds to $h(m)=0$, and the stationary distribution Eq. (8) boils down to an equilibrium distribution,

[TABLE]

The function $\varepsilon(m)$ thus appears as an effective local energy associated to a local density $m$. The function $h(m)$ describes the asymmetry of the dynamics. In the linear case $h(m)=h_{0}m$, having in mind local detailed balance, the term $2h_{0}(m_{i}-m_{i}^{\prime})$ that enters the ratio $\varphi(m_{i}^{\prime}|S_{i})/\varphi(m_{i}|S_{i}^{\prime})$ (with $S_{i}=S_{i}^{\prime}$) can be interpreted as the work done by a driving force $f=2h_{0}$ associated with a displaced mass $m_{i}-m_{i}^{\prime}$ on a unit distance (one lattice spacing). This case is thus physically meaningful, and we will focus on it when dealing with specific examples (keeping $f$ rather than $h_{0}$ as the driving parameter).

When $h(m)\neq 0$, the non-equilibrium steady-state distribution $P(\mathbf{m})$ given in Eq. (11) can be rewritten as

[TABLE]

It may be convenient to rewrite the distribution in a more compact as

[TABLE]

by introducing the global observables

[TABLE]

The presence of the hyperbolic cosine in Eq. (16) yields long-range correlations, that read to leading order in the driving force

[TABLE]

In more intuitive terms, these correlations are generated by the synchronous dynamics over two different partitions of the lattice. More details on the evaluation of the correlation function and on the expression of the pair and single mass distribution can be found in B.

In the following, the arbitrary inverse temperature scale $\beta$ is set to unity, unless stated otherwise.

4 Non-equilibrium free energy

4.1 Large deviation form of the partition function

It is natural to define from the partition function

[TABLE]

a nonequilibrium (intensive) free energy

[TABLE]

if this limit exists. This means that $Z_{N}(M)$ takes at large $N$ a large deviation form

[TABLE]

To evaluate $I(\rho)$, we follow the saddle-node method presented in [32] (note that the standard Gärtner-Ellis theorem cannot be applied in a straightforward way because $Z_{N}(M)$ is not a probability distribution). Plugging into Eq. (20) the Laplace representation of the delta function,

[TABLE]

with $a$ an arbitrary real number, we end up with

[TABLE]

where

[TABLE]

Note that the real part of $\zeta$ (equal to $a$) is chosen small enough for the integrals to converge. One then has, setting $M=N\rho$ with $\rho$ the average density,

[TABLE]

where we have introduced the function222We take here a determination of the logarithm in the complex plane such that the integration path does not cross the branch cut.

[TABLE]

The function $\lambda(\zeta)$ plays the same role as the scaled cumulant generating function in the Gärtner-Ellis theorem [32]. Assuming that $\lambda(\zeta)-\rho\zeta$ has a saddle-point $\zeta^{\ast}(\rho)$ defined by

[TABLE]

one can choose $a=\zeta^{\ast}$ in the integral on the rhs of Eq. (26), leading through a saddle-point evaluation of the integral to

[TABLE]

Note that the saddle-point $\zeta^{\ast}$ is necessarily real because $Z_{N}$ is real. The large deviation function $I(\rho)$ introduced in Eq. (21) is then given by

[TABLE]

Since the saddle-point $\zeta^{\ast}$ corresponds to a minimum of the function $\rho\zeta-\lambda(\zeta)$ along a line parallel to the imaginary axis, it also corresponds to a maximum of this function along the real axis, so that

[TABLE]

where $\mathcal{D}\subset\mathbb{R}$ is the domain of definition of $\lambda(\zeta)$ over the real axis.

As an example, we evaluate explicitly the large deviation function in the specific case of linear functions $\varepsilon(m)=\varepsilon_{0}m$ and $h(m)=\frac{1}{2}fm$, using the parameterization Eq. (13) of the functions $v(m)$ and $w(m)$. One obtains

[TABLE]

and

[TABLE]

The saddle-point $\zeta^{\ast}(\rho)$, as defined in Eq.(28), is given by

[TABLE]

so that the large deviation function reads, from Eq. (30),

[TABLE]

Note that here and in what follows, we emphasize the $f$-dependence of the large deviation function by denoting it as $I(\rho,f)$ when considering the specific case $h(m)=\frac{1}{2}fm$. At equilibrium, for $f=0$, the large deviation function reduces to the equilibrium free energy $I(\rho,0)=\epsilon_{0}\rho-1-\ln\rho$ (we recall that temperature is set to unity).

4.2 Pressure and chemical potential

We have seen above that the large deviation function $I(\rho)$ plays the role of a nonequilibrium free energy density. The associated extensive free energy simply reads

[TABLE]

From this non-equilibrium free energy, one can define, by analogy with equilibrium, a non-equilibrium thermodynamic pressure $p$ and a non-equilibrium chemical potential $\mu$ [28, 33]

[TABLE]

(note that $N$ plays here the role of the volume). From Eq. (36), one then obtains, using also Eqs. (28) and (30),

[TABLE]

From these definitions, it follows that the non-equilibrium (intensive) free energy $I(\rho)$ can be expressed as in equilibrium (Euler relation)

[TABLE]

In a non-equilibrium context, this relation was also postulated in [28]. In the example $\varepsilon(m)=\varepsilon_{0}m$ and $h(m)=\frac{1}{2}fm$, one finds from Eq. (34) and (35)

[TABLE]

5 Dynamical characterization of the ‘distance’ to equilibrium

We discuss here two different dynamical measures of how far the system is from equilibrium, namely, dynamical quantities that vanish at equilibrium. Note that the quantities we compute are not necessarily positive, but their absolute value might be interpreted as a ‘distance’ to equilibrium333We use here the term ‘distance’ in a loose sense, since the quantities considered do not satisfy mathematical properties (like symmetry under exchange) of a true distance.. The first quantity is the average mass flux $\Phi$, which is directly related to the bias in the redistribution probability. The second one is the entropy production rate $\sigma$, which is rather a measure of the breaking of detailed balance, or in other words, a global measure of probability fluxes in configuration space.

5.1 Stationary mass flux

We now evaluate the stationary mass flux between two sites $i$ and $i+1$ (which, due to mass conservation, is independent of $i$). During a given time step, a mass is transferred between $i$ and $i+1$ only if the link $(i,i+1)$ belongs to the chosen partition ($\mathcal{P}_{1}$ or $\mathcal{P}_{2}$) of the lattice; mass transfer thus occurs with probability $\frac{1}{2}$. The average flux $\Phi$ then reads

[TABLE]

where $m_{i}$ is the mass on site $i$ before a redistribution occurs on the link $(i,i+1)$, while $m_{i}^{\prime}$ is the mass on site $i$ after the redistribution. The masses $m_{i}$ and $m_{i+1}$ before redistribution are assumed to follow the steady-state distribution $P(m_{i},m_{i+1})$ given in B —see Eq. (72); one thus has $\langle m_{i}\rangle=\rho$. Note that the time step has been set to unity.

The average mass $\langle m_{i}^{\prime}\rangle$ after redistribution can be expressed as

[TABLE]

After some algebra, one finds

[TABLE]

The calculation can be carried out explicitly on the example $\varepsilon(m)=\varepsilon_{0}m$ and $h(m)=\frac{1}{2}fm$, yielding

[TABLE]

The average mass flux $\Phi$ then reads, using Eqs. (41) and (43),

[TABLE]

Also, using the explicit expression of $\mu(\rho)$ (using Eq. (41)

[TABLE]

Furthermore, one can notice that the flux, which can be interpreted as a response of the system to the driving force $f$ (when $h(m)=\frac{1}{2}fm$), is directly related to the free energy $I$ (21) as explicitly shown in C :

[TABLE]

5.2 Entropy production rate

An alternative dynamical measure of the degree of irreversibility is given by the entropy production rate. For a discrete time Markov process, the (time-dependent) entropy production rate (i.e., the entropy production per time step) is defined as [34]444Note that notations in [34] do not follow the same convention, as $P(\omega|\omega^{\prime})$ denotes there the probability of a transition from a configuration $\omega$ to a configuration $\omega^{\prime}$, while we use here a (somehow more standard) conditional probability notation where $T(\mathbf{m}|\mathbf{m}^{\prime})$ is the transition probability from $\mathbf{m}^{\prime}$ to $\mathbf{m}$.

[TABLE]

The advantage of this form is that the positivity of $\sigma$ is visible, as it involves products of factors of equal sign. In steady state, the entropy production rate is the opposite of the entropy flow, $\Delta_{\rm int}S=-\Delta_{\rm ext}S$, yielding the simpler expression [34]

[TABLE]

The entropy production rate $\Delta_{\rm int}S$ can be evaluated in the present model, yielding (technical details are reported in D):

[TABLE]

where $P_{j}(\mathbf{m})$ is defined in Eq. (9). One thus recovers, as expected, that $\sigma=0$ at equilibrium, when $P_{1}(\mathbf{m})=P_{2}(\mathbf{m})$. Eq. (52) can be rewritten in terms of the observables $E$ and $H$ defined in Eq. (18), as

[TABLE]

Since the entropy production rate is extensive with system size, it is convenient to define the density $\sigma=\Delta_{\rm int}S/N$ of entropy production rate, in the limit $N\to\infty$. A way to evaluate $\sigma$ in practice is to introduce the generalized partition function $Z_{N}(M,\theta)$, obtained by replacing $h(m)$ by $\theta h(m)$ where $\theta$ is a real parameter, yielding

[TABLE]

Assuming a large deviation form $Z_{N}(N\rho,\theta)\sim e^{-NJ(\rho,\theta)}$, one can then write

[TABLE]

The large deviation function $J(\rho,\theta)$ can be evaluated in the same way as $I(\rho)$, simply replacing $h(m)$ by $\theta h(m)$ in the calculation of $\lambda(\zeta)$ —see Eqs. (27) and (31).

In the specific case $h(m)=\frac{1}{2}fm$, one can also write the entropy production rate in terms of the non-equilibrium free energy $I(\rho,f)$ as

[TABLE]

Given that the flux $\Phi$ is equal to $-\partial I/\partial f$, the entropy production $\sigma$ reads

[TABLE]

using Eq. (49). One then recovers the usual expression of the local entropy production interpreted as the average local work injected in the system (times the inverse temperature that is equal to $1$ here). Note that if the inverse temperature $\beta\neq 1$, one finds $\sigma=\beta f\,\Phi$. This result is consistent with the local detailed balance interpretation of the dynamics briefly discussed in Sect. 3.3.

6 Static characterizations of the ‘distance’ to equilibrium

Having discussed dynamical characterizations of the distance to equilibrium, we turn in this section to static characterizations of this distance, namely, measures of the ’degree of non-equilibrium’ that are based only on the steady-state probability distribution $P(\mathbf{m})$, without any explicit reference to the dynamics.

6.1 Difference of Gibbs free energy functional

One possible such measure is the difference of Gibbs free energy functional between the nonequilibrium and equilibrium distributions, for the same temperature of the thermal bath. Note that for the sake of clarity, we explicitly take into account in this subsection the temperature $T=\beta^{-1}$ (previously set to $T=1$). For an arbitrary probability distribution $P(\mathbf{m})$ over the configuration space of the model, the Gibbs free energy functional $\mathcal{F}[P]$ is defined as

[TABLE]

Given that the equilibrium distribution $P_{\rm eq}(\mathbf{m})$ at temperature $T$ minimizes the functional $\mathcal{F}[P]$, the quantity

[TABLE]

satisfies $\Delta\mathcal{F}\geq 0$ for any distribution $P$ (note that we have introduced the factor $1/N$ to make $\Delta\mathcal{F}$ an intensive quantity). It is thus natural to interpret $\Delta\mathcal{F}$ as a measure of the distance to equilibrium. Note that $\beta N\Delta\mathcal{F}$ identifies with the Kullback-Leibler divergence

[TABLE]

In the present model, a straightforward calculation yields

[TABLE]

The last integral can be evaluated explicitly in the case $h(m)=\frac{1}{2}fm$, where one has

[TABLE]

where $\Psi(M^{\prime}|M)$ is the distribution of the total mass over even sites $M^{\prime}=\sum_{k=1}^{N^{\prime}}m_{2k}$, given the total mass $M$ in the system. By symmetry, the most probable value of $M^{\prime}$ is $M/2$, so that by a saddle-point argument, the last integral in Eq. (62) is equal to zero at order $N$, with only possible subextensive corrections. One thus finds from Eqs. (61) and (21), for $N\to\infty$,

[TABLE]

so that $\Delta\mathcal{F}$ also identifies in this case with the difference of free energy as defined by the large deviation function $I(\rho,f)$ of the partition function $Z_{N}(M)$ —a quantity a priori distinct from the Gibbs free energy functional, as seen from Eq. (61).

6.2 Non-equilibrium order parameter

A non-equilibrium order parameter $\Psi$ has been introduced by Sasa and Tasaki [28] as (the opposite of) the derivative of the non-equilibrium free energy with respect to the driving force. In the present model with $h(m)=\frac{1}{2}fm$, this definition leads to

[TABLE]

Several remarks are in order here. First, this definition is similar to the relation linking, at equilibrium, an order parameter like the magnetization to its conjugate field, hence the name ‘non-equilibrium order parameter’. Second, an alternative definition, involving the derivation with respect to the (mass or particle) flux, has also been proposed in [28]. Third, we use here an intensive order parameter instead of the extensive order parameter originally introduced in [28].

Since the non-equilibrium free-energy $I(\rho,f)$ is, from symmetry arguments, an even function of $f$, $\Psi(\rho,f)$ is an odd function of $f$, and thus vanishes for $f=0$, consistently with the interpretation of $\Psi$ as a non-equilibrium order parameter.

Using Eq. (49), the non-equilibrium order parameter $\Psi$ simply boils down to the mass flux,

[TABLE]

Although $\Psi$ turns out to be equal to $\Phi$, the two quantities differ in essence: $\Psi$ is a static order parameter, while the flux $\Phi$ is a dynamical quantity. Introducing explicitly a time step $\Delta t$ in the model (this time step is been set to $\Delta t=1$ up to now), we would have $\Phi=\Psi/\Delta t$, showing that both quantities have different dimensions.

7 Discussion and conclusion

In this paper, we have introduced a mass transport model with synchronous dynamics for which the steady-state distribution takes a simple non-factorized form, and can be determined explicitly. The knowledge of the steady-state distribution allows for a straightforward evaluation of local distributions of mass, and of a non-equilibrium free energy. The main advantages of this model are on the one hand the simplicity of calculations, and on the other hand the explicit dependence of the steady-state distribution on the driving field —at odds with, for instance, the ZRP and related mass transport models [17].

In addition, we have evaluated several quantities, either static or dynamic, that characterize the ‘degree of non-equilibrium’ of the steady state of the system. These include the mass flux $\Phi$, the entropy production rate per site $\sigma$, the difference $\Delta\mathcal{F}$ of Gibbs free energy functional (per site) between the non-equilibrium and equilibrium states, as well as the non-equilibrium order parameter $\Psi$ introduced by Sasa and Tasaki [28] as the derivative of the non-equilibrium free-energy with respect to the driving force. We have found that all these non-equilibrium parameters are closely related one to the other, and that (at least in the case of a density-independent driving force $f$) the non-equilibrium order parameter $\Psi$ may be seen as a key parameter from which the others can be evaluated. In particular, we have found that

[TABLE]

For a non-zero applied force $f$, all these parameters have a non-zero value. The present mass transport model may thus be considered as a genuine non-equilibrium model. This is to be contrasted, for instance, with more standard mass transport models [19, 17] (including the ZRP) which, in spite of the presence of a non-zero particle flux, have vanishing values of $\Psi$ and $\Delta\mathcal{F}$, because their steady-state distribution is independent of the driving.

Future work may consider possible extensions of the model with asynchronous dynamics, where more complicated forms of the steady-state distribution (involving, e.g., matrix-product states) are likely to be needed. Applications of the model to the field of glassy dynamics could also be considered, by including kinetic constraints in the spirit of the model introduced in [26].

Appendix A Evaluation of the integral terms in the master equation

Calculations of the integrals appearing in the steady-state master equation, as formulated in Eq. (11), are straightforward. We provide here the explicit calculation in the case $j=k=1$ [see Eq. (12)], using again the short notation $S_{i}\equiv m_{i}+m_{i+1}$ and $S_{i}^{\prime}\equiv m_{i}^{\prime}+m_{i+1}^{\prime}$:

[TABLE]

Given that

[TABLE]

one eventually obtains

[TABLE]

Calculations for other values of $j,k$ follow the same lines. For instance, for $k=1$ and $j=2$, $v$ and $w$ are exchanged in the l.h.s. of Eq. (68), but the result is the same since the convolution product is commutative.

Appendix B One- and two-site mass distributions in the thermodynamic limit

We derive in this appendix the one- and two-site mass distributions in the limit of an infinitely large system ($N\to\infty$), keeping the average density $\rho$ fixed.

B.1 Joint mass distribution on a pair of sites

The easiest distribution to compute is the joint distribution of masses $P(m_{i},m_{i+1})$ on neighboring sites. Integrating Eq. (8) over the $N-2$ remaining variables $m_{j}$ ($j\neq i,i+1$), one finds

[TABLE]

Using the large deviation form of $Z_{N}$, one finds

[TABLE]

Hence the distribution $P(m_{i},m_{i+1})$ can be written as

[TABLE]

where $C_{2}(\rho)$ is a normalization constant. It is convenient at this stage to introduce the auxiliary distributions $p_{v}(m)$ and $p_{w}(m)$ defined as

[TABLE]

where $c_{v}$ and $c_{w}$ are normalization constants. In this way, the distribution $P(m_{i},m_{i+1})$ given in Eq. (72) can be reformulated as

[TABLE]

The same calculation holds for the joint distribution $P_{j}(m_{i},m_{i+j})$ of the masses $m_{i}$ and $m_{i+j}$ on distant sites $i$ and $i+j$, as long as $j$ is odd. One thus has

[TABLE]

When $j$ is even, the calculation is slightly more complicated; one has

[TABLE]

with $N^{\prime}=N/2$ and where the quantity $Z_{N_{1},N_{2}}(M)$ is defined as

[TABLE]

However, in the limit $N^{\prime}\to\infty$, the two prefactors $Z_{N^{\prime}-2,N^{\prime}}/Z_{N}$ and $Z_{N^{\prime},N^{\prime}-2}/Z_{N}$ have the same limit, again given by Eq. (71). Hence the distribution reduces to

[TABLE]

Using the more physically meaningful parameterization in terms of the functions $\varepsilon(\rho)$ and $h(\rho)$, the distribution $P_{j}(m_{i},m_{i+j})$ can also be written for all $j>0$ in the form

[TABLE]

As an explicit example, $P_{j}(m_{i},m_{i+j})$ reads in the specific case $\varepsilon(m)=\varepsilon_{0}m$ and $h(m)=h_{0}m$

[TABLE]

where $\mu(\rho)$ is given by Eq. (41).

B.2 Two-point correlation

The two-point correlation function $G_{j}$ between the masses $m_{i}$ and $m_{i+j}$, defined as

[TABLE]

then takes a simple form. From Eqs. (75) and (78), one has for $k>0$

[TABLE]

where $\langle\dots\rangle_{v}$ and $\langle\dots\rangle_{w}$ are averages over the distributions $p_{v}(m)$ and $p_{w}(m)$ respectively. Obviously, $G_{j}$ is $2$-periodic for $j>0$. In the example $\varepsilon(m)=\varepsilon_{0}m$ and $h(m)=h_{0}m$, $G_{j}$ is given by

[TABLE]

In the limit where $f=2h_{0}$ is small, one can expand $G_{j}$ to leading order, yielding

[TABLE]

B.3 Single-site distribution

The single-site distribution $p(m)$ is obtained by integrating the two-site distribution over one of the masses. Using for instance Eq. (74), we get

[TABLE]

or equivalently, in terms of $\varepsilon(m)$ and $h(m)$,

[TABLE]

with $c(\rho)$ a normalization constant.

Appendix C Link between the flux $\phi$ and the nonequilibrium free energy $I$

When one goes from a configuration $\mathbf{m}$ to another one $\mathbf{m}^{\prime}$, the local instantaneous current $\Delta_{i,i+1}(\mathbf{m},\mathbf{m}^{\prime})$ that goes to the right on the link $(i,i+1)$ is

[TABLE]

Summing over all links, the total mass transfered during the transition $\mathbf{m}\to\mathbf{m}^{\prime}$, $\Delta(\mathbf{m},\mathbf{m}^{\prime})$, is given by

[TABLE]

On average,

[TABLE]

Since $\int\rmd\mathbf{m}^{\prime}\,T_{k}(\mathbf{m}^{\prime}|\mathbf{m})=1$ ($k=1,2$), the terms involving $\sum_{i=1}^{N}m_{i}$ cancel out. Using Eq. (10) and (12), one gets

[TABLE]

To go further, one needs to use the physical interpretation of the dynamics, given in Sect. 3.3. Indeed, using Eq. (10) and (13), one can notice that in the (linear) case where $h(m)=\frac{1}{2}fm$,

[TABLE]

Eventually, using Eq. (20), the total averaged mass transfered is equal to

[TABLE]

leading to the final expression of the current $\phi$ (mass transfered per site)

[TABLE]

thus proving the relation in Eq. (49).

Appendix D Evaluation of the entropy production rate

In this appendix, we evaluate the entropy production in the model defined in Sect. 2. From Eq. (5), the transition rate $T(\mathbf{m}^{\prime}|\mathbf{m})$ takes the form

[TABLE]

where $T_{1}(\mathbf{m}^{\prime}|\mathbf{m})$ and $T_{2}(\mathbf{m}^{\prime}|\mathbf{m})$ respectively describe redistributions over the partitions $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ of the lattice. For a given configuration $\mathbf{m}$, we define the sets $\mathcal{D}_{1}(\mathbf{m})$ and $\mathcal{D}_{2}(\mathbf{m})$ as the subsets of configurations $\mathbf{m}^{\prime}$ accessible from $\mathbf{m}$ through redistributions over the partitions $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$. More formally, one has for $j\in\{1,2\}$,

[TABLE]

Using the subsets $\mathcal{D}_{1}(\mathbf{m})$ and $\mathcal{D}_{2}(\mathbf{m})$, one can express the ratio of reciprocal, nonzero transition probabilities, so that the entropy production reads, in steady state,

[TABLE]

The ratios of transition rates can be expressed as

[TABLE]

The restriction of the integration domains to the subsets $\mathcal{D}_{1}(\mathbf{m})$ and $\mathcal{D}_{2}(\mathbf{m})$ in Eq. (97) was needed only to be able to properly define the ratio of reverse transition probabilities. Once Eq. (97) is rewritten in terms of the observables $E(\mathbf{m})$ and $H(\mathbf{m})$, the integration domains no longer need to be restricted to these subsets since the transition probabilities $T_{1}(\mathbf{m}^{\prime}|\mathbf{m})$ and $T_{2}(\mathbf{m}^{\prime}|\mathbf{m})$ appearing in the integrals vanish by definition outside the subsets $\mathcal{D}_{1}(\mathbf{m})$ and $\mathcal{D}_{2}(\mathbf{m})$. Hence one has

[TABLE]

The part of the integral involving $E$ is easily shown to vanish. Using the form $P(\mathbf{m})=\frac{1}{2}[P_{1}(\mathbf{m})+P_{2}(\mathbf{m})]$ of the probability distribution —see Eq. (9)— one has thanks to Eq. (12) that $\int\rmd\mathbf{m}\,T_{k}(\mathbf{m}^{\prime}|\mathbf{m})P_{j}(\mathbf{m})=P_{k}(\mathbf{m})$. The $H$-dependent part in Eq. (100) can then be simplified, after a straightforward calculation, to

[TABLE]

which is precisely Eq. (52).