Ehsan Roohi
email: e.roohi@ferdowsi.um.ac.ir

Stefan Stefanov

On Efficient Collision Schemes in the Direct Simulation Monte Carlo (DSMC): From micro/nano flows to hypersonic flows

Keywords: Rarefied flow, Boltzmann equation, Kac master equations, DSMC, collision models, simplified Bernoulli trials, micro/nano flows, hypersonic flow

The motivation of this review paper is to present a detailed summary of different collision models developed in the framework of the direct simulation Monte Carlo (DSMC) method. The emphasis is put on a newly developed collision model, i.e., the Simplified Bernoulli trial (SBT), which permits efficient low-memory simulation of rarefied gas flows. The paper starts with a brief review of the governing equations of the rarefied gas dynamics including Boltzmann and Kac master equations and reiterates that the linear Kac equation reduces to a non-linear Boltzmann equation under the assumption of molecular chaos. An introduction to the DSMC method is provided, and principles of collision algorithms in the DSMC are discussed. A distinction is made between those collision models that are based on classical kinetic theory (time counter, no time counter (NTC), and nearest neighbor (NN)) and the other class that could be derived mathematically from the Kac master equation (pseudo-Poisson process, ballot box, majorant frequency, Null collision, Bernoulli trial and its variants). To provide a deeper insight, the

Introduction

With the rapid progress in design and utilization of micro/nano-electro-mechanical-system (MEMS/NEMS), the requirement for detailed analysis of fluid and thermal characteristics in these systems increases. Flow in MEMS/NEMS reveals rarefied behavior as the characteristic length of these devices is comparable to the gas mean free path. Therefore, flow treatment in MEMS/NEMS requires the same numerical approach as the flow over high altitude flying objects. The scaling parameter determining the degree of flow rarefaction is the Knudsen number, which is defined as the ratio of the gas mean free path () to the characteristics dimension (L) of the geometry. That is Kn= /L. Different rarefaction regimes are defined according to the Knudsen number, i.e., continuum (Kn<0.001), slip (0.001<Kn<0.1), transition (0.1<Kn<10), and free molecular (Kn>10) [START_REF] Karniadakis | Microflows and nanoflows: fundamentals and simulation[END_REF]. However, it should be reminded that this classification is mainly based on data obtained from experiments and numerical studies of isothermal gaseous flows in long one-dimensional (1D) microchannels. For flows in 2D and 3D complex geometries, the range of above regimes is debatable and should be reconsidered for each studied problem separately [START_REF] Mohammadzadeh | Thermal and second-law analysis of a micro-or nano-cavity using direct-simulation Monte Carlo[END_REF]. The fundamental kinetic equations describing the rarefied gas flow are the Boltzmann nonlinear equation [START_REF] Kennard | Kinetic Theory of Gases: With an Introduction to Statistical Mechanics[END_REF] and its probabilistic alternative the linear Kac stochastic equation [START_REF] Kac | Probability and related topics in physical sciences[END_REF][START_REF] Kac | Probability and Adjacent Questions[END_REF].

Direct simulation Monte Carlo (DSMC) is a statistical approach widely employed for treating rarefied gas flows [START_REF] Bird | Monte Carlo simulation of gas flows[END_REF][START_REF] Bird | Molecular Gas Dynamics[END_REF][START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. The straightforwardness and strong physical logic of its procedures are some of the advantages of the DSMC method for studying rarefied gas flows. The basic DSMC is defined as a numerical particle representation of the Boltzmann equation. However, Kac stochastic equation was also utilized to derive alternative collision models [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF][START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF]. The first published paper on DSMC by Graeme Bird was concerned to the simulation of a 0-D relaxation problem [START_REF] Bird | Approach to translational equilibrium in a rigid sphere gas[END_REF]. In that paper, 30,000 collision events per hour were simulated using the hard sphere model of 500 simulator particles. With today's computers, more than 10 billion collision events can be simulated per hour for this relaxation problem. Progress in computer technology and development of new DSMC models has allowed for very complex physics to be modeled in the DSMC method, i.e., chemical reactions, evaporation, and condensation of substances, ionization and radiation [START_REF] Bird | Aspects of the Structure of Strong Shock Waves[END_REF][START_REF] Bird | The DSMC Method[END_REF]. In DSMC, the time evolution of the particle system within a small time interval (t) is split into two consecutive stepsfree motion of all particles and collisions of particles localized in the neighborhoods of given points in the space, i.e. in cells of a computational grid assuming fixed their coordinates and changing only their velocities in results of binary collisions. The collision algorithm plays the significant role in the DSMC method and calculates the most sophisticated term of the Boltzmann or Kac stochastic equation. In this review we are focused on the more complicated collision step, which defines the basic features of the DSMC method. The DSMC schemes could be categorized into two general classes/groups concerning the treatment of collisions. Considering the Boltzmann definition for the kinetic equation of a rarefied gas, the concept of the first -group of collision schemes is based on the principle of the maximum collision rate per time step. In this group, the "No Time Counter (NTC)" [START_REF] Bird | The Perception of Numerical Methods in Rarefied Gas Dynamics[END_REF] scheme, considers a superior number of maximum collision rate per time step which determines the number of randomly selected particle pairs that should be checked for accepted collision, while the other methods in this group, e.g. "Time Counter (TC)" [START_REF] Bird | Shock Wave Structure in a Rigid Sphere Gas[END_REF], "Null Collision (NC)" [START_REF] Koura | Null-collision technique in the direct-simulation Monte Carlo method[END_REF][START_REF] Nanbu | Direct Simulation Scheme Derived from the Boltzmann-Equation .I. Monocomponent Gases[END_REF], and "Majorant Frequency Scheme (MFS)" [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF] use a time-interval of 𝛿𝑡 𝑖 for each captured collision and a time-step interval of ∆𝑡 for the DSMC procedure, and continue the collision process until ∑ 𝛿𝑡 𝑖 > ∆𝑡 𝑖 . It should be noted that besides the inherent discretization problems of deterministic numerical approaches, DSMC calculations are accompanied by two extra problems: (a) the presence of statistical noise in output results, and (b) the dependence of results on the number of particles per cell and possibility of repeated collisions. It has been shown [START_REF] Baker | Variance reduction for Monte Carlo solutions of the Boltzmann equation[END_REF][START_REF] Homolle | Low-variance deviational simulation Monte Carlo[END_REF][START_REF] Homolle | A low-variance deviational simulation Monte Carlo for the Boltzmann equation[END_REF][START_REF] Dimarco | Fluid solver independent hybrid methods for multiscale kinetic equations[END_REF] that a modified variance reduced Monte Carlo simulation, which takes into account the asymptotic properties of near continuum low-speed regimes, is capable of overcoming the first problem. In response to the second problem, other type collision schemes, proposed by Belotserkovskii and Yanitskiy [START_REF] Belotserkovskii | The statistical particles-in-cells method for solving rarefied gas dynamics problems[END_REF] and Yanitskiy [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF], were constructed on the base of the Kac stochastic equation. Contrary to the former group of collision schemes, the latter group, in accordance with the Kac stochastic model, defines a collision probability function for each particle pair and check all pair combinations for collision occurrence. The Yanitskiy approach [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF][START_REF] Belotserkovskii | The statistical particles-in-cells method for solving rarefied gas dynamics problems[END_REF] then led to the introduction of the Bernoulli-Trials collision scheme (BT) which benefitted from the avoidance of the repeat collisions. Stefanov introduced a simplified variant of the Bernoulli Trials scheme entitled as 'SBT' [START_REF] Stefanov | Particle Monte Carlo algorithms with small number of particles in grid cells[END_REF][START_REF] Stefanov | On DSMC calculations of rarefied gas flows with small number of particles in cells[END_REF]. Unlike the former scheme (BT), which has a quadratic dependency of the computational cost on the particle number in cells, the latter one (SBT) has a linear dependency and a higher computational efficiency. As an evolution of the SBT scheme to a method which can intelligently prefer collisions for closer pairs, Goshayeshi et al. [START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF] introduced an intelligent variant of the SBT scheme entitled as "ISBT" which provides semicognition of distance for the collision scheme. This semi-cognition reduces approximately 25-32% of the overall mean collision separation distance (MCS) in collision cells. Understanding the recent notes of Gallis et al. [START_REF] Gallis | Convergence behavior of a new DSMC algorithm[END_REF][START_REF] Gallis | Effect of collision-partner selection schemes on the accuracy and efficiency of the direct simulation Monte Carlo method[END_REF] in the preference of choosing a near neighbor partner rather than the nearest neighbor one, which consequently leads to the saving of the collision scheme from losing some of its probable collisions during the advection phase of particles, the ISBT scheme also follows the same strategy of near neighbor pair-selection.

Considering the fact that smaller mean collision separation distances would cause more realistic collisions and prevent from the angular momentum reduction, the modifications have been proposed to both groups of collision schemes. Based on the logic of the first group, LeBeau et al. [START_REF] Lebeau | Virtual sub-cells for the direct simulation Monte Carlo method[END_REF] introduced the virtual sub-cell (VSC) method which performs an O(N 2 ) operation to sort all N simulators in a cell to find the nearest-neighbor to any simulator chosen for collision. Alternatively, Bird [START_REF] Bird | Visual DSMC Program for Two-Dimensional and Axially Symmetric Flows, The DS2V Program User's Guide[END_REF] proposed the transient adaptive sub-cell (TAS) scheme that subdivides cells into sub-cells, and collision pairs are selected from the same or neighboring sub-cells. In the second group of collision schemes, however, there were some constraints in the development of modifications to reduce MCS. For instance, because the collision process is conducted in a hierarchical order of indexed particles, it is not possible to directly search for the nearest pair in such a way like the VSC does, or search in the neighboring subcells in a way that NTC-TAS does. Therefore, in the case of TAS usage, it became only possible to have subcells adapted with 4-5 particles. Therefore, the selection of the nearest-neighbor in the second group of collision schemes is bound with the nearer neighbor selection [START_REF] Goshayeshi | DSMC Simulation of Hypersonic Flows Using an Improved SBT-TAS Technique[END_REF].

The aim of the current review is to present and evaluate a recent Kac-based collision model in the direct simulation Monte Carlo (DSMC) method, Simplified Bernoulli trial (SBT) and its variants, i.e., SBT on transient adaptive subcells (SBT-TAS) and intelligent SBT (ISBT) in the treatment of a wide spectrum of rarefied gas flows either at micro/nano scales or hypersonic flow regimes.

Accuracy, memory and CPU requirements, time consumption and rigorousness of this collision family will be discussed in details. The paper continues with a derivation of the mathematical relation between the Boltzmann and Kac master equations and then introduces the collision schemes based on classical kinetic theory (Boltzmann Equation) in the DSMC method followed by the derivation of Kac-based family collision models such as the Bernoulli trial (BT) and Simplified Bernoulli trial family schemes and presentation of validations of the SBT collision family.

Collision models based on the Boltzmann equation

The most detailed level of description of a system consisting of large number (N) of molecules is given by the Newton equation. The evolution equations of this system are given by [START_REF] Khlopkov | Development of Monte Carlo Methods in Hypersonic Aerodynamics[END_REF]:

𝑚 𝑗 𝑑 2 𝑥 ⃗ 𝑗 𝑑𝑡 2 = ∑ 𝑅 𝑖𝑗 𝑁 𝑖≠𝑗 (1)
, where R ij is the force between pair (i-j). However, the solution of such a system is difficult due to requirement of specifying of the initial coordinates and velocity of each molecule in the system, complexity of the force function, as well as solving a 3N differential-coupled equation.

Following the Gibbs formalism, rather than considering a single system, an ensemble of systems in the 6N-dimensional (3 space and 3 velocity coordinate, i.e. phase space, for each particle) is distributed according to the N-particle velocity (probability) distribution function (F N ). This ensemble system could be treated by the Liouville equation, which describes the time evolution of the phase space distribution function [START_REF] Khlopkov | Development of Monte Carlo Methods in Hypersonic Aerodynamics[END_REF]:

𝜕𝐹 𝑁 𝜕𝑡 + ∑ 𝑐 𝑖 𝜕𝐹 𝑁 𝜕𝑥 𝑖 𝑁 𝑖=1 + ∑ ∑ 𝑅 𝑖𝑗 𝑚 𝑁 𝑖=1 𝑁 𝑖≠𝑗 𝜕𝐹 𝑁 𝜕𝑐 𝑖 = 0 (2)
The Liouville equation and all the succeeding kinetic equations following from the BBGKY hierarchy (Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy, or sometimes called (𝑙) , 𝑐). The velocity distribution function quantifies the number of particles, f(t, 𝑥 (𝑙) , 𝑐) d𝑥 (𝑙) dc, which are located in an infinitesimal volume 𝑑𝑥 (𝑙) at position 𝑥 (𝑙) and whose velocities are in the infinitesimal interval 𝑑𝑐 around 𝑐. The Boltzmann equation is

given by [START_REF] Kogan | Rarefied Gas Dynamics[END_REF] (

, where n is number density, t is time, r is space vector, c is velocity space vector, F is external force per unit mass, c r is the relative velocity between a molecule of velocity class c and one with velocity class c 1 , σ T dΩ is the differential cross-section for the collision of a molecule of class c with another one having class c 1 such that their post-collision velocities are c* and c 1 *,

        4 2 * * 1 1 1 0 .. r nf c nf F nf n f f ff c d dc t r c               
 respectively, and functions f, f 1 , f* and f 1 * are the corresponding velocity distribution functions for the molecule and its collision partner before and after the collisions. In fact, f is the probability of finding a particle in the velocity space volume element 𝑑 3 𝑐around 𝑐, provided that the particle is located in the space volume element 𝑑 3 𝑥 around 𝑥. Collision in the Boltzmann equation is an instantaneous, random jump process in particles velocity localized in a point. Boltzmann equation is valid for dilute gasses so that only binary collisions between the molecules are possible. The sophisticated and nonlinear properties of the collision integral and its containing of seven variables create severe difficulties for the numerical analysis of the Boltzmann equation. High dimensionality and probabilistic nature of the kinetic processes as well as complex molecular collision models are among the main prerequisites for the application of the direct simulation Monte Carlo method to handle the Boltzmann equation. Direct simulation Monte Carlo, as developed by Bird [START_REF] Bird | Monte Carlo simulation of gas flows[END_REF], is a numerical simulation of the evolution of a collection of particles.

Historically, the application of Monte Carlo method for rarefied gas flows was initiated by the pioneering works of Kogan and Perepukhov for simulation of free-molecular flows around flying objects [START_REF] Kogan | Rarefied Gas Dynamics[END_REF]. In fact, Bird used principles of Monte Carlo method, such as acceptance-rejection, to perform numerical computations of the collision terms of the Boltzmann equation, as will be described in the following section. It should be noted that the first DSMC application by Bird was aimed for simple homogeneous relaxation problems, in which a non-equilibrium gas velocity distribution is driven toward equilibrium by intermolecular collisions [START_REF] Bird | Approach to translational equilibrium in a rigid sphere gas[END_REF]. A flowchart depicting different collision models employed in the DSMC method is shown in Fig. 2. In the following, details of these collision models are provided with sufficient discussions.

Time counter (TC) collision model

The time counter method is the first collision schemes introduced by Bird [START_REF] Bird | The Perception of Numerical Methods in Rarefied Gas Dynamics[END_REF]. In this method, the total number of collisions in a cell is assumed to be: Where overbar quantity means averaged valued, and 𝑛 = 𝑁/𝑉 𝑐 is number density. For every particle pair chosen at random, the following time increment is added to the time counter of the cell:

𝑁 𝐶𝑜𝑙𝑙 = Δ𝑡𝑁𝑛𝐹 𝑛𝑢𝑚 𝜎 𝑇 𝑐 𝑟 ̅̅̅̅̅̅ 2 (4)

Kac

Δ𝑡 𝑐𝑜𝑙𝑙 = 2 𝑁𝑛𝐹 𝑛𝑢𝑚 𝜎 𝑇 𝑐 𝑟 (5)

Fig. 3 TC collision procedure

The TC model assumes that each collision contributes equally in the temporal advancement of the collision numbers in the cell. The procedure is repeated until

∑ Δ𝑡 𝑐𝑜𝑙𝑙 ≤ Δ𝑡 𝑐𝑒𝑙𝑙 (6)
Flowchart of the TC collision model is shown in Fig. 3. The computational cost of TC algorithm is of the order O(N). It is mentioned that TC model needs around 20 particles per cell to produce correct results. However, the TC model could not produce correct solution under extreme nonequilibrium conditions such as the front of a very strong shock. In fact, the acceptance of an unlikely collision can advance the time by an interval that is much larger than the time step and Time advancm ent

• Advance time step of the collision using Eq. ( 5) pair (i,j)

• Random selection of a collision pair (i,j) from N l particle in cell l.

Velocity update

•change the particle velocities Ending condition

• Repeat the procesure until condition of Eq. ( 6) is satisfied the overall collision rate can be distorted. Lutisan [START_REF] Lutišan | The Treatment of Molecular Collisions in DSMC Methods[END_REF] showed that TC model could not reproduce

Poisson distribution function for probability distribution of number of collisions.

No-Time-Counter (NTC) collision model [8, 14]

According to the classical kinetic theory utilized in deriving the Boltzmann equation, the collision probability (P) between two particles over the time interval dt is equal to the ratio of volume swept by their total collision cross section moving at the relative speed (c r ) to the volume of the cell which contains particles (V c ), i.e., 𝑃 = 𝐹 𝑁 𝜎 𝑇 𝑐 𝑟 𝑑𝑡/𝑉 𝑐 [START_REF] Bird | Molecular Gas Dynamics[END_REF] The average number of simulated particles in a cell is N=nV c /F n , where n is the gas number density and F N is the ratio of number of real molecules to simulated particles. All possible collisions could be checked by choosing all N(N-1)/2 collision pairs in the cell and checking the collision probability P for all of them, i.e., collision is accepted if P>R nf . This collision algorithm is called Bernoulli trials (BT) [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF], and is inefficient because P is a small quantity and number of choices is proportional to the square of number of particles, O (N 2 ). To increase the efficiency of the collision procedure and reduce the computational costs to O(N), Bird suggested to repair the TC scheme by using the idea of maximum collision cross-section approach well known in Monte Carlo methods used for simulation of linear integral equations describing the transport of neutrons. The idea appears almost around the same time in the papers of Koura [START_REF] Koura | Null-collision technique in the direct-simulation Monte Carlo method[END_REF] and Bird [START_REF] Bird | The Perception of Numerical Methods in Rarefied Gas Dynamics[END_REF]. Then, the acceptance-rejection technique and no time counter method or null-collision methods were used for the first time in DSMC. The universal acceptance-rejection technique [START_REF] Von Neumann | Various techniques used in connection with random digits. Monte Carlo methods[END_REF] requires a superior estimation of the collision probability. It can be written as:

𝑃 𝑚𝑎𝑥 = 𝐹 𝑁 (𝜎 𝑇 𝑐 𝑟 ) 𝑚𝑎𝑥 𝑑𝑡/𝑉 𝑐 (8)
The number of pair selections per time step is calculated by multiplying this equation by N(N-1)/2. Therefore, in this procedure, that is called "No Time Counter (NTC) method, the maximum number of particle pairs checked for collision is computed as:

𝑁 𝑐𝑜𝑙𝑙 = 1/2𝑁(𝑁 -1)𝐹 𝑁 (𝜎 𝑇 𝑐 𝑟 ) 𝑚𝑎𝑥 𝑑𝑡/𝑉 𝑐 [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF] , and this number of pairs is selected at each cell per time step, and the collision of pair (𝑖, 𝑗) is accepted with the following probability given by Eq. [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF].

(𝜎 𝑇 𝑐 𝑟 ) 𝑖𝑗 (𝜎 𝑇 𝑐 𝑟 ) 𝑚𝑎𝑥 > 𝑅 𝑓 [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF] Note that initially, (𝜎 𝑇 𝑐 𝑟 ) 𝑚𝑎𝑥 is set to an appropriately chosen reference value in each cell, and then it is updated if the product of (𝜎 𝑇 𝑐 𝑟 ) of the chosen pair becomes larger than the reference value. This procedure keeps the computational expanse of NTC as O(N). The sequence of NTC collision procedure is shown in Fig. 4. First, an upper limit, i.e., N coll pairs to be checked is calculated; second, by using acceptance-rejection procedure the actual number of accepted collisions in a cell is defined. It should be reminded that in Bird's monograph [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF], and other references describing NTC, instead the term 𝑁(𝑁 -1) was written the product 𝑁𝑁 ̅ , where over bar quantity means averaged value.

The relation of the collision process to the Poisson distribution and its consequence, i.e., equality

N coll
• Compute maximum number of collision pairs (Eq. ( 9)) pair (i,j)

• Selection of a collision pair (i,j) randomly from N l particle in cell l.

if accepted

• checking the collision probability (Eq. ( 10))

Velocity update

•change the particle velocities if collision is accepted of 𝑁 ̅ (𝑁 -1) ̅̅̅̅̅̅̅̅̅̅ = 𝑁 2 ̅̅̅̅ , were pointed out by Yanitskiy [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF], then reemphasized and utilized, as given in Eq. [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF], by Stefanov and Cercignani [START_REF] Stefanov | Monte Carlo simulation of Bénard's instability in a rarefied gas[END_REF], and finally utilized by Bird in his sophisticated algorithm [START_REF] Bird | Sophisticated versus simple DSMC[END_REF][START_REF] Bird | Sophisticated DSMC: Notes prepared for a short course at the DSMC07 meeting[END_REF].

An insufficient number of particles per cell can be a source of stochastic errors in the NTC model.

However, if enough number of particles in cells is employed, NTC scheme is an efficient approach for modeling the intermolecular collisions. But, NTC scheme is CPU demanding in complex, 3D simulations where a huge number of cells and particles are required for proper flow simulation. In fact, NTC scheme requires enough number of particles per cell, i.e., N=10-20, to provide accurate results. However, due to random selection of collision pairs, one of the main shortcomings of the NTC collision is the possibility of repeated collisions, i.e., the same pair of particles is selected repeatedly for collision within one or several successive time steps without occurrence of events of collisions with other particles. This results in inaccuracy of the collision process in cells containing small number of particles. One remedy is to construct an array and keeping track of the record of collision pairs to avoid repeated collisions for a particle pair.

However, this introduces additional complexities. To use the acceptance-rejection procedure, the pairs should be chosen randomly from the whole available set of N particles to keep a correct probability for collision corresponding to the Boltzmann collision frequency. In this case, all probabilities and collision frequencies are calculated correctly. Direct avoiding of repeated collisions means that one introduces a condition that leads to a non-uniform distribution of probabilities, because for the first chosen pair one chooses randomly among N(N-1)/2 pair of particles but the second choice is performed from among (N)(N-1)/2-1 pairs of particles. The effect is stronger if one keeps this direct collision control condition for successive time steps.

Evidently, this error is minor for large number of particles but it is of considerable importance if there are a few particles per cell. There is another disadvantage in the NTC scheme, the use of remainders to keep an accurate number of collisions in time. The problem arises from the practical calculation of N coll by using eq. ( 9), which, in general case, is not an integer value and the integer part of N coll is used. When the number N coll is large it works fine but when N coll ~ 1 or less the reminder procedure introduces additional error related to the fact the probability distribution of the remainders gives different from the Poisson law variation of the number of collisions.

Nearest Neighbor (NN) scheme

LeBeau et al. [START_REF] Lebeau | Virtual sub-cells for the direct simulation Monte Carlo method[END_REF] introduced the "virtual subcell" (VSC) method as an improvement in the particle selection procedure of the NTC, in which the first particle is chosen randomly, while the second partner is selected from available neighbors of the selected particle, i.e., the nearest neighbor (NN) of a given particle is found and selected. The flowchart of NN collision model is shown in Fig. 5. Method VSC needs the collision separation distance to be calculated for all available particles in the cell to find the nearest neighbor one. To reduce computational costs to O(N), Bird proposed to avoid the calculation of all intermolecular distances at the start of collision procedure and restrict to the calculation of an intermolecular distance between the first randomly chosen particle, i, and the nearest among the other particles.

The NTC-NN algorithm was implemented in the new versions of DSMC codes distributed by Gream Bird [START_REF] Bird | The DSMC Method[END_REF]. Bird et al. [START_REF] Bird | Accuracy and efficiency of the sophisticated direct simulation Monte Carlo algorithm for simulating noncontinuum gas flows[END_REF] and Gallis et al. [START_REF] Gallis | Effect of collision-partner selection schemes on the accuracy and efficiency of the direct simulation Monte Carlo method[END_REF] compared the convergence behavior of the NN collision model compared with the standard NTC one, and showed that NN achieves high efficiency as it minimizes the mean collision separation between collision partners; however, the new algorithm is very sensitive to the selection of the time step and almost for all cell needs smaller time steps compared to the original NTC. They showed that NN significantly reduces the computational resources required for a DSMC simulation to achieve a particular level of accuracy if its computational set-up is chosen appropriately. As a modification of NN model, Macrossan [START_REF] Macrossan | Searching for a near neighbor particle in DSMC cells using pseudo-subcells[END_REF] suggested the 'pseudo-subcell' collision method, in which the search for a collision partner discontinues if a 'near-enough' particle is found, i.e. whenever another particle is found within the 'pseudo-subcell' of radius  centered on the first particle. Macrossan showed that in structured cells, the pseudo-subcell method gives a 5% increase in mean collision separation (MCS) compared to the NN method, and is up to 20% faster than the NN method.

As mentioned, NN needs lower time steps as it ignores a key point in pair selection: within a time step particles with higher velocities may reach and collide particles, which are further than their nearest neighbor; therefore, a small enough time step should be utilized. On the other hand, according to the Boltzmann equation, the working space is phase space with three coordinates (positions) and three velocities variables. From this view point, NN takes into account only the space of coordinates using the criterion of minimum distance for colliding pairs. This means that NN works well when the whole velocity distribution function is available at every coordinate point. To have this fulfilled, one should have enough particles in all small volumes. Otherwise, there will be negative effects and erroneous solutions. Bird [START_REF] Bird | Sophisticated DSMC: Notes prepared for a short course at the DSMC07 meeting[END_REF] showed that NN needs at least 7 particles per collision cell to predict equilibrium collision frequency correctly.

It should be noted that Bird states that lower values of separation of free paths (SOF), which is the ratio of mean collision separation (MCS) to the mean free path (SOF=MCS/) largely selfvalidates the accuracy of DSMC solutions [START_REF] Bird | The DSMC Method[END_REF]. However, it should be mentioned that this is a necessary but not a sufficient condition. So, it should be noted that low SOF did not always guarantee correct solution if it is not accompanied with enough particle per cells (PPC) and fine enough cell.

Fig. 5 NN collision procedure

A note on DSMC cells and subcells

The DSMC method uses the cell for the selection of collision partners and for the sampling of the macroscopic properties. To avoid statistical scatters and reduce inherent errors in the standard NTC method, the number of DSMC particles per cell should be as large as possible, generally around twenty. On the other hand, in the selection of collision pairs, it is desirable to choose close particles and reduce the mean collision separation of pairs and minimize the smearing of gradients. These conflicting requirements can be remedied by dividing the sampling cell into a set of sub-cells and perform pair selection and collision from inside the subcells [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. However, the first use of the subcells was suggested after the objection raised by Meiburg [START_REF] Meiburg | Comparison of the molecular dynamics method and the direct simulation technique for flows around simple geometries[END_REF]. As angular momentum is not conserved in a simulated collision, Meiburg [START_REF] Meiburg | Comparison of the molecular dynamics method and the direct simulation technique for flows around simple geometries[END_REF] concluded that if colliding particles are selected from opposite sides of a cell, the collision introduces significant error.

Meiburg pointed out this process in a flow with vorticity and noted that the angular velocity of the collision partners about the center of the collision is not necessarily conserved. However, Bird

N coll

• Compute number of collision pairs (Eq. ( 9)) pair (i,j)

• Select the first collision particle (i) randomly from N l particle in cell l. The second particle (j) is the closest one to particle (i) if accepted

• checking the collision probability (Eq. ( 10)) Velocity update

•change the particle velocities if collision is accepted mentioned that this is essentially a cell size effect, i.e., the cell size must be very small in comparison with the mean free path in regions with large gradients, such as Knudsen layers. As Bird mentioned, the problem with Meiburg's calculation was that the cells he had used had a linear dimension of 3 which was excessively large by a factor of at least ten [START_REF] Mareschal | Microscopic simulations of complex flows[END_REF]. However, to overcome this problem, Bird [START_REF] Bird | Direct simulation of high-vorticity gas flows[END_REF] introduced subcells within each cell, where a potential collision partner is selected from the same subcell or from a nearby subcell if no collision partner could be found in the same subcell. Thus, the use of subcells reduces the mean collision separation (MCS) to a fraction of the subcell size, rather than a fraction of the cell size. Stefanov et al. [START_REF] Stefanov | Monte Carlo analysis of macroscopic fluctuations in a rarefied Hypersonic flow around a cylinder[END_REF][START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. I. Convection regime[END_REF] suggested the idea of dynamic subdivision of subcells, or as called: collision cells. They adjusted the subcell such that the subcell size remains smaller than the mean free path (λ) at every time step. In order to provide a sufficient number of subcells that is compatible with the instantaneous number of particles inside each collision cell, Bird proposed the transient adaptive subcell (TAS) technique [START_REF] Bird | Visual DSMC program for two-dimensional and axially symmetric flows[END_REF][START_REF] Bird | The DS2V/3V program suite for DSMC calculations[END_REF]. TAS technique was successfully employed in combination with NTC collision model in structured [START_REF] Bird | Visual DSMC program for two-dimensional and axially symmetric flows[END_REF] and unstructured cells [START_REF] Su | Implementation of a transient adaptive sub-cell module for the parallel-DSMC code using unstructured grids[END_REF]. In fact, before reaching the steady state, the number of particles in each cell varies during the computational process; consequently, the number of transient subcells varies according to local density gradients, see Fig. 6. To determine the number of subcells in each direction in a 2-D cell, number of particles per subcell, PPSC, should be set as one of the inputs. For each cell, this transient layer of subcells is fabricated by a special number of divisions along x (D x ) and y (D y ) directions as follows:

𝐷 𝑥 = √ 𝑁 𝑐 𝐴𝑅 × 𝐸𝑛𝑣𝑒𝑙𝑜𝑝 × 𝑃𝑃𝑆𝐶 ( 11
)
𝐷 𝑦 = 𝐷 𝑥 × 𝐴𝑅 ; 𝐴𝑅 = ∆𝑦 𝑐𝑒𝑙𝑙 ∆𝑥 𝑐𝑒𝑙𝑙 ; 𝐸𝑛𝑣𝑒𝑙𝑜𝑝 = V 𝑐 ∆𝑦 𝑐𝑒𝑙𝑙 × ∆𝑥 𝑐𝑒𝑙𝑙 ( 12
)
, where 𝐴𝑅 is the aspect ratio, 𝑃𝑃𝑆𝐶 is the desired number of particles per subcell and 𝐸𝑛𝑣𝑒𝑙𝑜𝑝 is used to increase the number of subcells, in case the cell is not fully rectangular. Using TAS, the selection of collision pairs is performed in subcells level, therefore, the effect of cell-size on the accuracy of the solution is reduced and much coarser basic collision cells could be employed.

However, to apply correctly the TAS strategy one should complete it with a second element concerning the time step. The time step should be coupled with the smallest subcell size in order to fulfil the well-known Courant condition, which for the DSMC method could be replaced with a simple rule required on the stage of particle free motion i.e. the time step to be chosen small enough in order the probability for all particles to pass further than one neighboring subcell to tend to zero. As shown in section 3.9.2 this point was taken into account for controlling the time step.

In the sophisticated version of DSMC, the idea of using cells and subcells was upgraded to use separate sampling and collision cells. In the NN collision model, collision cells should be adapted to 7 particles and sampling cells are adapted to around 20 particles, respectively [START_REF] Bird | Sophisticated DSMC: Notes prepared for a short course at the DSMC07 meeting[END_REF].

Kac Master Equation

Introduction

In the homogenous case, i.e. a kinetic equation without streaming term, the time evolution of velocity distribution function of N particle system 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) due to inter-molecular collision could be described by the Kac stochastic model [START_REF] Kac | Probability and related topics in physical sciences[END_REF][START_REF] Kac | Probability and Adjacent Questions[END_REF][START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF]. The Kac master equation is a probabilistic analog of Liouville equation considering only binary collisions of the particles with given probability. Consider a system of { 𝑥 (𝑙) , 𝐶 𝑁 (𝑙) } = { 𝑥 𝑗 (𝑙) (𝑡 𝑘 ), 𝑐 𝑗 (𝑙) (𝑡 𝑘 ) }, 𝑗 = 1, . . . , 𝑁 (𝑙) particles in a cell (l) with volume 𝑉 (𝑙) . The Kac stochastic model can be described by the following direct Kolmogorov equation:

𝜕 𝜕𝑡 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) = 1 𝑉 𝑙 ∑ 𝑔 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙 ∫ [𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) 𝑖𝑗 ) -𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) )] 4𝜋 𝑑𝜎 𝑖𝑗 (13)
, where 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) 𝑖𝑗 ) is the N-particle probabilistic velocity distribution function after the collision of i-j particles, i.e., 𝑐 𝑁 (𝑙) 𝑖𝑗 = {𝑐 1 , … . . , 𝑐 𝑖-1 , 𝑐 𝑖 ′ , 𝑐 𝑖-1 , … . , 𝑐 𝑗-1 , 𝑐 𝑗 ′ , 𝑐 𝑗+1 , … . , 𝑐 𝑁 (𝑙) }, which denotes that velocities of pair i-j particles are updated to their post-collision velocities. The sum over 1 ≤ 𝑖 < 𝑗 ≤ 𝑁 𝑙 means summation over N(N-1)/2 collision pairs. Kac stochastic model is a linear integro-differential equation that describes the time behavior of the N-particle distribution function. In fact, The Kac stochastic model is a jump-like strictly Markovian process over the hypersphere Ω(𝑁, 𝐸, 𝑃) of 3N-4 dimensions in Euclidean space R 3N , where N is the number of particles, E is the kinetic energy and P is the momentum of velocity components. This hypersphere is formed by crossing of a hypersphere of total kinetic energy of particles 𝑐 1 2 + ⋯ . +𝑐 𝑛 2 = 𝐸 = Cte with three hyperplanes 𝑐 1 + ⋯ . +𝑐 𝑛 = 𝑃 = Cte. Strictly Markovian processes are described by a pair of adjoined equations, i.e., direct and reverse Kolmogorov equations [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF].

The one particle distribution function could be obtained from 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) 𝑖𝑗 ) as follows:

𝑓 1 (𝑛) (𝑐 1, 𝑥 (𝑙) , 𝑡) = ∫ 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) )

𝑐 2 2 +⋯+𝑐 𝑛 2 =𝑛-𝑐 2 𝑑𝜎 (14)
, i.e., the integration over all velocity classes, expect to that of c 1, results in a reduced distribution function of velocity class c 1 . Assuming number of particles approaches infinity and using the local assumption of molecular chaos, which permits to write two-particle distribution function in terms of one particle distribution function:

𝑓 2 (𝑐 1 , 𝑐 2 , 𝑥 (𝑙) , 𝑡) = 𝑓 1 (𝑐 1 , 𝑥 (𝑙) , 𝑡)𝑓 1 (𝑐 2 , 𝑥 (𝑙) , 𝑡)

, the Boltzmann equation without the streaming term could be recovered from the Kac stochastic model.

In general, the hypothesis of molecular chaos is one of the most discussed assumptions in kinetic theory and statistical physics. There are not well-defined criteria for its validity. However, it can be accepted as a good approximation that the molecular chaos hypothesis is valid if the gas is close to a local equilibrium state and the provided number of particles, N, becomes very large. It is obvious that the Kac master equation does not require the fulfillment of molecular chaos hypothesis although it is partly used for generating stochastic collision parameters. One can conclude that all collision schemes derived from the Kac master equation keep this property.

Following this line, it could be shown that most of collision schemes, except the Nanbu collision scheme [START_REF] Nanbu | Direct Simulation Scheme Derived from the Boltzmann-Equation .I. Monocomponent Gases[END_REF], have a more direct connection to the Kac master equation than the Boltzmann equation. In fact, similar to Liouville equation, Kac model contains all possible correlation between particles velocities. Reducing N-particle distribution function via BBGKY chain of equations one consequently neglects (N-1)-particle correlation functions. The last is all binary correlation functions, among them are self-correlations functions that means velocity and other kind of fluctuations. The last relation between velocity distribution functions is between twoparticle and one-particle distribution functions, which reads as:

𝑓 2 (𝑐 1 , 𝑐 2 , 𝑥 (𝑙) , 𝑡) = 𝑓 1 (𝑐 1 , 𝑥 (𝑙) , 𝑡)𝑓 1 (𝑐 2 , 𝑥 (𝑙) , 𝑡) + 𝐶𝑜𝑟𝑟(𝑐 1 , 𝑐 2 )

Where 𝐶𝑜𝑟𝑟(𝑐 1 , 𝑐 2 ) corresponds to correlations between distribution functions of 𝑐 1 and 𝑐 2 .

Molecular chaos assumption means 𝐶𝑜𝑟𝑟(𝑐 1 , 𝑐 2 ) = 0, which leads to equality [START_REF] Bird | Shock Wave Structure in a Rigid Sphere Gas[END_REF] employed in the derivation of the Boltzmann equation. Thus, the correlations do not exist in the Boltzmann equation solution. But DSMC can evaluate them realistically as shown in the Bird's monograph [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] and supported, for example, by numerical results shown in Ref. [START_REF] Stefanov | Monte Carlo analysis of macroscopic fluctuations in a rarefied Hypersonic flow around a cylinder[END_REF] obtained for velocity correlations and fluctuations showing that the DSMC method adequately evaluates binary correlation functions, which cannot be obtained from the Boltzmann equation. This circumstance determines the place of the DSMC method in flowchart of hierarchy of the governing equation shown in Fig. 1. However, it should be noted that the calculation of the maximum collision rate used in collision schemes based on TC and NTC collision algorithms is founded on the Eq. ( 4) of the classic Boltzmann kinetic theory, which reflects on the classification of the models given in Fig. 2.

In the operator form, Kac stochastic model, Eq. ( 13), could be written as follows [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF]:

𝜕 𝜕𝑡 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) = [ ∑ 𝑤 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙
(𝑇 𝑖𝑗 -𝐼)] 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) )

= 𝜐(𝑇 -𝐼)𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) )

, where

𝐼𝜓 = 𝜓, 𝑇 𝑖𝑗 = ∫ 𝜓(𝑐 𝑖𝑗 )𝐵(𝑔 𝑖𝑗 , 𝜃)𝑑Ω(𝜃), 4𝜋 𝑇𝜓 = ∑ 𝑤 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙 𝑇 𝑖𝑗 𝜓 (18) (19)
, and 𝐵(𝑔 𝑖𝑗 , 𝜃) is a scattering kernel [START_REF] Deshpande | Monte Carlo simulation for molecular gas dynamics[END_REF]. The transformation 𝑇 𝑖𝑗 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) results in probability density distribution of 𝑐 𝑁 (𝑙) 𝑖𝑗 over Ω(𝑁, 𝐸, 𝑃) after the collision of a pair (𝑐 𝑖 , 𝑐 𝑗 ) under the condition that before the collision point 𝑐 𝑁 (𝑙) is distributed over Ω(𝑁, 𝐸, 𝑃) with the probability density 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ).

In the Kac model, the probability of time interval between two successive collisions is distributed according to the exponential distribution of a Poisson process as follows:

𝑃𝑟𝑜𝑏(𝛿𝑡 > 𝑡) = 𝑒 -𝜐𝑡 𝜐 = ∑ 𝑤 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙 ; 𝑤 𝑖𝑗 = 𝜎 𝑖𝑗 𝑔 𝑖𝑗 𝑉 𝑙 (20)
, where collision probability of pair (i, j) as denoted by (𝑤 𝑖𝑗 ) is a function of σ ij , which is the collision cross-section, and g ij = |c i -c j |, which is the particles relative velocities.

In case of the given state at time t 0 , the operator form of the Kac master equation could be solved at time 𝑡 with the definition of the transition operator 𝐺(𝑡) in the following form:

F N (l) (t, x (l) , C N (l) ) = G(t) F N (l) (t 0 , x (l) , C N (l) ) (21)
Where the transition operator G(t) could be written as:

𝐺(𝑡) = 𝑒𝑥𝑝 [𝑡 ∑ 𝑤 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙 (𝑇 𝑖𝑗 -𝐼)] = exp [𝑡𝑣(𝑇 -𝐼)] (22)
In the next sections, following Yanitskiy [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF], a family of collision models derived from Kac stochastic model and built from different variations of the exponential transitional operator formula, Eq. ( 22), will be discussed. These collision models define a collision probability function for each particle pair and checks all pair combinations for collision. In fact, the modeling of collision is converted to a statistical realization of the evolution of Kac's model during a period of time.

A Pseudo-Poisson Process (Collision Frequency Technique)

Considering small time step, we could assume: 𝜐 = 𝜐 𝑐 + 𝑂(𝑡), neglecting the higher order terms, the operator G is approximated as follows [START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF]:

𝐺 1 (𝑡) =
The above-mentioned algorithm is exact for pseudo-Maxwellian molecules, where the collision probability does not depends on relative velocity, and it has an order of accuracy of O(t) in the general case. Frequency 𝜐 𝑐 is computed at each time step and order of computations is O (N 2 ).

Sampling of the number of collisions s(t) is performed according to probability distribution P(𝑡𝜐 𝑐 ). A pair of colliding particles is sampled with the probability of 𝑊 𝑖𝑗 .

For non-Maxwellian gases, the following algorithm is performed [START_REF] Koura | Transient Couette flow of rarefied binary gas mixtures[END_REF] (see Fig. 7):

1-Collision time for a pair of (i, j) with velocities (V, V 1 ) is sampled form the exponential distribution with the collision frequency 𝜐 𝑐 as follows:

𝑃𝑟𝑜𝑏(𝛿𝑡 𝑐 ) = 𝜐 𝑐 𝑒 -𝜐 𝑐 𝑡 (25)
In an inverse sampling method widely employed in the Monte Carlo method, the sampled time step is given by:

𝛿𝑡 𝑐 = -𝑙𝑛(𝑅 𝑛𝑓 ) 𝑣 𝑐 (26)
The above procedure is continued until ∑ 𝛿𝑡 𝑐 exceeds cell time step. The pseudo-Poisson process was applied for the first time by Koura in the collision model [START_REF] Koura | Transient Couette flow of rarefied binary gas mixtures[END_REF] for counting the collision of a non-Maxwellian gas, where the collision frequency was computed as follows:

𝜐 𝑐 ≈ 𝐹 𝑁 𝑉 𝑐 ∑ ∑ (𝜎 𝑖𝑗 𝑔 𝑖𝑗 ) 𝑗=1,≠𝑖 𝑁 𝑖=1 (27)
However, in this case one needs to choose a pair for collision with a given probability. For this purpose, Koura derived it from the Boltzmann assumption of molecular chaos (for details see [START_REF] Koura | Transient Couette flow of rarefied binary gas mixtures[END_REF]). This fact led us to link the Koura's approach more to the Boltzmann equation, than the Kac stochastic model.

It should be reminded that 𝜎 𝑖𝑗 is a function of relative velocity and other gas properties, i.e. for variable had sphere gases, it is given by:

𝜎 𝑖𝑗 = 𝜋𝑑 𝑟𝑒𝑓 2 1 Γ(2.5 -𝜔) ( 2𝑘𝑇 𝑟𝑒𝑓 𝑚 𝑟 𝑔 2 ) 𝜔-0.5 ( 28
)
Where 𝜔 is the viscosity-temperature exponent: it is equal to 0.5 for hard sphere gases and is equal to 1 for Maxwell molecules [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]. Main drawbacks of the collision frequency technique are the significant time required for computation of Eq. ( 27) and the assumption of constant 𝜐 𝑐 during the time step; which requires the condition of 𝜐 𝑐 /𝜐 𝑐 ≪ 1where 𝜐 𝑐 is the difference between 𝜐 𝑐 at two successive time steps. • choose a pair for collision with a given probability, i.e., see [START_REF] Koura | Transient Couette flow of rarefied binary gas mixtures[END_REF] Velocity update

•change the particle velocities and advance time step until cell time step is reached.

Majorant frequency scheme (MFS) collision scheme

Majorant frequency scheme (MFS) could be derived either from the Leontovich master kinetic equation (MKE), which describes the behavior of an N-particle gas model with binary collisions [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF][START_REF] Ivanov | Direct statistical simulation method and master kinetic equation[END_REF][START_REF] Ivanov | Comparative analysis of algorithms of direct statistical method in rarefied gas dynamics[END_REF], or from the Kac master equation [START_REF] Ivanov | Statistical simulation of reentry capsule aerodynamics in hypersonic near-continuum flows[END_REF]. Similar to Kac master equation, Leontovich equation may be transformed to the Boltzmann equation as N→ ∞ and the molecular chaos condition is satisfied. Full details of the derivation of the MFS scheme from the Kac master equation is provided in Ref. [START_REF] Ivanov | Statistical simulation of reentry capsule aerodynamics in hypersonic near-continuum flows[END_REF]. In the MFS scheme, the time between two eventual collisions is computed from a Poisson distribution with the maximum (majorant) probability value as follows:

𝑃𝑟𝑜𝑏(𝛿𝑡 𝑐 ) = 𝜐 𝑚𝑎𝑥 𝑒 -𝜐 𝑚𝑎𝑥 𝑡 ( 29
)
𝛿𝑡 𝑐 = -𝑙𝑛(𝑅 𝑛𝑓 ) 𝑣 𝑚𝑎𝑥 (30)
The collisional pair is uniformly chosen from 𝑁(𝑁 -1)/2 available pairs. The majorant frequency is defined as:

𝜐 𝑚𝑎𝑥 = 1/2𝑁(𝑁 -1)𝐹 𝑁 (𝜎 𝑇 𝑐 𝑟 ) 𝑚𝑎𝑥 /𝑉 𝑐 (31)
MFS scheme starts with calculating this frequency. MFS performs as many collisions in one time step until the summation of computed 𝛿𝑡 𝑐 's exceeds the specified cell time step. The procedure of MFS collision is shown in Fig. 8. More details of MFS scheme and its validation could be found in Refs. [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF][START_REF] Ivanov | Statistical simulation of reentry capsule aerodynamics in hypersonic near-continuum flows[END_REF]. It is worth noting that the choice of the collision frequency in form Eq. [START_REF] Khlopkov | Development of Monte Carlo Methods in Hypersonic Aerodynamics[END_REF], which contains an assumption for statistical independency of molecular velocities, and respectively, is based on the molecular chaos hypothesis, is important for understanding the stochastic properties of the MFS collision algorithm and its close relation with the Bird's NTC scheme. The MFS model improves the probabilistic characteristics of the process by including the Poisson process for determining the maximum number of pairs of particles. Unlike Bird's NTC, in MFS it is an integer random number obeying the Poisson law and its averaged value within a time step is exactly equal to the calculated in NTC maximum number of particle pairs to check for collision by using Eq. ( 9). Thus, the MFS algorithm avoids the use of reminders, included in NTC to balance the integer number of collisions during the simulation.

Fig. 8 The Majorant Frequency Scheme procedure per each cell.

Bernoulli Trials [9]

For small interval 𝑡, operator G(t) can be expanded with respect to 𝑡 degrees and terms of order equal or higher than O(t 2 ) are neglected. The transition operator G(t) corresponding to this approximation is obtained as:

𝐺 2 (𝑡) ≈ ∏ 𝑒𝑥𝑝 1≤𝑖<𝑗≤𝑁 𝑙 [𝑡(𝑇 𝑖𝑗 -𝐼)𝑤 𝑖𝑗 ] = ∏ ∏ exp[𝑡(𝑇 𝑖𝑗 -𝐼)𝑤 𝑖𝑗 ] 𝑁 𝑗=𝑖+1 𝑁-1 𝑖=1 (32)
One could replace the exponential term in every co-factor in Eq. ( 32) by a linear approximation over t, thus:

•Calculation of the majorant frequency from Eq. ( 31), [sum time step]0

•Calculating the time to the next eventual collision from the Poisson distribution using Eq. (30)

• If [sum time step]=[sum time step]+time step exceeds DSMC time step (dt), end the loop, otherwise

•choose collision pair (i,j) randomly from N l particle in cell l.

•checking the collision probability (Eq. ( 10)), Note: Employed random number should be different from that used in Eq. (30)

•change particle velocities if collision is accepted

𝐺 2 (t) = ∏ ∏ [(1 -𝑡𝑤 𝑖𝑗 )I + 𝑡𝑤 𝑖𝑗 T ij ] N (l) j=i+1 N (l) -1 i=1 = ∏ ∏ [(1 -W ij )I + W ij T ij ] N (l) j=i+1 N (l) -1 i=1 (33)
, where

W ij = 𝜎 𝑖𝑗 𝑔 𝑖𝑗 𝑡 𝑉 (𝑙) . (34)
Every co-factor in the right-hand-side of Eq. ( 33) transforms the probability distribution function 𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) over Ω(𝑁, 𝐸, 𝑃) into 𝐹 𝑁 (𝑙) ̅̅̅̅̅̅ (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) over Ω(𝑁, 𝐸, 𝑃) according to:

𝐹 𝑁 (𝑙) ̅̅̅̅̅̅ (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) = (1 -tw ij )𝐹 𝑁 (𝑙) (𝑡, 𝑥 (𝑙) , 𝑐 𝑁 (𝑙) ) + tw ij ∫ 𝜓(𝑐 𝑖𝑗 )𝐵(𝑔 𝑖𝑗 , 𝜃)𝑑Ω(𝜃) 4𝜋

If the time step t is considered small such that:

𝑝 𝑖𝑗 = 𝑡𝑤 𝑖𝑗 ≤ 1 (36)
For all possible values of 𝑐 𝑁 (𝑙) , Eq. ( 33) may have an apparent probability approach. The collision scheme based on Eq. ( 33) considers all N(N-1)/2 pairs of particles (𝑐 𝑖 , 𝑐 𝑗 ) and to consider the collision of each pair of (𝑐 𝑖 , 𝑐 𝑗 ) with probability 𝑝 𝑖𝑗 . This scheme is of first order accuracy in terms of time but its computational costs is of order N 2 . The condition given by Eq. ( 36) could be replaced with a weaker one as follows:

𝑃{𝑡𝑤 𝑖𝑗 > 1} ≪ 1 (37)
by choosing appropriate time step and cell size. The above procedure states that for all of the available particle pairs in the collision cell (l), the acceptance-rejection should be checked, i.e. the following inequality should be checked for all available particle pairs(i,j) {i < j = 1, . . . , N (l) } (before their velocities are changed to post collision values):

W ij = σ ij g ij t V (l) > 𝑅 𝑛𝑓 (38)
The procedure of SBT is depicted in Fig. 9.

Ballot Box collision scheme [9]

If the transition operator is approximated by a linear term in terms of t as:

𝐺 3 (𝑡) = 𝐼 + 𝑡 ∑ 𝑤 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙 (𝑇 𝑖𝑗 -𝐼) (39)
, the number of particle pairs is denoted by k= N(N-1)/2, and if one introduces the time step (t) inside the summation operator into the operator G 3 (t), replacing 𝑡𝑤 𝑖𝑗 by 1 k (𝑘𝑡𝑤 𝑖𝑗 ), then:

G 3 (t) = [1 -∑ 1 k (𝑘𝑡𝑤 𝑖𝑗 ) 1≤𝑖<𝑗≤𝑁 𝑙 ]I + ∑ 1 k (𝑘𝑡𝑤 𝑖𝑗 )T ij 1≤𝑖<𝑗≤𝑁 𝑙 (40)
This form of operator G 3 (t) can be interpreted in a straightforward way that lead to the Ballot Box collision model. The condition 𝑃 𝑖𝑗 {𝑘𝑡𝑤 𝑖𝑗 > 1} ≪ 1 should also be satisfied. In the Ballot Box scheme, the following procedure is performed one times for each cell (see Fig. 10):

Inexing and first particle selection

• Local cross referencing of particle is performed in the cell •choose first particle in sequense from list: i=1, ... ,N l -1 second particle selection

• Choose the second particle from the list of indexed particles, i.e., j=i+1, ..., N l if accepted

• checking the collision probability (Eq. ( 38))

Velocity update

•change the particle velocities if collision is accepted a) At each time step, only a single pair of molecules is randomly selected with the probability of 1/k from the whole range of particles.

b) The possibility of collision is evaluated with the probability of 𝑊 𝑖𝑗 = 𝑘𝑡𝑤 𝑖𝑗 .

W ij = N(N -1)σ ij g ij t 2V (l) > 𝑅 𝑛𝑓 (41)
If the collision was accepted, the velocities are upgraded to their post-collision values. The laboriousness of the scheme depends linearly on the number of particles, i.e., O(N). The problem with the Ballot Box scheme is that 𝑡 ≪ ∆𝑡 (∆𝑡 is the optimal time step of simulation) in order to fulfill condition {P(𝑘𝑡𝑤 𝑖𝑗 > 1)}<<1. The collision algorithm Ballot-Box must be repeated 𝑚 times in time step ∆𝑡, so that ∆𝑡 = 𝑚𝑡. As a result, within the time step ∆𝑡 the scheme cannot prevent from repeated collisions. • checking the collision probability (Eq. 41)

Velocity update

•change the particle velocities if collision is accepted

Null-collision (maximum frequency approach)

This method is first derived by Koura [START_REF] Koura | Null-collision technique in the direct-simulation Monte Carlo method[END_REF] for pseudo-Maxwellian gases and then extended and modified to arbitrary gases, (g) [START_REF] Ivanov | Theoretical analysis of traditional and modern schemes of the DSMC method[END_REF] (and called majorant frequency therein). In the Nullcollision scheme, the maximum evaluation of the 𝜐 in (N, E, P) is computed in terms of 𝜐 𝑚𝑎𝑥 (N, E, P). The collision operator T is replaced with T max while it keeps the condition of:

𝜐 0 (𝑇 𝑚𝑎𝑥 -𝐼) = 𝜐(𝑇 -𝐼), therefore, the transition operator is written as follows:

G 4 (t) = exp[𝑡𝜐 𝑚𝑎𝑥 (𝑇 𝑚𝑎𝑥 -𝐼)] 𝑇 𝑚𝑎𝑥 = (1 - 𝜐 𝜐 𝑚𝑎𝑥 ) 𝐼 + 𝜐 𝜐 𝑚𝑎𝑥 𝑇 (42) (43)
As the maximum value 𝜐 𝑚𝑎𝑥 is a constant, G 4 (t) defines a transition operator of a pseudo-Poisson process where the number of collisions are distributed according to

𝑃 𝑠 (𝜐 𝑚𝑎𝑥 𝑡) = (𝑡𝜐 𝑚𝑎𝑥 ) 𝑠 𝑠! exp (-𝑡𝜐 𝑚𝑎𝑥 ) (44)
Any pair selected for collision possibility has a collision chance of

𝜐 𝜐 𝑚𝑎𝑥

and collision rejection chance of (1 -

𝜐 𝜐 𝑚𝑎𝑥

). Similar to Ballot-Box scheme and to introduce an equal opportunity for selecting a pair for collision, the maximum frequency value, 𝑤 𝑚𝑎𝑥 is computed. Values of 𝜐 𝑚𝑎𝑥 and 𝑤 𝑚𝑎𝑥 are chosen such that: 𝜐 𝑚𝑎𝑥 = 𝑘𝑤 𝑚𝑎𝑥 , where k is the number of particle pairs, then 𝑇 𝑚𝑎𝑥 is transformed to:

𝑇 𝑚𝑎𝑥 = ∑ 1 𝑘 1≤𝑖≤𝑗≤𝑁 [(1 - 𝑤 𝑖𝑗 𝑤 𝑚𝑎𝑥 ) 𝐼 + 𝑤 𝑖𝑗 𝑤 𝑚𝑎𝑥 𝑇 𝑖𝑗 ] (45)
The numerical procedure of Koura's method is as follows:

Repeat the procedure below s-times [START_REF] Koura | Null-collision technique in the direct-simulation Monte Carlo method[END_REF]:

1-Compute maximum frequency as follows:

𝜐 𝑚𝑎𝑥 = ( 1 2 )𝑛𝑁(𝑆) 𝑚𝑎𝑥 (46)
𝑆 𝑚𝑎𝑥 = σ(𝑔) 𝑚𝑎𝑥 [START_REF] Bird | The DS2V/3V program suite for DSMC calculations[END_REF] 2-A random pair (i,j) is selected with the same probability from all particles.

3-The collision is accepted with the probability of pair (i,j)

•Select the particle pais (i, j) randomly from N l particle in cell l.

if accepted

• checking the collision probability (Eq. ( 10))

Velocity update

•change the particle velocities if collision is accepted

Simplified Bernoulli-trials

Bernoulli-trials (BT) scheme checks collisions with a given probability for all available collision pairs, i.e., N(N-1), which makes this algorithm quite time consuming, i.e. O(N 2 ). Recently, Stefanov [START_REF] Stefanov | Particle Monte Carlo algorithms with small number of particles in grid cells[END_REF][START_REF] Stefanov | On DSMC calculations of rarefied gas flows with small number of particles in cells[END_REF] showed that in small time-step intervals, it is possible to simplify the collision process and instead of checking all available possibilities, for the first selected as shown in BT particle, say i, one can choose randomly the second particle among the next sorted-particles in the list. Hence the number of selections would be reduced to (N-1) while theoretically the probability for collision of each pair remains unchanged based on the strict Markov process equal to the corresponding one in the BT method. Thus, in SBT each chosen particle pair is checked for collision only once within a time step and possibility of repeated collision is prevented while the computational costs are much reduced compared to the BT method.

As simplification, Stefanov showed that it is possible to extend the internal product in the right hand side of the G(t) (Eq. 22) in a series of j with respect to 𝑡 to reach to a new simplified transition operator G 5 (t) as follows [START_REF] Stefanov | Particle Monte Carlo algorithms with small number of particles in grid cells[END_REF]:

G 5 (t) = ∏ [(1 -∑ 1 k (𝑘𝑡𝑤 𝑖𝑗 ) N (l) j=i+1 ) I + ∑ 1 k ((𝑘𝑡𝑤 𝑖𝑗 )T ij ) N (l) j=i+1 ] N (l) -1 i=1 (49)
, where k = (N (l) -i). The algorithmic interpretation of operator G 5 (t) states that instead of checking N(N-1) pairs (Note that the inner product in Eq. ( 33) considering all possible collision pairs N(N-1) is replaced here by a summation), it is possible in each acceptance-rejection, select the second collision particle at random and reduce the number of collision checking to (N-1). The SBT algorithm permits simulations with far less number of particles per cell compared to NTC, < N > ~ 1-2, where < > means mean value, and with reduced computational costs compared to BT algorithm. The numerical procedure of SBT is as follows: Particles in the l th cell should be locally indexed to produce a particle list numbered as 1…N l . The first particle of the collision pair (i, j), say i, is selected in sequence from the following particle list: i = 1…N l -1. The second particle, say j, is chosen randomly with the probability of 1/k from k = N l -i particles taking place in the list after particle i.

𝑗 = (𝑖 + 1) + 𝑖𝑛𝑡(𝑘 × 𝑅 𝑛𝑓1 ) ( 50
)
Each pair is then checked for collision with the probability, which with taking into account factor 𝐹 𝑛𝑢𝑚 and time step ∆𝑡 reads

𝑊 𝑖𝑗 = (𝑁 𝑙 -𝑖 )𝐹 𝑛𝑢𝑚 ∆𝑡𝜎 𝑖𝑗 𝑔 𝑖𝑗 𝑉 𝑙 > 𝑅 𝑛𝑓2 . ( 51
)
It should be noted that the ∆𝑡 should be adjusted so that in great amount 𝑊 𝑖𝑗 does not exceed unity , say

𝑝𝑟𝑜𝑏{𝑊 𝑖𝑗 ≥ 1} → 0 (52)
It is worth noting that the probability for 𝑊 𝑖𝑗 > 1 should be kept always close to zero by choosing appropriate time step and cell size. The SBT procedure avoids the production of at least part of the eventually successively repeated collisions which occurs in the NTC scheme when it is applied with a small number of particles. More details of the theoretical background of this scheme are available in Refs. [START_REF] Stefanov | Particle Monte Carlo algorithms with small number of particles in grid cells[END_REF][START_REF] Stefanov | On DSMC calculations of rarefied gas flows with small number of particles in cells[END_REF]. The sequence of SBT collision procedure is shown in Fig. 12.

Fig. 12 SBT collision procedure per each cell

One of the benefits of SBT appears if one decreases cell size (dx) and time step (dt) without increasing the number of particles. Thus, this refines the simulation process when approaching closely the real process of molecular dynamic collision. The difference is only that SBT process is stochastic (random choice of the collision parameters). Simply, this process is localization of collisions more precisely with O(dx,dt) error. This process is not possible in NTC and NN methods because these algorithms always calculate a maximum number of collisions based on formula that estimate a gas in equilibrium and allows an estimation requiring a larger number of collisions. The situation with MFS scheme is better than NTC and NN, as there are not reminders in the estimation of real number of collisions with an integer, but all of these schemes assume a

Poisson distribution of the collision pairs to be checked for collision. However, this idea does not work well when there are small number of particles because these schemes must to define preliminary a number of pairs (random or not) to be checked. Here it is possible the situation that the Poisson distribution or N coll formula gives 5 pairs to be checked even when there are only 2 particles in the cell. Or Poisson distribution predicts one pair to be checked when there are 5 particles in the cell. This happens very often when there are small mean numbers of particle in the •checking the collision probability from Eq. ( 51) velocity update

•change particle velocities cell. This gives considerable fluctuations in the collision process which leads to considerable errors in results.

A note on repeated collisions in SBT and NTC

The SBT procedure, as described in Sec. 3.7, prevents the possibility of repeated collisions in each time step. However, for possible repeated collisions in successive time steps that can occur in SBT there are two possible, one mathematical and other physical, answers. The mathematical answer is based on the fact that the SBT algorithm was derived from the Kac collision model assuming a strict Markov collision process within a time step [START_REF] Amiri | DSMC Simulation of Micro/Nano Flows using SBT-TAS Technique[END_REF]. The collision is based on its initial state, i.e., set of particle velocities in a cell at the beginning of each time step and the transition collision probabilities determined for each particle pair in a cell. After completing the collisions within a time step, a new state, i.e. velocity distribution is formed that is employed in the initial state for the collisions within the next step. At this moment, the step of particle free motion that changes the set of particles in the considered cell is neglected. It is obvious that the past of the collision history in a strict Markovian process is neglected. That is, the error realized by repeated collisions is included in the error of the Markov approximation of the collision process. Practically, the second answer is that first the exchange of particles between cells during the particle free motion occurs and then the probability of a collision of a given pair depends on the time stepthe smaller time step leads to smaller probability for collision. The probability for a repeated collision is proportional to the product of pair probabilities. Therefore, it can be shown that the order of the events is O(t 2 ). Since the accuracy of the splitting scheme is O(t), one can neglect the repeated collisions as introducing an error of higher order for very small time steps.

Altogether with the particle free motion step this allows to neglect the effect of the repeated collisions realized in successive time steps.

In the NTC and NN schemes, on the other hand, the situation is different, i.e., the decrease of time step leads to decrease of the maximum collision number in a collision cell within a time step.

However, for a given particle pair the collision probability is determined by the acceptancerejection rule with the probability given by Eq. ( 10) that does not depend on time step. Thus, the probability for a repeated collision of a twice chosen pair in the NTC algorithm does not decrease with the decrease of time step. That is why; NTC codes must introduce a preventing repeated collision rule within more time steps. However, the direct prevention of repeated collisions introduces extra errors in all schemes using the maximum collision frequency approach (NTC, MFS, Null-collision) when the number of particles per each cell is small. These schemes are based on the combination of acceptance-rejections procedure and the random choice of the particle pairs. When the number of particles is small, to complete the time step the collision procedure may check many times the same pair that previously was chosen and rejected. Next time the routine can choose the same pair and check it until a collision was realized. Using the direct prevention procedure, a conditional probability with removing from consideration some pairs destroys the correct random choice of pairs and the non-uniformity of probability distribution becomes significant. This leads to deviation from the correct probability, given in the Kac stochastic model for collision of each particle pair. This explains while the NTC, NN and MFS schemes cannot work accurately with small number of particles in cells. This point will be supported with the results presented in the following sections.

SBT on dual grids [24]

One reason for errors in the DSMC method stems from particles that are at a collision separation distance and have the potential to be selected as a collision pair, but they are not checked for collision as they are located in two neighboring cells. This difficulty may be overcome by using a dual grid (staggered grid) that is created by translation of cells (or subcells). This strategy improves the collision accuracy specially when there are a few particles per cell. For an orthogonal grid, displaced grids are created by translating cells in each axis direction by a distance equal to a half-cell size. Figure 13 shows the formation of the displaced grid in structured cells. In this figure, solid lines show the main grid and dashed lines show the displaced grid that is obtained by displacing the main gird by a magnitude of Δ𝑥/2 in axial direction and Δ𝑦/2 in normal direction. It should be noted that the number of cells in the displaced grid is (𝑁𝐶𝑋 + 1)(𝑁𝐶𝑌 + 1), where 𝑁𝐶𝑋 and 𝑁𝐶𝑌 are the number of cells in the original gird in 𝑥 and 𝑦 directions, respectively. The volume of internal cells in the displaced grid is equal to the volume of internal cells in the main grid. Nevertheless, the volume of boundary cells must be adjusted. Figure 13 also depicts two particles that in spite of being so close, they are not selected for collision as they are in two different cells in the main grid. After the grid is being displaced, these particles belong to one cell and possibility of collision will be checked. It should be noted that the numerical results obtained by using the dual staggered grid showed a slight improvement with regard to the space resolution compared to a single grid SBT application. (2) The second group is those which are in flow but they are vacant, thus they are assumed inactive in the collision occurrence.

(3) The third group of subcells is the ones that are totally or partially in the flow and at the same time contain at least one particle. These subcells called as "in-flow active" subcells. They are the subcells that SBT could use them as the areas for performing a collision. In the probability function of the SBT scheme, Equation ( 51), V 𝑙 is the volume that is assigned to the particles.

Since the collision cell is considered as a homogeneous space, Goshayeshi et al. [START_REF] Goshayeshi | DSMC Simulation of Hypersonic Flows Using an Improved SBT-TAS Technique[END_REF] suggested that the cell volume is equally distributed among only in-flow active subcells (type 3). In regards of this definition, we can modify the probability function in the form of the Equation [START_REF] Ivanov | Direct statistical simulation method and master kinetic equation[END_REF]. For the l th subcell, it will be:

W ij 𝑙 = F num (N 𝑙 -i) (σg) ij ∆t V 𝑙 V 𝑙 = V c NS type 3 (53) (54)
, where w ij l is the collision probability of the pair (i, j) in the l th subcell, N l is the number of indexed particles for this subcell, NS type 3 is the number of type 3 subcells, and V l is calculated in a way that all the collision volume (V c ) is equally distributed among type 3 subcells.

Corresponding to this volume propagation, Fig 14-b demonstrates how the cell volume is heterogeneously divided among the type 3 subcells. Eq. ( 54) could be substituted the direct calculation of the subcell volume, which could be relatively time consuming. This equation preserves a consistency with probabilistic nature of the DSMC method and is not limited to boundary cells, i.e., it could be utilized for all subcells in the simulation domain.

SBT-TAS: Time step interval control [30]

Previous investigation showed that there is a requirement of time step interval control at the subcell level [START_REF] Goshayeshi | DSMC Simulation of Hypersonic Flows Using an Improved SBT-TAS Technique[END_REF]. There are some restricting rules on time step in the DSMC method, in both of motion and collision stages of particles. In this sense, any manipulation that violates these rules will result in errors. The SBT's features require that the TAS implementation to depend on the local parameters of the subcells. In fact, SBT cannot accept dependent TAS grids similar to NN or NTC does. This incapability comes from the nature of SBT, which force the scheme to determine the collision-frequency by investigating through binary considerations of particles, unlike other schemes that fulfill this task with one step estimation of the total number of collision events. Consequently, implementing TAS in SBT will produce the so called "subcells" as smaller collision cells, while the time step interval employed in these subcells are used from main cells. This time step interval is calculated such that it produces a consistency between the allowed distances for the movement (transit time) and the allowed collision-frequency based on the mean collision time which has been determined in the main cell. Of course, there is no necessity that the mean collision time of the main cell and each of its subcells be the same. This phenomenon will be intensified if the considered cell is placed in complex regions such as shock interactions. In these situations, implementing SBT-TAS using collision cell time step will lead to a time-error, and the violations of DSMC's movement and collision decoupling will occur.

In this state, those subcells that have a time step smaller than their actual value (S 2 -S 5 subcells), will under-predict number of collisions. On the other hand, subcells that have a time-step larger than the actual value, subcell S 1 , over-predict the number of collisions, while particle-transition paths are larger than the allowed amount. Consequently, violations of movement and collision decoupling will emerge in these subsets, see Fig. 16. , is a constant that almost preserves its magnitude, even when the time-step of the cell (𝑇 𝑀-𝐶 ) is being controlled in such a way that it becomes equal to the available smallest mean collision time among concurrent subcells of the collision cell. Therefore, any attempt to avoid this error by reducing the working time-step will be unsuitable and it will just decrease the speed of the numerical procedure, while the 𝑆 𝑣 still exists in the cell space as a time-scaling factor.

In the SBT scheme, the acceptance-rejection depends on the value of the probability function (𝑤 𝑖𝑗 𝑙 ) of the two colliding particle pairs (𝑖, 𝑗) in the 𝑙 𝑡ℎ subcell. It is observed that the cells that are probable of containing S 1 subcells are mostly placed in the shock or shock-shock interaction regions. The sensible characteristic of these cells is that they have a heterogeneous distribution of particles within their space, which is a key factor in generating S 1 subcells. With increasing the usage of the TAS grid on these cells, the number of S 1 subcells will also increase. Hence, the TAS grid should be fabricated over these cells on such an extent that there would be no S 1 subcell. We assume that 𝑊 𝑖𝑗 𝑙 > 1 is a sign of the presence of S 1 subcells. Therefore, the controlling variable for each collision cell is defined as this:

𝐶. 𝑉 = 𝑁 𝐸 𝑁 𝐶 (55)
, where N E is the times that w ij l > 1 occurs and N C is the number of collision occurrences. In fact, it is a ratio between the times that the probability function (w ij l ) exceeds the acceptable limit to the all times that collision happens in that area, and at this work, it is also called as P-exceed ratio or P_E_R. Therefore, the implemented solution suggested in Ref. [START_REF] Goshayeshi | DSMC Simulation of Hypersonic Flows Using an Improved SBT-TAS Technique[END_REF] is to enable the system to sense the degree of the generated error and secure it to a safe zone, where the threat of creation S 1 subcells is at the minimum state. At the same time, this solution should let the TAS system to work with its maximum capacity for the SBT scheme. This combination of SBT-TAS with this control mechanism is called as, "SBT-TAS C.M.".

Intelligent SBT (ISBT)

The SBT collision scheme is designed to obtain correct collision frequencies with the minimum number of particles per cell. A modified version of this scheme, called intelligent SBT (ISBT) has a semi-perception of inter-particle distances and is able to prioritize the acceptance of closer pairs, was suggested by Goshayeshi et al. [START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF]. The aim of developing of ISBT is to alter the particle sorting and acceptance-rejection process in the SBT scheme to bias the selections in support of closer partners. ISBT has two main stages. The aim of the first step is to define a smart way of locally indexing particles. The particle indexing should be modified in a way that a semiperception of distance appears in the indexing process, i.e. in a collision cell with 𝑁 𝑐 particles, it should be an indexing method that satisfies the condition of (𝛿 𝑖 < 𝛿 𝑗 , 0 < 𝑖 < 𝑗 ≤ 𝑁 𝑐 ), where 𝛿 𝑖 is the distance of the 𝑖 𝑡ℎ particle to a reference point ® (see Fig. 17-a). Investigations over the optimized location of the ® point exhibits that the indexed points, relative to their hierarchical order, would be closer to each other if the ® point is located on the corners of a rectangle that surrounds the particles (compare Figs. 17 (b-c)). Therefore, Goshayeshi et al. [START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF] suggested the random selection of cell/subcell corners to efficiently index particles quite close to each other.

However, since the 𝛿 distance does not consider the angular distances, there are still some particles that are indexed successively while they are far from each other, specifically those particles on the opposite sides of the diagonal (e.g. particles 4 and 5 in the Fig. 17-c). Statistically, these particles are not the majority and, on the other hand, any other way for particle indexing would unavoidably produce situations that some indexed particles would be far from each other.

It is crucial to note that while along the indexing direction (diagonal RB in Fig. 17-a) particles are prevented to collide in opposite corners (because they are indexed at the bottom and the top of the hierarchical order), along the other one (diagonal AC in Fig. 17-a) particles might be indexed as neighbors, and their collision, reverse the angular momentum associated with that collision pairs. Since these kinds of collisions might be inevitably supported in the ISBT scheme, as a solution, Goshayeshi et al. [START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF] suggested using elongated cell/subcells. In a quadrangular cell/subcell ('RABC' in Fig. 17-d) the maximum distance (MD) that two far-indexed particles could have is √2 times greater than the side length. However, using elongated cell/subcell obliges the time step to be reduced accordingly. Therefore, in cases of hypersonic shock interaction simulations, Goshayeshi et al. [START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF] suggested using elongated cells/subcells with elongation factor of two (𝑚 = 2), while the time step is reduced about one third. Therefore, ISBT technique could be implemented alongside with the transient adaptive sub-cell technique (TAS) as described in Sec.

M D Q M D R 1 M D R 2 L/m A B E F C D , 2 /  L MD Q m L MD R / 2 / 1  1 2 / 4 2 2    m m L MD R
for m>1, MD R2 <MD R1 <MD Q Fig. 18 ISBT collision procedure • choosing the other particle from Eq. ( 50).

if accepted

•checking the collision probability from Eq. ( 51) using the same random number used in Eq. ( 50) velocity update

•change particle velocities if collision is accepted

Fig. 20 SBT-TAS results for the equilibrium collision frequency ratio in the cavity using a 25×25 grid with 16×16 subcells per each cell and PPSC=2 [START_REF] Amiri | DSMC Simulation of Micro/Nano Flows using SBT-TAS Technique[END_REF].

The collision frequency of the ISBT scheme in an equilibrium state, for the argon flow in the cavity test case with U lid = 0 m/s, Kn=0.01 and with surface and fluid temperatures have the same value of 273 K is evaluated using DS2V code [START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF]. Fig. 21-a, demonstrates the ability of the ISBT scheme to correctly predict the number of collisions, and Fig. 21-b, states that the deviation of the 𝐶𝐹 𝑟𝑎𝑡𝑖𝑜 from the unity in the ISBT collision scheme, similar to the SBT scheme, is bounded with acceptable limit of 0.002. As the simulated cavity has diffuse reflecting walls, the eventual correlations cannot propagate at long distances The dimensionless number, 𝑆𝑂𝐹, which is the mean collision separation distance (𝑀𝐶𝑆) divided by the local mean free path (λ), is selected as a parameter for measuring the quality of collisions.

Fourier flow: comparison with theory

The convergence behavior and accuracy of the DSMC technique for Fourier problems for nearcontinuum conditions was investigated for the NTC method [START_REF] Gallis | Molecular gas dynamics observations of Chapman-Enskog behavior and departures therefrom in nonequilibrium gases[END_REF][START_REF] Rader | Direct simulation Monte Carlo convergence behavior of the hard-sphere-gas thermal conductivity for Fourier heat flow[END_REF]. A similar evaluation was made for the SBT method [START_REF] Taheri | Evaluation of the SBT collision model for near-continuum nano Fourier flows[END_REF]. The hard-sphere argon was simulated at the reference pressure for the conductivity ratio approach to unity as cell size reduces.

Cavity flow

The SBT solution was compared with the NTC and MFS methods for the cavity geometry at Kn=0.005 and a moving lid at the velocity of 𝑈 𝑙𝑖𝑑 = 100 𝑚/𝑠 and 10 𝑚/𝑠 [START_REF] Amiri | DSMC Simulation of Low Knudsen Micro/Nano Flows using Small Number of Particles per Cells[END_REF]. Cavity flow is considered as a classical test case for DSMC solvers [START_REF] John | Investigation of heat and mass transfer in a lid-driven cavity under nonequilibrium flow conditions[END_REF][START_REF] John | Effects of incomplete surface accommodation on nonequilibrium heat transfer in cavity flow: A parallel DSMC study[END_REF][START_REF] Akhlaghi | A New iterative wall heat flux specifying technique in DSMC for heating/cooling simulations of MEMS/NEMS[END_REF][START_REF] Mohammadzadeh | Thermal and second-law analysis of a micro-or nanocavity using direct-simulation Monte Carlo[END_REF][START_REF] Mohammadzadeh | A parallel DSMC investigation of monatomic/diatomic gas flows in micro/nano cavity[END_REF][START_REF] Akhlaghi | A Novel Algorithm for Implementing Specified Wall Hat Flux in DSMC: Application to Micro/Nano Flows and Hypersonic Flows[END_REF][START_REF] Mohammadzadeh | DSMC and R13 modeling of the adiabatic surface[END_REF]. Variable hard sphere (VHS) model of the monatomic argon, m=6.63×10 -26 Kg and d=4.17×10 -10 m was considered as the gaseous flow. Using a 400 × 400 grid cells, the minimum required particle per cell was investigated for the NTC scheme. As Fig. 23 depicts temperature on the centerline of the cavity, < 𝑁 > = 20 (< > means average) could reasonably provide the sufficient accuracy for the NTC scheme. Fig. 24 compares temperature in the vertical centerline from the NTC with the SBT and MFS solutions.

The figure indicates that the only scheme which provide accurate results with N=2 is the SBT scheme. For the case of Kn=0.005, T w =300 K and U lid = 100 m/s, 25×25 grid, PPC=512 and PPSC=2, solution of the SBT on dual grid (SBT-Dual) is compared with the SBT-TAS and NTC solutions, the latter here called correct solution and is obtained on a 200×200 grid with 2 fixed subcells in each direction and 20 particles per cell [START_REF] Amiri | DSMC Simulation of Micro/Nano Flows using SBT-TAS Technique[END_REF]. Since the employed grid is coarse, it is quite expected that if we do not employ TAS, incorrect SBT solution will be attained, as is shown in Fig. 25, which depicts velocity components and temperature profile on the horizontal axis of the cavity. The figure shows that using dual grid does not improve the accuracy of the SBT-TAS solution more. Fig. 26 compares SOF from the SBT-Dual, SBT-TAS and SBT-TAS-Dual. The figures shows that using TAS drastically reduced the SOF, and use of dual grid reduced the SOF by half, however, had not any impact of the accuracy because the respective SOF was almost below 0.3, the minimum value that Bird stated could guarantee accuracy of solutions [START_REF] Bird | The DSMC Method[END_REF].

Flat plate flow

Hypersonic flow over a small scale plat of length 90 nm with the surface temperature of 𝑇 𝑤𝑎𝑙𝑙 = 500𝐾, Kn= 0.01 and T in =300 K was considered as a test case for the SBT, NTC, MFS and SBT-TAS schemes [START_REF] Amiri | DSMC Simulation of Low Knudsen Micro/Nano Flows using Small Number of Particles per Cells[END_REF]. The effect of number of particles per each cell is depicted in Fig. 27, where the figure indicates that both of NTC and MFS schemes require around 20 particles per cell to predict particle independent results while the pure SBT scheme works well even with one particle per cell in this case, i.e., sensitivity of the SBT solution to the PPC is quite low.

Rarefied flow past an airfoil

Rarefied supersonic and subsonic gas flow around a small-scale NACA 0012 airfoil was simulated using SBT and collision scheme [START_REF] Shoja-Sani | Investigation of Rarefied Gas Flow around NACA 0012 Airfoils using DSMC and NS solvers[END_REF][START_REF] Le | Rarefied gas flow simulations of NACA 0012 airfoil and sharp 25-55 biconic subject to high order nonequilibrium boundary conditions[END_REF]. The small-scale airfoils are widely employed in micro-air-vehicles. Flow field is investigated at Mach=2, Kn=0.26 and α=45°. Fig. 30 shows distribution of pressure coefficient and normalized slip velocity on the airfoil surface obtained using different values of PPC using SBT. In this case, the figure well indicate that SBT could provide accurate solution using even PPC=1. In addition, SBT solution shows almost no variations, expect statistical noises, in results with the changes in PPC. Figures 30-c-d

Micro and nano Nozzle flow

Flow inside micro and nano nozzles has been of interest of various researchers, see [START_REF] Darbandi | Study of Supersonic-Subsonic Gas Flows through Micro-Nano Scale Nozzles Using Unstructured DSMC Solver[END_REF][START_REF] Sebastião | Numerical simulation of heat transfer and pressure distributions in micro nozzles with surface discontinuities on the divergent contour[END_REF][START_REF] Sebastião | Impact of surface discontinuities on flowfield structure of a micro nozzle array[END_REF][START_REF] Lijo | Analysis of Supersonic Micro nozzle Flows[END_REF].

Accuracy of the SBT technique in treating rarefied micro and nano scale nozzle was considered [START_REF] Saadati | Detailed Investigation of Flow and Thermal Field in Micro/Nano Nozzle using Simplified Bernouli Trial (SBT) Collision Scheme in DSMC[END_REF]. Nitrogen flow with an inlet Knudsen number of Kn inlet =7.781×10 -4 was considered. Fig. 31 shows that SBT technique using PPC=2 at the steady state condition could provide particle independent results. Fig. 32 compares Mach number contours inside the micronozzle using both of SBT and NTC schemes, where both solutions show a good agreement. It should be noted that NTC needs PPC=10 to provide accurate solutions for all flow properties inside the nozzle. Hypersonic nitrogen flow of Mach 15.6 passing over a 25°-55° degree biconic geometry is a quite complex test case for numerical schemes as it consists of a laminar recirculation zone, laminar expanding zone, high-speed low-density region, and low-speed high-density region [START_REF] Knight | RTO WG 10: Test Cases for CFD Validation of Hypersonic Flight[END_REF][START_REF] Holden | Measurements of Regions of Low Density Laminar Shock Wave/Boundary Layer Interactions in Hypersonic Flows and Comparison with Navier-Stokes Predictions[END_REF][START_REF] Moss | Direct simulation Monte Carlo simulations of Hypersonic flows with shock interactions[END_REF][START_REF] Harvey | A review of a validation exercise on the use of the DSMC method to compute viscous/inviscid interactions in hypersonic flow[END_REF][START_REF] Titov | Implications of Slip Boundary Conditions on Surface Properties in Hypersonic Flows[END_REF]. Fig. 34 shows the biconic and structure of different shock interactions occurred over this geometry.

Re-entry test case with chemical reactions [85]

The SBT-TAS and NN collision models are applied to treat a typical hypersonic atmospheric re-entry problem exposing to dissociation and recombination chain reactions of real-air. The Total Collision Energy (TCE) [START_REF] Bird | Simulation of multi-dimensional and chemically reacting flows[END_REF][START_REF] Gupta | A Review of Reaction Rates and Thermodynamic Transport Properties for an 11-Species Air Model for Chemical and Thermal Non Equilibrium Calculations to 30,000 K[END_REF] model, with reactions listed in Table 3, is used to simulate chemical reactions in both collision schemes of NN and SBT-TAS. Fig. 37-a shows the schematic of the re-entry geometry and also indicates investigated working conditions. Fig. 37-b shows temperature contours from both SBT-TAS and NN scheme, with suitable agreement expect small differences in the wake region. Figs. 37-c-d compares temperature jump and species distribution over the surface, with an almost reasonable agreement between both schemes.

Concluding remarks

The current paper provides a review of the DSMC collision schemes, originating either from the Boltzmann equations (i.e., TC, NTC, NN) or the Kac stochastic model (i.e., MFS, null collision, BT, SBT

  Bogoliubov hierarchy) including the last Boltzmann equation have a probabilistic nature because their main unknown is the probability distribution function. Though Liouville equation is simpler than the Newton equation, it considers collisions of N molecules and its solution is quite challenging. A less expensive description is attained by approximating the flow description using only n-particle distribution functions, which determine the probability to simultaneously find n particles independently of the state of the remaining (N-n) particles. Following the BBGKY hierarchy, a chain of linked equations could be obtained for reduced n-particle distribution functions, where the equation for the single-particle distribution function f(t, 𝑥(𝑙) , 𝑐) corresponds to the Boltzmann equation, which considers only binary collisions. The flowchart of hierarchy of the governing equations is shown in Fig. 1. Boltzmann equation does not describe the motion of particles, as Newton's equation does, but the time-dependent change of the probability distribution function. In fact, the Boltzmann equation is a nonlinear, seven-dimensional, integrodifferential equation governing the spatio-temporal evolution equation of the one-particle velocity distribution function f(t, 𝑥

Fig. 1 :

 1 Fig. 1: Flowchart of hierarchy of governing equations

Fig. 2 :

 2 Fig. 2: DSMC collision models derived from the Boltzmann equation and Kac master equation.

Fig. 4

 4 Fig. 4 NTC collision procedure, repeated N coll time for each DSMC collision cell.

Fig. 6 :

 6 Fig. 6: Schematic of the number of subcells in each cell using TAS, bold lines depict cell borders and narrow lines depict subcell borders [49].

Fig. 7

 7 Fig. 7 Pseudo-Poisson (collision frequency) collision procedure.

Fig. 9

 9 Fig.9BT collision procedure[START_REF] Yanitskiy | Operator approach to Direct Simulation Monte Carlo theory in rarefied gas dynamics[END_REF][START_REF] Stefanov | Particle Monte Carlo algorithms with small number of particles in grid cells[END_REF][START_REF] Stefanov | On DSMC calculations of rarefied gas flows with small number of particles in cells[END_REF]

Fig. 10

 10 Fig. 10 Ballot Box collision procedure per each cell [9].

4 -

 4 Collisions are ended if summation of collision time (sampled from the Poisson distribution) exceeds the cell time step. > 𝑡) = 𝜐 𝑚𝑎𝑥 𝑒 -𝜐 𝑚𝑎𝑥 𝑡 (48) If 𝜐 𝑚𝑎𝑥 =cte, this technique is one of the possible approaches to generate a Poisson random value 𝑠(𝑡). A classical algorithm for producing a random number from a Poisson distribution is developed by Knuth [56], as given in appendix 1.

Fig. 11

 11 Fig. 11 Null collision procedure per each cell.

•

  Compute maximum frequncy (Eq. 46) and sample time step from the Poisson distribution (For example, see appendix 1)

  cross-referencing of particles in cell l to number particles from 1 to N l first particle selection •choosing first particle in sequense from list: i = 1, ... ,N l -1 second particle selection • choosing the other particle from Eq. (50).if accepted

Fig. 13

 13 Fig.[START_REF] Bird | The DSMC Method[END_REF] Schematic of the displaced grid (dashed line): solid lines are the main grid, grey cell is a typical cell created after grid displacement. (copied with permission from Ref.[START_REF] Amiri | DSMC Simulation of Low Knudsen Micro/Nano Flows using Small Number of Particles per Cells[END_REF])

Fig. 14 :

 14 Fig. 14: Different subcell types produced around a circular cylinder (grey region); (a) rectangular subcell divisions for indexing particles, (b) subcell divisions based on the assigned volume in the estimation of probability function (divisions are tried to represent an equivalent volume of V 𝑙 ) .

Fig. 15 depictsFig. 15 :

 1515 Fig.15depicts this probable time variation in a series of subcells. In SBT-TAS implementation,

Fig. 16 :

 16 Fig. 16: S 1 subcell and time-step violation.

3. 9 .Fig. 17 :

 917 Fig. 17:Distribution of some particles in a cell/subcell; (a) a system of particle indexing based on the distance of particles to a definite center point ®, (b) the indexing results for a group of nine particles when the ® is located on the center of the cell, (c) when this point is on the left-down corner of the cell, and (d) the geometrical comparison between one quadrangle cell/subcell ("RABC") and its equal-area elongated cell/subcell ("RDEF"), while the calculated values for the maximum distance (MD) between two farindexed particles are shown; they denote that elongated cell/subcell has smaller MD value.

Fig. 19 :

 19 Fig. 19: Pseudo-circular subcells in the ISBT scheme that causes a biased selection for closer collision pairs. (Reprinted with permission from B. Goshayeshi, E. Roohi, and S. Stefanov, Phys. Fluids 27, 107104 (2015). Copyright 2015 American Institute of Physics)

Fig. 21 :

 21 Fig. 21: (a) contour of CF ratio in the cavity at the equilibrium state, using the 100×100 grid and 10 particles per cell, calculated by the ISBT scheme, and (b) CF ratio along the horizontal center line calculated by the NN, ISBT and SBT schemes. (Reprinted with permission from B. Goshayeshi, E. Roohi, and S. Stefanov, Phys. Fluids 27, 107104 (2015). Copyright 2015 American Institute of Physics)

Fig. 22 :

 22 Fig. 22: (a) Comparison of SBT and NTC solution methods for Fourier flow.

Fig. 22 -

 22 Fig. 22-c compares dependency of the SBT and NTC solutions for the wall heat flux on the

Fig. 23 Fig. 24

 2324 Fig.[START_REF] Stefanov | Particle Monte Carlo algorithms with small number of particles in grid cells[END_REF] Effect of number of particle per cell on the accuracy of the NTC collision scheme in prediction of temperature field in the horizontal centerline of the cavity. (reproduced with permission from[START_REF] Amiri | DSMC Simulation of Low Knudsen Micro/Nano Flows using Small Number of Particles per Cells[END_REF])

Fig. 25

 25 Fig.[START_REF] Goshayeshi | A Novel Bernoulli Trials Collision Scheme in the DSMC Method with Intelligence over Particle Distances[END_REF] Results of the TAS technique on the horizontal axis of the cavity obtained using 25×25 grid-PPC=512, PPSC=0.5[START_REF] Amiri | DSMC Simulation of Micro/Nano Flows using SBT-TAS Technique[END_REF]

Fig. 26

 26 Fig. 26 SOF (MCS/) in the cavity flow (a) SBT without applying TAS (b) SBT-TAS (c) SBT-TAS-Dual [49]

Fig. 27 Fig. 28 2 Fig 29 .

 2728229 Fig.[START_REF] Gallis | Effect of collision-partner selection schemes on the accuracy and efficiency of the direct simulation Monte Carlo method[END_REF] Effect of average number of particles per cell on the accuracy of the temperature profile at y/L=0.2 for the flat plate case. Profiles are shifted to right for MFS and SBT schemes to make differences more clear. (reproduced with permission from[START_REF] Amiri | DSMC Simulation of Low Knudsen Micro/Nano Flows using Small Number of Particles per Cells[END_REF])

Fig 30 .

 30 Fig 30. Flow field investigation around an airfoil using SBT collision technique [68].

Fig. 31

 31 Fig. 31 Particle dependency of the SBT scheme for (a) Mach number, (b) normalized temperature at the nozzle exit [76].

Fig. 32 Fig. 33 :

 3233 Fig. 32 Comparison of SBT and NTC techniques for Mach number contours inside the micronozzle [76].4.7 Hypersonic flow over cylinder and biconic geometry

Fig. 34 :Fig. 35 :Fig. 36 :

 343536 Fig. 34: Biconic geometry (left) and shock structures (right) around the geometry obtained by DSMC (Reprinted with permission from B. Goshayeshi, E. Roohi, and S. Stefanov, Phys. Fluids 27, 107104 (2015). Copyright 2015 American Institute of Physics)

Fig. 37 :

 37 Fig. 37: Comparison of SBT-TAS and NN solutions for re-entry flow with chemical reactions [85].

  𝑒𝑥𝑝[𝑡𝜐 𝑐 (𝑇 𝑖𝑗 -𝐼)] = exp (-𝑡𝜐 𝑐 ) ∑

			∞ 𝑠=0	(𝑡𝜐 𝑐 ) 𝑠 𝑠!	𝑇 𝑠
	where: 𝑇 = ∑ 𝑊 𝑖𝑗 1≤𝑖<𝑗≤𝑁 𝑙	𝑇 𝑖𝑗 ;	𝑊 𝑖𝑗 =	𝑤 𝑖𝑗 𝜐 𝑐

Table 1

 1 states that the ISBT, compared with the SBT collision scheme, reduced the 𝑆𝑂𝐹 value quite considerably.

	RATIO CF Ratio
	1.002
	1.0015
	1.001
	1.0005
	1
	0.9995
	0.999
	0.9985
	0.998

Table 1 :

 1 Comparison of the 𝑆𝑂𝐹 = (

		𝑀𝐶𝑆 𝜆	) value in the lid-driven cavity flow between NN, ISBT and SBT
	collision schemes. (Reprinted with permission from B. Goshayeshi, E. Roohi, and S. Stefanov, Phys. Fluids
		27, 107104 (2015). Copyright 2015 American Institute of Physics)
	Collision scheme	Division	𝑆𝑂𝐹 value	𝑆𝑂𝐹 normalized by the NN	𝑆𝑂𝐹 normalized by the
		size in		value	SBT value
		DS2V code
	NN	100-100	0.0477	1	0.182
	ISBT	100-100	0.1740	3.66	0.666
	SBT	100-100	0.2621	5.49	1

Table 2 .

 2 Relative computational time for cavity flow. (reproduced with permission from[START_REF] Amiri | DSMC Simulation of Low Knudsen Micro/Nano Flows using Small Number of Particles per Cells[END_REF])

		<N>	time
	SBT	2	0.79
	NTC	250	1.15
		125	1.07
		20	1
		5	0.98
		2	0.75

Table 3 :

 3 The employed chemical reactions in the investigated test case.

	No.	Reactions
	Dissociation/Recombination reactions
	1	O 2 +N↔O+O+N
	2	O 2 +NO↔O+O+NO
	3	O 2 +N 2 ↔O+O+N 2
	4	O 2 +O 2 ↔O+O+O 2
	5	O 2 +O↔O+O+O
	6	N 2 +O↔N+N+O
	7	N 2 +O 2 ↔N+N+O 2
	8	N 2 +NO↔N+N+NO
	9	N 2 +N 2 ↔N+N+N 2
	10	N 2 +N↔N+N+N
	11	NO+N 2 ↔N+O+N 2
	12	NO+O 2 ↔N+O+O 2
	13	NO+NO↔N+O+NO
	14	NO+O↔N+O+O
	15	NO+N↔N+O+N
		Exchange reactions
	16	NO+O↔O2+N
	17	N2+O↔NO+N

Acknowledgements

Acknowledgements

The first author would like to acknowledge financial supports of Ferdowsi University of Mashhad under grant No. 39035. The second author would like to acknowledge the financial support provided by the EU network program H2020 under grant MIGRATE No. 643095.

SBT on transient adaptive subcells (SBT-TAS) [49]

Previous works [START_REF] Stefanov | On DSMC calculations of rarefied gas flows with small number of particles in cells[END_REF] showed that computational efficiency of the SBT method decreases with increasing the number of particles in a collision cell. Using TAS technique, this problem may be avoided because the transient adaptive subcells are set dynamically in such a way that the number of particles always remains limited in every subcell. During the simulation, the particle distribution in cells is not uniform; therefore, there may be subcells which are empty or contain only one particle. In this case, no collisions are possible in such cells. In the NTC-TAS approach suggested by Su et al. [START_REF] Su | Implementation of a transient adaptive sub-cell module for the parallel-DSMC code using unstructured grids[END_REF], if such a case occurs, usually a search among the neighboring subcells is performed to find a partner for collision pair. It could be shown that there is no need to follow this procedure in the SBT-TAS algorithm. Instead, it is suggested to use a two-half-time-step collision algorithm and shift the transient adaptive subcell grid in the second half-time step in order to prevent the reduction of collision frequency in the cell. By using a dual subcell grid, the single particles in cells within the first half-time step will have the chance to find a collision partner within the second half-time step.

SBT-TAS: volume estimation on curved boundaries [30]

Unlike the NTC collision scheme, the collision probability function of the SBT scheme, Eq. 51, depends on subcell volumes, therefore, there is a necessity of having an accurate volume computation (approximation) of subcells in the SBT scheme. TAS divides the collision cell into equally spaced rectangular spaces, and the SBT procedure can independently operate within their space. By producing subcell divisions over a generally shaped collision cell, there are three types of subcell (Fig. 14-a):

(1) The first group is those which are out of the borders of the collision cell and they are not considered in the flow. The ISBT procedure could be implemented as follows:

The first particle, i.e. 𝑖, is sequentially selected from 1 to 𝑁 𝑐 , while its pair, i.e. 𝑗, at each sequent is randomly selected among the further indexed particles, according to the previous paragraph and following Eq. ( 50). The main difference with the SBT is that instead of regenerating the second random number (𝑅 𝑛𝑓2 ), which was used in the acceptance-rejection section of the SBT scheme in Eq. ( 51), ISBT does not produce any new random number and implements the first generated one (𝑅 𝑛𝑓1 ):

, where 𝑊 𝑖𝑗 is calculated through the Eq. ( 51). Figure 18 show the sequences of the ISBT procedures. As a result of this policy, pseudo-circular subcells (Fig. 19) will be created that allow an easier acceptance for closer available-pairs and a more difficult condition for further particles;

it is because that smaller number for 𝑗, as the second particle in the hierarchical list of indexed numbers, is only obtained if the 𝑅 𝑛𝑓1 is closer to zero, and conversely, greater number for 𝑗, as the second particle, is merely obtained if the 𝑅 𝑛𝑓1 is further from zero.

Validations of the SBT collision family

After presenting different collision algorithms in the framework of the DSMC method, this section presents a brief review of results obtained with the SBT collision family emphasizing the accuracy and suitability of this collision model for a wide spectrum of rarefied flow fields.

Collision frequency ratio from SBT-TAS and ISBT

The equilibrium collision frequency ratio of a cavity flow at 𝐾𝑛 = 0.005, T w =300 K, simulated on a 25×25 grid using SBT-TAS technique, is shown in Fig. 20 [START_REF] Amiri | DSMC Simulation of Micro/Nano Flows using SBT-TAS Technique[END_REF]. In this state, the theoretical equilibrium collision rate per molecule (𝐶𝐹 𝑡ℎ ) is given by [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF]:

, where 𝑛, 𝑑, 𝐾 𝐵 , 𝑇 𝑟𝑒𝑓 , 𝑚 𝑠 , and are number density, gas molecular diameter, Boltzmann constant, reference temperature, molecular mass, and viscosity-temperature exponent, respectively. 𝐶𝐹 𝑛𝑢𝑚 represents the numerical value of this theoretical equilibrium collision rate that is calculated by the division of the number of collisions in each cell (𝑁 𝑐𝑜𝑙𝑙 ) on the execution time (𝑇𝑖𝑚𝑒) and half of the mean particle numbers per cell (0.5 N p ):

𝐶𝐹 𝑟𝑎𝑡𝑖𝑜 is the ratio of the 𝐶𝐹 𝑛𝑢𝑚 to its theoretical value-given by Eq. ( 57)-it must have a magnitude close to unity. Using TAS and setting 2 particles per subcell, i.e., PPSC=2, 16×16 subcells were employed in each cell in average. As Fig. 20 shows, mean deviation of the equilibrium collision frequency ratio of one is in order of 10 -4 and its maximum is about 0.2%. This proves that number of collisions in SBT-TAS scheme is accurate and coincides with the theory even using small number of particles per subcells (collision cells), i.e. two particles. families). Collision models derived from the Boltzmann equation are derived according to the principles of the classic kinetic theory, i.e., considering the concept of a collision cylinder and the maximum number of collisions that could occur in the ratio of the collision cylinder volume to the cell volume. These models needs around 10-20 particles per cell and usually suffer from eventually repeated collisions. On other side, the straightforward attempt to prevent repeated collision makes the collision frequency incorrect for small number of particle in cell. The models based on the Kac stochastic model are mainly derived from it by using mathematical operations, i.e., simplification and re-arrangement of the transition operator (G(t)) that appears in the Kac stochastic equation. These approximations keep the basic properties of the Kac model within small enough time step at arbitrary mean number of particles in cells when the time step and cell size are coupled in an appropriate manner. The first Kac-based models, such as the BT and Pseudo-Poisson (collision frequency) models, suffered from high computational costs, while a newly developed collision class, called the Simplified Bernoulli trial (SBT), has comparable computational costs to the standard models based on the Boltzmann equation such as NTC. A wide spectrum of rarefied flow test cases, ranging from simple equilibrium collision frequency to complicated shock-shock interactions and chemically reacting flows, was simulated by using the SBT collision family models and it was shown that SBT collision family could provide accurate solutions using small particles per cell. The benefit of the SBT collision model will be more clearly pronounced in complex 3D simulations, where a large number of particles is required if one relies on the standard NTC scheme. A hybrid NTC-SBT scheme could be also utilized for these complex 3D cases [START_REF] Stefanov | Rayleigh-Bénard flow of a rarefied gas and its attractors. I. Convection regime[END_REF][START_REF] Stefanov | Rayleigh-Benard flow of a rarefied gas and its attractors. Part 2: Chaotic and periodic convective regimes[END_REF].

Appendix