Derivation of Fourier’s law in stochastic energy exchange systems

# On the derivation of Fourier’s law in stochastic energy exchange systems

Pierre Gaspard†, Thomas Gilbert‡ Center for Nonlinear Phenomena and Complex Systems,
Université Libre de Bruxelles, C. P. 231, Campus Plaine, B-1050 Brussels, Belgium
###### Abstract

We present a detailed derivation of Fourier’s law in a class of stochastic energy exchange systems that naturally characterize two-dimensional mechanical systems of locally confined particles in interaction. The stochastic systems consist of an array of energy variables which can be partially exchanged among nearest neighbours at variable rates. We provide two independent derivations of the thermal conductivity and prove this quantity is identical to the frequency of energy exchanges. The first derivation relies on the diffusion of the Helfand moment, which is determined solely by static averages. The second approach relies on a gradient expansion of the probability measure around a non-equilibrium stationary state. The linear part of the heat current is determined by local thermal equilibrium distributions which solve a Boltzmann-like equation. A numerical scheme is presented with computations of the conductivity along our two methods. The results are in excellent agreement with our theory.

###### pacs:
05.10.Gg, 05.20.Dd, 05.60.-k, 05.70.Ln
: J. Stat. Mech.

## 1 Introduction

A key to a comprehensive derivation of transport properties starting from a microscopic theory is to identify the conditions under which scales in a macroscopic system. This problem is especially challenging for interacting particle systems. To establish the conditions under which scales separate is to provide an understanding as to how a large number of particles whose microscopic motion is described by Newton’s equations organise themselves in irreversible flow patterns at the macroscopic level.

In the context of thermodynamics, macroscopic equations such as Fourier’s heat law derive from conservation laws supplemented by phenomenological ones. The former concern say the conservation of mass or energy, and reflect the existence of the corresponding conservation laws at the level of Newton’s equations. Phenomenological laws on the other hand provide linear relations between the currents associated with the flow of conserved quantities and thermodynamic forces in the form of gradients of the conserved quantities, thereby introducing a set of transport coefficients. Though the values of these coefficients can be precisely measured and are usually tabulated for the sake of their use in the framework of applied thermodynamics, they cannot be otherwise determined without knowledge of the underlying dynamics. A first-principles based computation of the transport coefficients consequently requires a deep knowledge and understanding of the dynamics, as well as their statistical properties.

For this purpose, a common procedure in non-equilibrium statistical physics –see e.g. Uhlenbeck’s discussion of Bogoliubov’s approach to the description of a gas in [3] or van Kampen’s views on the role of stochastic processes in physics [4]– is to apply the two-step programme which consists of (i) identifying an intermediate level of description –a mesoscopic scale– where the Newtonian dynamics can be consistently approximated by a set of stochastic equations, and (ii) subsequently analyzing the statistical properties of this stochastic system so as to compute its transport properties.

The first part of this programme was successfully completed in a recent set of papers [1, 2], in which we introduced a class of Hamiltonian dynamical systems describing the two-dimensional motion of locally confined hard-disc particles undergoing elastic collisions with each other. Whereas the confining mechanism prevents any mass transport in these systems, energy exchanges can take place through binary collisions among particles belonging to neighbouring cells. Under the assumption that binary collisions are rare compared to wall collisions, it was argued that, as a consequence of the rapid decay of statistical correlations of the confining dynamics, the global multi-particle probability distribution of the system typically reaches local equilibrium distributions at the kinetic energy of each individual particle before energy exchanges proceed. This mechanism naturally yields a stochastic description of the time evolution of the probability distribution of the local energies in terms of a master equation to be described below. An important property in that respect is that the accuracy of this stochastic reduction can be controlled to arbitrary precision by simply tuning the system’s parameters.

In [1, 2], arguments were presented supporting the result that the heat conductivity associated with this master equation has a simple analytical form, given by the frequency of energy exchanges. In a sense, this result is similar to the fact that, for instance, uniform random walks describing tracer dynamics have diffusion coefficients given in terms of jump probabilities of the walkers. However, contrary to such systems, the stochastic system at hand lacks a special property which facilitates the derivation of such results, namely the gradient condition. A system obeys the gradient condition if there exists a local function such that the current, whether of mass or energy, can be written as the difference of this function evaluated at separate coordinates [5]. In such cases, the diffusion coefficient is determined through a static average only, and is therefore easy to compute. Given that our stochastic system does not verify the gradient condition, it is in fact remarkable that the heat conductivity should take such a simple form.

In this paper, we complete the second part of the programme described above and provide a systematic derivation of the heat conductivity associated with heat transport in the class of stochastic systems derived in [1, 2]. We justify in particular why, in spite of the fact that our systems do not obey the gradient condition, the heat conductivity can still be computed from the Green-Kubo formula through a static average only. We achieve this by an alternative method which consists in considering the stationary heat flux produced by a temperature gradient across the system. We suppose the system is a two-dimensional slab. Along the first dimension, the system has finite extension, with the corresponding borders in contact with stochastic thermal baths at different temperatures. The system may be taken to periodic along the second dimension. As a result, a temperature gradient develops along the first dimension, which induces a heat flux from the warmer border to the colder one. We then set up a scheme to compute the resulting non-equilibrium stationary state and obtain the linear relation connecting the temperature gradient to the heat flux. This scheme proceeds by consistently solving the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy up to pair distributions through a type of Chapman-Enskog gradient expansion. As the heat current associated with our systems depends on neighbouring pairs only, the knowledge of the pair distribution function to first order in the temperature gradient makes possible the computation of the heat current in the linear regime and, consequently, that of the thermal conductivity, thereby completing the derivation of Fourier’s law.

Our theoretical results are furthermore supported by the results of direct numerical simulations of the stochastic system, with excellent agreement. Simulation methods of kinetic processes governed by a master equation have been extensively developed after Gillespie’s original work [6, 7]. Here we describe a method which accounts for continuous energy exchanges among all the pairs of neighbouring cells.

The paper is organised as follows. The master equation is described in section 2, with relevant definitions. Statistical objects are described in section 3, identifying equilibrium and non-equilibrium stationary states. Local temperatures are defined in section 4 and the energy conservation law derived, thereby identifying heat currents. Section 5 offers a first computation of the heat conductivity through the method of Helfand moments. An alternative computation is presented in section 6, solving the BBGKY hierarchy as sketched above. The numerical scheme for the computation of the heat conductivity is discussed in section 7, with a presentation of the results. Finally, conclusions are drawn in section 8.

## 2 Master equation

Consider a lattice of confining two-dimensional cells, each containing a single hard-disc particle, and such that particles in neighbouring cells may perform elastic collisions with each other, thereby exchanging energy. Such mechanical systems with hard-core confining mechanisms have recently been considered in [1, 2] (see figure 1). In these systems, one distinguishes the local dynamics from the interacting dynamics. On the one hand, the local dynamics are characterized by a wall-collision frequency, , which depends on the geometry of the confining cell as well as on the kinetic energy of the moving particles. On the other hand, the interacting dynamics are characterized by the frequency of binary collisions, i.e. the frequency of collisions between neighbouring particles, .

Under specific conditions, the scale separation, , is achieved, which is to say that individual particles typically perform many collisions with the walls of their confining cells, rattling about their cages at higher frequency than that of binary collisions. In this regime, Liouville’s equation governing the time evolution of phase-space densities reduces to a master equation for the time evolution of local energies. Moreover the validity of this reduction is controlled by the scale separation between the two collision frequencies and , and becomes exact in the limit of vanishing binary collision frequency. On the contrary, in the absence of a clear separation between binary and wall collision frequencies, the dynamics is not reducible to such a stochastic description and cannot be addressed in the context of this paper –see [2] for a discussion of that limit.

Throughout this article, we assume the validity of the scale separation between the wall and binary collision frequencies and focus on the stochastic reduction of the energy exchange dynamics, considering, as our starting point, the mesoscopic level description of the time evolution of probability densities as a stochastic evolution. The time evolution is thus specified by a master equation which accounts for the energy exchanges between neighbouring cells. Note that the dimensionality of the dynamics specifies the maximum dimension of the energy cell array we may consider. In the case of the dynamics depicted in figure 1, the array of energy variables would be two-dimensional. However, without loss of generality, we can simplify the description and consider instead a one-dimensional array according to which every cell has two neighbours instead of four, one on each side111Though it is of course possible to consider a two-dimensional array, as for instance shown in figure 1, we will focus our attention on one-dimensional arrays because, in the presence of a temperature gradient along, say, the horizontal direction, no conduction takes place in the vertical direction. The only relevant energy exchange processes happen along the direction of the temperature gradient.. A similar construction can be carried out in three dimensions, starting with systems of confined hard balls.

We can thus leave aside the underlying dynamics and consider instead a system of cells along a one-dimensional axis with energies and let denote the time-dependent energy distribution associated with this system. The time evolution of this object is specified according to

 to0.0pt$$∂tPN(ϵ1,…,ϵN,t)= (1) N∑a=1∫\/d\/η[W(ϵa+η,ϵa+1−η|ϵa,ϵa+1)PN(…,ϵa+η,ϵa+1−η,…,t) −W(ϵa,ϵa+1|ϵa−η,ϵa+1+η)PN(…,ϵa,ϵa+1,…,t)], where, for the sake of the argument, we assume periodic boundary conditions and identify cells and . The kernel specifies the energy transition rates, whereby a pair of energies exchange an amount of energy. In the case of systems such as shown in figure 1, it is given by  W(ϵa,ϵb|ϵa−η,ϵb+η)=2ρmm2(2π)2|Lρ,ρm(2)|∫\/d\/ϕ\/d\/\biR∫^\bieab⋅\bivab>0\/d\/\biva\/d\/\bivb (2) to0.0pt$$×^\bieab⋅\bivabδ(ϵa−m2v2a)δ(ϵb−m2v2b)δ(η−m2[(^\bieab⋅\biva)2−(^\bieab⋅\bivb)2]),

where denotes the unit vector joining particles and with respective velocities and , is the angle between the direction of this unit vector and a reference axis, denotes the position of the center of mass of the two particles, their masses, their radii, and is the configuration-space volume which they occupy, with a parameter characterizing the geometry of the confining cell. Thus is the amount of energy exchanged by the two particles of respective energies and in the collision process.

The expression (2) of extends beyond the mechanical systems described above. It applies to all mechanical models in which locally confined particles interact with their nearest neighbours through hard-core collisions. The dimensionality of the dynamics, whether two or three, is of course relevant to the specific form of the kernel, in particular as far the computation of the velocity integrals goes. However, whether the underlying dynamics is two- or three-dimensional, it can be shown that the transport properties of the energy exchange process to be established below are obtained in similar ways.

In the above expression, one shows that the spatial integral decouples from the velocity integral . The former quantity represents the volume that the center of mass of particles and can occupy given that the two particles are in contact, performing a collision. Now rescaling the time variable by a factor , which amounts to converting time to the units of the inverse of the square root of an energy, we can rid equation (2) of all geometric factors and write the expression of the kernel in terms of Jacobi elliptic functions of the first kind, denoted . Assuming , the kernel takes the expression

 W(ϵa,ϵb|ϵa−η,ϵb+η)=√2π3×⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩√1ϵaK(ϵb+ηϵa),−ϵb<η<ϵa−ϵb,√1ϵb+ηK(ϵaϵb+η),ϵa−ϵb<η<0,√1ϵbK(ϵa−ηϵb),0<η<ϵa. (3)

For , the kernel is defined in a similar way, using the symmetry of with respect to its arguments, i.e. exchanging and and changing the sign of .

The rate at which two neighbouring cells with energies and exchange energy is obtained by integrating over the range of possible values of the amount of energy exchanged,

 ν(ϵa,ϵb)=∫ϵa−ϵb\/d\/ηW(ϵa,ϵb|ϵa−η,ϵb+η), (4)

which, using equation (3), is easily computed in terms of first and second kind elliptic functions, and , to read

 ν(ϵa,ϵb)=√8ϵbπ3[2E(ϵaϵb)−(1−ϵaϵb)K(ϵaϵb)](ϵb≥ϵa). (5)

This is a symmetric function of , , i.e. if . It is homogeneous in the sense that, given a positive constant , we have .

Likewise, the average heat exchanged between two cells at respective energies and , is

 j(ϵa,ϵb) ≡ ∫ϵa−ϵb\/d\/ηηW(ϵa,ϵb|ϵa−η,ϵb+η), = 23(ϵa−ϵb)ν(ϵa,ϵb), = 2(ϵa−ϵb)3√8ϵbπ3[2E(ϵaϵb)−(1−ϵaϵb)K(ϵaϵb)](ϵb≥ϵa),

which is an antisymmetric function, i.e. defined by when .

## 3 Stationary states

A stationary state of (1) is a time-independent probability distribution such that, for all pairs ,

 ∫\/d\/η[W(ϵa+η,ϵb−η|ϵa,ϵb)P(ss)N(…,ϵa+η,…,ϵb−η,…) −W(ϵa,ϵb|ϵa−η,ϵb+η)P(ss)N(…,ϵa,…,ϵb,…)]=0. (7)

### 3.1 Equilibrium states

Equilibrium measures have densities which depend on the local energies through their sum only, e.g. the micro-canonical distribution and the canonical distribution . The former is associated with an isolated -cells system, the latter to a system in contact with a thermal bath at fixed inverse temperature . Note that we assume the temperature is measured in the units of the energy, or, equivalently, that the Boltzmann constant is unity, .

The equilibrium average of the rate of energy exchanges (4) is the collision frequency, which in the canonical ensemble at temperature , is given by

 νB=1T2∫\/d\/ϵa\/d\/ϵbν(ϵa,ϵb)exp[−(ϵa+ϵb)/T]=√T. (8)

Hence the choice of the time scale.

We further notice that, for finite size systems of cells and total energy , i.e. average energy per cell, the collision frequency is obtained from the micro-canonical average,

 to0.0pt$$νB =1(NT)2∫NT0\/d\/ϵa∫NT−ϵa0\/d\/ϵbν(ϵa,ϵb)(N−1)(N−2)(1−ϵa+ϵbNT)N−3, (9) =√T+O(N−1). For future sake, we also notice the identity  ∫\/d\/x\/d\/y(x−y)2ν(x,y)exp(−x−y)=3. (10) ### 3.2 Non-equilibrium states A non-equilibrium stationary state occurs e.g. when the two ends of a one-dimensional channel are put in contact with thermal baths at inverse temperatures and , with , so that a stationary heat flux will flow from the hot to the cold end. One might hope that the corresponding stationary state would have the local product structure  P(ss)N(ϵ1,…,ϵN)=β1⋯βNexp[−β1ϵ1−⋯−βNϵN], (11) with a temperature profile specified according to Fourier’s law . However such a distribution is not stationary since an energy exchange displaces an amount of energy between two cells at different local temperatures. We will however see in Section 6 that (11) is indeed a good approximation to the actual non-equilibrium stationary state and can be used to compute the linear response. ## 4 Conservation of energy In order to analyze heat transport in our system, we need to consider the time evolution of local temperatures  kBT(a,t)≡⟨ϵa⟩. (12) We have, using equation (1), and assuming we have a one-dimensional lattice of cells so that only terms and contribute,  ∂tkBT(a,t)=∫\/d\/ϵ1…\/d\/ϵNϵa∂tPN(ϵ1,…,ϵN,t), =12∑a′,b′∫\/d\/ϵ1…\/d\/ϵN\/d\/η[ϵaW(ϵa′+η,ϵb′−η|ϵa′,ϵb′)PN(…,ϵ′a+η,…,ϵb′−η,…,t) −ϵaW(ϵa′,ϵb′|ϵa′−η,ϵb′+η)PN(…,ϵa′,…,ϵb′,…,t)], =∫\/d\/ϵ1…\/d\/ϵN\/d\/η[ϵaW(ϵa+η,ϵa+1−η|ϵa,ϵa+1)PN(…,ϵa+η,ϵa+1−η,…,t) −ϵaW(ϵa,ϵa+1|ϵa−η,ϵa+1+η)PN(…,ϵa,ϵa+1,…,t) +ϵaW(ϵa−1+η,ϵa−η|ϵa−1,ϵa)PN(…,ϵa−1+η,ϵa−η,…,t) −ϵaW(ϵa−1,ϵa|ϵa−1−η,ϵa+η)PN(…,ϵa−1,ϵa,…,t)], =∫\/d\/ϵ1…\/d\/ϵN\/d\/η[−ηW(ϵa,ϵa+1|ϵa−η,ϵa+1+η)PN(…,ϵa,ϵa+1,…,t) +ηW(ϵa−1,ϵa|ϵa−1−η,ϵa+η)PN(…,ϵa−1,ϵa,…,t)]. (13) The two terms on the RHS of the last line denote the energy currents flowing between the neighbouring cells,  ∂tkBT(a,t)=Ja−1,a(t)−Ja,a+1(t), (14) where  Ja,b(t)≡∫\/d\/ϵ1…\/d\/ϵNj(ϵa,ϵb)PN(ϵ1,…,ϵN,t) (15) denotes the average energy flux from cell to cell with respect to the distribution . Thus equation (14) is a conservation equation for the energy. We remark here that we have dropped the dimensional factors, setting the cells lengths to unity. Fourier’s law relates the heat current to the local gradient of temperature through a linear law which involves the coefficient of heat conductivity. We now set out to derive Fourier’s law and thereby compute this quantity. ## 5 Helfand moment The heat conductivity which is associated with heat transport can be computed by considering the linear growth in time of the mean squared Helfand moment for that process, which undergoes a deterministic diffusion in phase space [8]. This is a generalisation of Einstein’s formula which relates the diffusion coefficient to the mean squared displacement. Assuming a one-dimensional lattice of cells with unit distance between neighbouring cells, we may write the Helfand moment associated with heat transport as  H(t)=N∑a=1aϵa(t). (16) Changes in this quantity occur whenever an energy exchange takes place between any two neighbouring cells. Thus considering a stochastic realisation of the energy exchange process, we have a sequence of times at which a given pair of cells performs an energy exchange , i.e. and . The corresponding Helfand moment therefore evolves according to  H(τn)=H(τn−1)+η(ϵkn,ϵkn+1), (17) and has overall displacement  ΔH(τn)≡H(τn)−H(τ0)=n∑i=1η(ϵki,ϵki+1). (18) The thermal conductivity can be computed in terms of the equilibrium average of the mean squared displacement of the Helfand moment according to  κ=limN→∞1N(E/N)2limn→∞⟨12τnΔH(τn)2⟩E/N, (19) where denotes a micro-canonical equilibrium average of the cells system at energy , which involves both energy exchanges and relaxation times between then. When is large, boils down to an average with respect to the canonical distribution at temperature . In equation (19), the time corresponding to events is, by the law of large numbers, the typical time it takes for collision events to occur in a system of cells, which is easily expressed in terms of the collision frequency (8),  limn→∞nτn=NνB. (20) Substituting from equation (18) into equation (19), we obtain the thermal conductivity,  κ=limN→∞1NT2limn→∞[n∑i=1⟨12τnη(ϵki,ϵki+1)2⟩T (21) +2n−1∑i=1n∑j=i+1⟨12τnη(ϵki,ϵki+1)η(ϵkj,ϵkj+1)⟩T], which is equivalent to the corresponding result using the Green-Kubo formula. The first of the two terms in the brackets on the RHS of (21) is a sum of identical static averages, which, for large , is approximated by the canonical equilibrium average  ⟨η(ϵa,ϵb)2⟩T=1νBT2∫\/d\/ϵa\/d\/ϵb\/d\/ηη2W(ϵa,ϵb|ϵa−η,ϵb+η)exp[−(ϵa+ϵb)/T], =2T2. (22) Notice here that, in order to compute the second moment of the average energy exchanged between the two cells, we divided by , the rate at which these exchanges occur. The second term on the RHS of equation (21), on the other hand, is a sum of dynamic averages which are generally difficult to compute. We contend that it goes to zero when . This is to say that cross-correlations of energy transfers exist only so long as the system size is finite. The reason for this is that, in large enough systems, new energies keep entering the dynamic averages in equation (21), as if the systems were in contact with stochastic reservoirs, which is enough to destroy these correlations. In the following section, we turn to the evaluation of the non-equilibrium stationary state and will provide more definitive arguments that concur with these heuristic reasoning. In section 7 we discuss the results of numerical computations of the Helfand moment which provide further insight into the decay of these dynamic averages as . Accepting the claim that dynamic averages vanish in equation (21) and thus that only static averages contribute to the heat conductivity, we conclude that the thermal conductivity associated with the process defined by (1) is equal to the collision frequency222The dimensions of these quantities differ by a length squared –the separation between the energy cells– which we have set to unity., which, for a general equilibrium temperature , is  κ=νB=√T. (23) ## 6 Chapman-Enskog gradient expansion Consider the marginals of the -cell distribution function, given by the one- and two-cell distribution functions,  P(1)a(ϵa,t)=∫∏i≠a\/d\/ϵiPN(ϵ1,…,…,ϵN,t),P(2)a,b(ϵa,ϵb,t)=∫∏i≠a,b\/d\/ϵiPN(ϵ1,…,…,ϵN,t). (24) The time evolution of the former can be written in terms of the latter,  ∂tP(1)a(ϵ,t) = ∫∏i≠a\/d\/ϵi∂tPN(ϵ1,…,ϵa−1,ϵ,ϵa+1,…,ϵN,t), (25) = ∫\/d\/ϵ′\/d\/η[W(ϵ+η,ϵ′−η|ϵ,ϵ′)P(2)a,a+1(ϵ+η,ϵ′−η) −W(ϵ,ϵ′|ϵ−η,ϵ′+η)P(2)a,a+1(ϵ,ϵ′) +W(ϵ+η,ϵ′−η|ϵ,ϵ′)P(2)a,a−1(ϵ+η,ϵ′−η) −W(ϵ,ϵ′|ϵ−η,ϵ′+η)P(2)a,a−1(ϵ,ϵ′)]. This equation is the first step of the BBGKY hierarchy. We want to solve this equation in the non-equilibrium stationary state which is generated by putting the system boundaries in contact with thermal reservoirs at two different temperatures, , with respective distributions and , at the corresponding inverse temperatures and . In the presence of such a temperature gradient, we expect that a uniform stationary heat current, as defined by equation (15) and henceforth denoted , will establish itself throughout the system. Provided the local temperature gradient is small, this current would be given, according to Fourier’s law, by the product of the local temperature gradient and heat conductivity :  JH=−κ(Ta+1−Ta)+O[(T+−T−)/N]2. (26) This equality should hold for any , and corresponding to the two thermal reservoirs. Notice that the heat conductivity is itself a function of the temperature and thus varies from site to site. In order to derive Fourier’s law (26), as well as the expression of the heat conductivity , and thus establish equation (23), we compute the heat current using equation (15), which, in terms of the two-cell distribution , becomes  JH≡∫\/d\/ϵ\/d\/ϵ′j(ϵ,ϵ′)P(2)a,a+1(ϵ,ϵ′). (27) As noticed above, on the LHS of this equation is defined to first order in the local temperature gradients. We therefore need to compute to that order as well, which is to say it must be a stationary state of equation (25) up to second order in the local temperature gradient. In order to find this approximate stationary state, we perform a cluster expansion of the two cell distribution function and start by assuming that it factorises into the product of two one-cell distribution functions,  P(2)a,b(ϵ,ϵ′)=P(1)a(ϵ)P(1)b(ϵ′). (28) This solution is as yet incomplete and must be understood as being only the first step in obtaining the consistent solution which is to be written down below. Substituting equation (28) into equation (25), it thus reduces to a Boltzmann-Kac equation, which, owing to the symmetries of the kernel, can be written as  ∂tP(1)a(ϵ,t) =∫\/d\/ϵ′\/d\/ηW(ϵ,ϵ′|ϵ−η,ϵ′+η)[P(1)a(ϵ+η)P(1)a+1(ϵ′−η) −P(1)a(ϵ)P(1)a+1(ϵ′)+P(1)a(ϵ+η)P(1)a−1(ϵ′−η)−P(1)a(ϵ)P(1)a−1(ϵ′)]. A first approximation to the stationary state of this equation is given by the local thermal equilibrium solution,  P(1)a(ϵ)=βae−βaϵ. (30) By definition of the local temperature, we have  ∫∞0\/d\/ϵϵP(1)a(ϵ)=β−1a=Ta. (31) One verifies that the form (30) is an approximate stationary solution of the Boltzmann-Kac equation (6). Indeed, let denote the equilibrium temperature, . We have  P(1)a(ϵ+η)P(1)a+1(ϵ′−η)−P(1)a(ϵ)P(1)a+1(ϵ′)=βaβa+1e−βaϵ−βa+1ϵ′[eη(βa+1−βa)−1], =β2e−β(ϵ+ϵ′)ηδβ+O(δβ2), (32) Likewise  P(1)a(ϵ+η)P(1)a−1(ϵ′−η)−P(1)a(ϵ)P(1)a−1(ϵ′)=−β2e−β(ϵ+ϵ′)ηδβ+O(δβ2). (33) Therefore, to , the sum of equations (32) and (33) vanishes. This is to say the local thermal equilibrium distributions (30) give approximate stationary states of the Boltzmann-Kac equation (6) to that order. Going back to the two-cell distribution, equation (28), we proceed to compute the next order of the BBGKY hierarchy. One might hope that local equilibrium solutions are also a good approximation to the two-cell distribution function. This is however not the case. We can see this by considering the time evolution of the two-cell distribution,  ∂tP(2)a,a+1(ϵ,ϵ′)=∫\/d\/ηW(ϵ,ϵ′|ϵ−η,ϵ′+η)[P(2)a,a+1(ϵ+η,ϵ′−η)−P(2)a,a+1(ϵ,ϵ′)] (34) +∫\/d\/ϵ′′\/d\/ηW(ϵ,ϵ′′|ϵ−η,ϵ′′+η)[P(3)a−1,a,a+1(ϵ′′−η,ϵ+η,ϵ′)−P(3)a−1,a,a+1(ϵ′′,ϵ,ϵ′)] +∫\/d\/ϵ′′\/d\/ηW(ϵ′,ϵ′′|ϵ′−η,ϵ′′+η)[P(3)a,a+1,a+2(ϵ,ϵ′+η,ϵ′′−η)−P(3)a,a+1,a+2(ϵ,ϵ′,ϵ′′)]. We momentarily suppose that the three-cell distributions can be consistently written as a product of one-cell distributions,  P(3)a,b,c(ϵ,ϵ′,ϵ′′)=P(1)a(ϵ)P(1)b(ϵ′)P(1)c(ϵ′′). (35) Plugging this form into equation (34) and assuming a stationary state, we arrive, after expanding all terms to , to the equation  0=β2e−β(ϵ+ϵ′)δβ{j(ϵ,ϵ′)−β∫\/d\/ϵ′′e−βϵ′′j(ϵ,ϵ′′)+β∫\/d\/ϵ′′e−βϵ′′j(ϵ′,ϵ′′)}. (36) This is however a contradiction since the current does not have the gradient form which (36) implies,  j(ϵ,ϵ′)≠β∫\/d\/ϵ′′e−βϵ′′j(ϵ,ϵ′′)−β∫\/d\/ϵ′′e−βϵ′′j(ϵ′,ϵ′′). (37) Indeed, the only possible solution of this equation is a current given in terms of the difference of two local functions, say , which, clearly, equation (2) does not satisfy. To work around this difficulty, we must go back to equation (28) and perform a cluster expansion to include corrections in the form  P(2)a,a±1(ϵ,ϵ′) =P(1)a(ϵ)P(1)a±1(ϵ′)+(βa±1−βa)Q(2)a,a±1(ϵ,ϵ′), (38) ≡β2e−β(ϵ+ϵ′)[1∓δβϵ′±δββq(2)a,a±1(βϵ,βϵ′)], where the first two terms in the second line come from the expansion of the one-cell distributions and is one plus an expression derived from . Notice that, for definiteness, we require that , implying  ∫\/d\/ϵ′Q(2)a,b(ϵ,ϵ′)=0. (39) Substituting the form (38) into the RHS of equation (25), we must have the cancellation of all terms, and thus  to0.0pt$$0=∫\/d\/ϵ′\/d\/ηe−βϵ′W(ϵ,ϵ′|ϵ−η,ϵ′+η){q(2)a,a+1[β(ϵ+η),β(ϵ′−η)] (40) −q(2)a,a+1(βϵ,βϵ′)−q(2)a,a−1[β(ϵ+η),β(ϵ′−η)]+q(2)a,a−1(βϵ,βϵ′)}.

Given that the system is large, we expect to converge to the forms

 q(2)a,a+1(x,y)=q(x,y)+O(δβ),q(2)a,a−1(x,y)=q(y,x)+O(δβ). (41)

If so, let us consider to be the sum of a symmetric function and an antisymmetric function , , with and . It is obvious that the terms between the brackets of equation (40) cancel each other with respect to the symmetric part . As far as the antisymmetric part, on the other hand, we must have

 ∫∞0\/d\/ye−y[∫x−y\/d\/zW(x,y|x−z,y+z)h(x+z,y−z)−ν(x,y)h(x,y)]=0. (42)

We argue that the only possible solution of equation (42) is a function which depends on the sum of its arguments, which is obviously not an antisymmetric function. Therefore is a symmetric function of its arguments,

 q(y,x)=q(x,y). (43)

With the two-cell distribution (38), we write the three-cell distribution, now including the second order terms of the cluster expansion, as

 P(3)a−1,a,a+1(ϵ,ϵ′,ϵ′′)=P(1)a−1(ϵ)P(1)a(ϵ′)P(1)a+1(ϵ′′) +δβ[P(1)a−1(ϵ)Q(2)a,a+1(ϵ′,ϵ′′)+P(1)a+1(ϵ′′)Q(2)a−1,a(ϵ,ϵ′)], =β3e−β(ϵ+ϵ′+ϵ′′){1+δβ(ϵ−ϵ′′)+δββ[q(βϵ,βϵ′)+q(βϵ′,βϵ′′)]}. (44)

Plugging this expression into equation (34), we obtain an equation containing the terms of equation (36), but this time with additional terms involving ’s. Though this equation is complicated and yields a priori no simple explicit solution for , it is consistent with the fact that does not have the gradient property.

Nonetheless, finding the exact expression of is not necessary for our sake as the symmetry of , equation (43), turns out to be enough information on the two-cell distribution function to compute the stationary heat current (27). Indeed, we may write the stationary current as the average

 JH=∫\/d\/x\/d\/ye−(x+y){1−δββ[y−q(x,y)]}j(xβ,yβ). (45)

Notice from equation (2) that . Thus, to , there are only two possible contributions to the heat current that need to be taken into account. Namely,

 ∫\/d\/x\/d\/ye−(x+y)yj(x,y),∫\/d\/x\/d\/ye−(x+y)q(x,y)j(x,y). (46)

However the second of these expressions identically vanishes since is an anti-symmetric function of its arguments and , as we saw in equation (43), is symmetric. The steady state current thus becomes

 JH=−δββ5/2∫\/d\/x\/d\/ye−(x+y)yj(x,y), (47)

which, by symmetry, can be rewritten as

 JH= 12δββ5/2∫\/d\/x\/d\/ye−(x+y)(x−y)j(x,y), (48) =13δββ5/2∫\/d\/x\/d\/ye−(x+y