# The Nosé–Hoover thermostated Lorentz gas *footnote **footnote *Note: Most of the figures are in poor quality output. The originals are many MB large. They can be obtained upon request.

## I introduction

In a system of particles under an external force a nonequilibrium steady state can be obtained by applying a thermostat . Deterministic and time reversible bulk thermostating is based on introducing a momentum dependent friction coefficient in the equations of motion. One type of this mechanism is the Nosé–Hoover thermostat . It creates a canonical ensemble in equilibrium and yields a stationary nonequilibrium distribution of velocities in nonequilibrium. Another version, the Gaussian isokinetic thermostat , leads to a microcanonical density for the velocity components in equilibrium and to a constant kinetic energy in nonequilibrium. Though the microscopic dynamics of these thermostated systems is time reversible the macroscopic dynamics is irreversible in nonequilibrium . This is related to a contraction onto a fractal attractor .

Characteristic features of thermostated many particle systems, like a nonequilibrium steady state and a fractal attractor, have been recovered for a specific one particle system, the Gaussian thermostated Lorentz gas . The periodic Lorentz gas consists of a particle that moves through a triangular lattice of hard disks and is elastically reflected at a collision with a diskA model almost identical to the driven periodic Lorentz gas, except for some geometric restrictions, is the Galton board, which has been invented in 1873 to study probability distributions .. It serves as a standard model in the field of chaos and transport, see e. g., . In contrast to many particle systems a one particle system reflects more strongly the properties of a thermostat. For the Gaussian thermostated Lorentz gas a complicated dependence of the attractor on the field strength results .

A second simple one particle system which has been much investigated concerning its chaotic properties is the Nosé–Hoover thermostated harmonic oscillator. In contrast to the Lorentz gas, the dynamics of this system is generically nonergodic . However, one can obtain an ergodic dynamics for this system by the additional control of the square of the kinetic energy .

Recently an alternative thermostating mechanism, called thermostating by deterministic scattering, has been introduced for the periodic Lorentz gas This mechanism has been applied later to a system of hard disk under a temperature gradient and shear.. This deterministic and time reversible mechanism is based on including energy transfer between moving particle and disk scatterer at a collision, instead of using a momentum dependent friction coefficient. It leads to a canonical probability density for the particle in equilibrium, and in nonequilibrium it keeps the energy of the particle on average constant. In Refs.  this model has been compared to the Gaussian thermostated Lorentz gas: In nonequilibrium one finds an attractor for this model which is similar to the fractal attractor of the Gaussian thermostated Lorentz gas, but in contrast to the Gaussian case the attractor is phase space filling even for high field strengths. For both models the conductivity is a nonlinear decreasing function with increasing field strength on a coarse scale.

Based on its construction the method of thermostating by deterministic scattering is in fact closer to the Nosé–Hoover thermostat. This motivates us to apply for the first time the Nosé–Hoover thermostat to the periodic Lorentz gas. In Section II we introduce the Nosé–Hoover thermostat and discuss some variations of it. In Section III we define the periodic Lorentz gas and the thermostats we will study. We investigate these models in equilibrium in section IV and in nonequilibrium in Section V. Conclusions are drawn in Section VI.

## II The Nosé–Hoover thermostat and some variations

In the following sections we consider a one particle system in two dimensions with the position coordinates $`\stackrel{}{q}=(q_x,q_y)`$ and the momentum coordinates $`\stackrel{}{p}=(p_x,p_y)`$. The mass of the particle has been set equal to unity. The equations of motion for the Nosé–Hoover thermostat are then given by 

$`\dot{\stackrel{}{q}}`$ $`=`$ $`\stackrel{}{p}`$ (1)   
$`\dot{\stackrel{}{p}}`$ $`=`$ $`\stackrel{}{\epsilon }\zeta \stackrel{}{p}`$ (2)   
$`\dot{\zeta }`$ $`=`$ $`({\displaystyle \frac{p^2}{2T}}1){\displaystyle \frac{1}{\tau ^2}}.`$ (3)   
The thermostat variable $`\zeta `$ couples the particle dynamics to a reservoir. It controls the kinetic energy of the particle $`p^2/2`$ such that $`<p^2>=2T`$. This holds even in nonequilibrium as induced by an electric field $`\stackrel{}{\epsilon }`$. $`\tau `$ is the response time of the thermostat. Performing the limit $`\tau 0`$ in Eqs. (3) approximates the Gaussian thermostat with $`\zeta =(\stackrel{}{\epsilon }\stackrel{}{p})/p^2`$. In the limit $`\tau \mathrm{}`$ the friction coefficient approaches a constant, $`\zeta =\zeta _c`$, and the equations of motion are not time reversible anymore. The dynamics of this dissipative limit has been investigated in .

A generalization of the Nosé–Hoover thermostat to control higher even moments of $`p`$ has been introduced by Hoover . The moments are fixed according to the momentum relations for the Gaussian distribution. By such a more detailed control of the nonequilibrium steady state statistical dynamical properties, like ergodicity, can be improved, as has been mentioned for the Nosé–Hoover thermostated harmonic oscillator in the introduction. In principle the method of the control of the even moments can straightforwardly be extended to the control of the odd moments. However, such a thermostat would involve an additional parameter to control the respective current of the subsystem. Thus, the corresponding reservoir would be more than a single thermal reservoir, which is physically not desirable, apart from the fact that the respective equations of motion would not be time reversible anymore.

We briefly note that there exist other formal generalizations , or modifications , of the Nosé–Hoover thermostat in the literature. They have been critically reviewed in Refs. .

## III Variations of the Nosé–Hoover thermostat for the periodic Lorentz gas

The basic themostating method we investigate in this paper is the Nosé–Hoover thermostat, Eqs. (3). In the following we introduce three variations of it which are time reversible and result in a dissipative dynamics in nonequilibrium.

The first variation is the Nosé–Hoover thermostat with a field dependent coupling to the reservoir,

$`\dot{\stackrel{}{q}}`$ $`=`$ $`\stackrel{}{p}`$ (4)   
$`\dot{p_x}`$ $`=`$ $`\epsilon _x(1+\epsilon _x)\zeta p_x`$ (5)   
$`\dot{p_y}`$ $`=`$ $`\epsilon _y(1+\epsilon _y)\zeta p_y`$ (6)   
$`\dot{\zeta }`$ $`=`$ $`({\displaystyle \frac{p^2}{2T}}1){\displaystyle \frac{1}{\tau ^2}},`$ (7)   
which is obtained by including the factors $`1+\epsilon _x`$, resp. $`1+\epsilon _y`$, in Eqs. (3). Alternatively, these equations can be written by defining two field dependent friction coefficients, $`\xi _x=(1+\epsilon _x)\zeta `$ and $`\xi _y=(1+\epsilon _y)\zeta `$, which are governed by $`\dot{\xi }_x=(p^2/2T1)(1+\epsilon _x)/\tau ^2`$ and $`\dot{\xi }_y=(p^2/2T1)(1+\epsilon _y)/\tau ^2`$, respectively. It then becomes clear that for each momentum component there exists a separate response time of the reservoir which reads $`\tau /\sqrt{1+\epsilon _x}`$, resp. $`\tau /\sqrt{1+\epsilon _y}`$. The basic advantage of this thermostat is that the response times are now adjusted to the corresponding component of the field strength such that with increasing field strength the response time decreases. The standard Nosé–Hoover thermostat Eqs. (3) is contained as a special case in equilibrium .

The second variation goes back to Hoover . It includes a control of $`<p^4>=8T^2`$,

$`\dot{\stackrel{}{q}}`$ $`=`$ $`\stackrel{}{p}`$ (8)   
$`\dot{\stackrel{}{p}}`$ $`=`$ $`\stackrel{}{\epsilon }\zeta _1\stackrel{}{p}\zeta _2{\displaystyle \frac{p^2}{2T}}\stackrel{}{p}`$ (9)   
$`\dot{\zeta _1}`$ $`=`$ $`({\displaystyle \frac{p^2}{2T}}1){\displaystyle \frac{1}{\tau ^2}}`$ (10)   
$`\dot{\zeta _2}`$ $`=`$ $`{\displaystyle \frac{p^2}{2T}}({\displaystyle \frac{p^2}{2T}}2){\displaystyle \frac{1}{\tau ^2}}.`$ (11)   
As mentioned in the previous sections this variation can improve statistical dynamical properties, like ergodicity.

The third variation controls $`p_x^2`$ and $`p_y^2`$ separately. This is performed by defining two independent reservoirs for the $`x`$– and $`y`$–direction,

$`\dot{\stackrel{}{q}}`$ $`=`$ $`\stackrel{}{p}`$ (12)   
$`\dot{p_x}`$ $`=`$ $`\epsilon _x\zeta _xp_x`$ (13)   
$`\dot{p_y}`$ $`=`$ $`\epsilon _y\zeta _yp_y`$ (14)   
$`\dot{\zeta _x}`$ $`=`$ $`({\displaystyle \frac{p_x^2}{T}}1){\displaystyle \frac{1}{\tau _x^2}}`$ (15)   
$`\dot{\zeta _y}`$ $`=`$ $`({\displaystyle \frac{p_y^2}{T}}1){\displaystyle \frac{1}{\tau _y^2}}.`$ (16)   
This variation more deeply intervenes in the microscopic dynamics by forcing the single components $`p_x^2`$ and $`p_y^2`$ separately towards canonical distributions. However, in contrast to the previous variations a curiosity is hidden in it: At a collision the thermostat $`\zeta _x,\zeta _y`$ is uncorrelated to the thermostated variables $`p_x`$ and $`p_y`$, because $`p_x`$ and $`p_y`$ change at a collision whereas $`\zeta _x,\zeta _y`$ remain the same. Therefore the thermostat does not work efficiently. But we have not found any reflection of this curiosity in the macroscopic behavior.

We study the dynamics of these models in one Lorentz gas cell with periodic boundaries, see Fig. 1(a). As the radius of the disk we take $`r=1`$. For the spacing between two neighboring disks we choose $`w0.2361`$, corresponding to a density equal to $`4/5`$ of the maximum packing density of the scatterers. A collisionless free flight of the particle is avoided for this parameter . The relevant variables of the dynamical system are defined in Fig. 1(b): $`\beta `$ is the angular coordinate of the point at which the particle elastically collides with the disk, and $`\gamma `$ is the angle of incidence at this point.

The equations of motion are integrated by a fourth order Runge Kutta algorithm with a step size of $`dt=0.005`$ between two collisions. The collision of the particle with the disk has been determined with a precision of $`10^7`$. Unless declared otherwise the temperature is chosen to $`T=0.5`$.

## IV Equilibrium

Inserting the initial condition $`(p_0,\zeta _0)=(\sqrt{2T},0)`$ in Eqs. (3) with $`\epsilon =0`$ the velocity of the particle becomes a constant, $`p=\sqrt{2T}`$, that means the Nosé–Hoover thermostat does not act in the Lorentz gas and the dynamics is microcanonical. For other initial conditions one observes in computer simulations that $`p^2`$ and $`\zeta `$ oscillate periodically and that the dynamics of this one particle system is nonergodic.

A stability analysis confirms the numerical results: The Nosé–Hoover thermostated equations of motion Eqs. (3) can be reduced for $`\epsilon =0`$ to

$`\dot{p^2}`$ $`=`$ $`2\zeta p^2`$ (17)   
$`\dot{\zeta }`$ $`=`$ $`({\displaystyle \frac{p^2}{2T}}1){\displaystyle \frac{1}{\tau ^2}}.`$ (18)   
Eqs. (18) are also valid at the moment of a collision, because $`p^2`$ and $`\zeta `$ are not changed by a collision. The fixed point of Eqs. (18) is $`(p^2,\zeta )=(2T,0)`$ with the eigenvalues $`\lambda _{1/2}=\pm \sqrt{4T/\tau ^2}`$ and is thus elliptic.

The additional control of $`<p^4>`$ destroys the microcanonical probability density but it is not sufficient to obtain an “exact” dynamics<sup>§</sup><sup>§</sup>§The dynamics is “exact” if every initial density of nonzero measure converges to the same stationary density  in the periodic Lorentz gas. Different initial conditions still lead to a different shape of the probability density $`\varrho (p_x)`$, as is shown in Fig. 2(a).

In contrast the separate control of $`p_x^2`$ and $`p_y^2`$ leads to an exact dynamics in equilibrium corresponding to the canonical probability density $`\varrho (p_x)`$, as shown in Fig. 2(b).

## V Nonequilibrium

We now apply an external electric field $`\stackrel{}{\epsilon }`$ parallel to the $`x`$–axis. The Nosé–Hoover thermostat and the related models then lead to well defined nonequilibrium steady states with constant average energy of the particle.

### A Probability density $`\varrho (p_x)`$

The probability density $`\varrho (p_x)`$ for the Nosé–Hoover thermostat for $`\epsilon =0.5`$ is presented in Fig. 3(a). For $`\tau ^2=0.01`$ the density shows some remains of the deformed microcanonical density of the Gaussian thermostat, whereas for $`\tau ^2=1`$ and $`\tau ^2=1000`$ the density becomes similar to the density of thermostating by deterministic scattering, which is related to a canonical distribution .

Fig. 3(b) shows $`\varrho (p_x)`$ for the three variations of the Nosé–Hoover thermostat. We have chosen here $`T=0.60029`$ which corresponds to the temperature in the bulk for thermostating by deterministic scattering at a parametric temperature of $`T=0.5`$. The density of the Nosé–Hoover thermostat with field dependent coupling to the reservoir is very close to the density of thermostating by deterministic scattering. The density of the Nosé–Hoover thermostat with separate control of $`p_x^2`$ and $`p_y^2`$ looks like a superposition of the densities of the Nosé–Hoover thermostat for small and for large $`\tau `$.

In all models the mean value of $`\varrho (p_x)`$ is positive, indicating a current parallel to the field direction.

### B Conductivity

The conductivity $`\sigma =<p_x>/\epsilon `$ for the Nosé–Hoover thermostat is shown in Fig. 4. For $`\tau ^2=0.01`$ the curve is very similar to the conductivity of the Gaussian thermostated Lorentz gas . For $`\tau ^2=1000`$ the curve is more stretched along the $`\epsilon `$–axis and globally not decreasing anymore. In contrast the conductivity as obtained from thermostating by deterministic scattering  is a globally decreasing function. According to the Einstein relation, in the limit $`\epsilon 0`$ $`\sigma `$ should approach the equilibrium diffusion coefficient $`D`$ of the periodic Lorentz gas, which for $`w=0.2361`$ has the value $`D0.21`$ . This is hard to see for $`\tau ^2=1000`$, because for $`\epsilon 0`$ the probability density changes drastically from a smooth, canonical like density to a non-smooth density. It is also difficult to see any linear response in computer simulations, as has already been discussed for the Gaussian thermostated Lorentz gas and for thermostating by deterministic scattering in Ref. .

### C Attractor

Fig. 5 shows the Poincaré section of $`(\beta ,\mathrm{sin}(\gamma ))`$Because the symbols for the angles vary in the literature we mention again that the angle $`\beta `$ gives the location of the collision relative to the field direction and $`\gamma `$ is the angle of incidence. at the moment of the collision for the Nosé–Hoover thermostat, for the three variations and for thermostating by deterministic scattering. Again, we have chosen $`T=0.60029`$ to compare the results with thermostating by deterministic scattering. For all models the structure of the attractor is qualitatively the same as the structure of the fractal attractor obtained for the Gaussian thermostated Lorentz gas . However, the fine structure varies with the models and with the response time $`\tau `$. For the Nosé–Hoover thermostat with $`\tau ^2=0.01`$, see Fig 5(a), the structure is most pronounced, whereas for the control of $`p_x^2`$ and $`p_y^2`$ separately, see Fig. 5(e), the structure is least visible.

### D Bifurcation diagram

The angle $`\beta `$ at the moment of the collision is presented as a function of the field strength for the Nosé–Hoover thermostat in Fig. 6. For all three values of $`\tau ^2`$ the attractor is phase space filling for small field strengths $`\epsilon <1.3`$ and contracts onto a periodic orbit with increasing field strength. For $`\tau ^2=0.01`$ the scenario is similar to the one of the Gaussian thermostated Lorentz gas In Fig. 10 of  the angle of flight after a collision is plotted and the field is parallel to the negative x-axis.. For $`\tau ^2=1`$ the scenario looses its richness, but it gets a bit more complicated again for $`\tau ^2=1000`$. In contrast to the Nosé–Hoover thermostat, the attractor of thermostating by deterministic scattering remains phase space filling even for large $`\epsilon `$ .

The bifurcation diagram in the dissipative limit $`\tau \mathrm{}`$ of Eqs. (3) with a constant friction coefficient $`\zeta _c`$ shows an inverse scenario to Fig. 6, as is presented in Fig. 7. For small $`\epsilon `$ the trajectory is a so–called creeping orbit , then it changes to a periodic orbit, and for large $`\epsilon `$ the attractor gets phase space filling. By increasing $`\zeta _c`$ the strength of the dissipation increases with the consequence that the onset of chaotic behavior starts at higher field strengths. With respect to the numerics we remark that related to the different shapes of the attractors by varying $`\zeta _c`$ in the dissipative limit, the duration of the transient behavior of the Nosé–Hoover thermostat grows drastically for $`\tau \mathrm{}`$.

Figs. 8-10 show the bifurcation diagrams for the three variations of the Nosé–Hoover thermostat. In general one observes that even for high field strengths chaotic regions appear, in contrast to the Nosé–Hoover thermostat.

For the variation with the additional control of $`p^4`$ the attractor covers a bounded $`\beta `$–interval for $`\epsilon 3`$. For these field strengths the trajectory is a creeping orbit.

Fig. 10(b) depicts the attractor for the variation with separate control of $`p_x^2`$ and $`p_y^2`$ under different response times for the $`x`$– and $`y`$–direction. Since the field acts in the $`x`$–direction we have chosen a strong coupling, $`\tau _x^2=0.1`$, for the $`x`$–direction and a weak coupling, $`\tau _y^2=1000`$, for the y–direction. Up to a field strength $`\epsilon 6.5`$ no periodic window has been found. The bifurcation diagram for these parameters most strongly deviates from the Nosé–Hoover and Gaussian thermostated Lorentz gas and is closest to the bifurcation diagram of thermostating by deterministic scattering. However, in contrast to thermostating by deterministic scattering the attractor is more concentrated around $`\beta \pi `$.

We detected a numerical problem for the second and the third variation at several values of $`\tau `$. One observes that sporadically after large time intervals there appears a creeping orbit with a very low velocity of the particle, which is difficult to handle numerically. Whether this creeping orbit is the stationary state could probably be clarified by calculating Lyapunov exponents .

### E Thermodynamic entropy production and phase space volume contraction

A characteristic property of the Nosé–Hoover thermostat as well as of the Gaussian isokinetic thermostat is that the thermodynamic entropy production is equal to the phase space volume contraction rate . This equality can as well easily be verified for the variation with the additional control of $`p^4`$ and for the variation with separate control of $`p_x^2`$ and $`p_y^2`$.

On the other hand it does not hold for the Nosé–Hoover thermostat with field dependent coupling to the reservoir. From Eqs. (7) one gets for the phase space contraction rate of this model $`<\mathrm{div}\dot{\mathrm{\Gamma }}>=(2+\epsilon )<\zeta >`$ where $`\dot{\mathrm{\Gamma }}=(\dot{\stackrel{}{q}},\dot{\stackrel{}{p}},\dot{\zeta })`$. The precise relation between thermodynamic entropy production $`\dot{S}_{TD}=\epsilon <p_x>/T`$ and $`<\mathrm{div}\dot{\mathrm{\Gamma }}>`$ for this variation is obtained by calculating the energy balance between subsystem and reservoir:

$$E=\frac{p^2}{2}+T\tau ^2\zeta ^2$$
(19)   
is the total energy which in a nonequilibrium steady state should on average be zero,

$$<\frac{dE}{dt}>=0.$$
(20)   
Inserting Eqs. (7) with $`\epsilon _y=0`$ in Eqs. (20) leads to

$$\frac{\epsilon _x<p_x>}{T}=\frac{\epsilon _x<p_x^2\zeta >}{T}+2<\zeta >.$$
(21)   
Numerical simulations have shown that $`p_x^2`$ and $`\zeta `$ are no independent quantities in nonequilibrium. If $`p_x^2`$ and $`\zeta `$ would be independent and equipartitioning would be fulfilled, i. e., $`<p_x^2>=T`$, which is only the case in equilibrium, then Eq. (21) would lead to an identity between thermodynamic entropy production and phase space volume contraction. The results for $`<\mathrm{div}\dot{\mathrm{\Gamma }}>`$ and for $`\dot{S}_{TD}`$ as obtained from computer simulations for this system are presented in Table I at different $`\tau `$ and $`\epsilon `$. More details of the entropy production in this variation and in a Gaussian thermostat with the same property are discussed in Ref. .

## VI conclusions

We have investigated the Nosé–Hoover thermostat and three variations of it for the periodic Lorentz gas. All models are time reversible and lead to well defined nonequilibrium steady states with a constant average kinetic energy of the moving particle.

As a typical characteristic of deterministic and time reversible thermostating mechanisms it has been confirmed that in nonequilibrium all these systems contract onto attractors similar to the fractal attractor of the Gaussian thermostated Lorentz gas.

In equilibrium only the variation of the Nosé–Hoover thermostat with separate control of $`p_x^2`$ and $`p_y^2`$ leads to an exact dynamics with $`\varrho (p_x)`$ being canonical just like the corresponding density of thermostating by deterministic scatterin g.

In nonequilibrium the attractor of the Nosé–Hoover thermostat contracts onto a periodic orbit for higher field strength, analogous to the Gaussian thermostated Lorentz gas. However, the detailed scenario depends on the value of the response time $`\tau `$. Concerning the probability density in equilibrium, the attractor, and the conductivity in nonequilibrium, the properties of the standard Nosé–Hoover thermostat in the periodic Lorentz gas are closer to the properties of the Gaussian thermostat than to the properties of thermostating by deterministic scattering, although the Nosé–Hoover thermostat and thermostating by deterministic scattering share the property of keeping the energy of the particle on average constant in nonequilibrium, and for both thermostats the probability densities for the velocity components are related to a canonical probability density.

Concerning the bifurcation diagrams we find that the dynamics of the three variations of the Nosé–Hoover thermostat are in general “more chaotic”. Even for higher field strengths there exist pronounced chaotic regions. The separate control of $`p_x^2`$ and $`p_y^2`$ leads to a phase space filling attractor up to high field strengths in nonequilibrium. The equilibrium properties and the bifurcation diagram in nonequilibrium of this model are qualitatively closest to the properties of thermostating by deterministic scattering in comparison to the other versions.

For the Nosé–Hoover thermostat with field dependent coupling to the reservoir the thermodynamic entropy production is generically not equal to the phase space volume contraction in contrast to the other models.

An important question would be to look for common properties of all deterministic and time reversible thermostats. So far, only the existence of a fractal attractor in nonequilibrium appears to be typical. A more detailed investigation of these thermostating mechanisms should in particular involve a quantitative comparison, as, e. g., by means of computing Lyapunov exponents. An answer to this question would be helpful for obtaining a general characterization of nonequilibrium steady states.

Acknowledgments: We dedicate this article to G. Nicolis, a champion of chaoticity, on occasion of his 60th birthday. K. R. and R. K. thank G. Nicolis for his continuous support and for the possibility to collaborate with him on problems of thermostating. This work has been started during the workshop “Nonequilibrium Statistical Mechanics” (Vienna, February 1999). The authors thank the organizers G. Gallavotti, H. Spohn and H. Posch for the invitation to this meeting. K. R. thanks the European Commission for a TMR grant under contract no. ERBFMBICT96-1193. R. K. acknowledges as well financial support from the European Commission.