Fourier's Law in a Generalized Piston Model

A simplified, but non trivial, mechanical model -- gas of $N$ particles of mass $m$ in a box partitioned by $n$ mobile adiabatic walls of mass $M$ -- interacting with two thermal baths at different temperatures, is discussed in the framework of kinetic theory. Following an approach due to Smoluchowski, from an analysis of the collisions particles/walls, we derive the values of the main thermodynamic quantities for the stationary non-equilibrium states. The results are compared with extensive numerical simulations; in the limit of large $n$, $mN/M\gg 1$ and $m/M \ll 1$, we find a good approximation of Fourier's law.

velocities v i . The particles interact with the pistons via elastic collisions: where primes denote post-collisional velocities. The two external walls in z = 0 and z = L act as thermal baths at temperatures T A and T B . The interaction of the thermostats with the particles is the following: when a particle collides with the wall, it is reinjected into the system with a new velocity drawn from the probability distribution [22] with + for the case A and − for B, and k B is the Boltzmann constant, Θ(x) being the Heaviside step function: For the following discussion, we define the temperature of the i−th piston as and the temperature of the particle gas in the compartment j as an average on the particles between the (j − 1)-th and the j-th piston

III. SIMPLE THERMODYNAMIC CONSIDERATIONS
We expect, and this is fully in agreement with the numerical computations, that, given a generic initial condition, after a certain transient, the system reaches a stationary state. The positions of the walls fluctuate around their mean values z j , in a similar way to T (p) i and T i . The first non-trivial problem of the non equilibrium statistical mechanics is to determine z j , T (p) i and T i as function of the parameters of the model, i.e., n, L, m, M, T A and T B . We first analyze the simplest case of a single piston, where thermodynamics is sufficient to provide a complete description of the stationary state. Then, we consider the more general case of a multiple piston; in such a situation, thermodynamic relations are not enough to univocally determine the steady state: it is necessary to adopt a statistical mechanics approach. Such an approach will rely on three main assumptions, discussed in more detail below: small mass ratio m/M 1, local thermodynamic equilibrium in each compartment, and independence of the collisions particles/pistons.

A. Single Piston
If n = 1, using the equation for the perfect gas in each compartment, we immediately get the equations where p is the pressure, yielding Therefore, in this case, thermodynamics univocally determines the stationary state of the system.

B. Multiple Piston
We now consider the generalized piston model with n > 1. An analysis of the case n = 2 is enough to understand how the relations obtained from thermodynamics can be not sufficient to fully characterize the non equilibrium steady state. Indeed, we have the following relations: which give the constraints Therefore, since we have three variables (z 1 , z 2 , T 1 ) and only two constraints, thermodynamics is not enough to determine the stationary state. The computation can be easily extended to an arbitrary value of n, leading to that give the conditions with m = 2, . . . , n − 1, namely, only n constraints for 2n − 1 variables.

IV. COARSE-GRAINED DESCRIPTION AND EFFECTIVE LANGEVIN EQUATIONS
In order to obtain a statistical description of the system, we now derive effective stochastic equations, governing the dynamics of the relevant variables. Previous theoretical studies based on a Boltzmann equation approach on similar systems were reported in [13][14][15]. In particular, a generalized piston model was considered in Reference [15]. In those works, theoretical results were not compared to numerical simulations, so that the underlying hypotheses and their range of validity remained unclear.
Here, we present a different analysis. We will assume that the evolution of the observables is described by Langevin equations. The idea, coming back to Smoluchowski, is to integrate out the fast degrees of freedom of the gas particles by computing conditional averages, knowing the macroscopic variables: position z and velocity V of each piston, and temperatures T of the gases. In order to simplify the notation, let us denote by · this average, meaning the conditional average ·|z, V, T . We will compute the average change of a generic observable X in a small time interval ∆t, due to the collisions between the particles of the gas and the pistons.
Let us assume that in the stationary state the particles of the gas, in each compartment, have uniform space distribution ρ and a Maxwell-Boltzmann distribution φ(v) at temperature T : where ∆z is the length of the box containing the gas. We can obtain the rate of the collisions by considering the following equivalent problem: piston at rest and a particle, which moves with the relative velocity v − V . The point particles which collide against the piston in x in the time interval dt are: respectively, for particles on the left and on the right with respect to the piston. The Heaviside function Θ is necessary to have a collision.
Let us now consider a generic observable X j , depending in general on the velocities of the gas particles and of the pistons. We want to write down a Langevin equation: where both the drift term D j (X) (X being the vector of all relevant macroscopic variables in the system) and the noise term are due to collisions with the particles of the gas. We have: where ∆X j R coll is the variation of the observable due to the collisions with the particles that have the pistons on the left, while ∆X j L coll denotes the variation due to the collisions with the particles which have the piston on the right. In the stationary state, the macroscopic variables T , V and z do not depend on time. It means that the time derivative of the conditional average of a generic observable X is zero, namely: By using a perturbation development in the small parameter = m/M , it is possible to derive the relations between the average positions and the temperatures of the pistons and the temperatures of the gas. In order to obtain this result, we need to average on the piston velocity, which is a stochastic variable. Of course, in the steady state, the odd momenta have to be zero: The idea is to compute the average force produced by the particles of the gas, which collide against the piston, by computing the average momentum exchanged in the collisions. The evolution of the velocity of the piston is governed by the following stochastic differential equation: where η is zero-average white noise, usually K is a constant, and F coll is the average force which acts on the piston. This force is due to the collisions with the gas on the left and on the right, and depends on the average size of the box on the right and on the left of the piston, ∆z L and ∆z R , respectively, and on the temperatures of the gas on the left and on the right, T L and T R . Now let us set ∆X = V − V in Equation (15). As detailed in the Appendix, computing the left and right contributions to F coll and expanding in the small mass ratio = m/M , we obtain at the lowest orders: In the steady state, when the time derivative vanishes, we have from the order O(1) the following relation between the temperatures of the gas and the average length of the boxes: This relation is nothing but the one that can be derived from thermodynamics. From Equation (22), we obtain the obvious result V = 0. Analogous computations, reported in the Appendix, can be carried out for the variance of the pistons and the gas temperatures. In particular, for the variable ∆X = M (V 2 − V 2 ), at lowest orders in , we obtain the relations: In the steady state, by using the thermodynamic relations, the order O(1) Equation (24) is identically zero. After integrating over all values of the piston velocity and by requiring that the order O( ) vanishes, we obtain a relation between the temperature of the piston T (p) and those of the gases on the left and on the right: Finally, a condition on the temperature of the gas as a function of the velocity variance of the pistons on its left and on its right can be obtained by considering the variable ∆X = m(v 2 − v 2 )/N in Equation (15). This gives (see Appendix) the relations: By integrating over the velocity of the piston, because V is zero, the order O(1) is identically zero. By requiring that, in the stationary state, the order O( ) vanishes, we eventually obtain: As expected from thermodynamics, one can easily show that the gases in contact with the thermostats reach the temperatures of the thermal baths (see Appendix for a detailed explanation).

V. NUMERICAL SIMULATIONS
In order to check the range of validity of the above approach, we have performed extensive molecular dynamics simulations (using an event driven algorithm) of the system, varying the relevant parameters of the model and comparing the analytical prediction of Section 3 with the actual numerical results. The system is initialized with the pistons placed at equidistant positions, with zero velocity, while the gas particles are randomly distributed in each box, with velocities drawn from a Gaussian distribution at temperature intermediate with respect to those of the thermostats. We have, however, checked that the stationary state reached by the system is independent of the initial conditions. All data presented in the following are measured in the steady state, after the transient relaxation dynamics from the initial state.  In Figure 2, we show the piston temperatures T In Figure 3, we consider a system with a large number of pistons (n = 22), where the temperature profile shows a linear behavior, in agreement with Fourier's law. We report two sets of data, differing in the value of the parameter R = mN/M . Notice that the linear behavior extends for a wider range of z in the case R = 10. Lines represent linear fits of the data. Let us note that, for large n and R 1, T (p) vs. z is in good agreement with a linear behavior in the whole space interval (even close to the thermal baths). The derivation of the Langevin equation (see Appendix) is mainly based on the assumptions: • m/M 1; • a local thermodynamic equilibrium, i.e. in each compartment between the piston i − 1 and the piston i, the spatial distribution of the particles is uniform and the velocity probability distribution is a Gaussian function, whose variance is given by the gas temperature; • the collisions particles/pistons are independent.
The first assumption is easily checked, while the other ones are more difficult. We can expect that a necessary condition for their validity is that N must be large. We have measured the velocity distribution of the gas particles in the numerical simulations, finding that the above hypothesis is verified. Moreover, we expect that the number of recollisions decreases with the ratio m/M , and we have numerically checked that their contribution is negligible (the fraction of recollisions on the total number of collisions is about 0.1%). Figure 4 shows the probability distribution functions of particles colliding from left, φ L (v), and from right, φ R (v), with a moving wall: the agreement with a Gaussian assumption is rather evident.
In addition, the numerical computations show that the left-and right-moving particles in the same compartment have the following probability distributions: where T − (T + ) denotes the temperature of the particles incident on the piston i (i − 1) from the left (right), with T − ≈ T + , resulting in a small but finite heat flow proportional to (T − − T + ), in agreement with [15]. Let us note that the above probability distributions are different from the distributions φ L (v) and φ R (v), which are computed for particles actually colliding with the piston; this is the origin of the presence of the factor v appearing in the expressions in the caption of Figure 4.  The distribution of v in the compartment between pistons i − 1 and i, has the form: where A and B are constants. However, since we have a Gaussian distribution for the particle gas in the compartment: Let us note that all our results are based on an expansion in powers of , neglecting O( 2 ). In other words, because the heat rate is very small, this heat flux perturbs the Gaussian form of the probability distribution φ in a negligible way.
In order to better understand the dependence of our theoretical results on the parameter R, we have performed numerical simulations of the system for different values of R. From the previous remark, we expect to have an improvement of the agreement between the numerical results and the analytical predictions by increasing R. In Figure 5, we compare the piston temperature in a 4-piston system with theoretical predictions from Equation (26), finding a very good agreement for large values of R. This shows that the total mass of the gas contained in a box has to be larger than the mass of the piston, in order for the kinetic theory presented in the previous section to be accurate.

VI. CONCLUSIONS
We have studied a generalized piston in contact with two thermal baths at different temperatures. This system represents a simple but interesting case where the emergence of Fourier's law from a microscopic mechanical model can be studied. We have presented a kinetic theory treatment inspired by an approachà la Smoluchowski, and we have investigated the range of validity of these results with molecular dynamics numerical simulations. We have found that, in order for the theory to be accurate, the ratio R = mN/M should be large enough, namely the total mass of the gas in each compartment should be greater than that of the single piston.
We have considered ideal gas in our model, but we do not expect that the introduction of short-range interactions among gas particles, at least in not too dense cases, leads to significant differences in the behavior of the system. An interesting, non-trivial, future line of research in this model is the study of the relaxation to the steady state and the dynamical properties of the system, focusing on correlation and response functions.

Appendix A
In this appendix, we will derive the relations that determine the stationary state of the multi-component piston model described in the main text.

Piston Position
The effective force acting on the piston due to the collisions with gas particles, appearing in Equation (20), is given by two contributions, F coll = F L coll + F R coll . By taking into account the elastic collisions rule, Equation (1), these terms can be computed as follows: where in the first line we have used Equation (12). In the same way, by putting the observable ∆X = M (V − V ) in the Equation (17), we have: By putting together the previous relations, we obtain: By expanding in the small ratio = m/M 1, and assuming that V /v ∼ , we have: In the steady state, when the time derivative vanishes, we have from the order O(1) the following relation between the temperatures of the gas and the average lengths of the boxes:

Piston Fluctuations
Let us now consider the observable ∆X = M (V 2 − V 2 ). From Equation (16), using Equation (1) for the elastic collisions, we obtain: Using Equation (12), we have: In the same way, from Equation (17), we obtain: Putting together the contributions from the relations (A8) and (A9), we have By using a Taylor expansion in the small parameter = m M 1, we obtain: In the steady state, by using the thermodynamic relations, the order O(1) Equation (A11) is identically zero. After integrating over all values of the piston velocity and by requiring that the order O( ) vanishes, we obtain a relation between the temperature of the right piston and that of the left one: where we have used also the thermodynamic relation Equation (A6).

Temperature of the Gas
Let us now compute the average temperature of the gas. We have to distinguish between two different cases: • Gas between two moving walls • Gas between a moving wall and a thermostat a. Gas between Two Moving Walls Let us consider the observable ∆X equal to the difference of the gas energy before and after the collision: ∆X = m(v 2 − v 2 )/N . By putting the observable into the relation (16) and by taking into account the Equation (1), we obtain: By integrating, we have: In the same way, we obtain: Putting together Equations (A15) and (A16), one has: and using a Taylor expansion around the small parameter = m M 1 yields: By integrating over the velocity of the piston, because V is zero, the order O(1) is identically zero. By requiring that in the stationary state the order O( ) vanishes, we obtain a relation between the temperature of the gas and those of the near pistons:

Gas between a Piston and a Thermostat
Consider the gas that is near a thermostat and piston. The variation of the temperature due to the piston is the same as before and is given by the Equations (A15) or (A16), respectively, for a piston on the right or on left, with respect to the considered gas. On the side of the thermostat, after the collision, the particle takes a velocity according to the distribution given by the Equation (2). For instance, if the thermostat is that one on the left at temperature T 0 , we have: By integrating, we obtain: In order to compute the variation of the temperature, we have to put together Equations (A22) and (A15): By solving perturbatively around the small parameter = m M 1, we obtain: In the steady state, from relation (A24) at order O(1/ ), we have T = T 0 . This means that a gas near a thermostat reaches the temperature of the thermal bath: this is the result one can obtain from thermodynamics. We can obtain exactly the same result for a thermostat on the right with respect to the considered gas.