J Vides

B Braconnier
email: benjamin.braconnier@ifpenergiesnouvelles.fr

E Audit

C Berthon

B Nkonga

A Godunov-Type Solver for the Numerical Approximation of Gravitational Flows

Keywords: Euler-Poisson equations, approximate Riemann solver, relaxation scheme, source terms, gravitational effects

We present a new numerical method to approximate the solutions of an Euler-Poisson model, which is inherent to astrophysical flows where gravity plays an important role. We propose a discretization of gravity which ensures adequate coupling of the Poisson and Euler equations, paying particular attention to the gravity source term involved in the latter equations. In order to approximate this source term, its discretization is introduced into the approximate Riemann solver used for the Euler equations. A relaxation scheme is involved and its robustness is established. The method has been implemented in the software HERACLES [29] and several numerical experiments involving gravitational flows for astrophysics highlight the scheme.

Introduction

The present paper is devoted to the numerical approximation of the Euler equations when gravitational effects are taken into account. The associated solutions are governed by an Euler-Poisson model, given by the following system of partial differential equations (PDEs):

       ∂ t ρ + ∇•(ρu) = 0, ∂ t (ρu) + ∇•(ρu⊗u+ p) = -ρ∇φ, ∂ t (ρE) + ∇•((ρE+ p)u) = -ρu•∇φ, ∆φ = 4πGρ, (1.1)
where ρ > 0 is the density, u ∈ R d the velocity, and E the specific energy. The integer d refers to the space dimension. The thermodynamic pressure p is assumed to be governed by an equation of state p := p(ρ,ǫ),

where ǫ = E-|u| 2 /2 represents the specific internal energy. As usual (see [START_REF] Godlewski | Hyperbolic systems of conservations laws[END_REF][START_REF] Weyl | Shock waves in arbitrary fluids[END_REF]), the pressure law is assumed to satisfy

∂ ρ p(ρ,ǫ)+ p(ρ,ǫ) ρ 2 ∂ ǫ p(ρ,ǫ) > 0.
Here, G ≈ 6.67×10 -11 m 3 kg -1 s -2 is the gravitational constant. The gravitational potential φ is always a smooth function since it is the solution of the Laplace equation. For the sake of notation simplicity, the system is rewritten in compact form as

∂ t W+∇•F(W)+B(W)∇φ = 0, ∆φ = 4πGρ, (1.2)
with W = (ρ,ρu,ρE) T , F(W) = (ρu,ρu⊗u+ p,(ρE+ p)u) T ,

B(W) = ρ    0 T d I d u T   , (1.3)
where 0 d = (0,...,0) T is the null vector in R d and I d is an identity matrix of dimension d.

As usual, W:R d ×R + →Ω d is the state vector and F:Ω d →R 2+d is the flux function. Here, B : Ω d →R 2+d ×R d represents the gravitational contribution when multiplied by ∇φ. The convex set Ω d of the admissible state vectors is defined by

Ω d = W ∈ R 2+d ; ρ > 0, u ∈ R d , ǫ = E- |u| 2 2 > 0 . (1.4) System (1.
2) is completed with appropriate initial and boundary conditions that depend on the problem being considered, as will be seen in Section 5. Nevertheless, we note that in the field of astrophysics, it is common to find isolated boundary conditions for the potential φ, i.e., the potential is then approximated by a multipolar development of the mass distribution at the boundary (see [START_REF] James | The solution of Poisson's equation for isolated source distributions[END_REF]). For the initial data, we define W(x,0) = W 0 (x), φ(x,0) = φ 0 (x), (1.5) where the functions W 0 (x) and φ 0 (x) are chosen according to the specific physical settings of the problem addressed.

The Euler-Poisson model is certainly one of the most commonly used fluid models in plasma physics. However, simulations of plasma flows are not the purpose of the present work and will not be discussed in detail. We briefly wish to note that an abundant literature is devoted to the Euler-Poisson problem in plasmas (for instance, see [20-25, 39, 40]). The equations are used to describe ions and electrons flows, which are usually highly dynamic. A particular issue resides in the quasi-neutral limit (equality of the ion and electron densities) where the Poisson equation becomes singular. The derivation of efficient numerical schemes in this limit is a complex task. In fact, the space and time scales of the solution tend toward zero, leading to severe consistency and stability criteria. In order to overcome this difficulty, Degond et al. [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF][START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | Mach-number uniform asymptotic-preserving Gauge schemes for compressible flows[END_REF][START_REF] Degond | Analysis of an asymptotic preserving scheme for the Euler-Poisson system in the quasineutral limit[END_REF] (see also [START_REF] Degond | An Asymptotic Preserving scheme for the Euler equations in a strong magnetic field[END_REF]) proposed a reformulation of the Poisson equation which is not singular in the quasi-neutral limit. The resulting numerical scheme resolves successfully all the plasma regimes, which is particularly attractive when different regimes are present in a single domain.

However, in the field of astrophysical flows and fluid dynamics, the Poisson equation never bears singularities. The main challenge consists in discretizing the gravitational effects governed by B(W)∇φ. In fact, in the limit of a steady flow, the first three equations of system (1.1) give

   ∇•(ρu) = 0, ∇•(ρu⊗u+ p) = -ρ∇φ, ∇•((ρE+ p)u) = -ρu•∇φ. (1.6)
The main difficulty then resides in the design of a numerical scheme that preserves this asymptotic regime accurately. To address this issue, particular attention must be paid to the gravitational effects B(W)∇φ.

We now focus on system (1.2). We use an operator splitting approach (Yanenko splitting, first-order, or Strang splitting, second-order accurate), and in the context of the firstorder decomposition, equation (1.2) is solved in two steps. First, we restrict our attention to the system

∂ t W+∇•F(W)+B(W)∇φ = 0, (1.7)
where W ∈ Ω d is the unknown vector and the gravitational potential φ is an a priori given function. Equation (1.7) is in nonconservative form and the method used to solve it will be discussed later. Once the solution has been found, the second step consists in using the first component of W, i.e., ρ, to solve the Poisson equation

∆φ = 4πGρ, (1.8)
by means of a classical second-order finite difference approach. The operator splitting approach allows to decompose (1.2) into two subproblems and treat equations (1.7) and (1.8) individually, which are hyperbolic and elliptic, respectively. Several numerical strategies can be adopted in view of solving (1.7). The first idea [START_REF] Leveque | Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm[END_REF] is to decompose the solution into a hydrostatic steady part governed by system (1.6), and an acoustic dynamic part governed by the Euler equations without source terms. Classical Godunov-type based numerical schemes can be used to approximate the acoustic dynamic part. Unfortunately, this strategy tends to become much too complex in 3D or when a complete implicit formulation is required because of system (1.6) resolution.

The second idea is to construct a Godunov-type solver that takes into account the gravitational effects. A full conservative reformulation of (1.7) proposed in another context [START_REF] Balbus | On the dynamical foundations of α disks[END_REF] (see also Chieze [START_REF] Chièze | Elements of Hydrodynamics Applied to the Interstellar Medium, Starbursts: Triggers, Nature, and Evolution[END_REF]) could be used. Indeed, the authors introduce the specific energy, denoted E φ = E+φ, to rewrite (1.7) in the following conservative form:

       ∂ t (ρ) + ∇•(ρu) = 0, ∂ t (ρu) + ∇• ρu⊗u+ p+ ∇φ⊗∇φ 8πG = 0, ∂ t ρE φ + |∇φ| 2 8πG + ∇• ρE φ u-∇φ∂ t φ 4πG = 0.
Despite the conservation form, the nonlinearities involved in the differential operators make this approach of little interest from a scientific computing point of view. Another approach comes from the work of Greenberg et al. [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Greenberg | Analysis and approximation of conservation laws with source terms[END_REF] where the potential form of the source term is used with some benefits to derive well-balanced schemes. Such a numerical procedure has been widely improved in the literature (for instance, see [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF][START_REF] Gallouet | Some approximate Godunov schemes to compute shallowwater equations with topography[END_REF][START_REF] Gallardo | On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas[END_REF][START_REF] Jin | A steady-state capturing method for hyperbolic systems with geometrical source terms[END_REF][START_REF] Noelle | Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows[END_REF] for the shallow-water equations, [START_REF] Berthon | An HLLC scheme to solve the M1 Model of radiative transfer in two space dimensions[END_REF] in the framework of radiative transfer or [START_REF] Aregba-Driollet | Numerical approximation of Kerr-Debye equations[END_REF] for extensions of the Kerr-Debye model), and constitutes a relevant alternative to discretize (1.7). However, the nonlinear Riemann problem involved by this approach is too sophisticated in the present framework to be directly considered. Hence, the extension of this technique to the present context is obtained by involving a relaxation scheme. This relaxation technique is adopted to introduce relevant linearizations and thus derive a Suliciu-type relaxation solver for the Euler-Poisson system. From now on, we note that the obtained scheme is not well-balanced in the sense of [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Greenberg | Analysis and approximation of conservation laws with source terms[END_REF], but the steady states will be approximated with a better accuracy than with standard fractional step splitting approaches.

This paper is organized as follows. In the next section, we exhibit a relevant relaxation model to approximate the solutions of the system under consideration. One of the benefits of the proposed relaxation model is that it ensures a strong coupling between the Euler equations and the gravitational potential. The details of the Riemann solution associated with the homogeneous relaxation model are presented in Section 3. Section 4 is then devoted to the derivation of the relaxation scheme. Several numerical tests are presented and discussed in Section 5. Finally, concluding remarks are given in Section 6.

1-D relaxation model

We now consider the numerical approximation of system (1.7). For simplicity, the developments are proposed considering the 1-D model. The system can thus be written as

   ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) + ρ∂ x φ = 0, ∂ t (ρE) + ∂ x ((ρE+ p)u) + ρu∂ x φ = 0, (2.1)
where φ is an a priori given function, solution to the Poisson equation (1.8). With the definition of W, F and B given by (1.3), we write the system under the form

∂ t W+∂ x F(W)+B(W)∂ x φ = 0, (2.2)
and the 1-D state vector belongs to the phase space Ω 1 defined by (1.4).

The main difficulty lies with the source term, which requires an adequate approximation. We propose to derive a relaxation scheme, a classical numerical approach widely spread across the literature considering fluid flow simulations [1, 7-9, 17, 38]. After the work of Jin-Xin [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] or Bouchut [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF] (see also [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]), we propose to approximate the weak solutions of (2.1) by the weak solutions of a suitable first order PDE system with singular perturbations: the relaxation system. In order to make this strategy attractive, the relaxation system is generally designed to render the Riemann problem easily solvable. However, according to the work of Coquel-Perthame [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF] (for instance, see also [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Berthon | A positive preserving hight order VFRoe scheme for shallow water equations: A class of relaxation schemes[END_REF]), most of the nonlinearities of the initial system, called relaxation equilibrium system, must be preserved by the relaxation system to enforce accuracy of the resulting numerical scheme.

After the work of Suliciu [START_REF] Suliciu | On modelling phase transitions by means of rate-type constitutive equations, shock wave structure[END_REF][START_REF] Suliciu | Some stability-instability problems in phase transitions modelled by piecewise linear elastic or viscoelastic constitutive equations[END_REF] (see also [1, 7-9, 15, 17, 38] for several extensions), we first propose to relax the thermodynamic pressure p. In this sense, we suggest substituting p by an approximation π governed by a suitable evolution law supplemented by a relaxation source term. Such an evolution law can be deduced from the pressure equation. After usual computations (for instance, see [START_REF] Godlewski | Hyperbolic systems of conservations laws[END_REF]), the pressure equation is

∂ t p+u∂ x p+ρc 2 ∂ x u = 0, (2.3)
where the sound speed c is defined as

c 2 = ∂ ρ p(ρ,ǫ)+ p(ρ,ǫ) ρ 2 ∂ ǫ p(ρ,ǫ) > 0.
Let us underline that relation (2.3) does not depend on the gravitational potential φ. From (2.3), we suggest the following evolution equation to govern the additional variable π:

∂ t π +u∂ x π + a 2 ρ ∂ x u = 1 δ (p-π). (2.4)
A sub-characteristic Whitham condition [START_REF] Whitham | Linear and Nonlinear Waves[END_REF] will be considered later to fix the relaxation parameter a. The choice of a will be done to satisfy some robustness properties of the scheme. The parameter δ must go to zero in order to enforce a relaxation equilibrium limit characterized by π = p.

In addition, we propose to introduce a relaxation procedure to approximate the potential φ. As a consequence, φ is replaced by the new variable ψ satisfying

∂ t ψ = 1 δ (φ-ψ). (2.5)
The suggested relaxation model is then given by

             ∂ t ρ+∂ x (ρu) = 0, ∂ t (ρu)+∂ x (ρu 2 +π)+ρ∂ x ψ = 0, ∂ t (ρE)+∂ x ((ρE+π)u)+ρu∂ x ψ = 0, ∂ t π +u∂ x π + a 2 ρ ∂ x u = 1 δ (p-π), ∂ t ψ = 1 δ (φ-ψ).
(2.6)

Let us emphasize that, in the limit of δ = 0, at least formally, we have π = p and ψ = φ and the evolution laws of the density, momentum and energy (ρ, ρu, and ρE, respectively) in (2.6) are thus equivalent to those of (2.1).

For the sake of simplicity in the forthcoming developments, we introduce the condensed form of (2.6), i.e., the compact form of the relaxation system corresponding to (2.2), as

∂ t W δ +∂ x F δ (W δ )+B δ (W δ )∂ x ψ = 1 δ R δ (W δ ), (2.7)
which is written in terms of the relaxed variables given by

W δ = (ρ,ρu,ρE,ρπ,ψ) T , F δ (W δ ) = ρu,ρu 2 +π,(ρE+π)u, ρπ +a 2 u,0 T , B δ (W δ ) = (0,ρ,ρu,0,0) T , R δ (W δ ) = (0,0,0,ρ(p-π),φ-ψ) T .
The set of admissible relaxation state vectors is defined by

Ω δ = W δ ∈ R 5 ; ρ > 0, u ∈ R, ǫ = E- u 2 2 > 0, π ∈ R, ψ ∈ R ,
where the density ρ and specific internal energy ǫ are positive and must remain positive at all times to be physically meaningful.

The Riemann solution W R δ

We present the solution of the Riemann problem for the relaxation model (2.6) with a vanishing relaxation source term, i.e., δ = ∞. Simultaneously, relevant choices of the parameter a will be specified to enforce a suitable positive preserving property. First, we consider the algebra of system (2.6). We propose to introduce a change of variables and use the state vector V δ defined as

V δ = (ρ,u,ǫ,π,ψ) T .
Then, omitting the source term, (2.6) can be reformulated as

∂ t V δ +A δ (V δ )∂ x V δ = 0, (3.1)
where

A δ (V δ ) =       u ρ 0 0 0 0 u 0 1/ρ 1 0 π/ρ u 0 0 0 a 2 /ρ 0 u 0 0 0 0 0 0       .
An easy computation gives µ 0 = 0, µ u 1 = µ u 2 = u and µ ± = u±a/ρ as eigenvalues of the above system. It is then clear that the relaxation system (3.1) involves a stationary wave. Following the pioneering work of Greenberg-LeRoux [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Greenberg | Analysis and approximation of conservation laws with source terms[END_REF], we introduce a source term in the associated Riemann solver which naturally comes with a stationary wave (see also [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF]). With clear notation, the eigenvectors of A δ (V δ ) are

r 0 = ρ 2 ,-uρ,π,a 2 ,(u 2 ρ 2 -a 2 )/ρ T , r u 1 = (1,0,0,0,0) T , r u 2 = (0,0,1,0,0) T , r ± = ρ 2 ,±a,π,a 2 ,0 T ,
where r u 1 and r u 2 are the eigenvectors associated to the double eigenvalue u, from now on denoted by µ u . Using simple calculations, we deduce that all the characteristics fields of (2.6) are linearly degenerate. We complete the algebra of system (2.6) by exhibiting the Riemann invariants. Let us recall that the Riemann invariants, denoted I, associated with the eigenvector r, are defined by ∇I •r = 0. After straightforward computations, with the setting of π = π +a 2 /ρ and ǫ = ǫ-π 2 /(2a 2 ), we find the Riemann Invariants across each wave.

• Across the µ 0 -wave we have

I 0 1 = ρu, I 0 2 = π, I 0 3 = ǫ, I 0 4 = ψ+ u 2 2 - a 2 2ρ 2 . ( 3.2)
• Across the µ u -wave we have

I u 1 = u, I u 2 = π, I u 3 = ψ. (3.3)
• Across the µ ± -waves we have

I ± 1 = u± a ρ , I ± 2 = π, I ± 3 = ǫ, I ± 4 = ψ. (3.4)
We introduce several notations for the sake of simplicity. We begin by recalling that the specific volume τ is the reciprocal of the density, namely τ = 1/ρ. Thus, several of the quantities previously defined can be written in terms of τ instead of ρ. Additionally, given two constants X L and X R , we define

[X] = X R -X L , (3.5) X = X R +X L 2 . ( 3.6)
Now, we fix two admissible states, W L δ and W R δ , in Ω δ to define

W 0 δ (x) = W L δ if x < 0, W R δ if x > 0,
the initial data of the Riemann problem for system (2.6) δ=∞ . Thus, we want to find the solution to

∂ t V δ +A δ (V δ )∂ x V δ = 0 with Riemann initial data V 0 δ = V L if x < 0, V R if x > 0.
Since all eigenvalues are linearly degenerate, we can compute the exact solution W R δ to the Riemann problem, which has four wave speeds and three intermediate states.

To evaluate the intermediate constant states, recall that the Riemann invariants associated with an eigenvector are continuous across the corresponding contact discontinuity. The speeds of the stationary, contact, slowest and fastest waves are defined by S 0 , S M , S L and S R , respectively. The characterization of the solution is conditioned by the knowledge of the wave ordering, which is not straightforward. By construction, the wave speeds are

S 0 = 0, S L = u L -aτ L , S R = u R +aτ R . (3.7)
Since τ L > 0 and τ R > 0, we find that the difference S R -S L = u R -u L +a(τ R +τ L ) is positive as long as a is large enough. From a numerical point of view, the choice of the parameter a is crucial since it governs (in a way to be defined) the numerical diffusion involved in the scheme. As a consequence, choosing an a too large will produce a very viscous scheme while choosing it too small will lead to an unstable scheme. As usual, the parameter a cannot be smaller than the exact impedance Z = ρc involved in the initial system (2.1) (see [START_REF] Berthon | Numerical approximations of the 10-moment Gaussian closure[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] Chalons | Relaxation approximation of the Euler equations[END_REF]).

Hence, since S R >S L for large enough values of a, the wave order problem comes from the position of S R , S L and S M compared to the stationary wave S 0 . In order to cover all possibilities, we set the wave speed

σ = ū- [π] 2a , (3.8)
which corresponds to the velocity in the intermediate region of the Suliciu relaxation model used to approximate the Riemann solution of the Euler equations (see Bouchut [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws, and Well-Balanced Schemes for Sources[END_REF], for instance). Thus, only four distinct cases are possible and they are studied below.

In all cases, one can easily recover the intermediate velocity σ by setting [ψ] = 0. Recall that we consider only small values of [ψ], making this a valid approach to segregate Cases 2 and 3.

Case 1: S L > 0

In this case, we suppose that the wave ordering is S 0 = 0 < S L < S M < S R . Thus, we search for the exact Riemann solution having the following structure:

V L S 0 V * L S L V * * L S M V * R S R V R V L = ( τ L , u L , πL , ǫL , ψ L ) T V * L = ( τ * L , u * L , π * L , ǫ * L , ψ * L ) T V * * L = ( τ * * L , u * * L , π * * L , ǫ * * L , ψ * * L ) T V * R = ( τ * R , u * R , π * R , ǫ * R , ψ * R ) T V R = ( τ R , u R , πR , ǫR , ψ R ) T
The unknowns that compose this solution are then determined by considering the system given by the continuity of the Riemann invariants (3.2)-(3.4) across each field. Since the velocity and the quantity π are constant across S M , we have

S M = u * * L = u * R , π * * L = π * R , (3.9)
respectively. Additionally, we can easily derive the equalities

ψ R = ψ * R = ψ * * L = ψ * L and πL = π * L = π * * L , πR = π * R , (3.10) ǫL = ǫ * L = ǫ * * L , ǫR = ǫ * R . (3.11)
Across the S 0 -wave, we have the relations

u L /τ L = u * L /τ * L , (3.12a
)

ψ L + (u L ) 2 2 - a 2 2 (τ L ) 2 = ψ * L + (u * L ) 2 2 - a 2 2 (τ * L ) 2 . (3.12b)
Similarly, across the wave with speed S L , we have

u * L -aτ * L = u * * L -aτ * * L , (3.13)
and, across the S R -wave, we get

u * R +aτ * R = u R +aτ R . (3.14)
From (3.12a), we obtain

u * L = u L τ * L τ L . (3.15)
We substitute relations (3.15) and

ψ * L = ψ R into equation (3.12b) to produce (τ * L ) 2 = (τ L ) 2 - 2[ψ] (u L /τ L ) 2 -a 2 .
Since the specific volume τ * L has to be positive, we derive the expression

τ * L = τ L 1- 2[ψ] (u L ) 2 -a 2 (τ L ) 2 .
(3.16)

Recall that π = π+a 2 τ. We express the right hand sides of (3.10) in terms of this quantity. Thus, we obtain πL

= π * * L +a 2 τ * * L and πR = π * R +a 2 τ * R .
By performing algebraic manipulations on the above two equations and considering the second equation in (3.9), we derive

τ * R = τ * * L + [ π] a 2 . (3.17) Substitution of u * L from (3.15) into (3.13) gives u * * L = (u L /τ L -a)τ * L +aτ * * L = S M , (3.18)
and use of the relation

u * R = u R +aτ R -aτ * R = S M derived from (3.14) leads to τ * * L = 1 a (u R +aτ R -aτ * R -(u L /τ L -a)τ * L ). (3.19)
Then, using (3.17) and (3.19), we deduce the values

τ * * L = 1 2a S R -S L 1- 2[ψ] (u L ) 2 -a 2 (τ L ) 2 - [ π] a , (3.20
)

τ * R = 1 2a S R -S L 1- 2[ψ] (u L ) 2 -a 2 (τ L ) 2 + [ π] a . ( 3.21)
We skip the computation for the other variables since they can be easily deduced from one another. However, from (3.16), we can derive the necessary condition

[ψ] < 1 2 (u L ) 2 -a 2 (τ L ) 2 , ( 3.22)
where

(u L ) 2 -a 2 (τ L ) 2 >0 since S L >0, that is, u L +aτ L >u L -aτ L >0.
Additionally, in order for τ * * L and τ * R to be positive, the following condition must be satisfied:

|[ π]| < a S R -S L 1- 2[ψ] (u L ) 2 -a 2 (τ L ) 2 .
If the previous condition is not satisfied, to avoid the appearance of vacuum, we decide to degenerate the S L -wave by choosing a larger value of a such that S L <0<S R , σ>0 and, therefore, switch to Case 2.

Case 2: S

L < 0 < S R , σ > 0
Assume that the exact solution follows the wave ordering given by S L < S 0 < S M < S R , which implies looking for the exact Riemann solution of the form:

V L S L V * L S 0 V * * L S M V * R S R V R V L = ( τ L , u L , πL , ǫL , ψ L ) T V * L = ( τ * L , u * L , π * L , ǫ * L , ψ * L ) T V * * L = ( τ * * L , u * * L , π * * L , ǫ * * L , ψ * * L ) T V * R = ( τ * R , u * R , π * R , ǫ * R , ψ * R ) T V R = ( τ R , u R , πR , ǫR , ψ R ) T
Once again, using the continuity of the Riemann invariants (3.2)-(3.4) across each field, we obtain

u L -aτ L = u * L -aτ * L , (3.23a) u * L /τ * L = u * * L /τ * * L , (3.23b
)

ψ * L + (u * L ) 2 2 - a 2 2 (τ * L ) 2 = ψ * * L + (u * * L ) 2 2 - a 2 2 (τ * * L ) 2 , (3.23c) u * R +aτ * R = u R +aτ R . (3.23d)
Since the velocity and π are constant across S M , we have

S M = u * * L = u * R , π * * L = π * R . (3.24)
Additionally, we can easily derive the equalities

ψ R = ψ * R = ψ * * L , ψ L = ψ * L , (3.25) πL = π * L = π * * L , πR = π * R , (3.26) ǫL = ǫ * L = ǫ * * L , ǫR = ǫ * R . (3.27)
By using equation (3.23b), we solve for the velocity u * L as

u * L = u * * L τ * L τ * * L . (3.28) Direct substitution of relations (3.28), ψ * L = ψ L and ψ * * L = ψ R into equation (3.23c) yields τ * L τ * * L 2 = 1+ 2[ψ] (u * * L ) 2 -a 2 (τ * * L ) 2 .
(3.29)

From (3.23a) and (3.28), we get

τ * L τ * * L = u L -aτ L u * * L -aτ * * L . (3.30)
By expanding the right hand sides of (3.26), recalling that π = π +a 2 τ, we obtain

πL = π * * L +a 2 τ * * L and πR = π * R +a 2 τ * R .
We perform algebraic manipulations on the above two equations and use the second equation in (3.24). Thus, we derive

τ * R = τ * * L + [ π] a 2 , (3.31)
which we can substitute into equation (3.23d) to obtain

u * * L = u * R = α-aτ * * L , (3.32)
with

α = S R - [ π] a .
Substitution of relations (3.32) and (3.30) into (3.29) gives the second order polynomial

(α-2aτ * * L ) 2 + 2[ψ] α (α-2aτ * * L )-(S L ) 2 = 0, (3.33)
which has two real roots. Thus, we find

τ * * L = 1 2a   α+ [ψ] α ± [ψ] α 2 +(S L ) 2   . (3.34)
The suitable root is chosen by considering physical criteria. Recall that the specific volume τ * * L has to be a positive quantity. Therefore, we select the solution capable of restoring the Riemann solution to the Suliciu model as [ψ] tends to zero. Hence,

τ * * L = 1 2a   α+ [ψ] α + [ψ] α 2 +(S L ) 2   , (3.35)
and

τ * R = 1 2a   S R + [ π] a + [ψ] α + [ψ] α 2 +(S L ) 2   . (3.36)
Direct substitution of τ * * L into (3.32) leads to

S M = 1 2   α- [ψ] α - [ψ] α 2 +(S L ) 2   . ( 3
(α-|S L |) = 1 2 u L +u R - [π] a = σ. (3.38)
By performing steps similar to the ones used to find τ * * L , we find the second order polynomial

(S L +2aτ * L ) 2 - 2[ψ] S L (S L +2aτ * L )-α 2 = 0, (3.39)
which is used to obtain an expression for τ * L . After some algebraic manipulations, we get

τ * L = 1 2a   -S L + [ψ] S L +sgn(α) [ψ] S L 2 +α 2   . (3.40)
Again, we omit the computation for the remaining unknowns since they can be easily deduced from one another.

Case 3:

S L < 0 < S R , σ < 0
In this case, we assume that the wave ordering is S L < S M < S 0 < S R , and we look for the exact Riemann solution having the structure:

V L S L V * L S M V * * R S 0 V * R S R V R V L = ( τ L , u L , πL , ǫL , ψ L ) T V * L = ( τ * L , u * L , π * L , ǫ * L , ψ * L ) T V * * R = ( τ * * R , u * * R , π * * R , ǫ * * R , ψ * * R ) T V * R = ( τ * R , u * R , π * R , ǫ * R , ψ * R ) T V R = ( τ R , u R , πR , ǫR , ψ R ) T
The continuity of the Riemann invariants (3.2)-(3.4) across the linearly degenerate fields reads

u L -aτ L = u * L -aτ * L , (3.41a) u * R /τ * R = u * * R /τ * * R , (3.41b
)

ψ * R + (u * R ) 2 2 - a 2 2 (τ * R ) 2 = ψ * * R + (u * * R ) 2 2 - a 2 2 (τ * * R ) 2 , (3.41c) u * R +aτ * R = u R +aτ R . (3.41d)
Across the wave with speed S M , we have

S M = u * * R = u * L , π * * R = π * L . (3.42)
In addition, we also derive the equalities

ψ L = ψ * L = ψ * * R , ψ R = ψ * R , (3.43) πR = π * R = π * * R , πL = π * L , (3.44) ǫR = ǫ * R = ǫ * * R , ǫL = ǫ * L . (3.45)
In order to find the solution for this case, we use a methodology similar to the one presented in Case 2 and derive the specific volumes

τ * R = 1 2a   S R + [ψ] S R -sgn(β) [ψ] S R 2 +β 2   , (3.46
)

τ * * R = 1 2a   -β+ [ψ] β + [ψ] β 2 +(S R ) 2   , (3.47
)

τ * L = 1 2a   -S L + [ π] a + [ψ] β + [ψ] β 2 +(S R ) 2   , (3.48)
with

β = S L - [ π] a .
We omit the computation for the remaining unknowns, which easily follows from the relations defined in (3.41)- (3.45). However, due to its importance, we give the solution for S M = u * * R = u * L , with u * L = S L +aτ * L . Thus,

S M = 1 2   β + [ψ] β + [ψ] β 2 +(S R ) 2   . (3.49) If β is large enough such that |[ψ]| << |β|, then the value of S M is approximately 1 2 (β+|S R |) = 1 2 u L +u R - [π] a , (3.50)
and consequently, one recovers the speed σ, as defined in equation (3.8).

Case 4: S R < 0

Suppose that the exact solution follows the wave ordering S L < S M < S R < S 0 = 0 and is composed of five constant states:

V L , V * L , V * * R , V * R , and V R . V L S L V * L S M V * * R S R V * R S 0 V R V L = ( τ L , u L , πL , ǫL , ψ L ) T V * L = ( τ * L , u * L , π * L , ǫ * L , ψ * L ) T V * * R = ( τ * * R , u * * R , π * * R , ǫ * * R , ψ * * R ) T V * R = ( τ * R , u * R , π * R , ǫ * R , ψ * R ) T V R = ( τ R , u R , πR , ǫR , ψ R ) T
By means of the trivial Riemann invariants, we search for the Riemann solution having the structure presented above. Across S M , we have

S M = u * * R = u * L , π * * R = π * L . (3.51)
We also derive the equalities

ψ L = ψ * L = ψ * * R = ψ * R and πR = π * R = π * * R , πL = π * L , (3.52) ǫR = ǫ * R = ǫ * * R , ǫL = ǫ * L . (3.53)
Additionally, we have the relations

u L -aτ L = u * L -aτ * L , (3.54a) u R /τ R = u * R /τ * R , (3.54b) ψ R + (u R ) 2 2 - a 2 2 (τ R ) 2 = ψ * R + (u * R ) 2 2 - a 2 2 (τ * R ) 2 , (3.54c) u * R +aτ * R = u * * R +aτ * * R . (3.54d)
Using an approach similar to the one presented in Case 1, we derive the values for the following specific volumes:

τ * R = τ R 1+ 2[ψ] (u R ) 2 -a 2 (τ R ) 2 ,
(3.55)

τ * * R = 1 2a -S L +S R 1+ 2[ψ] (u R ) 2 -a 2 (τ R ) 2 + [ π] a , (3.56) τ * L = 1 2a -S L +S R 1+ 2[ψ] (u R ) 2 -a 2 (τ R ) 2 - [ π] a . ( 3

.57)

Again, we decide to skip the computation for the other variables since they can be easily deduced from one another. However, from (3.55), we derive the necessary condition

[ψ] > - 1 2 (u R ) 2 -a 2 (τ R ) 2 . (3.58)
In this case, recall that S R < 0. Therefore, u R -aτ R < u R +aτ R < 0 and the quantity (u R ) 2 -a 2 (τ R ) 2 is strictly positive. Now, for τ * * R and τ * L to be positive specific volumes, we must satisfy the following condition:

|[ π]| < a -S L +S R 1+ 2[ψ] (u R ) 2 -a 2 (τ R ) 2 .
If the previous condition is not satisfied, in order to avoid the appearance of vacuum, we decide to degenerate the S R -wave by choosing a larger value of a such that we can switch to Case 3.

Remarks

Summarizing the necessary conditions from Cases 2 and 3, we obtain

|[ψ]| < min (|α|,|β|).
In addition, from Cases 1 and 4, we have

|[ψ]| < 1 2 min (u R ) 2 -a 2 (τ R ) 2 , (u L ) 2 -a 2 (τ L ) 2 .
Thus, we confirm that the solver is robust for small values of [ψ].

1-D relaxation scheme

Based on the relaxation model (2.6), we present the numerical procedure used to approximate the solution of (2.1). We consider uniform space meshes of size ∆x = x i+1/2 -x i-1/2 for i ∈ Z. The time increment is given by ∆t and we note t n+1 = t n +∆t for n in N.

As mentioned in Section 2, we suppose for the moment that for each time step t n , the approximate potential φ h is known. Hence,

φ h (x,t n ) = φ n i if x ∈ x i-1 2 ,x i+ 1 2 . ( 4.1)
At time t = 0, the initial data is given by

W 0 i = 1 ∆x x i+ 1 2 x i-1 2 W 0 (x)dx, (4.2)
where W 0 (x) is the initial data function defined in (1.5). Additionally, at time t n , we assume that a piecewise constant approximate solution, denoted by W h (x,t n ), is known and defined by

W h (x,t n ) = W n i if x ∈ x i-1 2 ,x i+ 1 2 . (4.3)

Two-step splitting technique

The obtained approximations are evolved in time using a two-step splitting technique, which will be described in detail subsequently.

First step: evolution in time

During the first step, the time evolution is obtained considering the relaxation model (2.6) but omitting the source term, i.e., δ = ∞. More specifically, for all 0 < t < ∆t, we look for the weak solutions W h δ (x,t n +t) of the Cauchy problem

∂ t W δ +∂ x F δ (W δ )+B δ (W δ )∂ x ψ = 0, (4.4) subject to initial data W h δ (x,t n ) given by W h δ (x,t n ) = (W δ ) n i if x ∈ x i-1 2 ,x i+ 1 2 (4.5) = (ρ n i , (ρu) n i , (ρE) n i , (ρπ) n i = ρ n i p n i , ψ n i = φ n i ). (4.6)
We note that this initial data coincides with a relaxation equilibrium state since we have enforced π n i = p n i and ψ n i = φ n i . The solution W h δ (x,t n +t) of the above Cauchy problem (4.4)-(4.5) is composed of the juxtaposition of the non-interacting Riemann problem solutions set at each interface x i+1/2 for i ∈ Z and for all 0 < t < ∆t, with ∆t small enough such that it satisfies the following CFL like condition:

∆t ∆x max i∈Z |µ - i+ 1 2 |,|µ + i-1 2 | ≤ 1 2 , ( 4.7)
where the pair (µ - i+1/2 ,µ + i+1/2 ) describes the slowest and fastest wave speeds, respectively. The description and details of these waves is given in Section 3.

For the sake of clarity in the notations, we denote by W R δ ( x t ;W L δ ,W R δ ) the Riemann solution associated with (4.4) where the initial data at time t n is given by two constant states W L δ and W R δ separated by a discontinuity located at x i+ 1 2

. Consequently, we have

W h δ (x,t n +t) = W R δ x-x i+ 1 2 t ;(W δ ) n i ,(W δ ) n i+1 , x ∈ (x i ,x i+1 ), 0 < t < ∆t,
where the explicit form of W R δ is described in Section 3. The projection of this solution over the piecewise constant functions reads

(W δ ) n+1,- i = 1 ∆x x i+ 1 2 x i-1 2 W h δ (x,t n +∆t)dx = 1 ∆x x i x i-1 2 W R δ x-x i-1 2 ∆t ;(W δ ) n i-1 ,(W δ ) n i dx + 1 ∆x x i+ 1 2 x i W R δ x-x i+ 1 2 ∆t ;(W δ ) n i ,(W δ ) n i+1 dx. (4.8)
Let us note that the standard conservative flux balance cannot be reached in the present work because of the potential source term. However, an easy extension of the well-known formalism introduced by Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyper-bolic conservation laws[END_REF] gives

(W δ ) n+1,- i = (W δ ) n i - ∆t ∆x F L δ ((W δ ) n i ,(W δ ) n i+1 )-F R δ ((W δ ) n i-1 ,(W δ ) n i ) , (4.9)
where

F L δ ((W δ ) n i ,(W δ ) n i+1 ) =F δ ((W δ ) n i ) - 1 ∆t x i+ 1 2 x i W R δ x-x i+ 1 2 ∆t ;(W δ ) n i ,(W δ ) n i+1 -(W δ ) n i dx, (4.10) F R δ ((W δ ) n i ,(W δ ) n i+1 ) =F δ ((W δ ) n i+1 ) + 1 ∆t x i+1 x i+ 1 2 W R δ x-x i+ 1 2 ∆t ;(W δ ) n i ,(W δ ) n i+1 -(W δ ) n i+1 dx.
(4.11)

Second step: relaxation

During the second step of the scheme, the following system is solved:

∂ t W δ = 1 δ R δ (W δ ), (4.12)
where the initial data is defined by the piecewise constant approximation (W δ ) n+1,- i . As δ tends to zero, the updated approximate equilibrium solution W n+1 i , ∀i ∈ Z, is given by

W n+1 i = ρ n+1,- i ,(ρu) n+1,- i ,(ρE) n+1,- i T , ( 4

Remarks and summary

The description of the numerical scheme to approximate the solutions to subsystem (2.1) is complete once the relaxation parameter a is correctly characterized. Concerning this parameter, introduced in equation (2.4), let us note that it can be chosen locally at the interfaces since at each interface x i+1/2 , the Riemann solution depends not only on the left and right states (W δ ) n i and (W δ ) n i+1 , but also on the parameter a. Therefore, this constant a := a i+1/2 must be fixed according to the robustness conditions described in Section 3. Now, involving the usual framework of the finite volume methods, the proposed relaxation scheme is summarized below:

W n+1 i = W n i - ∆t ∆x F L,n i+ 1 2 -F R,n i-1 2 , ( 4.14)
where

F L,n i+ 1 2 = F L δ (W δ (W n i ),W δ (W n i+1 )), (4.15) F R,n i+ 1 2 = F R δ (W δ (W n i ),W δ (W n i+1 )), (4.16)
with W δ (W n i ) = (W δ ) n i defined according to the relaxation equilibrium, i.e., π n i = p n i and ψ n i = φ n i . Indeed, considering the projection step and the relaxation system solution, the numerical fluxes can be written as functions of the relaxation equilibrium states

F L,n i+ 1 2 = F L,n i+ 1 2 (ρ n i , u n i , (ρE) n i , φ n i , ρ n i+1 , u n i+1 , (ρE) n i+1 , φ n i+1 ), (4.17)
F R,n i+ 1 2 = F R,n i+ 1 2 (ρ n i , u n i , (ρE) n i , φ n i , ρ n i+1 , u n i+1 , (ρE) n i+1 , φ n i+1 ). (4.18)
Let us remark that we kept the notation F α,n i+1/2 , with α = L or R, for the fluctuations since the numerical scheme is not conservative. In fact, the nonconservative operator, coming from B(W)∂ x φ in (2.2), is null except at the interfaces x i+1/2 (i ∈Z) because of the evolution equation (2.5) that governs the relaxation potential. In Section 3, we clearly see that the continuity of the flux function is lost across any interface.

Finally, we proceed to remark that the above presented scheme assumes the approximate potential at time t n to be known, as mentioned at the beginning of this section. This comes from the fact that, as explained in Section 1, we use an operator splitting approach to decompose system (1.2) into two subsystems (1.7) and (1.8). Hence, in the one-dimensional case being considered, once we obtain the approximate solution W n+1 i by means of (4.14), we use ρ n+1 i to solve the Poisson equation and thus obtain φ n+1 i . We recall that the one-dimensional Poisson equation is

∂ xx φ = 4πGρ, (4.19)
and it can be discretized by means of a classical second-order finite difference approach which yields a tridiagonal matrix. There are numerous ways of solving the resulting matrix equation that can be categorized into direct and iterative methods.

Numerical results

The purpose of this section is to illustrate the interest of the proposed relaxation solver for the Euler-Poisson model. In order to address numerical issues, several astrophysical tests are performed. In the first two, we consider a fluid subject to a constant external gravitational field. In the third test, the Lane-Emden equation [START_REF] Kippenhahn | Stellar Structure and Evolution[END_REF], which describes the hydrostatic equilibrium of a self-gravitating star.

For the problems involving hydrostatic equilibrium flows, the proposed relaxation scheme is compared with a standard fractional step splitting method. This approach also consists on treating subsystems (1.7) and (1.8) individually, with the sole difference being that for the former subsystem, one alternates between solving the Euler equations without source terms

∂ t W+∇•F(W) = 0, (5.1)
and the ordinary differential equation

∂ t W = -B(W)∇φ, (5.2)
for each time step. Hence, we avoid the necessity of directly incorporating the source terms inside the solver, making this approach quite simple. Hereafter, we will refer to this splitting method as the standard method.

1-D Equilibrium flow

Let us consider a 1-D equilibrium flow. In this case, we neglect the gravitational interactions between particles and suppose that they are subject to a uniform gravitational field. Therefore, we define the gravitational potential φ = gx, with constant g > 0. We suppose that the fluid is governed by an isentropic equation of state p = c 2 ρ, where c is a non-zero constant. Under this assumption, the hydrostatic equilibrium flow satisfies ∂ x p = -ρg and we deduce the density as a function of position

ρ eq (x) = ρ(0)e -gx/c 2 . ( 5.3)
From this relation, we can define the characteristic scales associated to this experiment:

the gravitational length L re f = c 2 /g and time T re f = L re f /g. In order to illustrate the accuracy of the proposed method, we propose to compute this equilibrium flow and perform an error analysis in the L 2 norm. In this context, we set c = 1m.s -1 , L = 1m and g = 10m.s -2 . The characteristic scales are then L re f = 0.1m and T re f ≈ 0.316s. Hence, we consider meshes with L ≥ L re f and run the experiment a sufficiently long time T = 50s > 100T re f . The CFL number is chosen as 0.5. At the initial time, the density is set according to (5.3) with ρ(0) = 10. Additionally, we have homogeneous Dirichlet boundary conditions for the velocity and non-homogeneous Neumann boundary conditions for the other variables, which are set to ensure the hydrostatic equilibrium.

Figure 1 displays the L 2 error of the estimated density compared to the analytical solution (5.3) for different mesh sizes ∆x in a log-log scale. The slope of the error gives the order of accuracy of the method. Using the least squares method, we obtain a linear fit y = 0.9556 x+0.8632 and thus, the method's order is approximately 0.96.

On the left in Figure 2, we plot the logarithm of the L 2 error versus time for 1000 grid points. We can observe that the L 2 error grows exponentially for a short time and then stabilizes itself, and this holds true for both methods being compared. However, the error for the relaxation method is considerably smaller than the one for the standard method. Additionally, the evolution in time of the L 2 norm of the velocity, in logarithmic scale, is shown on the right in Figure 2 and one can observe that the velocity field is far lower for the relaxation method. The relaxation scheme clearly enables a better approximation of the equilibrium solution (5.3). ) is used, the total number of time steps for both methods is of the same order of magnitude, specifically, 1.001×10 5 and 1.006×10 5 steps for the relaxation and standard methods, respectively.

Rayleigh Taylor instability

In this section, we propose to validate the proposed method with a Rayleigh Taylor instability occurring as a consequence of a heavy fluid driven into a lighter one under the acceleration of a gravitational field. Initially, an unstable interface separates the fluids with different densities.

We carry out this experiment in the two-dimensional plane approximation (translational invariance along the z-axis). The derived 1-D numerical method given by (4.14) is easily extended to consider 2-D simulations. As is usual, we suggest

W n+1 ij = W n ij - ∆t ∆x F L,n i+ 1 2 -F R,n i-1 2 - ∆t ∆y G L,n j+ 1 2 -G R,n j-1 2 ,
where G L,n j+1/2 and G R,n j+1/2 are the numerical flux functions in the y-direction. We emphasize that the 2D Euler-Poisson extension is standard and does not involve difficulties in its derivation.

The two-dimensional computational domain is set to L x = 4m and L y = 1m and the gravity in the x-direction to g = -1.0m.s -2 . The domain contains a gas governed by an ideal gas equation of state p = (γ-1)ρǫ, with γ = 1.4. The gas densities to the right and left of the interface are ρ L = 1kg.m -3 and ρ R = 2kg.m -3 , respectively. The shape of the discontinuity or interface is given by x = 1 2 L x 1-1 10 cos(( y L y -1 2 )π) . First, this unsteady problem with imposed reflecting boundary conditions is run, using a CFL factor of 0.5, on a coarse and on a fine grid, composed of 600×150 and 4000×1000 points, respectively. The results are displayed in Figure 3. A graphical comparison indicates that a more accurate solution is obtained using the finer grid. Additionally, given that approximately first-order scheme is too diffusive, we perform another computation with the relaxation scheme associated with a second-order space method based on a MUSCL reconstruction used in [START_REF] Van Leer | Towards the ultimate conservative difference scheme V. A second-order sequel to Godunov's method[END_REF][START_REF] Van Leer ; Engquist-Osher | On the relation between the upwind-differencing schemes of Godunov[END_REF] associated with four different limiters: minmod [START_REF] Roe | Characteristic-based schemes for the Euler equations[END_REF], van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second order scheme[END_REF], superbee [START_REF] Roe | Characteristic-based schemes for the Euler equations[END_REF], and the positivity preserving limiter [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] Suresh | Positivity-preserving schemes in multidimensions[END_REF]. We use the coarse mesh composed of 600×150 points to show the results in Figure 4. Visually, the positivity preserving limiter is the less diffusive for this particular experiment and clearly more accurate than the first-order scheme on a fine grid.

In all experiments, along the evolution, we obtain the formation of a mushroom cap and the developing of side rolls, as expected. In general, we observe that the Rayleigh Taylor instability has been successfully reproduced.

1-D Equilibrium: Lane-Emden equation

We propose to compute the state of a self gravitational fluid at equilibrium and governed by a polytropic equation of state. In this context, equations (1.1) are rewritten in spherical coordinates (r,θ, ϕ) and we assume rotational invariance around the axes e ϕ and e θ . The system reduces to

           ∂ t (ρ)+ 1 r 2 ∂ r (r 2 ρu r ) = 0, ∂ t (ρu r )+ 1 r 2 ∂ r (r 2 ρu 2 r )+∂ r (p)+ρ∂ r (φ) = 0, 1 r 2 ∂ r r 2 ∂ r φ = 4πGρ.
(5.4)

The polytropic equation of state that governs the gas is given by

p = κρ γ = κρ 1+ 1 n ,
where κ is the polytropic constant, γ the adiabatic exponent and n the polytropic index. In this context, the 1-D equilibrium flows have been characterized by the Lane-Emden equation [START_REF] Kippenhahn | Stellar Structure and Evolution[END_REF], which will be derived now. First, assuming the fluid to be at hydrostatic equilibrium (∂ t = 0 and u r = 0), system (5.4) becomes

∂ r (p) = -ρ∂ r (φ), 1 r 2 ∂ r r 2 ∂ r (φ) = 4πGρ.
Combining the momentum equation and the polytropic relation we obtain, after space integration, a relation that links the fluid density to the potential

ρ = -φ (n+1)κ n .
Using the Poisson equation, we find the following second order ordinary differential equation for the potential:

∂ 2 r (φ)+ 2 r ∂ r (φ) = 4πG -φ (n+1)κ n .
(5.5) Defining

z = Ar, w = φ φ c = ρ ρ c 1 n , A 2 = 4πG (n+1) n κ n (-φ c ) n-1 = 4πG (n+1)κ (ρ c ) n-1 n ,
with φ c and ρ c being the potential and the density at the center of the domain, respectively, the differential equation (5.5) becomes the Lane-Emden equation [START_REF] Kippenhahn | Stellar Structure and Evolution[END_REF] 1

z 2 ∂ z z 2 ∂ z w +w n = 0. ( 5.6)
We are only concerned with solutions to the ordinary differential equation (5.6) that are finite at the center z=0. Hence, we assume that w(z,n) is a solution that fulfils the central boundary conditions w(0,n) = 1 and ∂ z w(0,n) = 0. Additionally, for n = 0,1,5, we have

w(z,n = 0) = 1- z 2 6 , w(z,n = 1) = sin(z) z , w(z,n = 5) = 1 √ 1+z 2 /3 .
We propose to validate the relaxation method with a polytropic gas of index n = 1. In order to reproduce this experiment, the numerical codes must be extended to the spherical finite volume formulation. For this, we propose to integrate the system over the spacetime domain V i ×T n , where

V i = [r i-1/2 ,r i+1/2 ]×[0,π]×[0,2π] and T n = [t n ,t n+1 ].
Using the elementary volume expression dv = r 2 sinθdrdθdϕ, we integrate the time differential terms using a finite volume approach and obtain

V i T i ∂ t (ρ)dtdv = V i (ρ(x,t n+1 )-ρ(x,t n ))dv = |V i |(ρ n+1 i -ρ n i ), V i T i ∂ t (ρu r )dtdv = |V i | (ρu r ) n+1 i -(ρu r ) n i ,
where the volume of a cell is |V i |=4π(r 3 i+1/2 -r 3 i-1/2 )/3. The integration of the divergence terms gives

T i V i 1 r 2 ∂ r (r 2 ρu r )dvdt = 4π T i r i+1/2 r i-1/2 ∂ r (r 2 ρu r )drdt = 4π∆t r 2 i+1/2 (ρu r ) L i+1/2 -r 2 i-1/2 (ρu r ) R i-1/2 ,
and

T i V i 1 r 2 ∂ r (r 2 ρu 2 r )dvdt = 4π∆t r 2 i+1/2 (ρu 2 r ) L i+1/2 -r 2 i-1/2 (ρu 2 r ) R i-1/2 ,
with (ρu) L,R i+1/2 already defined for the relations (4.10)-(4.11). For the gradient terms, we note that ρ∂ r (φ) = 0 over the integration domain V i and we obtain (5.8) To illustrate the accuracy of the relaxation scheme, we solve the Poisson equation with a standard LU method at each time step since efficiency is not required at the present time.

T i V i ∂ r (p)+ρ∂ r (φ)dvdt = 4πr 2 i T i r i+1/2 r i-1/2
In order to make numerical comparisons, the standard method is introduced. The first step is devoted to the resolution of system (5.4) with ∂ r (φ)=0 by using the proposed relaxation scheme in the spherical formulation (5.7). In the second step, we solve the following system:      ∂ t (ρ) = 0, ∂ t (ρu r )+ρ∂ r (φ) = 0, 1 r 2 ∂ r r 2 ∂ r φ = 4πGρ.

As a numerical test, we consider a polytropic gas governed by p = κρ 1+1/n with the polytropic constant κ=1000 and index n=1. The size of the experimental domain is [0,L r ] with L r = 6.0×10 5 m. At the center of the domain, we impose ρ c = 10kg.m -3 . We set the number of points to 100 and the final time of the experiment to 1.0×10 7 s. The density approximation is displayed in Figure 5.

Additionally, we observe the parasitic currents evolution during the calculation. For this, we plot the evolution in time of the L 2 norm of the velocity, in logarithmic scale, (Figure 6). The velocity field is far lower for the relaxation method than for the standard method, i.e., for the relaxation scheme, we have an approximate residual velocity of 0.01 and for the standard scheme, we get approximately 0.3. The relaxation scheme enables a better approximation of equilibrium solutions because it generates very low parasitic currents.

Conclusions

In this paper, we have derived a Suliciu-type relaxation solver to approximate the solutions to the Euler-Poisson system. In view of astrophysical simulations, the asymptotic regime of self gravitational equilibrium flows must be preserved. Although several models have been developed for this specific problem [START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | An Asymptotic Preserving scheme for the Euler equations in a strong magnetic field[END_REF], they become far too complex in 3-D or when an implicit formulation is required. We have presented here a simpler numerical scheme based on a Godunov-type solver deduced from a relaxation system to resolve accurately steady equilibrium flows. We proved the accuracy of our relaxation method on an equilibrium flow under constant gravitational forces. Then, we proposed the numerical simulation of a Rayleigh Taylor instability, an unsteady problem, which was successfully reproduced by the proposed numerical scheme. Finally, we reconsidered our scheme in spherical coordinates in order to recover precisely the solution of the Lane-Emden equation. Future work will be devoted to a well-balanced time implicit formulation of the method and simulations of 3-D astrophysical problems.

  relaxation equilibrium is reached.

Figure 1 :

 1 Figure 1: Accuracy of the relaxation method in the case of an equilibrium flow under constant gravitational field and with a fluid governed by an isentropic equation of state.

Figure 2 :

 2 Figure2: Results in logarithmic scale of the L 2 error (left) and L 2 norm of the velocity (right), as functions of time t, for Test 5.1 with 1000 grid points. Solid line: relaxation method proposed in this paper; Dashed line: standard method. Since the same half CFL condition (4.7) is used, the total number of time steps for both methods is of the same order of magnitude, specifically, 1.001×10 5 and 1.006×10 5 steps for the relaxation and standard methods, respectively.

Figure 3 :

 3 Figure 3: Rayleigh Taylor instability computed with the relaxation scheme proposed in this paper in a 2D computational domain. The results are given at times t = 2.4s, 4.0s, 5.6s, 7.2s.

Figure 4 :

 4 Figure 4: Rayleigh Taylor instability computed with a MUSCL-type second-order extension of the relaxation scheme proposed in this paper using different limiters. The computations are performed with 600×150 points and the results are given at times t = 2.4s, 4.0s, 5.6s, 7.2s.

Figure 5 :

 5 Figure 5: Numerical densities for Test 5.3 obtained by the relaxation and standard methods, compared to the exact solution ρ(r) = 10 * sin(z)/(z) with z = Ar. The total number of time steps required for this simulation is approximately 1.065×10 7 steps for both methods.

Figure 6 :

 6 Figure 6: L 2 norm of the velocity for Test 5.3, in logarithmic scale, as a function of time t, with 100 grid points. Solid line: relaxation method proposed in this paper; Dashed line: standard method.

  .37) Now, we verify that if α is large enough such that |[ψ]|<<|α|, then the value of the speed S M is approximately 1 2

  ∂ r (p)drdt, For the Poisson equation, we propose an implicit scheme deduced from the integration over the volume V i for which we obtain

	while	4πr 2 i T i	r i+1/2 r i-1/2	∂ r (p)drdt = 4πr 2 i ∆t (p) L i+1/2 -(p) R i-1/2 .
	The numerical scheme associated with system (5.4) and with the relaxation method pro-
	posed in this paper is thus given by
	          	ρ n+1 i ∆t -ρ n i (ρu r ) n+1 i ∆t + -(ρu r ) n 3r 2 i+1/2 (ρu r ) L i+1/2 -3r 2 i-1/2 (ρu r ) R i-1/2 r 3 i-1/2 i+1/2 -r 3 i +3 r 2 i+1/2 (ρu 2 r ) L i+1/2 -r 2 i-1/2 (ρu 2 = 0, r ) R i-1/2 +r 2 i (p) L i+1/2 -(p) R i-1/2 r 3 i+1/2 -r 3 i-1/2	= 0. (5.7)
		3r 2 i+1/2 i+1/2 -r 3 (r 3 i-1/2 )∆r	φ n+1 i+1 -	3(r 2 i+1/2 +r 2 i-1/2 ) (r 3 i+1/2 -r 3 i-1/2 )∆r	φ n+1 i	+	3r 2 i-1/2 i+1/2 -r 3 (r 3 i-1/2 )∆r	φ n+1 i-1 = 4πGρ n i .

Acknowledgments

This research was partially supported by the A.N.R. (Agence Nationale de la Recherche) through the projects SiNeRGHY (ANR-06-CIS6-009-01) and Anemos (ANR-11-MONU-002). The authors would also like to thank Neil Vaytet for his careful reading of the manuscript.


Mathieu Salanne
email: mathieu.salanne@upmc.fr

Leonardo J A Siqueira

Ari P Seitsonen

Paul A Madden

Barbara Kirchner

From molten salts to room temperature ionic liquids: Simulation studies on chloroaluminate systems

An interaction potential including chloride anion polarization effects, constructed from firstprinciples calculations, is used to examine the structure and transport properties of a series of chloroaluminate melts. A particular emphasis was given to the study of the equimolar mixture of aluminium chloride with 1-ethyl-3-methylimidazolium chloride, which forms a room temperature ionic liquid EMI + -AlCl - 4 . The structure yielded by the classical simulations performed within the framework of the polarizable ion model is compared to the results obtained from entirely electronic structure-based simulations: An excellent agreement between the two flavors of molecular dynamics is observed. When changing the organic cation EMI + by an inorganic cation with a smaller ionic radius (Li + , Na + , K + ), the chloroaluminate speciation becomes more complex, with the formation of Al2Cl - 7 in small amounts. The calculated transport properties (diffusion coefficients, electrical conductivity and viscosity) of EMI + -AlCl - 4 are in good agreement with experimental data.

I. INTRODUCTION

When monovalent chloride salts are dissolved in liquid AlCl 3 , ionic liquids with low melting points are formed. When the added cation is inorganic, the system becomes a molten salt, like NaAlCl 4 with T m = 152 • C . 1 When it is organic, the system is likely to be liquid below 100 • C, thus potentially becoming a room temperature ionic liquid (RTIL); for example T m (EMIAlCl 4 ) = 7 • C 2 (where EMI + stands for 1-ethyl-3-methylimidazolium cation).

A number of investigations using molecular dynamics simulations have been performed for ionic liquids directed at either structure or dynamical properties. 3 Given the molecular nature of cations or anions forming ionic liquids, earlier molecular dynamics simulations were performed with pairwise additive potentials based on nonpolarizable force fields such as OPLS or AMBER. [4][5][6][7][8][9][10][11] These simulations have been able to reproduce the structural properties of ionic liquids, for instance showing that a somewhat unusual hydrogen bond can take place in imidazolium based ionic between the imidazolium ring proton and the anion species. 12 This observation is in agreement with IR and NMR spectroscopy studies. [13][14][15] At intermediate spatial range (0.5 to 0.8 Å-1 ), a chainlength dependent pre-peak appears either in the static structure factor calculated from molecular dynamics trajectories 16,17 or obtained by small wide angle X-ray scattering (SWAXS) and/or small angle neutron scattering (SANS). [18][19][20] This feature has been assigned to structural segregation or inhomogeneity, i.e. the imidazolium rings and the alkyl chain regions of cations respectively cluster to form polar and apolar domains in ionic liquids where the cations have long alkyl chains. [START_REF] Bhargava | Balasubramanian[END_REF] Despite their capability to reproduce the structure, the non-polarizable force fields often begin to fail when predicting dynamical properties, such as diffusion coefficient, electrical conductivity, and viscosity. 22 In general the dynamical properties calculated with non-polarizable models are slower than is observed in experiments: The diffusion coefficient and electrical conductivities are lower while the viscosity is greater than the measured values. 11,[23][24][25][26] Note that pairwise additive potentials are often parameterised to implicitly include some aspects of the polarization effects. For example, the simulated fluid may become less viscous if the magnitudes of the partial charges are reduced, 27,28 but sometimes at the price of a poorer representation of other properties. These trends are particularly welldocumented in the case of the inorganic systems. 29 Disregarding the polarization method used, molecular dynamics simulations of molten salts in which polarization effects are included have shown that the liquids become more mobile. 30,31 Moreover, molecular dynamics simulations of ionic liquids including polarization have shown better agreement of dynamical properties with experiment, that is, the simulated liquids are either as mobile as the experimental sample or more mobile than the liquids simulated with non-polarizable model. [32][33][34][35][36] Among the studies using a polarizable model for ionic liquids, we would cite the recent and comprehensive work by Yan et al, 35,36 where the authors carefully investigated the role-played by polarization in the ionic liquid EMI + NO 3

-. The authors found that when polarization effects are included, the short-range interactions are enhanced due to the additional charge-dipole and dipole-dipole interactions, while long-range electrostatic interactions become more screened. Accordingly, the imidazolium ring proton -NO - 3 hydrogen bond strength is increased upon switching the polarization on. Another consequence of the stronger polar short-range interactions achieved with the polarizable model for the structure is the enhancement of the tendency towards structural inhomogeneity, which acts by pulling the apolar chains away from the polar center and reinforces the assembly of ethyl chain into the apolar domains. The attenuation of long-range electrostatic interaction is responsible for the increase of diffusion, conductivity and structural relaxation and the reduction of viscosity.

In this work, we study a series of chloroaluminate liquids. A particular emphasis is given to the room temperature ionic liquid consisting of an equimolar mixture of EMICl with AlCl 3 . We employ two flavors of molecular dynamics simulations, which differ by the method involved for the calculation of the interactions: In the first, they are obtained from a full density functional theory (DFT) calculation, and in the second an analytical force field (including polarization effects) is used. Both methods have already been used successfully for the study of pure AlCl 3 . 37,38 In a first step we detail how the parameters of the interaction potential for EMICl-AlCl 3 have been obtained. The structure and dynamics of this melt are then analyzed and compared to a series of molten salts consisting of equimolar mixtures of AlCl 3 with LiCl, NaCl and KCl.

II. INTERACTION POTENTIALS

A. Functional form

In recent work on the simulation of inorganic ionic systems, interaction potentials that were used commonly included many-body effects such as polarization. [39][40][41][42][43] The most generally retained analytical expression 44

is V = i<j q i q j r ij + B ij e -a ij r ij -f ij 6 (r ij ) C ij 6 (r ij ) 6 (1) -f ij 8 (r ij ) C ij 8 (r ij ) 8 + V pol ,
where the three first terms are pairwise additive. The first term corresponds to the electrostatic interactions, and formal integer charges are used. The two following terms are expressed by a Born-Huggins-Mayer type potential; they consist of an exponentially decaying term accounting for the overlap repulsion of the electronic clouds at short distances, and a term which represents the dispersion effects. The C ij 6 and C ij 8 are the dipoledipole and dipole-quadrupole dispersion coefficients; f n are Tang-Toennies dispersion damping functions 45 describing the short-range penetration correction to the asymptotic multipole expansion of dispersion. These functions take the form

f ij n (r ij ) = 1 -e -b ij n r ij n k=0 (b ij n r ij ) k k! (2)
where the b ij n parameter sets the range of the damping effect. The polarization term reflects the distortion of the electronic density in response to the electric fields due to all the other ions. It is written

V pol = i,j q i µ j α g ij (r ij ) -q j µ i α g ji (r ij ) T ij α (3) -µ i α µ j β T ij αβ + i ( 1 2α i | µ i | 2 )
where T α and T αβ are the charge-dipole and dipoledipole interaction tensors and α i is the polarizability of ion i. Again, Tang-Toennies functions are included to account for the short-range damping effects:

g ij (r ij ) = 1 -c ij e -b ij r ij 4 k=0 (b ij r ij ) k k! (4)
These functions differ from the previously defined f ij n due to the presence of an additional parameter c ij , which measures the strength of the ion response to the shortrange effects.

The set of induced dipoles {µ i } i∈[1,N ] is treated as 3N additional degrees of freedom of the systems. The dipoles are determined at each time step by minimization of the total polarization energy and they depend on the positions of all the atoms at the corresponding time; therefore the polarization part of the potential is a many-body term.

In the case of room-temperature ionic liquids, molecular species are involved, and additional intramolecular interactions accounting for covalent bonding therefore have to be included. Moreover, the electrostatic interactions now involve partial charges for each atom and the Lennard-Jones type is more commonly used than the Born-Huggins-Mayer one for the remaining pairwise additive interactions. 5,6,[8][9][10][11] The repulsion and damping terms are then cast together in the following way:

V LJ = i<j 4 × ǫ ij σ ij r ij 12 - σ ij r ij 6 (5)

B. Parameterization

As soon as a functional form has been chosen for the interaction potential, all the atomic and atom pair parameters have to be carefully determined. In all cases, the overall charges on the Cl -, Al 3+ , and EMI + are fixed as -1, +3, and +1 respectively. Our approach consists in parameterizing the remaining terms on the basis of first-principles electronic structure calculations. 44 This is achieved by a "force-matching" method, which is generalized by using the first-principles calculated induced dipoles in addition to the forces in the fitting procedure. Amongst molten chlorides, this method has already been applied successfully to the case of LiCl, NaCl and KCl. 46 Here we extend this set of parameters in order to include the Al 3+ and 1-ethyl-3-methylimidazolium cations. In the present work, we have chosen to describe the EMI + cation as a non-polarizable species for two reasons. Firstly, polarization effects are dominated by the chloride anions (in systems like EMI-BF 4 involving less polarizable anions the situation is opposite, with the polarization of the cation playing the main role 47 ). Secondly, from a practical point of view, the computer time associated with the minimization procedure of the polarization term increases substantially with the number of polarizable species. For all the intramolecular/intermolecular interactions involving the EMI + cation only, the parameters were taken from the well-tested force field of Pádua and co-workers. 9 The topology of the molecule is shown on figure 1, and the partial charges and Lennard-Jones individual parameters are given in table I. The pair parameters involved in equation 5 are obtained from the usual Lorentz-Berthelot mixing rules:

ǫ ij = ǫ i × ǫ j (6)
σ ij = σ i + σ j 2 (7)
FIG. 2: Results of force matching for one EMICl-AlCl3 configuration. The three components of the force acting on each aluminium and chloride ions are compared (dots: DFT, line: PIM (dots: DFT, line: PIM).

We have then applied the force-matching method to fit all the other interactions, which involve the chloride anion on the one hand and either the aluminium or the EMI + cation atoms on the other hand. In order to ensure that the force-matching really refines the cation-anion terms, we have held fixed the parameters for the Cl-Cl interactions at the values previously determined for the pure alkali halides. 46 To this end, we extracted several configurations from some trajectories obtained from Car-Parrinello molecular dynamics simulations of the liquids EMICl and EMICl-AlCl 3 . Each of these configurations were then used again as an input to a plane wave density functional theory (DFT) calculation (the details of these calculations are given in the next section). We first performed an energy minimization in order to find the ground-state electronic structure. The forces on each chloride and aluminium (if present) ions were extracted at this stage. Then we transformed the Kohn-Sham orbitals into a set of maximally localized Wannier functions, 48 from which we could determine the induced dipole moment on the Cl -anions. The parameters in the polarizable potentials are then optimized by matching the dipole and forces from the potential to the DFT calculated ones. 44,49 Two distinct functions are minimized. The first one corresponds to the polarization part:

χ 2 D = 1 N Cl - N Cl - i=1 | µ i DFT -µ i PIM | 2 | µ i DFT | 2 (8)
where the subscript PIM stands for "Polarizable Ion Model". Since the chloride anion is the only polarizable species, for a given atom type α the two short-range damping parameters b Clα and c Clα were optimized with this function. The second function,

χ 2 F = 1 N Cl - N Cl - i=1 | F i DFT -F i PIM | 2 | F i DFT | 2 (9)
+ 1 N Al 3+ N Al 3+ i=1 | F i DFT -F i PIM | 2 | F i DFT | 2
, is used to parameterize the repulsion term of the interaction potential, i.e. for determining the various B ij and a ij values in equation 2. A well known limitation of the PBE functional that we used in our DFT calculation is that it does not take into account properly the dispersion interactions, 50 so that the C ij 6 and C ij 8 cannot be fitted by this procedure. Instead, we have used our previously determined Cl --Cl -dispersion coefficients 46 for calculating the various C ij n (n = 6,8) by using the Slater-Kirkwood relationship.

The quality of the representation of the forces is illustrated in figure 2 where we compare the three components of the forces acting on the aluminium (32 first atoms) and chloride ions calculated with our potential with the DFT calculated data, for one of the EMICl-AlCl 3 configurations. The agreement is seen to be uniformly good. For both RTIL systems, we obtained values of 0.13 for χ 2 F , which is similar to our previous work on molten LiCl, NaCl and KCl. Concerning the dipoles, the values of χ 2 D are respectively 0.39 and 0.03 on average for the EMICl and EMICl-AlCl 3 configurations. The poorer agreement in the first case is due to the fact that the electric field felt by the Cl -anions, which is only due to the EMI + cation, is only partially reproduced by the partial charges of the atoms. When aluminium cations are added to the melt, this electric field becomes dominated by the trivalent charges they carry, and this effect being much more easily captured by our model. The observation suggests that further refinement of the charges of the EMI + model would be a good target to further improve the predictive power of the simulations.

The parameters which were determined from the fitting procedure are summarized in tables II and III. No distinction has been made between the various H i , C i and N i atom types because it did not lead to any improvement of the fit quality. For the ab initio-molecular dynamics simulations we used density functional theory with the generalised gradient approximation of Perdew, Burke and Ernzerhof, PBE 51 as the exchange-correlation term in the Kohn-Sham equations, and we replaced the action of the core electrons on the valence orbitals with norm-conserving pseudo potentials of the Troullier-Martins type. 52,53 We expanded the wave functions in plane waves up to the cut-off energy of 70 Ry. We sampled the Brillouin zone at the Γ point, employing periodic boundary conditions.

B ij ×10 -6 a ij C ij 6 ×10 -6 C ij 8 ×10 -6 i -j (J mol -1 ) ( Å-1 ) (J mol -1 Å6 ) (J mol -1
We performed the simulations in the NVT ensemble, employing a Nosé-Hoover thermostat at a target temperature of 300 K and a characteristic frequency of 595 cm -1 , a stretching mode of the AlCl 3 molecules. We propagated the velocity Verlet equations of motion with a time step of 5 a.t.u. = 0.121 fs, and the fictitious mass in the Car-Parrinello molecular dynamics 54 for the electrons was chosen as 700 a.u.; we note that the fictitious dynamics increases the temperature seen by the ions somewhat. 55 A cubic simulation cell with a edge length of 22.577 Å containing both 32 anionic and cationic molecules, equalling to the experimental density of 1.293 g/cm 3 . The length of the trajectory was 20 ps. The CPMD simulation code was used. [START_REF]K developers group[END_REF]

B. Polarizable ion model-based molecular dynamics

For the inorganic molten salts (LiCl-AlCl 3 , NaCl-AlCl 3 and KCl-AlCl 3 ) the simulations were performed in the NVT ensemble, employing a Nosé-Hoover thermostat at a target temperature of 573 K, with a time step of 0.5 fs, for a total simulation time of 3 ns. For EMICl-AlCl 3 , several simulations were performed at the temperatures of 298, 433, 473, 523 and 573 K, with a time step of 1 fs. Note that the higher temperatures employed may be above the thermal decomposition temperature of the ionic liquid (no data could be found in the litterature); these simulations were performed in order to make a fair comparison of the structure and dynamics with the inorganic molten salts. The total simulation time was of 1.5 ns for each temperature except 298 K, for which it was of 3 ns. The FIST module of CP2K simulation package was used. 57 In all cases, cubic simulation cells were employed. The total number of atoms and the cell edge length, which were chosen according to the experimental densities (an extrapolation of the room temperature data was made for EMICl-AlCl 3 when necessary), 2,58 , are provided in table IV. The cut-off distance for the real space of the Ewald sum and the short range potential was set to 13 Å for EMICl-AlCl 3 and to half the size of the simulation cell otherwise. Unlike previous classical MD simulations of chloroaluminate-based ionic liquids, in which the AlCl - 4 ion was treated as a molecular species, 5,59 an advantage of the two simulation methods involved in the present study is that they both treat the chloride and aluminium ions as independent species. This is done intrinsically in electronic structure based simulations on one hand, and by construction of the polarizable ion model (PIM) on the other. This capability has been used, for example, in pure AlCl 3 liquid, in which structure is characterized by the formation of AlCl - 4 tetrahedra sharing an edge to form Al 2 Cl 6 dimers or larger clusters. MD simulation studies either employing similar interaction potentials as in the present study 60 or based on electronic structure calculations 38,61 were able to confirm quantitatively this picture, in agreement with X-ray diffraction experiments. 62 The formation of large chloroaluminate species was also investigated in simulations of a system consisiting of one EMICl pair dissolved in a AlCl 3 solvent. 63 In equimolar mixtures of AlCl 3 with monovalent chloride M Cl salts, the main "chemical reaction" involving chloroaluminate species which is expected to occur [64][65][66] is:

+ = Li + , Na + , K + or EMI + ). System T (K) N M + N Al 3+ N Cl -L ( Å) LiCl-AlCl3 573
AlCl - 4 + AlCl - 4 → Al 2 Cl - 7 + Cl - (10)
Here in the case of the equimolar EMICl-AlCl 3 the two simulation methods agree on the formation of isolated AlCl - 4 anions only, in agreement with experimental results. 65 Small differences are observed in the first neighbour distances, which are shorter in the simulation employing the polarizable model: The first peaks of the Al-Cl radial distribution functions respectively appear at 2.18 Å and 2.22 Å for PIM-MD and CPMD, while the corresponding Cl-Cl functions display some first peak maxima at 3.53 Å and 3.61 Å. These differences may be attributed to the choice of keeping the Cl --Cl -interaction potential parameters similar to the values which had been previously been determined for pure LiCl, NaCl and KCl molten salts. 46 Another source of difference could be the presence of attractive dispersion effects in the classical MD only. When M + is an inorganic cation (Li + , Na + , K + ), we observe a weakening of the Al-Cl links. This may be quantified from the corresponding radial distribution functions (RDF), for which we observe a decrease of the intensity of the first peak, which switches from a value of 47 in EMICl-AlCl 3 to one of 18 in LiCl-AlCl 3 (not shown). Some Cl -ions begin to jump from one aluminium coordination sphere to another. The signature of these exchanges can also be tracked in the RDFs which are given on figure 3: The minimum occurring after the first maximum no longer corresponds to a value of zero as in the case in EMICl-AlCl 3 . The value of this minimum also increases when the ionic radius of the M + cation decreases in the series. The position of the second peak of the RDF is slightly shifted toward larger distances when passing from Li + to Na + and then K + , but this distance increases by as much as ≈ 3 Å when passing from K + to EMI + (note that the peak also becomes much wider, with a pronounced shoulder). The absence of exchange of Cl - between AlCl - 4 coordination tetrahedra in EMICl-AlCl 3 can therefore easily be explained by the longer distance that has to be crossed during a jump. Going back to the chemical reaction defined in equation 10, we observe some signature of the formation of Al 2 Cl - 7 ions in our systems by examining the Al-Al RDFs (shown on figure 4). A first peak with a small intensity is obtained in all the systems except EMICl-AlCl 3 . The corresponding distance (ranging from 3.5 Å in LiCl-AlCl 3 to 3.7 Å in KCl-AlCl 3 ) is smaller than the length of two Al-Cl bonds, indicating pairs of Al 3+ ions linked by common Cl -anions. 60 To quantify the proportion of each chloroaluminate species present in the melt, we have followed the same procedure as in our previous work on LiF-BeF 2 mixtures, where we could provide a quantitative picture of the speciation in the melt for a wide range of composition. 67,68 Here we have chosen to define an Al-Cl-Al bond when two Al-Cl and the corresponding Al-Al distances are shorter than the corresponding RDFs minima. By extending the analysis of the linkages between ions, we can calculate the proportion of any different polynuclear ionic species present in the melt. The main species that we observed apart from the isolated AlCl - 4 anions were the Al 2 Cl - 7 and Al 3 Cl - 10 species; their proportions are provided in table V for a temperature of 573 K. Note that these numbers have to be taken cautiously: Due to the finite size of the simulation cells and to the relatively low number of chloride anions jumping or bridging events, some refinement may be necessary. We also observed some small proportions of transient chemical species such as AlCl -M + radial distribution functions at 573 K. In the case where M + is the EMI + cation, the position of the center of mass of the molecule was used to calculate the function. Small arrows indicate the position of the corresponding Al 3+ -Al 3+ RDF first maximum (occulting the first peak due to the formation of Al2Cl - 7 ions); the Na + one is not visible because it is confounded with the K + one.

From the analysis of speciation in the chloroaluminate species we can try to describe the equimolar mixtures of M Cl with AlCl 3 as generic M X molten salts (like e.g. alkali halides), where M + is the monovalent cation and X -is the AlCl - 4 anion. In simple M X molten salts, the structure around a given ion consists in alternating layers of opposite charges. As a result, the RDFs are characterized by strong oscillations, with the like-like RDFs maxima (minima) which are located at the same positions as the cation-anion RDF minimum (maximum). In the present work, we can compare the M + -Al 3+ RDF to the Al 3+ -Al 3+ one; the Al 3+ was chosen because it gives the AlCl - 4 ion center of mass; for the EMI + cation the M + -Al 3+ was calculated with respect to the imidazolium ion center of mass position. These functions are shown on figures 4 and 5. On the latter we also report the position of the Al 3+ -Al 3+ RDF first maximum (omitting the small peak due to the formation of Al 2 Cl - 7 ) with an arrow on the abcissa in order to facilitate the comparison with the position of the first minimum in the cation-anion term. We observe that these arrows are not located at the same position as the corresponding minima, but the difference remains small. This means that the structure deviates slightly from that of a simple M X molten salt. This conclusion still holds for EMICl-AlCl 3 at lower temperatures: The main features of the M + -Al 3+ RDF remain unchanged when passing from 573 K to 298 K.

From this figure we can observe somewhat different behaviour depending on the size of the M + cation. In the case of the LiCl-AlCl 3 melt, the M + -Al 3+ RDF minimum is located at a shorter distance than the corresponding Al 3+ -Al 3+ RDF first maximum; this is due to the small size of the Li + cation, which allows it to "penetrate" slightly into the AlCl - 4 anion van der Waals sphere. The difference becomes smaller for the NaCl-AlCl 3 and KCl-AlCl 3 melts for which the extrema of the two functions almost coincide, and the situation is reversed for EMICl-AlCl 3 where closer packing of the anions is allowed. The M + -Al 3+ RDF first peak is also broader in the case of EMICl-AlCl 3 compared to the alkali cations containing systems, but the difference is not dramatic. The two latter observations suggest that in the present system the EMI + cation behaves mostly like a monoatomic cation with a very large ionic radius. In the next section we will investigate its coordination structure on a more local scale, in order to understand if the charge delocalization as well as the presence of different atom types induce a particular ordering in the first solvation shell.

C. EMICl-AlCl3 equimolar mixture

In order to characterize further the intermolecular structure of the EMICl-AlCl 3 system, the radial distribution functions of the chloride anions and the atoms from the EMI + cations have been calculated. They are plotted on figures 6 (hydrogen atoms) and 7 (carbon atoms). The CPMD and PIM results obtained at room temperature are compared; an excellent agreement is observed for all the atoms. The small differences observed are more likely due to the poorer sampling in the CPMD simulation rather than to a deficiency of the classical interaction potential. In a previous attempt to obtain classical force field parameters through a force-matching procedure similar to ours on the ionic liquid dimethylimidazolium chloride (DMICl), Youngs et al. have been less successful in reproducing the structure obtained from CPMD. 69 Although they were able to reproduce the peak positions better than with previous classical MD potentials, the intensities of the first peaks where largely overestimated (especially in the case of the ring protons). The overstructuring obtained by these authors could partly be cancelled by multiplying all the ǫ i parameters by a factor of 2, at the price of a worse reproduction of the initial set of DFT forces (therefore leading to higher χ 2 F values). The better agreement observed here can be attributed to two factors. The main one is the use of a more complicated functional form for the analytical interaction potential, or in other words to the inclusion of anion polarization effects. The second is the fact that we focused on the chloride (and aluminium) ion forces in our force-fitting process; in Youngs et al. study all the parameters, including the intramolecular ones, were fitted during the procedure. The intramolecular contribution may have somewhat "masked" the intermolecular ones, thus leading to a smaller precision for the latter. It is worth noticing that the better agreement observed here was not obtained by fitting more parameters. Due to the choice of using the set of partial charges of Padua and coworkers for the EMI + ion atoms, only twelve parameters involving the Cl -and EMI + species had to be fitted in the parameterization process: The eight Cl-C i ,N i ,H i ,Al repulsion parameters (B ij and a ij ) and the four parameters b ij and c ij involved in the damping of the chloride ions induced dipoles by the hydrogen and aluminium atoms. A striking feature of the radial distribution functions involving the chloride ions and the EMI + atoms is that the most intense peaks are observed for the carbon atoms. The peak is more pronounced for the C 7 and C 8 , while its intensity is almost the same for the other carbons. This result contrasts sharply with some other RTILs, like the above-mentioned DMICl for example, for which the anions tend to stay in the vicinity of the imidazolium ring protons. 69,70 On the contrary, it resembles to the case of the 1-n-butyl-3-methylimidazolium hexafluorophosphate, for which the hydrogen bonds were found to be weaker than expected, as indicated by their short lifetimes. 71 In the present system the situation is due to the huge electrostatic interaction between Cl -anions and the (Lewis) acidic Al 3+ cation. To illustrate this the electrostatic potential mapped onto the isosurface of the electron density of the two individual ions is shown in figure 8. It provides insight into the charge distribution of each ion. For the AlCl - 4 we recognize in the red color (low electrostatic potential) the consequence of the nega-tive charge that is associated with this ion. The negative charge is distributed all over the chlorine atoms, with a higher negative charge on the inner side (the spherelike surface around the chloride is coloured in dark red on the interior and in lighter red at the exterior of the AlCl - 4 ion). Towards the center, close to the aluminium atom, we find a decreased negative charge shown as the blue color in the left panel of Fig. 6. The opposite is the case for the EMI + . Here the positive charge leads to the blue color (high electrostatic potential) around this ion. Upon closer inspection we find that the ring protons show the same blue color as most of the molecule. A slight decrease of the charges can be found in the methyl group protons and a stronger decrease (green color) can be found at the terminal ethyl protons. This is in accordance with the observation from the radial pair distribution functions and with chemical intuition that these protons are less acidic than the other protons of the cation.

In the PIM simulations, the chloride anions therefore have their induced dipoles oriented along the Al-Cl axis, which hinders the formation of hydrogen bonds. If we look in further detail at the chloride-proton radial distribution functions, we see that the H 2 -Cl -one shows the most pronounced peak, associated to the shortest separation. It is therefore clear that this proton is the most acidic coordination site. The functions involving the other two protons of the imidazolium ring (H 4 and H 5 ) also display some peaks at slightly larger distances. However, considering the protons from the ethyl and the methyl group, small pronounced peaks can also be observed. The methyl group hydrogen (H 8 ) function first maximum has the same height as the H 4 one; a similar observation is made for the terminal hydrogen atoms of the ethyl group (H 6 ) compared to the ring H 5 protons.

The arrangement of the chloride anions the imidazolium cations therefore seems to correspond more to an involved network including the alkyl chains as well as the ring of the imidazolium cation, rather than to individual ion pairs, confirming the results obtained for many imidazolium-based ionic liquids. 71,72 In order to confirm this picture, the three-dimensional density map was calculated for the chloride anions located in the first solvation shell of the EMI + . Views of these maps perpendicular and parallel to the imidazolium ring are provided on figure 9. Red, green and blue regions correspond respectively to high, intermediate and low probability of finding an anion around a given cation. By high probability we mean higher than 80 % of the maximum density of occurences in these figures, and low probability means in between 40 % and 60 % of the maximum density of occurences (probability lower than 40 % is not shown). The important features of the radial distribution functions, i.e. a preference for some regions close to the C 7 and C 8 carbons and an absence of important correlations with the protons, are recovered. In general, localization effects are much less important than in the DMICl case, and the formation of hydrogen bond is not the main influence governing the structure of EMICl-AlCl 3 . It is worth mentioning that in a study involving imidazolium cations with a series of halide anions of different size, the analysis of density maps showed that the larger the anion the lower is the probability of the existence of hydrogen bond between the anion and the cation protons. 8 The present results agree with this view when considering that the relevant anion in our system is the full AlCl - 4 entity. The structural analysis of the equimolar EMICl-AlCl 3 mixture shows that it can almost be viewed as a classical M X molten salt, where M + = EMI + and X -= AlCl - 4 , with a structure characterized by a competition between close packing and charge ordering effects occuring between the two species. Small deviations from the M X system are observed due to the existence of weak spatial correlations arising from the formation of hydrogen bonds between the ring protons and the chloride anions.

D. Transport

In the recent simulation work on the room temperature ionic liquids, it was shown that polarization effects are important for reproducing accurately the transport properties. [32][33][34]36,47 Here from the Einstein relation

D α = lim t→∞ 1 6t | δr i (t) | 2 (11)
where δr i (t) is the displacement of a typical ions of species α in time t, we obtained a diffusion coefficient of 1.04×10 -10 m 2 s -1 for the cation at 298 K. This compares very well with the experimental value measured by pulse-field gradient NMR spectroscopy, which is 0.95×10 -10 m 2 s -1 at 293 K. 73 For the anion we obtained a diffusion coefficient of 4.45×10 -10 m 2 s -1 at 298 K. From the mean squared displacement of the overall charge density it is also possible to calculate the electrical conductivity of a melt according to the following formula:

σ = e 2 k B T V lim t→∞ 1 6t | i q i δr i (t) | 2 , ( 12
)
where e is the elementary charge and q i the magnitude of the (partial or formal) charge on atom i. Notice that this expression involves the correlations between the displacements of different ions, which is essential in this case as the motion of the Al 3+ and Cl -ions within a given AlCl - 4 is strongly correlated. The viscosity is calculated by integration of the correlation function of the stress tensor:

η = 1 k B T V ∞ 0 σ αβ (0)σ αβ (τ ) dτ, (13)
where σ αβ is one of the components of the stress tensor. The two latter quantities, which are collective transport coefficients (compared to the diffusion coefficients, which are individual), require much longer simulation times to get converged values. Such quantities may therefore only be calculated with the analytic force-field model, they are well out of reach of full ab initio simulations. In the case of EMICl-AlCl 3 , we could extract their values for all the temperatures except at 298 K. They are shown on figure 10 in the form of an Arrhenius plot. On the same figure we also show the experimental values; 2 although those were obtained at lower temperatures, we observe a good match between the two sets of data, which confirms that our polarizable ion model yields the correct dynamics of the system. At the temperature of 573 K, which is above the melting point of the inorganic chloroaluminates, it is also possible to compare the diffusion coefficients of the various ions in all the systems. They are reported in table VI. For all the species, we observe an increase of the diffusion coefficient with the M + cation size. Note that the number density also changes when switching from one system to another, which is likely to be the main cause for this evolution. Still it is possible to come to the conclusions that in all systems the chloroaluminate anion diffuses more slowly than the corresponding M + cation. We can also conclude that the organic cation-based ionic liquid, EMICl-AlCl 3 , fits into the same pattern of behaviour as the inorganic cations-based compounds in regard to the variation of transport properties with cation size.

V. CONCLUSION

We have shown in this study that the polarizable ion model developed for simple molten salts 46 can successfully be transfered to chloroaluminate ionic liquids. The interaction potential was parameterized purely from firstprinciples through a force and dipole-matching procedure, and a common set of parameters could be obtained for the mixtures of AlCl 3 with several inorganic salts (LiCl, NaCl and KCl) as well as with the organic 1-ethyl-3-methylimidazolium chloride. An advantage of our model is that it allows us to quantify the chloroaluminate speciation due to the choice of treating Al 3+ and Cl -as independent ionic species.

A particular emphasis was given to the study of the equimolar mixture EMICl-AlCl 3 , for which the structure yielded by the classical simulations performed within the framework of the polarizable ion model is compared to the results obtained from entirely electronic structurebased simulations: An excellent agreement between the two flavors of molecular dynamics is obtained, which enables us to conclude that the structure of EMICl-AlCl 3 can almost be viewed as the one of a classical M X molten salt, where M + = EMI + and X -= AlCl - 4 . The structure is characterized by a competition between close packing and charge ordering effects occuring between the two species. Small deviations from the M X system are observed due to the existence of weak spatial correlations arising from the formation of hydrogen bonds between the ring protons and the chloride anions. When decreasing the cation ionic radius (EMI + → K + → Na + → Li + ) small proportions of Al 2 Cl - 7 due to the reaction

AlCl - 4 + AlCl - 4 → Al 2 Cl - 7 + Cl - (14)
start to be observed (Al 3 Cl - 10 is also observed in equimolar LiCl-AlCl 3 ). An interesting extension of this work would be to study the same series of systems under non-stoechiometric conditions. The influence of the M + cation chemical nature in the Lewis acidity of the melt can also be quantified using a recently developed method. 74 In future work we will try to perform such calculations in order to relate it to the chloroaluminate speciation determined in the present study.

The only many-body effect included in the interaction potential for EMICl-AlCl 3 is the anion polarization, which enabled us to perform long enough runs to get the dynamic properties of the system, which are out of reach of first-principles simulations. The calculated transport properties were compared to the experimental data, showing again a good agreement for both the individual (diffusion coefficients) and collective (electrical conductivities and viscosities) dynamics of the system.

FIG. 1 :

 1 FIG.1: Topology of the 1-ethyl-3-methylimidazolium molecular ion.

FIG. 3 :

 3 FIG.3: Al 3+ -Cl -radial distribution functions for the M Cl-AlCl3 mixtures at 573 K. Insert: Enlargment of the first minimum zone.

FIG. 4 :

 4 FIG.4: Al 3+ -Al 3+ radial distribution functions for the M Cl-AlCl3 mixtures at 573 K.

FIG. 6 :FIG. 7 :

 67 FIG.6: Radial distribution functions of the chloride anions and the protons from the EMI + cations. Top: CPMD results; bottom: PIM results. Left: Protons from the imidazolium ring. Right: Protons from the alkyl chains.

FIG. 8 :

 8 FIG. 8: The electrostatic potential mapped onto the electron density with a surface values of 0.067 e -. Å-3 . The colour scale ranges from -0.1 (red) to +0.1 (blue) atomic units. Left: AlCl - 4 , right: EMI + .

FIG. 9 :

 9 FIG.9: Probability density maps of location of nearestneighbor anions around the imidazolium ring in equimolar EMICl-AlCl3 mixture at 298 K. The left panel shows the view above the imidazolium ring plane, and the bottom panel shows the view on the ring plane. Red, green and blue regions correspond respectively to high, intermediate and low probability of finding an anion around a given cation. By high probability we mean higher than 80 % of the maximum density of occurences, and low probability means in between 40 % and 60 % of the maximum density of occurences (probability lower than 40 % is not shown).

FIG. 10 :

 10 FIG.10: Viscosities (top) and electrical conductivities (bottom) of EMICl-AlCl3 in the form of an Arrhenius plot. The simulated values, although obtained in a different temperature regime, compare well with the experimental ones.2

TABLE I :

 I Interaction potential parameters (partial charges, Lennard-Jones parameters) used for the EMI + cation.9 Note that the ǫ i and σ i parameters were used for the cation-cation van der Waals interactions only.

	Atom i	ǫ i	σ i	q i
		(kJ mol -1 ) ( Å) (e)
	N1	0.71128 3.25 0.15
	C2	0.29288 3.55 -0.11
	N3	0.71128 3.25 0.15
	C4	0.29288 3.55 -0.13
	C5	0.29288 3.55 -0.13
	C6	0.27614 3.50 -0.17
	C7	0.27614 3.50 -0.05
	C8	0.27614 3.50 -0.17
	H2	0.12552 2.42 0.21
	H4	0.12552 2.42 0.21
	H5	0.12552 2.42 0.21
	H6	0.12552 2.50 0.13
	H7	0.12552 2.50 0.06
	H8	0.12552 2.50 0.13

TABLE II :

 II Parameters for the pairwise additive terms of the interaction potential. b ij 6 = b ij 8 = 3.21 Å-1 for all pairs. Atom pair

TABLE III :

 III Parameters for the polarization part of the interaction potential. The polarizability of the Cl -, Na + and K + ions were respectively set to 2.96, 0.13 and 0.70 Å3 .

	Atom pair b ij	c ij	c ji
	i -j	( Å-1 ) (-) (-)
	Cl-Al	3.802 2.953 -
	Cl-Li	3.517 2.079 -
	Cl-Na	3.326 3.000 0.697
	Cl-K	3.084 3.000 0.917
	Cl-Hi	2.901 0.301 -
	III. SIMULATION DETAILS

A. Electronic structure-based molecular dynamics

TABLE IV :

 IV Number of atoms and cell edge length (L) employed in the simulations (M

TABLE V :

 V Proportions (given in percentage of aluminium ions) of the various chloroaluminate species at 573 K.

	System	[AlCl -4 ] [Al2Cl -7 ] [Al3Cl -10 ]
	LiCl-AlCl3	94.1	3.5	0.7
	NaCl-AlCl3	95.7	3.4	0.2
	KCl-AlCl3	98.2	1.6	0.0
	EMICl-AlCl3 100.0	0.0	0.0
	B. Effect of the M + cation nature on the
	intermolecular structure
	FIG. 5: AlCl -4

TABLE VI :

 VI Diffusion coefficients of the equimolar M Cl-AlCl3 systems at 573 K (M + = Li + , Na + , K + or EMI + ). All values are in 10 -10 m 2 .s -1 System D M + D Al 3+ D Cl -

	LiCl-AlCl3	16.4 6.19 6.28
	NaCl-AlCl3	20.6 6.63 6.65
	KCl-AlCl3	24.3 8.62 8.62
	EMICl-AlCl3 35.5 23.2 23.2

Acknowledgments

MS is grateful to the FAPESP and the Federal University of São Paulo for a travel grant which enabled this collaboration. LJAS also thanks FAPESP, proc. 08/08670-7.