Article Deterministic Thermal Reservoirs OPEN ACCESS

We explore the consequences of a deterministic microscopic thermostat-reservoir contact mechanism for hard disks where the collision rule at the boundary is modified. Numerical evidence and theoretical argument is given that suggests that an energy balance is achieved for a system of hard disks in contact with two reservoirs at equal temperatures. This system however produces entropy near the the system-reservoir boundaries and this entropy flows into the two reservoirs. Thus rather than producing an equilibrium state, the system is at a steady state with a steady entropy flow without any associated energy flux. The microscopic mechanisms associated with energy and entropy fluxes for this system are examined in detail.


Introduction
The first concept established in equilibrium statistical mechanics is the existence of an equilibrium state; that is, after a sufficient time macroscopic observables have steady time independent average values.The concepts of temperature and entropy as applied in equilibrium thermodynamics cannot apply to some transient system and they do not easily generalize to nonequilibrium systems.However, it is possible that nonequilibrium steady states may have a thermodynamic description, particularly those that are close to equilibrium.For an isolated equilibrium system where the entropy is a maximum, if the system is connected to one or more reservoirs and perturbed from equilibrium the system is said to produce entropy.The entropy production is transferred to a reservoir or surroundings and the intrinsic entropy of the nonequilibrium system is lower than at equilibrium.Thus, when the perturbation is removed, the relaxation is driven by an increase in the intrinsic entropy.Perhaps a more useful construct is to instead consider the entropy flow as a flow of energy through the system that maintains it in a steady state.This energy flow has been called the housekeeping heat [1] to distinguish it from other intrinsic heat flows.
In the canonical ensemble a system is in thermal contact with a reservoir which can exchange energy with the system but is sufficiently large so that its thermodynamic state is not effected.Here it is assumed that the system comes to an equilibrium state and the details of the thermal contact are not believed to be important, however, away from equilibrium the details of the thermal contact are central.We are only interested in exchanging energy to equilibrate a system, so we assign a fixed temperature or reservoir momentum to each reservoir.This is an external parameter that remains fixed, so effectively the reservoir's thermodynamic state cannot change through the process of interaction with the system.
From a macroscopic viewpoint, the concept of temperature at equilibrium is well established by the zeroth law of thermodynamics.This involves the idea of bringing together two systems and allowing them to exchange energy.However, as soon as this idea is considered on a microscopic level it is clear that a detailed mechanism that facilitates energy transfer is needed, and that involves the details of the microscopic thermal contact interaction at the boundary.For systems away from equilibrium the status of the zeroth law is not clear, and despite a number of studies of possible nonequilibrium temperatures, the question is not settled [2][3][4][5][6][7][8][9].There is a large literature on thermal conductivity in which reservoirs are considered as stochastic sources of thermalized particles [10] at some required temperature, but this approach precludes the use of modern dynamical system techniques.
Molecular dynamics simulations have proved to be a very effective means of testing theoretical approaches to the study of fluids both in equilibrium and nonequilibrium steady states [11].Given a particular atomic pair interaction, the results are free of approximations and have an accuracy limited only by statistical considerations.It is usual to use the equipartition theorem to define the kinetic temperature, so in a system of N particles in d spatial dimensions, the translational kinetic energy is kT 2 per degree of freedom, giving In this case we will extend the idea of a kinetic temperature to define the local kinetic temperature as the temperature of the single atom in a particular volume as T i = p 2 i /dm, and this will be used to determine the local temperature profile inside the system and to test equilibration.We identify separately the x and y components of the temperature as and we will use the difference between the average of the x and y components of the local kinetic temperature to give a measure of the deviation from local thermodynamic equilibrium.A change in local temperature near the system-reservoir boundary for a mechanical reservoir [12] termed a Kapitza resistance has been previously reported [10,13], but here we investigate the origin of this effect in our system.

The Generalized Gibbs Relation: Kinetic Theory
In nonequilibrium thermodynamics the kinetic-theory basis of the heat transport model in quasi-one-dimensional (QOD) systems has been discussed recently [14].The basic ingredient [15,16] is the Boltzmann entropy S(t) which is defined, up to a constant, to be where s(r, t) is the entropy density at position r at time t and we take Boltzmann's constant to equal one.The Boltzmann entropy is just the kinetic contribution and the configurational contribution written in terms of correlation functions is often referred to as the Green entropy [17].We assume that the system has reached a steady state so that the entropy, and all other thermodynamic properties, no longer depend on time.The distribution function f (r, v, t) = n(r, t)P (r, v, t) obeys the Boltzmann equation, which without external forces takes the following form where J[f ] is the collision integral and P (r, v, t) is the momentum probability density.The distribution function is normalized so that dvf (r, v, t) = n(r, t) where n(r, t) is the local number density of the system.The local entropy-balance equation is obtained by using the Boltzmann equation, Equation (4), to calculate the total time-derivative of the Boltzmann entropy, Equation (3).When the fluid is on average stationary, the result is expressed as where j s is the entropy flux and σ ent is the entropy production.The kinetic contribution to the entropy flux j s is given by j s (r) = − dv vf ln f The entropy-production per unit volume σ ent is given exactly by the first equality in Equation ( 8).However, the collision integral is too difficult to handle so we approximate it using the BGK approximation as ), thus we can estimate σ ent as where ν(r) is the local relaxation frequency proportional to ν 0 n(r) T (r) (ν 0 is a constant that depends on dimensionality).According to the H-theorem [18] σ ent (r) must be positive.All the local quantities associated with the entropy balance equation have been defined through the Boltzmann equation and their definitions are exact within its validity and are not restricted to the linear regime.In the Chapman-Enskog expansion solution to first order in the gradient (that is Navier-Stokes order [19]) the distribution function is given by f = f loc (1 + Φ), where f loc is the local-equilibrium distribution function, which is a Gaussian at the local temperature and number density.

The Model System
To specify the thermal contact between system and reservoir, it is necessary to move from a macroscopic picture to a microscopic one, describing the collision process that allows energy to flow at the system-reservoir boundary.We can separate simplified reservoir models into two classes; stochastic and deterministic.A typical stochastic model uses Gaussian reservoirs to provide a source of thermalized particles from a Gaussian distribution at a set temperature.These reservoirs cannot sense the temperature as they are simply sources of thermalized particles.
The QOD system introduced recently [20] can be modified to interact with an idealized mechanical heat reservoir in a deterministic and reversible way [21], to study both heat conduction in low dimensional systems and the Lyapunov mode structure.The deterministic reservoir allows the calculation of the usual dynamical systems properties as well as the thermodynamic properties.Similar systems have been used recently to study heat flow [22,23], Eckmann and Young [24] have developed Hamiltonian and stochastic models for heat transport in low dimensional systems.The Hamiltonian models consist of energy storage devices which couple to each other through the motion of tracer particles that carry the energy.Although these models principally store energy as rotational kinetic energy, they are similar in many ways to the QOD system considered here which store energy locally in the y degree of freedom.If we envisage taking the QOD system, expanding odd numbered disks to twice their diameter and shrinking the even numbered ones to points, we obtain a scatterer-tracer system but with tracers always confined between the same two scatterers.The chaotic scattering of hard disks in two-dimensions gives a mechanism to transfer energy between particles with different momentum components, and hence may be more realistic than one-dimensional models particularly for fluid systems.

System Dynamics
The QOD system contains N hard disks of diameter σ = 1 and mass m = 1 in a narrow channel L x × L y where the system width is less than twice the particle diameter L y < 2σ so that the disks remain ordered (see Figure 1).As the position of the i th particle is restricted to lie between particles i − 1 and i + 1 any property of particle i can be associated with the same local density at the average position x i of particle i.For example, the local number density n(r, t) is the inverse of the average volume occupied by the particle The equations of motion connecting the QOD system to the two reservoirs at x = 0 and L x define the thermal contact.When a particle collides with a reservoir wall the usual collision rule in the x direction is replaced by where p I = ± √ 2T I is the fixed value of the reservoir momentum for the I th reservoir (either p L for the left-hand side or p R for the right-hand side) for a given boundary temperature of T I .Note that the reservoir momentum is always directed from the reservoir into the system; so here p L > 0 and p R < 0. We use a system width of L y = 1.15σ, allowing L x to vary with the number of particles N to give the desired density ρ = N σ 2 /(L x L y ).The temperature profile is determined from the average components of the kinetic temperature of each particle [the average of Equation ( 2  The height L y is sufficiently small that the disks cannot pass one another.We choose the coordinate origin to be located at the bottom left corner of the system, and system boundaries are denoted by dashed lines.The boundaries at x = 0 and L x are hard walls and those at y = 0, L y periodic (that is, a (H,P) quasi-one-dimensional system).
The parameter represents the strength of the coupling between the reservoir and the system.If = 0 there is no interaction with the reservoir, while if = 1 the incoming momentum is completely replaced by the reservoir momentum.Here we use an intermediate value of = 0.5, which provides an effective mix of the incoming momentum with the reservoir momentum.The physical interpretation of Equation ( 9) is straightforward; when a collision occurs with a reservoir the reflected particle has the x component of its momentum changed to a linear combination of the reservoir momentum and its momentum before collision.This results in a strong peak in the reflected momentum near the value of the reservoir momentum which is broadened due to the randomness of the incoming particle momentum, shown in Figure 2.This produces a probability density P (r, v) which is quite different from Gaussian at the local temperature and density.Through interactions with other particles this reflected momentum distribution returns to a more Gaussian-like distribution after four or five particles, regardless of N .
Figure 2. The probability density divided by the local Gaussian probability density for the velocity of particles 1, N/2 and N respectively in an N hard-disk quasi-one-dimensional (QOD) system with reservoirs at the same temperature.Particles 1 and N are nearest the two reservoirs so their probability densities in v x are mirrored, while particle N/2 is in the middle of the QOD system.The horizontal axis is v x and the vertical axis is v y with the colour scale for the probability on the right-hand side.The probability density for particle 1 shows a strong injection peak around the reservoir velocity which is positive.The probability density for particle N/2 is essentially equal one everywhere except for v 2 x + v 2 y > 6 2 where the numerical probability density is zero (thus colored red).

Energy Balance
In a numerical simulation we can calculate the energy flux between system and reservoir by direct calculation or we can calculate it indirectly through the probability density of incident momentum and the collision frequency.Once obtained, the probability density of incident momentum p x will give all of the fluxes to and from the reservoirs per collision.The number of collisions per unit time n L with momentum −p x on the left-hand wall is given by the product of, the length of the collision cylinder |p x |δt for particles which collide within time δt, and the probability density that a particle has momentum −p x is given by P (−p x ).Thus the number of collisions for particles with momentum x )δt and the normalized incoming wall collision probability density is To calculate the energy balance at the left-hand reservoir we use the wall collision distribution and calculate the average energy transfer p 2 x − p 2 x .Using the collision rule from Equation ( 9) and taking the wall momentum as p L = √ 2T L , we look for the value of β for which there is on average no energy transfer in a collision.Therefore we require that The probability density for the incoming momenta gives p x = − π/2β and p 2 x = 2/β, where the incoming collisional temperature T = 1/β.We then obtain a quadratic equation for T /T L , Solving this equation gives Note that as − 2 is always negative when 1 > > 0, the term inside the square root is always positive.
Here we have chosen a temperature to associate with the incoming momentum distribution and found the value which produces an energy balance at the boundary.In this analysis the temperature T only depends on the value of , not on the system size N or the density ρ.As a concrete example, we consider a QOD system of 640 hard disks at density ρ = 0.8 connected to reservoirs with equal temperatures T L = T R = 2, which gives a fitted incoming momentum probability density P (p x ) = αp x exp(− 1 2 βp 2 x ) where α = 0.547 and β = 0.537, implying T = 1.862.This is not a good fit to the data and for the distribution to be properly normalised we require α = β.For = 1/2, there is only one physical solution to Equation ( 12), T = 1.78152, independent of both N and ρ.Whether this temperature should be associated with the particle closest to the boundary or some other bulk measure of the temperature is not immediately clear.In Table 1 we observe that the x and y temperatures of particle 1 converge quickly with increasing N , but to different values, and T x is greater than the boundary temperature of 2.
Table 1.The temperature of the left-hand boundary particle -particle 1 -for QOD systems as a function of system size with fixed reservoir temperatures of T L = T R = 2.The first two columns are for a density ρ = 0.03 while the last two columns are for a density ρ = 0.8.This data is also contained in Figure 3.The best predicted temperature T is very close to T y at ρ = 0.8, but different from the low density result.Numerically, we find that a two parameter form for the probability density gives a much better fit (see Figure 4), and the temperature is more consistent with the temperature at the centre of the system.When γ = 0, α = β and this reduces to the simpler distribution.This distribution can also be used to calculate the energy balance or flux across the boundary by considering the average per collision p 2 x − p 2 x and multiplying this by the collision frequency (for collisions with reservoir I).The parameters for the fits to the incoming momentum probabilities, shown in Table 2, are remarkably insensitive to system size N .
Figure 4.The probability density of incoming momentum (dots) and its functional fit (line) for both the left and right hand reservoirs for a QOD system of 80 hard-disks at a density ρ = 0.8 with = 0.5 and reservoir temperatures T = 2.Note that all incoming momenta for the left hand reservoir must be negative so P L (p x ) = 0 for p x > 0. Similarly, P R (p x ) = 0 for p x < 0. The functional form [Equation (13)] is an excellent fit to the data for both boundary distributions.To produce a nonequilibrium steady state it is sufficient to have reservoirs of different temperatures on each side of the QOD system.The energy entering the system from a boundary with reservoir momentum p I during a collision with a particle of incoming momentum p x is given by The time average of this quantity gives the flux of energy into the system so for a total system energy balance ∆e L must be equal in magnitude but opposite in sign to ∆e R .The energy flux is controlled by both the reservoir momentum p I , the value of , and the probability distribution for the incoming momentum, and goes to zero as → 0.

Temperature Profiles
The generic picture used to construct the canonical ensemble is a system in contact with a heat reservoir; a system in contact with two reservoirs at the same temperature is both canonical and in equilibrium with both reservoirs.In a mechanical dynamical system we expect that an arbitrary initial condition after some transient (a large number of collisions) will come to an equilibrium state, independent of the initial configuration, and remain there indefinitely.After this transient, there is no average energy transfer between system and reservoir.In the equilibrium state the components of the local kinetic temperature, T x and T y , would be expected to become equal, T x = T y .For the QOD system the temperature profiles can be presented as either functions of the average position of the particle, or as a function of particle number (as the order is fixed), but to be able to compare systems of different size we use a simple scaling function X = (i − 0.5)/N where i is the particle number.Figure 5 shows both x and y components of the temperature profile for QOD systems of 3 to 50 disks connected to two reservoirs with the same temperatures T L = T R = 1 as a function of x = (i − 0.5)/N .For these small systems the temperature profiles show very strong number dependence and there is no evidence of equilibration between x and y temperatures before the system size nears 50, but even here T x > T y everywhere in the system.The drop in T y for the two particles in contact with the reservoirs is a consistent feature in all profiles.Further, the temperature in the centre of the system is significantly lower than the temperature of the reservoirs and this varies with system size.The direct numerical calculation of the p x and p y -distributions gives good fits to Gaussians with temperatures equal to local kinetic x and y temperatures, away from the boundaries, thus the second moment calculated from the distribution is consistent with that calculated directly.
A systematic study of the temperature profile for larger systems shows that the same behaviour persists but becomes smaller with increasing N (see Figure 3).There is a strong dependence on system size as the bulk T is lowered with increasing N with no indication that the profiles will converge to some universal N → ∞ limit.The profiles seem almost completely independent of the density, as can be seen by comparing the two graphs in Figure 3.Although the temperature difference between the boundary and particle 1 is not given by the analysis above, we can construct an argument for the systematically lower temperatures in the centre of the QOD systems as N increases.As the difference in temperature between the boundary and T y (1) is independent of N and decreases monotonically with particle number with the same slope regardless of N , a larger system size will lead to a larger decrease and thus a lower centre temperature.These profiles are strongly dependent on system size and there are significant differences between T x and T y for the smaller systems.(b) The temperature profiles for a QOD system of 50 hard-disks showing that despite good equilibration of T x and T y locally, there remains a systematic difference between the two temperatures with T x > T y .(c) The temperature of the middle particle(s) in the QOD systems of different size plotted against 1/N .Whilst initially appearing to converge, for systems of more than 100 disks the change in temperature with system size increases.Notice that the temperature near the centre is more than 20% smaller than the nominal reservoir temperature T = 1.All graphs have the same vertical scale.The temperature differences T x − T y as a function of particle number are shown in Figure 6 for both low density and high density systems and we see in both cases that the effects are mostly visible in the first ten particles regardless of the system size and the density.We will refer to this behaviour as an order one effect as it is independent of system size N and always involves the same number of particles near a boundary.For 80 and 160 particles T x > T y so nowhere in the system is the temperature fully equilibrated.At the larger system sizes of 320 and 640 particles the results for T x − T y look like noise with a mean of zero, so the temperature has equilibrated.Thus the temperature equilibration in the bulk of the system is an order N effect, requiring somewhere between 160 and 320 particles, but near both boundaries there is a region of approximately 10 particles that never equilibrate.

Entropy
We consider a steady system with reservoirs of equal temperature at each end.As can be seen from Figures 3 and 5, even in the absence of an external temperature gradient both the x and y temperatures vary with position and are generally not equal.The Chapman-Enskog expansion solution, to first order in the gradient of the distribution function, is given by f = f loc (1+Φ) where f loc is the local-equilibrium distribution function, β x = 1/T x and β y = 1/T y depend upon x and are allowed to differ locally.Then f loc can be written as The right-hand side depends on x through n(x), β x and β y .This is an extension of the local equilibrium concept to a system which can have differing temperatures in the x and y directions as is commonly found for nonequilibrium steady states [25].The entropy density is obtained at the level of the local-equilibrium approximation as This entropy density obtained is reminiscent of the equilibrium Sackur-Tetrode equation, except here the hydrodynamic fields are local.
In the computer simulation of the QOD system a histogram of the velocity distribution for each particle is collected on a fixed grid size ζ.For a velocity distribution in two variables P (v x , v y ) the distribution at v x = v i and v y = v j is f (i, j) = nP (v x , v y ).The discrete distribution for each histogram cell f i,j is the integral of the probability distribution f (i, j) over the cell so The entropy can then be calculated for each particle as follows The convergence of the entropy density s(x) and the local equilibrium entropy density s loc (x) is quite different.The entropy density requires the convergence of the numerical velocity distribution itself while the local equilibrium entropy density only requires the convergence of the field variables n(x), β x (x) and β y (x).In Figure 7 we see that s loc is an upper bound for s and the two entropies are in good agreement with each other.The number dependence of the entropy is similar to the the number dependence of the temperature as S becomes smaller with increasing N .for ρ = 0.8.These profiles are strongly dependent on system size.They change with system size but the agreement between s loc (x) and s(x) is impressive in all cases.Note that s loc (x) > s(x) for all values of x.The largest difference between s(x) and s loc (x) is at the boundaries where the momentum distributions deviate most greatly from Gaussian.

Entropy Production and Flux Near Reservoirs
It is surprising to find that despite the absence of entropy production and entropy flux in the bulk of the system, it is evident near the reservoirs.There is entropy production in this steady state system and an associated entropy flux directed into each reservoir.In Figure 8, the entropy production calculated from Equation ( 8) is concentrated on the particles closest to each of the reservoirs.The entropy properties come directly from the BGK approximation and the velocity probability density for each particle.
The two-dimensional velocity probability density for each particle P (v) is calculated numerically on a grid of typically 301 × 301 with a resolution of 0.05.The entropy, entropy flux and entropy production are all calculated numerically by integrating this probability density using Equations ( 3), ( 7), ( 8) with an error in the order of the square of the resolution, that is 0.0025.The entropy flux and production are expected to be good representations at low density in the kinetic regime where ρ < 0.1, but would neglect potential contributions at higher density.In the steady state ∂s ∂t = 0 so ∇ • j s = σ ent .As both the entropy current and production are defined at each particle, we can associate these properties with the regions assigned to each particle and approximate the integral in Equation ( 18) as a sum.
Here we can choose the particle number m freely and ∆x m is related to the local density by ρ m = 1 ∆xmLy .We have separated the entropy production into the product of two terms; the first is the local frequency parameter of the BGK approximation ν m , and the other is the integral of velocity probability density.Thus we can write ν m in terms of known quantities calculated, using the BGK approximation only for a m , and ask what value of ν m is required to obtain consistent results.We have In Figure 8 we see that for a system in contact with two the equal reservoirs the entropy production is confined to a small number of particles nearest the wall and that the entropy production is largely independent of the system size (thus order one).The entropy production in a system with equal reservoirs is intensive, and only involves properties of the particles that are in contact with the reservoirs.The numerical convergence of these two quantities is quite different; the entropy production converges with the smoothness of the distribution increasing slowly and monotonically near the reservoirs and decreasing in the central part until a steady value is reached.Statistically, the probability density needs to be sampled sufficiently to obtain a reliable value over the whole two-dimensional space (v , v y ) and the convergence of the local distribution depends on the convergence to the local fields n(x), T x (x) and T y (x).
Using Equation (20) we calculate the BGK frequency that is required to give the simulation results that we observe in Figure 8.These results are plotted in Figure 9 and give a required frequency of about 0.02 at the lower density of 0.03, and a frequency of about 0.6 at the higher density of 0.8.The single particle collision frequency is independent of system size N and is about 2.5 times larger than the BGK frequency at ρ = 0.03.At the higher density ρ = 0.8 the single particle collision frequency is almost 20 times larger than the BGK frequency.This suggests that the BGK frequencies are not connected with the frequencies of collision events in the simulation but rather relate to the longer time scale associated with the relaxation of the distribution.Comparing the two BGK frequencies we see that they are consistent with a linear dependence on density but do not give particular support to that proposition.In Table 3 we give the number of collisions per particle for the different system sizes and densities we have considered.The numerical results obtained here need to be qualified in the sense that the process of heat conduction in this system and the relaxation to a steady state is generally very slow.
Although we base our interpretation on long simulations, we cannot be certain that other slower processes remain undetected.

Conclusions
We have investigated the consequences of a particular microscopic coupling model between a system of hard disks and a reservoir, which has the advantage of being relatively simple and deterministic.The modification to the collision rule at the system-reservoir boundary allows the system to come to a steady state where on average no energy flows between system and reservoir.We construct a qualitative argument for the energy balance that implies a lower temperature than the boundary temperature inside the system which is observed numerically.In this work we have considered reservoirs with equal temperatures to probe the microscopic mechanisms, but it is possible to have reservoirs of different temperatures to produce nonequilibrium steady states.Although the collision rule is deterministic, it does not correspond to deterministic thermostats often used in nonequilibrium steady states [26,27].
In all cases at fixed N the system comes to a steady state but the temperature profiles depend on the system size.For systems of 3 to 50 disks the steady state temperature profiles differ for the x and y directions, but as the size gets larger, somewhere between 160 and 320, local thermal equilibration occurs and T x T y .While these system size effects get smaller with increasing N it is still clear that the temperature in the centre of the system has not converged to a limit at N = 640.Indeed numerically there is no indication that such a limit exists.Evidence that the reservoir affects the temperature profile near the boundary in the same way regardless of system size, leading to the same temperature gradient, suggests that the centre temperature will be smaller for larger system.This is observed numerically, but again the argument is only qualitative.The major unresolved question is whether a large system limit exists for a QOD system with these reservoirs.
As well as systematic changes in the temperature there are also systematic changes in the local entropy that occur with system size.These lead to lower values of entropy in the centre of the system at both low and high densities.Although it can be shown that the entropy derived from the local equilibrium distribution is an upper bound on the local entropy calculated from the numerical distributions, there is remarkable agreement between the two entropies in all circumstances.
The deterministic coupling mechanism produces entropy in particles near the reservoir and that entropy is transported as an entropy flux density into the reservoir.The expected equilibrium configuration of a system between two reservoirs of the same temperature leads to a system that has a steady kinetic temperature which is lower than the temperatures of the reservoirs.By considering the probability distributions of the momentum of particles near the boundary, as well as the details of the reservoir dynamics, we provide a microscopic explanation of this effect.
Two effects that are order one (or independent of N ) are the temperature difference near the boundary T x −T y , and the entropy production and flow near a boundary.The temperature difference is independent of both N and ρ.The entropy production and flux is independent of N but does depend on the density.There is no energy flux or momentum flux associated with this entropy flux.These local effects that occur near the boundaries reinforce the picture of a reservoir with a strong effect on the system locally but sometimes propagates throughout the system, as we see in the profiles of temperature and entropy density.Such reservoir-system interactions are important in a wide variety of situations and here the analysis is made possible because the interaction is specified mechanically. )].

Figure 1 .
Figure 1.Schematic presentation of an N hard-disk quasi-one-dimensional (QOD) system.The height L y is sufficiently small that the disks cannot pass one another.We choose the coordinate origin to be located at the bottom left corner of the system, and system boundaries are denoted by dashed lines.The boundaries at x = 0 and L x are hard walls and those at y = 0, L y periodic (that is, a (H,P) quasi-one-dimensional system).

Figure 3 .
Figure3.The temperature profiles T x (red symbols) and T y (blue symbols) as a function of X = (i − 0.5)/N where i is the particle index for QOD systems of 80, 160, 320 and 640 hard-disks with T L = T R = 2. Plot (a) is the result for ρ = 0.03 and (b) is the result for ρ = 0.8.These profiles are strongly dependent on system size with the bulk temperature smaller with increasing N and there are local differences between T x and T y .The results are remarkably insensitive to density with only very slightly smaller values at the higher density.

Table 2 .
Parameters for the incoming momentum distributions at the left-hand side of QOD systems as a function of system size at density ρ = 0.8 with reservoir temperatures of T L = T R = 2.

Figure 5 .
Figure 5. (a) The temperature profiles T x (top) and T y (bottom) as a function of X for QOD systems of 3, 5, 7, 9 and 15 hard-disks at ρ = 0.8 as a function of X = (i − 0.5)/N .These profiles are strongly dependent on system size and there are significant differences between T x and T y for the smaller systems.(b) The temperature profiles for a QOD system of 50 hard-disks showing that despite good equilibration of T x and T y locally, there remains a systematic difference between the two temperatures with T x > T y .(c) The temperature of the middle particle(s) in the QOD systems of different size plotted against 1/N .Whilst initially appearing to converge, for systems of more than 100 disks the change in temperature with system size increases.Notice that the temperature near the centre is more than 20% smaller than the nominal reservoir temperature T = 1.All graphs have the same vertical scale.

Figure 6 .
Figure 6.The temperature difference T x − T y profiles as a function of the particle index for QOD systems of 80, 160, 320 and 640 hard-disks with T L = T R = 2.The vertical scale on both graphs is the same.(a) Results for ρ = 0.03 and (b) Results for ρ = 0.8.These profiles are remarkably independent of system size N and density ρ.

Figure 7 .
Figure 7.The entropy profiles s(x) (blue symbols) and s loc (x) (red line) as a function of X = (i − 0.5)/N where i is the particle index for QOD systems of 80, 160, 320 and 640 hard-disks with T L = T R = 2. (a) are the results for ρ = 0.03 and (b) are the resultsfor ρ = 0.8.These profiles are strongly dependent on system size.They change with system size but the agreement between s loc (x) and s(x) is impressive in all cases.Note that s loc (x) > s(x) for all values of x.The largest difference between s(x) and s loc (x) is at the boundaries where the momentum distributions deviate most greatly from Gaussian.

Figure 8 .
Figure 8.The entropy flux and entropy production plotted against particle number for an equilibrium QOD system of N hard-disks with T L = T R = 2. System sizes used are N = 80 (filled circles), 160 (open circles), 320 (open diamonds), 640 (crosses) and the results are remarkably insensitive to N .The entropy production is positive near the boundary and the entropy flux is directed into the boundary (reservoirs).Panel (a) is for and ρ = 0.03, (b) is for ρ = 0.8.

Figure 9 .
Figure 9.The recalculated values of the BGK relaxation frequency from the entropy production and entropy flux of QOD systems of 80 (red filled circles) and 160 disks (blue filled diamonds) at two different densities with reservoir temperatures of T L = T R = 2. System sizes used are N = 80, 160 and the results are insensitive to N .Panel (a) is for ρ = 0.03, panel (b) is for ρ = 0.8.

Table 3 .
Simulation lengths given as the number of collisions per particle as a function of system size at the two different densities with reservoir temperatures of T L = T R = 2.