Appendix A - The Direct Simulation Monte Carlo

The Direct Simulation Monte Carlo, also known as Bird algorithm, is a simulation scheme useful to solve the Boltzmann equation for gases. The Boltzmann equation is characterized by a linear partial differential operator plus the so-called collision integral (see paragraph 2.2.3):

$$Q(P,P)=\int_{\Re^3}\int_{S_+} (P'P_*'-PP_*)\vert\mathbf{V} \cdot \hat{\mathbf{n}}\vert d\mathbf{v}_*d\hat{\mathbf{n}}$$ (A.1)

It is a five-dimensional integral that makes difficult numerical approximations: if we want to approximate this integral by a quadrature formula, we need the evaluation of the integrand in $N^5$ points where $N$ is the number of points necessary to well approximate a one-dimensional integral. For example, if $N=20$, than we must evaluate the integrand in $3.2 \cdot 10^5$ points, at every time step of our integration algorithm. Particle simulation methods (like the Bird algorithm) avoid this problem by replacing the density distribution $P$ by a discrete measure of test particles, the so-called ``simulation gas''. The number of test particles can be anything between $100$ and several millions, depending on the capacity of the computer being used. The key idea of particle simulations stems from the consideration that the physical reality underlying the problem is a gas with, say, $10^{25}$ particles and that the Boltzmann equation is only a mathematical approximation to this reality: the use of this equation is a drastic reduction of the number of freedom degrees (even if it replaces a finite number of degrees with a six-dimensional infinite number of them, we can consider the Boltzmann equation the result of a compression of information), and therefore we can replace it with some other kind of approximation that reduces the number of freedom degrees. Therefore, in particle simulations we return to the particle level, but restrict the number of particles in order to make it tractable, and consequently change the interaction rules in order to make them reflect the influence of the collisions on the behavior of the gas.

The Bird scheme [29,30] has been successfully employed for the simulation of rarefied gases for decades (for example in the calculations for the reentry of spacecrafts). Recently non constructive [218] and constructive [185] proofs of its convergence to the solution of the Boltzmann equation have been given.

The scheme consists of the following ingredients:

• Time discretization: time is discretized, i.e. $t=j \cdot \Delta t$ and the free flow operator $S_t^0$ (see paragraph 2.2.1) is approximated by $S_{\Delta t}^0 \sim 1+L^0 \Delta t$ (Euler-type approximation).

• Splitting: the complete evolution of the gas in the time integration step $\Delta t$ is divided in two sub-steps:
  
  

  This means that the particles evolve independently for a time $\Delta t$ and then the p.d.f. $P$ is modified in order to take into account the collision, by evaluating the integral $Q$ (this is the difficult part).

• Local homogenization (over space cells): the space domain is divided in non-overlapping cells $C_k$; the intermediate step in Eq. (79a) is performed on a p.d.f. $\bar{P}$ which is the spatial average of $P$ on those cells, i.e.:

  $$\bar{P}(\mathbf{r},\mathbf{v},j\Delta t)= \frac{1}{V(C_k)} \int C_k P(y,\mathbf{v},j\Delta t)dy$$ (A.3)

  
  

  and therefore also $\tilde{P}$ is homogeneous on the cells $C_k$ of volume $V(C_k)$. The free flow step in Eq. (79b) will destroy the local homogeneity. We assume that all the cells $C_k$ have the same diameter and call it $r_B$ ($B$ stays for Bird but also for Boltzmann, so that $r_B$ is called Bird radius or Boltzmann radius indifferently). The local homogenization is important for both practical and conceptual reasons: it makes more easy (and rapid) to simulate the collisions among particles and, at the same time, fulfills the Boltzmann idea that the spatial variation of the particle density over a cell with volume $dx dy dz$ could be neglected: it is the velocity dependence of the particles in such a cell that would make a difference.

.

Here follows the algorithm, step by step:

1. Evaluate the new positions at time $(j+1)\Delta t$ using the free flow operator, i.e. $x_i(j+1) = x_i(j)+v_i \Delta t$.

2. In every cell $C_k$ (containing $N_C$ particles and therefore $\binom{N_C}{2}$ collision pairs) do the following:
   1. reset the ``collision clock'' $t_C=0$.

   2. Choose randomly a pair $i_1,i_2$ in the cell with probability $p_{i_1i_2}$.
      
      

      with $k$ a normalization constant.

   3. Update the velocities of the particles $i_1$ and $i_2$ with the collision rule (for example the one in Eq. (2.33)), using a random choice of the collision parameter $\hat{\mathbf{n}}$.

   4. Update the collision clock $t_C= t_C+ \frac{2 V(C_k)}{N_C^2 \sigma_{i_1i_2}}$.

   5. go to step 2b until $t_C > \Delta t$ (this guarantees an appropriate number of collisions, consistent with the collision integral)

3. repeat the step 2 for all the cells.

4. go to step 1 and advance the time $t=t+\Delta t$.

We have used a slight different scheme where the particles in the cells have a fixed a-priori probability $p_c$ of colliding during a time step $\Delta t$: every particle is checked for collision, throwing a random number in the range $[0,1)$; if the number is lesser than $p_c$ then the particle is admitted for collision (it is particle $i_1$) and a collision mate $i_2$ is chosen in the cell with probability $p_{i_1i_2}$. In a homogeneous situation with $r_B$ of the order of the mean free path, this guarantees that the mean free time is about $\Delta t/ p_c$.

To study more dense situations, alternative schemes have been proposed that take into account spatial density correlations in a way consistent with the Enskog correction to the Boltzmann equation [90,165,166].