Entropy Production in Random Billiards

TL;DR

This paper introduces random billiards as a model for nonequilibrium systems, deriving a general entropy production formula that links mathematical and thermodynamic entropy, with explicit examples involving boundary thermostats.

Contribution

It develops a general entropy production rate formula for random billiards and connects it to thermodynamic entropy, including explicit models with boundary thermostats.

Findings

Derived a universal entropy production rate formula for random billiards.

Established a relation between entropy production and thermodynamic entropy.

Proved ergodicity and computed stationary distributions for specific models.

Abstract

We introduce a class of random mechanical systems called random billiards to study the problem of quantifying the irreversibility of nonequilibrium macroscopic systems. In a random billiard model, a point particle evolves by free motion through the interior of a spatial domain, and reflects according to a reflection operator, specified in the model by a Markov transition kernel, upon collision with the boundary of the domain. We derive a formula for entropy production rate that applies to a general class of random billiard systems. This formula establishes a relation between the purely mathematical concept of entropy production rate and textbook thermodynamic entropy, recovering in particular Clausius' formulation of the second law of thermodynamics. We also study an explicit class of examples whose reflection operator, referred to as the Maxwell-Smolukowski thermostat, models systems with boundary thermostats kept at possibly different temperatures. We prove that, under certain mild regularity conditions, the class of models are uniformly ergodic Markov chains and derive formulas for the stationary distribution and entropy production rate in terms of geometric and thermodynamic parameters.

Equations

$$\pi :N=\{(q,v)\in TM:q\in \partial M\}\to M$$

$$N_q^{\pm }:=\{v\in N_q:\pm ⟨v,\text{mathbbm}{n}_q{⟩}_q\ge 0\}.$$

$$\mathcal{T}(x):={\phi }_{t(x)}(x).$$

$$(\nu P_q)(f)={\int }_{N_q^{-}}P(f)(x)d\nu (x).$$

$${\rho }_q(v)=2\pi {\left(\frac{\beta (q)m}{2\pi }\right)}^{\frac{n+1}{2}}|⟨v,\text{mathbbm}{n}_q⟩|\exp \left\{-\beta (q)\frac{m|v{|}_q^2}{2}\right\}$$

$$J:N\to N,J(q,v)=(q,J_{q}v)=(q,-v)$$

$${\mathcal{R}}_q(u,v)=(J_{q}v,J_{q}u)$$

$$d{\zeta }_q(u,v)=d{\mu }_q^{-}(u)d{P}_{(q,u)}(v),$$

$$\mathcal{D}:=\{(x,y)\in N^{+}\times N^{+}:\pi (y)=\pi (\mathcal{T}(x))\}.$$

$$(\mathcal{B}f)(x)={\mathcal{B}}_x(f)=({\delta }_x\mathcal{B})(f)=P_{\mathcal{T}(x)}(f)$$

$$(\mathcal{B}f)(x)=\mathbb{E}\left[f(X_{i+1})|X_i=x\right]$$

$${\mathcal{D}}_{[0,n]}=\{(x_0,\dots ,x_n):(x_i,x_{i+1})\in \mathcal{D},i=0,\dots ,n-1\}.$$

$$\mathcal{R}:(x_0,\dots ,x_n)\mapsto (\mathcal{J}{x}_n,\dots ,\mathcal{J}{x}_0).$$

$${\int }_{N^{+}}f(x)d\nu (x)={\int }_{N^{+}}f(y)d{P}_{\mathcal{T}(x)}(y)d\nu (x).$$

$$d{\mathbb{P}}_{[0,n]}(x_0,\dots ,x_n)=d\nu (x_0)d{\mathcal{B}}_{x_0}(x_1)\cdots d{\mathcal{B}}_{x_{n-1}}(x_n).$$

$$({\mathcal{B}}^{-}f)(x):=\mathbb{E}\left[f(Y_{i+1})|Y_i=x\right].$$

$$d{\mathbb{P}}_{[0,n]}^{-}(x_0,\dots ,x_n)=d\nu (x_0)d{\mathcal{B}}_{x_0}^{-}(x_1)\cdots d{\mathcal{B}}_{x_{n-1}}^{-}(x_n).$$

$$H({\mathbb{P}}_1,{\mathbb{P}}_2):=\{\begin{array}{ll}{\int }_{\mathcal{D}}\log \left(\frac{d{\mathbb{P}}_1}{d{\mathbb{P}}_2}\right)d{\mathbb{P}}_1& \text{ if }{\mathbb{P}}_1\ll {\mathbb{P}}_2\text{ and }\log \left(\frac{d{\mathbb{P}}_1}{d{\mathbb{P}}_2}\right)\in L^1(\mathcal{D},{\mathbb{P}}_1)\\ +\infty & \text{ otherwise.}\end{array}$$

$$e_p:=\underset{n\to \infty }{\lim }\frac{1}{n}H\left({\mathbb{P}}_{[0,n]},{\mathbb{P}}_{[0,n]}^{-}\right)$$

$$d\nu (x)=d{\nu }_q(u)d({\pi }_{\ast }\nu )(q).$$

$$d{\mu }^{\pm }(x)=d{\mu }_q^{\pm }(u)d\bar{A}(q)$$

$$d\sigma (q,u)=C⟨u,\text{mathbbm}{n}_q⟩d{V}_q^{\text{S}}(u)dA(q),$$

$$e_p=-{\text{bigintssss}}_{\partial M}\frac{{\nu }_q^{+}(E_0)-{\nu }_q^{-}(E_0)}{\kappa T(q)}d\bar{A}(q)-{\text{bigintssss}}_{N^{+}}\frac{\Phi (\pi (\mathcal{T}(x)))-\Phi (\pi (x))}{\kappa T(\pi (x))}d\nu (x)\ge 0.$$

$$P_x(f)=\alpha (q){\mu }_q^{+}(f)+(1-\alpha (q))f({\text{emRef}}_{q}x).$$

$$d\nu (q,v)={\rho }_q(v)d{V}_q(v)d\bar{A}(q)$$

$${\nu }^{\mathrm{s}}=(I-Q{)}^{-1}\pi ,$$

$$Q_{ij}=(1-{\alpha }_i)p_{ij}{A}_i/A_j,{\pi }_i={\bar{A}}_i{\alpha }_i{\mu }_i^{\mathrm{s}}$$

$$e_p=-\sum _{j=1}^N\frac{{\nu }_j^{+}(E_0)-{\nu }_j^{-}(E_0)}{\kappa T_j}=-\sum _{i,j=1}^N\frac{{\nu }_j^{+}(E_0)-{\nu }_i^{+}(E_0)}{\kappa T_j}{p}_{ij}\frac{A_i}{A_j}$$

$$\Vert \mathcal{B}f{\Vert }_2^2={\int }_{N^{+}}|P_{\mathcal{T}(x)}(f){|}^{2}d\nu (x)\le {\int }_{N^{+}}{\int }_{N_{\pi (\mathcal{T}(x))}^{+}}|f(y){|}^{2}d{P}_{\mathcal{T}(x)}(y)d\nu (x)=(\nu \mathcal{B})(|f{|}^2).$$

$$⟨{\mathcal{B}}^{-}f,g⟩={\int }_{N^{+}}({\mathcal{B}}^{-}f)(x)\overline{g(x)}d\nu (x)=\mathbb{E}\left[({\mathcal{B}}^{-}f)(X)\overline{g(X)}\right]=\mathbb{E}\left[\mathbb{E}\left[f(Y_{i+1})|Y_i=X\right]\overline{g(X)}\right].$$

Full text

We consider a class of random mechanical systems called random billiards to study the problem of quantifying the irreversibility of nonequilibrium macroscopic systems. In a random billiard model, a point particle evolves by free motion through the interior of a spatial domain, and reflects according to a reflection operator, specified in the model by a Markov transition kernel, upon collision with the boundary of the domain. We derive a formula for entropy production rate that applies to a general class of random billiard systems. This formula establishes a relation between the purely mathematical concept of entropy production rate and textbook thermodynamic entropy, recovering in particular Clausius’ formulation of the second law of thermodynamics. We also study an explicit class of examples whose reflection operator, referred to as the Maxwell-Smoluchowski thermostat, models systems with boundary thermostats kept at possibly different temperatures. We prove that, under certain mild regularity conditions, the class of models are uniformly ergodic Markov chains and derive formulas for the stationary distribution and entropy production rate in terms of geometric and thermodynamic parameters.

1 Introduction

Model and results. The overarching focus of this paper is to study the dynamics and thermodynamic properties of a class of random mechanical systems, referred to as random billiards, that serve as a concrete model for the rigorous, analytic study of nonequilibrium phenomena. Random billiards can be seen as a Markov chain model variation on mathematical billiards with a random reflection in place of the classical law of specular reflection. The random reflections are specified by a choice of Markov transition operator, which we call the random reflection operator, defined on the space of post-collision velocities.

The main results of the paper are summarized as follows, with more detailed and explicit statements to come in the next section. We use the notion entropy production rate, defined in terms of the Kullback-Leibler divergence of a general stochastic process, to give a characterization of the irreversibility of random billiard systems in terms of the thermodynamic and geometric parameters that characterize the Markov transition operator. In particular, we give an explicit formula for the entropy production rate in terms of the stationary distribution of our Markov chain model and show that the entropy production rate is always a non-negative quantity. The result, which can be seen as the second law of thermodynamics for random billiard models, lays the groundwork for a more detailed study of how certain aspects of the models, such as the geometry of the billiard table or the choice of random reflection operator, affect the irreversibility of a system. In the second half of the paper, we study random billiard systems with a random reflection operator called the Maxwell-Smoluchowski reflection model which models a boundary thermostat, or external heat bath, that keeps the system out of equilibrium. We show that random billiard systems with the Maxwell-Smoluchowski reflection model are uniformly ergodic Markov chains, give an explicit description of the stationary distribution, and compute the entropy production rate analytically and numerically for various examples.

Main technical issues. We aim to formalize problems in nonequilibrium statistical physics which are related to the second law of thermodynamics. In particular, a primary problem is to quantify the notion of irreversibility of macroscopic systems, which are defined by reversible microscopic behavior, using entropy production rate. This has been done for a large class of mesoscopic models, namely stochastic differential equations (SDE), and a large class of microscopic models, namely countable state Markov chains; see the book [17] for a comprehensive overview. However, as far as we are aware, a rigorous study of entropy production rate for continuous state Markov chains, particularly those with non-compact state space, has not been done. The present work addresses this gap in the literature in the specific context of random billiard Markov chains, but the techniques and statements of results should hold for more general non-compact state space Markov chain models. The main work in establishing the entropy production rate formula and subsequent second law of thermodynamics is in showing that the probability measure for the Markov chain in forward time is mutually absolutely continuous with respect to the time reversed measure, from which it follows through careful calculations that the entropy production rate is well defined and satisfies a formula that makes explicit its dependence on the steady state stationary measure of the system. The work here requires a bit more detailed analysis than the SDE case because we cannot use the Cameron-Martin-Girsanov formula that is available for diffusion processes.

In the second part of the paper, we use coupling in order to prove the uniform ergodicity of random billiard models with the Maxwell-Smoluchowski reflection model. Such techniques are now standard, but there are always context specific details that need to be addressed, particularly in the case of non-compact state spaces. It should also be emphasized here that our model is particularly amenable to exact analytic study. Beyond proving uniform ergodicity, we are able give an explicit formula for the stationary measure, even in the case when the boundary thermostats are kept at different temperatures. The formula is given both in terms of temperatures and the geometry of the billiard table of the system. We believe that the techniques here should also serve as a basis for the further study of random billiards with boundary thermostats modeled by other reflection laws.

Related work. The work in this paper fits into several established, overlapping lines of research—random billiard models, models of heat conduction, and the mathematical analysis of nonequilibrium phenomena through entropy production rate. These lines of research are united by the theme of studying models that are simple enough to be amenable to rigorous analysis but rich enough to demonstrate transport phenomena. While this philosophy is taken by much of the related work to be discussed below, the present work is distinguished in the way that explicit formulas linking entropy production rate, temperature, and microscopic parameters defining the system are attained.

In random billiard models, much work has been done in studying the ergodicity of models with independent, identically distributed (i.i.d.) reflections. In [10], uniform ergodicity is proved for random billiards with the Knudsen reflection model that fixes the speed of the particle on tables with $C^{1}$ boundary. This work was extended in [4] to more general, but fixed speed reflection models, along with more relaxed restrictions on the boundary of the billiard domain. The issue of ergodicity for non-i.i.d. models, derived from Markov operators that model more complex boundary microstructures is arguably more subtle. Along these lines, ergodicity and rate of approach to equilibrium have been studied for certain classes of boundary microstructure in [11, 12, 13]. For reflection models that model boundary thermostats less is known. In [5], it is shown that the Maxwell-Boltzmann distribution is the equilibrium distribution for a large class of examples with a single heat bath. Ergodicity was proved for reflection models with multiple temperatures in [18], but this is for a special class of billiard tables, using the geometry of dispersing Sinai billiards.

While there has been a resurgence of interest in the recent past for mechanical models of heat conduction and the study of nonequilibrium behavior in systems with thermostatted boundaries, there are few results that give explicit descriptions of the steady state distribution of the system and few results that give analytic, functional relationships between thermodynamic quantities like the entropy production rate to parameters of the system beyond thermostat temperatures. We should note that in the case of anharmonic oscillators driven by stochastic heat baths, there is in fact a fairly complete set of results. Under certain general assumptions, they can be modeled by stochastic differential equations. There are results on the existence and uniqueness of stationary distributions [8] and rate of convergence, as well as other statistical properties, are well understood [26, 25]. Moreover, techniques as in [17] can be used to characterize and study entropy production rate [27]. On the other hand, these are results at the mesoscopic scale, but things are much less complete for models defined at the microscopic scale, which include random billiards, from which mesoscopic models should arise through universal limiting laws. In addition to random billiard models, other microscopic models include those based on purely Hamiltonian mechanical systems coupled to stochastic heat baths [22, 3, 19, 22, 9, 29] and hybrid models consisting of mechanical systems with a stochastic component meant to make the systems more tractable [21]. While much of the work on these models has been in studying the stationary states and rate of convergence to stationarity, the study of entropy production is a bit more limited. In [6], fluctuations of entropy production rate are studied for perturbed Lorentz gas. Surveys such as [16, 7, 15, 24, 28] study entropy production for general classes of dynamical systems and Markov processes, but this approach is a bit different than the one we take, where we are concerned with how explicit, microscopic details of the model affect entropy production rate.

Organization of the paper. Section 2 introduces definitions and notation, and concludes with more precise statements of main results. Section 3 establishes basic facts about random billiard Markov chains and random reflection operators. Section 4 establishes the preliminary details for the Second Law of Thermodynamics for random billiards and gives a proof of this result. Section 5 introduces random billiards with the Maxwell-Smoluchowski thermostat model. There, we prove the uniform ergodicity of such Markov chains, giving an explicit formula for the stationary distribution, present an explicit formula for entropy production rate, and present some analytical and numerical results for some representative examples. We conclude with Section 6 where we give a final example which shows that random billiard systems can be used to produce work against an external force, creating what we call random billiard heat engines. We present a short numerical study of efficiency and work production for random billiard heat engines.

2 Definitions and Main Results

Let $M$ be a smooth manifold of dimension $n$ with boundary $\partial M$. The boundary may contain points of non-differentiability, where the tangent space to $\partial M$ is not defined. Points where the tangent space is defined will be called regular. The notion of a manifold with corners as defined in [20] is general enough to include all the examples of interest to us. As the issue of regularity is not central to the results of this paper, we simply assume throughout that all points under consideration are regular.

Let $\pi:TM\rightarrow M$ indicate the tangent bundle to $M$. The notations $x=(q,v)\in TM$ and $v\in T_{q}M$ will be used. We are mostly concerned with tangent vectors at boundary points. Thus it makes sense to introduce the bundle

[TABLE]

corresponding to the restriction of $\pi$ to the boundary. Let $N_{q}:=N\cap T_{q}M$. We equip $M$ with a Riemannian metric $\langle\cdot,\cdot\rangle$. The inward pointing unit normal vector to the boundary at a regular point $q$ will be written $\mathbbm{n}_{q}$. Vectors $v\in N_{q}$ pointing to the interior of $M$ constitute the set $N_{q}^{+}$ of post-collision velocities; the negative of vectors in $N_{q}^{+}$ constitute the set $N_{q}^{-}$ of pre-collision velocities. Thus

[TABLE]

The disjoint union of the $N_{q}^{\pm}$ over the regular boundary points defines $N^{\pm}$.

We suppose that $\partial M$ has finite $(n-1)$-dimensional volume and denote by $dA(q)$ the Riemannian volume element of $\partial M$ at $q$. The notation $\bar{A}=A/A(\partial M)$ will indicate the normalized volume when $\partial M$ has finite volume. Let $t\mapsto\varphi_{t}(x)$ denote the geodesic flow in $TM$, which is only defined for values of $t$ up to the moment when geodesics reach the boundary. Recall that the geodesic flow is the Hamiltonian flow (on the tangent bundle) for the free motion of a particle of mass $m$ with kinetic energy $E(q,v)=\frac{1}{2}m|v|_{q}^{2}$, where $|v|_{q}^{2}=\langle v,v\rangle_{q}$. For $x=(q,v)\in N^{+}$, let $t(x):=\inf\{t>0:\pi(\varphi_{t}(x))\in\partial M\}$. The return map (to the boundary) $\mathcal{T}:N^{+}\rightarrow N^{-}$ is defined as

[TABLE]

Upon reaching the boundary, the billiard trajectory (i.e., an orbit of the geodesic flow) is reflected back into the manifold by a choice of map from $N^{-}$ to $N^{+}$; in deterministic billiards, the standard choice is the specular reflection $v\mapsto v-2\langle v,\mathbbm{n}_{q}\rangle_{q}\mathbbm{n}_{q}$. (The theory of standard billiard dynamical systems largely deals with planar systems, as in [2], but see also [14].) For random billiards, this is done by means of a choice of reflection operator $P_{q}$ at each boundary point $q$, to be defined shortly.

It is possible to make this set-up more general by adding conservative forces, in which case the geodesic flow should be replaced with the Hamiltonian flow corresponding to a choice of potential function $\Phi:M\rightarrow\mathbb{R}$. Most of the paper will be concerned with free motion ($\Phi=0$) between collisions, except in Theorem 8 (the Second Law of Thermodynamics), where we assume the more general setting. We then indicate by $E_{0}$ the above kinetic energy function and use $E$ to denote the total energy $E=E_{0}+\Phi\circ\pi$. The introduction of potential forces does not affect the form of the reflection law, which is intended to model impulsive forces; that is, very strong forces acting on very short time intervals.

The definition of reflection operators given next is motivated by natural boundary conditions for the Boltzmann equation involving gas-surface interaction. See for example Chapter 1 in [1].

The space of Borel probability measures on a topological space $X$ will be indicated by $\mathcal{P}(X)$. Two probability measures will be called equivalent if they are mutually absolutely continuous. We often denote by $\mu(f):=\int_{X}f(x)\,d\mu(x)$ the integral of a function $f$ on $X$ with respect to a measure $\mu$. A measure $\mu_{\lambda}$ will be said to depend measurably on elements $\lambda$ of a measurable space if for any bounded continuous function $f$ on $X$ the map $\lambda\mapsto\mu_{\lambda}(f)$ is measurable.

Definition 1** (General reflection operator).**

A general reflection operator at a regular boundary point $q$ is a map $v\in N_{q}^{-}\mapsto P_{(q,v)}\in\mathcal{P}(N_{q}^{+})$. A general reflection operator $P$ on $M$ is the assignment of such an operator for each regular boundary point $q$. We suppose that $P$ depends measurably on $x=(q,v)$ in the sense that for any given bounded continuous function $f$ on $N^{+}$, the map $x\mapsto P(f)(x):=P_{x}(f)$ is measurable.

The operator $P$ on $M$ gives rise at each $q$ to a map $P_{q}$ from $\mathcal{P}(N_{q}^{-})$ to $\mathcal{P}(N_{q}^{+})$ defined by $\nu\mapsto\nu P_{q}$, where the integral of a test function $f$ on $N_{q}^{+}$ with respect to the latter is

[TABLE]

A reflection operator will be defined as a general reflection operator satisfying the condition of reciprocity, whose definition depends on the notion of a Maxwellian probability measure defined below. It is through the notion of reciprocity that the property of $\partial M$ having a given, fixed, temperature at a boundary point will be defined. Let $dV_{q}(v)$ denote the Riemannian volume measure on $T_{q}M$. A measure $\mu$ on $N_{q}$ is said to have density $\rho(v)$ if $d\mu(v)=\rho(v)\,dV_{q}(v)$.

Definition 2** (Maxwellian at temperature $T$).**

The Maxwellian (or Maxwell-Boltzmann probability distribution) at boundary point $q\in M$ and temperature $T(q)$ is the probability measure $\mu^{\pm}_{q}\in\mathcal{P}(N_{q}^{\pm})$ having density

[TABLE]

where $m$ is the mass of the billiard particle, $\beta(q)=1/\kappa T(q)$, and $\kappa$ is known as the Boltzmann constant.

At each regular boundary point $q$ consider the following linear involutions: the flip map

[TABLE]

and the time reversal map

[TABLE]

on $N_{q}^{-}\times N_{q}^{+}$. Given a general reflection operator $P$, define $\zeta_{q}\in\mathcal{P}(N^{-}_{q}\times N^{+}_{q})$ as

[TABLE]

where $\mu_{q}^{-}$ is a Maxwellian at $q$. If $F:X\rightarrow Y$ is a measurable map between measure spaces and $\mu$ is a probability measure on $X$, the push-forward of $\mu$ under $F$ is the probability measure $F_{*}\mu$ on $Y$ defined by $(F_{*}\mu)(E):=\mu(F^{-1}(E))$, where $E$ is a measurable subset of $Y$. If $F$ is a self-map of $X$, then $\mu$ is said to be invariant under $F$ if $F_{*}\mu=\mu$.

Definition 3** (Reciprocity).**

The general reflection operator $P$ has the property of reciprocity if at each regular boundary point $q$ the probability measure $\zeta_{q}$ (just defined above) is invariant under the time-reversal map $\mathcal{R}_{q}$. A general reflection operator satisfying reciprocity will be called simply a reflection operator.

The notion of reciprocity may be interpreted as a local detailed thermal equilibrium of the boundary at $q$. It says that if a particle of mass $m$ hits the boundary at $q$ with a random pre-collision velocity distributed according to the Maxwellian at temperature $T(q)$, then it is reflected with the same distribution (at the same temperature), and this random scattering process is time-reversible in the stochastic sense. A more general, and somewhat more technical, definition of reciprocity than that of Definition 3 will be formulated at the beginning of Section 4. Theorem 8, the second law of thermodynamics, will be proved for random billiards satisfying this more general notion.

The billiard map of a standard deterministic billiard system is the composition of the geodesic translation $\mathcal{T}:x\in N^{+}\mapsto\mathcal{T}(x)\in N^{-}$ defined above and specular reflection at $\pi(\mathcal{T}(x))$. In a random billiard system the second map is replaced with the random scattering post-collision velocity distributed according to $P_{\mathcal{T}(x)}$. See Figure 1. Define

[TABLE]

Then an orbit of the random billiard system is a sequence $\dots,x_{-1},x_{0},x_{1},\dots$ for which every pair $(x_{i},x_{i+1})$ belongs to $\mathcal{D}$. If the random billiard map sends $x$ to $y$, the probability distribution of $y$ is given by $\mathcal{B}_{x}:=P_{\mathcal{T}(x)}$. In what follows, it will be convenient to refer to the map $x\mapsto\mathcal{B}_{x}$ itself as the billiard map.

Definition 4** (Random billiard map).**

The map $\mathcal{B}:N^{+}\rightarrow\mathcal{P}(N^{+})$ defined by $\mathcal{B}_{x}=P_{\mathcal{T}(x)}$, where $P$ is a reflection operator, will be called the random billiard map associated to $P$.

We may use at different places the various notations

[TABLE]

where $\delta_{x}$ is the point mass supported at $x$ and $f$ is a test function on $N^{+}$.

Given a choice of initial probability distribution, we obtain from $\mathcal{B}$ a Markov chain $X_{0},X_{1},\dots$ on the state space $N^{+}$. Note that

[TABLE]

where $\mathbb{E}$ indicates (conditional) expectation. We define the space of finite chain segments of length $n+1$ by

[TABLE]

Note that $\mathcal{D}=\mathcal{D}_{[0,1]}$. The notion of entropy production rate, to be considered shortly, involves the time reversal of the process. Clearly, simple reversal of the Markov chain, in which a chain segment $(x_{0},\dots,x_{n})$ is mapped to $(x_{n},\cdots,x_{0})$, cannot correspond to the physical idea of time reversal—the direction of velocities must also be reversed. In order to define proper time reversal we introduce the map $\mathcal{J}=J\circ\mathcal{T}:N^{+}\rightarrow N^{+}$.

Definition 5** (Proper time reversal map).**

The proper time reversal map, or simply the reversal map, on chain segments is the map $\mathcal{R}:\mathcal{D}_{[0,n]}\rightarrow\mathcal{D}_{[0,n]}$ defined by

[TABLE]

A probability measure $\nu$ on $N^{+}$ is stationary for the random billiard process if $\nu\mathcal{B}=\nu$. Applied to a test function $f$ on $N^{+}$, this condition means that

[TABLE]

It is useful to also define $\nu^{-}:=\mathcal{T}_{*}\nu\in\mathcal{P}(N^{-})$. For emphasis we may on occasion write $\nu^{+}:=\nu$. The issue of existence, uniqueness, regularity, and ergodicity of stationary measures will be addressed in Section 5 for a general class of examples.

Let $X_{0},X_{1},\dots$ be a stationary Markov chain with state space $N^{+}$, transition operator $\mathcal{B}$, and stationary probability $\nu$. For a test function $f$ on $N^{+}$, $\mathbb{E}\left[f(X_{i})\right]=\nu(f).$ Finite chain segments in $\mathcal{D}_{[0,n]}$ are distributed according to the probability measure $\mathbb{P}_{[0,n]}$ defined by

[TABLE]

Given a stationary chain segment $X_{0},\dots,X_{n}$, let $(Y_{0},\dots,Y_{n})=\mathcal{R}(X_{0},\dots,X_{n})$ be its proper time reversal. We introduce the operator $\mathcal{B}^{-}$ defined, for now, on bounded continuous functions:

[TABLE]

We previously defined the probability measure $\mathbb{P}_{[0,n]}$ on the space of finite chain segments $\mathcal{D}_{[0,n]}$. Correspondingly, we define the probability measure $\mathbb{P}^{-}_{[0,n]}$ on the same space by

[TABLE]

Following [17] we make the definitions given next.

Definition 6** (Relative entropy).**

Suppose that $\mathbb{P}_{1}$ and $\mathbb{P}_{2}$ are two probability measures on a measurable space $(\mathcal{D},\mathcal{F})$. The relative entropy of $\mathbb{P}_{1}$ with respect to $\mathbb{P}_{2}$ is defined as

[TABLE]

Definition 7** (Entropy production rate).**

The entropy production rate of the stationary Markov chain $X_{0},X_{1},\dots$ defined by $\nu$ and $\mathcal{B}$ is defined by

[TABLE]

where $H\left(\mathbb{P}_{[0,n]},\mathbb{P}^{-}_{[0,n]}\right)$ is the relative entropy of $\mathbb{P}_{[0,n]}$ with respect to $\mathbb{P}^{-}_{[0,n]}$ restricted to the $\sigma$-algebra generated by $X_{0},\dots,X_{n}$.

Let $\nu_{q}$ be the probability measure on $N^{+}_{q}$ obtained by disintegrating $\nu$ with respect to $\pi:N^{+}\rightarrow\partial M$. So, if $x=(q,u)\in N^{+}$,

[TABLE]

Recall that $\nu^{-}=\mathcal{T}_{*}\nu$, and $\mu_{q}^{+}$ is the Maxwellian at boundary point $q$ for temperature $T(q)$. The following assumptions will be made concerning the stationary measure $\nu$:

$\pi_{*}\nu\ll A$ where $A$ is the $(n-1)$-dimensional Riemannian volume on $\partial M$; 2. 2.

The measures $\nu_{q}$, $\mu^{+}_{q}$, $J_{*}\nu_{q}^{-}$ are mutually equivalent.

We define the measures $\mu^{\pm}\in\mathcal{P}(N^{\pm})$ by

[TABLE]

where $\bar{A}=A/A(\partial M)$ is the normalized Riemannian volume measure on $\partial M$. Moreover, let $S\subset N$ denote the bundle of unit vectors in $N$, $S^{+}=S\cap N^{+}$, and $S^{+}_{q}=S^{+}\cap N_{q}$. Let $\sigma\in\mathcal{P}(S^{+})$ be given by

[TABLE]

where $C$ is a normalizing constant and $V_{q}^{\text{\tiny S}}$ is the Riemannian volume measure on $S_{q}^{+}$.

Theorem 8** (Second law of thermodynamics).**

Let $E_{0}(q,u)=\frac{1}{2}m|u|^{2}_{q}$ denote the kinetic energy function and let $\Phi(q)$ be a bounded measurable function defined on the boundary of $M$ which denotes the potential function of the system. Suppose that a stationary probability measure $\nu$ for the random billiard map $\mathcal{B}$ exists and satisfies assumptions 1 and 2 above. Then

[TABLE]

In words, the entropy production rate (per collision with the boundary) is the average of the energy gained at each iteration of the random billiard system divided by the temperature at the point of collision.

That the displayed quantity in the theorem is non-negative means that, in the absence of a potential energy function, the forward direction in time for the Markov chain is distinguished from the time-reversed process in the following sense: on average, energy is extracted from the regions of higher temperature of $\partial M$ and released into the regions of lower temperature.

With this theorem in hand, the importance of systematically understanding the stationary measure $\nu$ for random billiard systems is readily seen. The rest of our main results are centered around the study of the stationary measure for a class of examples whose reflection operator is of the Maxwell-Smoluchowski kind.

Definition 9** (Maxwell-Smoluchowski model).**

Let $\text{\em Ref}_{q}$ denote the specular reflection at the regular boundary point $q$, and fix $\alpha(q)\in(0,1]$. Let $\mu_{q}^{+}$ be the Mawellian at $q$. Define $P_{q}$ by its evaluation on a test function $f$ as follows:

[TABLE]

Thus the surface scattering process defined by this general reflection operator, known as the Maxwell-Smoluchowski model, has the effect of mapping a pre-collision velocity $v$ of an incoming particle at $q$ to the random post-collision velocity $V$ whose probability distribution is $\mu_{q}^{+}$ with probability $\alpha(q)$ and the specular reflection of $v$ with probability $1-\alpha(q)$.

The next theorem, informally stated below and stated precisely in Theorem 21 and Proposition 22 of Section 5, considers a random billiard system with Maxwell-Smoluchowski reflection operator whose boundary $\partial M$ is partitioned into $N$ components $\Gamma_{i}$ with temperatures $T(q)\equiv T_{i}$ and constants $\alpha(q)\equiv\alpha_{i}$, for $i=1,\ldots,N$.

Theorem 10**.**

Under mild regularity conditions on $\partial M$ (see Assumption 19), the random billiard Markov chain $(X_{n})_{n\geq 0}$ is uniformly ergodic. Moreover, when the boundary temperatures $T_{1},\ldots,T_{N}$ are equal, say to a constant $T_{0}>0$, the stationary distribution $\nu$ is given by

[TABLE]

where $\rho_{q}(v)$ is the Maxwellian density given in (1) with constant temperature $T(q)\equiv T_{0}$. When the boundary temperatures are not equal, the stationary distribution, expressed in polar coordinates as a measure on $S^{+}\times(0,\infty)$ is given by $d\nu(q,u,r)=d\nu^{\mathrm{s}}_{q}(r)d\sigma(q,u)$, where $\nu^{\mathrm{s}}_{q}$ is the stationary distribution of the speed after collision with boundary point $q$ and $\sigma$ is defined in Equation (2). The measure $\nu_{q}^{\mathrm{s}}$ is constant on components $\Gamma_{i}\subset\partial M$, so we let $\nu_{i}^{\mathrm{s}}:=\nu_{q}^{\mathrm{s}}$ for $q\in\Gamma_{i}$. Letting $\nu^{\mathrm{s}}$ be the $N$-dimensional vector with components $\nu_{i}^{\mathrm{s}}$,

[TABLE]

where $Q$ is an $N\times N$ matrix and $\pi$ an $N$-dimensional vector given by

[TABLE]

with $A_{i}=A(\Gamma_{i})$ and $\bar{A}_{i}=\bar{A}(\Gamma_{i})$ the volume measure and normalized volume measure of boundary component $\Gamma_{i}$ respectively. Finally, the entropy production rate is given by

[TABLE]

where $\nu^{+}=\nu$, $\nu^{-}=\mathcal{T}_{*}\nu$, $\nu_{j}^{\pm}$ is the restriction to $N_{j}^{\pm}=\{(q,u)\in N^{\pm}:q\in\Gamma_{i}\}$, and $p_{ij}$ is the conditional probability of the billiard particle colliding next with boundary component $\Gamma_{j}$ given that its current position is on boundary component $\Gamma_{i}$.

In Subsection 5.3, these results are used to compute the entropy production rate explicitly for a series of examples, emphasizing the influence of system parameters.

3 Basic facts about the Markov chain

Consider the Hilbert space $\mathcal{H}:=L^{2}(N^{+},\nu)$. No confusion should arise if we use the same notation $\langle\cdot,\cdot\rangle$ for the inner product on $\mathcal{H}$ as for the Riemannian metric on $M$. We now consider $\mathcal{B}$ as an operator on $\mathcal{H}$.

Proposition 11**.**

The billiard map $\mathcal{B}$, regarded as an operator on $\mathcal{H}$, has norm $\|\mathcal{B}\|_{2}=1$.

Proof.

This is an immediate consequence of the Cauchy-Schwarz inequality and the fact that $\nu$ and $P_{x}$ are probability measures:

[TABLE]

But stationarity implies $(\nu\mathcal{B})(|f|^{2})=\nu(|f|^{2})=\|f\|_{2}^{2}$, so that $\|\mathcal{B}\|_{2}\leq 1$. The observation $\mathcal{B}1=1$ concludes the proof. ∎

Proposition 12** (Transition operator for the time-reversed chain).**

Suppose that $\nu$ and $\mathcal{J}_{*}\nu$ are equivalent measures, so that the Radon-Nikodym derivative $\rho^{\nu}:=\frac{d(\mathcal{J}_{*}\nu)}{d\nu}$ is defined. Then $\mathcal{B}^{-}=\rho^{\nu}\mathcal{J}\mathcal{B}^{*}\mathcal{J}.$ Here $\mathcal{B}^{*}$ is the adjoint operator to $\mathcal{B}$, $\mathcal{J}$ is the composition operator $\mathcal{J}f:=f\circ\mathcal{J}$, and $\rho^{\nu}$ is identified with its multiplication operator.

Proof.

Let $X$ be a random variable with probability distribution $\nu$. Let the function $f$ on $N^{+}$ be in the domain of $\mathcal{B}^{-}$ and consider, for $g\in\mathcal{H}$, the inner product

[TABLE]

Now $Y_{i}=\mathcal{J}X_{n-i}$ and $Y_{i+1}=\mathcal{J}X_{n-i-1}$. Letting $l=n-i-1$, the rightmost term above becomes

[TABLE]

We then have

[TABLE]

Therefore $\mathcal{B}^{-}=\mathcal{J}^{*}\mathcal{B}^{*}\mathcal{J}$. Since $\left\langle\mathcal{J}f,g\right\rangle=\int_{N^{+}}(f\circ\mathcal{J})\overline{g}\,d\nu=\int_{N^{+}}f\overline{(g\circ\mathcal{J})}\rho^{\nu}\,d\nu=\left\langle f,\rho^{\nu}\mathcal{J}g\right\rangle$ we obtain $\mathcal{J}^{*}=\rho^{\nu}\mathcal{J}$, concluding the proof. ∎

A similar argument shows that the process corresponding to the simple reversal

[TABLE]

has transition operator $\mathcal{B}^{*}$. Moreover $\nu\mathcal{B}^{*}=\nu$ since $\nu\mathcal{B}^{*}f=\langle\mathcal{B}^{*}f,1\rangle=\langle f,\mathcal{B}1\rangle=\langle f,1\rangle=\nu(f).$ On the other hand, $\nu\mathcal{B}^{-}=\mathcal{J}_{*}\nu$ and it is not the case in general that $\nu=\mathcal{J}_{*}\nu$, so $\nu$ may not be stationary with respect to $\mathcal{B}^{-}$.

From the billiard map $\mathcal{B}$ and a stationary probability measure $\nu$ we define the probability measure $\eta\in\mathcal{P}(\mathcal{D})$ by $d\eta(x,y)=d\nu(x)\,d\mathcal{B}_{x}(y)$ and call $\eta$ the probability measure on forward pairs. The probability measure on backward pairs $\eta^{-}\in\mathcal{P}(\mathcal{D})$ is defined by $\eta^{-}=\mathcal{R}_{*}\eta$, where we recall that $\mathcal{R}(x,y)=(\mathcal{J}y,\mathcal{J}x)$. We assume that $\eta$ and $\eta^{-}$ are in the same measure class.

Proposition 13**.**

The measure $\eta^{-}$ satisfies $d\eta^{-}(x,y)=d\nu(x)\,d\mathcal{B}^{-}_{x}(y)$ for $(x,y)\in\mathcal{D}$.

Proof.

It suffices to show that the two measures, $\eta^{-}$ and the measure defined by the right-hand side of the equation, give the same integral on functions of the form $f\times g:(x,y)\mapsto f(x)g(y)$ where $f$ and $g$ are bounded continuous functions. In fact,

[TABLE]

The right-most term above is equal to

[TABLE]

Finally, the equality $\left\langle\mathcal{B}^{-}g,\overline{f}\right\rangle=\int_{\mathcal{D}}(f\times g)(x,y)\,d\mathcal{B}^{-}_{x}(y)\,d\nu(x)$ concludes the proof. ∎

Note that $\mathbb{P}_{[0,n]}$ and $\mathbb{P}^{-}_{[0,n]}$ are equivalent under the assumption that $\eta$ and $\eta^{-}$ are equivalent.

Observe that

[TABLE]

Also note that, for any function $f$ on $\mathcal{D}$

[TABLE]

From these observations we immediately obtain

[TABLE]

Proposition 14** (Entropy production rate).**

The entropy production rate for the random billiard system, under the assumption that the probabilities on pairs $\eta$ and $\eta^{-}$ are equivalent, takes the form

[TABLE]

In particular, this expression shows that $e_{p}\geq 0$.

Proof.

Due to Equation (3) we have $e_{p}=-\int_{\mathcal{D}}\log\left(\frac{d\eta^{-}}{d\eta\ }\right)\,d\eta$. Now observe that $\frac{d\eta^{-}}{d\eta\ }\circ\mathcal{R}=\frac{d\eta\ }{d\eta^{-}}.$ In fact, for any measurable set $E\subset\mathcal{D}$,

[TABLE]

Therefore,

[TABLE]

from which we conclude that $e_{p}=\frac{1}{2}\int_{\mathcal{D}}\left[d\eta-d\eta^{-}\right]\log\left(\frac{d\eta\ }{d\eta^{-}}\right)$ as claimed. It is apparent from this expression that $e_{p}\geq 0$. ∎

4 Second Law of Thermodynamics

The reciprocity property imposed on the reflection operator $P$, which is needed in order to make sense of the concept of boundary temperature, has not been used so far. Below, we rewrite the expression for $e_{p}$ obtained in the previous section, making use of this property. Before doing so, it is useful to introduce a generalized but natural notion of reciprocity as noted after Definition 3. This yields a more general notion of reflection operator that applies to manifolds having a local product structure, corresponding to billiard systems consisting of multiple rigid masses. In such cases we suppose the existence of a measurable family of subspaces $W_{q}\subset T_{q}M$ for each $q\in\partial M$ such that $\mathbbm{n}_{q}\in W_{q}$, and define the Maxwellian $\mu_{q}^{\pm}$ as in Equation (1), except that the dimension $n$ is now replaced with the dimension of $W_{q}$ and $dV_{q}$ is replaced with the volume measure on $W_{q}$. We then assume that the family of operators $P_{q}$ satisfies

$u\in W_{q}^{-}\mapsto\mathcal{P}(W_{q}^{+})$, where $W_{q}^{\pm}:=W_{q}\cap N^{\pm}_{q}$; 2. 2.

if $u\in T_{q}(\partial M)$ is perpendicular to $W_{q}$, then $P_{(q,u)}$ is the point mass at $u$; 3. 3.

reciprocity is defined for the measure $d\zeta_{q}(u,v)$ for $(u,v)\in W_{q}^{-}\times W_{q}^{+}$.

Definition 15**.**

The subspaces $W_{q}$ will be called directions of thermal contact or simply thermal directions. The boundary of $M$ is said to have temperature $T$ at $q\in\partial M$ if $P_{q}$ satisfies reciprocity with respect to the Maxwellian on $W_{q}^{-}$ having parameter $\beta=1/\kappa T$.

Thus if $u=u_{1}+u_{2}\in N_{q}^{-}$ is a pre-collision velocity decomposed into $u_{1}\in W_{q}^{-}$ and $u_{2}$ in the orthogonal complement $W_{q}^{\perp}$ of $W_{q}$, then the post-collision velocity is $U_{1}+u_{2}$ where $U_{1}$ is a random vector in $W^{+}_{q}$ distributed according to $P_{(q,u_{1})}$, and reciprocity holds with respect to a Maxwellian on the space of thermal directions.

A simple example will help to clarify the need for the above notion of thermal directions. Consider the system shown in Figure 2 describing two point masses $m_{1},m_{2}$ that can slide freely over an interval of length $l$. When the two masses collide with each other, their post-collision velocities are derived from the assumptions of conservation of kinetic energy and momentum and when they collide with the end-points of the interval, they reflect according to random reflection operators at temperatures $T_{1}$ and $T_{2}$.

The configuration space of the pair of masses is a right-triangle with the sides adjacent to the right angle having length $l$. A point $(x,y)$ represents the configuration in which $m_{1}$ is at $x$ and $m_{2}$ is at $y$. On the longer side are the configurations representing collisions of the two masses. Introducing new coordinates $x_{1}=\sqrt{\frac{m_{1}}{m}}x,\ \ \ x_{2}=\sqrt{\frac{m_{2}}{m}}y,$ where $m=m_{1}+m_{2}$, the total kinetic energy becomes a multiple of the square Euclidian norm, $|v|^{2}$ of the velocity vector $v=(\dot{x}_{1},\dot{x}_{2})$.

In this rescaled picture, the two-masses system becomes a random billiard system in which a point particle of mass $m=m_{1}+m_{2}$ moves freely inside the triangle and undergoes specular reflection on the hypotenuse while reflection on the shorter sides is random. On these sides the thermal direction is along the normal vector $\mathbbm{n}_{q}$ and their temperature is $T_{i}$, $i=1,2$.

As this example makes clear (see also the idealized heat engine given in Section 6), the moving particle in our definition of random billiards can be thought in general to represent the configuration of a system consisting of several moving rigid masses, possibly extended rigid bodies in dimension $3$. In such cases the configuration manifold can be a non-flat Riemannian manifold.

In preparation for the proof of Theorem 8 let us recall that the measures $\mu^{\pm}\in\mathcal{P}(N^{\pm})$ were defined as

[TABLE]

where $\bar{A}=A/A(\partial M)$ is the normalized Riemannian volume measure on $\partial M$. When the space $W_{q}$ of thermal directions is not all of $N_{q}$, we let $d\mu_{q}^{\pm}(u)=d\overline{\mu}_{q}(u_{1})\,dS^{\perp}_{q}(u_{2})$. This is the product measure of the Maxwellian along $W^{\pm}_{q}$ and the normalized volume measure on the hemisphere of radius $|u_{2}|_{q}$, where $u=u_{1}+u_{2}$ is the orthogonal decomposition of $u$ into its $W_{q}$ and $W_{q}^{\perp}$ components.

Proposition 16**.**

With the definitions from Section 3, and bringing into play the reciprocity property of the reflection operator $P$, we obtain

[TABLE]

where $x=(q,u)$.

Proof.

First observe that $\frac{d(\mathcal{J}_{*}\nu)}{d\mu^{+}}\circ\mathcal{J}=\frac{d\nu}{d(\mathcal{J}_{*}\mu^{+})}$. In fact, for any measurable subset $E\subset N^{+}$,

[TABLE]

The identity then follows from $\int_{E}d\nu=\int_{N^{+}}\frac{d\nu}{d(\mathcal{J}_{*}\mu^{+})}(x)\,d(\mathcal{J}_{*}\mu^{+})(x).$

Proceeding with the proof of Equation (4), we first recall that

[TABLE]

Reciprocity was defined by the relation

[TABLE]

(See Figure 4.) Using the measures $\mu^{\pm}$ on $N^{\pm}$ defined just prior to the statement of this proposition, we may rewrite the reciprocity property as

[TABLE]

Also notice that $d\mu^{\pm}(Jy)=d(J_{*}\mu^{\pm})(y)=d\mu^{\mp}(y)$ and $d\mu^{-}(\mathcal{T}(x))=d\mu^{+}(\mathcal{J}x)$. Therefore,

[TABLE]

This, in combination with the observation that began the proof, yields Equation (4). ∎

The factorization of the Radon-Nikodym derivative $d\eta/d\eta^{-}$ as a product of a function of $x$ and a function of $y$, as given in Proposition 16, allows us to express $e_{p}$ as an integral over $N^{+}$ rather than $\mathcal{D}$. This is indicated in the following proposition.

Proposition 17**.**

Define the function $\Lambda:N^{+}\rightarrow(0,\infty)$ given by $\Lambda(x)=\frac{d\nu}{d(\mathcal{J}_{*}\mu^{+})}(x)$. Then, given a stationary probability measure $\nu$ of the random billiard system,

[TABLE]

under the assumption that $\eta$ and $\eta^{-}$ are equivalent measures.

Proof.

Observe that for any $\nu$-integrable function $f$ on $N^{+}$,

[TABLE]

Using Proposition 16 and the general expression for the entropy production rate given in Proposition 14 we obtain

[TABLE]

where the term involving $\Lambda(x)$ vanished due to Equation (5). Also observe that

[TABLE]

and that

[TABLE]

Collecting the terms we obtain $e_{p}=\frac{1}{2}\nu\left(\log\left(\Lambda/\Lambda\circ\mathcal{J}\right)\right)$ as claimed. ∎

Let $m$ be the measure on $N^{+}$ defined by $dm(q,u)=dV_{q}(u)\,d\bar{A}(q)$ and $g(q,u)$ the density of $\nu$ with respect to $m$. Notice that

[TABLE]

where $E$ is the sum of the kinetic and potential energy functions: $E(q,u)=E_{0}(q,u)+\Phi(q)$, $E_{0}(q,u)=\frac{1}{2}m|u|^{2}_{q}$. Furthermore, invariance of $m$ under $\mathcal{J}$ gives

[TABLE]

where $h(x)=h(q,u)=C_{q}\exp\left\{-\beta(q)E(q,u)\right\}$. Define $l(x,y)=g(x)h(\mathcal{J}y)$ for $(x,y)\in\mathcal{D}$. These definitions give the expression

[TABLE]

where $l(x,y)=g(x)h(\mathcal{J}y)$. These definitions will be used in the proof of the Second Law.

Proof of Theorem 8.

The expression in the theorem statement could be derived by taking as starting point the expression for $e_{p}$ given in Proposition 17, but we prefer to begin with the general form for $e_{p}$ asserted in Proposition 14 and the function $l(x,y)$ appearing in (6). In fact, we start with the non-symmetrized expression for $e_{p}$ and use (in the fifth line) the property expressed in Equation (5):

[TABLE]

Equation (5), again, implies

[TABLE]

And stationarity of $\nu$ implies

[TABLE]

Therefore

[TABLE]

Now observe that $\mathcal{J}_{*}\nu=J_{*}\mathcal{T}_{*}\nu=J_{*}\nu^{-}$, which allows us to write

[TABLE]

Note that the kinetic energy function $E_{0}$ is invariant under $J$ and that

[TABLE]

The integral of the constant term $\log C_{q}$ against the difference of probability measures $\nu^{+}_{q}-J_{*}\nu^{-}_{q}$ gives [math], so we are left with

[TABLE]

Decomposing along the fibers of $\pi:N^{+}\rightarrow\partial M$ gives the desired expression. ∎

5 Multi-temperature Maxwell-Smoluchowski systems

The central purpose of this section is to study the entropy production rate for a system whose configuration space $M$ is a Riemannian manifold, with some mild regularity conditions, whose boundary is partitioned into components $\Gamma_{i}$ which are kept at temperatures $T_{i}$ respectively for $i=1,\ldots,N$. The notion of a thermostatted boundary is modeled using the Maxwell-Smoluchowski model, defined Definition 9. Recall that the model can be thought of as follows. Let constants $\alpha_{1},\ldots,\alpha_{N}\in(0,1]$ be given. Upon collision at any point of boundary component $i$, the post-collision velocity of the colliding particle is either chosen randomly according to the Maxwell-Boltzmann distribution with temperature $T_{i}$, with independent probability $\alpha_{i}$, or the particle reflects specularly with probability $1-\alpha_{i}$. When $\alpha_{i}$ is small, this model may be regarded as a random perturbation of an ordinary billiard system. More generally, it can be thought of as a model of thermalization, where the particle only takes on the temperature of the boundary thermostat after a geometrically distributed number of collisions.

The proof that the Maxwell-Smoluchowski model satisfies reciprocity amounts to the following elementary exercise. (In order to alleviate clutter, we omit the subscript $q$ from maps and measures.) For proving invariance of $\zeta$ under $\mathcal{R}$ it is sufficient to test the equality $\mathcal{R}_{*}\zeta=\zeta$ on functions of the form $(f\times g)(u,v)=f(u)g(v)$. For such functions,

[TABLE]

This last expression is seen to be [math] because $J_{*}\mu^{\pm}=\text{Ref}_{*}\mu^{\pm}=\mu^{\mp}$ and the flip and reflection maps commute. (Note that $\mu^{+}(f\circ J)=(J_{*}\mu^{+})(f)=\mu^{-}(f)$.)

5.1 Uniform ergodicity

Our present aim is to study the Markov chain $X_{n}$ on state space $N^{+}$ with Markov transition kernel $\mathcal{B}$. In the remainder of this section we restrict ourselves to the case in which $\mathcal{B}$ is induced by the Maxwell-Smoluchowski reflection operator, the billiard system is free of potential forces, and the bundle $W$ of thermal directions is all of $N$.

Before turning to the chain $X_{n}$, we introduce a related Markov chain derived from the projection of $X_{n}$ onto the bundle of unit vectors in $N^{+}$. Recall that $S\subset N$ denotes the bundle of unit vectors in $N$, $S^{+}=S\cap N^{+}$, and $S^{+}_{q}=S^{+}\cap N_{q}$. The hemisphere bundle $S^{+}$ is invariant under the standard billiard map, which is defined as the composition of the translation map $\mathcal{T}$ and specular reflection. This is clear since $\mathcal{T}$ and the specular reflection map preserve the Riemannian norm. The billiard measure $\sigma$ on $S^{+}$, introduced at the end of Section 2, is the measure invariant under the standard billiard map, obtained from the symplectic form as described, for example, in [5]. Recall that

[TABLE]

where $A$ is the Riemannian $(n-1)$-dimensional volume measure on $\partial M$ and $V_{q}^{\text{\tiny S}}$ is the Riemannian $(n-1)$-dimensional volume measure on $S_{q}^{+}$.

One property of the Maxwell-Smoluchowski reflection operator that should be highlighted is that it is projective according to the following definition. We denote by $\Pi_{q}:N_{q}^{+}\setminus\{0\}\rightarrow S_{q}^{+}$ the projection map $\Pi_{q}u=u/|u|_{q}$ and by $\Pi$ the corresponding projection map from $N^{+}$ minus the zero section onto $S^{+}$.

Definition 18** (Projective reflection operators).**

We say that the reflection operator $P$ is projective if for all nonzero $x=(q,u)\in N^{-}$, $\lambda>0$, and continuous $f:S_{q}^{+}\rightarrow\mathbb{R}$, the integral $P_{(q,\lambda u)}(f\circ\Pi_{q})$ does not depend on $\lambda$.

That the Maxwell-Smoluchowski model has this property is readily seen:

[TABLE]

but $\Pi_{q}\text{Ref}_{q}\lambda u=\Pi_{q}\text{Ref}_{q}u$, so the left-side of the equation does not depend on $\lambda$. Moreover, the associated billiard map $\mathcal{B}$ induces a map $\overline{\mathcal{B}}$ on $S^{+}$ as follows: given a continuous function $f$ on $S^{+}$,

[TABLE]

The operator $\overline{\mathcal{B}}$ thus acts as the Markov transition kernel for the Markov chain $\overline{X}_{n}:=\Pi(X_{n})$.

In [4], uniform ergodicity is established for a stochastic process related to $\overline{X}_{n}$, which informally can be described as follows. A particle moves with constant speed inside the domain $M$. Upon collision with the boundary, it is reflected in some random direction, not depending on the incoming direction, keeping the magnitude of its velocity constant. The distribution of the random direction is given by the so-called Knudsen cosine law, which is the velocity component of the billiard measure $\sigma$. The particle then continues on to the next collision point, where it is again reflected in a random direction, independent of the previous collision, and this continues ad infinitum. The sequence of collision points, which forms a Markov chain with state space $\partial M$, is referred to as the Knudsen random walk.

It is readily seen that the Knudsen random walk is precisely the process $\xi_{n}:=\pi(\overline{X}_{n})$ when $\alpha(q)\equiv 1$. The following conditions on $\partial M$ are adapted from [4].

Assumption 19**.**

Suppose that $\partial M$ is an $(n-1)$-dimensional, almost everywhere continuously differentiable surface satisfying the following Lipschitz condition. For each $q\in\partial M$, there exists $\epsilon_{q}>0$, an affine isometry $\mathcal{I}_{q}:M\to\mathbb{R}^{n}$, and a function $f_{q}:\mathbb{R}^{n-1}\to\mathbb{R}$ such that

• •

The function $f_{q}$ satisfies $f_{q}(0)=0$ and the Lipschitz condition. That is, there exists a constant $L_{q}>0$ such that $|f_{q}(p)-f_{q}(p^{\prime})|<L_{q}|p-p^{\prime}|$ for all $p,p^{\prime}\in\partial M$.

• •

The affine isometry satisfies $\mathcal{I}_{q}q=0$ and

[TABLE]

Theorem 20** (adapted from [4]).**

Suppose $\mathrm{diam}(M)<\infty$ and Assumption 19 holds. Then the normalized Riemannian volume $\bar{A}$ on $\partial M$ is the unique stationary distribution for the Knudsen random walk. Moreover, there exist constants $\beta_{0},\beta_{1}$, independent of the distribution of $\xi_{0}$, such that

[TABLE]

where $||\cdot||_{v}$ is the total variation norm.

In what follows, we extend the result above to show that the Markov chains $\overline{X}_{n}$ and $X_{n}$ are uniformly ergodic. Moreover, we give explicit expressions for the stationary measures of these chains. Before stating the theorem, we first establish some notation. Recall that we assume that $\partial M$ is partitioned into $N$ components $\Gamma_{i}$. Moreover $T(q)\equiv T_{i}$ and $\alpha(q)\equiv\alpha_{i}$ for all $q\in\Gamma_{i}$. Let

[TABLE]

be the one-step transition probability of the Knudsen random walk between components of $\partial M$. Let $A_{i}=A(\Gamma_{i})$ and $\bar{A}_{i}=\bar{A}(\Gamma_{i})$ be the volume measure and normalized volume measure of the boundary components $\Gamma_{i}$ respectively. Next, note that since the temperature is constant on boundary components, by identifying $N_{q}^{+}$ with the upper half space $\mathbb{R}^{n}_{+}=\{x\in\mathbb{R}^{n}:x\cdot e_{n}>0\}$, we can define $\mu_{i}:=\mu_{q}\in\mathcal{P}(N_{q}^{+})$ for $q\in\Gamma_{i}$ to be the Maxwellian associated to $\Gamma_{i}$. Moreover, it is readily apparent that $\mu_{i}$ can be disintegrated using polar coordinates as a product measure on $S_{q}^{+}\times(0,\infty)$. One can check that the component on $S_{q}^{+}$ is $d\sigma_{q}(u)=C_{q}\left\langle u,\mathbbm{n}_{q}\right\rangle\,dV_{q}^{\text{\tiny S}}(u)$, the Knudsen cosine law referred to above. We denote the component on $(0,\infty)$, the speed component, as $\mu_{i}^{\mathrm{s}}$:

[TABLE]

Finally, let $Q$ be the $N\times N$ matrix and $\pi$ the $N$-dimensional vector where

[TABLE]

Theorem 21**.**

Suppose $\mathrm{diam}(M)<\infty$ and that Assumption 19 holds. Then

The billiard measure $\sigma$ is stationary for $\overline{X}_{n}$. 2. 2.

There exists a unique stationary distribution $\nu$ for $X_{n}$. Moreover, there exist constants $b_{0},b_{1}$, independent of the distribution of $X_{0}$, such that

[TABLE]

where $||\cdot||_{v}$ is the total variation norm. That is, the chain $X_{n}$ is uniformly ergodic. 3. 3.

When the boundary temperatures $T_{1},\ldots,T_{N}$ are equal, say to a constant $T_{0}>0$, the stationary distribution $\nu$ is given by

[TABLE]

where $\rho_{q}(v)$ is the Maxwellian density given in (1) with constant temperature $T(q)\equiv T_{0}$. 4. 4.

When the boundary temperatures are not equal, the stationary distribution $\nu$, expressed in polar coordinates as a measure on $S^{+}\times(0,\infty)$ is given by $d\nu(q,u,r)=d\nu^{\mathrm{s}}_{q}(r)d\sigma(q,u)$, where $\nu^{\mathrm{s}}_{q}$ is the stationary distribution of the speed after collision with boundary point $q$. The measure $\nu_{q}^{\mathrm{s}}$ is constant on components $\Gamma_{i}\subset\partial M$, so we let $\nu_{i}^{\mathrm{s}}:=\nu_{q}^{\mathrm{s}}$ for $q\in\Gamma_{i}$. Letting $\nu^{\mathrm{s}}$ be the $N$-dimensional vector with components $\nu_{i}^{\mathrm{s}}$, we have that

[TABLE]

Proof.

Let $\mu$ be the measure derived from the Maxwellians $\mu_{q}$ as indicated just prior to the statement of Proposition 16. Note that $\overline{\mu}:=\Pi_{*}\mu=\sigma$. Moreover, it is a consequence of the reciprocity property of $P$ that the following time reversibility condition holds: for all $x,y\in N^{+}$,

[TABLE]

where $\overline{x}:=\Pi x$. It then follows that $\sigma\overline{\mathcal{B}}=\sigma$ and the first part of the theorem is proved. Moreover, the same argument can be used to show that the third part of the theorem holds.

For the second part of the theorem, we show that the chain $X_{n}$ is uniformly ergodic by showing that the state space $N^{+}$ is small in the following sense: there exists $m\in\mathbb{Z}_{+}$ and a nontrivial measure $\phi$ on $N^{+}$ (which is not necessarily a probability measure) such that for all $x\in N^{+}$ and measurable sets $A\subseteq N^{+}$, we have that $\mathcal{B}^{m}_{x}(A)\geq\phi(A)$ (see Theorem 16.0.2 of [23]). By Theorem 20, there exists $m\in\mathbb{Z}_{+}$ and a nontrivial measure $\overline{\phi}$ on $S^{+}$ so that the aforementioned condition for uniform ergodicity holds for the $\overline{X}_{n}$ chain. We can then construct the measure $\phi$ as follows:

[TABLE]

where $L_{x}(A):=\{|v|:\overline{x}=(q,v/|v|)\in\Pi(A)\}$ and $\rho_{q}^{\text{\tiny s}}(|v|)=\int_{S_{q}^{+}}|v|^{n-1}\rho_{q}(v)dV_{q}^{\text{\tiny S}}(v/|v|)$. Note that

[TABLE]

For the final part of the theorem, let $\nu_{0}$ be the initial distribution of the Markov chain and define $\nu_{k}=\nu_{0}\mathcal{B}^{k}$ to be the distribution of $X_{k}$. Further, note that the state space $N^{+}$ can be partitioned into components $N_{i}^{+}=\{x=(q,u)\in N^{+}:q\in\Gamma_{i}\}$ and we let $\nu_{k,i}$ denote the restriction of $\nu_{k}$ to $N_{i}^{+}$. Note that the stationary measure for $X_{n}$ must assign probability $\bar{A}_{i}$ to component $N_{i}^{+}$ since the normalized area measure on the boundary is stationary for $\overline{X}_{n}$. Moreover, the relation $\nu_{k+1}=\nu_{k}\mathcal{B}$ amounts to the system

[TABLE]

for $i=1,\ldots,N$. Expressed in matrix form, this yields $\nu_{k+1}=\pi+Q\nu_{k}$, which implies

[TABLE]

It is straightforward to check that this converges in total variation to $(I-Q)^{-1}\pi$.

∎

5.2 Entropy production rate formula

With an explicit expression for the stationary measure in hand, we are now ready to give a formula for the entropy production rate. It is clear that the hypotheses of Theorem 8 are satisfied by the stationary measure given in Theorem 21. The main work will be in computing $\nu_{q}^{+}(E_{0})-\nu_{q}^{-}(E_{0})$, where $E_{0}$ is the kinetic energy function.

Proposition 22**.**

The entropy production rate in the multi-temperature Maxwell-Smoluchowski system is given by

[TABLE]

where $\nu^{+}$ is the stationary measure $\nu$ given in Theorem 21, $\nu^{-}=\mathcal{T}_{*}\nu$, and $\nu_{j}^{\pm}$ is the restriction of the stationary measure $\nu^{\pm}$ to component $N_{j}^{\pm}=\{(q,u)\in N^{\pm}:q\in\Gamma_{i}\}$.

The factor $\nu_{j}^{+}(E_{0})-\nu_{i}^{+}(E_{0})$, which quantifies the stationary state energy exchange between wall $j$ and wall $i$, can be expressed more explicitly in terms the temperature gradient between wall $j$ and wall $i$, but we leave this to some explicit examples of tables below.

Proof.

Let $E_{0}(q,u)=\frac{1}{2}m|u|_{q}^{2}$. Observe that

[TABLE]

Next, note that $\nu_{j}^{-}(E_{0})=\sum_{i=1}^{N}p_{ij}A_{i}/A_{j}\nu_{i}^{+}(E_{0})$. Using the entropy production rate formula in Theorem 8 this gives

[TABLE]

where the final equality follows because $\sum_{i=1}^{N}p_{ij}A_{i}/A_{j}=P(\xi_{n}\in\Gamma_{i}\mid\xi_{n+1}\in\Gamma_{j})=1$ for each $j=1,\ldots,N$, when $\xi_{n}$ is the stationary Knudsen random walk. ∎

5.3 Examples

5.3.1 The two-plates system

We illustrate the formula giving the entropy production rate for the elementary system indicated in Figure 5. It consists of a particle that bounces back and forth between two parallel plates kept at temperatures $T_{1}$ and $T_{2}$. For the reflection operator we adopt the Maxwell-Smoluchowski model with parameters $\alpha_{1}$ and $\alpha_{2}$. Thus at any point $q$ of plate $i$, the post-collision velocity of the colliding particle has the Maxwell-Boltzmann distribution with temperature $T_{i}$ with probability $\alpha_{i}$, and with probability $1-\alpha_{i}$ the particle reflects specularly.

The manifold is taken to be the product of a flat torus and an interval, $M=\mathbb{T}^{2}\times[0,l]$, where $l$ is the distance between the two (torus) plates, which comprise the two connected components of the boundary $\partial M=\mathbb{T}^{2}\times\{0,l\}$. Let the index $1$ be associated with the left-side plate and $2$ with right-side plate. The phase space $N$ is the union of components $N^{\pm}_{i}$, where $i\in\{1,2\}$. Each $N^{\pm}_{i}$ is identified with $\mathbb{R}^{3}_{+}=\{u\in\mathbb{R}^{3}:u\cdot\mathbbm{n}>0\}$, where $\mathbbm{n}$ is the normal vector to plate $1$. We indicate by $\bar{i}$ the index opposite to $i$, so $\bar{1}=2$ and $\bar{2}=1$. Then the translation map assumes the form $\mathcal{T}(i,u)=(-\bar{i},u)$. We write $f_{i}(u)=f(i,u)$ where $u\in\mathbb{R}^{3}_{+}$

The billiard map applied to a function $f$ on $N$ has the form

[TABLE]

Because the billiard particle alternates between the two plates, a stationary measure for $\mathcal{B}$ must assign equal probabilities for each plate: $\nu(N_{i})=1/2$. We can then restrict attention to the two-step chain describing the sequence of returns to a plate. Let $\nu_{0}$ be an initial distribution for the Markov chain and define $\nu_{k}=\mathcal{B}^{k}\nu_{0}$. Let $\nu_{k,i}$ be the restriction of these measures to plate $i$. Then the equation $\nu_{k+1}=\nu_{k}\mathcal{B}$ amounts to the system

[TABLE]

Writing

[TABLE]

equation $\nu_{k+1}=\nu_{k}\mathcal{B}$ becomes $\nu_{k+1}=\nu_{k}Q+\pi$, from which we easily obtain

[TABLE]

where $\gamma=(1-\alpha_{1})(1-\alpha_{2})$. We conclude that the process converges to a stationary process having stationary distribution $\nu=\nu_{1}+\nu_{2}$ where

[TABLE]

where $c=1-\gamma=1-(1-\alpha_{1})(1-\alpha_{2})$. Recall that $\nu^{-}=\mathcal{T}_{*}\nu$. Thus $\nu_{i}^{-}=\nu_{\bar{i}}$:

[TABLE]

The stationary measure clearly satisfies the conditions of Theorem 8. Moreover

[TABLE]

A simple computation gives $\mu_{i}(E_{0})=kT_{i}$. The entropy formula given in Theorem 8 then yields the result:

[TABLE]

This expression is clearly non-negative. Let us denote by $Q_{i}$ the expected change in energy of the billiard particle at a collision with plate $i$ in the stationary regime. Then

[TABLE]

Thus if $T_{1}>T_{2}$, the particle takes, on average, an amount $Q:=\frac{\alpha_{1}\alpha_{2}}{2c}\left(\kappa T_{1}-\kappa T_{2}\right)$ of energy from plate $1$ and transfers it to plate $2$ per collision. Combining with the expression for $e_{p}$ results in

[TABLE]

It should be noticed that this expression is independent of the nature of the thermostat model assumed for the plates. On the other hand, the above expression showing that $Q$ is proportional to the temperature difference depends on the choice of model. In fact, the coefficient $\frac{\alpha_{1}\alpha_{2}}{2\left[1-(1-\alpha_{1})(1-\alpha_{2})\right]}$ tells how fast the system transfers energy from hot to cold plate, and may thus be interpreted as a heat conductivity parameter (measured per pair of collisions rather than time between collisions; the latter may also be calculated from the stationary measure).

5.3.2 A three-temperature system

Next, we wish to express the entropy production rate for systems with more than two temperatures. As a prototypical example we take $\partial M$ to be an equilateral triangle where boundary component $\Gamma_{i}\subset\partial M$ is equipped with parameters $T_{i}$ and $\alpha_{i}$ for $i=1,2,3$. Note that the calculations and formulas to be shown can be generalized to more general polygons, but we restrict our attention to the equilateral triangle for the sake of simplicity.

Following the notation in Subsection 5.1 and Equation (7), the side lengths $A_{i}$ are equal for all $i$ and the normalized side length $\bar{A}_{i}=1/3$. Moreover, by symmetry $p_{ij}=(1-\delta_{ij})/2$ for $i,j=1,2,3$, where $\delta_{ij}$ denotes the Kronecker delta. Using Theorem 21,

[TABLE]

and $Q_{jk}=(1-\alpha_{j})(1-\delta_{jk})/2$. As in the two-plates example, we let $Q_{i}$ be the expected change in energy of the billiard particle at a collision with boundary component $i$ in the stationary regime. Moreover, let $\overline{Q}_{ij}$ denote the average amount of energy from component $i$ which is transferred to component $j$ per collision.

A tedious but straightforward calculation gives

[TABLE]

for $i=1,2,3$, where $i^{\prime}=i+1\mod 3$ and $i^{\prime\prime}=i+2\mod 3$. Note that another simple computation yields that $\mu_{i}(E_{0})=\frac{3}{2^{3/2}}\kappa T_{i}$ for a two-dimensional billiard domain. Using the entropy formula, we have that

[TABLE]

We have expressed $Q_{i}$ as the sum of two differences $\overline{Q}_{ij}$ in order to emphasize that the expected change in energy involves a transfer of energy among pairs of boundary components. While the formulas are a bit more complicated in the case of multiple temperatures, they nevertheless capture the general qualitative property that the entropy production scales with the square of the temperature difference. On the other hand, the coefficients on the terms $\mu_{i}(E_{0})$ are specific to the model and express how fast the system transfers energy among boundary components.

5.3.3 A circular chamber system

We conclude this section with a numerical example to demonstrate how geometric features of the billiard table, as opposed to features of the collision model, can influence the entropy production rate.

The inset in Figure 6 shows the billiard table of interest. It consists of two overlapping discs of radius $r$ with centers at a distance $a$ apart. We call $a/2r$ the ratio parameter. The boundary of the table is the union of two symmetric arcs of circles kept at constant temperatures $T_{1}$ and $T_{2}$ and equipped with the Maxwell-Smoluchowski collision model. Note that when $a/2r=0$, the two discs coincide and the boundary segments are the right and left hemispheres; and when $a/2r=1$, the discs are in tangential contact. We are interested in the changes in entropy production rate $e_{p}$ due to varying the ratio parameter over the interval $(0,1)$.

Using the formula for the stationary measure in Theorem 21 and the formula for entropy production in Proposition 22, the rate $e_{p}$ is determined by the transition probabilities $p_{ij}$ between boundary segments of $\partial M$. These probabilities are then estimated through numerical simulation of the billiard dynamic by doing a large sampling from the phase space of initial conditions and keeping a tally of one-step transitions between boundary segments. That is, $p_{ij}$ is estimated to be the proportion of trajectories from the sample that start at boundary segment $i$, i.e. circular arc $i$, and then have a first collision somewhere along boundary segment $j$, where $i,j=1,2$.

The five graphs of Figure 6 correspond to $T_{2}-T_{1}=0,\ldots,4$. (These values are indicated on top of each graph.) Each curve in the graph gives values of $e_{p}$ for 40 values of the ratio parameter. An interesting feature exhibited by the graphs in Figure 6 is the sharp transition in the rate of decay of $e_{p}$ past a value of $a/2r$ roughly between $0.8$ and $1.0$ across different values of $T_{2}-T_{1}$.

6 Work production

We view this paper’s focus on entropy production and the second law as only a first step in a broader investigation of the stochastic thermodynamics for random billiards. It is natural to ask whether these systems can be used to explore a wider range of classical thermodynamic phenomena. In this final section we wish to show how the issue of work production can very naturally be brought into the general picture by exploring, mainly numerically, a random billiard model of a heat engine.

The system we use to illustrate this point is shown in Figure 7. It is only one extremely simple example of a general class of models that we call thermophoretic engines, to be considered more systematically in a future study. For now, this example will serve to show the possibilities offered by random billiard models in stochastic thermodynamics beyond the more restricted purview of the present paper.

In Figure 7, a point mass $m_{2}$ can move freely inside an equilateral triangular domain whose sides are of two types: two of them can thermally interact with the particle through the Maxwell-Smoluchowski thermostat model with temperatures $T_{h}$ and $T_{c}$; the third side at the top of the triangle can slide without friction as a kind of conveyor belt. This side has its own mass, $m_{1}$. A constant force $F$ is applied to $m_{1}$. Among the many possible choices for the model of interaction between the two masses we choose here the simplest, and not necessarily the most physically natural, that allows for an exchange of momentum: when $m_{2}$ collides with the top side of the triangle the perpendicular component of its velocity changes sign and the horizontal component, together with the velocity of the belt, change as in the two-masses example of Figure 2. That is, the interaction is as if $m_{1}$ and $m_{2}$ undergo a frontal elastic collision in dimension $1$.

We first describe how this system fits the definition of a random billiard. Let $L$ denote the length of the conveyor belt and let $\mathcal{S}$ be the interior of the triangle. The configuration manifold is then the $3$-dimensional space $M=\mathcal{S}\times\mathbb{R}/L\mathbb{Z}$. As in the two-masses system of Section 4, it is convenient to rescale the position coordinates of $L$ and $\mathcal{S}$ so as to make the total kinetic energy proportional to the square Euclidean norm of the combined velocity vector of the two masses. In other words, if $(x,y)$ are the coordinates on the plane of $\mathcal{S}$ and $z$ parametrizes the position of the conveyor belt, we set

[TABLE]

where $m=m_{1}+m_{2}$. Let us also denote $\gamma:=\sqrt{m_{2}/m_{1}}$. We introduce a positive orthonormal frame on the boundary of $M$ denoted $(e_{1},e_{2},e_{3})$ where $e_{1}$ points in the direction of the axis $x_{1}$, $e_{2}$ is tangent to the boundary of $M$, and $e_{3}$ is perpendicular to the boundary, pointing to the interior of $M$. Let $v$ and $V$ denote velocity vectors before and after a collision, respectively, expressed in the coordinate system given by $x=(x_{1},x_{2},x_{3})$. Collisions with the stationary sides of the triangle are assumed to satisfy the reciprocity condition such that at each $x\in\partial M$ the subspace of thermal directions, in the sense of Definition 15, is perpendicular to $e_{1}(x)$. Collisions with the sliding top side are deterministic. Conservation of energy and momentum implies that $V=Cv$ where the collision map $C\in SO(3)$ is the orthogonal involution given by

[TABLE]

Note that the restriction of $C$ to $T_{x}\partial M$ is a reflection. The orthogonal component of $v$ in the subspace spanned by $\gamma e_{1}-e_{2}$ and $e_{3}$ changes sign, whereas the component parallel to $e_{1}+\gamma e_{2}$ remains unchanged.

Observe how the particle of mass $m_{2}$ imparts momentum to the sliding side of the triangle as it mediates the energy transfer between the two thermal sides. In the absence of the force $F$, we should expect the conveyor belt to move with a steady drift conterclockwise if $T_{h}>T_{c}$, and move in the opposite direction when the temperatures are reversed. We should also expect, for some range of values of $F$ and $T_{h}-T_{c}$, to observe work being produced against the force. This is in fact what is obtained by the numerical simulation. Figure 8 shows the motion of the conveyor belt as a function of time for different values of $T_{h}-T_{c}$.

As expected, the stochastic motion of the conveyor belt exhibits a steady drift compatible with an energy flow from the hot side to cold. The middle graph among the seven described in Figure 8 corresponds to $T_{h}=T_{c}$ and it is also shown in a different scale in the inset figure.

If we now impose a constant external force $F$ on $m_{1}$, a fraction of the energy flow between the thermal sides is converted into work against $F$. More precisely, at each time $t>0$, let $Q_{h}(t)$ be the total amount of energy transferred to the system since $t=0$ due to collisions between the particle and the hot wall. Let $Q_{c}(t)$ be the energy transferred due to collisions with the cold wall. Let $x_{w}(t)$ be the signed displacement from the initial position $x_{w}(0)$ of the mass $m_{1}$ at time $t$ and let $W(t)=(x_{w}(t)-x_{w}(0))F$ be the work done by the system. The Carnot mean efficiency of the engine up to time $t$ is then

[TABLE]

Let us now take into account, in standard fashion, energy conservation. The internal energy of the system at time $t$ is $E(t)=E_{w}(t)+E_{p}(t)$, which consists of the kinetic energy $E_{w}(t)$ of mass $m_{1}$ and the kinetic energy $E_{p}(t)$ of the particle $m_{2}$. Then conservation of energy, or the First Law of Thermodynamics, gives at each $t$

[TABLE]

As $Q_{h}(t)$ is the accumulated energy that flows through the hot side up to time $t$, and since $E(t)$ can be expected to have a finite mean value over time, it should be the case (and is seen experimentally) that $(E(t)-E(0))/Q_{h}(t)\to 0$ for large values of $t$. This then yields the steady state expression of efficiency given by

[TABLE]

Figure 9 shows the characteristic curve of mean efficiency as a function of $F$. The simulation evaluates, for each value of the force, the efficiency over $5000$ sample runs, with each run of length corresponding to $1000$ collision events. The parameters are $T_{h}=50,T_{c}=1,m_{1}=1000,m_{2}=1.$ The vertical bars indicate $95\%$ confidence intervals.