Asha Kumari Meena

Harish Kumar
email: hkumar@maths.iitd.ac.in

Robust MUSCL Schemes for Ten-Moment Gaussian Closure Equations with Source Terms

Keywords: Ten-Moment equations, Finite Volume Methods, MUSCL Scheme, Positivity preserving schemes

In this article, we present positivity preserving second-order numerical schemes to approximate Ten-Moment Gaussian closure equations with source terms. The challenge here is to preserve the positivity of the density and the symmetric pressure tensor. We propose MUSCL type numerical schemes to overcome these difficulties. The principal components of the proposed schemes are a Strang splitting of the source terms, positivity preserving first order scheme and suitable linear reconstruction process which ensures the positivity of the reconstructed variables. To achieve positivity of reconstructed variables, we impose the additional restrictions on the slopes of the linear reconstructions. Additionally, the source is discretized using both explicit and implicit methods. In the case of explicit source discretization, we derive the appropriate condition on the time step for discretization to be positivity preserving. Implicit discretization of the source terms is shown to be unconditionally positivity preserving. Numerical examples are presented to demonstrate the superior robustness and stability of the proposed numerical schemes.

Introduction

Fluid components of the plasma flows are often modeled by Euler equations of compressible flows. These equations are derived by taking moments of Boltzmann equation with respect to the velocity. The resulting set of equations is then closed by assuming local thermodynamic equilibrium. However, for many applications (especially related to the plasma flows,( See [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF][START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF][START_REF] Brown | Numerical Solution of a 10-Moment Model for Nonequilibrium Gasdynamics[END_REF][START_REF] Dubroca | Magnetic field generation in plasmas due to anisotropic laser heating[END_REF][START_REF] Hakim | Extended MHD Modelling with Ten-Moment Equations[END_REF][START_REF] Johnson | Gaussian-Moment Relaxation Closures for Verifiable Numerical Simulation of Fast Magnetic Reconnection in Plasma[END_REF][START_REF] Morreeuw | Electron temperature anisotropy modeling and its effect on anisotropy-magnetic field coupling in an underdense laser heated plasma[END_REF][START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF][START_REF] Sangam | Anisotropic instability in a laser heated plasma[END_REF]), the local thermodynamic equilibrium assumption is not valid, and one need to take the anisotropic nature of the pressure into account. To achieve this, Levermore et al. proposed Ten-Moment equations model (See [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF][START_REF] Levermore | The Gaussian Moment Closure for Gas Dynamics[END_REF]), which is derived by the Gaussian closure of the kinetic model. This results in a hyperbolic system of conservation laws where the pressure is described using symmetric tensor.

Due to hyperbolicity of the system, we need to consider the weak solutions, as solutions may develop discontinuities like shock and contact waves. However, the class of weak solution is too large. An additional criterion in the form of entropy inequality is imposed to choose the physically relevant solution.

Numerical methods to discretize hyperbolic conservation laws are often based on finite volume schemes (or finite difference schemes) (See [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF], [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluids dynamics, A Practical Introduction[END_REF]). The initial solution is discretized using cell averages, which are then evolved by using numerical flux over the cell edges. The numerical flux is either based on the exact solution of the local Riemann problem or some approximation of it. The higher-order schemes are designed by reconstructing the solutions inside the cells using ENO, WENO or TVD limiters based procedures. For the time update, SSP-RK schemes are used.

A primary concern for any numerical scheme for hyperbolic conservation laws is to be robust, i.e., given the cell averages in the physically admissible set, we would like numerical schemes to produce updated cell averages in the same set. It is essential; otherwise, we may lose the hyperbolicity of the system. For the Euler equations of compressible flows, this is equivalent of preserving the positivity of the density and the pressure. There are several works on design of positivity preserving schemes for Euler equations (See [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF], [START_REF] Berthon | Robustness of MUSCL schemes for 2D unstructured meshes[END_REF], [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF], [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] and [START_REF] Zhang | On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF]).

For the Ten-Moment equations, the robustness translates into ensuring positivity of the density and the symmetric pressure tensor (unlike Euler equations, where the pressure is a scalar). In the case of first-order schemes, this is based on the construction of the suitable numerical flux. For Ten-Moment equations (1), a relaxation based scheme is proposed in [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF]. In addition to ensuring positivity, the scheme is also shown to be entropy stable. In [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF], a Ten-Moment equation based plasma flow model is considered, which is discretized using HLLC numerical flux solver and shown to be positivity preserving. For higher order schemes, a wave propagation based discretization is proposed for a plasma flow model based on Ten-Moment equation in [START_REF] Hakim | Extended MHD Modelling with Ten-Moment Equations[END_REF]. However, the scheme does not guarantee positivity. More recently, in [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF] a Ten-Moment based plasma flow model with source terms is discretized to simulate laser effects on matter. The discretization is based on an equivalent relaxation model, which also takes source terms into consideration. This results in a first-order scheme, which is shown to be positivity preserving and entropy stable.

In this article, we consider the Ten-Moment equations with source terms model considered in [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF]. We propose second-order discretizations which ensure the positiv-ity of the density and the pressure tensor. We achieve this as follows:

• For the numerical flux, we use positivity preserving approximated Riemann solvers based numerical fluxes. There are several examples of them, like Lax-Friedrichs, HLLE and HLLC solvers (See [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF]).

• To obtain second-order positivity preserving schemes, we follow [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF] and [START_REF] Waagan | A positive MUSCL-Hancock scheme for ideal magnetohydrodynamics[END_REF], and propose robust MUSCL (See [START_REF] Harten | On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws[END_REF]) reconstruction processes. We propose two slope limiters based on the reconstruction of the primitive variables, namely: Generalized slope limiter and Conservative slope limiters. We prescribe exact conditions for both the cases, for the schemes to be positivity preserving.

• The source is discretized using both explicit and implicit schemes. In the case of explicit discretization, a condition on time step is derived, which ensures positivity of the solution. The implicit treatment of the source is shown to be unconditionally positivity preserving. Furthermore, we demonstrate that we do not need to solve any system of algebraic equations to implement implicit source update.

• Both source and flux discretizations are then combined using Strang splitting, which ensures positivity of the whole scheme.

The rest of the article is organized as follows: In the following Section, we first present the HLLC solvers for Ten-Moment model without source terms, which is simplified from [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF]. In Section 4, we present the MUSCL based reconstruction process. In Section 5, we present the analysis of the source discretization, followed by time discretization in Section 6. Numerical simulations to demonstrate the superior robustness of the proposed schemes are presented in Section 7.

Ten-Moment Equations with Source Terms

Following [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF], we consider the Ten-Moment equations with source terms which model inhomogeneous heating of the electrons in dense plasmas using lasers. In two dimensions, these equations can be written as

∂ t ρ + ∇ • (ρv) = 0, ( 1a
)
∂ t (ρv) + ∇ • (ρv ⊗ v + p) = - 1 2 ρ∇W, (1b)
∂ t E + ∇ • ((E + p) ⊗ v) = - 1 4 ρ(∇W ⊗ v + v ⊗ ∇W ). (1c)
Here, ρ is the density, v = (v 1 , v 2 ) is the velocity vector and E is the symmetric energy tensor with components E 11 , E 12 and E 22 . The set of equations is closed by the equation of state,

E = 1 2 (p + ρv ⊗ v) , (2)
where p is the symmetric pressure tensor with components p 11 , p 12 and p 22 . The given function W (x, y, t) is the electron quiver energy due to laser light. The equation (1a) represent balance of mass, followed by equation (1b) for momentum balance and (1c) is the balance of energy tensor. The system (1) can also be written as,

∂ t u + ∂ x f x (u) + ∂ y f y (u) = s(u), (3)
with conservative state variable u = {ρ, ρv 1 , ρv 2 , E 11 , E 12 , E 22 } . The flux components are given by,

f x (u) =                ρv 1 ρv 2 1 + p 11 ρv 1 v 2 + p 12 (E 11 + p 11 )v 1 E 12 v 1 + 1 2 (p 11 v 2 + p 12 v 1 ) E 22 v 1 + p 12 v 2                , (4a)
and

f y (u) =                ρv 2 ρv 1 v 2 + p 12 ρv 2 2 + p 22 E 11 v 2 + p 12 v 1 E 12 v 2 + 1 2 (p 12 v 2 + p 22 v 1 ) (E 22 + p 22 )v 2                . ( 4b
)
The source terms can be written component wise as follows:

s(u) =                0 -1 2 ρ∂ x W 0 -1 2 ρv 1 ∂ x W -1 4 ρv 2 ∂ x W 0                +                0 0 -1 2 ρ∂ y W 0 -1 4 ρv 1 ∂ y W -1 2 ρv 2 ∂ y W                . (5)
For solutions to be physically meaningful, state variable u needs to be in the convex set of physically admissible states,

Ω = u ∈ R 6 | ρ > 0, and x p x > 0, ∀ x ∈ R 2 /{0} , (6)
i.e the density and the pressure tensor has to be positive. For the solution in Ω, we have the following result from [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF]:

Lemma 2.1 The system (1) without source terms, is hyperbolic for u ∈ Ω and admits the eigenvalues,

v.n, v.n ± 3(p.n).n ρ , v.n ± (p.n).n ρ ,
along the unitary vector n. The eigenvalue v.n has two order of multiplicity while other eigenvalues have one order of multiplicity. The eigenvalue v.n is associated to a linearly degenerate field. The eigenvalues v.n ± 3(p.n).n ρ are associated to a genuinely nonlinear field while eigenvalues v.n± (p.n).n ρ are associated to a linearly degenerate field.

In addition, we observe that the source terms are Ω-invariant, i.e.: Proof To simplify, we consider the case of W as a function of x and t only. Dependance on y can be considered in similar way. Then we have

dρ dt = 0, (7)
for the density component. For the momentum, using [START_REF] Dubroca | Magnetic field generation in plasmas due to anisotropic laser heating[END_REF] we get,

dρv 1 dt = ρ dv 1 dt = - 1 2 ρ∂ x W, which implies dv 1 dt = - 1 2 ∂ x W and dv 2 dt = 0. ( 8
)
Considering energy tensor, we get,

d dt p 11 + ρv 2 1 = -ρv 1 ∂ x W, d dt (p 12 + ρv 1 v 2 ) = - 1 2 ρv 1 ∂ x W, d dt p 22 + ρv 2 2 = 0,
So, for p 11 using (8),

dp 11 dt = -2ρv 1 dv 1 dt -ρv 1 ∂ x W = -2ρv 1 dv 1 dt + 1 2 ∂ x W = 0.
Similarly, we can also show that,

dp
So, the evolutions of density and the pressure tensor are not affected by source term. Hence source is Ω-invariant.

For the system of Ten-Moments equations (1), following [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF], we introduce entropy e and entropy flux q as follows:

s(p, ρ) = ln det p ρ 4
, e = -ρs and q = ev.

Then we can deduce the following entropy stability result from [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF].

International Journal on Finite Volumes Lemma 2.3 The weak solutions of (1) satisfies following entropy stability condition:

∂ t e + ∇ • q ≤ 0. ( 11
)
In addition, the model above can be shown to be rotational invariant (See [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF]). As we have seen that at the continuous level, solutions are in Ω, it is important to design numerical schemes which are Ω-invariant.

3 First Order Schemes for Homogeneous Model

We will consider the homogeneous equations in one dimensions to simplify the discussion i.e. we will consider the following equations:

∂ t u + ∂ x f x (u) = 0, (12)
instead of (3). The extension of the numerical schemes proposed here to the two dimensions is straight forward. Let us consider a computational mesh, with (x i ) i∈Z as the cell centres of the cells

I i = (x i-1 2 , x i+ 1 2 ) with x i+ 1 2 = x i + (x i+1 -x i )/2.
We will consider uniform mesh for simplicity i.e. we assume ∆x = x i+1 -x i for i ∈ Z. The cell averages of the state variable u at the time t n are defined as,

w n i = I i u(x, t n )dx. ( 13
)
We will assume that the cell averages (w n i ) i∈Z belongs to set of physically admissible states i.e. w n i ∈ Ω for all i ∈ Z. Our aim is to propose numerical schemes which ensures that the evolved cell average w n+1 i ∈ Ω. We define it formally as follows: Definition 3.1 Ω-invariance: An update U of the solution (w n i ) ∈ Ω, for all i ∈ Z, is called Ω-invariant (or positivity preserving), if updated solution

w n+1 i = U(w n i ) is also in Ω.
A first order finite volume scheme for the discretization of ( 12) can be written in the form,

w n+1 i = w n i - ∆t ∆x F x (w n i , w n i+1 ) -F x (w n i-1 , w n i ) , (14)
where F x is the numerical flux, which is conservative and consistent with continuous flux f x . For the first order schemes, the Ω-invariance of the scheme depends on the choice of numerical flux. Several numerical fluxes exists for the 10-moment equations, which are Ω invariant, e.g. HLLE, Rusanov, HLLC (See [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF]), relaxation solver (See [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF]). Here, we will now present HLLC solvers for Ten-Moment equations, which is simplified from the HLLC solver presented in (See [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF]) for a more general system.

HLLC Numerical Flux for Ten-Moment Equations

The HLLC numerical flux is based on the approximated Riemann solver consisting of two intermediate states, namely, w * l and w * r i.e we are looking at the solutions of the form:

w hllc =        w l , if s l ≥ x t w * l , if s l < x t ≤ v * 1 w * r , if v * 1 < x t ≤ s r w r , if s r ≤ x t . (15)
Here s l , s r are approximations of the fastest left and right wave, respectively (see [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluids dynamics, A Practical Introduction[END_REF]) and intermediate wave speed is v * 1 . The corresponding numerical flux is,

F hllc =        F l , if s l ≥ x t F * l , if s l < x t < v * 1 F * r , if v * 1 < x t < s r F r , if s r ≤ x t . (16)
Assuming that the expressions for s l and s r is known, we can determine states w * l and w * r , using Rankine-Hugoniot (RH) conditions across the three waves. In addition, following [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF], we also impose following conditions across the contacts:

v * 1 l = v * 1r = v * 1 , v * 2 l = v * 2r = v * 2 , (17a)
p * 11 l = p * 11r = p * 11 , p * 12 l = p * 12r = p * 12 . (17b)
Then we have the expressions for velocities as,

v * 1 = p 11 l -p 11r + ρ l v 1 l (v 1 l -s l ) -ρ r v 1r (v 1r -s r ) ρ l (v 1 l -s l ) -ρ r (v 1r -s r ) , (18a)
v * 2 = p 12 l -p 12r + ρ l v 2 l (v 1 l -s l ) -ρ r v 2r (v 1r -s r ) ρ l (v 1 l -s l ) -ρ r (v 1r -s r ) . (18b)
Expression for the densities are then derived using RH condition across s l and s r waves and we get,

ρ * l = ρ l (v 1 l -s l ) v * 1 -s l , and ρ * r = ρ r (v 1r -s r ) v * 1 -s r . ( 19
)
Now the pressures can be derived as follows:

p * 11 = p 11 l + ρ * l v * 1 (s l -v * 1 ) -ρ l v 1 l (s l -v 1 l ) (20a) = p 11r + ρ * r v * 1 (s r -v * 1 ) -ρ r v 1r (s r -v 1r ), p * 12 = p 12 l + ρ * l v * 2 (s l -v * 1 ) -ρ l v 2 l (s l -v 1 l ) (20b) = p 12r + ρ * r v * 2 (s r -v * 1 ) -ρ r v 2r (s r -v 1r ).
Now, we can derive expressions for the energy components. For left states we get,

E * 11 l = v 1 l p 11 l -v * 1 p * 11 + E 11 l (v1 l -s l ) v * 1 -s l , (21a)
Robust MUSCL Schemes for Ten-Moment Gaussian Closure Equations with Source Terms

E * 12 l = v 1 l p 12 l + v 2 l p 11 l -v * 1 p * 12 -v * 2 p * 11 2(v * 1 -s l ) + E 12 l (v 1 l -s l ) v * 1 -s l , (21b)
E * 22 l = v 1 l p 12 l -v * 1 p * 12 + E 22 l (v 1 l -s l ) v * 1 -s l . (21c)
Similarly, expression for E * 11r , E * 12r and E * 22r can also be derived. To prove positivity of the solver we need conditions on s l and s r . We have the following result from [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF].

Proposition 3.2 The HLLC Riemann solver for Ten-Moment equations is positivity preserving i.e. w hllc ∈ Ω for {w l , w r } ∈ Ω, if the fastest left and right wave speeds s l and s r , satisfies the following bounds:

s l ≤ min(v 1 l -ĉl , v 1 l -cl ), and s r ≥ max(v 1r + ĉr , v 1r + cr ), (22)
where,

ĉl = p 11 l ρ l , cl = p 12 l √ ρ l p 22 l , ĉr = p 11r ρ r
, and cr = p 12r √ ρ r p 22r .

There are several choices of wavespeeds s l and s r are possible (See [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluids dynamics, A Practical Introduction[END_REF]), which satisfies these bounds. We follow [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF] and consider following speeds:

s l = min(λ Roe min , v 1 l -c l , v 1 l -cl ), (23a)
s r = max(λ Roe max , v 1r + c r , v 1r + cr ), (23b)
where

c l = 3p 11 l ρ l , c r = 3p 11r ρ r , cl = p 12 l √ ρ l p 22 l
, and cr = p 12r √ ρ r p 22r .

and λ Roe min , λ Roe max are the left and right speed for Roe averaged state, i.e.:

λ Roe min = v1 -c, and

λ Roe max = v1 + c, with, c = H 11 l √ ρ l + H 11r √ ρ r √ ρ l + √ ρ r -v1 2 , v1 = v 1 l √ ρ l + v 1r √ ρ r √ ρ l + √ ρ r , and
H 11 l = v 2 1l + 3p 11 l ρ l , H 11r = v 2 1r + 3p 11r ρ r .

Second Order Schemes for Homogeneous Model

To achieve the second order of accuracy in space, we use MUSCL procedure. This is based on the linear reconstruction using the cell averages. Consider a scalar function u and its cell averages w i over the computational cells I i . A piecewise linear reconstruction based the cell averages is defined as,

p(x) = w i + Dw i (x -x i ), x ∈ I j , (24)
where the slope Dw i is defined using a TVD limiter (See [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluids dynamics, A Practical Introduction[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]). One such commonly used limiter is MinMod limiter given by,

Dw i = minmod w i -w i-1 ∆x , w i+1 -w i ∆x , (25)
where

minmod(a, b) = max(0, min(a, b)) + min(0, max(a, b)).
The traces on the reconstructed function p(x) on the edges of the cell I i are denoted as follows:

w n,± i = w n i ± ∆w i , ∆w i = ∆x 2 Dw i . ( 26
)
This reconstruction process is performed on primitive or conservative variables, component wise. In this work, we shall restrict ourself to reconstruction of the primitive variables, as this is the most commonly used reconstruction process in MUSCL scheme. However, note that this is not equal to direct reconstruction of the conservative variables. Let us assume that w n,± i is the reconstructed conservative variable obtained by reconstructed conservative variables component wise and converting them to conservative variable. The reconstruction process is called conservative if,

w n,+ i + w n,- i = 2w n i , (27)
hold. Note that this is not true for every reconstruction. Following the reconstruction, a second order scheme can be written as,

w n+1 i = w n i - ∆t ∆x F x (w n,+ i , w n,- i+1 ) -F x (w n,+ i-1 , w n,- i ) , (28)
The scheme is second order accurate in space, however there is no guarantee that it is Ω-invariant. We will now describe the process which ensures that. Following [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF], we introduce the intermediate state w n, * i , which satisfies,

w n i = αw n,- i + (1 -2α)w n, * i + αw n,+ i , (29)
for some α ∈ (0, 1 3 ]. Then we can rewrite (28) as follows:

w n+1 i = w n i - ∆t ∆x F (w n,+ i , w n,- i+1 ) -F (w n,+ i-1 , w n,- i ) = αw n,- i + (1 -2α)w n, * i + αw n,+ i - ∆t ∆x F (w n,+ i , w n,- i+1 ) -F (w n,+ i-1 , w n,- i ) = αw n,- i + (1 -2α)w n, * i + αw n,+ i - ∆t ∆x (F (w n,- i , w n, * i ) -F (w n,+ i-1 , w n,- i )) -(F (w n, * i , w n,+ i ) -F (w n,- i , w n, * i )) -(F (w n,+ i , w n,- i+1 ) -F (w n, * i , w n,+ i )) = α w n,- i - ∆t α∆x F (w n,- i , w n, * i ) -F (w n,+ i-1 , w n,- i ) +(1 -2α) w n, * i - ∆t (1 -2α)∆x F (w n, * i , w n,+ i ) -F (w n,- i , w n, * i ) +α w n,+ i - ∆t α∆x F (w n,+ i , w n,- i+1 ) -F (w n, * i , w n,+ i ) . ( 30
)
Robust MUSCL Schemes for Ten-Moment Gaussian Closure Equations with Source Terms

We observe that w n+1 i is a convex combination of first-order schemes with time steps ∆t α and ∆t (1-2α) . This results in following (See [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF]):

Theorem 4.1 Consider a first-order Ω-invariant (positivity preserving) scheme under a standard CFL condition with CFL number C and α ∈ (0, 1 3 ]. Then the MUSCL scheme (28) is Ω-invariant with CFL number αC if the following conditions hold:

a) w n i ∈ Ω, for all i ∈ Z. b) w n,± * i ∈ Ω, for all i ∈ Z.
The standard reconstruction process described above doesn't ensure second condition of the Theorem 4.1. We will now present two reconstruction processes which ensures these conditions. First limiter is a general slope limiter and second limiter is conservative slope limiter. Both are based on the reconstruction of the primitive variables.

General Slope limiting of Primitive Variable (GSPV)

Consider the primitive variables ŵ = (ρ, v 1 , v 2 , p 11 , p 12 , p 22 ) . We define cell edge values ŵn,± i , in cell I i as,

ρ n,± i = ρ n i ± ∆ρ i , v n,± 1,i = v n 1,i ± ∆v 1,i , v n,± 2,i = v n 2,i ± ∆v 2,i , p n,± 11,i = p n 11,i ± ∆p 11,i , p n,± 12,i = p n 12,i ± ∆p 12,i , p n,± 22,i = p n 22,i ± ∆p 22,i . (31)
To ensure the Ω-invariant of the scheme, we aim that the conservative variable w n,± i corresponding to ŵn,± i and the state w * are in Ω. This leads to the additional conditions on the slopes. We proceed as follows:

Positivity of ρ n,± i , p n,±
11,i and p n,± 22,i can be easily established if,

|∆ρ i | < ρ n i , |∆p 11,i | < p n 11,i , |∆p 22,i | < p 22,i . (32)
To ensure the positivity of the pressure tensor we need that the determinant of tensors, p n,± i is positive. This results in,

|p n 12,i + ∆p
We now consider the state w n, * i . We take α = 1 3 and observe that (See Equation (29)),

ρ n, * i = ρ n i . ( 34
)
Robust MUSCL Schemes for Ten-Moment Gaussian Closure Equations with Source Terms So, the density component of w n, * i is positive. The pressure components can be written as,

p n, * 11,i = p n 11,i -2 1 + 2 ∆ρ i ρ n i 2 ρ n i ∆v 2 1,i , (35a)
p n, * 12,i = p n 12,i -2 1 + 2 ∆ρ i ρ n i 2 ρ n i ∆v 1,i ∆v 2,i , (35b)
p n, * 22,i = p n 22,i -2 1 + 2 ∆ρ i ρ n i 2 ρ n i ∆v 2 2,i . (35c)
Also det(p n, * i ) can be written as,

p n, * 11,i p n, * 22,i -p n, * 2 12,i = p n 11,i p n 22,i -(p n 12,i ) 2 -2ρ n i 1 + 2 ∆ρ i ρ n i 2 p n 11,i ∆v 2 2,i + p n 22,i ∆v 2 1,i -2p n 12,i ∆v 1,i ∆v 2,i . (36)
The positivity of p n, * 11 and p n, * 22 is guaranteed with the following conditions on velocity slopes,

|∆v 1,i | < p n 11,i 2ρ n i 1 + 2 ∆ρ i ρ n I 2 , |∆v 2,i | < p n 22,i 2ρ n i 1 + 2 ∆ρ i ρ n i 2 . ( 37
)
To ensure the positivity of det(p n, * i ), we assume, ∆v Proposition 4.2 Assume that the initial data w n i is in the set Ω for all i ∈ Z. Then, the GSPV-reconstruction satisfies the sufficient conditions of the Theorem (4.1) with α = 1 3 if (32), (33), (37) and (38) hold. To implement conditions (32), (33), (37) and (38) numerically, we proceed as follows: First calculate the slopes of primitive variables using MinMod limiters. We can use any other limiter also. Then we update the slopes ∆ρ i , ∆p 11,i and ∆p 22,i to satisfy (32), by taking,

∆ρ i = max(-ρ n i , min(ρ n i , ∆ρ i )), (40a)
∆v 1,i = max(v l 1,i , min(v u 1,i , ∆v 1,i )), ( 42
) with v l 1,i = - p n 11,i 2ρ n i 1 + 2 ∆ρ i ρ n i 2 , v u 1,i = p n 11,i 2ρ n i 1 + 2 ∆ρ i ρ n i 2 ,
and ∆v

2,i = max(max(v l1 2,i , v l2 2,i ), min(min(v u1 2,i , v u2 2,i ), ∆v 2,i )), ( 43
) with v l1 2,i = - p n 22,i 2ρ n i 1 + 2 ∆ρ i ρ n i 2 , v u1 2,i = p n 22,i 2ρ n i 1 + 2 ∆ρ i ρ n i 2 , v l2
p n 11,i -2ρ n i ∆v 2 1,i 1 + 2 ∆ρ i ρ n i 2 2ρ n i p n 2 11,i 1 + 2 ∆ρ i ρ n i 2 .

Conservative Slope limiting of Primitive Variable (CSPV)

If the slopes are conservative we can rewrite Equation (30) as,

w n+1 i = 1 2 w n,- i - ∆t ∆x/2 F (w n,- i , w n,+ i ) -F (w n,+ i-1 , w n,- i ) + 1 2 w n,+ i - ∆t ∆x/2 F (w n,+ i , w n,- i+1 ) -F (w n,- i , w n,+ i ) .
This result in following result (See [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF]):

Theorem 4.
3 Consider a first-order domain invariant scheme under a standard CFL condition with CFL number C, then the MUSCL scheme is Ω-invariant with CFL number 1 2 C under the following conditions: a)

w n i ∈ Ω for all i ∈ Z, b) w n,± i ∈ Ω, for all i ∈ Z.
Let us consider the following cell edge values of primitive variables:

ρ n,± i = ρ n i ± ∆ρ i , v n,+ 1,i = v n 1,i + ρ n,- i ρ i ∆v 1,i , v n,- 1,i = v n 1,i - ρ n,+ i ρ i ∆v 1,i , v n,+ 2,i = v n 2,i + ρ n,- i ρ i ∆v 2,i , v n,- 2,i = v n 2,i - ρ n,+ i ρ i ∆v 2,i , p n,± 11,i = p n 11,i ± ∆p 11,i - ρ n,- i ρ n,+ i ρ i (∆v 1,i ) 2 , p n,± 12,i = p n 12,i ± ∆p 12,i - ρ n,- i ρ n,+ i ρ i ∆v 1,i ∆v 2,i , p n,± 22,i = p n 22,i ± ∆p 22,i - ρ n,- i ρ n,+ i ρ i (∆v 2,i ) 2 . ( 44
)
One can easily check that this reconstruction is conservative. Now we will derive the sufficient conditions for Theorem 4.3 by implementing restrictions on the slopes. For positivity of density we have,

|∆ρ i | < ρ n i . (45)
The positivity of p n,± 11,i can be achieved if,

|∆p 11,i | < p n 11,i , (∆v 1,i ) 2 < ρ n i (p n 11,i + ∆p 11,i ) ρ n,- i ρ n,+ i
, and (∆v

1,i ) 2 < ρ n i (p n 11,i -∆p 11,i ) ρ n,- i ρ n,+ i . (
ρ n,- i ρ n,+ i ρ n i ∆v 1,i ∆v 2,i < p n,+ 11,i p n,+ 22,i . (48b)
Robust MUSCL Schemes for Ten-Moment Gaussian Closure Equations with Source Terms Now we have the following result:

Proposition 4.4 Assume that the initial data w n i is in the set Ω for all i ∈ Z. Then, the CSPV-reconstruction satisfies the sufficient conditions of the Theorem (4.3) if the conditions (45)-(48) holds.

For the implementation of these conditions we proceeds as follows: First calculate the slopes of primitive variables using MinMod limiter. Again we can use any other limiter for reconstruction. We now modify the slope of density as

∆ρ i = max (-ρ n i , min ( ρ n i , ∆ρ i )) . (49)
Conditions on the slope of the velocity v 1,i is the following:

∆v 1,i = max(max(v l1 1,i , v l2 1,i ), min(min(v u1 1,i , v u2 1,i ), ∆v 1,i )), (52)
with,

v l1 1,i = - ρ n i (p n 11,i + ∆p 11,i ) ρ n,+ i ρ n,- i , and v u1 1,i = ρ n i (p n 11,i + ∆p 11,i ) ρ n,+ i ρ n,- i , v l2 1,i = - ρ n i (p n 11,i -∆p 11,i ) ρ n,+ i ρ n,- i , and v u2 1,i = ρ n i (p n 11,i -∆p 11,i ) ρ n,+ i ρ n,- i ,
Similarly, for v 2,i , we take,

∆v 2,i = max(max(v l1 2,i , v l2 2,i ), min(min(v u1 2,i , v u2 2,i ), ∆v 2,i )), (53)
with,

v l1 2,i = - ρ n i (p n 22,i + ∆p 22,i ) ρ n,+ i ρ n,- i , and v u1 2,i = ρ n i (p n 22,i + ∆p 22,i ) ρ n,+ i ρ n,- i , v l2
2,i = -

ρ n i (p n 22,i -∆p 22,i ) ρ n,+ i ρ n,- i ,

Discretization of the Source Terms

Let us consider the source ordinary differential equations,

dw i dt = s(w i ), with s(u) =         0 -1 2 ρ∂ x W 0 -1 2 ρv 1 ∂ x W -1 4 ρv 2 ∂ x W 0         . (55)
Here, we have considered the one dimensional case of the source term. Extension to higher dimensions is straight forward. Furthermore, let us assume that the given initial data w n i ∈ Ω. We want to discretize (55), so that w n+1 i ∈ Ω. We now present explicit and implicit Euler discretizations of (55).

Explicit Euler Source Update

An Euler explicit update of (55) is given by,

w n+1 i = w n i + ∆ts(w n i ). (56)
Let us denote this with w n+1 i = S e,1 ∆t (w n i ). As the source component corresponding to density is zero, we note that,

ρ n+1 i = ρ n i .
Similarly, we observe that v n+1

2,i = v n 2,i .
From the first momentum component we get,

ρ n+1 i v n+1 1,i = ρ n i v n 1,i - ∆t 2 ρ n i W n x,i ,
where W n x,i = ∂ x W (x i , t n ) for the given function W (x, t). This implies,

v n+1 1,i = v n 1,i - ∆t 2 W n x,i .
Considering energy tensor component E 11 , we observe that,

p n+1 11,i + ρ n+1 i (v n+1 1,i ) 2 = p n 11,i + ρ n i (v n 1,i ) 2 -∆tρ n i v n 1,i W n x,i , =⇒ p n+1 11,i -p n 11,i ρ n i = (v n 1,i ) 2 -(v n+1 1,i ) 2 -∆tv n 1,i W n x,i (57)
= - ∆t 2 2 W n x,i 2 .
So, finally we get,

p n+1 11,i = p n 11,i - ∆t 2 4 ρ n i W 2 x,i , which is positive if, ∆t ≤ p n 11,i ρ n i W n x,i 2 . ( 58
)
Similarly, we can show that,

p
ρ n i p n 22,i (W n x,i ) 2 . ( 59
)
Combining the discussion above, we have the following result:

Lemma 5.1 The explicit Euler source update S e,1 ∆t is Ω invariant under the time step restrictions (58) and (59).

Implicit Euler Source Update

We will now consider implicit discretization of the source terms. A first order Euler implicit scheme for the source ODE (55) is written by,

w n+1 i = w n i + ∆ts(w n+1 i ). ( 60
)
Let us denote this with w n+1 i = S i,1 ∆t (w n i ). Assuming w n i ∈ Ω, similar to the case of explicit discretization,

ρ n+1 i = ρ n i , v n+1
2,i = v n 2,i and p n+1 22,i = p n 22,i .

Note that we do not have to solve a system of equations to implement this update. We first update momentum by,

ρ n+1 i v n+1 1,i = ρ n i v n 1,i - ∆t 2 ρ n i W n+1 x,i .
which is then used to update energy tensor. Using this, from the momentum equation, we get,

v n+1 1,i = v n 1,i - ∆t 2 W n+1 x,i .
The evolution of energy component E 11 is,

p n+1 11,i + ρ n+1 i (v n+1 1,i ) 2 = p n 11,i + ρ n i (v n 1,i ) 2 -∆tρ n+1 i v n+1 1,i W n+1 x,i ,
Rearranging terms, we get,

p n+1 11,i -p n 11,i ρ n i = (v n 1,i ) 2 -(v n+1 1,i ) 2 -∆tv n+1 1,i W n+1 x,i = - ∆t 2 2 W n+1 x,i 2 + ∆t(v n 1,i -v n+1 1,i )W n+1 x,i = - ∆t 2 2 W n+1 x,i 2 + ∆t ∆t 2 W n+1 x,i 2 = ∆t 2 2 W n+1 x,i 2 .
So, finally we have,

p n+1 11,i = p n 11,i + ∆t 2 4 ρ n i W n+1 x,i 2 .
Hence, p n+1

Time Discretization

For the time discretization we use second order SSP-Runge Kutta methods (See [START_REF] Shu | TVD time discretizations[END_REF]). Let us denote flux update with second order SSP-RK scheme with H 2 ∆t , where each internal Euler update is given by (28). Similarly, let us denote S e,2 ∆t second order SSP-RK time update with each internal Euler update is given by (56) and S i,2

∆t second order SSP-RK time update with each internal update is given by (60).

Using Strang splitting, we propose two second order schemes:

• O2-exp w n+1 i = S e,2
w n i . (61)
• O2-imex

w n+1 i = S i,2 ∆t 2 H 2 ∆t S i,2 ∆t 2 w n i . ( 62
)
The source discretization consist of evaluation of derivatives W n x,i and W n+1 x,i of given function W . As in all our test cases the function W is smooth, we calculate the derivatives exactly and then evaluate them at the grid points. The exact solution of the Riemann problem consists of a shock wave and a rarefaction wave separated by contact discontinuity. So, numerical schemes are tested for the performance on all kinds of possible waves in solution.

Numerical results for the problem is presented in Figures 2. The solution is computed using 100 cells for MinMod, CSPV and GSPV slope limiter. We again observe that there is no difference in the performance of CSPV and GSPV limiters compare to the standard MinMod limiter. This holds for all the state variables and det p (See Figures 2(a)-2(f)). Also, all the waves are resolved by the three schemes. We note that the additional waves are present in v 2 and p 22 components (See Figure 2(c) and 2(e)). The simulation times for CSPV, GSPV and MinMod schemes were 2.8838, 2.8457 and 2.5258 seconds, respectively.

Two Shock Waves Riemann Problem

We consider the same domain as the previous Riemann problem with Riemann problem again centered at x = 0.0. The initial left and right states for the Two shock wave Riemann problem are given in Table 2. We assume outflow boundary conditions and solutions are computed till time t = 0.125. The exact solution of the problem is two shock waves moving away from each other and separated by a contact discontinuity. The numerical solutions are presented in Figures 3 which are computed using 100 cells. Solutions are compared for MinMod, CSPV and GSPV limiters. As we do not have any low density or low pressure areas, all the three schemes produce similar results and resolve all the waves with similar accuracy. Furthermore, additional waves are present in v 2 and p 22 components. The computational time of CSPV and GSPV schemes were 2.4524 and 2.3775 seconds, respectively. For the standard MinMod limiter it was 1.9167 seconds.

Two Rarefaction Waves Problem

We now consider the Riemann problem for which solution consists of two rarefaction waves separated by a contact. The domain and the point of the initial discontinuity are same as two above Riemann problems. The initial left and right states are given in Table 3. The solutions are computed till time t = 0.15 with outflow boundary conditions.

Computed solutions are presented in Figures 4. They are computed using 100 cells. Similar to previous Riemann problems we observe that the performance of the all the three schemes is similar. We also note that all the five waves are present in p 22 components. The computational time of the CSPV and GSPV schemes were 2.8194 and 2.4125 seconds. The MinMod scheme has the computational time of 2.3458 seconds.

Two Rarefaction Waves Problem with Near Vacuum State

To demonstrate the superior robustness of the presented schemes, we now consider a Riemann problem, for which the solutions contains low-density and low-pressure area (See [START_REF] Sangam | An HLLC scheme for Ten-Moments approximation coupled with magnetic field[END_REF]). The domain of the Riemann problem is the same as in the previous cases. The initial states are given in Table 4. The solution of the Riemann problem contains two Rarefaction waves moving from each other and creating a low density, low-pressure area in the center.

Computational results are presented at time t = 0.05 and using 100 cells in Figures 5. We note that standard MinMod limiter based scheme is not stable in this case and breaks down immediately. So, we are not able to present the results for the scheme. Significantly though, both CSPV and GSPV schemes are stable, and both schemes capture low-density and low-pressure areas. Also, the performance of both CSPV and GSPV schemes is comparable in this case. In addition, the computational time of CSPV scheme was 2.4978 seconds whereas GSPV took 2.3367 seconds.

Two Dimensional Near Vacuum Test Case

In this Section, we present a two-dimensional test case which contains low density and pressure areas. The test is the generalization of the one-dimensional case presented in Section 7.1.5. We consider the domain [-0.5, 0.5]×[-0.5×0.5] with outflow boundary conditions. This domain is filled with fluid at constant unit density with pressure p 11 = 2, p 12 = 0 and p 22 = 2. The velocity is taken to be 5 n, where n is the unit normal at the point directed outwards. So, the fluid is pushed outside the domain.

The numerical solutions are presented in Figures 6 at time t = 0.05 computed using 100 × 100 cells. In Figures 6(a), we have plotted the density for CSPV. We note that, as in the one-dimensional case, a low-density area has appeared in the center of the domain and both the limiters capture it. We also note that standard MinMod limiter based scheme fails in this case. Similarly, in Figure 6(b), we have plotted det(p) for CSPV limiters.

To compare both limiters more accurately, In Figures 6(c)-6(d) we have compared the cuts of two dimensional plots at y = 0 for density, pressure components p 11 , p 22 , and det(p). We observe that both schemes produce similar results. The simulation time of CSPV scheme was 617.37 seconds which was slightly more than the GSPV scheme which took 582.26 seconds.

Numerical Results: Non-Homogeneous Case

Two Rarefaction Waves with Gaussian Source Terms

To test the effect of the source terms we consider the one-dimensional test case from [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF]. We consider the domain [0, 4] with initial discontinuity at x = 2. The initial states are given in Table 5. The solution without source terms consists of two rarefaction waves leaving behind a low-density area in the middle, similar to the case of Section 7.1.5. For the one-dimensional test case with source terms, we considered Gaussian profile given by, W (x, t) = 25 exp(-200(x -2) 2 ).

The numerical results are presented using 500 cells at final time t = 0.1 in Figure [START_REF] Dubroca | Magnetic field generation in plasmas due to anisotropic laser heating[END_REF]. At this time without source term (Homogeneous Case) the low density area is not completely developed in the middle. So, proposed limiter and standard MinMod limiter produce similar results. However, when source terms are added, a near vacuum area has developed around point x = 2. So, one need a robust limiting scheme. we note that source term has created a low density area at the centre. Furthermore, we observe that the effects are anisotropic. In Figure (8(b)), we have plotted cut along line y = -x + 4 for both limiters (CSPV and GSPV ) and both schemes (O2-exp and O2-imex). We again note that solution for both IMEX and explicit schemes are similar and CSPV limiter provide more details compared to GSPV limiter. On the other hand, the computational time of GSPV limiter were slightly smaller than the CSPV limiter (See Table 8).

Realistic Simulation in Two Dimensions

In this example we consider a problem similar to the Example 7.4 from [START_REF] Berthon | An entropy preserving relaxation scheme for the Ten-Moments equations with source terms[END_REF]. We consider a domain of [0, 100] × [0, 100] filled with plasma with density 0.109885, initially at rest with pressure p 11 = p 22 = 1 and p 12 = 0. This is excited with source term only in x-direction with W (x, y) ≡ exp -x -50 10

Lemma 2 . 2

 22 The solution of the source ODE du dt = s(u) are in Ω if the initial conditions are in Ω.

Figure 1 :v 1 v 2 p 11 p 12 p 22

 122 Figure1: Convergence Rate: L1-error of density for CSPV, GSPV and MinMod solutions using HLLC solver. All the schemes achieve second order of accuracy.

  Determinant of pressure tensor

Figure 2 :

 2 Figure 2: Sod Shock Tube Problem: Numerical Solutions of MinMod, CSPV and GSPV limiters with HLLC flux using 100 cells at time t = 0.125. HLLC solver is used as numerical flux.

  Determinant of pressure tensor

Figure 3 :

 3 Figure 3: Two Shock Waves Problem: Numerical Solutions of MinMod, CSPV and GSPV limiters with HLLC flux using 100 cells. HLLC solver is used as numerical flux.

Figure 4 :

 4 Figure 4: Two Rarefaction Waves Riemann Problem: Numerical Solutions of Min-Mod, CSPV and GSPV limiters with HLLC flux using 100 cells. HLLC is used as numerical flux.

Figure 5 :

 5 Figure 5: Two Rarefaction Waves with Near Vacuum State: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 100 cells. HLLC solver is used as numerical flux.

Figure 6 :

 6 Figure 6: Two Dimensional Near Vacuum State: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 100 × 100 cells. HLLC solver is used as numerical flux.

Figure 7 :

 7 Figure 7: Two Rarefaction Waves with Gaussian Source Terms: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 500 cells. HLLC solver is used as numerical flux.

Figure 8 :

 8 Figure 8: Two Dimensional Near Vacuum State: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 100 × 100 cells. HLLC solver is used as numerical flux.

2 -

 2

Figure 9 :

 9 Figure 9: Two Dimensional Realistic Simulation with Absorption Coefficient v T = 0: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 100 × 100, 200 × 200 and 400 × 400 cells along line y = 50

Figure 10 :

 10 Figure 10: Two Dimensional Realistic Simulation with Absorption Coefficient v T = 1: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 100×100, 200 × 200 and 400 × 400 cells along line y = 50

Figure 11 :

 11 Figure 11: Two Dimensional Realistic Simulation Comparison for Absorption Coefficient v T = 0 and v T = 1: Numerical Solutions of CSPV and GSPV limiters with HLLC flux using 400 × 400 cells.

Figure 12 :

 12 Figure 12: Two Dimensional Realistic Simulation with Absorption Coefficient v T = 0 and v T = 1: Numerical Solutions of CSPV limiters with HLLC flux using 400 × 400 cells.

  12,i | < p n 11,i p n 22,i + p n 11,i ∆p 22,i + p n 22,i ∆p 11,i + ∆p 11,i ∆p 22,i ,

|p n 12,i -∆p 12,i | < p n 11,i p n 22,i -p n 11,i ∆p 22,i -p n 22,i ∆p 11,i + ∆p 11,i ∆p 22,i .

  The implicit Euler source update S i,1 ∆t is unconditionally Ω invariant.

	11,i is positive unconditionally. Similarly, we can check that p n+1 12,i = p n 12,i .
	This ensures the positivity of p n+1 11,i p n+1 22,i -p n+1 12,i	2	, unconditionally. We have the
	following result:
	Lemma 5.2

Table 1 :

 1 Initial Conditions for Sod's Shock Tube Riemann Problem problem are two constant states, given in the Table 1. Solutions are computed till time t = 0.125 and outflow boundary conditions are used.

Table 2 :

 2 Initial conditions for Two Shock Waves Riemann Problem for Homogenous Case ρ v 1 v 2 p 11 p 12 p 22

Table 3 :

 3 Initial conditions for Two Rarefaction Waves Riemann Problem ρ v 1 v 2 p 11 p 12 p 22

	Left 1 -5 0	2	0	2
	Right 1 5	0	2	0	2

Table 4 :

 4 Initial conditions for Two Rarefaction Waves Riemann Problem: Near Vacuum Test Case

Table 8 :

 8 Simulation time for Uniform Plasma State with Gaussian Source Terms

International Journal on Finite Volumes

We consider outflow boundary conditions. In addition, we also consider additional source term 2v T ρW for the energy part, two simulate inverse bremsstrahlung in International Journal on Finite Volumes

Acknowledgement: H. Kumar has been funded in part by SERB, DST grant with file no. YSS/2015/001663.

Both of the above schemes are Ω-invariant if each internal explicit time step satisfies corresponding time restriction. However, it is not possible to ensure this when calculating the time step using w n i . So, we reduce ∆t slightly using ∆t = ∆t -∆t 5 and check at each internal step if the corresponding stability condition is satisfied. If not, then we reduce initial time step more and repeat the process. In practice, we note that the solution will not change drastically over a time step. So, ∆t chosen above is sufficient.

Numerical Results

In this section, we will present the numerical experiments to exhibit the accuracy and robustness of the proposed algorithms. For the results presented here, we use HLLC flux as the numerical flux. This is because HLLC is more accurate (less dissipative) solvers compared to Lax-Friedrichs, Rusanov and HLLE solvers. Also, the focus of the work here is second-order reconstruction process, not first order solver. In the following Section 7.1, we present the test cases for the Ten-Moment equations without source terms. This is to show the accuracy and robustness of the proposed reconstruction procedure. In Section 7.2, we present computational results for the complete model, to show the robustness of source discretization. ]

Numerical Results for Homogeneous Case

Smooth Solutions: Rate of Convergence

To check the formal order of accuracy of the proposed limiters and their comparison with the standard MinMod limiters, we have designed a smooth solution of Ten Moment equations (1) without source terms, in one dimension case. We consider the domain [-0.5, 0.5] with initial density profile of ρ(x, 0) = 2 + sin(2πx). This is assumed to be moving with velocity v = (1, 0) . The pressure components are taken to be p 11 = p 22 = 1 and p 12 = 0. Assuming periodic boundary conditions, the exact solution is advection of density profile in x-direction, i.e. ρ(x, t) = 2 + sin(2π(x -t)). All other variables remains the same. We have plotted the L 1 -errors of density in Figure 1 for conservative (CSPV), general (GSPV) and standard MinMod slope limiters using HLLC solver. The errors are calculated using 20, 40, 80, 160, 320, 640 and 1280 cells. We have also plotted reference slope for second order convergence for comparison.

All the schemes converge with the second-order of accuracy. Furthermore, errors of the CSPV and GSPV solutions are same as of standard MinMod limiter. This is because the solution does not contain any low density or low-pressure area. So, CSPV and GSPV reduce to standard MinMod based reconstruction. Hence, all the schemes have almost same errors in this case.

Sod Shock Tube Problem

We consider the interval [-0.5, 0.5] to be domain and assume that the initial discontinuity is at x = 0.0. The initial conditions for the Sod's shock tube Riemann In Figure (7(a)), we have compared standard MinMod limiter for the homogeneous case with CSPV and GSPV limiter for the non-homogeneous case using explicit source discretization (O2-exp). Here, the solution using CSPV contains more details compared with GSPV. Similar observation is made for the case implicit source (O2-imex) in Figure (7(b)). In Figure (7(c)) we have compared GSPV limiter for the homogeneous case with GSPV in O2-exp and O2-imex. We find that there is no visible difference in explicit and IMEX scheme. This is because the source is not stiff and time step is governed by flux discretization. Similar observation can be made for CSPV limiter using Figure (7(d)). The simulation times of GSPV limiter were smaller than the CSPV limiter for both explicit and IMEX schemes (See Table 6).

Uniform Plasma State with Gaussian Source in Two Dimension

To demonstrate two dimensional effects of the source term we consider a uniform plasma state with initial conditions given in Table 7. in the domain [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF][START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF] × [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF][START_REF] Berthon | Numerical approximations of the 10-moment Gaussian-closure[END_REF] with outflow boundary conditions. The Gaussian source term considered is

Numerical results are presented using 100 × 100 cells at final time t = 0.1 in Figure [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]. In Figure (8(a)), we have plotted density using CPSV limiter with O2-exp ρ v 1 v 2 p 11 p 12 p 22 0.1 0 0 9 7 9

Table 7: Initial conditions for Uniform Plasma with Gaussian Source in Two Dimension plasma. We consider two cases, one with absorption coefficient v T to be 1 and another with no absorption i.e. v T = 0. In previous examples, we notice no visible difference in solutions for explicit and IMEX schemes. So, here we present results for IMEX schemes only. Solutions are evolved till time 0.5.

The numerical results are presented in Figures [START_REF] Hakim | Extended MHD Modelling with Ten-Moment Equations[END_REF], ( 10), ( 11) and [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. In Figure [START_REF] Hakim | Extended MHD Modelling with Ten-Moment Equations[END_REF], we ignore the absorption by setting coefficient v T = 0. We plot the density along the line y = 50 using 100 × 100, 200 × 200 and 400 × 400 cells for CSPV ( See Figure 9(a)) and GSPV ( See Figure 9(b)) limiters. In both the cases we observe convergence of the results when we refine the mesh.

To observe effect of the absorption coefficient we now consider case with v T = 1. In Figure 10 we plot density and p 11 for 100 × 100, 200 × 200 and 400 × 400 for both CSPV and GSPV along line y = 50. We observe that both density and pressure p 11 have converged. Furthermore, we have compared the results with v T = 0 and v T = 1 for 400 × 400 mesh in Figure [START_REF] Johnson | Gaussian-Moment Relaxation Closures for Verifiable Numerical Simulation of Fast Magnetic Reconnection in Plasma[END_REF]. We observe that density has decreased significantly at the centre of laser and pressure p 11 has changed its shape and now highest at that point.

In Figure 12 we have compared density contours with 400 × 400 mesh using CSPV limiter. Again we observed that density has decreased in the center when we take v T = 1.0.

Conclusion

In this work, we have presented positivity preserving second-order MUSCL scheme for Ten-Moment Gaussian closure model with source terms. This is achieved by enforcing suitable restrictions on the slopes of the reconstructed variables. We prove that under the presented restrictions on the slopes, schemes are positivity preserving. Furthermore, we have presented robust treatment of the source terms. We have presented numerical experiments for several test cases and compared the presented schemes with the standard second-order scheme. We note that the proposed restrictions on the slopes result in comparable results to the standard scheme for the cases without low density and pressure areas. For the cases, where we have low density or pressure areas, the presented schemes are shown to have superior robustness.