A Evaluation of the integral terms in the master equation

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

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

## 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 sites, labelled by , with periodic boundary conditions (); is assumed to be even, namely with integer. On each site , one defines a real positive mass . The model is endowed with a synchronous, discrete time dynamics1 The dynamics proceeds, at each time step , by parallel redistributions of mass between neighboring sites and on one of the two partitions and , randomly chosen with equal probability. Once a partition has been selected, all links belonging to the partition are simultaneously updated. To update a link , a new value of the mass on site is randomly drawn from the distribution

 φ(m′i|Si)=v(m′i)w(Si−m′i)v∗w(Si),Si≡mi+mi+1 (1)

where and are arbitrary positive functions, and is the convolution product of and ,

 v∗w(S)=∫S0\rmdmv(m)w(S−m). (2)

From mass conservation, the mass on site is, after redistribution, .

## 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 of all the masses in the system. The probability density evolves according to the discrete time master equation

 P(m′,t)=∫\rmdmT(m′|m)P(m,t) (3)

with , and where is the probability (density) to jump from configuration to configuration in a single time step. This transition probability is normalized according to

 ∫\rmdm′T(m′|m)=1. (4)

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

 T(m′|m)=12T1(m′|m)+12T2(m′|m) (5)

where

 T1(m′|m) = N′∏k=1φ(m′2k|S2k)δ(S′2k−S2k), (6) T2(m′|m) = N′∏k=1φ(m′2k+1|S2k+1)δ(S′2k+1−S2k+1), (7)

with the shorthand notations and .

In the following, we show that the distribution

 P(m)=1ZN(M)(N′∏k=1v(m2k)w(m2k+1)+N′∏k=1w(m2k)v(m2k+1))δ(N∑i=1mi−M) (8)

is a stationary solution of the master equation Eq. (3). In Eq. (8), is the (constant) total mass, and is a normalization factor. In some cases, it may be convenient to write in the form with, for

 Pj(m)=2ZN(M)Qj(m)δ(N∑i=1mi−M), (9)

having defined

 Q1(m)=N′∏k=1v(m2k)w(m2k+1),Q2(m)=N′∏k=1w(m2k)v(m2k+1). (10)

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

 Q1(m′)+Q2(m′) (11) =12∫\rmdm[T1(m′|m)+T2(m′|m)][Q1(m)+Q2(m)]

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 (see A)

 ∫\rmdmTk(m′|m)Qj(m′)=Qk(m). (12)

The sum of the four contributions appearing in the r.h.s. of Eq. (11) is thus equal to , so that Eq. (11) is satisfied. Hence the distribution 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 and as

 v(m)=e−βε(m)−βh(m),w(m)=e−βε(m)+βh(m) (13)

where we have defined

 e−βε(m)=√v(m)w(m),e−βh(m)=√v(m)w(m). (14)

The parameter , 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 , corresponds to , and the stationary distribution Eq. (8) boils down to an equilibrium distribution,

 P(m)=2ZN(M)e−β∑Ni=1ε(mi)δ(N∑i=1mi−M). (15)

The function thus appears as an effective local energy associated to a local density . The function describes the asymmetry of the dynamics. In the linear case , having in mind local detailed balance, the term that enters the ratio (with ) can be interpreted as the work done by a driving force associated with a displaced mass 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 rather than as the driving parameter).

When , the non-equilibrium steady-state distribution given in Eq. (11) can be rewritten as

 P(m)=2ZN(M)e−β∑Ni=1ε(mi)cosh(N∑i=1(−1)iβh(mi))δ(N∑i=1mi−M). (16)

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

 P(m)=2ZN(M)e−βE(m)cosh[βH(m)]δ(N∑i=1mi−M) (17)

by introducing the global observables

 E(m)=N∑i=1ε(mi),H(m)=N∑i=1(−1)ih(mi). (18)

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

 Gj=fρ≪1(2+(−1)j)(ε0−μ(ρ))2ρ2f2+O((ρf)4). (19)

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

 ZN(M)=∫\rmdm[Q1(m)+Q2(m)]δ(N∑i=1mi−M) (20)

a nonequilibrium (intensive) free energy

 I(ρ)=−limN→∞1NlnZ(Nρ), (21)

if this limit exists. This means that takes at large a large deviation form

 ZN(Nρ)∼e−NI(ρ). (22)

To evaluate , 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 is not a probability distribution). Plugging into Eq. (20) the Laplace representation of the delta function,

 δ(s)=12πi∫a+i∞a−i∞\rmdζeζs (23)

with an arbitrary real number, we end up with

 ZN(M)=12πi∫a+i∞a−i∞\rmdζe−ζM[^v(ζ)^w(ζ)]N/2 (24)

where

 ^v(ζ)=∫∞0\rmdmeζmv(m),^w(ζ)=∫∞0\rmdmeζmw(m). (25)

Note that the real part of (equal to ) is chosen small enough for the integrals to converge. One then has, setting with the average density,

 ZN(Nρ)=12πi∫a+i∞a−i∞\rmdζeN(λ(ζ)−ρζ) (26)

where we have introduced the function2

 λ(ζ)=12ln[^v(ζ)^w(ζ)] (27)

The function plays the same role as the scaled cumulant generating function in the Gärtner-Ellis theorem [32]. Assuming that has a saddle-point defined by

 dλdζ(ζ∗)=ρ, (28)

one can choose in the integral on the rhs of Eq. (26), leading through a saddle-point evaluation of the integral to

 ZN(Nρ)∼e−N[ρζ∗(ρ)−λ(ζ∗(ρ))]. (29)

Note that the saddle-point is necessarily real because is real. The large deviation function introduced in Eq. (21) is then given by

 I(ρ)=ρζ∗(ρ)−λ(ζ∗(ρ)). (30)

Since the saddle-point corresponds to a minimum of the function along a line parallel to the imaginary axis, it also corresponds to a maximum of this function along the real axis, so that

 I(ρ)=supζ∈D(ρζ−λ(ζ)). (31)

where is the domain of definition of over the real axis.

As an example, we evaluate explicitly the large deviation function in the specific case of linear functions and , using the parameterization Eq. (13) of the functions and . One obtains

 ^v(ζ)=1ε0−ζ+12f,^w(ζ)=1ε0−ζ−12f (32)

and

 λ(ζ)=−12ln((ε0−ζ)2−f24). (33)

The saddle-point , as defined in Eq.(28), is given by

 ζ∗(ρ)=ε0−1+√1+ρ2f22ρ, (34)

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

 I(ρ,f)=ε0ρ−1+√1+ρ2f22+12ln(1+√1+ρ2f22ρ2). (35)

Note that here and in what follows, we emphasize the -dependence of the large deviation function by denoting it as when considering the specific case . At equilibrium, for , the large deviation function reduces to the equilibrium free energy (we recall that temperature is set to unity).

### 4.2 Pressure and chemical potential

We have seen above that the large deviation function plays the role of a nonequilibrium free energy density. The associated extensive free energy simply reads

 F(M,N)=NI(MN). (36)

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

 p=−∂F∂N,μ=∂F∂M (37)

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

 μ(ρ) = I′(ρ)=ζ∗(ρ), (38) p(ρ) = −I(ρ)+ρI′(ρ). (39)

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

 I(ρ)=−p(ρ)+ρμ(ρ). (40)

In a non-equilibrium context, this relation was also postulated in [28]. In the example and , one finds from Eq. (34) and (35)

 μ(ρ) = ε0−1+√1+ρ2f22ρ (41) p(ρ) = −12ln(1+√1+ρ2f22ρ2). (42)

## 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 equilibrium3. The first quantity is the average mass flux , which is directly related to the bias in the redistribution probability. The second one is the entropy production rate , 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 and (which, due to mass conservation, is independent of ). During a given time step, a mass is transferred between and only if the link belongs to the chosen partition ( or ) of the lattice; mass transfer thus occurs with probability . The average flux then reads

 Φ=12(⟨mi⟩−⟨m′i⟩) (43)

where is the mass on site before a redistribution occurs on the link , while is the mass on site after the redistribution. The masses and before redistribution are assumed to follow the steady-state distribution given in B —see Eq. (72); one thus has . Note that the time step has been set to unity.

The average mass after redistribution can be expressed as

 ⟨m′i⟩=∫∞0\rmdmi∫∞0\rmdmi+1P(mi,mi+1)∫∞0\rmdm′im′iφ(m′i|mi+mi+1). (44)

After some algebra, one finds

 ⟨m′i⟩=2C2(ρ)∫∞0\rmdSe−μS∫S0\rmdm′m′v(m′)w(S−m′). (45)

The calculation can be carried out explicitly on the example and , yielding

 ⟨m′i⟩=1ε0−μ+f. (46)

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

 Φ=12(ρ−⟨m′i⟩)=f4(ε0−μ)2−f2. (47)

Also, using the explicit expression of (using Eq. (41)

 Φ=ρ2f2+2√1+ρ2f2. (48)

Furthermore, one can notice that the flux, which can be interpreted as a response of the system to the driving force (when ), is directly related to the free energy (21) as explicitly shown in C :

 Φ=−∂I(ρ,f)∂f (49)

### 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]4

 Extra open brace or missing close brace (50)

The advantage of this form is that the positivity of 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, , yielding the simpler expression [34]

 ΔintS=∫\rmdm\rmdm′T(m′|m)P(m)lnT(m′|m)T(m|m′). (51)

The entropy production rate can be evaluated in the present model, yielding (technical details are reported in D):

 ΔintS=12∫\rmdm[P1(m)−P2(m)]H(m). (52)

where is defined in Eq. (9). One thus recovers, as expected, that at equilibrium, when . Eq. (52) can be rewritten in terms of the observables and defined in Eq. (18), as

 ΔintS=1ZN(M)∫\rmdmH(m)e−E(m)sinh(H(m))δ(N∑i=1mi−M). (53)

Since the entropy production rate is extensive with system size, it is convenient to define the density of entropy production rate, in the limit . A way to evaluate in practice is to introduce the generalized partition function , obtained by replacing by where is a real parameter, yielding

 ZN(M,θ)=∫\rmdme−E(m)cosh(θH(m))δ(N∑i=1mi−M). (54)

Assuming a large deviation form , one can then write

 σ=−∂J∂θ(ρ,θ=1). (55)

The large deviation function can be evaluated in the same way as , simply replacing by in the calculation of —see Eqs. (27) and (31).

In the specific case , one can also write the entropy production rate in terms of the non-equilibrium free energy as

 σ=−f∂I∂f. (56)

Given that the flux is equal to , the entropy production reads

 σ=fΦ (57)

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 here). Note that if the inverse temperature , one finds . 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 , 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 (previously set to ). For an arbitrary probability distribution over the configuration space of the model, the Gibbs free energy functional is defined as

 Extra open brace or missing close brace (58)

Given that the equilibrium distribution at temperature minimizes the functional , the quantity

 Missing or unrecognized delimiter for \Big (59)

satisfies for any distribution (note that we have introduced the factor to make an intensive quantity). It is thus natural to interpret as a measure of the distance to equilibrium. Note that identifies with the Kullback-Leibler divergence

 D[P||Peq]=∫\rmdmP(m)lnP(m)Peq(m). (60)

In the present model, a straightforward calculation yields

 ΔF=1NlnZeqN(M)−1NlnZN(M)+1N∫\rmdmP(m)lncosh[βH(m)]. (61)

The last integral can be evaluated explicitly in the case , where one has

 ∫\rmdmP(m)lncosh[βH(m)] =∫\rmdM′∫dmP(m)δ(N′∑k=1m2k−M′)lncosh[βf(M−2M′)] =∫\rmdM′Ψ(M′|M)lncosh[βf(M−2M′)] (62)

where is the distribution of the total mass over even sites , given the total mass in the system. By symmetry, the most probable value of is , so that by a saddle-point argument, the last integral in Eq. (62) is equal to zero at order , with only possible subextensive corrections. One thus finds from Eqs. (61) and (21), for ,

 ΔF=I(ρ,f)−I(ρ,0) (63)

so that also identifies in this case with the difference of free energy as defined by the large deviation function of the partition function —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 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 , this definition leads to

 Ψ=−∂I∂f(ρ,f). (64)

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 is, from symmetry arguments, an even function of , is an odd function of , and thus vanishes for , consistently with the interpretation of as a non-equilibrium order parameter.

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

 Ψ(ρ,f)=Φ(ρ,f) (65)

Although turns out to be equal to , the two quantities differ in essence: is a static order parameter, while the flux is a dynamical quantity. Introducing explicitly a time step in the model (this time step is been set to up to now), we would have , 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 , the entropy production rate per site , the difference of Gibbs free energy functional (per site) between the non-equilibrium and equilibrium states, as well as the non-equilibrium order parameter 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 ) the non-equilibrium order parameter may be seen as a key parameter from which the others can be evaluated. In particular, we have found that

 Φ(ρ,f)=Ψ(ρ,f),σ=fΨ(ρ,f),ΔF(ρ,f)=∫f0\rmdf′Ψ(ρ,f′). (66)

For a non-zero applied force , 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 and , 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 [see Eq. (12)], using again the short notation and :

 ∫\rmdmT1(m′|m)Q1(m) (67) =1ZN′∏k=1∫∞0\rmdm2k∫∞0\rmdm2k+1φ(m′2k|S2k)v(m2k)w(m2k+1)δ(S′2k−S2k) =1ZN′∏k=1[v(m′2k)w(m′2k+1)v∗w(S′2k)∫∞0\rmdm2k∫∞0\rmdm2k+1v(m2k)w(m2k+1)δ(S′2k−S2k)].

Given that

 ∫∞0\rmdm2k∫∞0\rmdm2k+1v(m2k)w(m2k+1)δ(S′2k−S2k)=v∗w(S′2k) (68)

one eventually obtains

 ∫\rmdmT1(m′|m)Q1(m)=Q1(m). (69)

Calculations for other values of follow the same lines. For instance, for and , and 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 (), keeping the average density fixed.

### b.1 Joint mass distribution on a pair of sites

The easiest distribution to compute is the joint distribution of masses on neighboring sites. Integrating Eq. (8) over the remaining variables (), one finds

 P(mi,mi+1)=ZN−2(M−mi−mi+1)ZN(M)[v(mi)w(mi+1)+w(mi)v(mi+1)]. (70)

Using the large deviation form of , one finds

 limN→∞ZN−2(M−mi−mi+1)ZN(M)=exp[−2I(ρ)+μ(ρ)(mi+mi+1−2ρ)]. (71)

Hence the distribution can be written as

 P(mi,mi+1)=C2(ρ)eμ(ρ)(mi+mi+1)[v(mi)w(mi+1)+w(mi)v(mi+1)], (72)

where is a normalization constant. It is convenient at this stage to introduce the auxiliary distributions and defined as

 pv(m)=cv(ρ)eμ(ρ)mv(m),pw(m)=cw(ρ)eμ(ρ)mw(m), (73)

where and are normalization constants. In this way, the distribution given in Eq. (72) can be reformulated as

 P(mi,mi+1)=12[pv(mi)pw(mi+1)+pw(mi)pv(mi+1)]. (74)

The same calculation holds for the joint distribution of the masses and on distant sites and , as long as is odd. One thus has

 Pj(mi,mi+j)=12[pv(mi)pw(mi+j)+pw(mi)pv(mi+j)](j=2k−1,k>0). (75)

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

 Pj(mi,mi+j)=ZN′−2,N′(M−mi−mi+j)ZN(M)v(mi)v(mi+j) (76) +ZN′,N′−2(M−mi−mi+j)ZN(M)w(mi)w(mi+j)

with and where the quantity is defined as

 ZN1,N2(M)=∫N1+N2∏i=1\rmdmiN1∏i=1v(mi)N2∏i=N1+1w(mi)δ(N1+N2∑i=1mi−M). (77)

However, in the limit , the two prefactors and have the same limit, again given by Eq. (71). Hence the distribution reduces to

 Pj(mi,mi+j)=12[pv(mi)pv(mi+j)+pw(mi)pw(mi+j)](j=2k,k>0). (78)

Using the more physically meaningful parameterization in terms of the functions and , the distribution can also be written for all in the form

 P(mi,mi+1)=2C2(ρ)e−ε(mi)−ε(mi+j)+μ(ρ)(mi+mi+j)cosh[h(mi)+(−1)jh(mi+j)]. (79)

As an explicit example, reads in the specific case and

 P(mi,mi+j)=[(ε0−μ(ρ))2−h20]2(ε0−μ(ρ))2+(−1)jh20e−(ε0−μ(ρ))(mi+mi+j) (80) ×cosh(h0mi+(−1)jh0mi+j)

where is given by Eq. (41).

### b.2 Two-point correlation

The two-point correlation function between the masses and , defined as

 Gj=⟨mimi+j⟩−ρ2 (81)

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

 G2k−1 = ⟨m⟩v⟨m⟩w−ρ2 (82) G2k = 12(⟨m⟩2v+⟨m⟩2w)−ρ2 (83)

where and are averages over the distributions and respectively. Obviously, is -periodic for . In the example and , is given by

 Missing or unrecognized delimiter for \big (84)

In the limit where is small, one can expand to leading order, yielding

 Gj=fρ≪1(2+(−1)j)(ε0−μ(ρ))2ρ2f2+O((ρf)4)