Unified Gas-kinetic Wave-Particle Method IV: Multi-species Gas Mixture and Plasma Transport

In this paper, we extend the unified gas-kinetic wave-particle (UGKWP) method to the multi-species gas mixture and multiscale plasma transport. The construction of the scheme is based on the direct modeling on the mesh size and time step scales, and the local cell's Knudsen number determines the flow physics. The proposed scheme has the multiscale and asymptotic complexity diminishing properties. The multiscale property means that according to cell's Knudsen number the scheme can capture the non-equilibrium flow physics in the rarefied flow regime, and preserve the asymptotic Euler, Navier-Stokes, and magnetohydrodynamics limit in the continuum regime. The asymptotic complexity diminishing property means that the total degree of freedom of the scheme automatically decreases as cell's Knudsen number decreases. In the continuum regime, the scheme automatically degenerates from a kinetic solver to a hydrodynamic solver. In UGKWP, the evolution of microscopic velocity distribution is coupled with the evolution of macroscopic variables, and the particle evolution as well as the macroscopic fluxes are modeled from the time accumulating solution up to a time step scale from the kinetic model equation. For plasma transport, current scheme provides a smooth transition from particle in cell (PIC) method in the rarefied regime to the magnetohydrodynamic (MHD) solver in the continuum regime. In the continuum limit, the cell size and time step of the UGKWP method is not restricted to be less than the mean free path and mean collision time. In the highly magnetized regime, the cell size and time step are not restricted by the Debye length and plasma cyclotron period. The multiscale and asymptotic complexity diminishing properties of the scheme are verified by numerical tests in multiple flow regimes.


Introduction
Gas mixture and plasma widely exit in the universe and are extensively applied in the industry of aerospace, chemical, and nuclear engineering. Both gas mixture and plasma transport have multiscale flow dynamics. For the gas mixture, the flow regime various from the rarefied to continuum regime according to the Knudsen number. In the rarefied regime, the fundamental governing equation is the multi-species Boltzmann equation [1], which resolves the physics on the mean free path and mean collision time scale. The complex five-fold integral collision operator makes the Boltzmann equation difficult for both mathematical analysis and numerical simulation. Therefore, many kinetic models have been proposed, for example the McCormack model [2] that linearizes the nonlinear collision term with the assumption for the distribution function to slightly deviate from equilibrium; the Andries-Aoki-Perthame model [3] in which the collision in modeled by a single relaxation term; and other modified models that can recover transport coefficients correctly [4,5]. Although the kinetic equation resolves small scale flow physics, the high dimension of the equation puts barrier in practical 3D engineering applications.
The hydrodynamic model, namely the Euler or Navier-Stokes (NS) equations are mostly used in the continuum regime. For the plasma transport, the flow regime varies from rarefied regime to continuum regime according to the Knudsen number, and varies from the two fluid regime to magnetohydrodynamic (MHD) regime according to the normalized Larmor radius and Debye length. In the rarefied flow regime with large Knudsen number, the plasma flow physics is described by the kinetic Fokker-Planck-Landau equation coupled with the Maxwell equation [6]. In the hydrodynamic regime at small Knudsen number, the two-fluid hydrodynamic system coupled with Maxwell equation can describe the plasma flow dynamics in a more effective way, which takes into account the Hall effect, electron inertia effect, resistive effect, etc. [7]. In the highly magnetized flow regime where the normalized Larmor radius approaches zero and Debye length on the order of the reciprocal of the speed of light, a single fluid ideal MHD can be used to approximate the large scale plasma flow dynamics [8]. For both multiscale gas mixture and plasma transport, the hydrodynamic models are more effective, but limited in the continuum regime; while the kinetic models capture small scale physics, but have complex form and high dimension. Therefore, in order to capture flow physics in different regime in the corresponding most efficient way, the construction of multiscale model for gas mixture and plasma transport is highly demanded.
In general, the numerical methods for gas mixture and plasma transport can be categorized into the deterministic method and stochastic method. The deterministic discrete ordinate method (DVM) has great advantage in the simulation of low speed flow as it does not suffer from any statistical noise [9]. In the last several years, many deterministic numerical methods have been developed for multi-species gas mixture [10][11][12][13], as well as plasma transport [14,15]. On the other hand, when dealing with the high speed flow and 3D flow, the stochastic particle method shows advantage in term of computation efficiency. The direct simulation Monte Carlo (DSMC) method has been extended to gas mixture and chemical reaction [16]. For the simulation of plasma transport, the particle in cell (PIC) method has been developed and widely applied in industry [17]. For the traditional DVM, DSMC, and PIC methods, the numerical cell size is usually required to be less than the mean free path and Deybye length, and the time step is required to be less than the mean collision time. The cell size and time step constraints reduce the computational efficiency of the traditional DVM, DSMC, and PIC methods and it becomes impossible to use them in the continuum regime. In order to remove the constraints, the asymptotic preserving schemes have been proposed that can preserve the flow dynamics in the collisionless and Euler limiting regime [10,18].
The unified gas-kinetic scheme (UGKS) proposed by Xu et al. is a multiscale numerical numerical method for the simulation of gas flow [19,20]. In the last decade, the UGKS has been well developed and extended to the field of multiscale photon transport [21], plasma transport [15], gas-particle multiphase flow [22], neutron transport [23], etc. The two important ingredients of UGKS are: firstly, the evolution of velocity distribution function is coupled with the evolution of the macroscopic conservative variables; secondly, the numerical flux of UGKS is constructed from the integral solution of the kinetic equation which takes into account both particle free stream and collision effects. The UGKS has been proved to be a second order unified preserving scheme that can accurately capture the NS solution with cell size and time step being much larger than the mean free path and mean collision time [24], the same as traditional NS solvers in discretizing the macroscopic equations directly. To improve the efficiency of UGKS in the simulation of high speed flow, the unified gas-kinetic wave-particle (UGKWP) method has been proposed and applied in the simulation of multiscale gas dynamics and photon transport [24][25][26]. The construction of UGKWP method follows the direct modeling methodology of UGKS: the evolution of microscopic simulation particle is coupled with the evolution macroscopic variables, and the multiscale particle evolution equation is derived from the integral solution of the kinetic equation. The propose of this work is to extend the UGKWP method to the field of multiscale gas mixture and plasma transport.
The rest of the paper is organized as following. The governing equations for gas mixture and plasma transport are discussed in Section 2. In Section 3, the UGKWP methods for gas mixture and plasma transport are proposed. The unified preserving and asymptotic complexity diminishing properties of UGKWP are discussed in Section 4. The numerical examples are shown in Section 5, and the last Section 6 gives the conclusion.

Governing equations for multi-species gas mixture and plasma transport
This section is to present the governing equations based on which the scheme is constructed.
The multi-species Boltzmann equation is first reviewed, and then the kinetic model equation proposed by Andries, et al. [3] will be discussed, including its asymptotic behavior in continuum regime. The two fluid kinetic-Maxwell system, as well as the Hall-MHD equations will be presented as well.

Multi-species Boltzmann equation
A gas mixture composed of m species can be modeled by the multi-species Boltzmann equations, where f α (t, x, v) is the velocity distribution function of species α, and the collision between species α and k follows the integral operator where µ αk = mαm k mα+m k is the reduced mass, and n is the unit vector joining the centers of the two colliding spheres. The collision kernel B α,k depends on relative velocity v r = v − v * .
The macroscopic density ρ α , velocity U α , temperature T α , and energy E α of species α can be calculated by taking the moments of the velocity distribution f α , where m α and n α are the molecular mass and number density of species α. The total density ρ, total number density n, total momentum ρ U , and total energy E satisfy Boltzmann equation is a fundamental equation that describes the mean free path level flow physics, however, the five-fold collision operator is costly numerically. Simplified kinetic model equations are developed in the literature [27][28][29][30], including a relaxation-type kinetic model proposed by Andries, et al. [3]. Andries' model will be introduced in the next section, based on which the numerical schemes for multi-species gas mixture and plasma are constructed.

Kinetic model equation for multi-species gas mixture
The relaxation-type kinetic model equation that originally proposed by Gross and Krook [31] has been widely used in the numerical simulation of rarefied gas dynamics due to its simple formulation. Such BGK-type operator has been extended to model the multi-species collision by Andries, Aoki, and Perthame [3], which can be written as where the post collision distribution function is a Maxwellian distribuion and the parameters T * α and U * α are chosen to recover the exchanging relations for Maxwell molecule, which takes the form For Maxwell molecules, the interaction coefficient χ and relaxation parameter τ satisfy where a αk is the constant of proportionality in the intermolecular force law [32]. The advantage of Andries' kinetic model is that it satisfies the indifferentiability principle, entropy condition, and can recover the exchanging relationship of Maxwell molecules with such a simple relaxation form [3].
Based on Andries' model, Liu et al. proposed a BGK-Maxwell system for fully ionized plasma transport [15], which can be written as where the velocity distribution f α (t, x, v) of species α (α = i for ion and α = e for electron) is governed by a kinetic equation that coupled with the Maxwell equations for electromagnetic wave. In the Maxwell equation, E and B are the electric and magnetic field, c is speed of light, and 0 is the vacuum permittivity. In the kinetic equation, the Lorenz acceleration a α takes the form where e is electric charge, and m α is the particle mass of species α. The post collision distribution g α takes the same form of the Andries' model as given in Eq. (5), however the interspecies interaction coefficient χ ie is determined by the plasma electrical conductivity σ p [15] The hydrodynamic equations such as the Navier-Stokes equations, the Euler equations, and magnetohydrodynamic equations can be derived in the continuum regime. The asymptotic behavior of above kinetic model Eq.(4) and Eq.(7) will be briefly discussed in the next subsection.

Asymptotic behavior of the kinetic system
In this section, the asymptotic analysis is applied to give the corresponding hydrodynamic limits of the Andries' and BGK-Maxwell equations. Given the reference variables length L ∞ , temperature T ∞ , mass m ∞ , number density n ∞ , and magnetic field strength B ∞ , the following reference variables can be deduced, which are the reference velocity, time, density, electric field, acceleration, and velocity distribution respectively. Based on above reference variables, the Andries' kinetic model can be re-scaled as and the BGK-Maxwell system can be re-scaled as where the variables with a tilde stand for the re-scaled variables, and especiallyr andλ D are the normalized Larmor radius and Debye length, For the sake of simplicity, the tilde is omitted in the following parts of the paper.
In the continuum regime, the Andries' kinetic equation recovers the Navier-Stokes and Euler equations as τ → 0, According to the Chapman-Enskog theory [1], the distribution of Andries' kinetic model can be expanded as whereḡ is the Maxwellian distribution of the averaged quantities of all species that are evaluated from Eq.(3). The zero-th order expansion with respect to τ α gives the Euler equations [3], and the first order expansion gives the Navier-Stokes equations [3], The mass diffusion flux J α is and the shear stress σ and heat flux q satisfy  The BGK-Maxwell system converges to the two-fluid system and magnetohydrodynamics system as τ → 0 and r → 0, In the continuum regime with τ α 1, and χ ie ∼ 1, the distribution of BGK-Maxwell system can be expanded as where g α is the Maxwellian distribution of the macroscopic quantities of species α. The first order expansion gives hydrodynamic two-fluid equations where the strain rate tensor σ( U ) is The viscosity µ α and the thermal conductivity κ α can be expressed by the relaxation parameter τ α as In above two-fluid system, S i = −S e and Q i = −Q e are the corresponding momentum and energy exchange between electron and ion, In the magnetohydrodynamic regime with m e m i , the first order with respect to r, the zeroth order with respect of τ α and m e /m i of the two-fluid system give the Hall-MHD equations, where j = en i U i − en e U e is the current density, and σ is the electrical conductivity that relates to the interspecies interaction coefficient χ ie as given in Eq.(8). In the limit where λ D = c −1 and r → 0, one gets the ideal MHD equations, The asymptotic behavior of the Andries' and BGK-Maxwell system is given in above discussion.
In the next section, the unified gas-kinetic wave-particle method for gas mixture and plasma transport will be proposed.

UGKWP method for multi-species gas mixture
The unified gas-kinetic wave-particle method is a multiscale numerical method that preserves the asymptotic limits of the Andries' kinetic equations. The UGKWP method couples the evolution of the velocity distribution f α and the macroscopic quantities W α . The evolution of microscopic distribution and macroscopic variables will be given in the following subsections.

The evolution of microscopic velocity distribution function
Similar to the UGKWP method for single species gas [24], in current scheme the velocity distribution function is partially represented by an analytical distribution g +,c α and partially represented by stochastic particles P αk = (m αk , x αk , v αk ), which is as shown in Fig. 1. Here m αk is the mass of simulation particle P αk , which represents a cluster of real gas particles of species α, and x αk , v αk stand for the position and velocity of simulation particle P αk . The evolution of the microscopic velocity distribution function follows the integral solution of the kinetic equation (4).
can be written as where the equilibrium distribution is integrated along the characteristics Substituting the second order Taylor expansion of equilibrium into the integral solution, the numerical multiscale evolution solution for simulation particle can be obtained, where A physical interpretation of Eq.(17) is that a particle has a probability e −t/τ to free stream in a time period [0, t], and has a probability 1 − e −t/τ to interact with other particles and reach a velocity distribution g + α . The free stream particles are kept and the collisional particles get re-sampled from distribution g + α . The cumulative distribution function of the free streaming time t f is from which t f can be sampled as t f = −τ ln(η) with η a uniform distribution η ∼ U (0, 1). For a time step ∆t, the particles with t f ≥ ∆t will be collisionless particles, and the particles with t f < ∆t will be collisional particles. The procedure of updating particles in UGKWP method is Step 1: Sample free streaming time t f,αk for each particle P αk , and stream particle P αk for a time period of min(∆t, t f,αk ); Step 2: Keep collisionless particles, and remove collisional particles. Calculate the total conservative quantities of collisional particles W h i,α from the updated conservative quantities Step 3: Rebuild the microscopic velocity distribution. Calculate the analytical distribution g +,c α and re-sample collisional particles from distribution g +,f α .
In above particle updating procedure, total conservative quantities of collisionless particles in cell Ω i is denoted as W p i,α , and the total conservative quantities of collisional particle in cell Ω i is denoted as W h i,α . In the distribution rebuilding process, the which coresponding to the collisional and collisionless particles in the next time step from t n to t n+1 . The distribution g +,c α is recorded analytically, and the distribution g +,f α is re-sampled into stachastic particles. Above discussion gives the evolution of particles, and in the next subsection we will give the evolution of the conservative variables.

The evolution of macroscopic quantities
The evolution of macroscopic quantities is under the framework of finite volume scheme. The The finite volume scheme of W i,α follows where l s ∈ ∂Ω i is the cell interface with center x s and outer unit normal vector n s . The numerical flux of conservative variables F s,α at x s can be written as is the conservative moments of distribution function with ξ the internal degree of freedom. The time dependent distribution function f α ( x s , t, v) at cell interface is constructed from the integral solution of kinetic equation as given in Eq. (16). The above UGKWP flux for conservative variables can be split into the equilibrium flux and the free streaming flux First, we consider the equilibrium flux F eq s which can be calculated directly form the macroscopic flow field. Assume x s = 0 and t n = 0, the equilibrium g can be expanded as where g 0,α = g α (0, 0, v). The initial equilibrium g 0,α and its spatial and time derivatives can be obtained from the micro-macro consistency and compatible condition where g l α and g r α are the equilibrium distributions according to the reconstructed left and right side conservative variables at cell interface W l α , W r α , and ∇ x W α is the reconstructed spatial derivative of conservative variables at cell interface. In this paper, the van Leer limiter is used to achieve a second order accurate space reconstruction. Substitute the reconstructed equilibrium distribution Eq.(23) into the equilibrium flux Eq.(21), and we have where the time integration coefficients are Next we consider the free stream flux F f s,α . As stated in the last subsection, the initial distribution is represented partially by an analytical distribution g +,c α , and partially by particles, and 13 therefore the free stream flux F f s,α is also calculated partially from the reconstructed analytical distribution as F f,w s,α , and partially from particles as F f,p s,α . The initial analytical distribution g +,c α is reconstructed as which gives where the time integration coefficients are The net particle flux F f,p s,α is calculated as is the index set of the particles stream out cell Ω i during a time step, and P ∂Ω + i ,α is the index set of the particles stream in cell Ω i . Finally, the finite volume scheme for conservative variables is To solve W n+1 i,α from Eq. (27), the following two linear system needs to be solved. The first is the m × m linear system for m species velocity vector V n+1 The second m × m linear system is for m species internal energies e n+1 i = (e n+1 i,1 , e n+1 i,2 , ...e n+1 i,m ) Under the assumption of non-vacuum solutions (ρ n i,α > 0), each system admits a unique solution. The evolution of the microscopic velocity distribution and macroscopic quantities compose the UGKWP method for multi-species gas mixture.

UGKWP method for plasma transport
In this subsection, the UGKWP method for plasma transport will be proposed, which is the UGKWP method for multi-species coupled with the electromagnetic field. We split the BGK-Maxwell equations into the transport equations and the interaction equations. The transport equations including the electron ion transport and electromagnetic wave transport can be written as and the interaction equations are In the next two subsections, the numerical evolution equations for the transport equations and interaction equations will be presented respectively.

Evolution equations for the transport equations
In the transport equations, the electron and ion transport is decoupled from the electromagnetic wave transport. The numerical evolution equation for the electron and ion transport is the UGKWP method presented in Section 3. The Yee-grid based Crank-Nicolson scheme proposed by Yang el al. is used as the evolution equation for the electromagnetic wave transport [33].
The semi-implicit discretization of transverse electric wave equation on Yee mesh can be written And the semi-implicit discretization of magnetic wave equation is Substituting Eq.(29) and Eq.(30) into Eq.(31), an implicit equation for B z can be derived as which can be effectively solved by Douglas-Gunn algorithm. The advantage of the Yee-grid based Crank-Nicolson scheme is that the divergence constraint of the Maxwell equation is numerically preserved; the dispersion and dissipation error is lower than the FDTD method; and the scheme is unconditionally stable, which removes the CFL constraint on time step.

Evolution equations for the interaction equations
Taking conservative moments on Eq.(28), one gets the macroscopic interaction equations The implicit discretization of the macroscopic interaction equations gives the following linear system, from which the electromagnetic field and macroscopic flow variables are updated to t n+1 , and the velocity of the simulation particles is updated by The evolutions of the transport equations and interaction equations compose of the UGKWP method for plasma transport.

Unified preserving and asymptotic complexity diminishing properties of UGKWP method
In this section, the multiscale property of UGKWP method will be discussed, and the computational complexity will be estimated. Guo et al. proposes the unified preserving property which assesses the accuracy of a kinetic scheme in continuum regime [34]. Crestetto et al. proposes the asymptotic complexity diminishing property of a kinetic scheme which assesses the computational complexity of a kinetic scheme in continuum regime [35]. In the following proposition, we show that the UGKWP method is a second order UP scheme and an asymptotic complexity diminishing scheme. Proposition 4.1. Holding the mesh size and time step, the UGKWP method satisfies: 2. The scheme becomes a second order scheme for Navier-Stokes equations τ → 0. 3. The total degree of freedom of the scheme N f → N h f as τ → 0, where N h f is the freedom of the hydrodynamic equations.

Proof.
1. In the collisionless limit, we have Therefore, all particles will be streamed for min(∆t, t f,α ) = ∆t. And the UGKWP method solves collisionless Boltzmann equation in collisionless regime. 2. In the continuum regime when τ → 0, we have The analytic flux F an α of macroscopic variables, namely the equilibrium flux and free streaming flux by analytic distribution function satisfies The sampled particle mass in UGKWP method is e −∆t/τ ρ h α Ω x and therefore the net free streaming flow contributed by particles passing through the cell interface, F f,p s,α ∼ O(e −∆t/τ ), diminishes. As τ → 0, Eq.(27) exponentially converges to It can be observed that the numerical flux of conservative variables is consistent with the Navier-Stokes flux given by first order Chapman-Enskog expansion Eq. (9). Therefore in the continuum regime, the UGKWP method converges to Eq.(34), which is a second order gas-kinetic Navier-Stokes solver [36], i.e., the same as the direct macroscopic NS solver in smooth flow region. 3. As τ → 0, the total mass of simulation particle e −∆t/τ ρ h Ω x → 0, and therefore the number of simulation particles N p → 0 in continuum regime. As τ → 0, the total degree of freedom N f = N h f + N p → N h f , and the UGKWP method is an asymptotic complexity diminishing scheme.

Asymptotic preserving property of UGKWP method for plasma transport
Property 4.1 states that the UGKWP method for plasma transport preserves the two fluid model in the hydrodynamic regime. In this subsection, the behavior of the UGKWP method in the highly magnetized regime is discussed. Proposition 4.2. In the highly magnetized regime as r → 0, λ D = c −1 , the linear system Eq.(32) is consistent to the magnetohydrodynamic equations.
Proof. The Crank-Nicolson scheme for electromagnetic wave propagation gives The implicit discretization of the macroscopic interaction equations Eq.(32) gives and therefore the total momentum equation gives which converges to a consistent MHD scheme.

Numerical tests
Five numerical tests are carried out in this section to verify the performance of the UGKWP method in various flow regime, including three 1D and two 2D tests. Firstly, the shock structure of binary gas mixture is calculated to show the capability of the UGKWP method in capturing the flow non-equilibrium in the rarefied regime. The second test is the Landau damping and two steam instability, showing that the scheme can accurately capture the interaction between plasma and electromagnetic field. The Brio-Wu and Orszag-Tang tests verifies the performance of the UGKWP method in different flow regimes. In the last, the scheme is applied to the magnetic reconnection problem to study how the electron ion collision affects the reconnection rate.

Shock structure of binary gas mixture
Normal shock structure is a standard test that verifies the ability of the scheme in capturing the non-equilibrium effect in rarefied regime. In this test, the Mach number is set as M = 1.5, the mass ratio of gas mixture is m B /m A = 0.5, diameter ratio d B /d A = 1, and the component concentration of B is χ B = 0.1. The hard sphere model is used and the reference mean free path is defined by For each component, the upstream and downstream conditions are related through Rankine-Hugoniot condition. The cell size is chosen to be ∆x = 0.5λ ∞ , and CFL number is 0.95. The mass of simulation particle is m p,α = 10 −2 , which corresponds to around a hundred of simulation particles per cell. The normalized density, velocity and temperature are compared to the reference UGKS solution [37], as shown in Fig.2. The UGKWP results well agree with the UGKS solution, which shows the capability of UGKWP in capturing the non-equilibrium flow physics.

Landau damping and two steam instability
The Landau damping and two steam instability are two classical phenomenons that have been well studied theoretically, and therefore we choose these two cases to test the accuracy of UGKWP method in capturing the interaction between plasma and electromagnetic field. First we consider the Landau damping. Consider a Vlasov-Poisson system that perturbed by a weak signal. The linear theory of Landau damping can be applied to predict the linear decay of electric energy with time [6]. The initial condition of linear Landau damping is with α = 0.01. The length of the domain in the x direction is L = 2π/k. The background ion distribution function is fixed, uniformly chosen so that the total net charge density for the system is zero. When perturbation parameter α = 0.01 is small enough, the Vlasov-Poisson system can be approximated by linearization around the Maxwellian equilibrium. The analytical damping rate of electric field can be derived accordingly. Numerical cell number in physical space is N x = 128, and the particle number in each cell is N p = 1000. We test our 20 scheme with different wave numbers and compare the numerical damping rates with theoretical values. For wave numbers k = 0.3 and k = 0.4, the evolution of the L 2 norm electric field is plotted in Fig. 3. It can be observed that the decay rates and oscillating frequencies ω = 1.16, 1.29 agree well with theoretical data.
Once a larger perturbation α = 0.5 and k = 0.5 is applied, the linear theory breaks down, and the nonlinear phenomenon occurs. The evolution of the electric energy calculated by UGKWP method is shown in Fig. 4 (a), The linear decay rate of electric energy is approximately equal to γ 1 = −0.287, which agrees well to the values obtained by Heath et al. [38]. The growth rate predicted by UGKWP method is approximately γ 2 = 0.078, which is between the value of 0.0815 computed by Rossmanith and Seal and 0.0746 by Heath et al. [39].
Next we consider the linear two stream instability problem with initial distribution function: with α = 0.001 and k = 0.2. The length of the domain in the x direction is L = 2π k . The background ion distribution function is fixed, uniformly to balance the charge density of electron.
After an initial transition, a linear growth rate of electric field can be theoretically predicted [6].
We apply the UGKWP method to calculate this two stream instability problem with physical cell number N x = 512 and simulation particle number N p = 1000 per cell. The electric energy result is shown in Fig. 4 (b), and good agreement between UGKWP solution and theoretical value can be observed. The phase space contour at t = 70 is shown in Fig. 5. Compared to the UGKS solution, the UGKWP solution provides more detailed flow structure.

Brio-Wu shock tube
The Brio-Wu shock tube is originally designed for MHD solvers in continuum regime. Here we calculate the Brio-Wu problem in rarefied (Kn=1), transitional (Kn=10 −2 ), and continuum (Kn=10 −4 ) regimes. The same initial condition as the Brio-Wu one is shown in Fig. 6. The ion to electron mass ratio is set to be 1836, and the ionic charge state is set to be unity. The normalized Debye length is λ D = 0.001, the normalized Larmor radius is r = 0.001, and the normalized speed of light is 1000. The grid points in physical space are N x = 1000. And the simulation particle mass is set as m p = 10 −5 . The UGKWP solutions in rarefied and transitional flow regimes are shown in Fig. 7 and 8, and compared to the reference UGKS solution. In Fig. 9, the UGKWP solution in continuum regime is compared to the reference UGKS solution.
The UGKWP solutions have good agreement with reference solutions in different flow regime.
Especially, it can be observed that the statistical noise significantly reduces as Knudsen number decreases thanks to the asymptotic complexity diminishing property of UGKWP. In the MHD regime, the UGKWP solution is shown in Fig. 10 and compared to the reference ideal-MHD solution.

Orszag-Tang Vortex
The Orszag-Tang Vortex problem was originally designed to study the MHD turbulence [40]. In this work, the problem is calculated in rarefied (Kn=1) and continuum (Kn=10 −4 ) regimes to verify the multiscale and asymptotic complexity diminishing property of UGKWP. The initial data for the current study is m i /m e = 25, n i = n e = γ 2 , P i = P e = γ, B y = sin(2x), in rarefied regime is shown in Fig. 11 compared to the reference UGKS solution, and the UGKWP result in continuum regime is shown in Fig. 12. A better agreement and lower noise can be observed in the continuum regime due to the asymptotic preserving and the asymptotic complexity diminishing property of UGKWP. In the MHD regime, the UGKWP solution with Kn=10 −4 and r = 0 are shown in Fig. 13-15, where in Fig. 15(b) the UGKWP pressure distribution along y = 0.625π is compared to the MHD solution [41].

Magnetic reconnection
Magnetic reconnection is an important phenomenon that transfers magnetic energy into flow energy by topological change of magnetic lines. In this test case, the UGKWP method is used to study the reconnection phenomenon in different flow regime, and study how the particle collision affects the collision rate as well as the topology of magnetic line. The simulation uses the same initial conditions as the GEM challenge problem [42]. The initial magnetic field is given by and a corresponding current sheet is carried by the electrons The initial number densities of electron and ion are n e = n i = 1/5 + sech 2 (y/λ).
The electron and ion pressures are set to be  Fig. 18. In the continuum regime, the topology of the magnetic field is symmetric, while in the transitional regime a magnetic island appears in the middle region at ω p t = 15 and merges into the big right island at ω p t = 30. Due to the magnetic island, two x-shape reconnection points form and the reconnection rate in the transitional regime is significantly increased 15 < ω p t < 30.
After ω p t = 30 when the middle magnetic island merges with the right one, the reconnection rate slows down to the same reconnection intensity as in the continuum regime, which is shown in Fig. 18(b).

Conclusion
In this work, we extend the unified gas-kinetic wave-particle method to the field of multi-species gas mixture and multiscale plasma transport. The construction of numerical scheme for multiscale transport is based on the direct modeling methodology [20], where the flow physics is In conclusion, the UGKWP method has great potential to solve multiscale transport problems in rarefied flow [24,25], radiative transfer [26,43], and plasma physics.

Sample collisional particles
Analytical distribution