Nonequilibrium Steady States for Certain
Hamiltonian Models

In this subsection we review the scheme proposed in [16] for deducing approximate mean disk energies, temperatures, and other profiles in chains that are out of equilibrium. This scheme is based on the following ideas:

The rationale for the first assumption is that when the length of the opening $|\gamma|$ is small, a particle will remain inside each cell it visits for a long time, eventually losing memory of past events. The second assumption is based on the idea of local thermal equilibrium (LTE) [22, 25], which, if true for these models, will imply the convergence of $\mu_{N,[xN]}$ to $\mu^{T(x),\rho(x)}$ as $N\to\infty$. The idea in [16] is to treat the approximations in Assumptions 1 and 2 as though they were exact, and to derive under these simplifying assumptions mean energy and particle density profiles.

In [16], simulations were performed for models in which $\Gamma_{0}$ is given by the area enclosed by four overlapping disks with small openings introduced, and results were reported for 1D-chains with various parameters. It was found that even for relatively small $N$, e.g., $N=30$, the predictions above compare surprisingly well against empirical data. This suggests that for the models tested, (i) memory effects leading to directional bias in particle movements are indeed negligible, and (ii) local equilibrium is, to a good approximation, achieved at relatively small $N$, i.e., finite size effects are insignificant.

Our first example is a 2D rectangular array on a rectangular-shaped domain with Dirichlet boundary conditions (to be explained). Fixed disks of radius $R\in(0,1)$ are placed at even lattice points $\{0,2,\cdots,2M\}\times\{0,2,\cdots,2N\}~.$ We think of $R$ as being close to $1$, so that each square configuration of 4 fixed disks bounds an open cell with 4 narrow gaps, with gap size $2(1-R)$. In the center of each cell, i.e., at odd lattice points $\in\{1,\cdots,2M-1\}\times\{1,\cdots,2N-1\}$, we place a small rotating disk of radius $r$. Heat baths are viewed as occupying the strip of cells just outside the system, injecting and removing particles through the associated openings. For simplicity, we consider only configurations where all baths along the same edge have the same parameters $(T,\rho)$. For example, if the bath along the left edge has parameters $(T_{L},\rho_{L})$, then particles with energies distributed according to Eq. (1) with $T=T_{L}$ are injected at rate $\rho_{L}$ into each cell in the leftmost column, through its opening on the left. We refer to the boundary condition described above as a “Dirichlet” boundary condition.

Returning to the finite system described in the first paragraph, let $J_{k,\ell}$ and $Q_{k,\ell}$ denote the particle and energy exit rate from the $(k,\ell)$th cell. The particle balance equation is the discrete Laplace equation

with boundary conditions $J=4\rho_{\rm L}$ on the left edge, $J=4\rho_{\rm R}$ on the right edge, etc. The $Q_{k,\ell}$ are described by the discrete Laplace equation as well, with boundary conditions $Q=4\rho_{\rm L}\cdot\frac{3}{2}T_{\rm L}$ on the left, $Q=4\rho_{\rm R}\cdot\frac{3}{2}T_{\rm R}$ on the right, and so on. In the continuum limit, these equations tend to the Laplace equation with the given boundary values. Solutions for $J$ and $Q$ are then used to compute the mean disk energy and mean particle count through

These formulas are virtually identical to those for 1D-chains given in Sect. 1.2, differing only in some constants, and can be deduced the same way. The differences in constants are due to the fact that each cell here is in contact with 4 (rather than 2) cells/baths.

Fig. 3 shows profiles for the jump rate $J$, the energy exit rate $Q$, the mean disk energy $D$, and the particle count $\kappa$. Boundary values have been filled in for ease of visualization. We have found that the discrepancy between simulation results and predicted profiles are greatest in the “L-shaped” regions at corners along the boundary of the domain. Away from these regions, the discrepancy is $\leq 1.5\%$ . That the discrepancy is larger here is clearly due to discontinuous boundary values and finite-size effects.

Before leaving this example, let us check our intuition against the plotted profiles. Observe first that panels (a) and (c), i.e., plots of jump rates $J$ and disk energies $D$, reflect the boundary conditions in a straightforward way. The saddle-shaped profile of $Q$ (panel (b)) is also easy to explain: it is because along the 4 edges clockwise from left, $Q=4\rho\times\frac{3}{2}T$ has ratios 4 : 3 : 5 : 1 . A factor that influences particle count is that all else being equal, particles tend to gather in low-temperature regions. The reason is that in a steady state, the total particle flux into and out of any region must be equal. At the lower temperature end, particles move slowly, hence more of them are needed to produce the same flux. For the boundary conditions considered here, the accumulation of particles on the left due to colder temperatures is exaggerated by the higher injection rate on the same side. This has resulted in the huge spread for $\kappa$ in panel (d), from $>120$ on the left to $<5$ on the right. If we change the boundary conditions to $\rho_{L}<\rho_{R}$, the profile for $\kappa$ would be quite a bit flatter. This is both predicted and confirmed in simulations.

The main novel feature of this example is the use of thermostats in the present context.¹¹1See, e.g., [8, 17, 18, 20, 24, 33, 39] for more examples and thorough discussions of thermostats in nonequilibrium statistical mechanics. A second feature of note is the use of exit-only boundaries, which for convenience we will refer to as open boundaries.

We consider in this example a cooling system on a finite cylindrical domain with an open boundary at one end. More precisely, the domain is defined by a rectangular array of fixed disks with a periodic boundary condition in the vertical direction. The right end of the cylinder is put in contact with heat baths, with which it has the same type of interaction as in the previous example. The left end of the cylinder has an open boundary, meaning particles can leave but no new particles enter through it. This configuration leads to a particle flux across the system, from right to left. As the relatively hot particles move through the system, they are cooled (either directly or indirectly) by a row of thermostated cells. For simplicity, we assume the thermostated cells occupy exactly one row which extends the full length of the cylinder, and all thermostats in these cells are set to a fixed temperature $T_{*}$ considerably lower than that of the heat bath.

First, one solves for $J$ as before, the only difference being that the boundary conditions for $J$ are periodic in the $y$ direction. The open boundary for $x=0$ corresponds to setting $J(0,y)\equiv 0$ for all $y$. Solving the boundary value problem readily yields $J$.

To find $Q$, the left and right boundaries are treated in a similar fashion as for $J$. To handle thermostats, we assume that the array is mapped to a domain $[0,a]\times[0,b]$ with $\{y=0\}$ and $\{y=b\}$ identified, and that the thermostats are at $y=0$. We then set, for each $(x,y)$ with $y\in\{0,b\}$, $Q(x,y)=\frac{3}{2}J(x,y)T_{*}$ , where $J(x,y)$ is as in the last paragraph. This is again a boundary value problem, which can be solved to give $Q$.

Once $J$ and $Q$ have been computed, profiles for $D$ and $\kappa$ are given by the same formulas as in the previous example.

We point out two reasons why one may expect greater discrepancies between predicted and simulated results in the present model than in the Dirichlet case: (1) In setting the “boundary values” for $Q$ along the row of thermostated cells, we have assumed that the outgoing particles have energy $\frac{3}{2}T_{*}$. In reality, a highly-energetic particle that enters a thermostated cell is likely to leave it with somewhat more energy than $\frac{3}{2}T_{*}$ : it interacts with the rotating disk only a finite number of times, keeping the normal component of its velocity each time. (2) At the open boundary, $J$ must necessarily go to zero, as must therefore particle count. Since a healthy number of particles is needed in each cell to ensure proper mixing, the situation near this boundary can be anomalous.

Also, with few particles, a practical concern is inadequate updating during the course of the simulation.

The effects highlighted in panel (d), however, seem confined to the thermostated cells. Away from them, the discrepancy between predicted and simulation profiles are about 1-4% for the mean disk energy and $\leq 7\%$ for the particle count. The latter shows larger discrepancies near the open (left) end of the domain, due presumably to $\kappa\approx 0$ near there.

Our third example is a hexagonal array with point sources and open boundaries. The array we use consists of fixed disks of radius $R$ placed on the vertices of a triangular mesh; each edge of the mesh has length 1. Each triangular array of 3 fixed disks bounds a cell. At the center of each cell, we place a small rotating disk of radius $r$. See Fig. 5. The system has open boundaries with no bath injections. Instead, it is driven by a finite number of point sources which pump particles directly into the interior of the array. We say a cell contains a point source with parameters $(T_{*},\rho_{*})$ if instead of being an ordinary cell containing a rotating disk, it emits particles having the energy distribution in (1) with $T=T_{*}$ at a rate of $\rho_{*}$  irrespective of what enters the cell. These particles are injected into each of the 3 neighboring cells with equal probability, i.e. at a rate of $\rho_{*}/3$ each. Notice that by this definition, point sources are “sinks” as well, in that they also absorb particles.

Nevertheless, one can still apply the prediction scheme for finite-size systems: we simply solve the discrete $J$ and $Q$ balance equations without taking continuum limits.²²2Another possibility is to obtain a nontrivial continuum limit by allowing the injection rate $\rho$ to scale with domain size. We do not consider this possibility further. These equations take the form

where $\sim$ is the neighbor relation on the hexagonal grid and the equation holds for all cells $c$ that do not contain a point source. As before, open boundaries are implemented by imagining a strip of cells $c$ just outside the array and setting $J_{c}=Q_{c}=0$ for these cells. Finally, a cell $c$ containing a point source is implemented by defining the values of $J_{c}$ and $Q_{c}$ for that cell by

where $(T_{*},\rho_{*})$ are the parameters for the point source in cell $c$ .

With the exception of $\kappa$, which depends explicitly on the geometry of cells, all profiles can be derived from $J$ and $Q$ in exactly the same way as before. The expression for $\kappa$ here is easily found to be

These results are as expected: the disk energy and particle count profiles clearly reflect parameters of the point sources. One also sees clearly that the number of particles tends to zero as one approaches the boundary of the domain, as it should for open boundaries. The sparsity of particles explains the relative noisiness near the boundaries of Fig. 6(a) (as discussed in relation to the left edge in Example 2).

We have extended the recipe for profile prediction proposed in [16] to a variety of situations and tested it in many concrete examples. In each of the three examples presented (and the many others not shown), we have found that (i) the predicted profiles are very simple to compute, and (ii) the discrepancies between predicted and simulated values are no more than a few percentage points. Given that our parameter choices take the systems tested very far out of equilibrium and the system sizes considered are only moderate, we find the agreement between predicted and simulated results strong. We conclude that the scheme proposed is effective for the class of models studied.

The first part of this subsection has some overlap with [15], which contains a phenomenological study of the same bias questions. Our findings are consistent with those in that earlier study, though we take a slightly different approach, and the bias-correction equations we derive (which are simpler than those in [15]) will be tested in models with significant bias.

The 1D-chain used here is exactly one row of the 2D-rectangular configuration used in the previous section, with the top and bottom gaps of each cell sealed, i.e., particles reflect off these two straight segments elastically. In Parts (A) and (B) below, we consider systems with fixed disk radius $R=0.95$ (leading to gap size $=0.1$) and rotating disk radius $r=0.2$, the same parameters as those used in Examples 1 and 2 of the previous section. These parameters will be varied, with consequences, in Part (C).

The analysis below applies only to situations when there is a sufficient number of particles in the cell. We assume a mean occupation number of at least about 5.

When every particle spends a very long time in a cell before moving on, one may expect the particles and the rotating disk in the cell to come close to equilibrating. Thus one may think that small gaps (or passageways) between cells speed up the convergence to local equilibrium. In [16] and [15], these gaps were taken to be very small: on average, a particle hits the rotating disk more than $12$ times on each visit to a cell. In Examples 1–3 in the present paper, these conditions are relaxed to 3-5 hits per visit. But are small gaps a necessary condition? Put differently, suppose the gaps are large (as in Fig. 7(c)), so that each particle moves through every cell quickly, hitting the rotating disk only once or twice (or even less) on each pass. Would such incomplete equilibration necessarily lead to anomalous behavior?

We now attempt to address this question. The short answer is: A particle’s failure to completely equilibrate with its cell environment on each pass is, in and of itself, not an obstruction to normal behavior. However, large gaps can amplify finite-size effects.

For lack of a better term, let us call a cell boundary “porous” if a particle passing through the cell is not likely to collide with the rotating disk too many times before exiting. For the type of models discussed earlier, porosity is equivalent to larger gap size. We propose below a connection between porosity and finite-size effects, and provide numerical confirmation of our thinking.

Building on previous ideas, our results in Sects. 2,  3.1, and 3.2 point to a very rich class of out-of-equilibrium Hamiltonian models for which we have a simple and effective scheme for predicting approximate energy and particle density profiles. We recapitulate here some of the ideas, and discuss further implications.

Retaining the particle-disk interaction introduced in [35] and [27], we have shown that the ideas first put forward in [16] are applicable to a very flexible collection of models:

First, the results are not limited to 1D; simulations have been carried out in 1D and 2D, and we have seen no reason why the type of prediction scheme studied in this paper cannot be applied in three and higher dimensions, provided the models have suitable interactions between moving particles and pinned-down objects. The type of lattices over which to place the rotating disks appears not to be an issue either; we have illustrated it for rectangular and hexagonal lattices. We have tested a number of different ways to maintain the system out of equilibrium: be it boundary driven (we have tested a variety of boundary conditions), or driven from the interior, using point sources or rows of thermostats – in all cases our prediction scheme has produced easy-to-compute profiles that show strong agreement with simulations.

To understand why memory and finite-size effects are so mild in the models tested in [16] and in Sect. 2, we have proposed an explanation in terms of geometric designs: In these models, the cells have mostly concave (scattering) walls, are connected by narrow passageways, and have rotating disks at good distances away from the walls. We have shown in Sect. 3.1, as have the authors of [15] earlier, that such designs give rise to very low degrees of directional bias. We have also explained in Sect. 3.2 that the small passageways, which facilitate equilibration within cells, also ensure the efficient absorption of bath and thermostat information, reducing finite-size effects. This explains to some degree the rather surprising fact that in [16], for example, the prediction scheme gives remarkably good results in chains with a mere 30 cells.

Going beyond the models discussed in the last paragraph, we have found bias to be small even as we relax considerably the parameters there, i.e., Assumption 1 in Sect. 1.1 applies to a quite large collection of models. However, there are equally natural geometries that give rise a non-negligible amount of bias. For these geometries, we have derived a correction scheme which is relatively straightforward to implement, and shown it to give accurate results (see Sect. 3.1). In a different direction, we have produced numerical evidence in support of our ideas on porosity of cell boundaries and finite-size effects (Sect. 3.2). If these heuristic ideas could be further developed and validated, the implications would be that the results in [16] and Sect. 2 of this paper will remain valid for systems with large gap sizes (perhaps no walls at all) – provided we go to larger system sizes.

We finish with two examples: Example 6(A), along with Example 5, represents one of the greatest departures from the models studied in [16] for which our prediction scheme remains effective; Example 6(B), a close cousin of Example 6(A), makes a connection to the models studied in [27].

We consider a 1D $N$-chain as before, with cells $i=1$ and $i=N$ being in contact with heat baths. Particles are replaced immediately upon exit. To be definite, we assume that when a particle exits the system at location $x\in\gamma$, where $\gamma$ is either the left opening of cell $1$ or the right opening of cell $N$, independently of its exit velocity it is replaced by a new particle which appears in the same position $x$ and has a new velocity $v^{\prime}$ sampled from the bath distribution in Eq. (1). Equivalently, one can think of the gaps leading outside the chain as being sealed off by two “thermostated walls,” and when a particle hits one of these walls, it is reflected but with a new velocity. With these particle-conserving dynamics, the total number of particles $M$ remains constant in time. Taking a continuum limit means letting $M,N\to\infty$ while keeping the density $M/N$ constant.

What is different here is that because particles are replaced upon exit, the boundary conditions for $J$ are $J_{0}=J_{1}$ and $J_{N}=J_{N+1}$ . These (discrete) homogeneous Neumann boundary conditions lead, assuming zero bias, to $J(x)\equiv J_{0}$, i.e., the $J$ profile is constant in space. The $Q$ profile has the same boundary conditions as before and is linear. This leads to a linear temperature profile. The solutions to the $J$ and $Q$ equations are unique only up to multiplicative constants: they are scaled to match the total particle number in the chain.

The discrepancy in the $J$ profile is, at first, a little puzzling: it is somewhat larger than in similar non-particle-conserving chains, and is not improved by doubling the total number of particles to $4000$. At the same time, estimates along the lines discussed in Sect. 3.1 show that the amount of bias here ($\Delta\rho_{out}/J\leq 0.0025$) is relatively small. A possible explanation is that in the present setting, bias can and does lead to a net gradient in $J$ across the length of the chain, whereas in non-particle-conserving systems the values of $J$ at the two ends are fixed by the baths.

Fig. 11(b) shows bias-corrected predictions for $J$ using the same numerical values of the bias coefficients $\alpha_{\rho}$, $\beta_{\rho}$, etc., as in Sect. 3. Agreement with simulation data is good, and lends further support to the idea that bias is a “local” phenomenon, governed mainly by the mean values of thermodynamic quantities and their gradients.