S T Éphane Brull
email: stephane.brull@math.u-bordeaux1.fr

Jacques Schneider
email: jacques.schneider@univ-tln.fr

DERIVATION OF A BGK MODEL FOR REACTING GAS MIXTURES

Keywords: Kinetic theory, gas mixtures, BGK models, entropy minimization, hydrodynamic limit subject classifications 35Q20, 35Q35

In this paper we derive a new relaxation model for reacting gas mixtures. We prove that this model satisfies the fundamental properties (equilibrium states, conservation laws, H-theorem, ...). We also consider the slow reaction regime. In this case a rigorous Chapman-Enskog procedure is performed and Navier-Stokes equations are derived.

Introduction

The modeling of reacting gas mixtures is essentially based on the microscopic phenomena that occur when particles are colliding. Part of those collisions are non reactive while the others lead to different species and feature the chemical reactions. The same distinction is done at the kinetic level by adding to the classical collision operators for gas mixtures collision operators that only feature chemical reactions. A quite general setting of those equations can be found in [START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. It is however interesting to focus on the case where only one chemical reaction is considered. A simplified model has been proposed by Rossani and Spiga [START_REF] Rossani | Spiga A note on the kinetic theory of chemically reacting gases[END_REF] where each species has only one degree of internal energy and where elastic collisions occur for the non reactive part of the collision operators. A generalization of this model to different degrees of internal energies was proposed in [START_REF] Groppi | Spiga Kinetic approach to chemical reactions and inelastic transitions in a rarefied gas[END_REF]. In particular they have considered a four species gas mixture with a reversible chemical reaction

A 1 + A 2 ↔ A 3 + A 4 .
(1.1) (Here A i stands for the species i).

The main motivation of this paper is the derivation of a relaxation model based on the simplified model of Rossani and Spiga. Our concerns are to include as much physical properties as possible (positivity, equilibrium states, conservation laws, H theorem) as well as minimizing computational cost. Such a model is the well known BGK model [START_REF] Bhatnagar | A model for collision processes in gases[END_REF] for monoatomic gas. But its drawback is that it cannot lead to the proper Prandtl number at the Navier-Stokes level. This problem has been overcome by Holway [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF] who introduced the so-called ES-BGK model. When inert gas mixtures are considered cross effects such as "mass" diffusion (Fick law) or thermo diffusion (Soret law) occur. The derivation of BGK models giving the right transport coefficients at the hydrodynamic limit is much more complicated. For example Kosuge [START_REF] Kosuge | Model Boltzmann equation for gas mixtures: Construction and numerical comparison[END_REF] has derived a BGK operator that is able to approximate those coefficients but which guarantees neither the non negativeness of the distribution functions nor the entropy decay. Andries & co-authors have successfully addressed the later problems [START_REF] Andries | Perthame A consistent BGK-type model for gas mixtures[END_REF] but their model is valid only for Maxwellian molecules and thus does not respond to the above requirements. The authors of the present paper have introduced a new concept of relaxation coefficients on non-conserved moments which allows to construct BGK type models. A first achievement was to recover the ES-BGK model in the monoatomic and polyatomic cases [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model Cont[END_REF][START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF]. But this idea extends to more general cases and we have derived in [START_REF] Brull | Derivation of BGK models for mixtures[END_REF] a relaxation model allowing both to recover the correct Fick law at the hydrodynamic and enjoying the required mathematical (and physical) properties.

Coming back to reactive gas mixtures many BGK models were constructed by Groppi, Spiga and co-authors (see for example [START_REF] Bisi | Spiga On a kinetic BGK model for slow chemical reactions[END_REF] for a complete bibliography) basing upon the "simplified" Boltzmann equation of Rossani and Spiga [START_REF] Rossani | Spiga A note on the kinetic theory of chemically reacting gases[END_REF]. Those models always feature only one relaxation operator per species that includes both mechanical and chemical reactions. For the whole set of models the mechanical part is based on [START_REF] Andries | Perthame A consistent BGK-type model for gas mixtures[END_REF] where Maxwellian molecules are considered. Hence the hydrodynamic limits in the case of slow reactions lead to the Navier-Stokes equations with uncorrect transport coefficients. The chemical contributions depend on the approximation of the chemical collision operators. There are essentially two types of models. The first ones are constructed using estimations of the exchanges of mass, momentum and (kinetic plus internal energy) per species [START_REF] Groppi | A Bhatnagar-Gross-Krook approach for chemically reacting mixtures[END_REF] and [START_REF] Bisi | Spiga On a kinetic BGK model for slow chemical reactions[END_REF]. So in some sense those models mimic the one by Andries & al [START_REF] Andries | Perthame A consistent BGK-type model for gas mixtures[END_REF]. However those values cannot be computed in general even under the assumption of Maxwellian molecules and strong restrictions must be made -distribution functions at mechanical equilibrium (Maxwellian distributions)to overcome this problem. This means that those models are essentially adapted to the case of slow reactions where the chemical reactions occur when the distribution of molecules are closed to mechanical equilibrium. The second model may be called "the" BGK model for mechanical plus chemical reactions [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]. The macroscopic parameters of the Maxwellian attractors are given under the requirements of number atoms conservation and conservation laws and bound together with some type of mass action law. This model enjoys the classical H-theorem for chemically reacting gas mixtures contrarily to the first class of models for which it is very unlikely that this property holds.

We intend in this paper to take advantage on one side on the modeling for inert gas mixtures and on the other side on the modeling for chemical reactions. Our methodology is based on a splitting between a mechanical operator and a chemical one as for the full Boltzmann equation. This procedure is particulary adapted to the case of slow reaction regime. More precisely we simply add to the model derived in [START_REF] Brull | Derivation of BGK models for mixtures[END_REF] a chemical relaxation operator directly derived from [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]. However remark that while the chemical BGK operator seems to be more adapted to the case of fast chemical reactions, our choice relies upon the possibility to prove a H-theorem. In doing so the present paper must be considered more as an example of the splitting methodology than an attempt to match the best approximation of the chemical contribution in the case of slow reaction regime.

The paper is organized as follows. In section 3 we make a short review of different Boltzmann equations for reacting mixtures. In the next section we show that the construction of the mechanical model can be achieved either by using experimental values (Fick law) or by algorithms developed by Ern and Giovangigli [START_REF] Ern | Multicomponent Transport Algorithms[END_REF]. Special attention is given as to the definition of the relaxation coefficients for the chemical model. Then we show that the whole relaxation model satisfies the H-theorem and that equilibrium states are Maxwellian functions which densities and temperature are bound together with the mass action law. Finally section 5 is devoted to the derivation of the Euler and the Navier-Stokes system in the slow reacting regime. The Chapman-Enskog procedure yields a first formulation of the Navier-Stokes. We focus on the calculation of the chemical terms at the first order. The second formulation of those equations (section 5.3.2) is obtained by using the work of Kurochkin, Makarenko and Tirskii [START_REF] Kurochkin | Tirskii Transport coefficients and the Onsager relations in the kinetic theroy of dense gas mixtures[END_REF] on the transport coefficients and the Onsager relations. Thus the Navier-Stokes equations are recast in the classical framework of polyatomic gases (see for example [START_REF] Giovangigli | Multicomponent flow modeling[END_REF]). This makes possible the comparison between the different coefficients (species multicomponent diffusion, thermal diffusion, etc) that are obtained on one side with our model and on the other side with the usual theory. This also allows to determine entirely the definition of our model as already mentioned above.

2. Notations Consider a gas mixtures with p components.

f i (t,x,v) (or for short f i ,i ∈ [1,p] with f := (f 1 ,••• ,f p )) represents the distribution function of a given species i.
We denote with n i , ρ i , u i , E i and T i the macroscopic quantities representing respectively the density and mass per unit volume, average velocity, internal energy per unit volume and finally temperature of a given species i. They are defined by the following relations:

n i = R 3 f i dv, ρ i = m i n i , n i u i = R 3 vf i dv, E i = n i 3 2 k B T i + E i = R 3 m i 2 v -u i 2 + E i f i dv,
where k B is the Boltzmann constant and is the Euclidean norm in R 3 . E i represents the internal energy of species i due to chemical links. In the same way macroscopic quantities for the mixture are defined by

n = p k=1 n k , ρ = p k=1 ρ k , ρu = p k=1 ρ k u k , ρ 2 u 2 + E = 4 i=1 1 2 ρ i u i 2 + n i E i , E = 4 i=1 n i 3 2 k B T + E i .
(2.1)

The Boltzmann equation for reactive gas mixtures

In this section we recall some backgrounds about the Boltzmann equations for reacting mixtures. In the first subsection we recall the Boltzmann equation for inert gas mixtures which will be considered here. Next in subsection 3.2 we review some collision operators for reacting mixtures.

Inert monoatomic gas mixtures.

For the sake of simplicity we are going to consider the mechanical Boltzmann equations for inert gas mixtures of p components, that is when collision between molecules are elastic. It reads [START_REF] Aoki | layer for a gas mixture[END_REF]:

∀i ∈ [1,p], ∂ t f i + v • ∇ x f i = k=p k=1 Q ki (f k ,f i ) := Q i (f,f), where Q ki (f k ,f i ) = R 3 ×S 2 (f k (w * ki )f i (v * ki ) -f k (w)f i (v)) σ ik (ω ω ω.V, V ) V dwdω ω ω.
Here Q ki is the Boltzmann collision operator between molecules of species i and k and σ ik = σ ki is the differential cross section which depend on the interaction potential between species i and k. Finally V = wv is the relative velocity. The post collisional velocities are given by

v * ki = v -2 m k m i + m k ((v -w) •ω ω ω)ω ω ω, w * ki = w + 2 m i m i + m k ((v -w) •ω ω ω)ω ω ω,
where m i represents the particle mass of species i. These equations satisfy the conservation of mass (per species) momentum and energy at microscopic level. The set of collisional invariants K is spanned by the set of functions (φ l ) l∈{1;8} defined by

    1 0 0 0     ,     0 1 0 0     ,     0 0 1 0     ,     0 0 0 1     ,     m 1 v α m 2 v α m 3 v α m 4 v α     ,     1 2 m 1 v 2 + E 1 1 2 m 2 v 2 + E 2 1 2 m 3 v 2 + E 3 1 2 m 4 v 2 + E 4     , (3.1)
for α = x,y,z.

Collision equilibria are given by Maxwellian distributions

f i = M i = n i m i 2πk B T 3 2 exp - m i 2k B T (v -u) 2 , (3.2)
for any densities n i ≥ 0, mean velocity u ∈ R 3 and positive temperature T . k B is the Boltzmann constant.

Next we need to introduce the space L 2 (M) equipped with its natural dot product:

Ψ,Φ = i=p i=1 R 3 ψ i φ i M i dv. (3.3)
The i th component of the linearized Boltzmann operator reads [START_REF] Aoki | layer for a gas mixture[END_REF] L

Bi (g) = 1 M i   p j=1 Q ji (M j ,M i g i ) + Q ji (M j g j ,M i )   . (3.4)
Then Ker (L B ) = K, L B is continuous, invertible and self adjoint negative on K ⊥ .

Reacting gas mixtures.

In this subsection we recall the three main different approaches leading to collision operators for gas mixtures.

In [START_REF] Rossani | Spiga A note on the kinetic theory of chemically reacting gases[END_REF] the authors consider a mixture of four gases undergoing a reversible bimolecular chemical reaction together with mechanical binary collisions. Then the kinetic equation writes for the i th species

∂f i ∂t + v • ∇f i = Q i + J i , i ∈ {1;4}, (3.5)
where Q i is the elastic collision operator and J i represents the chemical collision term. Let us describe briefly J i . Let m i (resp. E i ) be the particle mass (resp. the energy of chemical link) for a given species i. Then J 1 is defined by

J 1 (f ) = R 3 S 2 U (V -ε 12 )V I 34 12 (V,Ω.Ω )(( m 12 m 34 ) 3 f 3 (v )f 4 (w ) -f 1 (v)f 2 (w)) (3.6) dwdΩ ,
where (v ,w ) are the post-collisional velocities and m ij = mimj mi+mj is the reduced mass. The quantity ε ij is defined by

ε ij = 2 ∆E m ij with ∆E = E 3 + E 4 -E 1 -E 2 .
U denotes the unit step function and is a threshold for endothermic reactions. This term is linked to the energy of activation which must be given to the system in order that the endothermic reaction can happen. I 34 12 (V,Ω.Ω ) is the cross section where

V = w -v , Ω = (w -v)/V, Ω = (w -v )/ w -v .
The other chemical terms J i for i ∈ {2;3;4} can be obtained after suitable permutations.

Densities are not conserved but only the total mass. However one must also consider the conservation of atoms number during a reaction. In the case of bimolecular interactions this may be expressed by the following equalities (see [START_REF] Rossani | Spiga A note on the kinetic theory of chemically reacting gases[END_REF])

J i (f )dv = -J j (f )dv, (i,j) = (1,3), (1,4), (2,4), (3.7)
from which one also has

J 1 (f )dv = J 2 (f )dv, J 3 (f )dv = J 4 (f )dv. ( 3.8)
Finally the conservation of momentum and total energy (φ l ) l∈{5;8} must be added to this set of invariants. A generalization of this operator for many levels of energy per particles ( [START_REF] Groppi | Spiga Kinetic approach to chemical reactions and inelastic transitions in a rarefied gas[END_REF]) can be considered in the framework of Giovangigli and Ern ([20,[START_REF] Giovangigli | Multicomponent flow modeling[END_REF]) and references therein. Each distribution function reads f i (t,x,v,I) where I is the index of the quantum internal energy states of the i th . Then the chemical operator for f 1 reads

J i (f ) = J,K,L f 3 f 4 β 3K β 4L β 1I β 2J -f i f j W IJKL 1234 dwdv dw . (3.9)
J,K,L are all possible internal quantum energy states of species 2,3,4. In the case of [START_REF] Groppi | Spiga Kinetic approach to chemical reactions and inelastic transitions in a rarefied gas[END_REF] the ratio β 3K β 4L /(β 1I β 2J ) does not depend upon I,J,K,L and simply takes the value (m 12 /m 34 ) 3 while in [START_REF] Giovangigli | Multicomponent flow modeling[END_REF]

β iI = a iI m 3 i .
is the Plank constant and a iI is the degeneracy of the internal energy state I for species i. Within our framework we have on one side a iI = 1 and β iI = /m 3 i . However does not appear in (3.6) since we are considering bimolecular reaction. For the sake of simplicity we will set β iI = 1/m 3 i . Remark 3.1. The main difference between (3.6) and (3.9) is the interpretation of W IJKL 1234 which are probabilities of transition [START_REF] Ern | Giovangigli The kinetic chemical equilibrium regime[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. However those values as well as differential cross sections I 34 12 (V,Ω.Ω ) are unknown in general. But while angular deviation of the velocities are not accessible from experiments quite accurate evaluations of cross sections (as functions of the translational energy) have been carried out (see e.g [START_REF] Mccourt | Kuscer Nonequilibrium Phenomena in Polyatomic Gases[END_REF]). They allow to compute the concentrations of the different species at any quantum energy states (vibrational, rotational) as soon as those quantum scates can be evaluated for the molecules to be considered. This means that even under strong deviation from equilibrium (except for the translational energies) the chemical processes are quite well evaluated. Remark 3.2. One may consider the problem of finding relevant differential cross sections such as to obtain a correct hydrodynamical limit. In this vein one may cite the work of Desvillettes & co-authors [START_REF] Desvillettes | Salvarani A Kinetic Model Allowing to Obtain the Energy Law of Polytropic Gases in the Presence of Chemical Reactions[END_REF][START_REF] Bourgat | Perthame Microreversible collisions for polyatomic gases and Boltzmann's theorem[END_REF] who proposed a model for both mechanical and chemical reactions that allows to recover the energy law of polytropic gases at the Euler limit. The distribution function

f i (t,x,v) (or for short f i ,i ∈ [1,4] with f := (f 1 ,f 2 ,f 3 ,f 4
)) of a given species i evolves according to the kinetic equation:

∀i ∈ [1,4], ∂ t f i + v • ∇ x f i = R M E i (f ) + R CE i (f ) := R i (f ), (4.1)
where

R M E i (f ) (resp. R CE i (f ))
represents the mechanical (resp. chemical) part of the relaxation operator.

The mechanical BGK operators

The mechanical part of the model is the one constructed in [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]. For the sake of clarity we recall how it is derived. Let us now consider the formalism of the thermodynamic of irreversible processes. Mass flux for the i th species reads

J i = j=4 j=1 L ij ∇ G tr i + L i5 ∇ - 1 T (4.2)
where the specific Gibbs functions G tr i are defined by

G tr i = k B T m i ln n i Q tr i , with Q tr i = (2πm i k B T ) 3 2 . ( 4.3)
We may as well consider the full specific Gibbs functions by including partition functions for internal energies. But the present definition is more suited to the derivation of the hydrodynamic limit of our model. We will later on use the full specific Gibbs functions to recast the Navier-Stokes equations in a classical way (section 5.3.2). Remark here that L ij and L i5 are of opposite sign of the traditional thermodynamic coefficients. The reason is that they are directly derived from the Boltzmann equation through a Chapman-Enskog expansion. Notice that the density fluxes are generally rather expressed in terms of the phenomenological coefficients in which case the dependance of the above fluxes on gradients of the densities is called the Fick law. Whatever is the expression the corresponding coefficients depend on temperature T and also on the number densities and molecular masses of each component.

Then in order to recover the matrix (L ij ) i,j=1,...,4 one performs the following steps: 1. Define the symmetric non positive matrix L * whose elements are

L * ij := L ij √ ρ i ρ j T . (4.4)
In the case where the exact Fick coefficients F ij can be measured from experiments one sets

L * ij = n j k B T F ij (m i m j ) 3 2 √ n i n j . (4.5)
It is of course necessary that those coefficients are such that the matrix L * remains symmetric non positive. L * is also related to the diffusion matrix (D ij ) i,j=1,...,4 which definition is recalled in section . Then one also has

L * ij = - √ ρ i ρ j nk B T D ij , i,j = 1,...,4. (4.6)
This second formula makes possible the computation of L * by using the algorithms developed by Ern and Giovangigli [START_REF] Ern | Multicomponent Transport Algorithms[END_REF]. 2. Remark that L * always diagonalizes in an orthonormal basis:

L * = W T K * W
and that up to some permutation in W and K * the corresponding eigenvalues (k * r ) r are non null for r = 1,...,3 while k * 4 = 0 and corresponds to the conservation of the total mass ρ. Then set

λ r = -k * r -1 for r = 1,2,3 and λ 4 = 0.
3. Define the vector of velocities U = (u 1 ,...,u 4 ) T with the relation

U -U = N -1 W T I - 1 ν M ∆ W N U -U , (4.7)
where U = (u,u,u,u) T , U = (u 1 ,...,u p ) T and N and ∆ are the diagonal matrix which diagonal terms are respectively ( √ ρ 1 ,..., √ ρ 4 ) and (λ 1 ,...,λ 4 ).

Set

T = T - 1 3nk B W T I - 1 ν M ∆ W N U -U 2 . (4.8)
Definition 4.1. The i th component of the mechanical operator reads

R M E i (f ) = ν M (G i -f i ), i ∈ [1,4], (4.9
)

with ∀i ∈ [1,p], G i = n i (2πk B T /m i ) 3/2 exp - m i (v -u i ) 2 2k B T , (4.10)
where ν M > 0 represents the relaxation coefficient. This coefficient must be chosen with the constraint ν M ≥ max r λ r /2 in order to ensure the positivity of T * . 5. Definition of ν M : In our case ν M must be considered as a mean frequency of return to thermodynamic equilibria of the whole set of molecules. Assume that the shear viscosity η is given by some formula or approximated by the algorithms of Ern and Giovangigli. Then if η satisfies the condition

η ≤ nk B T max r λ r (4.11)
we set

ν M = nk B T /η. (4.12)
This definition allows to recover the viscosity η at the hydrodynamic limit.

It is also classical for the BGK model with one component. However we also have and the definition (4.12) is subjected to the relative speed between the shear stress and diffusion phenomena. For want of anything better we set

λ r = nk B T d * r , r = 1,
ν M = max nk B T η ,max r λ r , (4.13)
in such a way that the model is always well-posed.

4.1.2. The chemical BGK operators R CE i
represents the chemical part of the relaxation operator constructed in [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]. Definition 4.2. The i th component of the chemical operator reads

R CE i (f ) = ν C i ( Mi -f i ), (4.14)
where (ν C i ) i are the chemical collision frequencies and

Mi = ñi m i 2πk B T 3 2 exp - m i 2k B T (v -ũ) 2 . (4.15)
It is clear that within the above framework the collision frequency of e.g the first species should read (see [START_REF] Groppi | A Bhatnagar-Gross-Krook approach for chemically reacting mixtures[END_REF] or [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF])

1 n 1 R 3 S 2 U (V -ε 12 )V I 34 12 (V,Ω.Ω )f 1 (v)f 2 (w))dwdΩ , (4.16)
(if the chemical collision integral is (3.6)). However neither the differential cross section nor the functions f are known. Recall that the characteristic time of mechanical collisions is supposed to be much smaller than that of the chemical ones. In this case a possible approximation of the chemical terms (J i (f)) i (3.6) consists in replacing f with the set of functions at thermodynamical equilibrium M. Further simplification -Maxwellian molecules [START_REF] Groppi | A Bhatnagar-Gross-Krook approach for chemically reacting mixtures[END_REF] -allows to evaluate the collision frequencies which then read

ν C 1 = ν 34 12 2 √ 2π Γ 3 2 , ∆E k B T n 2 , ( 4
.17)

ν C 3 = ν 34 12 2 √ 2π Γ 3 2 , ∆E k B T m 12 m 34 3 2 exp ∆E k B T n 4 . (4.18)
Here ν 34 12 is a physical constant depending on the species involved in the chemical reaction but not on the specified internal energy states (E i ) i=1,...,4 which are considered here. Γ(α,x) is the uncomplete Gamma function defined by

Γ 3 2 , ∆E k B T = +∞ ∆E k B T exp -t t 1 2 dt = ∆E k B T 1 2 exp -∆E k B T +erfc ∆E k B T .
Thus N ∆E is somehow the activation energy of the reaction (1.1) when only one internal energy state per molecule is considered (N is the Avogadro number). Notation 4.1. The set of functions M = ( M1 , M2 , M3 , M4 ) is denoted M(f) and respectively Mi (f ) for each i. We will omit the dependance on f when it is self evident.

We recall how the parameters ñi , ũ and T are computed. The BGK models must satisfy the same laws than the chemical operators (3.6). That is the conservation of atoms number (3.7) together with the momentum and energy conservation. [START_REF] Andries | Perthame A consistent BGK-type model for gas mixtures[END_REF][START_REF] Bhatnagar | A model for collision processes in gases[END_REF], [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF][START_REF] Bhatnagar | A model for collision processes in gases[END_REF],

(R CE i (f ) + R CE j (f ))dv = 0, (i,j) = (1,3),
4 i=1 m i vR CE i (f )dv = 0, 4 i=1 1 2 m i v 2 + E i R CE i (f )dv = 0. (4.19)
There are not enough equations to define all unknown densities (ñ i ) i but it is convenient to write them in the following form:

ñi = n i + Λ i ν C 1 ν C i (ñ 1 -n 1 ), i = 2,3,4, (4.20)
where (Λ i ) i is the string of stoichiometric coefficients (1, 1, -1, -1). Thus if ñ1 is known all other densities are defined as well. Then the conservation of the mean velocity together with (4.20) allow to express ũ in term of the macroscopic fields:

ũ = 4 i=1 ν C i m i n i u i 4 i=1 ν C i m i n i . (4.21)
Finally the conservation of the total energy together with (4.20), (4.21) gives the following equation for the temperature T T =

4 i=1 ν C i n i 1 2 m i (u i ) 2 -ũ2 + 3 2 k B T i + ν C 1 ∆E(ñ 1 -n 1 ) 3 2 k B 4 i=1 ν C i n i . (4.22)
Thus T is also defined as soon as ñ1 is. This last quantity is obtained by solving the implicit equation

ν C 3 ν C 4 ñ1 (ν C 2 n 2 + ν C 1 (ñ 1 -n 1 )) ν C 2 (ν C 3 n 3 -ν C 1 (ñ 1 -n 1 ))(ν C 4 n 4 -ν C 1 (ñ 1 -n 1 )) exp(- ∆E k B T (ñ 1 ) ) = m 12 m 34 3 2 . (4.23)
As mentioned by the authors "the left hand side is a monotonically increasing function of ñ1 ranging from 0 to +∞ when ñ1 varies on its domain defined by the constraint of positivity of density and temperature fields. This guarantees existence and uniqueness of the solution to this equation" [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF].

Properties of the model

It is clear that those BGK models share the same collision invariants as the reacting Boltzmann equation itself. Now we show that the model satisfies an entropy dissipation property together with the H theorem. Let H(f ) be defined by

H(f ) = 4 i=1 R 3 f i ln f i m 3 i dv. ( 4
p i=1 R 3 R M E i (f ) ln(f i /m 3 i )dv ≤ 0, p i=1 R 3 R CE i (f ) ln(f i /m 3 i )dv ≤ 0. (4.25)
As a consequence H defined in (4.24) is a Lyapounov function for the whole model. Proof. On one side side we have R M E i (f )dv = 0 for all i so that

4 i=1 R 3 R M E i (f ) ln(f i /m 3 i )dv = 4 i=1 R 3 R M E i (f ) lnf i dv ≤ 0
as proved in [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]. On the other side [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF] it holds that

4 i=1 R 3 R CE i (f ) ln(f i /m 3 i )dv ≤ 0.
So Proposition 4.3 follows Theorem 4.4. The following assertions are equivalent: i) the entropy production rate is equal to 0

4 i=1 R 3 (R M E i (f ) + R CE i (f )) ln(f i /m 3 i )dv = 0. (4.26)
ii) For all i,

R M E i (f ) + R CE i (f ) = 0, (4.27)
iii) the distribution functions f i are at mechanical and chemical equilibrium:

∀i ∈ [1,4],f i = M i = n i (2πk B T /m i ) 3 2 exp - m i (v -u) 2 2k B T (4.28)
with the mass action law

n 1 n 2 n 3 n 4 = m 12 m 34 3/2 exp ∆E KT . ( 4

.29)

Proof. i) ⇒ iii) According to the above proposition 4.3 i) implies that

4 i=1 R 3 R M E i (f ) ln(f i /m 3 i )dv = 0.
This holds [START_REF] Brull | Derivation of BGK models for mixtures[END_REF] if and only if there exists macroscopic values n 1 ,...,n 4 ,u,T such that ∀i ∈ {1,2,3,4},

f i = M i = n i (2πk B T /m i ) 3 2 exp - m i (v -u) 2 2k B T .
However this condition is necessary but not sufficient for the entropy production rate to vanish. Plugging the above value of the distribution functions into the definition 4.2 of R CE i , we obtain

4 i=1 R 3 ν C i Mi -M i ln(M i /m 3 i )dv = 0.
A direct computation (see [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]) together with (4.23) gives

4 i=1 R 3 ν C i Mi -f i ln( Mi /m 3 i )dv = ν 1 (ñ 1 -n 1 ) ln ñ1 ñ2 ñ3 ñ4 ( m 3 m 4 m 1 m 2 ) 3 2 - ∆E k B T = 0, ∀f ≥ 0. (4.30)
Subtracting the two previous equations with

f i = M i gives 4 i=1 R 3 ν C i Mi -M i ln Mi M i dv = 0.
The convexity of the function x → xlnx yields f i = M i = Mi for all i ∈ {1;4} which satisfy the mass action law (4.23).

ii) ⇒ i) is trivial. iii) ⇒ ii). Assume that (f i ) i = (M i ) i such that density and temperature fields satisfy the mass action law

n 1 n 2 n 3 n 4 exp - ∆E k B T = m 12 m 34 3/2 . ( 4
= T + ν C 1 ∆E(ñ 1 -n 1 ) 3 2 k B 4 i=1 ν C i n i . (4.32)
One obtains the same equation as (4.31) by setting ñ1 = n 1 into (4.23). Thus n 1 is the unique solution of (4.23) and ñi = n i , ∀i according to (4.20). One also has T = T from (4.32) so that Mi = M i , ∀i and R CE i (f ) = 0.

Hydrodynamic limit for slow reactions

The slow reactions regime corresponds to the situation where the chemistry characteristic time is much smaller than the mechanical ones. This regime has been studied in ( [START_REF] Ern | Multicomponent Transport Algorithms[END_REF]) by estimating the ratio between the characteristic time of mechanical collisions and the chemistry characteristic time. This ratio is shown to be less than 0.16 for temperature varying between 300 and 1000 K. The hydrodynamic limit of our model is then obtained by studying the equations form

∀i ∈ [1,4], ∂ t f ε i + v • ∇ x f ε i = 1 ε R M E i (f ) + R CE i (f ), (5.1)
as ε → 0.

Asymptotic expansion The classical Chapman-Enskog expansion consists in setting

∀i ∈ [START_REF] Andries | Perthame A consistent BGK-type model for gas mixtures[END_REF][START_REF] Bhatnagar | A model for collision processes in gases[END_REF],

f ε i = f 0 i + εf 1 i + εf 2 i .... (5.2)
and plugging those expressions in (5.1). Then we get at order -1 in ε R M i (f ) = 0 for all i. So f 0 i = M i (see [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]). Setting f 1 i = M i g i the equations (5.1) at order 0 read:

∂ t M i + v • ∇ x M i = M i L M E i (g) + R CE i (M), ∀i, (5.3)
where L M E i is the i th component of the linearized operator of R M E around the Maxwellian distributions M. Its i th component is defined as usual by the formula:

DR M E i (M).Mg = M i L M E i (g) = lim τ →0 R M E i (M(1 + τ g)) -R M E i (M) τ , (5.4)
where DR M E (M).Mg represents the differential of R M E at M in the direction Mg. After computing ∂ t M i from the Euler system (5.9, 5.10, 5.11) in function of the space derivatives, the equations (5.3) will in turn be particularly easy to solve thanks to the exact formulation of the (pseudo-)inverse of L M E . In order to obtain an asymptotic expansion of (5.1) up to order 1 in ε, the term R M E i (f ) has to be expanded up to order 2. Set f 2 i = M i h i then the equation (5.1) up to the order 1 is derived from

∂ t (M i (1 + εg i )) + v • ∇ x (M i (1 + εg i )) = 1 ε R M E i M(1 + εg + ε 2 h) + R CE i (M(1 + εg)).
(

The Taylor expansion of

R M E i M(1 + εg + ε 2 h) gives R M E i M(1 + εg + ε 2 h) = R M E i (M) + εL M E i (g) + ε 2 L M E i (h) + 1 2 D 2 R M E i (M).(Mg,Mg) + O(ε 3 ). (5.6)
Next one must expands R CE i (M(1 + εg)) as well. The linearized operator L CE (g) around M is defined with

DR CE i (M).Mg = L CE i (g) = lim τ →0 1 τ R CE i (M(1 + τ g)) -R CE i (M) = lim τ →0 1 τ Mi (M(1 + τ g) -Mi (M) -τ g i M i .
(5.7)

Remark that we have not factorized M i in front of L CE i in this definition. The reason is that such a factor should not appear as suggests the last line of the above expression. Hence the expansion up to order 1 in ε of (5.5) gives

∂ t (M i (1 + εg i )) + v • ∇ x (M i (1 + εg i )) = L M E i (g) + R CE i (M) + ε(L CE i (g) + L M E i (h) + 1 2 D 2 R M E i (M).(Mg,Mg)) + O(ε 2 ).
(5.8)

Remark 5.1. According to their definition, the second order terms in ε of (5.6) conserve mass, total momentum and energy. Therefore they will neither contribute to the Euler system nor to the Navier-Stokes system.

The Euler system

The Euler equations are classically obtained by integration of the 8 first moments of (5.3). They read:

∂ t ρ i + ∇ x • (ρ i u) = m i Λ i ν C 1 (ñ 1 -n 1 ), (5.9)
∂ t (ρu) + ∇ x • (ρu ⊗ u + pI) = 0, (5.10) .11) where the total energy and the pressure are

∂ t E tot + ∇ • (E tot u + Pu) = 0, ( 5
E tot = 1 2 ρu 2 + E, P = nk B T I
and I is the identity matrix.

5.3. The Navier-Stokes system.

Setting of the problem

In order to obtain the Navier-Stokes system we have to perform the following steps:

1. compute g solution to the equations (5.3) by using the Euler system, 2. compute the linearized operator L CE by using formula (5.7), 3. integrate (5.8) with respect to the elements of K (3.1) and keep all terms up to the order 1 in ε. Those steps are performed in the following sections and the final result is given in Proposition 5.1 Proposition 5.1. The Navier-Stokes system for slow chemical reactions of (5.1) reads:

∂ t ρ i + ∇ • (ρ i u) + ε∇ • J i = m i Λ i ν C 1 (ñ 1 -n 1 ) + εω 1 i , (5.12)
∂ t (ρu) + ∇ • (ρu ⊗ u + P) + ε∇ • J u ) = 0, (5.13)

∂ t E tot + ∇ • (E tot u + Pu) + ε∇ • (J u [u]) + ε∇ • J q = 0, (5.14)
where the fluxes are defined by

J i = j=4 j=1 L ij ∇ G tr i T , J u = - nk B T ν M D(u), J q = -κ∇T + i E i J i , (5.15)
and the partial thermal conductivity for monoatomic gases

κ = 5k 2 B T 2ν M 4 i=1 n i m i . ( 5

.16)

The matrix L ij is defined in (4.2) and the Specific Gibbs functions (G tr i ) i in (4.3). D(u) is the traceless part of the deformation tensor

D(u) = 1 2 ∇u + (∇u) T - 1 3 (∇ • u) I. (m i Λ i ν C 1 (ñ 1 -n 1
)) i are the zero order chemical terms where ñ1 is the unique solution to the equation (4.23). Finally the first order perturbation of the chemical term reads

ω 1 i = m i Λ i ν C 1 τ (ñ 1 )h(g), τ (ñ 1 ) = m 12 m 34 3 2 ∆E k B T 2 / ξ (ñ 1 ) + ξ(ñ 1 ) ν 1 (∆E) 2 3 2 k 2 B T 2 4 i=1 ν C i n i , (5.17) h(g) = 1 ν M ∆E(ñ 1 -n 1 ) 4 i=1 (ν C i ) 2 ρ i 4 i=1 ν C i n i ν C 1 ν C i n i - 1 n .
(5.18)

Remark 5.2. When all chemical relaxation frequencies are equal the first order chemical perturbation ω 1 i vanishes. Moreover for fixed (t,x) the densities (ñ i ) i are solutions to the local mass action law

ñ1 ñ2 ñ3 ñ4 exp(- ∆E k B T (ñ 1 ) ) = m 12 m 34 3 2
together with the relations

ñi -n i = Λ i (ñ 1 -n 1 ), i = 2,3,4, T = T + ∆E(ñ 1 -n 1 ) 3 2 k B n .
In a thermodynamic framework -that is considering the space homogeneous problem related to (5.12) and (5.14) with initial conditions n i (0) = n i , ∀i and T (0) = T -those solutions correspond to the "true" chemical equilibrium states (when t → +∞) of the chemical reaction (1.1). Remark 5.3. According to the expression of the fluxes given in Proposition 5.1, the shear viscosity η and partial thermal conductivity for monoatomic species κ are equal to

η = nk B T ν M , κ = 5k 2 B T 2ν M 4 i=1 n i m i .
Define the mean molar mass of the mixture m as

m = N ρ 4 i=1 ρ i mi = N ρ n .
Then the Prandtl number P r in the monoatomic setting can be computed as

P r = 5 2 R m η κ = n 2 ρ 4 i=1 n i mi ,
where R = N k B is the constant of perfect gases. Hence in the situation of indifferentiability -that is when all molecules are the same -the Prandtl number is equal to 1 which is a drawback of the classical BGK model for a single monoatomic species.

Other formulation of the Navier-Stokes equations

It is clear that the present model does not pretend to reproduce all phenomena that occur for polyatomic gas mixtures. However it is possible to recast the above equations in a more familiar and useful framework such as that of the book of Giovangigli [START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. We first introduce the specific enthalpy and specific Gibbs free energy for each species. They are defined by

h i = 5 2 k B T m i + E i m i , G i = k B T m i ln n i Q i , Q i = (2πm i k B T ) 3 2 exp - E i k B T ,
where we have introduced the partition function for internal energy

Q int i = exp(-E i /k B T
) of the species i. Those definitions are similar to those introduced in [START_REF] Giovangigli | Multicomponent flow modeling[END_REF] when the molecules have multiple internal energy states. From this the mass fluxes can be conveniently rewritten in the following form

J i = 4 j=1 Lij ∇ G j T + Li5 ∇ - 1 T (5.19)
where J i is defined in (4.2) and

Lij = L ij , i,j = 1,...,4, Li5 = 4 j=1 L ij E j m j , i = 1,...,4. (5.20)
Then the heat flux J q (5.15) may be written in a "symmetric" form

J q = 4 i=1 L5i ∇ G j T + L55 ∇ - 1 T , (5.21)
where L5i = Li5 and

L55 = 4 i,j=1 E i m i L ij E j m j -κT 2 < 0, (5.22)
(κ is given in (5.16)). The null space of L is R(

T and L is symmetric non positive as shows the following decomposition L = P t LP + M with

M =       0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 L55       , P =        1 0 0 0 E1 m1 0 1 0 0 E1 m2 0 0 1 0 E3 m3 0 0 0 1 E4 m4 0 0 0 0 1        , L =       L11 L12 L13 L14 0 L21 L22 L23 L24 0 L31 L32 L33 L34 0 L41 L42 L43 L44 0 0 0 0 0 1      
.

We are now going to use the work of Kurochkin & co-authors [START_REF] Kurochkin | Tirskii Transport coefficients and the Onsager relations in the kinetic theroy of dense gas mixtures[END_REF] that bridges the TIP formulations of the fluxes with more classical and practical ones. This comparison is still possible in our case since the thermodynamic functions are clearly defined. Denote with A ∈ R 5×5 the matrix

A =       1 0 0 0 -h 1 0 1 0 0 -h 2 0 0 1 0 -h 3 0 0 0 1 -h 4 0 0 0 0 1       .
Next we define with L = A T LA. Then it is easily checked that the mass and heat fluxes read also

J i = - p j=1 ρ i D ij ∇p j p -ρ i θ i ∇T T , ∀i ∈ {1;4} (5.23)
J q = -λ∇T -p 4 i=1 θ i ∇p i p + 4 i=1 h i J i . (5.24)
Here (p i = n i k B T ) i are the partial pressures and the diffusion matrix (D ij ) i,j=1,...,4 , thermal diffusion coefficients (θ i ) i=1,...,4 and partial thermal conductivity λ are respectively found to be

D ij = - nk B L ij ρ i ρ j , θ i = Li5 ρ i T , λ = L55 T 2 .
(5.25)

Hence we recover the structural form of the fluxes in terms of phenomenological coefficients. Remark that the relation between the diffusion matrix and (L ij ) i,j=1,...,4 allows to approximate the later by using the algorithm developed by Ern and Giovangigli ( [START_REF] Ern | Multicomponent Transport Algorithms[END_REF]). As a consequence the mechanical model R M E (section 4.1.1) is made explicit once the shear viscosity η is known. This is also given by the above approximation.

Nevertheless the thermal diffusion coefficients and partial thermal conductivity are uncorrect as it can be seen from the formulae (5.20) and (5.22). Indeed the values of ( Li5 ) i and ( L55 ) i are not the "real" value coming from the thermodynamics of irreversible processes since our mechanical model just allows to recover the matrix (L ij ) i,j=1,...,4 at the hydrodynamic limit.

Let us now compare J u in (5.15) with the general form that holds for a mixture of polyatomic gases with chemical reactions which writes

J u = p reac I -α∇ • uI -ηD(u),
where p reac is the chemical pressure, α is the volume viscosity and η the shear viscosity.

In our case, since α = 0 we are only considering elastic collisions. But it is possible to recover the exact shear viscosity under suitable conditions (see section 4.1.1). Finally there is no chemical pressure in our case. This is due to the (extremely) limited number of possible quantum energy states per molecules. Finally the zero-order molecular production rate of the i th species (5.12) is quite different from the Maxwellian production rate of the chemical operators (3.6). That is

ν C i (ñ i -n i ) = Λ i R 3 J 1 (M)dv.
Only the negative parts can be compared if the frequencies are chosen such as in (5.2).

Computation of g

Before we set the first result we need to introduce some notations. Definition 5.2. C i is the vector which i th component is vu and others are 0. P K is the orthogonal projection on K with respect to the scalar product (3.3) and I is the identity operator.

Lemma 5.3. The perturbation g reads

g = (L M E ) -1 (Ψ M E ) - 1 ν M Ψ CE = g M E + g CE , (5.26)
where

Ψ M E i =   j=p j=1 k -1 B (I -P K )C j • ∇ µ j k B T + A : D(u) + B • ∇ 1 k B T   i , (5.27) Ψ CE i = -ν C i Mi -M i M i + ν C i ñi -n i n i + 2 3nk 2 B T 2 ( 1 2 m i (v -u) 2 - 3 2 k B T )ν C 1 (ñ 1 -n 1 )∆E.
(5.28)

A and B are the list of tensors defined by

(A) i = m i k B T (v -u) ⊗ (v -u) - 1 3 (v -u) 2 I , (B) i = (v -u) 5 2 k B T - 1 2 m i (v -u) 2 .
I is the identity and P K is the operator of projection on K with respect to the scalar product (3.3). Finally C i is the vector in R 3×4 which i th line is vu and the others are 0.

Proof. Using the Euler equations (5.9, 5.10, 5.11) the time derivatives of the macroscopic fields in (5.3) can be expressed in term of spatial derivatives. Then we use the formalism of [START_REF] Brull | Derivation of BGK models for mixtures[END_REF] which can be easily extended to the case of slow chemical reactions. This reads:

L M E (g) = j=p j=1 k -1 B (I -P K )m i (C j ) • ∇ G tr j + A : D(u) + B • ∇ 1 k B T -(I -P K )R CE (M).
(5.29)

Let us first give more details on

R CE i (M) = ν C i ( Mi -M i ),
which corresponds to the chemical contribution. ñi , ũ and T are defined by the relations (4.20, 4.21, 4.22) corresponding to the fields of M. Hence we have

ñi = n i + Λ i ν C 1 ν C i (ñ 1 -n 1 ), i = 2,3,4, (5.30)
ũ = u, T = T + ν C 1 ∆E(ñ 1 -n 1 ) 3 2 k B 4 i=1 ν C i n i , (5.31)
where ñ1 is the unique solution of the implicit equation (4.23). It is then an easy task to compute de i th line of (I -P K )R CE (M) by using the definition of K ⊥ and the result reads:

Ψ CE i = -ν C i Mi -M i M i + ν C i ñi -n i n i + 2 3nk 2 B T 2 ( 1 2 m i (v -u) 2 - 3 2 k B T ) ×ν C 1 (ñ 1 -n 1 )∆E.
We now need to compute

L M E -1 (Ψ M E i + Ψ CE i ). Recall that ∀φ ∈ K ⊥ , (L M E ) -1 (φ) = 1 ν M ((R -I C ) -1 • P C + (P C -I C ))(φ), (5.32)
(Lemma 8 in [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]) where C is the space spanned by the vectors

(I -P K )(C i ), i = 1,...,4
and R is a linear operator acting on C. Then one has

g M E = (L M E ) -1 (Ψ M E i
) from which one can compute the "elastic" part of the fluxes (see [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]). Next recall that

(I -P K )C i = j a ij C i , (a ij ) i,j ∈ R 4 × R 4 ,
and one can easily check that Ψ CE ∈ C ⊥ . So Lemma 5.3 follows. Remark 5.4. The decomposition of the perturbation g into a mechanical and a chemical part corresponds to the usual decomposition that is valid for the reactive Boltzmann equation given in [START_REF] Giovangigli | Multicomponent flow modeling[END_REF].

5.4. Computation of L CE . In order to compute L CE (g) from formula (5.7) we need a preliminary result. Lemma 5.4. Let Mg i be the "attractors" of the chemical term

R CE i (M(1 + εg)) then Mg i (v) = ñg i m i 2πk B Tg 3 2 exp -m i (v -ũg ) 2 2k B Tg , (5.33)
with

ñg i = ñi + εΛ i ν C 1 ν C i τ (ñ 1 )h(g) + O(ε 2 ), (5.34)
ũg = u + ε 4 i=1 ν C i m i n i u g i 4 i=1 ν C i m i n i , n i u g i = R 3 M i g i vdv, u g = (u g 1 ,u g 2 ,u g 3 ,u g 4 ), (5.35) Tg = T + εh(g) ν C 1 ∆E τ (ñ 1 ) 3 2 k B 4 i=1 ν C i n i + 1 + O(ε 2 ), (5.36)
where

T = T + ν C 1 ∆E(ñ 1 -n 1 ) 3 2 k B 4 i=1 ν C i n i , (5.37)
and

h(g) = 1 ν M ∆E(ñ 1 -n 1 ) 4 i=1 (ν C i ) 2 ρ i 4 i=1 ν C i n i ν C 1 ν C i n i - 1 n .
This tedious and lengthy proof is left to appendix. Then expanding Mg i around ε and using formula (5.7) we obtain

L CE i (g) = ν C i M i Λ i ν C 1 ν C i τ (ñ 1 ) ñi h(g) - h(g) 2πk B T + m i 2πk B T 2 v -u,u g + (v -u) 2 2k B T h(g) -ν C i M i g i .
(5.38) Remark 5.5. h(g) = 0 (5.18) when the collision frequencies (ν C i ) are equal. Hence in this case there are no contribution of order ε to the chemical term in the Navier-Stokes equations for densities evolutions. Remark 5.6. The expression (5.38) of L CE i (g) does not depend on the mechanical perturbation g M E (5.26). Then contrarily to [START_REF] Giovangigli | Multicomponent flow modeling[END_REF] there are no derivatives of the macroscopic fields in the first order production chemical terms. Proof. (Proposition 5.1). In order to obtain a Navier-Stokes system we have to consider for any i ∈ {1,2,3} the moments of M i (1 + εg i ). Writing g i = g M E i + g CE i we remark that fluxes of order ε depend as expected only on g M E i . Hence we can directly infer that the left hand sides of the Navier-Stokes equations are those obtained in ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]). Denote that one important point is the computation of matrix (L ij ) appearing in the definition of the mass flux J i . More precisely, using (5.32) we obtain L ij through the expression of the Ohnsager matrix from the Chapmann-Enskog expansion. Next one has to compute the chemical terms up to order ε. As concerns the equations for densities those terms read

R 3 (R CE i (M) + εL CE i (g))dv = Λ i ν C 1 (ñ 1 -n 1 ) + εΛ i ν C 1 τ (ñ 1 )h(g).
Recall that

R CE i (M(1 + εg)) = R CE i (M) + εL CE i (g) + O(ε 2 )
and that R CE (f ) satisfies equations (4.19) for any function f = (f 1 ,f 2 ,f 3 ,f 4 ). Then we have on one side

i m i R 3 (R CE i (M) + εL CE i (g))vdv = O(ε 2 ),
and on the other side

i R 3 ( 1 2 m i v 2 + E i )(R CE i (M) + εL CE i (g))dv = O(ε 2 ),
which end the proof.

Discussion and conclusion

The BGK model which is proposed here features the same structure as the simplified Boltzmann equation for reacting gas mixture proposed by Rossani and Spiga [START_REF] Rossani | Spiga A note on the kinetic theory of chemically reacting gases[END_REF] or of the more general Boltzmann equation with different level of internal energies [START_REF] Giovangigli | Multicomponent flow modeling[END_REF]. That is a sum of two different terms for the nonreactive and the reactive "collisions". Let us recall that our model is designed here to satisfy good mathematical properties (H-theorem, ...) and to present a splitting method that seems naturally adapted to the case of slow chemical reactions. The Chapman-Enskog expansion clearly outlines the respective contributions of the first order mechanical and chemical perturbations (5.26). As expected all transport coefficients (thermal diffusion, viscosity, Fourier law, ...) depend only on the modeling of the mechanical part of the BGK model. Its construction is based on the knowledge of the Fick law (4.5) or of the diffusion matrix (4.6) together with the shearn viscosity. In some cases the Fick law may be obtained from experiments but the corresponding matrix of coefficients must be symmetrical non positive. This property is satisfied if the diffusion matrix is approximated by algorithms developed by Ern and Giovangigli [START_REF] Ern | Multicomponent Transport Algorithms[END_REF]. The hydrodynamic limit of the model features good structural agreement (section 5.3.2) with the usual Navier-Stokes for polyatomic gases [START_REF] Giovangigli | Multicomponent flow modeling[END_REF] (section 5.3.2) features a good structural agreement. The correct diffusion matrix is recovered by construction. This holds also true for the shear viscosity under suitable condition (4.11). But the thermal diffusion coefficients and partial thermal conductivity are uncorrect. There is no volume viscosity as expected since we are only dealing with elastic collisions and molecules with only one fixed internal energy state. As a consequence there is also no chemical pressure. The zero order chemistry source terms (5.12) are not equal to the Maxwellian production rates of the real chemical operators contrarily to what is obtained from the theoretical BGK model of Bisi and Spiga [START_REF] Bisi | Spiga On a kinetic BGK model for slow chemical reactions[END_REF] and for the Grad approximation in [START_REF] Bisi | Spiga Grad's closure in the kinetic theory of a chemically reacting gas[END_REF]. But positiveness of temperature of the attractors and H-theorem are guaranteed at the kinetic level contrarily to the model in [START_REF] Bisi | Spiga On a kinetic BGK model for slow chemical reactions[END_REF]. Remark that perturbations of order ε are free from any derivatives of the spatial gradients as in [START_REF] Bisi | Spiga Grad's closure in the kinetic theory of a chemically reacting gas[END_REF]. One should also point out that those perturbations vanish when the chemical frequencies are equal. Finally a generalization of the mechanical operator as well as a construction of a chemical operator allowing to recover the Maxwellian production rates in the general case of polyatomic gases are in progress. ν C i m i n i .

Hence by setting

T g i = m i 3k B n i R 3 (v -u -εu g i ) 2 M i (1 + εg i )dv, (7.3)
we finally obtain the following expression for Tg ,

Tg = 4 i=1 1 2 ν C i m i n i (u + εu g i ) 2 -(u + εf (u g )) 2 + 4 i=1 ν C i 3 2 k B T g i + ν C 1 ∆E(ñ g 1 -n 1 ) 3 2 k B 4 i=1 ν C i n i . ( 7.4)
In order to simplify (7.4) that

4 i=1 1 2 ν C i m i n i (u + εu g i ) 2 -(u + εf (u g )) 2 = 4 i=1 ν C i m i n i 2εu.u g i -2εu.f (u g ) + ε 2 u g i 2 -ε 2 f (u g ) 2 = 2εu. 4 i=1 ν C i m i n i u g i -f (u g ) 4 i=1 ν C i m i n i + ε 2 4 i=1 ν C i m i n i u g i 2 -f (u g ) 2 = ε 2 4 i=1 ν C i m i n i u g i 2 -f (u g ) 2
according to the definition of f (u g ) (7.2). The relation (7.3) leads to

T g i = m i 3k B n i R 3 (v -u) 2 M i (1 + εg i )dv -2εu g i . R 3 (v -u)M i (1 + εg i )dv + ε 2 (u g i ) 2 R 3 M i (1 + εg i )dv .
Hence by setting Next by introducing τ (ñ 1 ) defined in (5.17), ñg 1 writes ñg 1 = ñ1 + ετ (ñ 1 )h(g) + O(ε 2 ). Hence according to (7.1) ñg i satisfies for all i ∈ {2;3;4} the relation (5.30). Moreover by plugging the expression of ñg 1 into (7.6) it holds that Tg can be defined from relation (5.36) Computation of h(g). Firstly as

h g i = m i 3k B n i R 3 (v -u) 2 M i g i dv, (7.5
R 3 (v -u)(v -u) 2 M i dv = 0, R 3 A i (v -u)(v -u) 2 M i dv = 0, R 3 B i (v -u)(v -u) 2 M i dv = 0,
the mechanical part g M E i of g i can be removed in (7.5) and h g i reduces to

h g i = m i 2k B n i R 3 M i (v -u) 2 g CE i dv.
That is

h g i = m i 2ν M k B n i R 3 (v -u) 2 ν C i ( M i -M i ) -ν C i M i (ñ i -n i ) - ν CE 1 ∆E 3nk 2 B T 2 (ñ 1 -n 1 )( 1 2 m i (v -u) 2 - 3 2 k B T ) M i dv = 3ν C i m i 2ν M n i ñi ( TM (ñ 1 ) -T ) - ν C i ν M nk B ∆E (ñ 1 -n 1 ) = ν C i ν M k B ∆E (ñ 1 -n 1 ) ν C 1 ν C i n i - 1 n
where we have used (5.37) in the last equation. Finally we obtain h(g) = ∆E(ñ 1 -n 1 )

4 i=1 (ν C i ) 2 ρ i ν M 4 i=1 ν C i n i ν C 1 ν C i n i - 1 n .

4. The relaxation model 4 . 1 .

 41 Definition of the model.

  2,3, and λ 4 = 0 where (d * r ) r are the eigenvalues (up to some permutation) of the matrix ( √ ρ i ρ j D ij ) ij (4.6). Then the condition (4.11) also reads η ≤ min r=1,2,3 d * r

. 31 )

 31 One has firstly G i = M i from (4.7) and (4.8) so that R M E i (M) = 0, ∀i (see Definition 4.1). Secondly equations (4.21) and(4.22) give ũ = u and T

5. 5 .

 5 Computation of the right-hand side of Navier Stokes system. This subsection is devoted to the end of the proof of Proposition 5.1.

m i n i u g i 4 i=1

 4

  .[START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]

	Proposition 4.3. The relaxation operators R M E i	(f ) and R CE i (f ) satisfy the following
	dissipative properties

We first want to establish a transcendental equation as ( [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]) whose ñg 1 is a solution. By using the mass conservation for i ∈ {2;3;4}

Next the conservation of the impulsion for the chemical operator gives

So defining u g i by

it comes that

Finally, (7.4) leads to

with

Next we aim to compute ñg 1 . Let us write the implicit equation (4.23) in the following way

where

The computation of ñg 1 is performed through the implicit function theorem at ε = 0. If ε = 0, ñg 1 = ñ1 and ).

Then by using the implicit function theorem ñg 1 can be expanded in ε.


PhD Magali Brunet
email: mbrunet@laas.fr

; Hicham

Mafhoz Kotb

Laurent Bouscayrol

Emmanuel Scheid

Michel Andrieux

Corinne Legros

Sylvie Schamm-Chardon

Patrick Desjardins

H Mafhoz Kotb

Nanocrystallized tetragonal metastable ZrO2 thin films deposited by metalorganic chemical vapor deposition for 3D capacitors

Keywords: Thin Film Devices, Sensors, and Actuators MOCVD, ZrO2, cubic/tetragonal, high-k dielectric, 3D capacitors, FTIR, TEM MOCVD, ZrO 2, cubic/tetragonal, high-k dielectric, 3D capacitors, FTIR, TEM

ZrO2 is a potential candidate for the realization of 3D capacitors on silicon for future Systems-on-Chips. The paper reports the deposition of ZrO2 thin films by metal-organic chemical vapour deposition on planar and 3D structures. Physico-chemical as well as electrical properties of the films are investigated. It is shown that the change of phase and microstructure of the film due to annealing at 900°C under O2 impacts directly on the electrical performance of the capacitors. Capacitance densities are 2 nF/mm² for planar capacitors and reach 8 nF/mm² for capacitors with pores etched in silicon with 4:1 aspect ratio.

Introduction

Today, portable electronic equipment is required to be lighter, more compact, multifunctional, rapid and highly reliable. High performance Systems-on-Chip (SoC) appear to be a logical answer to these challenges. Integrated passive devices (IPD) are particularly interesting for industry as far as the numerous passive components (capacitors, inductors and resistors) required in electronic systems are concerned. Integrated capacitors with a metal-insulator-metal (MIM) structure showing high capacitance per surface area are expected to find many applications such as decoupling, filtering output signals in DC-DC converters, energy storage or A/D conversion [START_REF] Ulrich | Integrated passive component technology[END_REF]. To reach high capacitance densities, two complementary ways must be explored: 3D architectures and integration of high-k dielectrics. In 3D architectures, very high capacitance densities can be reached through the increase of the surface to volume ratio of the electrodes. The technique of high-aspect ratio pores network etched in silicon by deep reactive ion etching (DRIE) has previously been developed [START_REF] Roozeboom | International Symposium on Microelectronics[END_REF]. An alternative technique of n-type silicon electrochemical etching with backside illumination in HF-based electrolytes has also been reported [3]. A capacitance density of up to 650 nF/mm² could be reached with the latter technique and conventional dielectrics (SiO 2 and Si 3 N 4 ) [4]. Introducing high-k dielectrics such as HfO 2 or ZrO 2 or very high-k dielectrics such as (Pb,Zr)TiO 3 or (Ba,Sr)TiO 3 perovskites instead of SiO 2 in these high-aspect ratio structures should allow for capacitance densities higher than 1 µF/mm². For instance, atomic layer deposition (ALD) has been used to deposit Al 2 O 3 dielectric layers into 30 µm deep pores in silicon [5] leading to capacitance densities of 440 nF/mm² or inside anodic aluminium oxide nanopores leading to capacitance densities of 1 µF/mm² [6]. However, the ALD rate is quite low: 0.1-1 Å/cycle with a cycle time ranging from 1 to 10 sec. This limits its use for very thin dielectric layers. For applications such as output filters of DC-DC converters, relatively high breakdown voltages are required (10-15 V) driving the need for thicker layers than for metal-oxide-semiconductor (MOS) type applications. Metal-organic chemical vapor deposition (MOCVD) allows deposition of quite thick layers of binary and ternary oxides with a high deposition rate (0.1-1000 Å/s).

As a high-k dielectric, ZrO 2 is a promising candidate for high density capacitors as it combines a high dielectric constant and relatively high breakdown voltages: 3-6 MV/cm [START_REF] Ulrich | Integrated passive component technology[END_REF]7,8,9]. The dielectric static response of ZrO 2 was calculated by Zhao et al. [10,[START_REF] Zhao | Novel Materials and Processes for Advanced CMOS[END_REF]. It is much higher for the cubic and tetragonal phase ((ε r(c) = 37 and ε r(t) = 38, respectively) than for the monoclinic phase (ε r(m) = 20). Therefore, tetragonal/cubic phases must be present if high capacitance densities are to be reached. We have already deposited ZrO 2 thin films with a tetragonal phase by direct liquid injection MOCVD (DLI-MOCVD) on planar and tridimensional surfaces (high-aspect ratio pores etched in silicon) [START_REF] Galicka-Fau | [END_REF][START_REF] Galicka-Fau | Proceedings of the International Symposium CVD XVII and EUROCVD 17[END_REF]. The influence of the deposition parameters on the physico-chemical properties of the films was investigated. The purpose of this paper is to address, for one set of deposition parameters, the effect of annealing on the morphology, the crystalline structure and the corresponding capacitance of the films for 2D and 3D MIM structures.

Fabrication and experimental methods

Planar and 3D capacitors were fabricated as follows: on n-type two inch silicon wafers, arrays of pores were etched by DRIE with the Bosch process using a photoresist as an etch mask [START_REF] Brunet | [END_REF]. The pores are 4 µm wide with an aspect ratio of 4:1. The wafers were dipped for 1 min in HF to remove the native oxide then in a bath containing tetramethyl ammonium hydroxide (22 wt%) and isopropyl alcohol (17%) heated up at 75°C to smooth the scalloping in the pores produced with the Bosch process. The silicon substrate was highly doped with a phosphorus diffusion step (POCl 3 in N 2 at 1050°C for 35 min) in order to obtain a MIM structure (without depletion in the semiconductor) and to reduce the series resistance of the capacitor. After a RCA cleaning of the substrates, a thin silicon dioxide layer (12 nm) was then grown thermally through dry oxidation at 900°C for 35 min in total (5 min under O 2 /N 2 ). This oxide layer was introduced in the stack for electrical characterization purposes. Indeed, it is known that leakage current through as-deposited high-k films is mainly due to as-grown defect states at the interface with silicon [15,16].

The inclusion of a thin SiO 2 layer, grown on purpose, is a solution to reduce these defect states and act as a barrier for electronic transport [17].

ZrO 2 thin films with a thickness between 50 and 60 nm were deposited by DLI-MOCVD using the following precursor: Zr 2 (O i Pr) 6 (thd) 2 dissolved in cyclohexane and evaporated at 240°C. Several deposition parameters such as substrate temperature (400-900°C), injection frequency (0.5-2 Hz), total pressure in the reactive chamber (10-1000 Pa) and evaporation temperature of the precursor (200-230°C) were investigated in a previous study in order to determine the optimal conditions for the growth of a tetragonal zirconia film [START_REF] Galicka-Fau | [END_REF]. All these parameters were investigated and their effects analyzed. It was found that peculiar deposition conditions provided a uniform film with no cracks and a good crystallization in the expected tetragonal structure [START_REF] Galicka-Fau | [END_REF][START_REF] Galicka-Fau | Proceedings of the International Symposium CVD XVII and EUROCVD 17[END_REF]. The substrate temperature has especially a strong influence on the crystal structure: a substrate temperature of 525°C was selected as it is a good compromise between amorphous deposition (below 450°C) and very poor conformal deposition on 3D structures (above 550°C). Compared to the samples prepared in reference [START_REF] Galicka-Fau | [END_REF], the injection frequency was lowered from 2 Hz to 1 Hz in order to favor the (111) t phase which was proven in the previous study to be more stable during further annealing.

Finally, the oxygen and nitrogen gas flow were set at 0.05 L/min. Several batches were deposited in such optimized conditions; each of them containing up to ten substrates (2" wafer and several 10*10 mm 2 samples) to check them using various characterization techniques. A post-annealing step was carried out on some of the substrates under oxygen gas flow at reduced pressure (100 Pa) to check its chemical stability.

Capacitors were finalized with a boron-doped polysilicon layer as the top electrode etched in a square shape (100 µm x 100 µm and 250 µm x 250 µm). The final step was the deposition of a metallic contact (Ti/Au) on the polysilicon layer as well as on top of the films for planar capacitors. As for the bottom electrode of the capacitors, an AuSb layer was deposited on the back side of the silicon substrate.

Concerning the material characterization, the crystal structure of ZrO 2 in the thin films was checked by grazing incidence (2°) X-ray diffraction (XRD) performed on a PANALYTICAL X'pert pro MRD diffractometer using Cu K radiation. All the patterns were realized with a detector angle (2) varying from 27 to 37°, a step of 0.05° and 144 s counting rate. Nanocrystallized tetragonal zirconia ZrO 2(t) , cubic zirconia ZrO 2(c) and monoclinic zirconia ZrO 2(m) were identified using JCPDS files 17-923, 27-997 and 37-1484 respectively. Fourier transform infra-red (FT-IR) spectroscopy was performed on samples using a Perkin Elmer 2000 FT-IR spectrophotometer in the range 4000-220 cm -1 as it is a sensitive technique to identify the polymorphs of zirconia, even if films are nanocrystallized. Each spectrum was obtained by averaging nine or more scans in order to obtain a good S/N ratio. Selected area electron diffraction pattern (SAED) was also performed on the planar films.

The thin film coverage in 3D structures was controlled on cross-sections with a Hitachi S-4800 scanning electron microscope (SEM). The local structural and chemical properties around the film/substrate interface and in the film were analyzed by transmission electron microscope (TEM)/high resolution (HRTEM) and electron energy loss spectroscopy with a field emission TEM (FEI Tecnai™ F20) operated at 200 kV and equipped with a spherical aberration corrector. TEM samples were prepared by classical procedures cutting the samples either perpendicular to the film/substrate interface (cross-section) or parallel to the film/substrate interface (plan-view). Electrical characterizations were conducted on the fabricated capacitors by performing impedance versus frequency measurements under point probes with an impedance analyser (Agilent 4294A) and current versus voltage measurements on a Keithley 4200 Semiconductor Characterization System.

Results.

Thin film characterization

Figure 1a shows the XRD results for a ZrO 2 /SiO 2 /Si planar stack before and after annealing under O 2 at 900°C. The as-deposited film crystallized in the tetragonal/cubic phase, and after annealing the (111) m line of the monoclinic phase appeared. It is difficult to clearly distinguish the difference between the ZrO 2(c) and ZrO 2(t) phases due to the similarity of their patterns. It was shown by Tsoutsou et al. [18] in their study on Ge-doped zirconia films grown by atomic oxygen beam deposition that it was possible to identify the tetragonal crystalline state through the splitting of both lines at 50° or 60°. In our patterns, the large broadening of the diffraction lines due to a very small crystallite size (lower than 10 nm) hinders the discrimination of both (c) and (t) phases.

Therefore, in order to identify without any doubt the presence of (c) or (t) phases in the films, FT-IR experiments were conducted. The FT-IR spectra presented in figure 1b were compared with reported stabilized cubic, tetragonal and monoclinic zirconia spectra [19,20]. Group theory predicts only one IR active mode for cubic zirconia that signs in figure 1b as a weak broad band between 600 and 350 cm -1 and an extra band at 589 cm -1 due to a secondary oscillator of the cubic phase. The three active modes reported for tetragonal ZrO 2 correspond to the weak bands that emerged from the broad band between 350 and 600 cm -1 . The monoclinic phase can be recognized when sharpen bands are added to the spectra. They correspond to the fifteen vibrational modes of the monoclinic phase of zirconia. The position of these bands is red shifted compared to those obtained theoretically on bulk zirconia. Moreover, they are broadened and flattened. These differences in the broadening and position of the absorption band could be attributed to various parameters such as the reduction of the particle size (nanostructuration), the shape of particles and a possible disordering of the oxygen sublattice [21][22][23]. However, despite the differences in the position of the bands, the presence of well-defined modes in the spectra supports the identification of the tetragonal and monoclinic phases rather than the cubic phase either in the as-deposited and annealed films. In the annealed films, the monoclinic character is emphasized.

From TEM analyses carried out on the ZrO 2 /SiO 2 /Si planar stack before and after annealing at 900°C under O 2 , the layer thickness was measured to be 64 and 53 nm on the as-deposited and the annealed samples, respectively (figure 2 a and b). It is worthy to note that the thickness difference between the as-deposited and the annealed sample is probably due to the annealing process that is able to densify and/or to eliminate the carboxylates impurities adsorbed during the CVD process. As a consequence, it is possible to observe a modification of the microstructure of the film after the heat treatment. Indeed, as shown in figure 2, the as-deposited ZrO 2 layer is polycrystalline, with crystallite sizes between 5 and 10 nm. The SAED patterns of the as-deposited ZrO 2 layer provide evidence of mainly the cubic and or tetragonal ZrO 2 structure (figure 2c, insert). As for XRD patterns, the two structures are difficult to separate. However, the ring corresponding to the 0.213 nm distance can only be attributed to the low level diffracting (102) plane of the ZrO 2(t) phase supporting the FT-IR conclusions of a major contribution of the tetragonal phase with respect to the cubic one. After annealing, a columnar organization is revealed with an increase in crystallite size due to the grain growth mechanism. The monoclinic structure appears there in addition to the cubic/tetragonal one (figure 2d, insert). In addition, the thermal SiO 2 layer, which was 12 nm, did not change after annealing (figure 2a et b).

In tridimensional structures, the coverage of the deposit was observed by SEM (figure 3).

The deposition is not fully conformal: there is a loss of about 54% between the top thickness and the bottom thickness. However, it was verified with TEM analysis (figure 4a to d), thanks to a cross-section inside a pore, perpendicular to its long axis, that the ZrO 2 layer is deposited homogeneously on the vertical pore walls. The same conclusions for the stack in planar configuration can be drawn concerning the structure of the deposit: the asdeposited ZrO 2 layer is constituted of very small crystallites of cubic/tetragonal phase with a major contribution of the tetragonal phase (SAED not shown) that grow after annealing and transform to the monoclinic phase (SAED not shown) (figure 4-e and f).

Electrical characterization

Figure 5 shows the capacitance density versus frequency of planar and 3D capacitors. The capacitance was 8 nF/mm² for 3D capacitors with as-deposited ZrO 2 , which is about three times better than planar capacitors. In addition, annealed samples always have lower capacitance than as deposited ones (5.4 nF/mm² for 3D capacitors, 1.9 nF/mm² for planar capacitors).

The measured capacitance is an equivalent one, which contains the capacitance of the SiO 2 /ZrO 2 stack with parasitic incorporated due to electrodes and interfaces. Models were thus developed on planar capacitors to extract the pure capacitance of the ZrO 2 layer and therefore, the permittivity. The model curves are superposed on the measured curves in figure 5 for planar capacitors with as-deposited ZrO 2 and annealed ZrO 2 . The equivalent circuits used for modeling are shown in figure 6. In these models, identified contributions are: C ox (the capacitance of the SiO 2 layer), Rp ox and Rp ZrO2 (the parallel resistance of each layer representative of the leakage current), and Rs (the series resistance due to electrodes).

In the model for the as-deposited film, the branches C1 R1, C2 R2, C3 R3 and C4 R4 are representative of interface states distributed in energy [24]. Only one C1 R1 interface branch was necessary for modeling the annealed film. With these models, the capacitance of the ZrO 2 layer (C ZrO2 ) was extracted as 190 pF for both the as-deposited film and the annealed film. Using a thickness of 64 nm and 53 nm for the as-deposited and annealed ZrO 2 films, respectively, the permittivities of theZrO 2 layers are calculated to be 27 and 22.7.

In figure 7, for as-deposited and annealed samples (planar capacitors), the current density is plotted versus the electrical field of the stack capacitor where E = V/t stack ; t stack = t SiO2 + (ε SiO2 /ε ZrO2 )*t ZrO2 [25] In both case, the current at low bias (below 1-2 MV/cm) followed the same behavior and was of the order of 1 nA/cm². At higher bias, i.e. above 2.5 MV/cm for the as-deposited layer and above 1.5 MV/cm for the annealed layer, the current mechanism changed: a sharp increase of the leakage current was observed. The breakdown took place at a lower bias (around 2 MV/cm) for the annealed sample than for the as-deposited sample.

Discussion

4.1.

On the material structure.

It is well known that preparation methods (physical or chemical processes), preparation parameters (temperature of growth and O 2 partial pressure) [7,26,27] and film thicknesses [27][28][29][30] greatly influence the crystalline structure of the deposited material. In particular, for the zirconia system, general trends can be drawn for ZrO 2 polymorph. The monoclinic phase is favoured by high growth temperature and thick films. Furthermore, the crystallographic structure is dependent on the size of the crystallites [31]. This point is well discussed in the works of Bernard, Filipovitch and Garvie [27,32,33]. The general principle discussed is that the tetragonal phase is stabilized in thin films because of the small size of the crystallites having a lower surface energy compared to the crystallites of the monoclinic phase.

It was shown in figure 1b that the FT-IR spectra of as-deposited and annealed films are clearly structured with an increase in the number of modes going from as-deposited to annealed films. This supports the results observed by XRD (figure 1a) showing the growing of the monoclinic phase with annealing. In the as-deposited films, the monoclinic phase is also present, though with a minor contribution, together with the cubic phase (secondary oscillator band at 589 cm -1 ). However, two observations, clearly visible in the as-deposited films, support the presence of the tetragonal phase in the films. First, if we consider that the cubic or tetragonal phase of ZrO 2 is mainly present in the as-deposited films as revealed in the XRD patterns, the presence of structured modes in the corresponding FT-IR spectra can only support the presence of mainly the tetragonal phase.

Second, we showed in a previous work [START_REF] Galicka-Fau | Proceedings of the International Symposium CVD XVII and EUROCVD 17[END_REF] that the cubic zirconia phase transforms more easily into the monoclinic phase than the tetragonal phase during annealing under oxygen flow. As a consequence, in the present study, if the cubic phase amount was very large in the as-deposited films, they should easily transform during annealing into the monoclinic phase and a large number of monoclinic absorption bands should be observed on the spectrum. As shown in figure 1b, the monoclinic absorption bands do not evolve strongly during annealing. This confirms that the as-deposited films crystallize mainly in the tetragonal structure.

In this work, for a film thickness of around 60 nm, it has been verified with three different methods (XRD, FT-IR and SAED) that the deposited ZrO 2 crystallizes mainly in its tetragonal structure with crystallite sizes of the order of 5 nm. In addition, the monoclinic phase appears in the corresponding annealed film with crystallite sizes of around 15 nm.

This work confirms our previous results that for thicknesses in the range of 50-100 nm, it is possible through particular deposition conditions to stabilize the tetragonal zirconia phase [START_REF] Galicka-Fau | [END_REF]. As a consequence, initial crystallite size associated with grain growth during deposition and/or annealing may be the key parameter to control the existence of the tetragonal metastable phase.

4.2.

On the electrical performance.

The measured permittivities are 27 and 22.7 for the as-deposited layer and the annealed layer, respectively. This difference is attributed to the crystalline phase of the ZrO 2 layer: the contribution of the monoclinic phase appearing after annealing presents a lower permittivity, which decreases the overall permittivity. The measured ε r for (c)/(t) phase is lower than the theoretical value (ε r = 38) calculated by Zhao et. al [10,[START_REF] Zhao | Novel Materials and Processes for Advanced CMOS[END_REF], but higher than permittivities usually reported in the literature: between 16 and 25 [7,9,29,34,35]. It was also observed that the electrical model for the annealed sample is much simpler than for the as-deposited one. This verifies the common knowledge that annealing reduces the defect states at the interfaces [15,16,17]. For 3D capacitors, the same decrease in capacitance was observed: it is 8 nF/mm² for as-deposited films and 5.4 nF/mm² after annealing. The gain in capacitance density from planar to 3D capacitors is about 3 for asdeposited layer or annealed layer. A simple geometry calculation based on relation [START_REF] Ulrich | Integrated passive component technology[END_REF] indicates that the enhancement factor (3D electrode surface (S 3D )/ Planar electrode surface (S plan )) should be 4.7 for the realized pores. with t stack the dielectric thickness, n the number of pores etched in the substrate, w the pore width and h the pore depth.

The difference between the measured enhancement and the calculated one is due to the ZrO 2 layer which is not conformal on the 3D structure: the layer at the bottom is thinner (25 nm) than at the top (55 nm) reducing the equivalent permittivity of the SiO 2 /ZrO 2 stack. The lack of conformity in high-aspect ratio structures can be explained by the fact that the deposition is diffusion limited because of a high surface reaction rate. Although low temperature and low pressure should reduce surface reaction and promote diffusion of reactive species in 3D structures against a surface reaction, the MOCVD equipment was not designed to operate below 100 Pa. Furthermore, lower temperature would lead to a non stoechiometric zirconium oxide, which would enhance leakage current.

The measured leakage currents are representative of the silicon oxide layer because in dielectrics stack, the oxide of best quality, SiO 2 in this case, drives the current mechanism.

It can be noted first that the leakage current density is much lower than for ZrO 2 layers deposited by MOCVD in the same conditions directly on silicon [36]. For bias above 1 MV/cm, the current density is increasing drastically: this behavior was already reported in the literature [37]. It is believed that in this configuration (bi-layer SiO 2 /ZrO 2 ), electrons first tunnel through the SiO 2 layer (direct tunneling) and then tunnel through traps located in the high-k layer [START_REF] Autran | Fundamental and Technological Aspects of High-k Gate Dielectrics[END_REF], which is representative of Poole-Frenkel emission.

Similarly to the leakage current, the breakdown voltage plotted in figure 7, is representative of the SiO 2 layer. The breakdown voltage is about 6 MV/cm for the asdeposited layer and 2 MV/cm after annealing. These values are much lower than what is the usually reported breakdown voltage for thermal SiO 2, i.e. around 13 MV/cm. The phosphorus diffusion step before the thermal oxidation is probably one cause of a lower breakdown voltage as doping impurities may diffuse in the oxide during its growth. On the other hand, the annealing step under O 2 is another cause of the degraded breakdown voltage, since the microstructure of the film is changing. Indeed, after annealing the grain size is increasing and a preferred orientation is favored in the layer: the (200) t peak is slightly increased in the XRD pattern (fig. 1a) and in the TEM analysis is associated with a columnar microstructure that is clearly visible (fig. 2b). This columnar structure leads to more conduction paths through grain boundaries, eventually causing premature breakdown of the stack.

Conclusions

Planar and 3D MIM-like capacitors with ZrO 2 as the high-k dielectric were prepared and characterized. The tetragonal phase of ZrO 2 , in films with a thickness ranging from 50 to 60 nm, was stabilized through peculiar MOCVD conditions. It was demonstrated that the change of phase and microstructure of the film due to annealing at 900°C under O 2 impacts directly on the electrical performance of the capacitors. Intrinsically, the permittivity of ZrO 2 is lowered with the annealing step due to the partial tetragonal/monoclinic phase transformation. Extrinsically, although the quality of the interfaces of the overall stack is improved, the change in microstructure (larger grains and columnar structure) is responsible for the higher leakage current and lower breakdown voltage.

3D capacitors showed four times larger capacitance density (8 nF/mm²) compared to planar ones (2 nF/mm²). Although the deposition parameters giving conformal deposition were not found for this precursor in this particular MOCVD equipment, it is believed that conformal deposition of ZrO 2 in very high-aspect ratio structures should lead to very high capacitance densities. Figure 3. SEM images of ZrO 2 /SiO 2 /Si coverage inside pores (4.9 µm*4.9 µm*11.9 µm).

Thickness of 55 nm and 25 nm at the top and at the bottom of the pore respectively.

Figure captures:Figure 1 .

 1 Figure captures:

Figure 2 .

 2 Figure 2. TEM images of ZrO 2 /SiO 2 /Si stack a) and c) as-deposited, b) and d) after an

Figure 4 .

 4 Figure 4. TEM images of ZrO 2 /SiO 2 /Si inside pores (4.9 µm*4.9 µm*11.9 µm) in cross-

Figure 5 .

 5 Figure 5. Measured specific capacitance versus frequency for planar and 3D capacitors, as-

Figure 6 .

 6 Figure 6. Equivalent circuits for planar capacitors with SiO 2 /ZrO 2 layer a) as-deposited, b)

Figure 7 .

 7 Figure 7. Characteristics J(E) of as-deposited and annealed planar capacitors.

  Figure 1 Click here to download Figures (if any): fig 1a et 1b ter.doc

  Figure 3 Click here to download Figures (if any): Figure3.pptx

  Figure 6 Click here to download Figures (if any): Figure_6a-6b.doc

Figure 7

 7 Figure 7 Click here to download Figures (if any): Figure7.docx

Acknowledgment

This work was funded by French National Research Agency (project ANR-06-JCJC-0081).

The authors are also grateful to I Gallet, M. Herbst-Ghysel and K. Galicka-Fau for their participation in the experimental process.