Heat conduction and the nonequilibrium stationary states of stochastic energy exchange processes

TL;DR

This paper provides a systematic analysis of non-equilibrium stationary states in a solvable heat conduction model, deriving Fourier's law and correlations without dual process analysis, and extends results to generalized gradient models.

Contribution

It offers a new systematic characterization of non-equilibrium states in the KMP model and its generalizations, avoiding dual process analysis.

Findings

Derivation of Fourier's law in the model

Explicit static correlations and covariance matrices

Extension to generalized gradient models

Abstract

I revisit the exactly solvable Kipnis--Marchioro--Presutti model of heat conduction [J. Stat. Phys. 27 65 (1982)] and describe, for one-dimensional systems of arbitrary sizes whose ends are in contact with thermal baths at different temperatures, a systematic characterization of their non-equilibrium stationary states. These arguments avoid resorting to the analysis of a dual process and yield a straightforward derivation of Fourier's law, as well as higher-order static correlations, such as the covariant matrix. The transposition of these results to families of gradient models generalizing the KMP model is established and specific cases are examined.

Figures

Tables

Table 1. Table 1 : Coefficients c r , n subscript 𝑐 𝑟 𝑛 c_{r,n} , equation ( 3.31 ), for degrees r = 2 , … , 5 𝑟 2 … 5 r=2,\dots,5 , obtained by solving equation ( 3.29 ).

$r\backslash n$	$0$	$1$	$2$	$3$	$4$	$5$
$2$	$\frac{19}{234}$	$\frac{5}{234}$	$\frac{19}{234}$
$3$	$\frac{67}{4914}$	$\frac{77}{9126}$	$\frac{-77}{9126}$	$\frac{-67}{4914}$
$4$	$\frac{16231}{1149876}$	$\frac{1394}{287469}$	$\frac{781}{127764}$	$\frac{1394}{287469}$	$\frac{16231}{1149876}$
$5$	$\frac{3824}{729729}$	$\frac{104123}{39670722}$	$\frac{25271}{100039212}$	$\frac{-25271}{100039212}$	$\frac{-104123}{39670722}$	$\frac{-3824}{729729}$

Equations

$$K({\xi }_a,{\xi }_b\to {\xi }_a^{\prime },{\xi }_b^{\prime })=\frac{\nu }{{\xi }_a+{\xi }_b}{\delta }_{\text{textsc}d}({\xi }_a+{\xi }_b-{\xi }_a^{\prime }-{\xi }_b^{\prime }){\Theta }_{\text{textsc}h}({\xi }_a^{\prime }){\Theta }_{\text{textsc}h}({\xi }_b^{\prime }){\delta }_{|a-b|,1},$$

$$\begin{array}{rl}{\xi }_n^{\prime }& =p({\xi }_n+{\xi }_{n+1}),\\ {\xi }_{n+1}^{\prime }& =(1-p)({\xi }_n+{\xi }_{n+1}).\end{array}$$

$$P_{\text{textsc}eq}^{(N)}(\dots ,{\xi }_a,{\xi }_b,\dots )K({\xi }_a,{\xi }_b\to {\xi }_a^{\prime },{\xi }_b^{\prime })=P_{\text{textsc}eq}^{(N)}(\dots ,{\xi }_a^{\prime },{\xi }_b^{\prime },\dots )K({\xi }_a^{\prime },{\xi }_b^{\prime }\to {\xi }_a,{\xi }_b),$$

$$p_{\beta }(\xi )=\beta {\text{rme}}^{-\beta \xi }.$$

$$f({\xi }_1,{\xi }_2)=\int \text{rmd}{\xi }_1^{\prime }\text{rmd}{\xi }_2^{\prime }K({\xi }_1,{\xi }_2\to {\xi }_1^{\prime },{\xi }_2^{\prime })=\nu ,$$

$$\int \text{rmd}{\xi }_1^{\prime }\text{rmd}{\xi }_2^{\prime }({\xi }_1-{\xi }_1^{\prime }{)}^{n}K({\xi }_1,{\xi }_2\to {\xi }_1^{\prime },{\xi }_2^{\prime })=\frac{\nu }{n+1}\frac{{\xi }_1^{n+1}+(-{)}^n{\xi }_2^{n+1}}{{\xi }_1+{\xi }_2}.$$

$$j({\xi }_1,{\xi }_2)=\frac{\nu }{2}({\xi }_1-{\xi }_2),$$

$$\kappa =\frac{\nu }{2}$$

$${\partial }_t{P}_t^{(N)}({\underline{\xi }}_N)=\sum _{a={}^{-N}/{}_2}^{{}^N/{}_2-1}{\hat{\mathcal{L}}}_{a,a+1}{P}_t^{(N)}({\underline{\xi }}_N)+{\hat{\mathcal{L}}}_{{}^{-N}/{}_2}{P}_t^{(N)}({\underline{\xi }}_N)+{\hat{\mathcal{L}}}_{{}^N/{}_2}{P}_t^{(N)}({\underline{\xi }}_N),$$

$${\hat{\mathcal{L}}}_{a,a+1}{P}_t^{(N)}(\{\dots ,{\xi }_a,{\xi }_{a+1},\dots \})=\text{cr}\int \text{rmd}{\xi }_a^{\prime }\text{rmd}{\xi }_{a+1}^{\prime }K({\xi }_a^{\prime },{\xi }_{a+1}^{\prime }\to {\xi }_a,{\xi }_{a+1})P_t^{(N)}(\{\dots ,{\xi }_a^{\prime },{\xi }_{a+1}^{\prime },\dots \})-\nu P_t^{(N)}({\underline{\xi }}_N),$$

$${\hat{\mathcal{L}}}_{{}^{-N}/{}_2}{P}_t^{(N)}(\{{\xi }_{{}^{-N}/{}_2},\dots \})=\text{cr}\int \text{rmd}{\xi }_{-}\text{rmd}{\xi }_{-}^{\prime }\text{rmd}{\xi }_{{}^{-N}/{}_2}^{\prime }K({\xi }_{-}^{\prime },{\xi }_{{}^{-N}/{}_2}^{\prime }\to {\xi }_{-},{\xi }_{{}^{-N}/{}_2})p_{{\beta }_{-}}({\xi }_{-}^{\prime })P_t^{(N)}(\{{\xi }_{{}^{-N}/{}_2}^{\prime },\dots \})-\nu P_t^{(N)}({\underline{\xi }}_N),$$

$${\hat{\mathcal{L}}}_{{}^N/{}_2}{P}_t^{(N)}(\{\dots ,{\xi }_{{}^N/{}_2}\})=\text{cr}\int \text{rmd}{\xi }_{+}\text{rmd}{\xi }_{+}^{\prime }\text{rmd}{\xi }_{{}^N/{}_2}^{\prime }K({\xi }_{{}^N/{}_2}^{\prime },{\xi }_{+}^{\prime }\to {\xi }_{{}^N/{}_2},{\xi }_{+})p_{{\beta }_{+}}({\xi }_{+}^{\prime })P_t^{(N)}(\{\dots ,{\xi }_{{}^N/{}_2}^{\prime }\})-\nu P_t^{(N)}({\underline{\xi }}_N).$$

$$P_{\text{textsc}ne}^{(N)}({\underline{\xi }}_N)=\prod _{a={}^{-N}/{}_2}^{{}^N/{}_2}{p}_{{\beta }_a}({\xi }_a)\sum _{n_{{}^{-N}/{}_2},\dots ,n_{{}^N/{}_2}=0}^{\infty }{\gamma }_{n_{{}^{-N}/{}_2}\dots n_{{}^N/{}_2}}{L}_{n_{{}^{-N}/{}_2}}({\beta }_{{}^{-N}/{}_2}{\xi }_{{}^{-N}/{}_2})\dots L_{n_{{}^N/{}_2}}({\beta }_{{}^N/{}_2}{\xi }_{{}^N/{}_2}).$$

$${\gamma }_{0,\dots ,0}=1.$$

$${\gamma }_{0,\dots ,0,\underset{n}{\underbracket{1}},0,\dots ,0}=0.$$

$$⟨j({\xi }_n,{\xi }_{n+1}){⟩}_{\text{textsc}ne}=\frac{\nu }{2}({\beta }_n^{-1}-{\beta }_{n+1}^{-1}).$$

$$⟨j({\xi }_n,{\xi }_{n+1}){⟩}_{\text{textsc}ne}=-\frac{\nu }{2}\frac{{\beta }_{+}^{-1}-{\beta }_{-}^{-1}}{N+2}.$$

$${\partial }_t⟨{\xi }_n{⟩}_{\text{textsc}ne}=0,$$

$$\int \text{rmd}{\xi }_n\text{rmd}{\xi }_{n\pm 1}{\beta }_n{\beta }_{n\pm 1}\frac{{\xi }_n}{{\xi }_n+{\xi }_{n\pm 1}}{\int }_{-{\xi }_n}^{{\xi }_{n\pm 1}}\text{rmd}\eta {\text{rme}}^{-{\beta }_n({\xi }_n+\eta )-{\beta }_{n\pm 1}({\xi }_{n\pm 1}-\eta )}\text{cr}\times L_p({\beta }_n({\xi }_n+\eta ))L_q({\beta }_{n\pm 1}({\xi }_{n\pm 1}-\eta ))\text{cr}=\frac{1}{2}({\beta }_n^{-1}+{\beta }_{n\pm 1}^{-1}){\delta }_{p,0}{\delta }_{q,0}-\frac{1}{2}{\beta }_n^{-1}{\delta }_{p,1}{\delta }_{q,0}-\frac{1}{2}{\beta }_{n\pm 1}^{-1}{\delta }_{p,0}{\delta }_{q,1};$$

$${\beta }_n^{-1}=\frac{1}{2}({\beta }_{n+1}^{-1}+{\beta }_{n-1}^{-1}).$$

$${\beta }_n^{-1}=\frac{1}{2}({\beta }_{+}^{-1}+{\beta }_{-}^{-1})+\frac{n}{N+2}({\beta }_{+}^{-1}-{\beta }_{-}^{-1}),$$

$${\beta }_m{\beta }_n⟨{\xi }_m{\xi }_n{⟩}_{\text{textsc}ne}-1-{\delta }_{m,n}=\{\begin{array}{ll}{\gamma }_{m:1,n:1},& m<n,\\ 2{\gamma }_{n:2},& m=n,\\ {\gamma }_{n:1,m:1},& m>n,\end{array}$$

$${\gamma }_{\pm {}^N/{}_2\pm 1:2}=0,$$

$$\begin{array}{rl}& {\gamma }_{-{}^N/{}_2-1:1,n:1}=0,\\ & {\gamma }_{n:1,{}^N/{}_2+1:1}=0,\end{array}$$

$${\beta }_{n-1}^{-2}{\gamma }_{n-1:2}+{\beta }_{n+1}^{-2}{\gamma }_{n+1:2}-4{\beta }_n^{-2}{\gamma }_{n:2}+{\beta }_{n-1}^{-1}{\beta }_n^{-1}{\gamma }_{n-1:1,n:1}+{\beta }_{n+1}^{-1}{\beta }_n^{-1}{\gamma }_{n:1,n+1:1}\text{cr}=-2\frac{({\beta }_{+}^{-1}-{\beta }_{-}^{-1}{)}^2}{(N+2{)}^2};$$

$$\frac{1}{2}({\beta }_{n-1}^{-1}{\beta }_{n+1}^{-1}{\gamma }_{n-1:1,n+1:1}+{\beta }_n^{-1}{\beta }_{n+2}^{-1}{\gamma }_{n:1,n+2:1})-\frac{5}{3}{\beta }_n^{-1}{\beta }_{n+1}^{-1}{\gamma }_{n:1,n+1:1}\text{cr}+\frac{1}{3}({\beta }_n^{-2}{\gamma }_{n:2}+{\beta }_{n+1}^{-2}{\gamma }_{n+1:2})=\frac{1}{6}\frac{({\beta }_{+}^{-1}-{\beta }_{-}^{-1}{)}^2}{(N+2{)}^2},$$

$${\beta }_{m-1}^{-1}{\beta }_n^{-1}{\gamma }_{m-1:1,n:1}+{\beta }_m^{-1}{\beta }_{n-1}^{-1}{\gamma }_{m:1,n-1:1}+{\beta }_{m+1}^{-1}{\beta }_n^{-1}{\gamma }_{m+1:1,n:1}\text{cr}+{\beta }_m^{-1}{\beta }_{n+1}^{-1}{\gamma }_{m:1,n+1:1}-4{\beta }_m^{-1}{\beta }_n^{-1}{\gamma }_{m:1,n:1}=0.$$

$${\beta }_m^{-1}{\beta }_n^{-1}{\gamma }_{m:1,n:1}\propto \left(\frac{1}{2}+\frac{m}{N+2}\right)\left(\frac{1}{2}-\frac{n}{N+2}\right),$$

$${\gamma }_{m:1,n:1}={\beta }_m{\beta }_n\frac{({\beta }_{+}^{-1}-{\beta }_{-}^{-1}{)}^2}{N+3}\left(\frac{1}{2}+\frac{m}{N+2}\right)\left(\frac{1}{2}-\frac{n}{N+2}\right),$$

$${\gamma }_{n:2}={\beta }_n^2\frac{({\beta }_{+}^{-1}-{\beta }_{-}^{-1}{)}^2}{N+3}\left[\frac{1}{4}-\frac{n^2}{(N+2{)}^2}+\frac{1}{2(N+2)}\right].$$

$$\frac{1}{4(N+2)}\sum_{n_{1},n_{2}=1}^{2N+3}\big{[}(-1)^{\tfrac{1}{2}(n_{1}+n_{2})}+(-1)^{\tfrac{3}{2}(n_{1}+n_{2})}\big{]}M_{n_{1},n_{2}}\cos\pi(n_{1}x+n_{2}y),$$

Full text

1 Introduction

The Kipnis–Marchioro–Presutti (KMP) model of heat conduction [1] consists of a one dimensional chain of harmonic oscillators which exchange energy among nearest neighbours through stochastic interactions. Considering a finite-size chain in contact with thermal reservoirs at different temperatures, the existence and unicity of the stationary measure was established, proving, in the infinite system-size limit, the convergence of the distribution of energies to a product measure of exponential Gibbs distributions whose temperatures interpolate linearly between the two baths’ temperatures. Moreover, the associated heat flux is proportional to the temperature gradient, with uniform coefficient specifying the heat conductivity.

In many respects, the KMP model is a prototypical example of a system of interacting particles on a lattice which exhibits normal heat conduction and is amenable to an exact solution. In particular, it is one among a few models whose large deviation functional is known, providing a characterisation of macroscopic energy fluctuations about the non-equilibrium stationary state [2].

A potent tool which has proven very useful in establishing the properties of the KMP model is that of duality [3], which allows to reduce the analysis of the model under consideration to that of a lattice gas of particles randomly hopping and mixing among neighbouring sites until they are absorbed at the boundaries; see also reference [2]. The notion that the $r$-point correlation functions in the nonequilibrium steady state can be obtained from absorption probabilities of $r$ dual particles has in the recent years led to a number of other fruitful applications of duality in the context of interacting particle systems [4, 5, 6, 7, 8, 9].

The objective pursued in this paper is to show that the nonequilibrium stationary state of the KMP process can in fact be fully characterised rather easily and without resorting to duality. The technique is based on a straightforward expansion of the stationary state in terms of orthogonal polynomials whose coefficients specify, for a given degree $r$, the $r$-point correlation functions. In so doing, we recover the results already obtained by Bertini et al. [10] for the two-point correlation function. Our technique is, however, more general and yields, for systems of any size, a systematic derivation of the coefficients of arbitrary degree, $r$. They are indeed determined through a closed set of linear equations which involve coefficients of degree $r^{\prime}\leq r$. This is of course consistent with the applicability of duality and is indeed key to its usefulness. Yet the simplicity and straightforwardness of the approach described here seems not to have been duly appreciated.

A remarkable property of the KMP model is that the deterministic part of the current associated with two neighbouring cells with fixed energies, i.e. the first moment of the KMP stochastic kernel associated with its generator, is proportional to the difference of their energies and is therefore similar in form to its thermodynamic counterpart, given by Fourier’s law. This property, which is known as the gradient property [11], is central to the model’s simplicity and the fact that its transport coefficient is given in terms of the current’s average value with respect to local thermal equilibria. Further models of heat transport similar to the KMP process in that they share the gradient property have been considered in recent years; see references [12, 5, 13, 14, 15, 16]. Of particular interest for our sake are so-called Brownian energy processes (BEP) [5], which have been extensively studied in the framework of duality [17]; see also references [18] and [19].

Here, we consider a family of models derived from the instantaneous thermalisation regime of BEP [5, 20], whose local equilibrium distributions are specified in terms of Gamma distributions with arbitrary shape parameters, $\alpha>0$ (the local temperatures specify the scale parameters). Gamma distributions with half integer shape parameters, $\alpha=n/2$, are typically encountered when considering the energy distributions of mechanical systems of $n$ particles [21, Chapter 1]. Our analysis is, however, not constrained by such restrictions; we consider positive real valued shape parameters and, in general, let them take different values in different cells. We show the nonequilibrium stationary states of such systems can in fact be characterised in the same way as that of the KMP process, irrespective of the configuration of shape parameters, wether uniform, periodic or disordered. The consequences of configurational disorder on heat conduction are usually investigated in the framework of harmonic chains [22]. Our models provide a different class and are amenable to analytic results.

The problem of characterising the nonequilibrium stationary states of such systems may however be quite complex since, in general, the temperature profile does not have a simple linear form. On the one hand, periodic shape parameter configurations are quite easy to treat, at least as far as the computation of the heat conductivity goes: it is proportional to the harmonic mean of the shape parameters. Disordered configurations of shape parameters may, on the other hand, lead to pathological cases, in particular, when the shape parameters can have arbitrarily small values. We thus focus more specifically on the study of systems with uniform and alternating shape parameters, which all have linear temperature profiles. An example of such a system was recently considered in reference [15], corresponding to shape parameters alternating between $1$ and $1/2$. There, its dual was shown to be a symmetric simple exclusion process with alternating jump rates [23, 24]; see also references [25, 26] for disordered cases. The linear temperature profile of the energy exchange processes allows us to obtain their two-point correlation functions. In the case of uniform shape parameters, we briefly discuss two limiting cases, corresponding respectively to small and large shape parameter values. In the former case, the limiting regime is such that indivisible energy packets perform random walks in the vicinity of empty (zero energy) cells and merge whenever two energy packets cross each other, which is similar to diffusion-aggregation models [27, 28]. The latter limits to a process where the total energy of the interacting cells is halved evenly among them.

The paper is organised as follows. The general features and properties of the KMP model are recalled in section 2. In section 3, we describe the construction of its stationary state by a polynomial expansion in the energy variables, focusing in sections 3.1 and 3.2 respectively on first and second degree contributions for which general solutions are easily inferred for any system size. In section 3.3, we consider the special case of a single-cell system and compare the results of our approach to the exact solution found in reference [10]. The case of a two-cell system is then considered in section 3.4 which allows to show that explicit solutions can easily be found for arbitrary degrees of the polynomial expansion of the stationary state. The extension of the KMP model to families of such models specified by a configuration of positive real shape parameters is discussed in section 4. In section 4.1 general results are obtained for the first degree terms of the polynomial expansion of the stationary states of such models and different models are discussed. The second degree coefficients are derived in section 4.2 and section 4.3 for uniform and alternating shape parameters respectively. Conclusions are drawn in section 5. The appendices provide a number of technical details pertaining to polynomial expansions of the stationary states.

2 The Kipnis–Marchioro–Presutti model

In the original KMP model [1] the state of a system of $N+1$ cells on a one-dimensional lattice is specified by a collection $\underline{\xi}_{N}\equiv\{\xi_{{}^{-N}\!\!/\!{}_{2}},\dots,\xi_{{}^{N}\!\!/\!{}_{2}}\}$ of positive real variables, $\xi_{i}\in\mathbb{R}_{+}$, interpreted as energies, which are let to interact pairwise through the stochastic kernel

[TABLE]

i.e. such that the combined energy $\xi_{a}+\xi_{b}$ of the interacting nearest neighbouring cells $a$ and $b$ is uniformly redistributed among themselves.

This amounts to picking a pair $\{n,n+1\}$ with uniform rate $\nu$ and drawing a uniformly distributed random number $p\in(0,1)$ such that

[TABLE]

Furthermore, the kernel (2.1) satisfies the detailed balance condition,

[TABLE]

with microcanonical equilibrium distribution specified by the condition $\xi_{{}^{-N}\!\!/\!{}_{2}}+\dots+\xi_{{}^{N}\!\!/\!{}_{2}}=\beta^{-1}(N+1)$, which, as $N\to\infty$, tends to the product of canonical distributions with inverse temperature $\beta$, $P_{\textsc{eq}}^{(N)}(\underline{\xi}_{N})=\prod_{a}p_{\beta}(\xi_{a})$,

[TABLE]

Beyond the uniform rate, identified with the zeroth moment of the kernel (2.1),

[TABLE]

the higher moments are

[TABLE]

In particular, the first moment defines the current,

[TABLE]

which is the gradient of the local energies, a simple instance of the gradient property; see reference [11]. As a consequence, the coefficient of heat conductivity is simply

[TABLE]

In other words, the infinite system-size limit of the heat current is linear in the local temperature gradient and, up to a sign, the coefficient of proportionality is $\kappa$, which is uniform throughout.

This property actually holds also for systems of finite sizes in a nonequilibrium stationary state associated with a temperature gradient. In the following section, we consider systems of arbitrary sizes whose ends are in contact with thermal baths at different temperatures and obtain a systematic characterisation of their stationary state in terms of the energy moments. In particular, we obtain equation (2.8) from the degree-$1$ contributions.

3 Non-equilibrium stationary state

To drive the system away from equilibrium, one lets the two ends of the system be coupled to thermal baths with respective inverse temperatures $\beta_{\pm}$ and distributions $p_{\beta_{\pm}}$, as specified by equation (2.4). Other choices may, however, be more convenient; see equation (3.24) below.

The time evolution of the distribution $P_{t}^{(N)}(\underline{\xi}_{N})$ thus satisfies the following master equation,

[TABLE]

where the local exchange operators are self-adjoint operators, defined by

[TABLE]

and the thermal boundary conditions are specified, on the left-hand boundary, by the operator

[TABLE]

and similarly for right-hand boundary,

[TABLE]

In these expressions the dots in the arguments of the distributions on the right-hand side stand for the components of $\underline{\xi}_{N}$ identical to those appearing on the left-hand side.

The non-equilibrium stationary state $P_{\textsc{ne}}^{(N)}$ of this process can be conveniently expanded in terms of Laguerre polynomials, which form an orthonormal set with respect to the weight function (2.4) and are the natural choice for this model:

[TABLE]

We are thus faced with the problem of determining the parameters $\beta_{n}$ and coefficients $\gamma_{n_{{}^{-N}\!\!/\!{}_{2}}\dots n_{{}^{N}\!\!/\!{}_{2}}}$ in equation (3.5). Given a set of indices $\{n_{{}^{-N}\!\!/\!{}_{2}},\dots,n_{{}^{N}\!\!/\!{}_{2}}\}$, we refer to the sum of the indices $r=\sum_{i}n_{i}$ as the degree of the coefficients. We show below that these coefficients can be determined by sets of equations which are closed, degree by degree in the sense that terms of degree $r$ are determined through a set of linear equations involving coefficients of degree less than or equal to $r$. Coefficients of degree $r$ can therefore be determined exactly, without having to invoke a closure approximation. With hindsight, one realises this property is key to the applicability of duality to the KMP process.

Prior to this, however, we note that, by normalisation of the stationary distribution, we must have:

[TABLE]

Moreover, requiring that the parameters $\beta_{n}$ are the inverse temperatures111Here and in the sequel, we denote the integration with respect to stationary state (3.5) by $\langle\cdot\rangle_{\textsc{ne}}$. , i.e. $\langle\xi_{n}\rangle_{\textsc{ne}}=\beta_{n}^{-1}$, imposes, for all $n$,

[TABLE]

Of course, this does not tell us what are the values of the parameters $\beta_{n}$, which we turn to below.

For ease of notation, we will denote the above elements by $\gamma_{n:1}$, meaning that the $n$th index is $1$ and all the others are zero. This notation extends to other combinations of indices in a self-explanatory way so that, e.g. degree-$2$ coefficients correspond to all combinations of $\gamma_{n:2}$ and $\gamma_{n:1,m:1}$, $n<m$.

3.1 Degree-$1$ contributions: heat

current

It is an immediate consequence of equations (3.6)-(3.7) that the stationary expectation value of the current (2.7) is

[TABLE]

Provided the parameters $\beta_{n}$ correspond to inverse local temperatures (and their pairwise differences to local temperature gradients), the above equation establishes Fourier’s law of heat conduction with uniform heat conductivity (2.8). Indeed, let $J_{\textsc{h}}(N+2)=\sum_{n}\langle j(\xi_{n},\xi_{n+1})\rangle_{\textsc{ne}}$ denote the total of the average currents through the $N+2$ energy pairs (including thermal baths). Equation (3.8) implies $J_{\textsc{h}}(N+2)=-\tfrac{1}{2}\nu(\beta_{+}^{-1}-\beta_{-}^{-1})$. Furthermore, since the current is constant throughout the system, we must therefore have, for every pair of cells,

[TABLE]

It thus remains to show the local temperature gradients are uniform and given by the ratio of the temperature difference of the thermal baths, $\beta_{+}^{-1}-\beta_{-}^{-1}$, by the total number of energy pairs, $N+2$.

To determine the relationship between the parameters $\beta_{n}$ and the baths inverse temperatures $\beta_{\pm}$, we consider the stationarity of the first energy moments,

[TABLE]

which is a short-hand notation for the multiplication of the right-hand side of equation (3.1) by $\xi_{n}$ and integrated over all variables $\xi_{{}^{-N}\!\!/\!{}_{2}},\dots,\xi_{{}^{-N}\!\!/\!{}_{2}}$. Its contributions arise from only two terms in equation (3.1), namely those representing the interaction of cell $n$ with cell $n-1$ ($a=n-1$) and cell $n+1$ ($a=n$). The contributions to each of these two terms come about from terms in equation (3.5) with coefficients $\gamma_{n;p,n+1:q}$ and $\gamma_{n-1:q,n;p}$, i.e. such that every index is zero except for the $n$th and $n\pm 1$th indices. For each such term, we have

[TABLE]

see A for details. Therefore only the terms of degrees [math] and $1$ in the stationary state (3.5) contribute to the first energy moments. Moreover, from equations (3.6) and (3.7), we see that only the first of the three terms on the right-hand side of equation (3.11) is actually relevant. Adding the contributions from the two pairs of cells and subtracting $2\beta_{n}^{-1}$, which arises from the loss term in (3.2), we obtain

[TABLE]

In terms of the baths’ inverse temperatures, this implies

[TABLE]

i.e. the stationary state is characterised by a linear profile of temperatures interpolating between the temperatures of the two baths.

The temperature gradient is therefore uniform across the system, so that equation (3.9) expresses the direct proportionality of the local current and local temperature gradient, with constant $\nu/2$ identified as the thermal conductivity. In particular, for $N\to\infty$, we recover the linear temperature profile obtained in reference [1] by other methods. Correspondingly, we expect the stationary state to converge to a product measure so that all coefficients $\gamma$ in (3.5) must vanish in this limit, at the exception of the zeroth degree coefficient (3.6). In the next section, we turn to the computation of the coefficients of the degree $2$ for finite $N$, recovering the results found in reference [10, Section 2.4].

3.2 Degree-$2$ contributions: covariant

matrix

Thus far, and thanks to the simplicity of equation (3.11), we have not had to cope with the computation of coefficients $\gamma$ involving terms of degree larger or equal to $2$ in the expansion of the stationary state (3.5). The arguments leading to equation (3.11) are in fact more general and extend to arbitrary order; see B. The identification of the degree-$2$ contributions to the stationary state (3.5) occurs from the stationarity of the elements of the covariant matrix,

[TABLE]

with the boundary conditions $\langle\xi_{\pm{}^{N}\!\!/\!{}_{2}\pm 1}\xi_{n}\rangle_{\textsc{ne}}=\beta_{\pm}^{-1}\beta_{n}^{-1}$ and $\langle\xi_{\pm{}^{N}\!\!/\!{}_{2}\pm 1}^{2}\rangle_{\textsc{ne}}=2\beta_{\pm}^{-2}$, i.e.,

[TABLE]

where $-{}^{N}\!\!/\!{}_{2}\leq n\leq{}^{N}\!\!/\!{}_{2}$.

Considering the stationarity of the quadratic energy moments $\langle\xi_{m}\xi_{n}\rangle_{\textsc{ne}}=0\,$, $-\tfrac{N}{2}\leq m\leq n\leq\tfrac{N}{2}$, we obtain results equivalent to those described in reference [10]:

for the diagonal elements,

[TABLE]

see equation (B.11); 2. 2.

the off-diagonal elements, $m=n+1$,

[TABLE]

which is the combination of a contribution of degree $2$ from the pair $(n,n+1)$, as in equation (B.11), and two contributions of degree $1$ from the pairs $(n-1,n)$ and $(n+1,n+2)$, such as in (B.10); 3. 3.

the off-diagonal elements with $|m-n|>1$ are the combinations of contributions degree $1$ (B.10),

[TABLE]

By symmetry of the elements of the covariant matrix, equations (3.16)-(3.18) provide a closed set of $N(N+1)/2$ equations for the coefficients $\gamma_{n:2}$ and $\gamma_{m:1,n:1}$, with $m<n$, which must be solved for the inverse temperatures (3.12) and boundary conditions (3.15).

To infer a solution of this set of equations, notice that equation (3.18) is actually a discrete form of the Poisson equation and admits among its non-trivial solutions bilinear functions of the form

[TABLE]

matching the boundary conditions (3.15b). Remarkably, it is easy to identify a set of solutions of the system (3.16)-(3.18) based on the form (3.19). Matching this form to equations (3.16)-(3.17), it is readily seen that bilinear solutions to equation (3.18),

[TABLE]

also solve equations (3.16)-(3.17), provided the diagonal elements of the covariant matrix are given by

[TABLE]

This solution, however, violates the boundary condition (3.15a). To account for this correction, Fourier modes must in general be added to equations (3.20)-(3.21),

[TABLE]

with amplitudes $M_{n_{1},n_{2}}$ which we will not explicitly compute.

As noted by Bertini et al. [10], the solution (3.21) becomes exact if alternatively one changes the boundary condition (3.15a) to

[TABLE]

which amounts to changing the energy distribution of the baths from the canonical distribution (2.4) to

[TABLE]

Although this is not the canonical distribution expected of a thermal bath, it can be interpreted as the distribution at the interface between the thermal bath and the system. The temperature at the interface is that of the thermal bath but a degree-$2$ correction proportional to the square of the system’s local temperature gradient modifies its energy distribution. Even then, correlations with the energy distribution in the system are absent, as they should.

Pairwise correlations thus take a very simple form in the large system-size limit, decaying with its inverse and falling off according to a quadratic function of the locations of the two cells. These results are similar to those discussed within an analogous framework in reference [29]. They also arise in lattice gases submitted to a density gradient [30].

A numerical implementation of the KMP model subject to boundary conditions (3.24) is straightforward. A comparison between the diagonal elements $\gamma_{n:2}$ of the covariant matrix (3.21) and numerical computations of these quantities is shown in figure 3.1 for systems of different sizes ranging from $N=0$ to $N=63$. The decay of the covariant matrix elements in direct proportion to the size of the system, consistent with their disappearing in the infinite system-size limit, as implied in reference [1], makes them difficult to measure accurately for large system sizes.

As to higher-degree contributions to the non-equilibrium stationary state, the results of B allow their computation from sets of closed linear algebraic equations which can in principle be carried out for any system size. Simple analytic solutions similar to equations (3.20)-(3.21) are, however, more difficult to obtain and it is not immediately clear whether conditions exist on the energy distribution of the baths under which a full-fledge expansion in terms of Fourier modes in the forme of equation (3.22) can be avoided. Their numerical computation is also more difficult than second degree contributions due to their faster decay; see figure 3.2.

3.3 Single-cell system

For $N=0$, i.e. a single cell interacting with two thermostats, it was shown in reference [10] that the stationary state has energy distribution

[TABLE]

where the last line follows by expanding $\rme^{-\beta\,\xi_{0}}$ about $\rme^{-\beta_{0}\,\xi_{0}}$, where $\beta_{0}^{-1}=\tfrac{1}{2}(\beta_{+}^{-1}+\beta_{-}^{-1})$ is the cell’s temperature, whose value is prescribed by the requirement that $P_{\textsc{ne}}^{(0)}$ be normalised. This last form is especially convenient for the sake of evaluating the moments of the distribution.

We may use the results of B to confirm the validity of the solution (3.25). Considering the stationarity of the $n$th energy moment, whose contributions arise from equations (3.3)-(3.4), the relevant terms are obtained from equation (B.8) with $c=n$ and $d=0$. These terms multiply the coefficients $\gamma_{0:p}$ ($q=0$) of the nonequilibrium stationary state (3.5). Their sum for $p=0,\dots,n$, combined with minus twice the $n$th energy moment, must be equal to zero and thus yields the value of $\gamma_{0:n}$ in terms of $\gamma_{0:p}$, $p=0,\dots,n-1$ and parameter $\beta_{0}$,

[TABLE]

Letting $\gamma_{0:0}=1$ and $\gamma_{0:1}=0$, one first computes $\beta_{0}$ from the $n=1$ term, then $\gamma_{0:n}$ for $n\geq 2$. The resulting series (3.25) follows from the identity

[TABLE]

where ${}_{2}{F}_{1}$ denotes the hypergeometric function [31, §15.2].

3.4 Two-cell system

The difficulty of extending the general result (3.25) of Bertini et al. [10] to larger system sizes is already apparent for $N=1$, i.e. two cells, each of which is in contact with a thermostat. Let $r\geq 0$ denote the degree. For each $n=0,\dots,r$ and starting from the bottom, the coefficients of degree $r$ are found by solving the linear system of equations, derived through application of equation (B.8) to operators (3.2)-(3.4) acting on $\xi_{{}^{-1}\!\!/\!{}_{2}}^{n}\,\xi_{{}^{1}\!\!/\!{}_{2}}^{r-n}$,

[TABLE]

Letting $\gamma_{{}^{-1}\!\!/\!{}_{2}:0,\,{}^{1}\!\!/\!{}_{2}:0}=1$ and $\gamma_{{}^{-1}\!\!/\!{}_{2}:1,\,{}^{1}\!\!/\!{}_{2}:0}=\gamma_{{}^{-1}\!\!/\!{}_{2}:0,\,{}^{1}\!\!/\!{}_{2}:1}=0$, one first computes, for $r=1$,

[TABLE]

Considering $r=2,\,3,\dots$, one then finds, for $n=0,\dots,\,r$, the coefficients $\gamma_{{}^{-1}\!\!/\!{}_{2}:n,\,{}^{1}\!\!/\!{}_{2}:r-n}$. Letting

[TABLE]

we obtain the results reported in table 1 and figure 3.3 where values are displayed for $2\leq r\leq 20$, which involves $119$ different coefficients (up to a sign for $r$ odd).

4 Parameter-dependent KMP models

Generalisations of the KMP kernel (2.1) have been obtained as instantaneous thermalisation limits of so-called Brownian energy processes [5]. Here we consider the kernels

[TABLE]

where $\underline{\alpha}=\{\alpha_{-{}^{N}\!\!/\!{}_{2}},\dots,\alpha_{{}^{N}\!\!/\!{}_{2}}\}$ is a sequence of positive parameters, $\alpha_{i}>0$, whose values may depend on the cell index. The KMP kernel (2.1) is a particular case, retrieved for the parameter values $\alpha_{i}=1$, for all $i$. More generally, when $\alpha_{i}$ has half-integer value, it is interpreted has half the number of degrees of freedom involved in the local Brownian energy process. Here we need not make such restrictions and consider $\alpha_{i}\in\mathbb{R}_{+}$.

The kernels (4.1) share important properties of the KMP kernel. They satisfy the detailed balance condition (2.3) with equilibrium distribution now specified by Dirichlet distributions which, for $N$ large, tend to the product of gamma distributions with shape parameter $\alpha_{i}$ and scale parameter given by the temperature $\beta^{-1}$,

[TABLE]

The canonical distribution (2.4) is recovered when $\alpha=1$. Furthermore, although the kernels exhibit explicit dependence on the outgoing energies $\xi^{\prime}_{a}$ and $\xi^{\prime}_{b}$, their zeroth moment is actually uniform,

[TABLE]

and consequently identical to equation (2.5). Irrespective of the actual configuration $\underline{\alpha}$ of shape parameters, all energy pairs thus exchange energy at the same rate $\nu$.

Moreover, the current, given by the first moment of the kernels (4.1),

[TABLE]

exhibits a dependence on the parameters whenever $\alpha_{1}\neq\alpha_{2}$. Thus it might appear that, unless the factor $\alpha_{1}\,\alpha_{2}/(\alpha_{1}+\alpha_{2})$ is independent of the pair of cells through which the measurement is taking place, the current (4.4) is not expressible as the difference of a local function evaluated at the two cells and does not have the gradient property [11]. Nevertheless, more general cases can be considered and their conductivity characterised, as we shall see shortly.

By analogy with equation (3.5), the non-equilibrium stationary state $P_{\textsc{ne}}^{(N)}$ of this process can be expanded about the local equilibria (4.2) in terms of polynomials,

[TABLE]

where the polynomials $J^{\alpha}_{n}$,

[TABLE]

define a complete set derived from the generalised Laguerre polynomials of parameter $\alpha-1$, whose elements are orthonormal with respect to the weight function $p_{\alpha,1}(x)$, as specified by equation (4.2).

4.1 Heat current and temperature profile

Taking the average of the current (4.4) with respect to the stationary state (4.5), we have

[TABLE]

Considering the total current, $J_{\textsc{h}}(N+2)=\sum_{n}\langle j(\xi_{n},\xi_{n+1})\rangle_{\textsc{ne}}$, since we must have $\langle j(\xi_{n},\xi_{n+1})\rangle_{\textsc{ne}}=J_{\textsc{h}}(N+2)/(N+2)$, we find

[TABLE]

Letting $N\to\infty$, provided the overall temperature gradient of the nonequilibrium stationary state is linear, which may depend on the configuration of shape parameters $\underline{\alpha}$, the heat conductivity is given by $\nu/2$ times the harmonic mean of the shape parameters,

[TABLE]

This is remindful of the problem of diffusion in one-dimensional disordered lattices [32].

To determine the local temperatures and verify the linearity of the overall temperature profile, we note that the results obtained for the characterisation of the non-equilibrium stationary state of the KMP model transpose to the kernels (4.1); see C. In particular, considering the stationarity of the first energy moments, we obtain the following identity governing the temperature profile,

[TABLE]

There are two simple cases such that this equation reduces to (3.12):

if all $\alpha_{n}$ are identical, or 2. 2.

if they alternate between two different values, say $\alpha_{0}$ and $\alpha_{1}$, depending on the parity of $n$.

In both cases, the temperature profile is strictly linear and given by (3.13). The current (4.4) has the gradient property and the heat conductivity of the resulting process (4.9) can also be inferred in terms of the static correlations only,

[TABLE]

Perhaps interestingly, this value is equal to the exchange frequency (4.3) when $\alpha_{1}\alpha_{2}=\alpha_{1}+\alpha_{2}$, which has positive solutions for $\alpha_{1},\,\alpha_{2}>1$. Such solutions include, in particular, $\alpha_{1}=\alpha_{2}=2$ and $\alpha_{1}=\tfrac{3}{2}$, $\alpha_{2}=3$, for which this identity may therefore be given a mechanical interpretation.

By extension, a linear temperature profile is also observed

if the sequence of shape parameters is periodic, i.e. such that $\alpha_{n}=\alpha_{n\,\mathrm{mod}\,p}$ where $p\geq 3$ is the period.

However, in such cases, and assuming $N+2$ is a multiple of $p$, periodic excursions of lengths $p$ are observed about the linear temperature profile (3.13),

[TABLE]

Letting $\gamma^{(+)}_{n}=\alpha_{n+1}(\alpha_{n}+\alpha_{n-1})$ and $\gamma^{(-)}_{n}=\alpha_{n-1}(\alpha_{n}+\alpha_{n+1})$, with $n=0,\dots,p-1$, the coefficients $\theta_{0}$, …, $\theta_{p-1}$ are determined through the set of $p$ equations,

[TABLE]

Other cases of interest are those of a random sequences of shape parameters, for which the temperature profiles have the form (4.12)-(4.13), but without the periodicity. Among such models, different classes can be distinguished, in particular:

if the shape parameters are drawn randomly from a finite set, 2. 5.

if the shape parameters are drawn randomly from the (countable) set of positive half integers, or 3. 6.

if the shape parameters are drawn randomly over the real positive numbers.

The last category (6) is a priori problematic, in particular when the values of the shape parameters can be arbitrarily small. Loosely speaking, in such cases, the temperature profile typically exhibits a step-like structure with points of discontinuity at sites where the shape parameters are small. Such a situation is excluded in cases (4)-(5). In case (4), fluctuations about the linear temperature profile are expected to be finite. Case (5) is interesting since, in the framework of Brownian energy processes, it amounts to randomly selecting the number of degrees of freedom at every site (with respect to some probability distribution on the set). Generally speaking, one expects equation (4.9) to hold provided a linear temperature profile is recovered upon local averaging at some intermediate scale, that is, provided $\theta_{n}$ vanishes when locally averaged on this intermediate scale.

We will not dwell further on such considerations which touch upon the broader problem of random walks in random environments [33, 34] and bears similarities with the problem of conduction in random one-dimensional chains [35], as well as, in higher dimensions, the percolation threshold [36]; it deserves a separate study. Rather, we focus below on models (1) and (2) and obtain the second-degree contributions to their non-equilibrium stationary states, from which the contributions to the covariant matrix are deduced.

4.2 Uniform $\alpha$

Letting $\alpha_{i}\equiv\alpha$, we have the temperature profile (3.13) and stationary current,

[TABLE]

with heat conductivity

[TABLE]

Considering the second degree contributions to the nonequilibrium stationary state, one finds that equation (3.18) remains unchanged. To find the second degree contributions, we observe that equations (3.16)-(3.17) transpose to:

[TABLE]

and

[TABLE]

Equations (3.20)-(3.21) thus become to

[TABLE]

and

[TABLE]

As with $\alpha=1$, the boundary condition (3.15a) has to be modified for equation (4.19) to become an exact solution. We must therefore impose

[TABLE]

which amounts to changing the energy distribution of the baths to

[TABLE]

The limit $\alpha\to 0$ is of particular interest. The integrated kernel (4.1) corresponding to this regime,

[TABLE]

is the cumulative distribution function of a regularised Beta distribution with vanishing shape parameter. Its value thus tends to $1/2$ for $-\xi_{b}<\eta<\xi_{a}$. The limiting energy exchange process thus induces a complete transfer of energy to either of the two interacting cells, $\eta=-\xi_{b}$ or $\eta=\xi_{a}$, with probability rate $\nu/2$, viz.

[TABLE]

In other words $p$ in equation (2.2) takes its values in the discrete set $\{0,1\}$, each with probability $1/2$. In an empty neighbourhood, energy packets thus perform a random walk, coalescing whenever two energy packets overlap.

The fact that the heat conductivity (4.15) goes to zero in this limit is a feature of the anomalously slow kinetics of heat transfer which results from the coupling of diffusion and aggregation [37]. As a simple illustration, let $\alpha\equiv\alpha_{N}$,

[TABLE]

with $\alpha_{0}>0$, in which case the stationary current (4.14) scales with $(N+2)^{-2}$.

With this choice of parametrisation, Equations (4.18)-(4.19) yield respectively

[TABLE]

and

[TABLE]

The off-diagonal elements thus behave similarly to equation (4.18), but the diagonal elements decay with the square root of the system size, slower than equation (4.19). At the same time, the distribution (4.21) converges to the usual gamma distribution faster than the coefficient (4.20) when $\alpha$ is fixed.

The opposite limit, $\alpha\to\infty$, is such that the sum of the energies of the interacting cells is exactly halved among them,

[TABLE]

The sources of randomness are thus restricted to the interaction times. As of energy transfers from the thermal boundaries, notice that, for large shape parameters, the energy distributions of the baths are sharply peaked about $\alpha\,\beta_{\pm}^{-1}$. Energies at the boundaries are therefore fixed and, in the bulk, fluctuations of the energies rescaled by the shape parameter, $\alpha^{-1}\xi_{n}$, about the local temperature $\beta_{n}^{-1}$ are inversely proportional to the system size.

Correlations in this regime can be studied, for instance, by assuming

[TABLE]

for which the current (4.14) is independent of $N$ and, similarly to conduction in harmonic chains [38], the heat conductivity (4.15) is infinite.

Whereas the off-diagonal elements of the second-degree contributions to the stationary state (4.18) become independent of $\alpha_{0}$, the diagonal elements (4.19) still decay with $(N+2)^{-1}$ but pick up a uniform contribution proportional to $\alpha_{0}$,

[TABLE]

The second-degree contribution to the energy distribution of the baths (4.20) must therefore scale with $N^{-1}$ rather than $N^{-2}$ for fixed $\alpha$.

4.3 Alternating $\alpha_{0}$ and $\alpha_{1}$

For the sake of the argument, let us think of $N$ as even. Let $i\in\{-\tfrac{N}{2},\dots,\tfrac{N}{2}\}$ and assume $\alpha_{i}=\alpha_{0}$ if $i$ is even and $\alpha_{i}=\alpha_{1}$ if $i$ is odd and denote by $\sigma(i)\in\{0,\,1\}$ the parity of $i$222For $N$ odd, the parity of $i$ should be interpreted as that of $i+\tfrac{1}{2}$. ,

[TABLE]

Equations (3.16)-(3.18) take on the expressions, for $m=n$,

[TABLE]

for $m=n+1$,

[TABLE]

and, for $|m-n|>1$,

[TABLE]

The solutions to equation (4.33) which match the boundary conditions (3.15b) take forms similar to that of equation (3.19), but with an extra factor $\sqrt{\alpha_{0}/\alpha_{1}}$ if both $m$ and $n$ are even, or $\sqrt{\alpha_{1}/\alpha_{0}}$ if both $m$ and $n$ are odd. By matching them to equations (4.31) and (4.32), one obtains the expressions:

[TABLE]

and

[TABLE]

These solutions are exact provided the thermal baths have energy distributions

[TABLE]

5 Concluding remarks

The Kipnis–Marchioro–Presutti model of heat conduction [1] belongs to a larger class of stochastic energy exchange Markov jump processes derived from the instantaneous thermalisation limit of Brownian energy processes [5]. On the one hand, the defining common feature of these models is that every pair of neighbouring cells exchanges energy among them at uniform rate. Their distinctive feature, on the other hand, is that the detailed balance condition is obeyed with respect to different canonical equilibrium energy distributions, which are identified by the sequence of shape parameters associated with every model.

The non-equilibrium stationary states resulting from the application of a temperature gradient at the system’s boundaries are amenable to analytic treatment. Indeed, the simple structure of the stochastic kernel allows for the determination of the stationary states in terms of the products of local canonical equilibrium distributions identified by their shape parameters and a multinomial expansion specified by the orthogonal polynomials associated with them. The collection of all coefficients of degree $r$ in this expansion is obtained by solving sets of linear equations derived by invoking the stationarity of $r$-point correlation functions.

The temperature profile and covariant matrix elements are thus obtained by considering the first and second degree terms in this expansion. Among the models we considered, those specified by sequences of alternating shape parameters yield a linear temperature profile under thermal boundary conditions at different temperatures. Explicit expressions of their covariant matrix elements were obtained upon the condition that the energy distributions associated with the thermal baths have second degree contributions tailored so as to eliminate Fourier components of non-zero wavelengths. This provides a generalisation of the results obtained by Bertini et al. in the context of the KMP model [10].

Furthermore, the strict linearity of the temperature profile, which is linked to the gradient property, is lost for more general sequences of shape parameters. Provided linearity is recovered at some intermediate scale, however, one recovers a simple expression of the heat conductivity in terms of the harmonic mean of the sequence of shape parameters, similar to the diffusion coefficient of disordered lattices [35].

Appendix A Computation of the degree $1$ contributions

To derive equation (3.11), replace the argument of the Laguerre polynomials by a differential operator:

[TABLE]

and likewise for $L_{q}(\beta_{n\pm 1}(\xi_{n\pm 1}-\eta))$. Now compute:

[TABLE]

Equation (3.11) thus becomes

[TABLE]

Evaluating the derivatives and using the identities

[TABLE]

we obtain the announced result.

Appendix B Computation of higher degree contributions

To compute contributions of degree higher than one, we must consider the transposition of equation (A.2) to factors of energy variables raised to arbitrary integer powers,

[TABLE]

To determine the action of the product of Laguerre operators on this expression, let $x=a\beta_{n}$ and $y=b\beta_{n\pm 1}$ and evaluate the results of derivatives at $x=\beta_{n}$ and $y=\beta_{n\pm 1}$. It is sufficient to consider the $x$-part only, for which we have

[TABLE]

where the binomial factor is $(1-s)\dots(p-s)/p!$ and thus vanishes for $p\geq s$.

Finally, the action of the product of the two Laguerre operators on equation (B.1) becomes

[TABLE]

which vanishes identically when $p+q>c+d$ and justifies our assertion that the coefficients in the expansion of the stationary state (3.5) are determined through sets of equations which are closed degreewise.

We proceed to evaluate this expression for the first few degrees.

B.1 Degree-[math] term

Letting $c=d=0$ in equation (B.8), we get back the conservation of probability,

[TABLE]

B.2 Degree-$1$ terms

For $c=1$ and $d=0$ (or equivalently $c=0$ and $d=1$), we retrieve the two contributions to the right-hand side of equation (3.11), which, up to factor $1/2$, are

[TABLE]

B.3 Degree-$2$ terms

For $c=2$ and $d=0$ (or $c=0$ and $d=2$), we have, apart from a factor $2/3$, the contributions

[TABLE]

whose sum yields the right-hand side of equation (3.16).

The contributions corresponding to $c=1$ and $d=1$, given by (B.11) multiplied by $1/3$, bring about one of three contributions to equation (3.17), namely the $(n,n+1)$ interaction, the two others, $(n-1,n)$ and $(n+1,n+2)$, involving corrections of degree $1$ only, given by equation (B.10).

B.4 Degree-$3$ terms

For $c=3$ and $d=0$ (or $c=0$ and $d=3$), we have, up to a factor $3/2$ the contributions

[TABLE]

The same contributions hold for $c=2$ and $d=1$ or $c=1$ and $d=2$, up to a factor one half. Altogether these contributions provide the means to extend our analysis of the stationary state (3.5) to degree $3$.

Appendix C Extension to parameter-dependent models

For the model (4.1), one finds that the right-hand side of equation (B.1) transposes to

[TABLE]

Acting on this expression with the generalised Laguerre operators $J^{\alpha_{0}}_{p}(-\partial_{a})$ and $J^{\alpha_{0}}_{q}(-\partial_{b})$ expanded in closed form, one obtains

[TABLE]

For $i$ and $j$ integers, observe that $\Gamma(\alpha+i+j)/\Gamma(\alpha+j)=(\alpha+j)_{i}$, the Pochhammer symbol, is a polynomial in $\alpha+j$ of degree $i$. It follows from an algebraic identity [39] that the corresponding binomial series in (C.2) vanish whenever the indices $k<c+d-i$ or $l<i$. As a consequence, the right-hand side of (C.2) is zero when $p+q>c+d$. Moreover, if $p\leq i$, the identities (A.4) generalise to

[TABLE]

Thus, let

[TABLE]

with the coefficients $C_{n}(\alpha,i)=0$ if $n>i$. After substituting these expressions in (C.2), we have

[TABLE]

Below, we proceed to evaluate this expression for terms of degrees $1$ and $2$.

C.1 Degree-$1$ terms

For $c=1$ and $d=0$, up to factor $\nu\alpha_{n}/(\alpha_{n}+\alpha_{n\pm 1})$, we have the three contributions

[TABLE]

The addition of the $p=0$ and $q=0$ contributions for the pairs $\{n-1,n\}$ and $\{n,n+1\}$ minus twice $\alpha_{n}\beta_{n}^{-1}$, which arises from the lost term in (3.1), yield equation (4.10).

C.2 Degree-$2$ terms

Considering $c=2$ and $d=0$, up to a common factor

[TABLE]

we have the six contributions

[TABLE]

The same expression (C.8) is obtained for $c=0$ and $d=2$, except for the factor (C.7), which is replaced by

[TABLE]

as well as for $c=1$ and $d=1$, with the factor (C.7) replaced by

[TABLE]

Combinations involving these terms lead to the expressions (4.16)-(4.17) and (4.31)-(4.32).

The author wishes to acknowledge the hospitality of the Erwin Schrödinger Institute, Vienna, on the occasion of the conference Hyperbolic Dynamics and Statistical Physics held in May 2016, where part of this work was presented. He receives financial support from the (Belgian) FRS-FNRS.