{'fraction_non_alphanumeric': 0.04765716191429048, 'fraction_numerical': 0.015441053860263464, 'mean_word_length': 4.302564102564102, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, '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}}
Uncertainty quantification for the BGK model of the Boltzmann equation using multilevel variance reduced Monte Carlo methods

We propose a control variate multilevel Monte Carlo method for the kinetic BGK model of the Boltzmann equation subject to random inputs. The method combines a multilevel Monte Carlo technique with the computation of the optimal control variate multipliers derived from local or global variance minimization problems. Consistency and convergence analysis for the method equipped with a second-order positivity-preserving and asymptotic-preserving scheme in space and time is also performed. Various numerical examples confirm that the optimized multilevel Monte Carlo method outperforms the classical multilevel Monte Carlo method especially for problems with discontinuities.

1. Introduction. Kinetic theory, from a statistical physics viewpoint [5], represents an essential tool to model the non-equilibrium dynamics in a variety of fields including rarefied gases, semiconductors, plasmas, and even large particle systems in biological and social sciences [4,28,21]. The most fundamental kinetic equation, the Boltzmann equation, describes the statistical behavior of a thermodynamic system by taking into account particle transport and binary collisions [3]: where f = f (x, v, t) is the phase space distribution function of position x, velocity v and time t. The collision term Q(f, f ) is a high-dimensional, quadratic integral operator. ε is the Knudsen number, defined as the ratio of the mean free path and the typical length scale. In most applications, ε varies from O(1), the kinetic regime, to ε ≪ 1, the fluid regime. Albeit wide applicability of the Boltzmann equation, the complexity of the collision operator Q(f, f ) makes both analysis and computation of the equation extremely challenging. Hence many simplified collisional models have been introduced to mimic the properties of the full Boltzmann operator. Among these, the Bhatnagar-Gross-Krook (BGK) model [1], which assumes a simple relaxation to equilibrium, has been widely used. The model reads as follows: where M [f ] is the so-called Maxwellian equilibrium function given by where ρ(x, t), U (x, t), T (x, t) are the density, bulk velocity, and temperature defined through the moments of f : Let Φ(v) = 1, v, 1 2 |v| 2 T , then one has the following conservation property: where m is the momentum and E is the total energy. Using (1.5), if we multiply (1.2) by Φ(v) and integrate over v, we obtain the following local conservation law: When ε → 0, formally we have f → M [f ] from (1.2). Replacing f by M [f ] in the above local conservation law yields the compressible Euler equations: In the last decades, research activities in kinetic theory have focused mainly on deterministic kinetic equations, both theoretically and numerically [4,28,9], ignoring the presence of uncertain/random inputs. In reality, uncertainties may arise in initial/boundary conditions and other parameters, like the details of the microscopic interaction, because of incomplete knowledge or imprecise measurement. Recently, there has been a significant interest to study the impact of these random inputs in kinetic equations, see [14] and the whole collection [16] for an overview. In order to quantify the above uncertainties, the construction of numerical methods for kinetic equations has been mostly oriented on stochastic Galerkin approximation based on generalized Polynomial Chaos expansion (gPC-sG), already successfully applied to many physical and engineering problems [13,29]. We mention that recently gPC-sG methods have been successfully applied also to direct simulation Monte Carlo methods for the Boltzmann equation [22]. Despite the fact that gPC-sG methods have been able to show spectral accuracy for smooth solutions, they suffer the drawback of the curse of dimensionality and their highly intrusive nature. On one hand, existing codes for simulating the deterministic kinetic problems need to be completely reconfigured to implement the gPC-sG method. On the other hand, intrusiveness can induce some non-physical approximations even when the deterministic numerical solvers possess the correct physical properties. For example, due to the gPC expansion, the methods may induce approximations with non positive density; furthermore, close to fluid regimes, it is well-known that the gPC-sG system may lose hyperbolicity and lead to spurious solutions [7].
Another class of methods for uncertainty quantification is based on statistical Monte Carlo (MC) sampling, where the random space is sampled and the underlying deterministic PDE is solved for each sample. The non-intrusiveness of the method enables the approximated solutions to inherit properties, like positivity preservation, of the existing deterministic solvers and makes the parallel computing feasible for implementation. However, the asymptotic convergence rate is non-improvable by the central limit theorem and accelerated algorithms are obtained through variance reduction techniques [2]. In this context, multifidelity methods for kinetic equations have been recently introduced in [10,11,17], see also the recent survey [23] for an introduction to the topic. These methods are capable to provide a significant speedup of the convergence properties of the Monte Carlo solver using as control variates simplified surrogate models that are cheaper to solve than the full model. A related line of research is based on the use of multilevel Monte Carlo (MLMC) methods (see for instance [20,19] and references therein for these methods applied to hyperbolic conservation laws), where the approximation of statistical expectation breaks up into telescopic sums of expectations of consecutive mesh sizes. These methods are closely related to multifidelity methods, since they essentially use in a recursive way the solution of the full model with various coarser meshes as surrogate models.
In this manuscript, following the above analogy we develop multilevel Monte Carlo methods in a control variate setting for the multiscale kinetic equations. Therefore, in our MLMC method each level in the telescopic sum depends on an additional parameter which is computed in order to minimize the variance of the solver. We will perform this strategy, both locally between two different levels, as well as globally among all levels. As a prototype kinetic equation to design our methodology we consider the BGK model (1.2) of the Boltzmann equation subject to random inputs. Following the well-posedness results in [24,25], we provide a direct analogue of the former to the BGK equation with random parameters. Due to the non-intrusiveness of MC type methods, approximations of the statistical moments can preserve properties from the deterministic solvers. We adopt the Implicit-Explicit Runge-Kutta (IMEX-RK) scheme from [15] to construct a second-order positivity-preserving (the distribution f is positive for all ε) and asymptotic-preserving (the scheme becomes a solver for the This manuscript is for review purposes only. limiting Euler sytem (1.7) when ε goes to zero) scheme for time and spatial discretizations.
Various numerical examples confirm the good performance of MLMC methods compared to MC methods and that the control variate MLMC method outperforms the classical MLMC method especially for problems with discontinuities.
The rest of this paper is organized as follows. In the next section, we introduce the BGK equation with random inputs and establish the well-posedness of the equation. The Monte Carlo methods and analysis are presented in section 3, whereas in section 4 we discuss their multilevel extension in a standard and control variate setting. In section 5 we show the numerical results obtained with standard MC, MLMC and control variate MLMC methods. Finally some conclusions are drawn in section 6. In a separate Appendix A we report the details of the dimension reduction method and the numerical scheme adopted to solve the deterministic BGK equation.
2. The BGK equation with random inputs. In this section we formulate systematically the BGK equation with random inputs and establish the well-posedness of the equation by extending the results in [24,25].
2.1. Setup of the problem. In the BGK equation, due to the uncertain parameters coming from the initial condition or boundary condition, the resulting solution f would be a random variable taking values in the functional space, where the solution of the BGK equation (1.2) lies in. In most circumstances, it is the physical observables or macroscopic quantities (such as ρ, U , T ) at certain time that are of interest, hence we will mainly consider the random variables taking values in L 1 (D), where D is the physical domain. Following the discussion in [19], we first present some basic concepts from probability theory and functional analysis.
Let (Ω, F , P) be a probability space with Ω being the set of elementary events, F the corresponding σ-algebra, and P the probability measure mapping Ω into [0, 1] such that P(Ω) = 1. A random variable taking values in L 1 (D), a separable Banach space, is defined to be any mapping X : Ω → L 1 (D) such that the set {w ∈ Ω : X(w) ∈ A} = {X ∈ A} ∈ F for any A ∈ G , where (L 1 (D), G ) is a measurable space.
To define the expectation and variance of random variables in L 1 (D), we need the concept of Bochner integral by extending the Lebesgue integral theory. The strong measurable mapping X : Ω → L 1 (D) is Bochner integrable if, for any probability measure P on the measurable space (Ω, F ), Moreover, any Bochner integrable random variable X : Ω → L 1 (D) can be approximated by a sequence of simple random variables {X n } n∈N defined as follows, To get moments like expectation or central moments like variance, similar as the derivation of the Lebesgue integral, the Bochner integral is defined by taking the limit of sequences of This manuscript is for review purposes only.
simple random variables {X n (w)}, for example, the k-th order moments is defined as n→∞ Ω X k n (w) dP(w), and the variance is defined as For the error analysis, we need to introduce the Banach space L p (Ω, F , P; L 1 (D)) with the norm and L ∞ (Ω, F , P; L 1 (D)) with the norm The BGK equation with random inputs hence reads The initial condition is given as For the boundary condition, we consider one of the following: • periodic boundary: f (w; x + a, v, t) = f (w; x, v, t) for x ∈ ∂D and some a ∈ R 3 ; • Dirichlet boundary: f (w; x, v, t) = g(w; x, v, t) for x ∈ ∂D; • purely diffusive Maxwell boundary: for x ∈ ∂D, where n is outward normal of ∂D and M w is given by where T w (w; x, t) is the wall temperature and ρ w (w; x, t) is chosen such that This manuscript is for review purposes only.

2.2.
Well-posedness of the equation and some estimates of the macroscopic quantities. In the following, we establish the well-posedness of the BGK equation (2.7) with random inputs. We also obtain some estimates for the macroscopic quantities ρ, U and T . For simplicity, we assume the periodic boundary condition and consider the uncertainty only arising in the initial condition f 0 .
First of all, some general estimates on the macroscopic quantities can be obtained pointwise in w following [25] for the deterministic BGK equation.
according to (2.9). Moreover, set (2.14) N Then the following estimates hold: where C is a positive constant depending only on q.
Based on the above estimates, one can obtain the existence and uniqueness of the solution to (2.7) also following [25] in a point-wise manner in w.

Theorem 2.2 ([25]). Set
Suppose the initial condition f 0 (w; x, v) ≥ 0 and that for some q > 5, (2.18) then, for fixed Knudsen number ε > 0, there exists a unique mild solution of the initial-value problem (2.7)-(2.10) with periodic boundary condition. Moreover, for all t > 0, the following bounds hold: This manuscript is for review purposes only.
where C is a constant depending only on q.
As a direct consequence of Proposition 2.1 and Theorem 2.2, we have the following corollary on the upper bounds of the macroscopic quantities. Corollary 2.3. Suppose the conditions in Theorem 2.2 hold. We also assume the Knudsen number ε ≥ ε 0 > 0. Then for all t > 0, the following bounds hold: where C 1 and C 2 are positive constants depending only on A 0 , C 0 and q.
Proof. By (2.16), (2.20) and (2.21), we have By (2.15), (2.20) and (2.21), we have Again using (2.15), we have Finally, by (2.25) and (2.20), we have 3. Standard Monte Carlo method. In this section, we describe the basic Monte Carlo sampling method to solve the BGK equation (2.7) and establish some error estimates. For simplicity, we will consider that the uncertainty only comes from the initial condition. The case for the random boundary condition is similar.
3.1. Monte Carlo method. Suppose we generate M independent and identically distributed (i.i.d.) random samples f i 0 , i = 1, . . . , M, according to the random initial condition f 0 (w; x, v). Then each f i 0 (w; x, v) will yield a unique analytical solution to (2.7) at time t, then U i and T i are given by Since it is the macroscopic quantities we are interested in, in the following, without further notice we will use a single variable q to denote ρ, |U | or T . Given the samples q i , i = 1, . . . , M , the MC estimate of the expectation E[q(w; x, t)] is given by To estimate the error between E[q(w; x, t)] and E M [q(w; x, t)], we need the following lemma.
Lemma 3.1. For every finite sequence {Y j } M j=1 of independent random variables with zero mean in L 2 (Ω; L 2 (D)), We have the following consistency theorem.
. As a result, the corresponding copies of macroscopic quantities .
Using the boundedness of domain D, .
Taking the expectation, noting that ∆q i are independent and using Lemma 3.1, we have As a direct result of Theorem 3.2 and Corollary 2.3, we have the following convergence theorem.

Theorem 3.3. Under assumptions of Theorem 2.2 and Corollary
Furthermore, for any M ∈ N + , there holds the error bound Proof. It only needs to note that

Monte
Carlo method with fully discrete scheme. To complete the error analysis, we need to consider the Monte Carlo method coupled with the fully discrete scheme for the BGK equation, which includes discretization in time, physical space and velocity space. The details are given in the Appendix A. Simply speaking, we are using the Gauss quadrature in the velocity space, second order IMEX-RK scheme for time discretization, and second order MUSCL finite volume scheme for spatial discretization (under the hyperbolic CFL condition ∆t ≤ C∆x). Overall, this leads to a second order positivity-preserving and asymptoticpreserving scheme for the deterministic BGK equation. In the following, we assume that the velocity discretization is accurate enough and ignore the work and error in velocity space. It is then reasonable to assume the numerical solution q ∆x,∆t (w; x, t 1 ), computed with mesh size ∆x and time step ∆t corresponding to initial data f 0 (w; x, v) up to time t 1 , satisfies the following error estimate point-wise in w: where the constant C(w) has an upper bound C w .
The MC estimate of the expectation E[q(w; x, t)] is now given by We have It is enough to apply Theorem 3.2 and Assumption 3.1.
The following corollary is a direct result of Theorem 3.4.
Furthermore, for any M ∈ N + , there holds the error bound: 4. Control variate multilevel Monte Carlo method. In this section we first introduce the multilevel Monte Carlo method and then following [11] we discuss the use of control variate techniques to optimize its variance reduction properties locally using two subsequent levels or globally among all levels.
This manuscript is for review purposes only.

Multilevel Monte Carlo method.
The MLMC method is defined as a multilevel discretization in x and t with a level l dependent number of samples M l . Suppose we have a nested triangulation {T l } L l=1 of the spatial domain D (L ∈ N + is the number of levels) such that the mesh size ∆x l at level l satisfies Set q i ∆x 0 ,∆t 0 (w; x, t) := 0, then given a target level L of spatial resolution, the MLMC estimate of the expectation E[q(w; x, t)] is given as follows Hence what we really sample is the difference of solutions at two consecutive levels. At each level l, we generate M l i.i.d. samples f i 0 , i = 1, . . . , M l , of the initial data f 0 on meshes ∆x l and ∆x l−1 respectively, and then use the fully discrete scheme for the BGK equation (2.7) to advance solutions q i ∆x l ,∆t l and q i ∆x l−1 ,∆t l−1 to certain time t. To simplify the notation, we set q ∆x 0 ,∆t 0 (w; x, t) := 0 and define the random variable Y l := q ∆x l ,∆t l (w; x, t) − q ∆x l−1 ,∆t l−1 (w; x, t), and the specific samples Y i l := q i ∆x l ,∆t l (w; x, t) − q i ∆x l−1 ,∆t l−1 (w; x, t). We have the following consistency and convergence results for the estimator (4.2). Proof.
This manuscript is for review purposes only.
For part I, Assumption 3.1 yields, For part II, using Lemma 3.1, Proof. From Theorem 4.1 for l = 1, and similarly for l ≥ 2,

4.2.
Quasi-optimal and optimal multilevel Monte Carlo method. In this section we generalize the previous MLMC method following [10]. To start with, take the 2 level MLMC method for example. Suppose we have a low fidelity (coarse mesh) approximation q 1 and a high fidelity (fine mesh) approximation q 2 , then the 2 level MLMC method with control variate (assume M 2 ≪ M 1 ) reads as follows This manuscript is for review purposes only.
where the multiplier λ has to be determined in order to minimize the variance. It can be shown that for independent samples the optimal value of λ is given by and in the above expression the covariance and variance are estimated directly from the Monte Carlo samples. Generally, suppose we have L levels of solutions {q ∆x i ,∆t i } i=1,...,L , from coarsest level q ∆x 1 ,∆t 1 to finest level q ∆x L ,∆t L . Then the MLMC method with control variates is given by Note that {λ l } L l=1 here are the coefficients to be determined and λ L = 1. If we only consider the variance reduction for each pair of consecutive levels, then we can easily get the analogy of (4.11) to estimate {λ l }, which we refer to as the quasi-optimal MLMC method: (4.13) However, if we focus on minimizing the overall variance of the estimator (4.12) and assume that the levels are independent, then denoting (4.14)λ l = L i=l λ i , l = 1, . . . , L, the optimality conditions yield a tridiagonal system forλ l : Cov[q ∆x l−1 ,∆t l−1 , q ∆x l ,∆t l ] = 0, l = 1, . . . , L − 1, where we assumedλ 0 = 0,λ L = 1 and q ∆x 0 ,∆t 0 = 0. A practical way to solve the above tridiagonal system is to rewrite (4.15) in terms of original λ i . For simplicity, we denote V[q ∆x l ,∆t l ] by V l and Cov[q ∆x l+1 ,∆t l+1 , q ∆x l ,∆t l ] by Cov l to get (4.16) which can be easily solved by recursive substitution. This is what we refer to as the optimal MLMC method.
Denote the correlation coefficient of q ∆x l ,∆t l and q ∆x l+1 ,∆t l+1 by we can prove the following consistency and convergence results for the estimator (4.12): then at time t = t 1 , Proof. The proof is similar to Theorem 4.1. All we need is to note that when λ is quasioptimal, we have for l ≥ 2, This manuscript is for review purposes only.
with the error bound Remark 4.5. We emphasize that the computational cost for quasi-optimal and optimal MLMC is the same as the standard MLMC method. In fact, we can use the data from MLMC to estimate the λ l using (4.13) or (4.16).

Numerical results.
In this section, we present several numerical examples for the BGK equation (2.7) with random initial condition or random boundary condition. The details of the deterministic solver are provided in Appendix A. Simply speaking, we are solving a reduced system (A.6) and (A.7), which is equivalent to the full BGK equation in one spatial dimension. We use the IMEX-RK scheme for time discretization and finite volume scheme for spatial discretization so that the overall method is second order in both time and space. We choose x ∈ [0, 1] and v ∈ [−5, 5], where 40 Legendre-Gauss quadrature points are used in the velocity space to ensure that the error in velocity is negligible. The CFL condition is fixed as ∆t = 0.1∆x.

Error evaluation.
In the following, we assume the uncertainties come from either the initial condition or boundary condition. Since the solution is a random field, the numerical error is a random quantity as well. For error analysis, we therefore compute a statistical estimator by averaging numerical errors from several independent experiments.
More precisely, for each method we perform K = 40 experiments, and get the corresponding approximations {q (j) (x, t)}, j = 1, . . . , K, where q can be ρ, U or T . We approximate the overall error in norm · L 2 (Ω;L 1 (D)) via where q ref (x, t) is the reference solution obtained using the stochastic collocation method [29] with 120 Legendre-Gauss collocation points and N x = 1280 spatial points. We are also interested in the error at each spatial point: Sometimes to better evaluate the error from the random domain, we would like to ignore the error induced by spatial discretization. To achieve so, we consider another kind of reference solution, q rel (x, t), obtained again using the stochastic collocation with 120 collocation points, while in the spatial domain we use the same finest mesh ∆x L as that in the corresponding MLMC method to obtain q (j) (x, t). Therefore, we can assess the error as

Test 1:
Smooth random initial condition. We first consider the BGK equation subject to random initial condition: and the random variable z obeys the uniform distribution on [−1, 1]. The periodic boundary condition is used and the Knudsen number ε = 1.
To determine the number of samples needed in MC and MLMC methods, we proceed as follows.
In the MC method, we consider a series of spatial discretizations: N = 10, 20, 30, 40, and for each case, we vary the sample size as M = 5, 10, 15, ... The results are shown in Figure 1 Figure 1 (right). Roughly we can see that M 1 ≈ 10240 gives the smallest error (the error saturates when the sample size further increases).
In Figure 2 we combine all the previous MC and MLMC results under the scale of workload. Since we are essentially solving 1D BGK problem, the workload for one deterministic run up to certain time with N spatial points is O(N 2 ). Then for the MC method with M samples, the This manuscript is for review purposes only.  quasi-optimal MLMC, and optimal MLMC. The results are shown in Figure 3, from which we clearly see the better accuracy of MLMC methods compared to standard MC. On the other hand, the difference of three MLMC methods are not obvious in this example. Next we examine the errors of the three MLMC methods as defined in (5.2),(5.3). The results are gathered in Figure 4. We can see that the three methods perform equally well in this test (at some points of the domain, the errors of quasi-optimal MLMC and optimal MLMC are slightly better than standard MLMC), largely because the solution is smooth.
To better understand this, we plot the values of λ 1 and λ 2 in the quasi-optimal and optimal MLMC methods in Figure 5. We can see that almost all values are not far from 1, which means the methods are not far from the standard MLMC.

Test 2: Shock tube problem.
In this test, we consider two kinds of shock tube problems with random initial condition. The first one has uncertainty in the interface location: The random variable z obeys the uniform distribution on [−1, 1]. We set the Knudsen number ε = 1.0e − 6 so that the problem is close to the fluid regime. Similarly as the previous example, we need to perform a series of tests to determine the number of samples needed in MC and MLMC methods. Figure 6 shows the analogous tests as those in Figure 1. The main difference from the previous example is that the errors saturate much quicker as the number of samples increases. This is due to the low regularity of the solution so that the error from spatial/temporal discretization dominants easily. In Figure 7 we combine both MC and MLMC results under the scale of workload. Similarly as we observed in Figure 2, with the same workload, the MLMC method can achieve the best accuracy.
This manuscript is for review purposes only.  From the right plot in Figure 6, we also see that M 1 ≈ 320 is the best number of samples at level 1 of MLMC method. Therefore, we fix the set of parameters in the following tests of MLMC methods: mesh sizes N 1 = 10, N 2 = 20, N 3 = 40, and sample sizes M 1 = 320, M 2 = 80, M 3 = 20. In Figures 8-11, we report the results obtained using standard MLMC, quasi-optimal MLMC, and optimal MLMC methods. We examine the approximations to both expectation E[q] and variance V[q] of the macroscopic quantities ρ, U and T . Note that the MLMC methods are based on the linearity of the expectation operator, not the variance operator. Hence to approximate the variance, we approximate separately two different expectations E[q 2 ] and E[q] and use them to obtain The results clearly show that both control variate MLMC methods outperform the standard MLMC in regions where the solution presents strong variations, namely close to the shock position. On the other hand, we did not observe any relevant gain using the optimal MLMC method over the quasi-optimal MLMC.  To better see the difference of three MLMC methods, we plot the values of λ 1 and λ 2 in the quasi-optimal and optimal MLMC methods for both tests in Figure 12 and 13. It is clear that for these problems with shocks/discontinuities the values are far from one in various regions of the computational domain. This is particularly true for the temperature and velocity in agreement with the corresponding errors observed in the previous figures.

Test 3: Sudden heating problem.
In the last test, we consider a problem with random boundary condition. The gas is initially in a constant state with ρ 0 = 1, U 0 = (0, 0, 0),  6. Conclusions. We have introduced a control variate multilevel Monte Carlo method for the BGK model of the Boltzmann equation with uncertainty. Well-posedness of the BGK equation with random parameters, consistency and convergence analysis for various MC type methods are established. Extensive numerical results confirm that the MLMC methods perform much better than the standard MC, and the control variate MLMC is capable to provide further improvement over the conventional MLMC, in particular for problems close to fluid regimes and in presence of discontinuities, where the fidelity degree of the various levels is reduced and traditional gPC-SG based methods may fail. On the other hand, we haven't observed any gain using a global optimal strategy in the variance reduction with respect to a local quasi-optimal strategy based on two subsequent levels, which is subject to future investigation. The approach here developed naturally extends to other kinetic equations of Boltzmann type which combines deterministic discretizations in the phase space with Monte Carlo sampling in the random space. Appendix A. Dimension reduction method and deterministic solver for the BGK equation.
In this Appendix, we briefly describe the dimension reduction method adopted to reduce the computational complexity of the BGK equation and the details of the numerical methods used to discretize time, physical space and velocity space. Since the Monte Carlo methods are non-intrusive, our discussion will be based on the deterministic equation (1.2) for simplicity.
A.1. The Chu reduction method. The BGK equation (1.2) lives in six-dimensional phase space whose computation can be extremely expensive. Under certain homogeneity assumptions, one can reduce the dimension using the so-called Chu reduction [6].
, and U = (U 1 , U 2 , U 3 ). If the physical domain D ∈ R 3 is effective only in one dimension and homogeneous in the other two dimensions (e.g., a tube), it is reasonable to assume the following: Then the equation (1.2) becomes 2T (x 1 , t) . The Chu reduction proceeds by introducing two distribution functions: It is then easy to derive that φ and ψ satisfy the following system: This manuscript is for review purposes only. Denoting R · dv 1 = · , it is easy to see the following relation holds Now our task is to solve the reduced 1D BGK system (A.6)-(A.7).
A.2. The fully discrete scheme. The fully discrete scheme used to solve (A.6)-(A.7) consists of three components: velocity discretization, time discretization, and spatial discretization.
Velocity discretization. In the velocity space, we follow the discrete velocity method (see Section 4.1.1 in [12] or [18] for example), which satisfies a discrete entropy decay property.
We first truncate the infinite velocity domain into a bounded interval [−R, R] and then discretize it using N v -point Gauss quadrature with (ξ k , w k ), k = 1, 2, . . . , N v as abscissae and weights. To obtain M φ , M ψ from φ and ψ, normally one could use the relation in (A.10), where the continuous integral is replaced by the Gauss quadrature. However, due to the domain truncation error, the resulting moments are not sufficiently accurate. To remove this  Time discretization. Due to the possibly stiff collision term, we use the implicit-explicit Runge-Kutta (IMEX-RK) scheme [8,26] for the time discretization. In particular, we employ the second-order IMEX-RK scheme proposed in [15], which is positivity preserving and asymptotic preserving (preserving the Euler limit without ∆t resolving ε).
This manuscript is for review purposes only. where the values of the coefficientsã ij , a ij , α are given in Section 2.6.1 of [15]. To implement the above scheme explicitly, we first solve the moment system (A.14) which is obtained by taking the moments of (A.13) and using (A.12). Hence we can obtain This manuscript is for review purposes only.