DSMC: A Statistical Mechanics Perspective

Three scientific eras dawned in the middle of the last century: the atomic age, the space age, and the age of the computer. Numerical algorithms were developed for the study of fluids under extreme conditions, such as nuclear explosions and orbital reentry, that were beyond the realm of conventional laboratory experiments. An early example is the Metropolis Monte Carlo method, which in 1953 showed that a particle representation could be used to map phase space and using statistical mechanics extract a fluid’s thermodynamic properties, such as the equation of state.[39] It was closely followed by molecular dynamics (MD), introduced in 1959 by Alder and Wainwright, which uses a similar particle representation but includes dynamics, allowing it to measure transport properties such as viscosity.[2]

In 1963, Graeme Bird submitted to Physics of Fluids a paper entitled “Approach to Translational Equilibrium in a Rigid Sphere Gas.” In this article he introduced a stochastic collision algorithm that would become the cornerstone of Direct Simulation Monte Carlo (DSMC). The paper was rejected (G.A. Bird, personal communication). The referee had failed to recognize that DSMC differed significantly from MD; after revising it to explain these differences, Bird’s paper was accepted.[12]

DSMC went on to become the dominant numerical method in aerospace engineering for rarefied gas dynamics while MD established popularity in the field of statistical mechanics for the simulation of microscopic systems. For many years, these communities remained separate, yet recently they have found common interests in topics such as nanofluidics and granular flows. This paper presents several short stories in which DSMC is used by researchers, such as myself, whose background is in statistical mechanics and molecular dynamics. As we shall see, by approaching DSMC from an unconventional perspective, a number of new avenues are opened. This narrative is intentionally personal, since Graeme Bird, to whom this paper is dedicated, was a colleague and friend for nearly 30 years. In the interest of space, a summary of the DSMC algorithm is omitted; see [4] for an introduction and [14, 15, 20] for full treatments.

In 1989 I joined the faculty at San Jose State University, located in Silicon Valley. About that time I met Berni Alder and started working with him at the nearby Livermore Lab. Alder took an immediate interest in DSMC and we eventually published 16 papers together, most of them on ways to extend the DSMC algorithm. Yet when I first explained DSMC to Berni, he was puzzled and asked “But why does it give the ideal gas law?” Graeme Bird did not address this question in his books on DSMC [13, 14, 15] since he always viewed the algorithm as numerically solving the Boltzmann equation, a perspective that was reinforced by a proof by Wagner [48].

In their seminal MD simulations Alder and Wainwright [3] measured the equation of state for hard spheres (mass $m$, diameter $d$) as

where $p$, $n$, $T$ are pressure, number density, and temperature; $k$ is the Boltzmann constant. The collision rate is $\Gamma$ and

is the virial. For a collision between particles $\alpha$ and $\beta$, $\Delta\vec{v}_{\alpha}$ is the change in velocity for particle $\alpha$ and $\vec{r}_{\alpha\beta}=\vec{r}_{\alpha}-\vec{r}_{\beta}$ is their separation. The virial is an average over collisions and depends on the particle diameter since $|\vec{r}_{\alpha\beta}|=d$. In (1) the pressure is the sum of translational and collisional contributions to momentum transfer.

In DSMC collisions occur within cells and the probability of a collision is independent of the positions of the particles, so $\langle\vec{r}_{\alpha\beta}\rangle=0$ and $\langle\Delta\vec{v}_{\alpha}\cdot\vec{r}_{\alpha\beta}\rangle=\langle\Delta\vec{v}_{\alpha}\rangle\cdot\langle\vec{r}_{\alpha\beta}\rangle=0$ so $\Theta=0$. Since the virial contribution is zero, the equation of state for DSMC is the ideal gas law. Note that Wagner’s proof is for infinitesimal collision cells yet we have the ideal gas law even for finite cell size due to symmetry in the selection of collision partners in DSMC.

Alder found this unsatisfactory, so, together with Frank Alexander, we formulated the Consistent Boltzmann Algorithm (CBA).[5] It is identical to DSMC except that collisions change both the velocities and the positions of colliding particles. Specifically, after computing the post-collision velocities the positions are shifted as

where $\vec{\hat{a}}_{\alpha}=\Delta\vec{v}_{\alpha}/|\Delta\vec{v}_{\alpha}|$ is a unit vector in the direction of the apse line (line between centers) for a hard sphere collision with this change in velocity. Simulations verified that CBA gives the correct pressure and reasonably accurate transport coefficients up to moderately high densities.

The Consistent Universal Boltzmann Algorithm (CUBA) makes the post-collision displacement a function of density and temperature, which allows one to choose the desired equation of state.[6, 34] For example, by choosing the van der Waals equation, we were able to simulate the condensation of vapor into liquid droplets. However, for simulating liquids, we found that CUBA was not competitive with MD in terms of computational cost.

Today, the most popular variants of DSMC for dense gases are based on the Enskog equation.[27, 40] Instead of using a CBA displacement, the correct virial $\Theta$ can be obtained by using a collision probability that depends on both the relative velocity and the relative position for a collision pair.[25] In other words, the Enskog perspective is that particle positions affect collisions, while in CBA, collisions affect particle positions.

In his first book, Graeme Bird said of DSMC, “the results are remarkably insensitive…and no deleterious effects are generally present…if the cell size approaches the mean free path or if the time step approaches the mean collision time.”[13] Nevertheless, the demonstration program in this book for the Rayleigh problem sets the collision cell size to half a mean free path and the time step to a fifth of a mean free collision time.

A common way of quantifying the error due to cell size and time step is to measure viscosity or thermal conductivity and, at low Knudsen number, compare with Chapman-Enskog theory. These transport coefficients can be obtained by measuring the momentum (or heat) flux in a DSMC simulation with a linear velocity (or temperature) gradient.[19] From the perspective of non-equilibrium thermodynamics, the transport coefficient is an Onsager coefficient that relates the thermodynamic flux to the thermodynamic force due to the gradient.[24, 41]

Yet in the early days of molecular dynamics, such non-equilibrium measurements were not feasible since simulations were limited to a few hundred particles. Instead, transport properties were measured in simulations of equilibrium systems using results from statistical mechanics.

The basic idea is to apply the fluctuation-dissipation theorem, which connects statistical fluctuations to the rate of dissipation. The best known example is Brownian motion in which the mobility, $\mu$, of a particle is obtained from measuring the variance of its position in time as

By contrast, a non-equilibrium measurement obtains the mobility $\mu=v_{\mathrm{T}}/F_{\mathrm{ext}}$ from the terminal velocity, $v_{\mathrm{T}}$, for the particle under an external force of magnitude $F_{\mathrm{ext}}$.

Using Einstein-Helfand theory one can measure the viscosity of a fluid from molecular trajectories as [1]

where $N$ is the number of molecules, $V$ is the system volume, $x_{i}$, $y_{i}$ are components of $\vec{r}_{i}$, and $\dot{x}\equiv dx/dt$. This allows us to measure the viscosity (and other transport coefficients) from equilibrium simulations. A related approach uses the off-diagonal stress tensor component,

where $\vec{F}_{ij}$ is the force between molecules. Our previous result now takes the more familiar form from Green-Kubo theory,

Notice from (6) that $J$ is the sum of a kinetic (ballistic) contribution and a contribution due to the transfer of momentum by intermolecular forces.

For hard sphere collisions Wainwright wrote [49]

where $t_{c}$ is the time of a collision. He was able to evaluate (7) and express it as the sum of three contributions, $\eta=\eta^{\mathrm{K}}+\eta^{\mathrm{K}\times\mathrm{C}}+\eta^{\mathrm{C}}.$ The kinetic contribution $\eta^{\mathrm{K}}$ is precisely the Chapman-Enskog viscosity. The cross term $\eta^{\mathrm{K}\times\mathrm{C}}$ and the collision term $\eta^{\mathrm{C}}$ are dense gas corrections.

In DSMC, the collision term $\eta^{\mathrm{C}}$ is not zero due to the finite distance between the collision partners (i.e., the finite collision cell size). Unlike the equation of state contribution due to the virial (see (2)), this contribution to the viscosity does not average to zero over DSMC collisions within a cell. In fact, following a similar analysis as Wainwright one finds

For cubic collision cells of size $\ell$ the average $\left\langle(y_{\alpha}-y_{\beta})^{2}\right\rangle=\ell^{2}/6$; at equilibrium $\left\langle(\Delta\dot{x}_{\alpha}-\Delta\dot{x}_{\beta})^{2}\right\rangle=4kT/3m$. In DSMC the cross term $\eta^{\mathrm{K}\times\mathrm{C}}$ is zero and so the hard sphere viscosity is [7]

where $\lambda$ is the mean free path. Note that the error goes as $\ell^{2}$; for a cell size of one mean free path the viscosity deviates from its Chapmann-Enskog value by 11 percent. A similar analysis (with a similar result) may be applied to obtain the error in thermal conductivity due to finite cell size. Finally, this approach also gives the error in the transport coefficients due to a finite time step $\tau$; in DSMC this error goes as $\tau^{2}$.[32, 33]

In his last book Bird wrote “(In) early DSMC programs… the mean spacing between collision pairs was too large in comparison with the mean free path and this led to the introduction of sub-cells…(and) an option to select the collision pairs from the nearest-neighbour pairs within the cell.”[15] He also mentions that “The time step is set to a specified small fraction of the sampled mean collision time.” A challenge for modern DSMC codes is to minimize the distance between colliding particles and to minimize the time between ballistic movement steps while avoiding unphysical bias errors (see Section 5).

My introduction to DSMC came in 1983 as I was finishing my doctoral studies at The University of Texas at Austin in what is now The Ilya Prigogine Center for Studies in Statistical Mechanics and Complex Systems. At that time there were theoretical results, confirmed by light scattering experiments, indicating that hydrodynamic fluctuations had weak, long-ranged correlations in non-equilibrium systems (e.g., under a temperature gradient).

In his first book, Bird describes fluctuations in DSMC as an unavoidable annoyance but that with a sufficiently large sample size the desired mean value can be measured to the desired accuracy.[13] One commonly held misconception was that these fluctuations came from the Monte Carlo nature of the collision algorithm in DSMC. Yet the very same fluctuations are also found in molecular dynamics, which has deterministic dynamics.

Due to computer memory constraints, the number of particles (or “simulators”) in a DSMC calculation is typically a small fraction of the number of physical molecules. For example, at standard temperature and pressure, the mean free path in air is roughly 60 nanometers and a cubic mean free path contains about 6000 molecules. In contrast, a typical DSMC simulation uses 20 particles per cubic mean free path, that is, each DSMC particle represents $F_{N}=300$ physical molecules. This number scales as the Knudsen number so it is much higher for rarefied flows.

Statistical mechanics gives us the standard deviation of hydrodynamic variables, allowing us to estimate the expected statistical standard error in particle simulations.[35] For example, in a cell with $N_{\mathrm{c}}$ particles, the variance in each component of fluid velocity is $\langle\delta u_{x}^{2}\rangle=kT/mN_{\mathrm{c}}$.[36] For both MD and DSMC the fractional error in estimating $u_{x}$ from $S$ independent statistical samples is

where $\mathrm{Ma}$ is the Mach number of the flow. This result illustrates why nanofluidic flows, in which $\mathrm{Ma}\ll 1$, are much harder to measure in DSMC compared to hypersonic flows. In DSMC $N_{\mathrm{c}}$ is smaller than the number of physical molecules by a factor of $F_{N}$, the variance in fluid velocity increases by the same factor.

After accounting for this amplification factor, we find that hydrodynamic fluctuations in DSMC are physically correct [31] and yield hydrodynamic information. For example, from the time correlation of density fluctuations one may measure the sound speed and viscosity from the location and width of the Brillouin peak, mimicking the technique used by light scattering.[10, 21] This dynamic structure factor is a useful way to compare DSMC models (e.g. internal energy relaxation) with laboratory measurements and theoretical predictions.[22]

The first numerical observations of long-ranged correlations of hydrodynamic fluctuations were made using DSMC for systems with a temperature gradient [37] and a velocity gradient [30]. The “giant fluctuation” phenomenon, first reported in diffusive mixing experiments [46], is due to such non-equilibrium correlations for fluctuations of concentration and fluid velocity, as confirmed by DSMC [26]. Most recently DSMC simulations have been used to study the impact of thermal fluctuations on turbulence.[28, 38] They confirm theoretical predictions [9] and numerical fluctuating hydrodynamic observations [11] that molecular fluctuations dominate the turbulent energy cascade at the Kolmogorov length scale.

Bird wrote in his last book, “While the fluctuations are unphysical when $F_{N}$ is large, they are physically realistic… (with) a one-to-one correspondence between real and simulated molecules. This is another instance of DSMC going beyond the Boltzmann equation because fluctuations are neglected in the Boltzmann model.” His last three papers were on Brownian motion.[16, 17, 18]

In 1887 Maxwell presented a thought experiment, “… conceive of a being whose faculties are so sharpened that he can follow every molecule in its course … so as to allow only the swifter molecules to pass from (chamber) A to B, and only the slower molecules to pass from B to A. He will thus, without expenditure of work, raise the temperature of B and lower that of A, in contradiction to the second law of thermodynamics.” The DSMC algorithm has full information and control of the particle’s dynamics and modern implementations have many complex stages. How can we be sure that we do not have a Maxwell demon hiding in our simulations? In other words, how do we test that our simulations satisfy the Second Law of Thermodynamics? Boltzmann’s H-theorem is not helpful, since, due to fluctuations, it only applies in the limit $N\rightarrow\infty$.

Consider an isolated system with microstate dynamics[44]

where $p_{i}$ is the probability of state $i$ and $W_{ij}$ is the transition rate from $i$ to $j$. At equilibrium $\dot{p}_{i}=0$. If $W_{ij}=W_{ji}$, then we have microscopic reversibility. At equilibrium we have detailed balance ($p_{j}W_{ji}=p_{i}W_{ij}$) and ergodicity ($p_{i}=p_{j}$) iff we have microscopic reversibility. Finally, the dynamics satisfies the Second Law ($\dot{S}=-k\sum_{j}\dot{p}_{j}\ln p_{j}\geq 0$) iff we have ergodicity.[44]

Detailed balance is a necessary condition for microscopic reversibility, which is a sufficient condition for the dynamics to obey the Second Law. Yet, it is not uncommon to find DSMC implementations that do not satisfy detailed balance. For example, some of the early surface scattering models were ad hoc empirical fits to experimental data, which led Cercignani and Lampis to formulate a model that imposed detailed balance.[23] Many implementations of inflow/outflow boundaries also fail to satisfy detailed balance.[45] If machine learning methods are used to handle some elements in a DSMC calculation, such as chemical reactions and internal degrees of freedom, we may not know if detailed balance is satisfied.[8] And even when the algorithm in a DSMC code is theoretically sound, there is always the possibility that the implementation has a ’bug’ (coding error).

In testing for the presence of Maxwell demons in a DSMC code, equilibrium hydrodynamic fluctuations can serve as a coal mine canary. For example, at equilibrium, the number of particles in each cell is independently Poisson distributed, and this result is independent of the collision model, boundary conditions, etc. Statistical mechanics provides further relations for the variances and correlations of fluid velocity, temperature, internal energy, and other variables.[29, 36] While this is a necessarily but not sufficient statistical test, it is a sensitive one that can reveal malevolent problems. Ensuring the correctness of simulation data is increasingly important as we enter the age of artificial intelligence. Given the AI systems’ voracious appetite for such data, we must avoid having machines learn the wrong lessons.

Graeme Bird and Berni Alder met for the first time at the 22nd Rarefied Gas Dynamics Symposium in Sydney, Australia, nearly 40 years after their original publications on DSMC and MD. In the 25 years since that meeting, we find increasing cooperation between the two algorithms’ communities. One example is the study of rotational relaxation using all-atom molecular dynamics simulations combined with DSMC.[47] Yet perhaps the best example is the shared heritage of LAMPPS [43] and SPARTA [42], two of the most widely used codes for molecular dynamics and direct simulation Monte Carlo. Although separated at birth, MD and DSMC are now the patriarchs overseeing a large family of particle-based fluid algorithms.

The author acknowledges support from the US Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under contract no. DE-AC02-05CH11231.