A stochastic kinetic scheme for multi-scale plasma transport with uncertainty quantification
AA stochastic kinetic scheme for multi-scale plasma transport withuncertainty quantification
Tianbai Xiao a, ∗ , Martin Frank a a Karlsruhe Institute of Technology, Karlsruhe, Germany
Abstract
Plasmas present a diverse set of behaviors in different regimes. Given the intrinsic multi-scale nature of plasma dynamics, classical theoretical and numerical methods are often em-ployed at separate scales with corresponding assumptions and approximations. Clearly, thecoarse-grained modeling may introduce considerable uncertainties between the field solutionsof flow and electromagnetic variables, and the real plasma physics. To study the emergence,propagation and evolution of randomness from gyrations of charged particles to magnetohy-drodynamics poses great opportunities and challenges to develop both sound theories andreliable numerical algorithms. In this paper, a physics-oriented stochastic kinetic schemewill be developed that includes random inputs from both flow and electromagnetic fields viaa hybridization of stochastic Galerkin and collocation methods. Based on the BGK-typerelaxation model of the multi-component Boltzmann equation, a scale-dependent kineticcentral-upwind flux function is designed in both physical and particle velocity space, andthe governing equations in the discrete temporal-spatial-random domain are constructed.By solving Maxwell’s equations with the wave-propagation method, the evolutions of ions,electrons and electromagnetic field are coupled throughout the simulation. We prove thatthe scheme is formally asymptotic-preserving in the Vlasov, magnetohydrodynamical, andneutral Euler regimes with the inclusion of random variables. Therefore, it can be used forthe study of multi-scale and multi-physics plasma system under the effects of uncertainties,and provide scale-adaptive physical solutions under different ratios among numerical cellsize, particle mean free path and gyroradius (or time step, local particle collision time andplasma period). Numerical experiments including one-dimensional Landau Damping, thetwo-stream instability and the Brio-Wu shock tube problem with one- to three-dimensionalvelocity settings, and each under stochastic initial conditions with one-dimensional uncer-tainty, will be presented to validate the scheme.
Keywords:
Boltzmann equation, plasma physics, kinetic theory, uncertaintyquantification, multi-scale methods ∗ Corresponding author
Email addresses: [email protected] (Tianbai Xiao), [email protected] (Martin Frank)
Preprint submitted to Elsevier June 8, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J un . Introduction Plasma applications cover an extremely wide range of density n from 10 to 10 m − and temperature k B T from 10 − to 10 eV [1]. As F. Chen wrote in his famous monograph[2], ”What makes plasmas particularly difficult to analyze is the fact that the densities fallin an intermediate range. They behave sometimes like fluids, and sometimes like a collectionof individual particles.” A qualitative classification of typical plasma regimes with respectto density and temperature is presented in Fig. 1. In this paper, we confine ourselves tothe study of classical, non-relativistic and weakly coupled plasmas, e.g. the magnetosphere,whose behaviors can be well described by the kinetic theory of gases. Figure 1: An illustrative demonstration of different plasma regimes and related phenomena [1].
Given the intrinsic multi-scale nature in plasma physics, classical theories are devoted todifferent governing equations at different hierarchies. For example, at a particle mean freepath which is larger than the Debye length λ D = ( ε k B T ele /ne ) / , and a mean collisiontime which is larger than the reciprocal of the plasma frequency ω p = ( ne /ε m ion ) / , themotions of ions and electrons can be depicted statistically through kinetic equations, e.g.the Vlasov equation and the Fokker-Planck-Landau equation. On the other hand, at themacroscopic level with intensive intermolecular collisions, the fluid dynamic equations areroutinely used to model the collective behaviors of charged particles, i.e. the magnetohy-drodynamics (MHD) equations.Since the 1950s, rapid development has been made in deterministic numerical meth-ods for plasma simulations [3–10]. However, given the coarse-grained approximation in the2eld theories of plasmas and errors inherited from numerical simulations, considerable un-certainties may be introduced inevitably. One typical example are the uncertain inputs ofthe initial/boundary value problem. Furthermore, for the evaluation of collision kernel inthe kinetic equations, the phenomenological model parameters often need to be calibratedby experiments to reproduce correct transport coefficients, which introduce errors into thesimulations. Another example goes to the vacuum permeability employed in the Maxwell’sequations for electromagnetic fields. In the SI system which has gone into force in 2019 [11],this value is measured experimentally as µ = 8 . × − F · m − , with a relativestandard deviation being 1 . × − .Uncertainty quantification (UQ) is a thriving subject that quantifies one’s lack of knowl-edge concerning a physical reality. It applies itself to answer the challenging questions, e.g.how predictive are the simulation results from the idealized models, and how can one explic-itly assess the effects of uncertainties on the quality of model predictions. Depending on themethodology to model the random variables, the methods for UQ study can be classifiedinto intrusive and non-intrusive ones. In the former case, a series of realizations of ran-dom inputs are generated based on a prescribed probability distribution. Each realizationis solved by a deterministic solver, and then a post-processing is employed to estimated theuncertainties. In contrast, intrusive methods work in a way such that we reformulate theoriginal deterministic system.One prevalent intrusive strategy is the stochastic Galerkin (SG) method [12], in whichthe solutions are expressed into orthogonal polynomials of the input random parameters. Itpromises spectral convergence in random space when the solution depends smoothly on thestochastic parameters [13]. In a nonlinear Galerkin system, all the expansion coefficients areessentially coupled, which becomes cumbersome in massive computations. The stochasticcollocation (SC) method [14], although a non-intrusive method, can be seen as a middle way.It combines the strengths of non-intrusive sampling and SG by evaluating the generalizedpolynomial chaos (gPC) expansions on quadrature points in random space. As a result,a set of decoupled equations can be derived and solved with deterministic solvers on eachquadrature point. Provided the solutions posses sufficient smoothness over random space,the SC methods maintain similar convergence as SG, but suffers from aliasing errors due tolimited number of quadrature points.Although the UQ field has undergone rapid development over the past few years, itsapplications on plasma physics mainly focus on the two limits of Vlasov [15–18] and MHD [19,20] with standard stochastic settings. To the best of the authors’ knowledge, limited work hasbeen conducted on the evolutionary process of uncertainty in multi-scale physics. Given thenonlinear system including intermolecular collisions, initial inputs, fluid-surface interactionsand geometric complexities, uncertainties may emerge from molecular-level nature, developupwards, affect macroscopic collective behaviors, and vice versa. To study the emergence,propagation and evolution of uncertainty poses great opportunities and challenges to developboth sound theories and reliable multi-scale numerical algorithms.It is noticeable that tracking the evolution of stochastic variables with either polynomialchaos or quadrature rules is similar in spirit to solving kinetic equations in phase spacewith moment or discrete velocity methods. The advantages of SG and SC methods can be3ombined when the integrals that are necessary for SG inside the algorithm are computednumerically using SC. In this paper, we follow the strategy proposed in [21] and developa stochastic kinetic scheme for multi-scale plasma transport and we couple it to Maxwell’sequations. Based on the Bhatnagar-Gross-Krook (BGK) type relaxation model for multi-component plasmas, a scale-dependent central-upwind flux function is constructed in bothphysical and particle velocity space, which considers simultaneously the individual particletransports and their collective behaviors. The update of source terms of plasma and elec-tromagnetic fields are solved in a coupled way implicitly. We thus combine the advantagesof SG and SC methods with the construction principle of kinetic schemes, and obtain anefficient and accurate scheme for cross-scale BGK-Maxwell system with uncertainties. Therandomly initial inputs of both flow and electromagnetic fields are considered.The rest of this paper is organized as follows. Sec. 2 is a brief introduction of kinetictheory of plasma and its stochastic formulation. Sec. 3 presents the numerical implemen-tation of the current scheme and detailed solution algorithm. Sec. 4 includes numericalexperiments to demonstrate the performance of the current scheme and analyze some newphysical observations. The last section is the conclusion.
2. Deterministic and stochastic theories
The gas kinetic theory describes the time-space evolution of particle distribution function.With a separate modeling of particle transport and collision processes, the evolution equationof monatomic plasmas writes as ∂f α ∂t + u · ∇ x f α + q α m α ( E + u × B ) · ∇ u f α = Q α , (1)where α = ion, ele denotes a specific ion or electron, ( q α , m α ) are particle charge and mass,and ( E , B ) are electric and magnetic fields respectively. For the Coulomb collisions betweencharged particles, the limiting case of Boltzmann collision integral leads to the Fokker-Planck-Landau operator, Q α = N (cid:88) β (cid:26) ∇ u · (cid:90) Φ ( u − u (cid:48) ) [ f β ( u (cid:48) ) ∇ u f α ( u ) − f α ( u ) ∇ u f β ( u (cid:48) )] d u (cid:48) (cid:27) , (2)where ( u , u (cid:48) ) are the velocities of two classes of particles and Φ ( u ) = ( | u | I − u ⊗ u ) / | u | is a 3 × Q α = ν α ( M α − f α ) , (3)4here ν α is the collision frequency. The equilibrium distribution is defined based on thelocal modified macroscopic variables, i.e., M α = n α (cid:18) m α πk B ¯ T α (cid:19) exp (cid:18) − m α k B ¯ T α ( u − ¯U α ) (cid:19) , (4)where n α is number density, m α is molecular mass and k B is the Boltzmann constant. Thedesign of modified temperature ¯ T α and velocity ¯U α is based on the idea that the macroscopictransfer rates in the moments equations derived from BGK model should be consistent withthat from multi-component Boltzmann equation. For elastic scattering, the evaluation ofmodified variables for Maxwell and hard sphere molecules can be written as, ¯U α = U α + τ α (cid:88) r m r m α + m r ν αr ( U r − U α ) , k B ¯ T α = 32 k B T α − m α ¯U α − U α ) + τ α (cid:88) r m α m r ( m α + m r ) ν αr (cid:20) k B T r − k B T α + m r U r − U α ) (cid:21) , (5)where ν αr is the frequency of intermolecular interactions which can be derived throughspecific molecule models [23], and it determines the relaxation time by τ α = 1 / (cid:80) Nr ν αr .Here we adopt the hard-sphere molecules, i.e., ν αr = 4 √ πn r (cid:18) k B T α m α + 2 k B T r m r (cid:19) / (cid:18) d α + d r (cid:19) , (6)where d is the kinetic molecule diameter.Macroscopic conservative flow variables are related to the moments of the particle dis-tribution function, W α = ρ α ρ α U α ρ α E α = (cid:90) m α f α (cid:36)d u , (7)where E α = ( U α ) / k B T α / m α is total energy density, and (cid:36) = (cid:0) , u , u (cid:1) T is thevector of collision invariants. Hence, macroscopic transport equations can be derived bytaking moments of the kinetic equation with respect to the collision invariants, i.e., ∂n α ∂t + ∇ x · ( n α U α ) = 0 ,∂ ( ρ α U α ) ∂t + ∇ x · ( ρ α U α U α ) = ∇ x · P α + n α q α ( E + U α × B ) + R α ,∂ ( ρ α E α ) ∂t + ∇ x · ( ρ α E α U α ) = ∇ x · ( P α U α ) − ∇ x · q α + n α q α U α · E + H α , (8)5here the source terms R α and H α in the balance laws come from the moments of collisionterm respectively, R α = (cid:90) u m α ν α ( M α − f α ) d u = (cid:88) r m α m r m α + m r n α v αr ( U r − U α ) ,H α = (cid:90)
12 ( u − U ) m α ν α ( M α − f α ) d u = (cid:88) r m α m r ( m α + m r ) n α v αr (cid:20) k B T r − k B T α + m r U r − U α ) (cid:21) . (9)Eq.(8) is consistent with the well-known Braginskii’s two-fluid model [24]. Several sources of uncertainty can be considered in the BGK equation. Here we con-sider uncertain initial and boundary conditions, which turn the deterministic system intostochastic case. We employ the generalized polynomial chaos (gPC) expansion of particledistribution with degree N , i.e., f α ( t, x , u , z ) (cid:39) f αN = N (cid:88) | i | =0 ˆ f αi ( t, x , u )Φ i ( z ) = ˆ f Tα Φ , (10)where i could be a scalar or a K -dimensional vector i = ( i , i , · · · , i K ) with | i | = i + i + · · · + i K . The ˆ f α i is the coefficient of i -th polynomial chaos expansion, and the basis functionsused are orthogonal polynomials { Φ i ( z ) } satisfying the following constraints, E [Φ j ( z )Φ k ( z )] = γ k δ jk , ≤ | j | , | k | ≤ N, (11)where γ k = E [Φ k ( z )] , ≤ | k | ≤ N, (12)are the normalization factors. The expectation value defines a scalar product, E [Φ j ( z )Φ k ( z )] = (cid:90) I z Φ j ( z )Φ k ( z ) (cid:37) ( z ) d z , (13)where (cid:37) ( z ) is the probability density function. In practice, it can be evaluated theoreticallyor with numerical quadrature rule, i.e., E [Φ j ( z )Φ k ( z )] = (cid:88) i Φ j ( z i )Φ k ( z i ) w ( z i ) , (14)where w ( z i ) is the corresponding quadrature weight function in random space. In the fol-lowing we adopt a uniform notation (cid:104) Φ j Φ k (cid:105) to denote the integrals over random space fromEq.(13) and (14). 6iven the correspondence between macroscopic and mesoscopic variables, from Eq.(7)we can derive, W α (cid:39) (cid:90) f αN (cid:36)d u = (cid:90) N (cid:88) i ˆ f αi ( t, x , u )Φ i ( z ) (cid:36)d u = (cid:88) i (cid:18)(cid:90) ˆ f αi (cid:36)d u (cid:19) Φ i (cid:39) W αN = N (cid:88) i ˆ w αi Φ i = ˆ w Tα Φ . (15)After substituting the Eq.(10) into the kinetic equation (1) and (3), and performing aGalerkin projection, we then obtain ∂ ˆ f α ∂t + u · ∇ x ˆ f α + ˆ G α = ˆ Q α = ν α ( ˆ m α − ˆ f α ) , (16)where ˆ Q α is the gPC coefficient vector of the projection from collision operator to thepolynomial basis, ˆ Q α = ν α ( ˆ m α − ˆ f α ) , (17)with ˆ m α being the vector of gPC coefficients of Maxwellian distribution (which dependsimplicitly on the ˆ f α ) and ν α being a deterministic collision frequency. The electromagneticforcing term ˆ G α isˆ G α = ˆ G Tα Φ = N (cid:88) i ˆ G i Φ i , ˆ G i = q α m α (cid:80) Nj (cid:80) Nk (cid:16) ˆ E j + u × ˆ B j (cid:17) ∇ u ˆ f αk (cid:104) Φ j Φ k Φ i (cid:105)(cid:104) Φ i (cid:105) . (18)Notice that the gPC coefficients of both ˆ G α and ˆ Q α are nonlinear functions of the statevariables (cf. Eq.(5)). For a deterministic system, the evaluation of the Maxwellian distribution given in Eq.(4)is straight-forward. However, given a generalized polynomial chaos (gPC) system, the mul-tiplication and division can’t be operated directly on the stochastic moments without modi-fying the orthogonal basis. Starting from a known particle distribution function in Eq.(10),here we draw a brief outline to approximately evaluate the Maxwellian distribution functionin the gPC expansion.1. Derive the macroscopic conservative variables from particle distribution function withgPC expansion, W αN = ρ αN ( ρ α U α ) N ( ρ α E α ) N = N (cid:88) i (cid:18)(cid:90) ˆ f αi ψd u (cid:19) Φ i ; (19)2. Locate conservative variables on quadrature points z j of random space and calculate7rimitive variables, e.g. flow velocity U α ( z j ) = ( ρ α U α ) N ( z j ) ρ αN ( z j ) , (20)and T α ( z j ) = ( ρ α E α ) N ( z j ) − ( ρ α U α ) N ( z j ) / ρ αN ( z j )3 k B ρ αN / m α , (21)and then calculate the modified velocity and temperature via Eq.(5).3. Calculate Maxwellian distribution on quadrature points M α ( z j ) = n α ( z j ) (cid:18) m α πk B ¯ T α ( z j ) (cid:19) e − ¯ λ α ( z j )( u − ¯ U α ( z j )) , (22)and decompose it into a gPC expansion M αN = N (cid:88) i ˆ m αi Φ i , (23)with each coefficient in the expansion being given by a quadrature ruleˆ m αi = (cid:104)M α , Φ i (cid:105)(cid:104) Φ i (cid:105) = (cid:80) j M α ( z j )Φ i ( z j ) (cid:37) ( z j ) (cid:82) I z (Φ i ( z )) (cid:37) ( z ) d z . (24) For the self-consistent problem of plasma dynamics, the evolutions of electric and mag-netic fields ( E , B ) are coupled with the motions of charged particles, which can be describedby the linear Maxwell’s equations in vacuum, ∂ E ∂t − c ∇ x × B = − ε J ,∂ B ∂t + ∇ x × E = 0 , ∇ x · E = σε , ∇ x · B = 0 . (25)Here σ = e ( n i − n e ) is the net charge density, J is the current, and the speed of light isrelated to the permeability and permittivity of vacuum with c = ( µ ε ) − / . To ensurethe divergence constraints in numerical simulations, some techniques can be used in solvingMaxwell’s equations. Here we employ the perfectly hyperbolic Maxwell’s equations (PHM)825], ∂ E ∂t − c ∇ x × B + χc ∇ x φ = − ε J ,∂ B ∂t + ∇ x × E + γ ∇ x ψ = 0 , χ ∂φ∂t + ∇ x · E = σε ,ε µ γ ∂ψ∂t + ∇ x · B = 0 , (26)where φ, ψ are two additional correction potentials, and the propagation speed of errorsfor the divergence of magnetic and electric fields are γc and χc correspondingly. With thestochastic Galerkin formulation, the PHM system can be rewritten as ∂ ˆ E ∂t − c ∇ x × ˆ B + χc ∇ x ˆ φ = − ( ˆ J /ε ) ,∂ ˆ B ∂t + ∇ x × ˆ E + γ ∇ x ˆ ψ = 0 , χ ∂ ˆ φ ∂t + ∇ x · ˆ E = ( ˆ σ/ε ) ,ε µ γ ∂ ˆ ψ ∂t + ∇ x · ˆ B = 0 . (27)The parameters in Maxwell’s equations are assumed to be deterministic.
3. Solution algorithm
The current numerical algorithm is implemented within the finite volume framework.We adopt the notation of cell averaged macroscopic conservative variables and particle dis-tribution function in a control volume, W α ( t n , x i , z k ) = ( W α ) ni,k = 1Ω i ( x )Ω k ( z ) (cid:90) Ω i (cid:90) Ω k W α ( t n , x , z ) d x d z ,f α ( t n , x i , u j , z k ) = ( f α ) ni,j,k = 1Ω i ( x )Ω j ( u )Ω k ( z ) (cid:90) Ω i (cid:90) Ω j (cid:90) Ω k f α ( t n , x , u , z ) d x d u d z , along with the coefficient vector in the gPC expansions,ˆ W α ( t n , x i ) = ( ˆ W α ) ni = 1Ω i ( x ) (cid:90) Ω i ˆ W α ( t n , x ) d x , ˆ f α ( t n , x i , u j ) = ( ˆ f α ) ni,j = 1Ω i ( x )Ω j ( u ) (cid:90) Ω i (cid:90) Ω j ˆ f α ( t n , x , u ) d x d u , i , Ω j and Ω k are the cell area in the discretized physical, velocity and random space.The update of the stochastic Galerkin coefficients for particle distribution function canbe written as, ( ˆ f α ) n +1 i,j =( ˆ f α ) ni,j + 1Ω i (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω i S r ( ˆ F α ) fr,j dt + 1Ω j (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω j S r ( ˆ F α ) fi,r dt + (cid:90) t n +1 t n ( ˆ Q α ) fi,j dt. (28)where ( ˆ F α ) fr is the time-dependent fluxes for distribution function at interface r in physicaland velocity space, S r is the interface area, and ˆ Q fα is the collision term. Taking velocitymoments of Eq.(28), we obtain the corresponding macroscopic system,( ˆ W α ) n +1 i =( ˆ W α ) ni + 1Ω i (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω i S r · ( ˆ F α ) Wr dt + (cid:90) t n +1 t n ( ˆ G α ) Wi dt + (cid:90) t n +1 t n ( ˆ Q α ) Wi dt, (29)where ( ˆ F α ) Wr is the flux functions for macroscopic conservative variables, S r = n S r is theinterface area vector, and ˆ G Wα is the external force terms related to Eq.(8).Notice that the gPC coefficients of different orders are coupled through the externalforce term ˆ G α (17) and the collision term ˆ Q α (18). Instead of solving this large nonlinearsystem, we can locate the gPC system onto the collocation points in random space. It resultsa decoupled system that can be solved efficiently. Therefore, we combine the advantagesof stochastic Galerkin and collcation methods, which is one of the novelties of this paper.To make use of it, in the solution algorithm, we first update the gPC coefficients to anintermediate step t ∗ , ( ˆ W α ) ∗ i = ( ˆ W α ) ni + 1Ω i (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω i S r · ( ˆ F α ) Wr dt, ( ˆ f α ) ∗ i,j = ( ˆ f α ) ni,j + 1Ω i (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω i S r ( ˆ F α ) fr,j dt, (30)which are then evaluated on random quadrature cell Ω k , W ∗ i,k = W ∗ Ni ( z k ) = N (cid:88) m ˆ w ∗ i,m ( z k ) Φ m ( z k ) ,f ∗ i,j,k = f ∗ Ni,j ( z k ) = N (cid:88) m ˆ f ∗ i,j,m ( z k ) Φ m ( z k ) . (31)10fterwards, the collision and forcing term are evaluated on collocation points via( W α ) n +1 i,k = ( W α ) ∗ i,k + (cid:90) t n +1 t n ( G α ) Wi,k dt + (cid:90) t n +1 t n ( Q α ) Wi,k dt, (32)and ( f α ) n +1 i,j,k = ( f α ) ∗ i,j,k + 1Ω j (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω j S r ( F α ) fi,r,k dt + (cid:90) t n +1 t n ( Q α ) fi,j,k dt, (33)where ( F α ) fr is the numerical flux at interface r in velocity space. In the solution algorithmloop, Eq.(32) can be solved first, and then the updated variables at t n +1 can be employedto evaluate the Maxwellian distribution in Eq.(33) implicitly.For plasma transport, the evolution of electromagnetic field should be solved in a coupledway with flow field. The hybrid Galerkin-collocation method is employed as well, where thegPC coefficients are stepped to the intermediate state first,ˆ M ∗ i = ˆ M ni + 1Ω i (cid:90) t n +1 t n (cid:88) r ∆ S r · ˆ F Mr dt, (34)where ˆ F M is the flux functions for the electromagnetic fields ( ˆ E , ˆ B , ˆ φ , ˆ ψ ). Then the sourceterms are solved via, E n +1 i = E ∗ i − eε (cid:90) t n +1 t n (( n U ) ion − ( n U ) ele ) dt,φ n +1 i = φ ∗ i + eε (cid:90) t n +1 t n ( n ion − n ele ) dt. (35) Based on the finite volume framework, a scale-dependent interface flux function is neededin multi-scale modeling and simulation. Different from a purely upwind flux which losesefficiency in the collisional limit, we here develop a kinetic central-upwind flux functionbased on an integral solution of the kinetic model equation. The integral solution originatesfrom Kogan’s monograph on rarefied gas dynamics [26] and has been employed by a seriesof gas-kinetic schemes [27–31]. Let us rewrite the stochastic BGK equation (16) for the gPCcoefficients vector along the characteristics, D ˆ f α Dt + ν α ˆ f α = ν α ˆ m α . (36)We assume that the collision frequency ν α is kept fixed at the value that can be computedfrom the macroscopic variables at the previous timestep. Then the following integral solution11olds along the characteristics,ˆ f α ( t, x , u ) = ν α (cid:90) t ˆ m α ( t (cid:48) , x (cid:48) , u (cid:48) ) e − ν α ( t − t (cid:48) ) dt (cid:48) + e − ν α t ˆ f α (0 , x , u ) , (37)where x (cid:48) = x − u (cid:48) ( t − t (cid:48) ) − a ( t − t (cid:48) ) is the particle trajectory, and x = x − u (cid:48) t − a t is thelocation at initial time t = 0. The above solution indicates a self-conditioned mechanismfor multi-scale gas dynamics. For example, when the evolving time t is much less than themean collision time τ = 1 /ν α , the latter term in Eq.(37) dominates and describes the freetransport of particles. And if t is much larger than τ , the second term approaches to zero,and then the distribution function will be an accumulation of equilibrium state along thecharacteristic lines, which provides the underlying wave-propagation mechanism for hydro-dynamic solutions. In the following, we present a detailed strategy for the construction ofnumerical fluxes. Since the electromagnetic force term a is a stochastic variable, we do anoperator splitting to evaluate fluxes in physical and velocity space.With the simplified notations of physical cell interface x i +1 / = 0 and initial time t n = 0,Eq.(37) along physical trajectories of particles can be rewritten into the following form,ˆ f α ( t, , u j ) = ν α (cid:90) t ˆ m α ( t (cid:48) , x (cid:48) , u j ) e − ν α ( t − t (cid:48) ) dt (cid:48) + e − ν α t ˆ f α (0 , − u j t, u j ) , (38)where ˆ f α (0 , − u j t, u j ) is the initial distribution at each time step.In the numerical scheme, the initial distribution function at cell interface can be obtainedthrough reconstruction, i.e.,( ˆ f α )(0 , ± , u j ) = (cid:40) ( ˆ f α ) Li +1 / ,j , x = 0 − , ( ˆ f α ) Ri +1 / ,j , x = 0 + . (39)The initial distributions ( ˆ f α ) L,Ri +1 / ,j at the left and right hand sides of a cell interface areobtained through the van-Leer limiter.The macroscopic conservative variables in the gPC expansions at the interface can beevaluated by taking moments over velocity space,ˆ w = (cid:88) u j > ˆ f Li +1 / ,j ψ ∆ u j + (cid:88) u j < ˆ f Ri +1 / ,j ψ ∆ u j . The equilibrium distribution at interface can be determined as illustrated in Sec. 2.3.After all coefficients are obtained, the interface distribution function becomesˆ f α ( t, , u j ) = (cid:0) − e − ν α t (cid:1) ( ˆ m α ) j + e − ν α t (cid:104) ( ˆ f α ) Li +1 / ,j H [ u j ] + ( ˆ f α ) Ri +1 / ,j (1 − H [ u j ]) (cid:105) , (40)where H ( u ) is the heaviside step function. The above interface distribution function can be12egarded as a combination of central difference and upwind methods. With the variationof the ratio between evolving time t (i.e., the time step in the computation) and collisiontime τ α = 1 /ν α , the factor e − ν α t plays as a modulator and provides a self-adaptive solutionbetween equilibrium and non-equilibrium physics.After the coefficients of distribution function at all orders are determined, the corre-sponding gPC expansion can be expressed as, f αN ( t, , u j ) = N (cid:88) m =0 ˆ f m ( t, , u j )Φ m ( z ) , (41)and the corresponding fluxes of particle distribution function and conservative flow variablescan be evaluated via F fαN ( t, , u j , z ) = u j f αN ( t, , u j , z ) , F WαN ( t, , z ) = (cid:90) u f αN ( t, , u , z ) (cid:36)d u (cid:39) (cid:88) w j u j f αN ( t, , u j , z ) (cid:36) j , (42)where u j denotes quadrature points in particle velocity space, and w j is its integral weightin velocity space, and the time-integrated fluxes in Eq.(29) and (28) can be evaluated withrespect to time in Eq.(40). Besides the flow variables, the numerical fluxes of electromagnetic fields in Eq.(27) arecalculated by the wave-propagation method developed by Hakim et al. [32].
Besides the construction of the interface flux, the source terms need to be evaluatedinside each control volume within each time step. In this part, we show the detailed updatealgorithm for the collision and external force term with stochastic collocation method.
Given the intermediate macroscopic variables in the discrete cell (Ω i , Ω j , Ω k ), the macro-scopic system writes as follows,( ρ α U α ) n +1 i,k = ( ρ α U α ) ∗ i,k + ∆ tq α ( n α ( E + U α × B )) n +1 i,k + (cid:90) t n +1 t n (cid:16) ν α ρ α ( U ∗ α − U ∗ α ) (cid:17) i,k dt, ( ρ α E α ) n +1 i,k = ( ρ α E α ) ∗ i,k + ∆ tq α ( n α E · U α ) n +1 i,k + (cid:90) t n +1 t n (cid:16) ν α ( ρ α E ∗ α − ρ α E ∗ α ) (cid:17) i,k dt, E n +1 i,k = E ∗ i,k − e ∆ tε (( n U ) ion − ( n U ) ele ) n +1 i,k ,φ n +1 i,k = φ ∗ i,k + e ∆ tε ( n ion − n ele ) n +1 i,k . (43)13e split the above system into two parts. First, the flow variables are evolved with respectto mixture source terms,( ρ α U α ) ∗∗ i,k = ( ρ α U α ) ∗ i,k + (cid:90) t n +1 t n (cid:16) ν α ρ n +1 α ( U ∗ α − U ∗ α ) (cid:17) i,k dt, ( ρ α E α ) ∗∗ i,k = ( ρ α E α ) ∗ i,k + (cid:90) t n +1 t n (cid:16) ν α ρ n +1 α (cid:16) E ∗ α − E ∗ α (cid:17)(cid:17) i,k dt. (44)The relaxation integrals are evaluated through the Rosenbrock method [33] to overcomepossible stiffness.Then we solve the electromagnetic sources implicitly, in which a linear system can besolved, i.e., ( ρ α U α ) n +1 i,k = ( ρ α U α ) ∗∗ i,k + ∆ tq α ( n α ( E + U α × B )) n +1 i,k , ( ρ α E α ) n +1 i,k = ( ρ α E α ) ∗∗ i,k + ∆ tq α ( n α E · U α ) n +1 i,k , E n +1 i,k = E ∗ i,k − e ∆ tε (( n U ) ion − ( n U ) ele ) n +1 i,k ,φ n +1 i,k = φ ∗ i,k + e ∆ tε ( n ion − n ele ) n +1 i,k . (45) The updated macroscopic variables can be used to step the particle distribution functionsof ion and electron implicitly. Let us consider the kinetic equation, ∂f α ∂t + a α · ∇ u f α = ν α ( M α − f α ) , (46)where the electromagnetic force is a α = q α m α (cid:0) E n +1 + u α × B n +1 (cid:1) . Making the simplified notations of velocity cell interface u j +1 / = 0 and initial time t n = 0 again, we write the integral solution of Eq.(46) as, f α ( t, x i , , z k ) = ν α (cid:90) t M ∗ α ( t (cid:48) , x i , u (cid:48) , z k ) e − ν α ( t − t (cid:48) ) dt (cid:48) + e − ν α t f ∗ α (0 , x i , − a α t, z k )= M ∗ α (0 , x , u − a α t, z k ) (cid:0) − e − ν α t (cid:1) + f ∗ α (0 , x i , u − a α t, z k ) e − ν α t . (47)Similar to physical space, the above interface distribution function can also be regarded asa combination of central difference and upwind methods in particle velocity space. Withthe variation of the ratio between evolving time t (i.e., the time step in the computation)and collision time τ α = 1 /ν α , it provides a self-adaptive solution from equilibrium to non-equilibrium. Based the above solution, the interface flux in particle velocity space can be14onstructed as F fα ( t, x i , , z k ) = a α f α ( t, x i , , z k ) , (48)and the time-integrated flux can be evaluated directly with respect to t .For the collision term, with the updated flow variables at t n +1 , an implicit update can bearranged. Therefore, the update algorithm of particle distribution function can be writtenas, ( f α ) n +1 i,j,k = ( f α ) ∗ i,j,k + 1Ω j (cid:90) t n +1 t n (cid:88) S r ∈ ∂ Ω j S r ( F α ) fi,r,k dt + ∆ tν n +1 α (cid:0) ( M α ) n +1 i,j,k − ( f α ) n +1 i,j,k (cid:1) . (49) In the current scheme, the time step is determined by the Courant-Friedrichs-Lewy con-dition in phase space,∆ t = C min (cid:18) min( | ∆ x | )max( | u | ) + max( | U | ) , min( | ∆ x | ) c , min( | ∆ u | )max( | a | ) (cid:19) , (50)where C is the CFL number, u = ( u ion , u ele ) is particle velocity, U = ( U ion , U ele ) is fluidvelocity, and c is speed of light. In the computation, c usually takes a pseudo value whichis less than 3 . × m/s but a few orders larger than particle velocity. In this part we will present theoretical analysis on the current kinetic model and numer-ical algorithm, with special focus on their asymptotic limits.
Let us return to the BGK equation (1). From Eq.(6), we see that the collision frequency ispositively correlated with plasma density and temperature. When the plasma gets rarefied,the intensity of collision term decreases correspondingly. As the collision frequency ν α goesto zero, the BGK equation automatically reduces to the Vlasov equation, i.e., ∂f α ∂t + u · ∇ x f α + q α m α ( E + u α × B ) · ∇ u f α = 0 . (51)On the other hand, as the collision frequency increases, then the two-fluid system (8) isequivalent with the BGK equation, which describes the motions of plasma as fluid. Underthe quasineutral assumption n i (cid:39) n e with λ D (cid:28) L ( L is the characteristic length ofsystem), m i (cid:29) m e and ν i , ν e (cid:29)
0, it can be further degenerated to the single-fluid Hall-15HD equations, ∂ρ∂t + ∇ x · ( ρ U ) = 0 ,∂ ( ρ U ) ∂t + ∇ x · ( ρ UU ) = −∇ x p + J × B ,∂ ( ρ E ) ∂t + ∇ x · ( ρ E U ) = ∇ x · ( p U ) + J · E , E + U × B = η J + 1 en ( J × B + ∇ x p ) ,∂σ∂t + ∇ x · J = 0 . (52)where σ is electric charge density and J is net current. Together with Maxwells equations,this set is able to describe the equilibrium state of the plasma.The specific resistivity η is proportional to the interaction frequency between ions andelectrons, i.e., η = m e ne ν ie . (53)Therefore, given the fully conductive conditions with ν ie → ∂ρ∂t + ∇ x · ( ρ U ) = 0 ,∂ ( ρ U ) ∂t + ∇ x · ( ρ UU ) = −∇ x p + ( B · ∇ x ) B µ − ∇ x (cid:18) B µ (cid:19) ,∂ ( ρ E ) ∂t + ∇ x · ( ρ E U ) = ∇ x · ( p U ) + 1 µ ρ U · ( ∇ x × B × B ) ,∂ B ∂t + ∇ x × ( U × B ) = 0 . (54)In contrast to the ideal MHD regime, we can also consider the non-conductive limitwhere the interspecies molecular interaction is intensive with ν ie → ∞ . As a result, theplasma now behaves like dielectric material, and the two-fluid system deduces to the Eulerequations, ∂ρ∂t + ∇ x · ( ρ U ) = 0 ,∂ ( ρ U ) ∂t + ∇ x · ( ρ UU ) = −∇ x p,∂ ( ρ E ) ∂t + ∇ x · ( ρ E U ) = ∇ x · ( p U ) . (55) Besides the theoretical modeling, the asymptotic preserving property for the limitingsolutions is also a preferred nature for the kinetic scheme. Let us consider the interfacefluxes constructed in Sec. 3.2 and 3.3 first. In the collisionless limit where ν α approaches16ero, the relation ν α ∆ t (cid:28) f α ) i +1 / ,j,k = ( f α ) Li +1 / ,j H ( u j ) + ( f α ) Ri +1 / ,j (1 − H ( u j )) , ( f α ) i,j +1 / ,k = f ∗ α ( x i , u − a α t, z k ) , (56)which constitute a fully upwind scheme for the Vlasov equation.On the other hand, in the hydrodynamic regime with intensive collisions, the interfacedistribution becomes, ( f α ) i +1 / ,j,k = ( M α ) i +1 / ,j , ( f α ) i,j +1 / ,k = M ∗ α ( x i , u − a α t, z k ) , (57)leading to near-equilibrium MHD or Euler solutions. Besides the flux functions, the treatment of source terms plays an important role for theasymptotic property of scheme. Let us rewrite the electromagnetic sources in Eq.(45) asfollows, ( ρ ion U ion ) n +1 i,k = ( ρ ion U ion ) ∗∗ i,k + ∆ te ( n ion ( E + U ion × B )) n +1 i,k , ( ρ ele U ele ) n +1 i,k = ( ρ ele U ele ) ∗∗ i,k − ∆ te ( n ele ( E + U ele × B )) n +1 i,k , ( ρ ion E ion ) n +1 i,k = ( ρ ion E ion ) ∗∗ i,k + ∆ te ( n ion E · U ion ) n +1 i,k , ( ρ ele E ele ) n +1 i,k = ( ρ ele E ele ) ∗∗ i,k − ∆ te ( n ele E · U ele ) n +1 i,k , E n +1 i,k = E ∗ i,k − e ∆ tε ( n ion U ion − n ele U ele ) n +1 i,k ,φ n +1 i,k = φ ∗ i,k + e ∆ tε ( n ion − n ele ) n +1 i,k . (58)Adding the momentum and energy equations of ion and electron, we get( ρ U ) n +1 i,k − ( ρ U ) ∗∗ i,k = ∆ te (( n ion − n ele ) E ) n +1 i,k + ∆ t (( J ion + J ele ) × B ) n +1 i,k , ( ρ E ) n +1 i,k − ( ρ E ) ∗∗ i,k = ∆ t (( J ion + J ele ) · E ) n +1 i,k . (59)With the quasineutral assumption, it can be simplified as( ρ U ) n +1 i,k − ( ρ U ) ∗∗ i,k = ∆ t ( J × B ) n +1 i,k , ( ρ E ) n +1 i,k − ( ρ E ) ∗∗ i,k = ∆ t ( J · E ) n +1 i,k . (60)At the same time, let us multiply the first equation in Eq.(58) by m ele and second by m ion and subtract the latter from the former, which results, m ion m ele ( n ( U ion − U ele )) n +1 i,k − m ion m ele ( n ( U ion − U ele )) ∗∗ i,k = ∆ te ( m ion + m ele ) ( n E ) n +1 i,k + ∆ te ( n ( m ele U ion + m ion U ele ) × B ) n +1 i,k . (61)17iven m ion (cid:29) m ele , the above the equation can be degenerated to E n +1 i,k + U n +1 i,k × B n +1 i,k = 1 eρ n +1 i,k (cid:16) m ion m ele e ( J n +1 i,k − J ∗∗ i,k ) + m ion J n +1 i,k × B n +1 i,k (cid:17) , (62)which is the generalized Ohm’s law. At a time scale where the inertial effects (i.e., cyclotronfrequency) are unimportant, the time difference of current is ignorable. The second termis the Hall current term, which becomes unimportant in the ideal MHD limit. Let ν ie → E n +1 i,k + U n +1 i,k × B n +1 i,k = 0 . (63) The flowchart of the current solution algorithm is summarized in Fig. 2.
4. Numerical experiments
In this section, we will present some numerical results. The goal of numerical experimentsis not simply to validate the performance of the current scheme, but also to present andanalyze new physical observations. In order to demonstrate the multi-scale nature of the al-gorithm, simulations from Vlasov to magnetohydrodynamics (MHD) regimes are presented.The following dimensionless flow variables are introduced in the simulations,˜ x = x L , ˜ t = tL /U , ˜ m = mm ion , ˜ n = nn , ˜ U = U U , ˜ T = TT , ˜ P = P m ion n U , ˜ q = q m ion n U , ˜ f = fn U , ˜ u = u U , ˜ B = B B , ˜ E = E B U , ˜ σ = σen , ˜ J = J en U , ˜ λ D = λ D r g , where U = (cid:112) k B T /m ion is the thermal velocity of ions, λ D = (cid:112) ε k B T /n e is the Debyelength, and r g = m ion U /eB is the gyroradius of ion in the reference state. For brevity,the tilde notation for dimensionless variables will be removed henceforth. Therefore, the18 tart initialize plasma and eletromanetic fieldcalculate time stepreconstruct distribution function and electromagnetic fieldscalculate flow fluxes F W and F f update W *, f*, and M * by fluxesupdate W ** by interspecies collisionupdate W n+1 and M n+1 from electromagnetic source term timeout / converge?end calculate electromagnetic fluxes F M calculate electromagnetic force a n+1 and flux F f in velocity spacecalculate Maxwellian M n+1 update f n+1 from velocity flux and BGK term Figure 2: Flowchart of solution algorithm. dimensionless BGK-Maxwell system becomes, ∂f α ∂t + u · ∇ x f α + 1 r g m α ( E + u × B ) · ∇ u f α = ν α ( M α − f α ) ,∂ E ∂t − c ∇ x × B = − λ D r g J ,∂ B ∂t + ∇ x × E = 0 , χ ∂φ∂t + ∇ x · E = σλ D r g , c γ ∂ψ∂t + ∇ x · B = 0 . (64)19 .1. Landau damping The Landau damping is a physical phenomenon first predicted by Landau [34] basedon theoretical derivation, namely the exponential decay effect of electromagnetic waves incollisionless plasmas. Here, we use the example of Landau damping to verify the performanceof the current stochastic scheme in the Vlasov limit. For brevity, we consider the one-dimensional case first, and thus the Maxwell’s equations for the electrostatic field degeneratesinto the Poisson equation for the electric potential.
Consider the following macroscopic system n ion n ele UT t=0 =
11 + αξ cos( kx )01 , where α is the amplitude of electromagnetic wave, k is the wave number, and ξ is a randomparameter. The corresponding particle distribution functions at t = 0 are the Maxwellian, f α = n α (cid:16) m α πT (cid:17) / exp (cid:16) − m α T ( u − U ) (cid:17) . The computational setup is listed in Table 1. Given the minor amplitude α , the Vlasov-Poisson system can be regarded as the Maxwellian plus lineazised perturbation near equi-librium. With the large mass ratio, we fix the slow motions of ions as background, and solvethe evolution of electrons. The standard non-intrusive stochastic collocation method is alsoemployed to provide reference solution. Besides, the damping rate of electric field energycan be derived by linear theory [5]. We hereby combine α and ξ into a new perturbationamplitude β , and derive the following relation as benchmark, E ( x, t ) (cid:39) β × . e − . t sin( kx ) cos(1 . t − . . (65) Table 1: Computational setup of linear Landau damping. t x N x m ion /m ele u ele N u α (0 ,
40] [ − π/k, π/k ] 128 1836 [ − ,
5] 128 0.01 k ξ
Polynomial N z (gPC) N z (quad) Kn cfl0.5 U (0 ,
1) Legendre 5 9 100 0.2Fig. 3 presents the time evolution of electric energy. As is shown, the current stochastickinetic scheme provides equivalent numerical solution as the standard collocation method.20he expected energy damping rate is consistent with the theoretical damping rate − . ω = 1 . x = 0. From the zoom-in expected distribution function around u = 0, we see that thelow-speed particles resonate with the electromagnetic wave and gradually absorb energy,resulting in decreased value of particle distribution. In spite of the minor variation ofexpected value around the Maxwellian for the linear damping case, the standard deviationof solutions shows an increasingly symmetric oscillation in the velocity space during theresonance process. It indicates a more significant sensitivity compared to the mean field, andprovides a quantitative description of the non-equilibrium effects triggered by the particle-wave resonance. The same initial conditions and computational setups as linear case are followed, exceptfor the enhanced amplitude α = 0 .
5. The evolution of electric energy is presented in Fig.5. As the amplitude of the electromagnetic wave increases, nonlinear effects would emergecorrespondingly, resulting in a rise in energy after the initial damping.Fig. 6 shows the time evolution of expected particle distribution over the phase space( x, u ), and Fig. 7 picks out the expectation and variance of particle distribution at physicaldomain center x = 0. Given the increasing intensity of radio field, here the particle distri-bution function is deformed in phase space. Consistent with the linear case, the change ofvariance here is more significant than expectation value, indicating a stronger sensibility. The two-stream instability is another typical phenomenon in for collisionless plasmas.To a certain extent, it can also be regarded as an inverse phenomenon of Landau damping.
Consider the following initial system, n ion n ele UTf ion f ele t=0 = (cid:104) αξ (cid:16) cos(2 kx )+cos(3 kx )1 . + cos( kx ) (cid:17)(cid:105) n ion (cid:0) m ion πT (cid:1) / exp (cid:0) − m ion T ( u − U ) (cid:1) n ele (cid:0) m ele πT (cid:1) / (1 + 5 u ) exp (cid:0) − m ele T ( u − U ) (cid:1) , Table 2: Computational setup of linear two-stream instability. t x N x m ion /m ele u ele N u α (0 ,
70] [0 , π/k ] 128 1836 [ − ,
5] 128 0.001 k ξ
Polynomial N z (gPC) N z (quad) Kn cfl0.5 U (0 ,
1) Legendre 5 9 100 0.2The evolution of electric field energy is shown in Fig. 8, from which it can be seen thatthe numerical and theoretical solutions [2] fit well, and the uncertainties of electric energypropagate in the similar pattern as mean field. Fig. 9 provides the particle distributionfunction over the phase space ( x, u ) at t = 70. The swirling pattern is clearly identifiedin both expectation and variance values. More fine structures can be seen in the standarddeviation, which provides a clearer way to quantify the stochastic evolution of particles. In the following let us consider the nonlinear case. The initial system is given as, n ion n ele UTf ion f ele t=0 = (1 + αξ ∗ cos( kx ))00 . n ion (cid:0) m ion πT (cid:1) / exp (cid:0) − m ion T ( u − U ) (cid:1) n ele (cid:0) m ele πT (cid:1) / (cid:0) exp (cid:0) − m ele T ( u − . (cid:1) + exp (cid:0) − m ele T ( u + 0 . (cid:1)(cid:1) . The computational setups adopted here is shown in Table 3.
Table 3: Computational setup of nonlinear two-stream instability. t x N x m ion /m ele u ele N u α (0 ,
70] [0 , π/k ] 256 1836 [ − ,
5] 256 0.001 k ξ
Polynomial N z (gPC) N z (quad) Kn cfl0.5 U (0 ,
1) Legendre 5 9 100 0.2Fig. 10 shows the time evolution of expected particle distribution over the phase space( x, u ). Due to the increasing kinetic energy from two particle streams, the particle dis-tribution function is stretched and warped in the phase space. The pattern of standarddeviation is similar as mean field, yet presenting more fine-scale structures, indicating astronger sensitivity with respect to stochastic parameters.22 .3. Brio-Wu shock tube
After the validation of current scheme in the Vlasov limit, let us turn to another lim-iting regime, i.e. the magnetohydrodynamic (MHD) transport problem under continuumassumption. In this case, we employ the Brio-Wu shock tube [7] as benchmark case for idealMHD solutions.In the simulation, the initial background of deterministic solutions is consistent with theSod problem, with an additional magnetic discontinuity, i.e., n ion n ele UpE x E y E z B x B y B z φψ L = . , for the left half, and n ion n ele UpE x E y E z B x B y B z φψ R = . . . . − , for the right half. The computational setup is listed in Table 4. To recover the ideal MHDsolutions, here the interspecies coefficients ν ie is set as zero. First, we consider the uncertainties of plasma density in the left half tube, with n ion,L = n ele,L = ξ . The numerical expected solutions of macroscopic variables under different ref-erence gyroradius at t = 0 . r g is small. At r g = 0 . able 4: Computational setup of Brio-Wu shock tube. t x N x m ion /m ele u ion u ele cfl(0 , .
1] [0 ,
1] 400 1836 [ − ,
5] [ − , × (cid:113) m ion m ele ξ Polynomial N z (gPC) N z (quad) Kn r g λ D U (0 . , .
05) Legendre 5 9 1 . × − [0 . , r g increases, the electric and Lorentz force decreases accordingly. Theeffects of charge separation are enhanced, and the omitted terms, e.g. the Hall current term,become important. The plasma gradually manifests itself as dielectric material. As a result,the magnetic diffusion is enhanced, and the wave structures from discontinuous magneticfield B y are flattened. When r g = 100, at this point the gyroradius is much larger thanthe characteristic length of the shock tube, and thus plasma behaves as neutral gas, andthe Brio-Wu shock tube degenerates into a standard Sod case. In Fig. 11(m) and (n), wesee the current scheme provides equivalent solutions as Euler equations, indicating the APproperty of the current scheme in the Euler limit.Fig. 12 presents the standard deviations at the same output instant. Generally, theuncertainties travel along with the wave structure of expectation values and present sim-ilar propagating patterns. Typical wave structures serve as sources of local maximums ofvariance. Compared with the expected value, its variance is more sensitive to physical dis-continuities and holds finer-scale structures due to the spectrum formulation in the randomspace. As a result, the overshoots near contact discontinuity and shock, as well as theoscillations around central compound wave are observed.Fig. 13 presents the expectations and standard deviations of particle distribution func-tion at different reference gyrorarius. As is shown, the discontinuities in macroscopic ex-pectations and overshoots in standard deviations come from the uncertainties contained inthe particle distribution function near the center of velocity space. From MHD to Eulerregimes, the randomness on particles get reduced and smoothed, resulting in gentle profilesof macroscopic quantities. In the second case, we turn to the uncertainties from magnetic field in the left half tube,with B y,L = ξ . The numerical solutions of macroscopic variables under different referencegyroradius at t = 0 . r g increases, the wave structures around tube center becomes evenmore complicated in the transition regime r g = 0 .
01, and simplify again when it comes to r g = 0 .
1, which demonstrates the complex nonlinearity from Hall current, two-fluid effects,etc. The increasing gyroradius doesn’t lead to a simple process of monotonic variation.This test case clearly shows the consistency and distinction of propagation modes be-tween expectation value and variance. It also illustrates the capacity of current scheme tosimulate multi-scale and multi-physics plasma transports, and capture the propagation ofuncertainties in different regimes.
5. Conclusion
Plasma dynamics is associated with an intrinsic multi-scale nature due to the largevariations of particle density and temperature, as well as the characteristic scales of thelocal structures. Based on the multi-component BGK model and Maxwell’s equations, astochastic kinetic scheme with hybrid Galerkin-collocation strategy has been constructed inthis paper, which allows for a unified numerical simulation for multi-scale plasma physics.Based on the cross-scale modeling, the solution algorithm is able to capture both equilibriummagnetohyrodynamics and non-equilibrium gyrations of charged particles simultaneously,and recover the scale-dependent plasma physics along with the emergence, propagation, andevolution of randomness. The asymptotic-preserving property of the scheme is validatedthrough theoretical analysis and numerical tests.
Acknowledgement
The current research is funded by the Alexander von Humboldt Foundation.25 eferences [1] National Research Council, Plasma Science Committee, Plasma 2010 Committee, et al.
Plasma science:advancing knowledge in the national interest , volume 3. National Academies Press, 2008.[2] Francis F Chen.
Introduction to plasma physics . Springer Science & Business Media, 2012.[3] Vahid Vahedi and Maheswaran Surendra. A monte carlo collision model for the particle-in-cell method:applications to argon and oxygen discharges.
Computer Physics Communications , 87(1-2):179–198,1995.[4] Francis Filbet, Eric Sonnendr¨ucker, and Pierre Bertrand. Conservative numerical schemes for the vlasovequation.
Journal of Computational Physics , 172(1):166–187, 2001.[5] Nicolas Crouseilles, Thomas Respaud, and Eric Sonnendr¨ucker. A forward semi-lagrangian method forthe numerical solution of the vlasov equation.
Computer Physics Communications , 180(10):1730–1745,2009.[6] Jing-Mei Qiu and Andrew Christlieb. A conservative high order semi-lagrangian weno method for thevlasov equation.
Journal of Computational Physics , 229(4):1130–1149, 2010.[7] Moysey Brio and Cheng Chin Wu. An upwind differencing scheme for the equations of ideal magneto-hydrodynamics.
Journal of computational physics , 75(2):400–422, 1988.[8] Kenneth G Powell, Philip L Roe, Timur J Linde, Tamas I Gombosi, and Darren L De Zeeuw. Asolution-adaptive upwind scheme for ideal magnetohydrodynamics.
Journal of Computational Physics ,154(2):284–309, 1999.[9] Pierre Degond, Fabrice Deluzet, Laurent Navoret, An-Bang Sun, and Marie-H´el`ene Vignal. Asymptotic-preserving particle-in-cell method for the vlasov–poisson system near quasineutrality.
Journal of Com-putational Physics , 229(16):5630–5652, 2010.[10] Pierre Degond and Fabrice Deluzet. Asymptotic-preserving methods and multiscale models for plasmaphysics.
Journal of Computational Physics , 336:429–457, 2017.[11] National Institute of Standards and Technology.
The NIST Reference on Constants, Units, and Un-certainty , 5 2019.[12] Dongbin Xiu.
Numerical methods for stochastic computations: a spectral method approach . Princetonuniversity press, 2010.[13] Claudio Canuto and Alfio Quarteroni. Approximation results for orthogonal polynomials in sobolevspaces.
Mathematics of Computation , 38(157):67–86, 1982.[14] Dongbin Xiu and Jan S Hesthaven. High-order collocation methods for differential equations withrandom inputs.
SIAM Journal on Scientific Computing , 27(3):1118–1139, 2005.[15] Jingwei Hu, Shi Jin, and Ruiwen Shu. A stochastic galerkin method for the fokker–planck–landauequation with random uncertainties. In
XVI International Conference on Hyperbolic Problems: Theory,Numerics, Applications , pages 1–19. Springer, 2016.[16] Jingwei Hu and Shi Jin. Uncertainty quantification for kinetic equations. In
Uncertainty Quantificationfor Hyperbolic and Kinetic Equations , pages 193–229. Springer, 2017.[17] Shi Jin and Yuhua Zhu. Hypocoercivity and uniform regularity for the vlasov–poisson–fokker–plancksystem with uncertainty and multiple scales.
SIAM Journal on Mathematical Analysis , 50(2):1790–1816, 2018.[18] Zhiyan Ding and Shi Jin. Random regularity of a nonlinear landau damping solution for the vlasov-poisson equations with random inputs.
International Journal for Uncertainty Quantification , 9(2),2019.[19] Edward G Phillips and Howard C Elman. A stochastic approach to uncertainty in the equations ofmhd kinematics.
Journal of Computational Physics , 284:334–350, 2015.[20] Kazuo Yamazaki. Stochastic hall-magneto-hydrodynamics system in three and two and a half dimen-sions.
Journal of Statistical Physics , 166(2):368–397, 2017.[21] Tianbai Xiao and Martin Frank. A stochastic kinetic scheme for multi-scale flow transport with uncer-tainty quantification. arXiv preprint arXiv:2002.00277 , 2020.[22] Pierre Andries, Kazuo Aoki, and Benoit Perthame. A consistent bgk-type model for gas mixtures.
Journal of Statistical Physics , 106(5):993–1018, 2002.
23] TF Morse. Energy and momentum exchange between nonequipartition gases.
The Physics of Fluids ,6(10):1420–1427, 1963.[24] SI Braginskii. Transport processes in a plasma.
Reviews of plasma physics , 1, 1965.[25] C-D Munz, Pascal Omnes, Rudolf Schneider, Eric Sonnendr¨ucker, and Ursula Voss. Divergence correc-tion techniques for maxwell solvers based on a hyperbolic model.
Journal of Computational Physics ,161(2):484–511, 2000.[26] Mikhail N Kogan.
Rarefied Gas Dynamics . Plenum Press, 1969.[27] Kun Xu and Juan-Chen Huang. A unified gas-kinetic scheme for continuum and rarefied flows.
Journalof Computational Physics , 229(20):7747–7764, 2010.[28] Tianbai Xiao, Qingdong Cai, and Kun Xu. A well-balanced unified gas-kinetic scheme for multiscaleflow transport under gravitational field.
Journal of Computational Physics , 332:475 – 491, 2017.[29] Chang Liu and Kun Xu. A unified gas kinetic scheme for continuum and rarefied flows v: Multiscaleand multi-component plasma transport.
Communications in Computational Physics , 22(5):1175–1223,2017.[30] Tianbai Xiao, Kun Xu, and Qingdong Cai. A unified gas-kinetic scheme for multiscale and multicom-ponent flow transport.
Applied Mathematics and Mechanics , pages 1–18, 2019.[31] Tianbai Xiao, Chang Liu, Kun Xu, and Qingdong Cai. A velocity-space adaptive unified gas kineticscheme for continuum and rarefied flows.
Journal of Computational Physics , page 109535, 2020.[32] Ammar Hakim, John Loverich, and Uri Shumlak. A high resolution wave propagation scheme for idealtwo-fluid plasma equations.
Journal of Computational Physics , 219(1):418–442, 2006.[33] Lawrence F Shampine. Implementation of rosenbrock methods.
ACM Transactions on MathematicalSoftware (TOMS) , 8(2):93–113, 1982.[34] L. Landau. On the vibration of the electronic plasma.
Journal of Physics , 10(1):25 – 34, 1946.
10 20 30 40-10-8-6-4 t E x pe c t a t i on CurrentCollocationLinear theory (a) Expectation value t S t anda r d de v i a t i on CurrentCollocation (b) Standard deviationFigure 3: Expectation value and standard deviation of electric field energy (logarithmic) in linear landaudamping. (a) Expectation value (b) Standard deviationFigure 4: Time evolved expectation value and standard deviation of particle distribution function at x = 0in linear landau damping. able 5: Nomenclature of stochastic kinetic scheme. k B Boltzmann constant ε vacuum permittivity µ vacuum permeability e elementary charge c speed of light q α charge of species α , with q α = ± e ( α = ion, ele ) m α particle mass of species αn α number density of species αρ α density of species α U α macroscopic velocity of species αT α temperature of species α P α stress tensor of species αp α pressure of species α q α heat flux of species α W α macroscopic conservative variables of species α (density, momentum, energy) f α particle distribution function of species α u particle velocity a α electromagnetic force acting on species αQ α kinetic collision operator of species α M α equilibrium distribution function of species α E electric field B magnetic field φ correction potential for E in perfectly hyperbolic Maxwell’s equations ψ correction potential for B in perfectly hyperbolic Maxwell’s equations σ charge density J current density M abbreviation of electromagnetic variables ( E , B , φ, ψ ) (cid:36) vector of collision invariants ν α collision frequency of species αν ie interaction frequency between ion and electron τ α Collision time of species α with τ α = 1 /ν α F Wα flux for macroscopic conservative variables of species α F fα flux for particle distribution function of species α F M flux for electromagnetic variables L N generalized polynomial chaos expansion of any stochastic variable L with degree N ˆ l vector of generalized polynomial chaos coefficients of any stochastic variable L Φ orthogonal polynomials in random space z (cid:37) probability density function of random variable z Kn Knudsen number r L gyroradius λ D Debye length H [ x ] Heaviside step function 29
10 20 30 40 50-8-6-4-20 t E x pe c t a t i on CurrentCollocation (a) Expectation value t S t anda r d de v i a t i on CurrentCollocation (b) Standard deviationFigure 5: Expectation value and standard deviation of electric field energy (logarithmic) in nonlinear landaudamping. -5.0 -2.5 0.0 2.5 5.0-4-2024 x u (a) t = 0 -5.0 -2.5 0.0 2.5 5.0-4-2024 x u (b) t = 5 -5.0 -2.5 0.0 2.5 5.0-4-2024 x u (c) t = 10 -5.0 -2.5 0.0 2.5 5.0-4-2024 x u (d) t = 20 -5.0 -2.5 0.0 2.5 5.0-4-2024 x u (e) t = 30 -5.0 -2.5 0.0 2.5 5.0-4-2024 x u (f) t = 50Figure 6: Time evolved expectation value of particle distribution function over phase space ( x, u ) in nonlinearlandau damping. a) Expectation value (b) Standard deviationFigure 7: Time evolved expectation value and standard deviation of particle distribution function at x = 0in nonlinear landau damping. t E x pe c t a t i on CurrentCollocationLinear theory (a) Expectation value t S t anda r d de v i a t i on CurrentCollocation (b) Standard deviationFigure 8: Expectation value and standard deviation of electric field energy (logarithmic) in linear two-streaminstability. x u (a) Expectation value x u (b) Standard deviationFigure 9: Expectation value and standard deviation of particle distribution function over phase space ( x, u )in linear two-stream instability. x u (a) Expectation value
10 20 30 40-4-2024 x u (b) Standard deviationFigure 10: Expectation value and standard deviation of particle distribution function over phase space ( x, u )in nonlinear two-stream instability. .00 0.25 0.50 0.75 1.000.20.40.60.81.0 x E x pe c t a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation)N (MHD) (a) E ( N ) x E x pe c t a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation)U (MHD) (b) E ( U ) x E x pe c t a t i on By (Current)By (Collocation)By (MHD) (c) E ( B y ) x E x pe c t a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation) (d) E ( N ) x E x pe c t a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation) (e) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (f) E ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (g) E ( N ) x E x pe c t a t i on U (Current)U (Collocation) (h) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (i) E ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (j) E ( N ) x E x pe c t a t i on U (Current)U (Collocation) (k) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (l) E ( B y ) x E x pe c t a t i on N (Current)N (Collocation)N (Neutral) (m) E ( N ) x E x pe c t a t i on U (Current)U (Collocation)U (Neutral) (n) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (o) E ( B y )Figure 11: Expectation values of N , U and B y in Brio-Wu shock tube with density uncertainty at t = 0 . r g = 0 . r g = 0 .
01, row 3: r g = 0 .
1, row 4: r g = 1, row 5: r g = 100).= 100).
1, row 4: r g = 1, row 5: r g = 100).= 100). .00 0.25 0.50 0.75 1.000.000.010.020.03 x S t anda r d de v i a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation)N (MHD) (a) S ( N ) x S t anda r d de v i a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation)U (MHD) (b) S ( U ) x S t anda r d de v i a t i on By (Current)By (Collocation)By (MHD) (c) S ( B y ) x E x pe c t a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation) (d) S ( N ) x E x pe c t a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation) (e) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (f) S ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (g) S ( N ) x E x pe c t a t i on U (Current)U (Collocation) (h) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (i) S ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (j) S ( N ) x E x pe c t a t i on U (Current)U (Collocation) (k) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (l) S ( B y ) x E x pe c t a t i on N (Current)N (Collocation)N (Neutral) (m) S ( N ) x E x pe c t a t i on U (Current)U (Collocation)U (Neutral) (n) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (o) S ( B y )Figure 12: Standard deviations of N , U and B y in Brio-Wu shock tube with density uncertainty at t = 0 . r g = 0 . r g = 0 .
01, row 3: r g = 0 .
1, row 4: r g = 1, row 5: r g = 100).= 100).
1, row 4: r g = 1, row 5: r g = 100).= 100). .10.20.30.40.5 (a) E ( h ion ) (b) S ( h ion ) (c) E ( h ion ) (d) S ( h ion ) (e) E ( h ion ) (f) S ( h ion ) (g) E ( h ion ) (h) S ( h ion )Figure 13: Expectations and variances of reduced ion distribution h ion , in Brio-Wu shock tube with densityuncertainty at t = 0 . r g = 0 . r g = 0 .
01, row 3: r g = 0 .
1, row 4: r g = 1). .00 0.25 0.50 0.75 1.000.20.40.60.81.0 x E x pe c t a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation)N (MHD) (a) E ( N ) x E x pe c t a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation)U (MHD) (b) E ( U ) x E x pe c t a t i on By (Current)By (Collocation)By (MHD) (c) E ( B y ) x E x pe c t a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation) (d) E ( N ) x E x pe c t a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation) (e) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (f) E ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (g) E ( N ) x E x pe c t a t i on U (Current)U (Collocation) (h) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (i) E ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (j) E ( N ) x E x pe c t a t i on U (Current)U (Collocation) (k) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (l) E ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (m) E ( N ) x E x pe c t a t i on U (Current)U (Collocation) (n) E ( U ) x E x pe c t a t i on By (Current)By (Collocation) (o) E ( B y )Figure 14: Expectation values of N , U and B y in Brio-Wu shock tube with magnetic uncertainty at t = 0 . r g = 0 . r g = 0 .
01, row 3: r g = 0 .
1, row 4: r g = 1, row 5: r g = 100).= 100).
1, row 4: r g = 1, row 5: r g = 100).= 100). .00 0.25 0.50 0.75 1.000.000.010.020.030.040.050.06 x S t anda r d de v i a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation)N (MHD) (a) S ( N ) x S t anda r d de v i a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation)U (MHD) (b) S ( U ) x S t anda r d de v i a t i on By (Current)By (Collocation)By (MHD) (c) S ( B y ) x E x pe c t a t i on Ni (Current)Ne (Current)Ni (Collocation)Ne (Collocation) (d) S ( N ) x E x pe c t a t i on Ui (Current)Ue (Current)Ui (Collocation)Ue (Collocation) (e) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (f) S ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (g) S ( N ) x E x pe c t a t i on U (Current)U (Collocation) (h) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (i) S ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (j) S ( N ) x E x pe c t a t i on U (Current)U (Collocation) (k) S ( U ) x E x pe c t a t i on By (Current)By (Collocation) (l) S ( B y ) x E x pe c t a t i on N (Current)N (Collocation) (m) S ( N ) x E x pe c t a t i on U (Current)U (Collocation) (n) S ( U ) -2 -2 -2 -2 -2 -2 x E x pe c t a t i on By (Current)By (Collocation) (o) S ( B y )Figure 15: Standard deviations of N , U and B y in Brio-Wu shock tube with magnetic uncertainty at t = 0 . r g = 0 . r g = 0 .
01, row 3: r g = 0 .
1, row 4: r g = 1, row 5: r g = 100).= 100).