{'fraction_non_alphanumeric': 0.041163556531284305, 'fraction_numerical': 0.01181851445298207, 'mean_word_length': 4.643402849473467, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exponential BGK Integrator for Multi-Scale Particle-Based Kinetic Simulations

Despite the development of an extensive toolbox of multi-scale rarefied flow simulators, such simulations remain challenging due to the significant disparity of collisional and macroscopic spatio-temporal scales. Our study offers a novel and consistent numerical scheme for a coupled treatment of particles advection and collision governed by the BGK evolution, honouring positivity of the velocity distribution. Our method shares its framework, in spirit, with the unified gas kinetic class of multi-scale schemes. Yet it provides attractive features for particle-based stochastic simulations, readily implementable to existing direct simulation Monte-Carlo codes. We demonstrate accuracy and performance of the devised scheme for prototypic gas flows, over a wide-range of rarefaction parameters. Due to the resulting robustness and flexibility of the devised exponential BGK integrator, the scheme paves the way towards more affordable simulations of largeand multi-scale rarefied gas phenomena.


I. INTRODUCTION
The direct simulation Monte Carlo (DSMC) method has become the standard approach for simulations of rarefied nonequilibrium gas flows [1]. It is built on the main idea of employing discrete particle collisions in order to statistically mimic the molecular collisions. An important success element of DSMC rests on the particle-based treatment of the phasespace, which allows rather straightforward integration of the inner energy modes, chemical reactions as well as boundary interactions. Despite the established accuracy of DSMC across the whole rarefaction regimes, it comes with a drawback of dense operations near the continuum as the mean free path and the collision frequency have to be resolved. As a consequence, it is still challenging to simulate flows with large variation of the Knudsen number, covering both continuum and rarefied phenomena. To overcome this central issue, different approaches have been proposed. The most direct solution is to couple continuum solvers based on the Navier-Stokes equations with DSMC [2,3]. However since the two solvers operate at different levels of the flow description, i.e. macroscopic variables in the Navier-Stokes and particle-based probabilities in DSMC, various problems arise which limit the generality and robustness of such approaches.
Another possibility is offered by multi-scale methods based on the discretisation of the phase-space, as provided e.g. in Discrete-Velocity-Methods (DVMs) [4]. Along DVM, recently there have been several developments such as UGKS [5,6] and DUKGS [7], where Bhatnagar-Gross-Krook (BGK) approximations of the collision term are computed efficiently. Due to the deterministic treatment of the kinetic problem, these methods provide noise-free solutions, which make them well suited for low-Mach flows. However this comes with the price of the velocity-space discretisation, which might render these methods inefficient for high-Mach non-equilibrium flows with significantly extended velocity domains.
An interesting remedy to circumvent this problem is provided by the UGKWP method [8]. Here the non-equilibrium part of the velocity distribution is represented by particles, to reduce the corresponding cost of the velocity discretisation. Furthermore, such hybrid representation allows for noise-free solutions in the hydrodynamic limit, as particles are only employed in the non-equilibrium portions of the flow. Nevertheless, a mixed particle-DVM treatment of the distribution leads to implementation challenges, especially if the approach is to be integrated into the existing mature DSMC solvers.
At a different front, quite a few particle-based BGK [9][10][11] and Fokker-Planck [12][13][14] (FP) methods came to the fore in recent years, as they could simply be coupled with DSMC while also being efficient at moderate/low Knudsen flows. On the BGK side, especially the particle-based ellipsoidal statistical BGK [15,16] (ES-BGK) and the Shakhov BGK methods were investigated [9,17]. For the FP models, the main work was done on the Entropic FP model (EFP) [13], cubic-FP model [18], the ellipsoidal statistical FP (ESFP) model [14] and the model of Bogomolov [19]. However while the proposed schemes might reduce the cost of dense collisions in the continuum regime, they still need fine resolutions. This is due to the typical first-order treatment of particles evolution which is implied by splitting between free flight and collision/relaxation sub-time-steps.
In this study, we devise a particle-based multi-scale BGK solver. We remain in the realm of pure particle representation of the distribution with identical positive statistical weights. This avoids implementation overheads of integrating the scheme with existing DSMC solvers. Moreover, we anticipate straightforward extensions to more complex non-equilibrium phenomena (e.g. chemical reactions). The multi-scale capability of the scheme is achieved by coupled position-velocity integration of the particles evolution, which is discussed in the following, after brief review of the govern-ing theory.

II. THEORY
The Boltzmann equation describes a monatomic gas flow with the corresponding distribution function f = f (x, v,t) at position x and velocity v where external forces are neglected and δ f /δt| coll follows the Boltzmann collision integral Here, S 2 ⊂ R 3 is the unit sphere, n is the unit vector of the scattered velocities, B is the collision kernel and the superscript denotes post-collision velocities. The multiple integration of this non-linear collision term, besides highdimensionality of the solution domain make the Boltzmann collision integral computationally complex. For this reason, the DSMC method reduces the collision integral to the Monte-Carlo sampling of collision events among random particles representing the gas flow. Either through direct solvers or by random particles, the mean free path as well as the collision frequency must be resolved, leading to a significant computational effort for low Knudsen number flows.

A. BGK Approximation
The BGK model approximates the collision term by a nonlinear relaxation form, where the distribution f relaxes towards a target f t with a certain relaxation frequency ν. It is assumed, in the original BGK model, that the target distribution is the Maxwellian with the number density n, the molecular mass m, temperature T , the Boltzmann constant k B and the thermal velocity c = v − u with the average flow velocity u [20]. The relaxation frequency gives rise to the viscosity The Maxwellian distribution, as the target, leads to the fixed Prandtl number of Pr = µc P /K = 1, whereas the Prandtl number for monatomic gases is close to 2/3 [21]. To overcome this problem, several extensions of the BGK model were introduced in the past. Some of these models transform the target distribution function e.g. the ellipsoidal statistical BGK model [22] or the Shakhov BGK model [23], while others modify the relaxation frequency from a constant to a function of the microscopic velocities as described by Struchtrup [24].

B. Particle based BGK solver
While the BGK relaxation has a much simpler construct compared to the Boltzmann collision integral, still stiff relaxation might be encountered if explicit time integrations are applied. Furthermore, to obtain asymptotical convergence to the Navier-Stokes equations, the advection (particle movement) and the relaxation term should be treated in a coupled way, as discussed by Xu and Huang [5], Guo et al. [7].
In the stochastic particle BGK method (SP-BGK), the time integration of Eq. (3) is typically done by assuming a constant target distribution [15,16] The advection part follows a first-order operator-splitting which means that the particles of f (v, x,t + ∆t) are moved along their trajectory for a full time step ∆t to reach f (v, x + v∆t,t + ∆t). This first-order time integration has some advantages for particle methods. For ν∆t 1, the time integration Eq. (6) recovers the forward Euler method Yet in contrast to the forward Euler Eq. (7), the prefactors of f (v, x,t) and f t are always positive in Eq. (6), also when ν∆t > 1. In the stochastic particle method context, Eq. (6) can be easily realised if each particle within a cell gets a new velocity sampled from f t with the well-posed probability (1 − e −ν∆t ). For the Euler forward method, however, similar realisations only work as long as ν∆t remains below 1. Otherwise, particles with negative weights have to be introduced in order to construct the desired distribution function, as the prefactor of f (v, x,t) becomes negative in Eq. (7). While ν∆t > 1 may not be relevant for the first-order explicit time integration anyhow, the unconditional positivity of the distribution function reveals the advantage of Eq. (6) over Eq. (7) for particle methods, both in terms of the robustness as well as the implementation.
To improve the explicit time integration to a secondorder method coupling advection and relaxation, we follow UGKS [5,6] and DUGKS [7]. The Crank-Nicolson method is used for the time integration Next, let us introduce two additional distributions After inserting the above-introduced distributions into Eq. (8), one sees thatf (v, x + v∆t,t + ∆t) is nothing but the particles off (v, x,t) which are moved along their trajectories for ∆t, i.e.f This approach provides an implicit integration of the coupled advection and relaxation, suitable in theory also for stochastic particle methods. Note that, in the case of DVM (including UGKS and DUGKS), some additional work has to be done for the flux reconstruction. Nevertheless, for stochastic particle methods the problem of negative pre-factors for large ∆t arises again, so that additional particles with partly negative weighting factors would have to be added to construct the distribution functions. A practical solution was provided by Fei et al. [25,26], where an additional collision term was inserted in which the current distribution function is approximated by a Grad-13 approximation. This allows the advection and the relaxation process to be solved together. The choice of the Grad-13 distribution also ensures that the Navier-Stokes limit is asymptotically preserved. Finally, a multi-scale parameter was defined, allowing a smooth transition between the normal SP-BGK method and the one subject to the Grad-13 approximation.
Even though a good performance of the SP-BGK method regularized by Grad-13 approximation has been demonstrated, it is still desirable to construct a rigorous SP-BGK algorithm with overarching properties of the DSMC method, without the need of introducing auxiliary approximations to deal with the negative pre-factors. In the following, we construct such second-order particle method without making additional approximations on the distribution functions. Thus, the introduction of an additional multi-scale parameter is omitted and the construction would come consistent with the BGK equation, in a fashion similar to UGKS and DUGKS.

C. Exponential Time Differencing of the BGK equation
The exponential time differencing has advantages especially in dealing with ordinary differential equations that contain a stiff linear term [27] where c is a constant that leads to the stiff system, and F is the non-linear part. The main idea is to integrate the stiff term in the exact form, and then approximate the nonlinear term numerically (typically explicitly), as following which is exact up to the point where the occurring integral is approximated. If we now look at the relaxation term of the BGK equation in a form corresponding to Eq. (12) we face an additional problem that the non-linear part ν f t is also stiff (since it scales with ν). Following Eq. (13), we get e νs ν f t (t + s) ds. (15) Now in order to treat the stiffness of the non-linear term, an implicit integration, such as the Crank-Nicolson method similar to the DUGKS method, can be carried out. In the following, we carry out this procedure and it is shown that it has various advantages for particle methods, since the prefactors mentioned always remain positive.

Linear approximation
Let f t t+s and f t t represent the numerical approximations of the target distribution at time steps t + s and t, respectively. To construct a second-order scheme, we adopt the following linear approximation for the target distribution leading to a Crank-Nicolson type scheme. Using this approximation, the integration of Eq. (15) yields to where f t+∆t and f t denote the numerical approximations of f (t + ∆t) and f (t), respectively. After a few rearrangements and taking particle movement into account, Eq. (17) leads to the following integration scheme and Equation (18) already contains the simple form of the first order SP-BGK method (6). Furthermore the prefactors A ∈ [0, 0.5] and B ∈ [0.5, 1] are non-negative. Following Xu and Huang [5], Guo et al. [7], two additional distributions can be introducedf whereby once morẽ is honoured. Therefore, by moving the particles off (v, x,t) along the trajectoriesf (v, x + v∆t,t + ∆t), the advection and relaxation processes are coupled and the time integration remains second-order accurate. Starting from f t (v, x),f (v, x,t) can be constructed from Eq. (21). Afterwards,f (v, x + v∆t,t + ∆t) is realised by moving the particles fromf (v, x,t), and finally f t+∆t (v, x + v∆t) is constructed via While all prefactors remain non-negative, a re-normalisation has to be introduced as the prefactors of Eqs. (21) and (24) do not sum up to unity. To avoid deleting and creating particles, the introduced distribution functions are normalized by the factor and To further simplify the time integration, one can also omit the intermediate step via f and only track the additional distribution functionsf * andf * . By substituting Eq. (26) into (27), This is a convenient formulation as similar to the original SP-BGK, particles have to be sampled from the target distribution with the probability (1 − e −ν∆t ) to constructf * fromf * . The two most common target distributions are the ellipsoidal statistical BGK (ES-BGK) and the Shakhov model (S-BGK). The ES-BGK distribution is given as with the matrix Here, δ i j is the Kronecker delta, Pr is the Prandtl number, ρ is the mass density and p i j is the trace-less pressure tensor In the case of the S-BGK model, the target distribution is defined as with the heat flux vector Therefore to evaluate Eq. (28), either knowledge about p i j ( f ) or q( f ) become necessary. Due to the normalization of the distributions, as already extensively discussed by Xu and Huang [5], Guo et al. [7], the mass, momentum and energy are conserved by the collision operator and can be thus determined directly from the additional distributions and ρε = 3 2 with ρε being the internal energy. Neither p i j ( f ) nor q( f ) would be preserved during the relaxation process, yet can be calculated directly fromf * and Following the ES-BGK target distribution, we can summarize the algorithmic steps of the devised second-order time integration, whose computatinoal complexity is almost identical to the first-order conventional SP-BGK methods: 1. Initialize the particles in the simulation domain (similar to DSMC) method to obtain f t .
3. Move the particles off * in the physical domain (and apply the boundary conditions similar to DSMC) to con- 5. Repeat steps 3 and 4 for the entire duration of the simulation.

Exponential approximation
The devised scheme can be further improved by applying exponential relaxation instead of the linear approximation of Eq. (16). The two additional distributions are then modified with the prefactorss resulting in a modified γ (exp) for the normalization Thus the normalized additional distributions for the exponential approximation are given aŝ Interestingly, this exponential approximation just swaps the pre-factors in Eqs. (43) and (44) compared to the Eqs. (26) and (27) of the linear approach . Therefore, this results again inf * ,(exp) = e −ν∆tf * ,(exp) The only difference that arises is the calculation of the pressure tensor and the heat flux fromf * ,(exp) and q( f t ) = ν∆tPr e −ν∆tPr 1 − e −ν∆tPr q(f * ,(exp) ). (47)

Properties of the ED-BGK method
The behaviour of the introduced particle scheme with exponential differencing (ED-BGK) is analyzed in the limits of free-molecular and continuum regimes. In the following and without loss of generality, the discussion is limited to the linear approximation. In the free molecular regime (τ ∆t), Eqs. (26) and (27) simplify tô which is consistent with the collision-less limit implying On the other hand, for the continuum limit (τ ∆t), the distribution relaxes after a time step ∆t based Eqs. (26) and (27) gives (shown in Appendix A) which recovers the Chapman-Enskog approximation of the Navier-Stokes solution as discussed by Guo et al. [7].

III. IMPLEMENTATION
The presented BGK particle method is implemented in the PIC-DSMC code PICLas [28,29] as described in details by Pfeiffer [9]. The main concept of the particle BGK method especially the energy and momentum conservation, is based on the works of Gallis and Torczynski [15,16]. The particles are moved in the physical space sorted in a computational mesh, collide with boundaries and their macroscopic properties are sampled, all similar to DSMC. However, in contrast to the DSMC method, the collision step is replaced by a jump relaxation with the probability according to Eq. (28), towards the target distribution.f * is constructed out of f by using Eq. (26) or (43) for the linear or the exponential approximation, respectively. Next,f * is constructed by moving the particles fromf * in the next time steps. The difference to the common SP-BGK method is that the pressure tensor or heat flux vector must be relaxed according to equations (37) and (38) or (46) and (47) depending whether the linear or the exponential approximation is used before sampling from the target distribution, since these are from f and not fromf * . For the ellipsoidal statistical target distribution, only the pressure tensor, and for the Shakhov target distribution, only the heat flux vector is necessary.
To determine the correct relaxation frequency ν, the well known temperature dependency of the viscosity µ is used, with the reference temperature T re f , and µ re f the reference dynamic viscosity at T re f [30], and ω V HS as the parameter of the variable hard sphere model (VHS). For a VHS gas, the reference dynamic viscosity can be calculated with the VHS reference diameter d re f of the particles The sampling process itself as well as a detailed discussion of the possible energy and momentum conservation schemes for the particle BGK method can be found in Pfeiffer [9]. For momentum and energy conservation, the final velocities of the particles are corrected according to where u = ∑ N i=1 v i /N is the bulk flow velocity before the relaxation, v i the particle velocity after the relaxation (uncorrected), and u = ∑ N i=1 v i /N. Note that, v i = v i , if no relaxation occurs for particle i. Due to Equation (54) ensures the momentum conservation. The energy conservation is achieved by choosing α as with T the temperature before relaxation, and T after the relaxation process (in the absence of energy correction).
In order to achieve second-order accuracy in space, the values that enter the BGK operator are linearly interpolated, instead of being constant values per computational cell. A more detailed discussion, including an alternative to the linear interpolation, especially for particle methods, can be found in Fei et al. [25,26]. In this work, a conventional linear interpolation was used, similar to Particle-In-Cell codes, as have already been described for the PICLas code by Fasoulas et al. [29], among others.

A. Homogeneous relaxation
The first test case, similar to Fei et al. [25], considers a spatially homogeneous relaxation (in a stationary frame of reference) to examine the accuracy of the time integration scheme. In a single adiabatic cell, argon with the particle density n = 2.7 · 10 25 m −3 and temperature T = 273 K is simulated, subject to the initial condition of the 13-moment Grads distribution [31] f (57) This allows to initialize the gas in a non-equilibrium state with stresses and heat fluxes. Each deviatoric entry of the pressure tensor was chosen according to p i j = 3nk B T , whereas the heat flux is set based on q i = n(3k B T ) 3/2 m −1/2 . Relaxation of the pressure tensor and the heat flux vector is analyzed for different time step sizes. The normal SP-BGK method, the proposed exponential differencing SP-BGK (ED-SP-BGK) with the linear approximation (Eq. (16)), and the exponential approximation (Eq. (39)) are employed. A finely resolved DSMC calculation is deployed as the reference. As shown in Fig. 1, the error in the pressure tensor for the standard SP-BGK is significant, as the target distribution is approximated as a constant per time step. Comparing the results depicted in Fig. 2 and Fig. 3, it can be seen that the relaxation is also correctly reproduced for time steps greater than the collision time by the devised ED-SP-BGK model. In the case of the linear approximation, the error in the relaxation of the pressure tensor is greater than the one resulting from the exponential approximation, while the heat flux is reproduced slightly better with the linear approximation.

B. Sod shock tube
The well-known Sod test case is considered in two different regimes. Argon was again used for the simulations. The physical domain has a length of 1 m. The starting temperature in the entire simulation domain is 273 K. However, the particle density differs in the left n l , x ∈ [0, 0.5] and the right n r , x ∈ [0.5, 1] sub-domains. For the rarefied case, the densities n l = 1.508E19 m −3 and n r = 1.885E18 m −3 were chosen. This results in a Knudsen number of about Kn ≈ 0.1. In the continuum case, the imposed densities are n l = 1.508E21 m −3 and n r = 1.885E20 m −3 resulting in a Knudsen number Kn ≈ 0.001. The reference solutions were generated with DSMC.
In the rarefied case, a computational mesh with 60 grid cells was used for both DSMC and ED-SP-BGK simulations. Also the time step ∆t = 3E − 5 s was the same for both methods, resulting in ν∆t ≈ 0.035 for the ED-SP-BGK method. As depicted in Fig. 4, the DSMC result is very well matched with the ED-SP-BGK method.
In the case of the continuum flow, while spatio-temporal discretization of DSMC was adjusted to resolve the collisional scales, the mesh similar to the rarefied setting was used for ED-SP-BGK. Therefore, an adaptive subcell method besides the time step size of ∆t DSMC = 2E − 7 s were used to resolve the mean free path and the mean collision time for the DSMC method. The ED-SP-BGK simulations were performed with the same time step as in the rarefied simulation. Thus, a fac- tor of ν∆t ≈ 3.5 was obtained. As can be seen in Fig. 5, the ED-SP-BGK method achieves similar number density, mean velocity and temperature with significantly coarser discretization. Both the linear and exponential approximations were tested, yielding negligible differences. For the next test cases, therefore, only the linear approximation was used.

C. Couette flow
The next test case is a steady planar Couette flow of argon, confined between the two parallel plates which move in opposite directions with the velocity v wall = ±500 m/s. The initial temperature of the gas as well as the wall temperature are T = 273 K leading to a supersonic Couette flow. The Couette flow was carried out in three different Knudsen number ranges. The conditions are shown in Table I. In the table, the reference time step is also given, for which the solution of the standard SP-BGK method becomes identical to the proposed ED-SP-BGK method. In order to investigate solely the influence of the time integration on the performance of the two methods, identical spatial in- terpolation was performed in both SP-BGK and ED-SP-BGK.
An important feature of wall bounded gas flows is the gas-boundary interaction. In stochastic particle methods, the time and location where a particle hits the boundary are estimated using linear interpolation. This has been shown to reset the scheme back to the first order [32,33], close to the boundaries. Consistent high-order treatment of boundaries normally involve extensive modifications in the underlying stochastic process (e.g. see walk on sphere scheme for the case of random walks [34,35]) Such treatments go beyond the scope of this study, yet it should be noted that a certain reduction in the order of time integration will be observed in following scenarios where boundary interactions become dominant.

Kn=0.1
The first Couette flow is in the rarefied range with Kn = 0.1.  maximum ν∆t ≈ 0.8 is obtained in the simulation region. We observe that the ED-SP-BGK method achieves significantly better results compared to the SP-BGK method, even in the rarefied range where ν∆ < t. In Fig. 7

Kn=0.01
In the transition regime with Kn = 0.01, the performance of the SP-BGK method becomes worse for larger time step sizes as depicted in Fig. 8. The non-dimensional time step reaches ν∆t ≈ 2 for the largest time step size of 16∆t. The ED-SP-BGK method still shows good agreements with the reference solution. The time step is bounded by the CFL condition. The velocity plot depicted in Fig. 9 again shows good agreement for both methods.

Kn=0.001
The last test case enters the continuum range with Kn = 0.001. Here the required time step size to accurately resolve the solution becomes prohibitively small for the SP-BGK method, as shown in Fig. 10. We get ν∆t ≈ 5 corresponding to the time step size of 32∆t, where again the ED-SP-BGK scheme performs reasonably well. For the velocity plots Fig. 11, the picture is the same as before. Both methods can represent the velocity very well.

IV. CONCLUSION
Particle methods comprise an attractive class of algorithms, which are typically employed for rarefied gas simulations. Despite their popularity and robustness, their numerical accuracy and convergence-order hardly go beyond first-order. Multi-scale flow phenomena with realistic rarefaction regimes, however, require high-order schemes and thus entail fundamental improvements on numerical aspects of the stochastic particle methods. This study enhances SP-BGK schemes by introducing a consistent second-order treatment of the BGK relaxation, using the exponential differencing method. The structure of the devised ED-SP-BGK method allows for a straightforward implementation and minimal overhead with respect to the traditional SP-BGK method, as the particle weights remain non-negative irrespective of the spatio-temporal scales. Besides correct asymptotic limits of the proposed ED-SP-BGK scheme in both free-molecular and continuum limits, its performance was investigated in a few standard settings covering a large variation of the Knudsen number. It is demonstrated that the devised scheme performs better than the SP-BGK method in all considered scenarios including both rarefied and continuum regimes. Extension of the scheme for more complex gas flow scenarios besides im-