# Thermostating by Deterministic Scattering: Heat and Shear Flow

## I Introduction

Driving macroscopic systems out of equilibrium requires external forces. Now, the very existence of a nonequilibrium steady state implies that the temperature of the system must remain time independent. One way to prevent the system from heating up indefinitely in nonequilibrium is the introduction of a thermostating algorithm . Starting from molecular dynamics simulations Evans, Hoover, Nosé and others proposed deterministic thermostats to model equilibrium and nonequilibrium fluids . In this formalism the (average) internal energy of the dynamical system is kept constant by subjecting the particles to fictitious frictional forces, thus leading to microcanonical or canonical distributions in phase space . The major feature of this mechanism is its deterministic and time-reversible character, which is in contrast to stochastic thermostats . This allows to elaborate on the connection between microscopic reversibility and macroscopic irreversibility and has led to interesting new links between statistical physics and dynamical systems theory , especially to relations between transport coefficients and Lyapunov exponents , and between entropy production and phase-space contraction . Although used almost exclusively in the context of nonequilibrium systems, the abovementioned thermostating mechanism presents the drawback that the dynamical equations themselves are altered, even in equilibrium. This raises the question of whether some of the results are due to the special nature of this thermostating formalism or are of general validity .

Recently, an alternative mechanism in which thermalization is achieved in a deterministic and time-reversible way has been put forward by Klages et al.  and has been applied to a periodic Lorentz gas under an external field. In the present paper we apply this thermostating method to an interacting many-particle system subjected to nonequilibrium boundary conditions giving rise to thermal conduction and to shear flow. The model is closely related to that of Chernov and Lebowitz , who study a hard disk fluid driven out of equilibrium into a steady state shear flow by applying special scattering rules at the boundaries in which the particle velocity is kept constant. Our model is introduced in Section II where the thermalization mechanism is tested under equilibrium conditions. In Section III we move on to the case of an imposed temperature gradient and a velocity field by adapting the scattering rules, and we compute the respective transport coefficients. Having a deterministic and time-reversible system at hand we proceed in Section IV to investigate the relation between thermodynamic entropy production and phase-space contraction rate in nonequilibrium stationary states. The main conclusions are drawn in Section V.

## II Equilibrium state

Consider a two-dimensional system of hard discs confined in a square box of length $`L`$ with periodic boundary conditions along the x-axis, i.e., the left and right sides at $`x=\pm L/2`$ are identified. At the top and bottom sides of the box, $`y=\pm L/2`$, we introduce rigid walls where the discs are reflected according to certain rules to be defined later. The discs interact among themselves via impulsive hard collisions so that the bulk dynamics is purely conservative. In the following and in all the numerical computations we use reduced units by setting the particle mass $`m`$, the disk diameter $`\sigma `$ and the Boltzmann constant $`k_B`$ equal to one.

Before proceeding to the nonequilibrium case we define the disc-wall collision rules in equilibrium and check whether the system is well-behaved. Now, in equilibrium the bulk distribution is Gaussian with a temperature $`T`$, and the in- and outgoing fluxes at the top and bottom wall have the form (see )

$$\mathrm{\Phi }(v_x,v_y)=(2\pi T^3)^{1/2}|v_y|\mathrm{exp}\left(\frac{v_x^2+v_y^2}{2T}\right),$$
(1)   
with $`v_y<0`$ for the bottom wall and $`v_y>0`$ for the top wall. Imposing stochastic boundary conditions on the systems in this setting would mean that for every incoming particle the outgoing velocities are chosen randomly according to Eq. (1). In practice, this is usually done by drawing numbers from two independent uniformly distributed random generators $`\zeta ,\xi [0,1]`$ and then transforming these numbers with the invertible map $`𝓣\text{ }^1:[0,1]\times [0,1][0,\mathrm{})\times [0,\mathrm{})`$ as

$$(v_x,v_y)=𝓣\text{ }^1(\zeta ,\xi )=\sqrt{2T}(\text{erf}^1(\zeta ),\sqrt{\mathrm{ln}(\xi )}),$$
(2)   
which amounts to transforming the uniform densities $`\rho (\zeta )=1`$ and $`\rho (\xi )=1`$ onto $`\mathrm{\Phi }(v_x,v_y)`$ according to

$$\rho (\zeta )\rho (\xi )\left|\frac{d\zeta d\xi }{dv_xdv_y}\right|=\left|\frac{𝓣\text{ }(v_x,v_y)}{v_xv_y}\right|=(2/\pi T^3)^{1/2}|v_y|\mathrm{exp}\left(\frac{v_x^2+v_y^2}{2T}\right).$$
(3)   
Note that so far we have restricted Eqs. (2),(3) to positive velocities $`v_x,v_y[0,\mathrm{})`$, which implies a normalization factor in Eq. (3) being different to the one of Eq. (1). In analogy to stochastic boundaries we now define the deterministic scattering at the walls as follows. First, take the incoming velocities $`v_x`$, $`v_y`$ and transform them via $`𝓣\text{ }(v_x,v_y)=(\zeta ,\xi )`$ onto the unit square. Second, use a two-dimensional, invertible, phase-space conserving chaotic map $`:[0,1]\times [0,1][0,1]\times [0,1]`$ to obtain $`(\zeta ^{},\xi ^{})=(\zeta ,\xi )`$. Finally, transform back to the outgoing velocities via $`(v_x^{},v_y^{})=𝓣\text{ }^1(\zeta ^{},\xi ^{})`$. In order to render the collision process time-reversible, we also have to distinguish between particles with positive and negative tangential velocities by using $``$ and $`^1`$, respectively. Thus, particles going in with positive (negative) velocities have to go out with positive (negative) velocities and the full collision rules read

$`(v_x^{},v_y^{})`$ $`=`$ $`𝓣\text{ }^1𝓣\text{ }(v_x,v_y),v_x0`$ (4)   
$`(v_x^{},v_y^{})`$ $`=`$ $`𝓣\text{ }^1^1𝓣\text{ }(v_x,v_y),v_x<0,`$ (5)   
where $`𝓣`$  is meant to be applied to the modulus of the velocities . Since both the positive and the negative side of the tangential velocity distribution of Eq. (1) are normalized to $`1/2`$, this normalization factor has to be incorporated in Eq.(3) to render the full desired flux $`\mathrm{\Phi }`$ equivalent to the one of Eq. (1). Rewriting Eq. (3) in polar coordinates yields precisely the transformation used in , in the limiting case where it mimics a reservoir with infinitely many degrees of freedom. It should also be realized that for obtaining the transformation $`𝓣`$  the total number of degrees of freedom of the reservoir has been projected out onto a single velocity variable, which couples the bulk to the reservoir. Eckmann et al.  used a similar idea to go from a Hamiltonian reservoir with infinitely many degrees of freedom to a reduced description when modeling heat transfer via a finite chain of nonlinear oscillators. It remains to assign the form of the chaotic map $``$ and we shall first adopt the choice of a baker map, as in Refs. ,

$$(\zeta ^{},\xi ^{})=(\zeta ,\xi )=\{\begin{array}{cc}(2\zeta ,\xi /2);\hfill & 0\zeta 1/2\hfill \\ (2\zeta 1,(\xi +1)/2);\hfill & 1/2<\zeta 1\hfill \end{array}.$$
(6)   
Later on we will investigate the consequences of choosing other mappings like the standard map (see, e.g., ). Since in equilibrium the in- and outgoing fluxes have the same form as in Eq.(1), the baker map yields a uniform density and our scattering prescription in Eq.(4) can be viewed as a deterministic and time-reversible counterpart of stochastic boundary conditions. $``$ being chaotic, the initial and final momentum and energy of any single particle are certainly different, but both quantities should be conserved on the average. The latter is confirmed by numerical experiments in equilibrium, where, as usual in hard disk simulations, we follow a collision-to-collision approach . Keeping the volume fraction occupied by $`N=100`$ hard disks equal to $`\rho =0.1`$ sets the length of the box equal to $`L=28.0`$. After some transient behavior which depends on the temperature of the initial configuration the bulk distribution is Gaussian with zero mean and mean kinetic energy $`T/2`$ in each directions. The in- and outgoing fluxes at the walls are correctly equipartioned with $`T`$ as well and have the desired form of Eq.(1), so the system reproduces the correct statistic properties. We close this section by a remark on how we measure the temperature of a flux to or from the boundaries. As the temperature of the tangential component we use the variance of the velocity distribution, $`T_x:=(v_xv_x_x)^2_x`$, where $`<>_x`$ denotes an average over the density $`\rho (v_x)`$. On the other hand, since in the normal direction we actually measure a flux, the appropriate prescription to measure the temperature of this component is $`T_y:=\left[v_y\right]_y/\left[v_y^1\right]_y`$, where $`[]_y`$ represents an average over the flux $`\mathrm{\Phi }`$ and the denominator serves as a normalization. The temperatures of the in- and outgoing fluxes at the wall are then defined as $`T_{i/o}:=(T_x+T_y)/2`$, and $`T_w:=(T_i+T_o)/2`$.

## III Nonequilibrium steady state

### A Heat flow

#### 1 The Model

In the following, we explicitly indicate the dependence of the transformation $`𝓣`$  on the parameter $`T`$ by writing $`𝓣\text{ }_T`$. This immediately indicates how we may drive our system to thermal nonequilibrium: We just have to use different values of this parameter for the upper ($`T^u`$) and the lower wall ($`T^d`$). We deliberately avoid to use the word ’temperature’ for this parameter, since in contrast to stochastic boundary conditions we have generally no idea how a different $`T`$ affects the actual temperature of the wall in the sense of the definition given above. In a nonequilibrium situation the temperature of the ingoing flux $`\mathrm{\Phi }_i`$ generally does not match exactly the parameter $`T`$. Therefore, we do not transform onto the uniform invariant density of the baker map anymore, and consequently, the outgoing flux might have all kinds of shapes or temperatures. Nevertheless, the hope is that the mapping $``$ will be chaotic enough to smooth out most of the differences and to produce a reasonable outgoing flux $`\mathrm{\Phi }_o`$ such that the system is correctly thermostated. And this is indeed what we find in the numerical experiments.

#### 2 Numerical Results

We set $`T^u=2`$, $`T^d=1`$ and $`\rho =0.1`$ and average over about 40000 particle-particle collisions per particle and about 6000 particle-wall collisions per particle. We divide the available vertical height $`L1`$ into 20 equally spaced horizontal layers and calculate the time averages of the number density $`n(y)`$ of the particles, the mean velocities $`u_x(y)=v_x`$, $`u_y(y)=v_y`$, and the variances $`(v_xu_x)^2`$, $`(v_yu_y)^2`$. Furthermore, we record the time average of the kinetic energy transfer and measure the temperatures of the in- and outgoing fluxes of both walls as described in the preceding section. The temperatures at the walls are then defined as the mean value of the in- and outgoing temperatures, $`T_w^{u/d}:=(T_i^{u/d}+T_o^{u/d})/2`$. Time series plots of these quantities confirm the existence of a nonequilibrium stationary state (NSS) induced by the temperature gradient.

Fig. 1 shows the temperature profile between the upper and the lower wall. Apart from boundary effects it is approximately linear, and the respective kinetic energy is equipartitioned between the two degrees of freedom. The parameters $`T^u`$ and $`T^d`$ are represented as (\*), and we find a ’temperature’ jump, whereas the measured temperatures $`T_w^{u/d}`$ (+) at the walls seem to continue the bulk profile reasonably well. The profile of the number density $`n=4/\pi \rho `$ is depicted in Fig. 2. Note again the boundary effects. The densities of the in-coming particles at the upper wall are Gaussian shaped (Figs. 3a,c), whereas the outgoing densities (Figs. 3b,d) show cusps due to the folding property of the baker map. Nevertheless, the baker map produces a reasonable outgoing flux which generates a NSS.

In order to examine the bulk behavior we now compute the thermal conductivity in our computer experiment and compare it to the theoretical value. For this purpose, we measure the heat flux $`Q`$ across the boundaries and estimate the temperature gradient $`dT(y)/dy`$ by a linear least square fit to the experimental profile. To discard boundary effects we use only data in the bulk of the system, namely from layer 3 to layer 18, i.e., excluding the top two and the bottom two layers. The experimental heat conductivity is then defined as

$$\lambda _{exp}=Q\left(\frac{dy}{dT}\right),$$
(7)   
whereas the theoretical expression for the conductivity of a gas of hard disks with unit mass and unit diameter as predicted by Enskog’s theory reads 

$$\lambda _l=1.0292\sqrt{\frac{T}{\pi }}\left[\frac{1}{\chi }+\frac{3}{2}bn+0.8718(bn)^2\chi \right].$$
(8)   
Here, $`b`$ is the second virial coefficient, $`b=\pi /2`$, and $`\chi `$ is the Enskog scaling factor, which is just the pair correlation function in contact,

$$\chi =\frac{1\frac{7}{16}\frac{\pi }{4}n}{(1\frac{\pi }{4}n)^2}.$$
(9)   
Since (8) depends on local values of $`T`$ and $`n`$ we define the theoretical effective conductivity $`\lambda _{th}`$ as the harmonic mean over the layers ,

$$\lambda _{th}=\left(1/N_{layers}\underset{l=1}{\overset{N_{layers}}{}}1/\lambda _l\right)^1.$$
(10)   
Table I compares $`\lambda _{exp}`$ to $`\lambda _{th}`$ by showing the ratio of the experimental to the theoretical conductivity for different particle numbers and temperature differences. The agreement is quite good, so our thermostating mechanism produces a NSS which is in agreement with hydrodynamics.

Furthermore, going into the hydrodynamic limit by increasing the number of particles we observe that the discontinuity in the outgoing flux of Figs. 3(b,d) diminishes, as expected, since both the in- and outgoing flux come closer to local equilibrium.

### B Shear Flow

Inspired by the recently proposed model of Chernov and Lebowitz  for a boundary driven planar Couette-flow in a nonequilibrium steady state, we now proceed to check whether it is possible to combine our thermostating mechanism with a positive (negative) drift imposed onto the upper (lower) wall, respectively. Chernov and Lebowitz chose a purely Hamiltonian bulk and simulated the drift at the boundaries by rotating the angle of the particle velocity at the moment of the scattering event with the wall while keeping the absolute value of the velocity constant. This setting could be formulated in a time-reversible way and keeps the total energy of the system generically constant. Here we separate the thermostating mechanism and the drift of the walls by introducing the map

$$𝒮_d(v_x,v_y)=(v_x+d,v_y),$$
(11)   
and by applying this shift to the ’thermostated’ velocities. Time-reversibility forces us to do the same before thermostating. Thus, the full particle-wall interaction reads

(Model I) (12)   
$`(v_x^{},v_y^{})`$ $`=`$ $`𝒮_d𝓣\text{ }_T^1𝓣\text{ }_T𝒮_d(v_x,v_y),v_xd`$ (13)   
$`(v_x^{},v_y^{})`$ $`=`$ $`𝒮_d𝓣\text{ }_T^1^1𝓣\text{ }_T𝒮_d(v_x,v_y),v_x<d,`$ (14)   
where shifts of different sign are used for the upper (lower) wall to let the walls move into opposite directions. Other prescriptions to impose a shear will be investigated in more detail in the following section.

In the simulations we set $`d=\pm 0.05`$, $`T=T^u=T^d=1.0`$, $`N=100`$ and $`\rho =0.1`$. As we expected, we find a NSS with a linear shear profile along the x-direction (Fig. 4), where the drift velocity of the wall $`u_w`$ (\*) is defined as the average between the in- and outgoing tangential velocities. The temperature profile is shown in Fig. 5, with the wall temperatures $`T_w`$ (+) defined as above. As can be seen in the plots, none of these values correspond to the parameters $`T`$ (\*) or $`d`$. Nevertheless, we obtain a linear shear profile $`u_x(y)`$ and an almost quadratic temperature profile $`T(y)`$, as predicted by hydrodynamics .

For a comparison of experimental and theoretical viscosity we follow the same procedure as above, i.e., we estimate the experimental shear rate $`du_x(y)/dy`$ by a linear least square fit $`u_x(y)=\gamma y`$, again discarding the outermost layers. Through the measured momentum transfer from wall to wall $`\mathrm{\Pi }`$ the experimental viscosity is given as

$$\eta _{exp}=\mathrm{\Pi }/\gamma ,$$
(15)   
and the theoretical value as calculated by the Enskog theory has the form 

$$\eta _l=1.022\frac{1}{2}\sqrt{\frac{T}{\pi }}\left[\frac{1}{\chi }+bn+0.8729(bn)^2\chi \right].$$
(16)   
Again $`b=\pi /2`$ denotes the second virial coefficient, $`\chi `$ is given by Eq.(9) and we use the arithmetic mean of the viscosity over layers 3 - 18 to compute the theoretical viscosity, i.e., $`\eta _{th}=1/N_{layers}\eta _l`$. Table II shows good agreement between the experimental and the theoretical values for different numbers of particles, so again the thermostating mechanism leads to the correct macroscopic behavior. However, the discontinuities in the $`v_x`$-velocity distribution of the scattered particles should be noticed (see Fig.6). This time these discontinuities do not diminish or disappear in the hydrodynamic limit.

## IV Entropy production and Phase-Space Contraction

Having a deterministic and reversible dynamics at hand we can now turn to properties beyond the usual hydrodynamic ones and, in particular, investigate the conjectured identity between phase space contraction rate and thermodynamic entropy production in the light of our formalism . In an isolated macroscopic system the entropy is a thermodynamic potential and therefore plays the central role in determining the time evolution and the final equilibrium state. Yet, its microscopic interpretation out of equilibrium remains controversial (see, e.g., ) and the situation is even much less clear in NSS. For the class of models where thermostating is ensured by friction coefficients an exact equality between entropy production and phase space contraction rate in NSS has been inferred on the basis of a global balance between the system and the reservoir . For the Chernov-Lebowitz model an approximate equality has also been found . Still, it is not clear under which circumstances this relation holds in general .

### A Equilibrium State

We begin with the simplest case of equilibrium described in section II, where the thermodynamic entropy production $`\overline{R}_{eq}`$ vanishes. The bulk dynamics being Hamiltonian, phase-space contraction can only occur during collisions with a wall. Since these collisions take place ’ instantaneously’, we ignore the bulk particles and restrict ourselves to the compression related to a single collision during the time interval $`dt`$. The phase space contraction is then given by the ratio of the one-particle phase-space volume after the collision $`(dx^{}dy^{}dv_x^{}dv_y^{})`$ to the one before the collision, $`(dxdydv_xdv_y)`$, and can thus be obtained from the Jacobi determinant of the scattering process. One easily sees that $`dx^{}=dx`$ and $`|dy^{}/dy|=|v_y^{}dt/v_ydt|`$ . Furthermore,

$`\left|{\displaystyle \frac{dv_x^{}dv_y^{}}{dv_xdv_y}}\right|`$ $`=`$ $`\left|{\displaystyle \frac{𝓣\text{ }}{v_xv_y}}{\displaystyle \frac{^{(1)}}{xy}}{\displaystyle \frac{𝓣\text{ }^1}{v_x^{}v_y^{}}}\right|=\left|{\displaystyle \frac{𝓣\text{ }}{v_xv_y}}\left[{\displaystyle \frac{𝓣\text{ }}{v_x^{}v_y^{}}}\right]^1\right|`$ (17)   
$`=`$ $`\left|{\displaystyle \frac{v_y}{v_y^{}}}\right|\mathrm{exp}\left({\displaystyle \frac{v_x^2+v_y^2v_x^2v_y^2}{2T}}\right),`$ (18)   
where step two follows from the phase-space conservation of $`/^1`$, and the last line is obtained from Eqs.(2)(3). Hence, in a particle-wall collision the phase-space volume is changed by a factor of

$$\left|\frac{dv_x^{}dv_y^{}dx^{}dy^{}}{dv_xdv_ydxdy}\right|=\mathrm{exp}\left(\frac{v_x^2+v_y^2v_x^2v_y^2}{2T}\right).$$
(19)   
The mean exponential rate of compression of the phase space volume per unit time is thus given by

$$\overline{P}=<\mathrm{ln}\left|\frac{dv_x^{}dv_y^{}dx^{}dy^{}}{dv_xdv_ydxdy}\right|>=(v_x^2+v_y^2v_x^2v_y^2)/2T,$$
(20)   
where the brackets $`<>`$ denote a time average over all collisions at the top and the bottom walls. In equilibrium the in- and outgoing fluxes associated to these collisions have the same statistical properties, so $`\overline{P}_{eq}`$ sums up to zero and

$$\overline{P}_{eq}=\overline{R}_{eq},$$
(21)   
This is fully confirmed by the simulations.

### B NSS

The thermodynamic entropy production $`\sigma `$ per unit volume of our system in NSS is given by the Onsager form 

$$\sigma (y)=\frac{\mathrm{\Pi }}{T}\frac{du_x}{dy}+J(y)\frac{d}{dy}\left(\frac{1}{T}\right),$$
(22)   
where $`\mathrm{\Pi }`$ is the x-momentum flux in the negative y-direction, and $`J(y)`$ is the heat flux in the positive y-direction.

#### 1 Heat Flow

Imposing only a temperature gradient on our system like in section III A, the first term in Eq. (22) is identical to zero and the total entropy production $`\overline{R}`$ in the steady state is then

$$\overline{R}=_{Volume}\sigma 𝑑𝐫=_{Surface}J/T𝑑s=J_w^u/T_w^u+J_w^d/T_w^d$$
(23)   
The right hand side of Eq.(23) is the outward entropy flux $`J_w/T_w`$ across the walls of the container. Note that there is no temperature slip at the walls with respect to the correctly defined temperature values, as indicated by the simulation results in Figs.1 and 5.

On the other hand, the exponential phase-space contraction rate now reads

$$\overline{P}=(v_x^2+v_y^2v_x^2v_y^2)/2T_u_u+(v_x^2+v_y^2v_x^2v_y^2)/2T_d_d,$$
(24)   
where we averaged over the upper and the lower wall separately. Since in NSS $`J_w^u=(v_x^2+v_y^2v_x^2v_y^2)/2_u=J_w^d=(v_x^2+v_y^2v_x^2v_y^2)/2_d`$, the ratios of entropy production to exponential phase space contraction rate reduce to

$$\frac{\overline{R}^{u/d}}{\overline{P}^{u/d}}=\frac{T^{u/d}}{T_w^{u/d}}.$$
(25)   
In the hydrodynamic limit the in- and the outgoing fluxes approach local equilibrium, implying $`T_i^{u/d}T_o^{u/d}T_w^{u/d}T^{u/d}`$ for both walls. Therefore, the ratios in Eq.(25) should go to unity. The numerical results in Table III confirm this expectation, leading to a good agreement between entropy production and exponential phase-space contraction rate.

#### 2 Shear Flow

We follow the same procedure as in the preceding section. For a stationary shear flow the hydrodynamic entropy production $`\sigma `$ per unit volume in Eq. (22) can be written as 

$$\sigma (y)=\frac{\mathrm{\Pi }}{T}\frac{du_x}{dy}+J(y)\frac{d}{dy}\left(\frac{1}{T}\right)=\mathrm{\Pi }\frac{d}{dy}\left(\frac{u_x}{T}\right).$$
(26)   
The second step in Eq.(26) follows from the fact that in NSS $`\lambda dT/dy=J(y)=\mathrm{\Pi }u_x(y)`$. The total entropy production $`\overline{R}`$ in the steady shear flow state is then 

$`\overline{R}={\displaystyle _{Volume}}\sigma 𝑑𝐫={\displaystyle _{Surface}}\mathrm{\Pi }u/T𝑑s`$ (27)   
$`=J_w/T_w=2L^2\mathrm{\Pi }\left(u_w/L\right)/T_w=L^2\mathrm{\Pi }\gamma /T_w.`$ (28)   
In the macroscopic formulation of irreversible thermodynamics Eq.(27) is interpreted as an equality, in the stationary state, between the entropy produced in the interior and the entropy flow carried across the walls. Our shift map $`𝒮_d`$ mimics moving walls with drift velocities $`\pm u_w`$. The work performed at these walls is converted by the viscous bulk into heat and then again absorbed by the walls which now act as infinite thermal reservoirs. By imagining that the walls act as an ’equilibrium’ thermal bath at temperature $`T_w`$, $`\overline{R}`$ can be interpreted as their entropy increase rate.

For Model I (Eqs.(13)(14)) the mean exponential phase-space contraction rate takes the form ($`𝒮_d/v_xv_y1`$)

$$\overline{P}=[v_x^2+v_y^2v_x^2v_y^22d(v_x^{}+v_x)]/2T,$$
(29)   
whereas the entropy production is given by

$$J_w/T_w=[v_x^2+v_y^2v_x^2v_y^2<v_x^{}>^2+<v_x>^2]/2T_w.$$
(30)   
In Table IV the ratios of $`L^2\mathrm{\Pi }\gamma /J_w`$ as obtained from the simulations are reported and the relation of phase-space contraction rate to entropy production is subsequently checked. Whereas the equality between entropy production and entropy flow via heat transfer is confirmed, we observe a significant difference between entropy production and phase-space contraction which subsists in the hydrodynamic limit. This mismatch is perhaps not so unexpected in view of the distorted outgoing fluxes (see Fig.6). Nevertheless, one could argue that the fine structure of these distributions may depend on the specific characteristics of the baker map chosen to model the collision process. We therefore also considered the standard map

$$\stackrel{~}{}:\{\begin{array}{c}\xi ^{}=\xi \frac{k}{2\pi }\mathrm{sin}(2\pi \zeta ),\hfill \\ \zeta ^{}=\zeta +\xi ^{},\hfill \end{array}$$
(31)   
with the parameter $`k=100`$ to ensure that we are in the hyperbolic regime . We found that the discrepancy between entropy production and phase-space contraction rate remains: although the deterministic map now seems to be chaotic enough to smooth out the fine structure of the outgoing densities, the discontinuity at $`d`$ survives. Actually, as long as Model I is adopted it becomes clear that in NSS there will always be more out- than ingoing particles with $`v_xd`$ at the upper wall (and with $`v_xd`$ at the lower wall). Thus, the Gaussian halves in Fig.7 (b) will never match to a full Gaussian even in the hydrodynamic limit, and $`\mathrm{\Phi }_i`$ and $`\mathrm{\Phi }_o`$ will never come close to a local equilibrium. To circumvent this problem we modify the map $`𝓣`$  and investigate the following model:

(Model II) (32)   
$`𝓣\text{ }_\pm (v_x,v_y)`$ $`=`$ $`({\displaystyle \frac{\text{erf}\left[(|v_x|d)/\sqrt{2T}\right]\pm \text{erf}(d/\sqrt{2T})}{1\pm \text{erf}(d/\sqrt{2T})}},\mathrm{exp}(v_y^2/2T))`$ (33)   
$`(v_x^{},v_y^{})`$ $`=`$ $`𝓣\text{ }_+^1\stackrel{~}{}𝓣\text{ }_{}(v_x,v_y),v_x0`$ (35)   
$`(v_x^{},v_y^{})`$ $`=`$ $`𝓣\text{ }_{}^1\stackrel{~}{}^1𝓣\text{ }_+(v_x,v_y),v_x<0.`$ (36)   
This model is also time-reversible, but in contrast to the former one no particle changes its tangential direction during the scattering. There is still a gap in the outgoing distribution of Fig.(8) (b), however, simulations show that this gap disappears in the hydrodynamic limit thus bringing the in- and the outgoing distributions close to local equilibrium. Furthermore, we note that whereas we were not able to give a relation between the parameter $`d`$ and the actual wall velocity $`u_w`$ for Model I, in case of Model II $`u_w`$ converges to $`d`$ in the hydrodynamic limit. For this reason we chose $`d=0.5`$ in the following, since this value yields the same order of the wall velocity as $`d=0.1`$ for Model I. Proceeding now to the phase-space contraction rate we find that it takes the form

$$\overline{P}^{u/d}=(n_+n_{})\mathrm{ln}\frac{1\pm \text{erf}(d/\sqrt{2T})}{1\text{erf}(d/\sqrt{2T})}[v_x^2+v_y^2v_x^2v_y^22d(v_x^{}+v_x)]/2T_{u/d},$$
(37)   
where $`n_\pm `$ are the collision rates for positive and negative tangential velocities, and the additional term (cf. Eq.(29)) results from the different denominators in Eq.(33). Note that one has to average over the upper and the lower wall separately. Again we compare $`\overline{R}`$ and $`\overline{P}`$ (Table V), but although the outgoing flux approaches now a Gaussian in the hydrodynamic limit the two quantities still do not match. This result can be understood in more detail by rearranging the terms in Eq.(37) as

$`\overline{P}^{u/d}`$ $`=`$ $`[v_x^2+v_y^2v_x^2v_y^2<v_x^{}>^2+<v_x>^2]/2T_{u/d}`$ (39)   
$`[<v_x^{}>^2<v_x>^22d(v_x^{}+v_x)]/2T_{u/d}(n_+n_{})\mathrm{ln}{\displaystyle \frac{1\pm \text{erf}(d/\sqrt{2T})}{1\text{erf}(d/\sqrt{2T})}}.`$   
Since $`TT_w`$ in the hydrodynamic limit, the first term clearly corresponds to the entropy production Eq.(30). However, the second and the third terms provide additional contributions. For $`u_wd`$ and $`d0`$ they are both of order $`d^2`$ and can be interpreted as a phase space contraction due to a friction parallel to the walls . These two terms apparently depend on the specific modeling of the collision process at the wall. They may physically be interpreted as representing certain properties of a wall, like a roughness, or an anisotropy. Actually, the second term already appeared in Model I, see Eq.(29). The price we had to pay in Model II for the fluxes getting close to local equilibrium is the additional third term in Eq.(39), which does not compensate the second one.

The foregoing analysis shows clearly what to do to get rid of the additional term in Eq.(37): we have to use the same forward and backward transformations $`𝓣\text{ }_\pm `$ in Eqs.(33,35). If one still wants to transform onto a full Gaussian in the hydrodynamic limit time-reversibility has to be given up. This leads us to propose the model

(Model III) (40)   
$`𝓣\text{ }_{}(v_x,v_y)`$ $`=`$ $`({\displaystyle \frac{\text{erf}\left[(v_xd)/\sqrt{2T}\right]+1}{2}},\mathrm{exp}(v_y^2/2T))`$ (41)   
$`(v_x^{},v_y^{})`$ $`=`$ $`𝓣\text{ }_{}^1\stackrel{~}{}𝓣\text{ }_{}(v_x,v_y).`$ (43)   
which is still deterministic, but no longer time-reversible. The phase-space contraction is now given as

$$\overline{P}=[v_x^2+v_y^2v_x^2v_y^22d(v_x^{}v_x)]/2T.$$
(44)   
Fig. 9 shows that the in- and the outgoing fluxes are getting close to local equilibrium, implying that the velocity of the wall $`u_w`$ goes to $`d`$ and the wall temperature $`T_w`$ goes to $`T`$. Consequently, Eq.(44) should converge to the correct thermodynamic entropy production of Eq.(30) in the hydrodynamic limit, and this is indeed what we observe in Table V. This implies that time-reversibility does not appear to be an essential ingredient for having a relation between phase space contraction and entropy production, as was already stated in Refs. . We remark that we consider the lack of time-reversibility in Model III rather as a technical difficulty of how we define our scattering rules than a fundamental property of this model.

## V Conclusion

We have applied a novel thermostating mechanism to an interacting many-particle system. Under this formalism the system is thermalized through scattering at the boundaries while the bulk is left Hamiltonian. We have shown how this deterministic and time-reversible thermostating mechanism is related to conventional stochastic boundary conditions. For a two-dimensional system of hard disks, this thermostat yields a stationary nonequilibrium heat or shear flow state. Transport coefficients obtained from computer simulations, such as thermal conductivity and viscosity, agree with the values obtained from Enskog’s theory.

Having a time-reversible and deterministic system we also examined the relation between microscopic reversibility and macroscopic irreversibility in terms of entropy production. We find that entropy production and exponential phase-space contraction rate do in general not agree. When the NSS is created by a temperature gradient both quantities converge in the hydrodynamic limit. By subjecting the system to a shear we examined three different versions of scattering rules, of which one (Model III) produced an agreement.

Our results indicate that neither time-reversibility nor the existence of a local thermodynamic equilibrium at the walls are sufficient conditions for obtaining an identity between phase space contraction and entropy production. A class of systems where such an identity is guaranteed by default are the ones thermostated by velocity-dependent friction coefficients . We suggest that in general, that is, by using other ways of deterministic and time-reversible thermostating, such an identity may not necessarily exist. We would expect the same to hold for any system where the interaction between bulk and reservoir depends on the details of the microscopic scattering rules.

As a next step it would be important to compute the spectrum of Lyapunov exponents for the models presented in this paper. This would enable to check, for example, the validity of formulas which express transport coefficients in terms of sums of Lyapunov exponents , and the existence of a so-called conjugate pairing rule of Lyapunov exponents . Moreover, it would be interesting to verify the fluctuation theorem, as it has been done recently for the Chernov-Lebowitz model .

Acknowledgments   
Helpful discussions with P.Gaspard, M.Mareschal and K. Rateitschak are gratefully acknowledged. R.K. wants to thank the Deutsche Forschungsgemeinschaft (DFG) for financial support. This work is supported, in part, by the Interuniversity Attraction Pole program of the Belgian Federal Office of Scientific, Technical and Cultural Affairs and by the Training and Mobility Program of the European Commission.