Heat conduction and Fourier's law in a class of many particle dispersing billiards

We consider the motion of many confined billiard balls in interaction and discuss their transport and chaotic properties. In spite of the absence of mass transport, due to confinement, energy transport can take place through binary collisions between neighbouring particles. We explore the conditions under which relaxation to local equilibrium occurs on time scales much shorter than that of binary collisions, which characterize the transport of energy, and subsequent relaxation to local thermal equilibrium. Starting from the pseudo-Liouville equation for the time evolution of phase-space distributions, we derive a master equation which governs the energy exchange between the system constituents. We thus obtain analytical results relating the transport coefficient of thermal conductivity to the frequency of collision events and compute these quantities. We also provide estimates of the Lyapunov exponents and Kolmogorov-Sinai entropy under the assumption of scale separation. The validity of our results is confirmed by extensive numerical studies.


Introduction
Understanding the dynamical origin of the mechanisms which underly the phenomenology of heat conduction has remained one of the major open problems of statistical mechanics ever since Fourier's seminal work [1]. Fourier himself actually warned his reader that the effects of heat conduction "make up a special range of phenomena which cannot be explained by the principles of motion and equilibrium," thus seemingly rejecting the possibility of a fundamental level of description. Of course, with the subsequent developments of the molecular kinetic theory of heat, starting with the works of the founding-fathers of statistical mechanics, Boltzmann, Gibbs and Maxwell, and later with the definite triumph of atomism, thanks to Perrin's 1908 measurement of the Avogadro number in relation to Einstein's work on Brownian motion, Fourier's earlier perception was soon discarded as it became clear that heat transport was indeed the effect of mechanical causes. Yet, after well over a century of hard labor, the community is still actively hunting for a first principles based derivation of Fourier's law, some authors going so far as to promise "a bottle of very good wine to anyone who provides [a satisfactory answer to this challenge]" [2]. Thus the challenge is, starting from the Hamiltonian time evolution of a system of interacting particles which models a fluid or a crystal, to derive the conditions under which the heat flux and temperature gradients are linearly related by the coefficient of heat conductivity. This is the embodiment of Fourier's law.
As outlined in [2], it is necessary, in order to achieve this, that the model statistical properties be fully determined in terms of the local temperature, a notion which involves that of local thermalization. To visualize this, we imagine that the system is divided up into a large number of small volume elements, each large enough to contain a number of particles whose statistics is accurately described by the equilibrium statistics at the local temperature associated to the volume element under consideration. To establish this property, one must ensure that time scales separate, which implies that the volume elements settle to local thermal equilibrium on time scales much shorter than the ones which characterize the transport of heat at the macroscopic scale. Local thermal equilibrium therefore relies on very strong ergodic properties of the model.
The natural framework to apply this programme is that of chaotic billiards. Indeed non-interacting particle billiards, where individual tracers move and make specular collisions among a periodic array of fixed convex scatterers (the periodic Lorentz gases) are the only known Hamiltonian systems for which mass transport [3] and shear and bulk viscosities [4] have been established in a rigorous way. Here, rather than local thermal equilibrium, it is a relaxation to local equilibrium, occurring on the constant energy sheet of the individual tracer particles, which allows to model the mass transport by a random walk and yields an analytic estimate of the diffusion coefficient [5]. Arguably Lorentz gases, whether periodic or disordered, in or out of equilibrium, have played, over more than a century, a privileged and most important role in the development of transport and kinetic theories [6]. However, in the absence of interaction among the tracer particles, there is no mechanism for energy exchange and therefore no process of thermalization before heat is conducted. If, on the other hand, one adds some interaction between the tracers, say, as suggested in [2], by assigning them a size, the resulting system will typically have transport equations for the diffusion of both mass and energy. Though these systems may be conceptually simple, they are usually mathematically too difficult to handle with the appropriate level or rigor.
An example of billiard with interacting tracer particles is the modified Lorentz gas with rotating discs proposed by Mejía-Monasterio et al. [7]. This system exhibits normal transport with non-trivial coupled equations for heat and mass transport. Though this model has attracted much attention among mathematicians [8,9,10,11,12,13], the rigorous proof of Fourier's law for such a system of arbitrary size and number of tracers relies on assumptions which, though they are plausible, are themselves not resolved. Moreover the equations which govern the collisions and the discs rotations do not, as far as we know, derive from a Hamiltonian description.
Yet a remedy to such limitations with respect to the number of particles in interaction that a rigorous treatment will handle may have been just around the corner [14], and, this, within a Hamiltonian framework. Indeed in [15], Bunimovich et al. introduced a class of dispersing billiard tables with particles that are in geometric confinement -i. e. trapped within cells-but that can nevertheless interact among particles belonging to neighbouring cells. The authors proved ergodicity and strong chaotic properties of such systems with arbitrary number of particles.
More recently, the idea that energy transport can be modeled as a slow diffusion process resulting from the coupling of fast energy-conserving dynamics has led to proofs of central limit theorems in the context of models of random walks and coupled maps which describe the diffusion of energy in a strongly chaotic, fast changing environment [16,17]. Although the extension of these results to symplectic coupled maps, let alone Hamiltonian flows, is not yet on the horizon, it is our belief that such systems as the dispersing billiard tables introduced in [15] will, if any, lend themselves to a fully rigorous treatment of heat transport within a Hamiltonian framework. As we announced in [18], the reason why we should be so hopeful is that the particles confinement has two important consequences : first, relaxation to local thermal equilibrium is preceded by a relaxation of individual particles to local equilibrium ‡, which occurs at constant energy within each cell, and has strong ergodic properties that guarantee the rapid decay of statistical correlations; and, second, heat transport, unlike in the rotating discs model, can be controlled by the mere geometry of the billiard, which also controls the absence of mass transport. As of the first property, the relaxation to a local equilibrium before energy exchanges take place is characterized by a fast time scale, much faster than that of relaxation to local thermal equilibrium among neighbouring cells, which is itself much faster than the hydrodynamic relaxation scale. There is therefore a hierarchy of three separate time scales in this system, the first accounting for relaxation at the microscopic scale of individual cells, the second one at the mesoscopic scale of neighbouring cells, and the third one at the macroscopic scale of the whole system.
In this paper, we achieve two important milestones towards a complete first principles derivation of the transport properties of such models. Having defined the model, we establish the conditions for separation of time scales and relaxation to local equilibrium, identifying a critical geometry where binary collisions become impossible. Assuming relaxation to local equilibrium holds, we go on to considering the time evolution of phase-space densities and derive, from it, a master equation which governs the exchange of energy in the system [19], thus going from a microscopic scale description of the Liouville equation to the mesoscopic scale at which energy transport takes place. We regard this as the first milestone, namely identifying the conditions under which one can rigorously reduce the level of descrition from the deterministic dynamics at the microscopic level to a stochastic process described by a master equation at the mesoscopic level of energy exchanges.
This master equation is then used to compute the frequency of binary collisions and to derive Fourier's law and the macroscopic heat equation, which results from the application of a small temperature gradient between neighbouring cells. This is our second milestone : an analytic formula for heat conductivity, exact for the stochastic system, and thus valid for the determinisitc system at the critical geometry limit. These results are then checked against numerical computations of these quantities, with outstanding agreement, and shown to extend beyond the critical geometry, with very good accuracy, to a wide range of parameter values.
We further characterize the chaotic properties of the model and offer arguments to account for the spectrum of Lyapunov exponents of the system, as well as the Kolmogorov-Sinai entropy, expressions which are exact at the critical geometry. Again, these results are very nicely confirmed by our numerical computations.
The paper is organized as follows. The models, which we coin lattice billiards, are introduced in section 2. Their main geometric properties are established, distinguishing transitions between insulating and conducting regimes under the tuning of a single parameter. The same parameter controls the time scales separation responsible for local equilibrium. Section 3 provides the derivation of the master equation which, under the assumption of local equilibrium, governs energy transport. The main observables are computed and their scaling properties discussed. In section 4 we review the properties of the model and assess the validity of the results of section 3 under the scope of our numerical computations. The chaotic properties of the models are discussed in section 5. We use simple theoretical arguments to predict some of these properties and compare them to numerical computations. Finally, conclusions are drawn in section 6.

Lattice billiards
To introduce our model, we start by considering the uniform motion of a point particle about a dispersing billiard table, B ρ , defined by the domain exterior to four overlapping discs of radii ρ, centered at the four corners of a square of sides l. The radius is thus restricted to the interval l/2 ≤ ρ < l/ √ 2, where the lower bound is the overlap (or bounded horizon) condition, and the upper bound is reached when B ρ is empty. Figure 1. Two equivalent representations of a dispersing billiard table : (left) a point particle moves uniformly inside the domain B ρ and performs specular collisions with its boundary; (right) a disc of radius ρ m moves uniformly inside the domain B ρ and performs specular collisions with fixed discs of radii ρ f = ρ − ρ m . The radius ρ m is a free parameter which is allowed to take any value between 0 and ρ.
As illustrated in figure 1, the motion of a point particle in this environment is a limiting case of a class of equivalent dispersing billiards, whereby the point particle becomes a moving disc with radius ρ m , 0 ≤ ρ m ≤ ρ, and bounces off fixed discs of radii ρ f = ρ − ρ m . In all these cases, the motion of the center of the moving disc is equivalent to that of the point particle in B ρ . The border of the domain B ρ constitutes the walls of the billiards.
In the absence of cusps (which occur at ρ = l/2), the ergodic and hyperbolic properties of these billiards are well established [20]. In particular, the long term statistics of the billiard map, which takes the particle from one collision event to the next, preserves the measure cos φ dr dφ, where φ denotes the angle that the particle post-collisional velocity makes with respect to the normal vector to the boundary. A direct consequence of this invariance is a general formula which relates the mean free path, ℓ, to the billiard table area, |B ρ |, and perimeter |∂B ρ |, ℓ = π|B ρ |/|∂B ρ |.
The ratio between the speed of the particle, which we denote v, and the mean free path gives the wall collision frequency §, ν c = v/ℓ, whose computation is shown in figure  2. With this quantity, one can relate the billiard map iterations to the time-continuous dynamics of the flow. In particular, the billiard map has two Lyapunov exponents, opposite in signs and equal in magnitudes, which, multiplied by the wall collision § Anticipating the more general definition of the wall collision frequency for interacting particle billiard cells, we adopt the subscript "c" in reference to "critical" for reasons to be clarified below.
frequency, correspond to the two non-zero Lyapunov exponents of the flow (there are two additional zero exponents related to the direction of the flow and conservation of energy). The results of numerical computations of the positive one, denoted λ + , are shown in figure 2 for different values of the parameters ρ. Thus let Q(i, j) denote the rhombus of sides l centered at point The rhombic billiard cell illustrated in figure 1 becomes a domain centered at point (c ij , d ij ), defined according to where δ[., .] denotes the usual Euclidean distance between two points, and (p k , q k ) = (±l/ √ 2, 0), (0, ±l/ √ 2) are the coordinates of the d ≡ 4 discs at the corners of the rhombus. Now consider a number of copies {B ρ (i, j)} i,j which tessellate a two-dimensional domain. We define a lattice billiard as a collection of billiard cells Each individual cell of this billiard table possesses a single moving particle of radius ρ m , 0 ≤ ρ m ≤ ρ and unit mass. All the moving particles are assumed to have independent initial coordinates within their respective cells, with the proviso that no overlap can occur between any pair of moving particles. The system energy, E = Nk B T (with N = n 1 × n 2 , the number of moving particles, T the system temperature, and k B Boltzmann's constant) is constant and assumed to be initially randomly divided among the kinetic energies of the moving particles, E = i,j ǫ ij , ǫ ij = mv 2 ij /2, where v ij denotes the speed of particle (i, j).
Energy exchanges occur when two moving particles located in neighbouring cells collide. Such events can take place provided the radii of the moving particles ρ m is large enough compared to ρ. Indeed the value of the critical radius, below which binary collisions do not occur, is determined by half the separation between the corners of two neighbouring cells, For the sake of illustration, the unlikely occurrence of a binary collision event at the critical radius ρ m = ρ c is shown in figure 3.
All the collisions are elastic and conserve energy, so that the dynamical system is Hamiltonian with 2N degrees of freedom. Its phase space of positions and velocities is 4N-dimensional. Accordingly, the sensitivity to initial conditions of the dynamics is characterized by 4N Lyapunov exponents, {λ i } 4N i=1 , obeying the pairing rule of symplectic systems, λ 4N −i+1 = −λ i , i = 1, . . . , 2N.
Collision events between two moving particles are referred to as binary collision events and will be distinguished from wall collision events, which occur between the moving particles and the walls of their respective confining cells. The occurrences of the former are characterized by a binary collision frequency, ν b , and the latter by a wall collision frequency, ν w . Both frequencies depend on the difference ρ m − ρ c , separating the moving particles radii from the critical radius, equation (4). By definition of ρ c , the binary collision frequency vanishes at ρ m = ρ c , ν b | ρm=ρc = 0, and, correspondingly, the wall collision frequency at the critical radius is the collision frequency of the singlecell billiard, ν w | ρm=ρc = ν c . We will assume from now, unless otherwise stated, that the system is globally isolated and apply periodic boundary conditions at the borders, thereby identifying B ρ (i + kn 1 , j + ln 2 ) with B ρ (i, j) for any k, l ∈ Z, Examples of such billiards are displayed in figure 4. Obviously the quincunx rhombic lattice structure, which is generated by the rhombic cells, is but one among . Examples of lattice billiards with triangular (top), rhombic (middle) and hexagonal (bottom) tilings. The coloured particles move among an array of fixed black discs. The radii of both fixed and mobile discs are chosen so that (i) every moving particle is geometrically confined to its own billiard cell (identified as the area delimited by the exterior intersection of the black circles around the fixed discs), but (ii) can nevertheless exchange energy with the moving particles in the neighbouring cells through binary collisions. The solid broken lines show the trajectories of the moving particles centers about their respective cells. The colours are coded according to the particles kinetic temperatures (from blue to red with increasing temperature). different possible structures. Triangular, upright square, or hexagonal cells can be used as alternative periodic structures. One might also cover the plane with random or quasi-crystalline tessellations. The only relevant assumptions in what follows is that the moving particles must be confined to their (dispersing) billiard cells and that binary interactions between neighbouring cells can be turned on and off by tuning the system parameters.
The two important features of such lattice billiards is that (i) there is no mass transport across the billiard cells since the moving particles are confined to their respective cells, and (ii) energy transport can occur through binary collision events which take place when the particles of two neighbouring cells come into contact. In periodic structures such as the quincunx rhombic lattice, the possibility of such collisions is controlled by tuning the parameter ρ m about the critical radius ρ c , keeping ρ fixed.
We can therefore distinguish two separate regimes : • Insulating billiard cells: 0 ≤ ρ m < ρ c Absence of interaction between the moving discs. No transport process across the individual cells can happen; • Conducting billiard cells: ρ c < ρ m < ρ Binary collision events are possible. Energy transport across the individual cells takes place.
The case ρ m = ρ c is singular. We will refer to the critical geometry as the limit ρ m > → ρ c . In the insulating regime, there is no interaction among moving particles so that the billiard cells are decoupled. The moving particles are independent and their kinetic energies are individually conserved, resulting in 2N zero Lyapunov exponents. The equilibrium measure in turn has a product structure and phase-space distributions are locally uniform with respect to the particles positions and velocity directions. The N positive Lyapunov exponents of the system are all equal to the positive Lyapunov exponent of the single-cell dispersing billiard, up to a factor corresponding to the particles speeds, v ij = 2ǫ ij /m : λ ij + = v ij λ + , where λ + is the Lyapunov exponent of the single-cell billiard measured per unit length.
When particles are allowed to interact, on the other hand, local energies are exchanged through collision events. Thus only the total energy is conserved in the conducting regime. The ergodicity of such systems of geometrically confined particles in interaction was proven by Bunimovich et al. [15]. The resulting dynamical system, whose equilibrium measure is the microcanonical one (taking into consideration that particles are otherwise uniformly distributed within their respective cells), enjoys the K-property. This implies ergodicity, mixing, and strong chaotic properties, including the positivity of the Kolmogorov-Sinai entropy. Two Lyapunov exponents are zero, one associated to the conservation of energy, the other to the direction of the flow. The 2(2N − 1) remaining Lyapunov exponents form non-vanishing pairs of exponents with opposite signs, The regime of interest to us is that corresponding to particles interacting rarely, which is to say, in analogy with a solid, that particles mostly vibrate inside their cells, ignorant of each other, and only seldom making collisions with their neighbours, thereby exchanging energy. As we turn on the interaction and let ρ m ρ c , binary collisions, though they can occur, will remain unlikely. This is to say that the binary collision frequency, ν b , will, in this regime, remain small with respect to the wall collision frequency, ν w , which in the absence of interaction and, in particular, at the critical geometry, we recall is equal to the wall collision frequency of the single-cell billiard, ν w | ρm=ρc = ν c . When ρ m ρ c , we therefore expect ν w ≫ ν b , as well as ν w ≃ ν c . In words : time scales separate. The consequence is that relaxation to local equilibrium -i. e. uniformization of the distribution of the particles positions and velocity directions at fixed speeds-occurs typically much faster than the energy exchange which drives the relaxation to the global equilibrium. This mechanism justifies resorting to kinetic theory in order to compute the transport properties of the model.

From Liouville's equation to the master equation
The phase-space probability density is specified by the N-particles distribution function where r a and v a , a = 1, . . . , N, denote the ath particle position and velocity vectors. The index a stands for the label (i, j) of the cells defined by equation (2). For our system, as is customary for hard sphere dynamics, this distribution satisfies a pseudo-Liouville equation [21], which is well defined despite the singularity of the hard-core interactions. This equation, which describes the time evolution of p N is composed of three types of terms: (i) the advection terms, which account for the displacement of the moving particles within their respective billiard cells; (ii) the wall collision terms, which account for the wall collision events, between the moving particles and the d fixed scattering discs which form the cells walls; and (iii) the binary collision terms, which account for binary collision events, between moving particles belonging to neighbouring billiard cells : Each wall collision term involves a single moving particle with index a and one of the d fixed discs in the corresponding cell, with index k and position R k . Let r ak = r a − R k denote their relative position. Following [22], we have whereê denotes the normal unit vector to the fixed disc k in the cell of particle a.
Likewise the binary collision operator, written in terms of the relative positions r ab and velocities v ab of particles a and b, and the unit vectorê ab that connects them, is We notice that only the terms B (a,b) corresponding to first neighbours are non-vanishing and contribute to the double sum on the RHS of equation (5).
Provided we have a separation of time scales between wall and binary collisions, the advection and wall collision terms on the RHS of equation (5) will typically dominate the dynamics on the short time ∼ 1/ν w , which follows every binary collision event, thus ensuring, thanks to the mixing within individual billiard cells, the relaxation of the phase-space distribution p N to local equilibrium well before the occurrence of the next binary event, whose time scale is quickly relaxes to a locally uniform distribution, which depends only on the local energies, justifying the introduction of where v a ≡ |v a |. On the time scale of binary collision events, this distribution subsequently relaxes to the global microcanonical equilibrium distribution. This process accounts for the transport of energy, and can be characterized by the master equation [19] ∂ t P where W (ǫ a , ǫ b |ǫ a − η, ǫ b + η) denotes the probability that an energy η be transferred from particle a to particle b as the result of a binary collision event between them. This equation is a closure for the local equilibrium distribution P (leq) N , obtained from equation (5) under the assumption that ν w ≫ ν b . The first two terms on the RHS of equation (5) are eliminated because they leave invariant the local distribution N a=1 δ(ǫ a − mv 2 a /2). There remain the contributions (7) from the binary collisions, which, under the assumption that the local distibutions are uniform with respect to the positions and velocity directions, yield the following expression of W : where the first integration is performed over the positions of the center of mass, R ≡ (r a + r b )/2, between the two particles a and b, given that they are in contact and both located in their respective cells, and over the angle φ of the unit vector connecting a and b,ê ab = (cos φ, sin φ). The normalizing factor |L ρ,ρm (2)| denotes the 4-volume of the billiard corresponding to two neighbouring cells a and b, which, with the assumption that ρ m ρ c , can be approximated by |L ρ,ρm (2)| ≃ |B ρ | 2 . This substitution amounts to neglecting the overlap between the two particles; see equation (23) for a refinement of that approximation. We point out that the position and velocity integrations in equation (10) can be formally decoupled; in this way, we can prove that the transition rate W is given in terms of Jacobian elliptic functions, see Appendix A.

Geometric factor
As we show in Appendix A, an important property of the master equation (9) is that the factor which accounts for the geometry of collision events factorizes from the part of the kernel that accounts for energy exchanges. Therefore, as ρ m → ρ c , the critical value of the radius at which binary collision events become impossible, which is the regime where the billiard properties are accurately described in terms of the master equation above, the geometric factor dφdR encloses the scaling properties of observables with respect to the billiard geometry. We now compute this quantity.
A binary collision occurs when particles a and b come to a distance 2ρ m of each other, with r a ∈ B ρ (a) and r b ∈ B ρ (b). Let R = (x, y) be the center of mass coordinates and φ be the angle between the particles relative position and the axis connecting the center of the cells. Taking a reference frame centered between the cells, we may write where σ i = ±1 and i = a or b. The integral to be evaluated is the volume of the triplets (x, y, φ) about the origin so that As illustrated in figure 5, for different orientations φ of the vector connecting the two particles, this is a region bounded by four arc-circles, which we denote by where σ, τ = ±1. As seen from figure 5, the area is connected for −φ T ≤ φ ≤ φ T , where φ T is the angle φ at which opposite arcs intersect, Beyond that value, the area splits into two triangular areas. These areas shrink to zero at the angle φ given by Figure 5. Possible positions of the center of mass (x, y) for different values of φ, see equation (11), at a binary collision event. Here ρ = 13/25l and ρ m = 13/50l.
We denote by x 1 < x 2 < x 3 < x 4 the four corners of the rectangular domain, For φ ≥ φ T , the points at which the opposite arcs intersect are given by Combining equations (13)- (17) together, we can make use of the symmetry φ → −φ and write the integral to be computed as where which is the area bounded by the four arcs y σ,τ , and is the area of the overlapping opposite arcs y −1,+1 and y +1,−1 , which occurs when The computation of these expressions is easily performed numerically, and the result shown in figure 6.
Near the critical geometry, we expand the quantity (18) in powers of the difference ρ m − ρ c , The first two coefficients vanish, so that the leading term corresponds to n = 3. The first few coefficients are derived in Appendix B and given by : We further note that the computation of the two-cell 4-volume, |L ρ,ρm (2)|, which, as noticed above, is approximated by the square of the single cell area |B ρ | 2 , can be improved using equation (21). Indeed it is easily seen that d|L ρ,ρm (2)|/dρ m = −2α(ρ, ρ m ), which implies It is immediate to check that corrections |L ρ,ρm

Binary collision frequency
Having computed the transition rates of the master equation with respect to the billiard geometry, we now turn to the computation of observables. The first quantity of interest, which can be readily computed from equation (10), is the binary collision frequency ν b . This is an equilibrium quantity which, in the (global) microcanonical ensemble with energy E = ǫ 1 + . . . + ǫ N , involves the two-particle energy distribution, and can be written as Taking the large N limit and letting E = Nk B T ≡ N/β, we can write P (eq) 2 ]. Substituting this expression into the above equation and inserting the expression of W from equation (10), we obtain, after some calculations, This expression involves the geometric factor (18) in the first bracket. The leading term in the second bracket is the canonical expression of the binary collision frequency. The second term is a positive finite N correction, which is useful in that it shows that the binary collision frequency decreases to its asymptotic value as N → ∞.

Rescaled master equation
The binary collision frequency (26) defines a natural dimensionless time scale for the stochastic process described by the master equation (9) with transition rates (10). The master equation can thus be converted to dimensionless form by rescaling the energies by a reference thermal energy and time by the corresponding asymptotic (N → ∞) value of the binary collision frequency (26). Introducing the variables with the transition rates w(e a , e b |e a − h, e b + h) = 2 π 3 2 ) This master equation shows that all the properties of heat conduction are rescaled by the binary collision frequency and temperature in its limit of validity where the collision frequency vanishes. In particular, this shows that the coefficient of heat conduction is proportional to the binary collision frequency in this limit, as explained in the next subsection.

Thermal conductivity
Starting from the master equation (9), we derive an equation for the evolution of the kinetic energy of each moving particle, ǫ a ≡ k B T a , which defines the local temperature, where . denotes an average with respect to the energy distributions. By the structure of equation (9), such an equation can be expressed in terms of the transfer of energy due to the binary collisions between neighbouring cells. The time evolution of the average local kinetic energy is given by with the energy flux defined as This expresses the local conservation of energy. Over long time scales, the probability distribution becomes controlled by this local conservation of energy, the slowest variables being the local kinetic energies ǫ a or, equivalently, the local temperatures T a as defined above. This holds even though statistical correlations develop between the local energies in the probability distributions. These statistical correlations are well known for transport processes ruled by master equations such as Eq. (9) [19,23] and are observed in the present system as well.
To be specific, we consider a one-dimensional chain, extending along the x-axis, formed with a succession of pairs of rhombic billiard cells arranged in quincunx, similarly to the middle panel of figure 4, except the vertical height is here only one unit of length. The unit of horizontal and vertical lengths is thus l √ 2 and there are two cells per each unit of length. Similar results hold for different choices of geometry modulo straightforward adaptations.
We imagine that the system is in a non-equilibrium state, with a small temperature difference δT about an average temperature T between neighbouring cells, and consider the average heat transfered from cell a at inverse temperature β a = β + δβ/2 to cell b at inverse temperature β b = β − δβ/2, δβ = −δT /(k B T 2 ), both cells being assumed to be in thermal equilibrium at their respective temperatures. The statistical correlations we observe in the present system are of the order of δβ 2 , as it is the case in other systems [19]. In the non-equilibrium state, these statistical correlations are controlled in the long-time limit by the local temperatures. Since the process is here ruled by a Markovian master equation in the limit of small binary collision frequency, we get the equation of heat for the temperature where the local temperature is here written as T (x, t) = ǫ ij /k B , x = √ 2l(i + j/2), i = 1, . . . , N/2, j = 0, 1 [see equation (1)].
According to the rescaling property of the master equation discussed in section 3.4, the heat conductivity is proportional to the binary collision frequency : with a dimensionless constant A. An analytical estimation can be obtained by transforming the master equation into a hierarchy of equations for all the moments of the probability distribution: ǫ a , ǫ a ǫ b , ǫ a ǫ b ǫ c ,... The evolution equations of these moments are coupled.
Truncating the hierarchy at the equations for the averages ǫ a , we get the approximate heat conductivity : So that A = 1 in this approximation. The same result holds if we include the equations for the moments ǫ a ǫ b with |a − b| ≤ 1. Though the approximate result (36) above does not rule out possible corrections to A = 1, we argue that A = 1 is an exact property of the master equation (9) in the infinite system limit. This claim is borne out by extensive studies of the stochastic process described by equation (9) that will be reported in a separate publication [24]. The focus is here on the billiard systems whose conductivity may however bear corrections to this identity, due we believe to lack of sufficient separation of time scales between wall and binary collision events. In the following, the results of numerical computations are presented which support these claims.

Numerical results
The above formulae for the binary collision frequency, equation (26), and the thermal conductivity, equation (36), together with the expressions of the geometric factors, equations (22) and (23), provide a detailed picture of the mechanism which governs the transport of heat in our model. Numerical computations of these quantities further add to this picture and provide strong evidence of the validity of our theoretical approach.

Binary and wall collision frequencies
For the sake of computing the binary collision frequency, we simulate the quasi-one dimensional channel of N cells with rhombic shapes and apply periodic boundary conditions at the horizontal ends of the channel. Thus each cell has two neighbouring We notice that the heat capacity per particle is equal to c V = (1/N )∂E/∂T = k B , so that the thermal conductivity is also equal to the thermal diffusivity in units where k B = 1.

cells, left and right, and interactions between any two neighbouring cells can occur through both top and bottom corners ¶.
In figure 7, we compare the computations of ν b to ν w for a system of N = 10 cells at unit temperature. The parameters are taken to be ρ = 9/25 and ρ m = 3/25, . . . , 17/50 by steps of 1/50. Both collision frequencies are computed in the units of the microcanonical average velocity, v N ≡ 2 N N/2(N − 1)!/(2N − 1)!!, which, as N → ∞, converges to the canonical average velocity v ∞ = π/2. The wall collision frequency is compared to the collision frequency of the isolated cells ν c , itself measured in the units of the single particle velocity. As expected, ν b /ν w ≪ 1 and ν w ≃ ν c for ρ m ρ c . The crossover ν b ≃ ν w occurs at ρ m ≃ 11/50 for this value of ρ.  The thermal conductivity can be obtained by computing the heat flux in a nonequilibrium stationary state. Such stationary states occur when the two ends of the channel are put in contact with heat baths at separate temperatures, say T − and T + , T − < T + . Let a system of N cells be in contact with two thermostated cells at respective temperatures T ± , and let these cell indices be n = ±(N + 1)/2 (we take N odd for the sake of definiteness). Provided the difference between the bath temperatures is small, T + −T − ≪ (T + +T − )/2, a linear gradient of temperature establishes through the system, with local temperatures (37) ¶ We mention that this set-up allows for re-collision between two particles (under stringent conditions), due to the vertical periodic boundary conditions. This conflicts with the assumptions of [15], but does not seem to affect the results as far as numerics are concerned.
Under these conditions, the nonequilibrium stationary state is expected to be locally well approximated by a canonical equilibrium at temperature T n .
This can be checked numerically. In fact the local thermal equilibrium is verified (by comparing the moments of ǫ n to their Gaussian expectation values) under the weaker property of small local temperature gradients, i. e. (T + − T − )/N ≪ (T + + T − )/2, for which the temperature profile is generally not linear since the thermal conductivity depends on the temperature. Indeed, since κ ∝ T 1/2 , we expect in that case, according to Fourier's law, the profile The thermal conductivity can therefore be computed from the heat exchanges of the chain in contact with the two cells at the ends of the chain, respectively thermalized at temperatures T ± δT /2, δT ≪ T . The thermalization of the end cells is achieved by randomizing the velocities of the two particles at every collision they make with their cells walls, according to the usual thermalization procedure of particles colliding with thermalized walls [28].
First, we consider a chain containing a single cell in order to test the validity of the master equation. In this case, the procedure amounts to simulating a single particle confined to its cell and performing random collisions with stochastic particles which penetrate the cell corners according to the statistics of binary collisions. For a chain with a single cell, the heat conductance is given by Eq. (36) with A = 1, as there are no correlations with the stochastic particles. Figure 8 shows the results of the computations of the heat conductance + with this method and provides a comparison with the binary collision frequency ν b on the one hand (left panel), as well as with the results of our kinetic theory predictions on the other hand (right panel).
The agreement between the data and equations (26), (36), (35) and the computation of the integrals (18)- (20), especially as ρ m → ρ c , demonstrates the validity of the stochastic description of the billiard system, equation (9).
Next, we increase the size of the chain in order to reach the thermal conductivity in the limit of an arbitrarily large chain. The results are that statistical correlations appear between the kinetic energies along the chain. As we show below, their influence on the computed value of the conductivity diminishes as ρ m → ρ c .
Fix T − = 0.5 and T + = 1.5 to be the baths temperatures, and let the size of the system increase from N = 1 to N = 20 as ρ m is progressively decreased from ρ m = 11/50 to ρ m = 13/100, with fixed ρ = 9/25. As one can see from figure 8, this range of values of ρ m crosses over from a regime where the separation of time scales is not effective + Here and in the sequel, the thermal conductance or conductivity are further divided by l 2 √ T , where l = 1/ √ 2 is the rhombic cell size, so as to eliminate its length and temperature dependences, thus defining the reduced thermal conductivity κ * = κ/(l 2 √ T ). In these expressions and from here on, we further set k B ≡ 1.  to one where it appears to be and where the stochastic model should therefore be a reasonable approximation to the process of energy transport in the billiard.
For all values of ρ m , we measured the temperature profile and heat fluxes throughout the system and inferred the value of the reduced thermal conductivity by linearly extrapolating the ratios between the average heat flux and local temperature gradient divided by the square root of the local temperature as functions of 1/N to the vertical axis intercept, corresponding to N = ∞. Given the parameter values, every realisation was carried out over a time corresponding to 1,000 interactions between the system and baths and repeated over 10 4 realisations. For N up to 20, this time provided satisfying stationary statistics, with temperature profiles verifying equation (38) and statistically constant heat fluxes.
The results of the linear regression used to compute the value of κ/ν b for selective values of ρ m are shown in figures 9 and 10.
Our results thus make it plausible that the ratio between the thermal conductivity and binary collision frequency approaches unity as the parameter ρ m decreases towards the critical value ρ c . As will be show in a separate publication [24], this is indeed a property of the stochastic model described by equation (9) and has been verified numerically by direct simulation of the master equation within an accuracy of 4 digits. As of the billiard system, it is unfortunately difficult to improve the results beyond those presented here as the CPU times necessary to either increase N or decrease ρ m quickly become prohibitive. Nevertheless, the data displayed in figures 9 and 10 offer convincing evidence that the thermal conductivity is well approximated by the binary collision frequency so long as the separation of time scales between wall and binary collision events is effective.

Helfand moment.
The computation of the thermal conductivity can also be performed in the global equilibrium microcanonical ensemble using the method of Helfand moments [25,26,27].
The Helfand moment has expression H(t) = a x a (t)ǫ a (t), where x a (t) denotes the horizontal position of particle a at time t and ǫ a (t) = |v a (t)| 2 /2 its kinetic energy (the masses are taken to be unity). The computation of the time evolution of this quantity proceeds by discrete steps, integrating the Helfand moment from one collision event to the next, whether between a particle and the walls of its cell, or between two particles. Let {τ n } n∈Z denote the times at successive collision events. In the absence of binary collisions, the energies are locally conserved and the Helfand moment changes according to H(τ n ) = H(τ n−1 ) + a [x a (τ n ) − x a (τ n−1 )]ǫ a (τ n−1 ). If, on the other hand, a binary collision occurs between particles k and l, the Helfand moment changes by an additional Computing the time average of the squared Helfand moment, we obtain an expression of the thermal conductivity according to where L = N/2 is the horizontal length of the system. Figure 11 shows the results of a computation of the thermal conductivity through equation (39) for different system sizes. Though the actual values of κ vary wildly with N, it is clear that a finite asymptotic value is reached for N ≃ 10 2 . In this case, the constant of proportionality in equation (35) takes the value A = 0.98 ± 0.08, close to 1. Similar results were obtained for other parameter values, and other cell geometries as well.

Lyapunov spectrum
A key aspect of our model, which justifies the assumption of local equilibrium, is that it is strongly chaotic. This property can be illustrated through the computation of the Lyapunov spectrum and Kolmogorov-Sinai entropy in equilibrium conditions.
As mentioned earlier, in the absence of interaction between the cells, ρ m < ρ c , The Lyapunov spectrum of a system of N cells has N positive and N negative Lyapunov exponents, which, if divided by the average speed of the particle to which they are attached, are all equal in absolute value. This reference value we denote by λ + . The 2N remaining Lyapunov exponents vanish.
As we increase ρ m and let the particles interact, we expect that, in the regime 0 < ρ m − ρ c ≪ 1, where binary collision events are rare, the Lyapunov exponents will essentially be determined by λ + multiplied by a factor which is specified by the particle velocities. The exchange of velocities thus produces an ordering of the exponents which can be computed as shown below. We note that the other half of the spectrum, which remains zero in this approximation, will only pick up positive values as a result of the interactions.
Assume for the sake of the argument that N is large. The probability that a given particle with velocity v has exponent λ = vλ + less than a value λ i = v i λ + can be approximated by the probability that the particle velocity be less than v i , which, if we assume a canonical form of the equilibrium distribution, is But this probability is simply (N − i + 1/2)/N. Therefore the half of the positive Lyapunov exponent spectrum, which is associated to the isolated motion of particles within their cells, becomes, in the presence of rare collision events, with ordering λ 1 > λ 2 > . . . > λ N . In particular, the largest exponent λ 1 grows like √ ln N . The Kolmogorov-Sinai entropy on the other hand is extensive : h KS = Nλ + π/(2mβ).
Refined expressions can be computed by taking the microcanonical distribution associated to a finite N. In particular, the expressions of the Lyapunov exponents become We mention in passing that similar arguments are relevant and can be used to approximate the Lyapunov spectrum (actually half of it) of other models, such as a mixture of light and heavy particles [30]. Let us again insist that equations (41) and (42) account for only N of the 2N − 1 positive Lyapunov exponents. N − 1 are vanishing within this approximation. For all these Lyapunov exponents, we expect the corrections to vanish with the binary collision frequency as we go to the critical geometry. Figures 12 and 13 show the results of a numerical computation of the whole spectrum of Lyapunov exponents and corresponding Kolmogorov-Sinai entropy for the one-dimensional channel of rhombic cells. The agreement with equation (42) as ρ m → ρ c is very good, at the exception perhaps of the last few among the first N exponents, whose convergence to the asymptotic value (42) appears to be slower. Interestingly, the largest exponents have a minimum which occurs at about the value of ρ m for which the binary and wall collision frequencies have ratio unity, see figure 7. Indeed, for larger radii ρ m , the spectrum is similar to that of a channel of hard discs (without obstacles) [29]. Note that the same holds for the ratio between the thermal conductivity and binary collision frequency, as seen from figure 8. Thus one can interpret the occurrence of a minimum of the largest exponent as evidence of a crossover from a near close-packing solid-like phase (ρ m ρ) to a gaseous-like phase trapped in a rigid structure (ρ m ρ c ).

Conclusions
To summarize, lattice billiards form the simplest class of Hamiltonian models for which one can observe normal transport of energy, consistent with Fourier's law. Geometric confinement restricts the transport properties of this system to heat conductivity alone, thereby avoiding the complications of coupling mass and heat transports, which are common to other many particle billiards.
The strong chaotic properties of the isolated billiard cells warrant, in a parametric regime where interactions among moving particles are seldom, the property of local equilibrium. This is to say, assuming that wall collision frequencies are order-ofmagnitude larger than binary collision frequencies, that one is entitled to making a Markovian approximation according to which phase-space distributions are spread over individual cells. We are thus allowed to ignore the details of the distribution at the level of individual cells and coarse-grain the phase-space distributions to a manyparticle energy distribution function, thereby going from the pseudo-Liouville equation, governing the microscopic statistical evolution, to the master equation, which accounts for the energy exchanges at a mesoscopic cell-size scale and local thermalization. The energy exchange process further drives the relaxation of the whole system to global equilibrium.
This separation of scales, from the cell scale dynamics, corresponding to the microscopic level, to the energy exchange among neighbouring cells at the mesoscopic level, and to the relaxation of the system to thermal equilibrium at the macroscopic level, is characterized by three different rates. The process of relaxation to local equilibrium has a rate given by the wall collision frequency, much larger than the rate of binary collisions, which characterizes the rate of energy exchanges which accompany the relaxation to local thermal equilibrium, itself much larger than the hydrodynamic relaxation rate, given by the binary collision rate divided by the square of the macroscopic length of the system.
On this basis, having reduced the deterministic dynamics of the many particles motions to a stochastic process of energy exchanges between neighboring cells, we are able to derive Fourier's law and the macroscopic heat equation.
The energy transport master equation can be solved with the result, to be presented elsewhere [24], that the binary collision frequency and heat conductivity are equal. Under the assumption that the wall and binary collision time scales of the billiard are well separated, the transposition of this result to the billiard dynamics is that : (i) The heat conductivity of the mechanical model is proportional to the binary collision frequency, i. e. the rate of collisions among neighbouring particles, κ with a constant A that is exactly 1 at the critical geometry, ρ m → ρ c , where the evolution of probability densities is rigorously described by the master equation, and remains close to unity over a large range of parameter values for which we conclude the time scale separation is effective and the master equation therefore gives a good approximation to the energy transport process of the billiard ; (ii) The heat conductivity and the binary collision frequency both vanish in the limit of insulating system, ρ m → ρ c , with (ρ m − ρ c ) 3 , where the coefficient c 3 depends on the specific geometry of the binary collisions.
Though both results are exact strictly speaking only in the limit ρ m → ρ c , the first one, according to our numerical computations, is robust and holds throughout the range of parameters for which the binary collision frequency is much less than the wall collision frequency, ν b ≪ ν w . The deviations from A = 1 which we observed at intermediary values of ρ m , where the separation of time scales is less effective, are interpreted as actual deviations of the energy transport process of the billiard from that described by the master equation, where correlations between the motions of neighboring particles must be accounted for.
Under the conditions of local equilibrium, the Lyapunov spectrum has a simple structure, half of it being determined according to random velocity distributions within the microcanonical ensemble, while the other half remains close to zero. The analytic expression of the Lyapunov spectrum that we obtained is thus exact at the critical geometry. The Kolmogorov-Sinai entropy is equal to the sum of the positive Lyapunov exponents and thus determined by half of them. It is extensive in the number of particles in the system, whereas the largest Lyapunov exponent grows like the square root of the logarithm of that number. As we mentioned earlier, we believe our method is relevant to the computation of the Lyapunov spectrum of other models of interacting particles [30]. The computation of the full spectrum, particularly regarding the effect of binary collisions on the exponents, remains an open problem. angle integration is performed over the domain D + such that √ ǫ a cos θ a > √ ǫ b cos θ b .
With the above expression (A.1), the explicit φ dependence has disappeared so that we have effectively decoupled the velocity integration from the integration over the direction of the relative position between the two colliding particles. We can further transform this expression in terms of Jacobian elliptic functions as follows.
Let x i = cos θ i , i = a, b, in equation (A.1), which becomes where θ(.) is the Heaviside step function. We thus have to perform the x a and x b integrations along the line defined by the argument of the delta function, and that satisfies the condition To carry out this computation, we have to consider the following alternatives : where K denotes the Jabobian elliptic function of the first kind, Thus the kernel is, in this case, This case is similar to case (i), with equation (A.5) replaced by and −1 < x a < +1. The expression of the kernel corresponding to this case is therefore This case is similar to case (ii), with − (ǫ b + η)/ǫ a < x a < + (ǫ b + η)/ǫ a . In this case, the expression of the kernel is given by W (ǫ a , ǫ b |ǫ a − η, ǫ b + η) = 2ρ m π 2 |L ρ,ρm (2)| 2 mǫ a K ǫ b + η ǫ a dφdR. (A.11) The cases with ǫ a > ǫ b are obtained from the cases above with the roles of a and b interchanged and η → −η.