Stéphane Brull
email: stephane.brull@math.u-bordeaux1.fr.

Vincent Pavan
email: vincent.pavan@polytech.univ-mrs.fr.

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

Derivation of BGK models for mixtures

Keywords: Kinetic theory, gas mixtures, BGK models, entropy minimization, hydrodynamic limit

) on the introduction of relaxation coefficients and a principle of minimization of the entropy under constraints of moments. These free parameters are compared with the free parameters introduced in the Thermodynamics of Irreversible Processes approach of the Navier-Stokes system. At the end the BGK model is proved to satisfy Fick and Newton laws.

1 Introduction.

The BGK equation [START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF] is a seminal simplified model of the nonlinear Boltzmann equation of gas dynamics. Its features are to replace the complicated integral of collisions with a relaxation model while keeping some important physical and mathematical properties of the interaction term (conservation laws, H-theorem, equilibrium states,...). Thus it is very interesting from a numerical point of view. Its drawbacks rely upon a somewhat bad description of the distribution function far from equilibrium and closed to equilibrium in the sense that the Navier-Stokes limit might be not correct (Prandt number). Though a correction can be brought to obtain this important ratio between viscosity and heat conductivity (see the Ellipsodal Statistical Model by Holway [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF]). In the case of gas mixtures things are less clear. The Boltzmann equation(s) is still considered to describe the true behavior of the distribution function but BGK models and their extensions meet fundamental difficulties. The first one is that conservations laws stand only for the whole set of particles; exchanges of momentum and energy occur between each species. The second is that the hydrodynamical limit is much more complicated. Phenomena such as diffusion or thermal diffusion must be considered. Many of those problems have been widely studied by Japanese researchers (see for example [START_REF] Sone | Doi Kinetic theory analysis of gas flows condensing on a plane condensed phase: Case of a mixture of a vapor and noncondensable gas[END_REF][START_REF] Aoki | The behaviour of a vapor-gas mixture in the continuum limit: Asymptotic analysis based on the Boltzman equation[END_REF][START_REF] Aoki | Vapor flows caused by evaporation and condensation on two parallel plane surfaces: Effect of the presence of a noncondensable gas[END_REF][START_REF] Aoki | Vapor flows with evaporation and condensation in the continuum limit: effect of a trace of non condensable gas[END_REF] and show great difference with the usual monatomic case. Finally theoretical aspects such as existence theorems [7,[START_REF] Brull | The Boltzmann equation for a two component gas in the slab for soft forces[END_REF][START_REF] Brull | The stationary Boltzmann equation for a two component gas in the slab for different masses[END_REF] or study of a binary mixture closed to a local equilibrium [START_REF] Brull | Problem of evaporation-condensation for a two component gas in the slab[END_REF] confirm the specificity of multi-component case. Coming back to modeling there exists however a great variety of BGK models which traces back to the work of Gross and Krook [START_REF] Gross-M | Krook Model for collision processes in gases: smallamplitude oscillations of charged two-component systems[END_REF] to the most recent model by Kosuge [START_REF] Kosuge | Model Boltzmann equation for gas mixtures: Construction and numerical comparison[END_REF]. A first family of model has been proposed [START_REF] Gross-M | Krook Model for collision processes in gases: smallamplitude oscillations of charged two-component systems[END_REF][START_REF] Sirovich | Kinetic modeling of gas mixtures[END_REF][START_REF] Hamel | Kinetic model for binary gas mixtures[END_REF][START_REF]Morse Kinetic model equations for a gas mixture[END_REF] (often in the context of Maxwellian molecules). Those BGK models intend to mimic the Boltzmann equations for self and cross interaction and some of them give the right exchange of momentum and energy between each species [START_REF] Hamel | Kinetic model for binary gas mixtures[END_REF][START_REF]Morse Kinetic model equations for a gas mixture[END_REF]. A second family adds to the previous properties the so-called indifferentiability principle that reduces the whole family of relaxation models to a single one when all masses and cross sections are equal [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF][START_REF] Aoki | Perthame A consistent BGK-type model for gas mixtures[END_REF]. Finally Kosuge has generalized the model of Garzo & al [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF] to relaxation terms that allow to approximate exchanges of momentum, energy, cross moments of order two and heat flux between each species. Though numerical results may be quite good for some models or may enjoy nice mathematical properties for others it is quite surprising that none of them has attempted to reach the right hydrodynamical limit. In general their authors rather study the hydrodynamical limit of their model and eventually compared it with the right one [START_REF] Hamel | Two-fluid hydrodynamic equations for a neutral, disparatemass, binary mixture[END_REF]. Our approach goes the other way. We first consider the right hydrodynamical limit and then construct a model that allows to recover the kinetic coefficients -in the present case the Fick law and the viscosity.

The paper is constructed as follows. In section 2 we recall the main properties of the Boltzmann equation and the linearized operator which lead us to propose a list of necessary properties that should satisfy any approximation of the interaction terms. We also introduce a space of moments denoted by C (see [START_REF] Aoki | Perthame A consistent BGK-type model for gas mixtures[END_REF] that features the non-conservation of momentum per species. Section 3 is the cornerstone of our paper. After recalling the Navier-Stokes equation for gas mixtures and how the thermodynamic of irreversible processes allow to define the fluxes and the different laws (Fick, Soret, Dufour) that do not exist in the case of a single specie, we revisit the way those fluxes can be computed starting from the kinetic theory. Though this connection has already been done in the celebrated book by De Groot and Mazur [START_REF] De Groot | Nonequilibrium Thermodynamics[END_REF] and in a certain sense by Chapman and Cowling [START_REF] Chapman | The mathematical theory of non uniform gases[END_REF], we present a synthetic and clear computation of the kinetic coefficients that is based on one hand on the mathematical properties of the linearized Boltzmann operator and on the other hand on the space C . Section 4 is devoted to the construction of our model. It leans on two ideas that were introduced by two of the authors [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF][START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF]: the operator must relax certain moments of the distribution functions so as to obtain the correct hydrodynamical limit and and the model is constructed via a principle of minimization of entropy under moment constraints. In section 5 we prove that the whole set of properties that were introduced in section 2 are satisfied and we easely derive from them and from the formalism of section 3 the hydrodynamical limit. We also address the properties of indifferentiability and correct hydrodynamical limit for BGK models. Finally section 6 is devoted to a discussion implying physics, modeling and mathematical aspects.

2 Notations and recall about the Boltzmann and other general kinetic equation for gas mixtures.

2.1 The Boltzmann operator for mixtures.

Detailed expression.

Let us consider a gas mixtures with p components. The functions

f i (t, x, v) (or for short f i , i ∈ [1, p] with f := (f 1 , • • • , f p )
) associated to a given species i evolves according to the Boltzmann equation:

∀i ∈ [1, p] , ∂ t f i + v • ∇ x f i = k=p k=1 Q ki (f k , f i ) := Q i (f, f) , (1)
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 wich depend on the interaction potential of the 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) • ω ω ω) ω ω ω.
Those equations satisfy the conservation of momentum and energy at a microscopic level

m i v + m k w = m i v * ki + m k w * ki , m i v 2 + m k w 2 = m i v * ki 2 + m k v * ki 2 .

Macroscopic quantities for the mixture.

We denote with n i , ρ i , u i , E i , T i and E i the macroscopic quantities representing respectively the number density, density, average velocity, energy per unit volume, energy per particle and finally temperature of a given specie 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 v f i dv, E i = 1 2 ρ i u i 2 + n i E i , E i = 3 2 k B T i = m i 2n i R 3 v -u i 2 f i dv,
where k B is the Boltzmann constant. In the same way, we can introduce hydrodynamic quantities for the mixture

n = p k=1 n k , ρ = p k=1 ρ k , ρu = p k=1 ρ k u k , ( 2
)
nE + ρ 2 u 2 = E = p k=1 E k , E = 3 2 k B T. (3)
Given a mixture of p species featuring the parameters n i , u, T , an important list of functions are the Maxwellians of equilibrium reading as:

∀i ∈ [1, p] , M i = n i (2πk B T /m i ) 3 2 exp - m i (v -u) 2 2k B T . (4)
In the sequel we will denote by M := (M

H (f) := i=p i=1 R 3 (f i ln (f i ) -f i ) dv.

Spaces of interest.

Using the above notations we note as L 2 (M) the set of measurable functions Ψ = (ψ 1 , • • • , ψ p ) such that:

Ψ 2 := i=p i=1 R 3 ψ 2 i M i < +∞.
This space is equipped using its natural dot product:

Ψ, Φ = i=p i=1 R 3 ψ i φ i M i dv.
In the sequel, we will often use the dot product notation even if what lies in the bracket is not a list of p scalar functions but is a list of p tensorial functions. For instance, if Ψ, Φ are two lists of p (symmetrical) tensorial functions ψ i , φ i , i ∈ [1, p], then the dot product notation ψ, φ should be understood as

ψ, φ := i=p i=1 R 3 ψ i ⊗ φ i M i dv,
where ⊗ denotes the usual tensorial product. Note also the following convention: assume that V is a list of p vectors and T a list of p square matrixes. Then a notation like αV, αT (where α is a scalar) means that the scalar α is distributed on each line of the list of vectors or matrixes. Besides, if β β β is a vector, then a notation like β β β • V means that we distribute the dot product by β β β on each line of the vector V. The same way, if γ γ γ is a square matrix, then a notation like γ γ γ : T means that we distribute the (total) dot product (between matrix) by γ γ γ on each line of the line of T. Finally, if s is a tensor and S a list of p tensors, then a notation like s ⊗ S means that we distribute the tensor product by s on the left on each line of the list S.

Introduce now the following subspace K of L 2 (M ) of dimension p + 4 generated by the following list of functions:

     1 0 . . . 0      , • • • ,      0 0 . . . 1      ,      m 1 v x m 2 v x . . . m p v x      ,      m 1 v y m 2 v y . . . m p v y      ,      m 1 v z m 2 v z . . . m p v z      ,      m 1 v 2 m 2 v 2 . . . m p v 2      .
This space is of dimension p + 4 and the above list of functions is noted φ φ φ l , l ∈ [1, p + 4]. K is the natural space of collisional invariants but contrarily to the case of monatomic gas there exists a "complementary" space C of moments of degree 1 in velocity which is not conserved. Our purpose is now to exhibit a basis of C. Definition 1. Let C i be the vector having 0 everywhere except on its i th line where it lies vu. Note by P K the orthogonal projection on K and I the identity operator. Then we define C as the space generated by the vectors

(I -P K ) (C i ) , i ∈ [1, p].
We have the following lemma Lemma 1. The family (I -P K ) (C i ) , i ∈ [1, p] is composed of p-1 independant "vectors" and as a consequence the dimension of C is 3 (p -1).

We postpone the proof of this lemma in appendix.

2.2 Properties of the Boltzmann operator.

2.2.1 Collisional invariants, H-theorem and equilibrium states.

1. Collisional invariants: for all f = (f 1 , • • • , f p ) the following equation holds i=p i=1 j=p j=1 R 3 Q ji (f j , f i ) φ i dv = 0 iff φ φ φ ∈ K.

H-theorem: for any list of nonnegative functions

f = (f 1 , • • • , f p ) it holds that i=p i=1 j=p j=1 R 3 Q j,i (f j , f i ) ln f i dv ≤ 0. (5)
3. Equilibrium states: the equality holds in the above equation if and only if f is at thermodynamical equilibrium i.e there exists hydrodynamic values

n 1 , • • • , n p , u, T such that ∀i ∈ [1, p] , f i = M i
In such a case we denote f = M. Moreover the equality in (5) holds if and only if

i=p i=1 j=p j=1 R 3 Q j,i (f j , f i )dv = 0 that is iff f i = M i .

Linearized operator.

Define now the linearized Boltzmann operator

L B := (L B,1 , • • • , L B,p ) operating on g = (g 1 , • • • , g p ) ∈ L 2 (M) as ∀i ∈ [1, p] , L B,i (g) := lim ǫ→0 1 ǫM i j=p j=1 Q ji (M j (1 + ǫg j ) , M i (1 + ǫg i )) .
Basically, such an operator may be computed as:

L B,i (g) = 1 M i   p j=1 Q ji (M j , M i g i ) + Q ji (M j g j , M i )   . (6)
Property 1. The operator L B is self-adjoint and non-negative and its kernel is exactly K.

Another important property of the linearized operator is the following: Property 2. L B is the sum of a diagonal operator f → νf , such that (νf ) i = ν i (v) f i and a compact operator K (for regular enough differential cross sections) such that for any i ∈ {1, ..., p},

(Kf ) i = p j=1 R 3 ×S 2 M j f i (v * ji )σ ji (V • ω ω ω, V ) V dω ω ωdw + 1 M j Q ji (M j f j , M i ) .
As a consequence L B is a Fredholm operator and L B : K ⊥ → K ⊥ is continuous and invertible. For the sequel, we note as L -1 B the pseudo inverse of L B . It is an obvious self adjoint positive operator on K ⊥ .

Proof. This proof is given in details in [START_REF] Aoki | layer for a gas mixture[END_REF] in the context of a two component gas and can be easily generalized to a given number of components.

2.3

The Boltzmann operator as a model for other kinetic operators: minimum properties for a BGK operator.

In many practical situations the computation of the full Boltzmann operator is expansive. This is the reason why one looks for BGK ones which will denoted as R := (R 1 , • • • , R p ) in the sequel. For mono-specie system such operators exist already and an abundant literature may be found on the subject. However for multi-species systems things are less clear even if there exist already different models. Our framework is that of the classical Boltzmann operator, that is R must satisfy the following list of properties. The first ones read

R i (f) = 0 ⇔ ∃n i , u, T s.t. ∀i ∈ [1, p] , f i = M i (7)
while the other are those of section 2.2.1

∀f, f i ≥ 0, ∀φ φ φ, i=p i=1 R 3 R i (f) φ i dv = 0 ⇔ φ φ φ ∈ K, (8)
∀f,

f i ≥ 0, i=p i=1 R 3 R i (f) ln (f i ) dv ≤ 0, (9)
R 3 R i (f) ln (f i ) dv = 0 ⇔ ∃n i , u, T s.t. ∀i ∈ [1, p] , f i = M i . (10)
If these properties are satisfied we define the linearized BGK operator

L := (L 1 , • • • , L p ) acting on g ∈ L 2 (M) with ∀i ∈ [1, p] , L i (g) = lim ǫ→0 1 ǫM i R i (M 1 (1 + ǫg 1 ) , • • • , M p (1 + ǫg 1 )) .
Then the second list of properties concerns the linear operator.

Ker (L) = K and L is Fredholm on K ⊥ , (11)
L is self adjoint negative on K ⊥ . ( 12
)
As a consequence of the Fredholm property the operator L is continuous and one to one on K ⊥ as well as its pseudo inverse L -1 : K ⊥ → K ⊥ which is also self adjoint and negative on K ⊥ .

Definition 2. Any kinetic operator satisfying all above properties (7,[START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF][START_REF] Brull | Problem of evaporation-condensation for a two component gas in the slab[END_REF][START_REF] Brull | The stationary Boltzmann equation for a two component gas in the slab for different masses[END_REF][START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]) is said to be properly defined.

3 The Navier-Stokes equations for mixtures and the Boltzmann expansion: importance of the linearized operator.

In this section we aim firstly at introducing the hydrodynamic equations that lead the behaviour of gas mixtures. We recall how the "new" terms appearing in the Navier-Stokes equations are interpreted from a physical point of view using the thermodynamics of irreversibility processes (TIP). Yet the measurements of the kinetic or hydrodynamic coefficients is a hard task and it is of interest to look at the hydrodynamical limit of the Boltzmann equation to approximate them. Following the path of Chapman and Cowling for binary gas mixtures we secondly perform a Chapman-Enskog expansion of the Boltzmann equation in the case of more than two species. Though the previous physically correct Navier-Stokes equations are necessary to make some conclusions concerning the linearized Boltzmann operator. This comparison between the two theories also allows to derive some conditions on the hydrodynamical coefficients. In particular we are going to see that the computation of those coefficients and their properties only depends on the linearized kinetic operator L.

The Navier-Stokes equations.

The Navier-Stokes system for a mixture of p components reads

∀i ∈ [1, p] , ∂ t n i + ∇ • (n i u + J i ) = 0, (13)
∂ t (ρ u) + ∇ • (P + ρ u ⊗ u + J u ) = 0, ( 14
) ∂ t E + ∇ • (Eu + P [u] + J u [u] + J q ) = 0, (15)
where J i , J u J q are respectively the mass, momentum and heat fluxes. Such fluxes have been observed experimentally since a long time. From a phenomenological point of view they have been linked to the gradients of macroscopic variables like density, velocity and temperature through coefficients or tensors. One of the challenge of the transport theory consists in finding models to compute these coefficients. Technically this can be achieved from at least two ways: one by using the general construction of the thermodynamic of irreversible processes (TIP), the other by considering the hydrodynamical limit of the Boltzmann equation. The second method will be detailed in this article. But we first need to give few words about the first one.

The Thermodynamic of Irreversible Processes.

The classical reference concerning the thermodynamics of irreversibility precocesses (TIP) for gas mixtures is the book of De Groot and Mazur [START_REF] De Groot | Nonequilibrium Thermodynamics[END_REF]. In this book it is shown how to define the different fluxes J i , J u and J q basing on the second principle of thermodynamic. Though the structure of those fluxes are expressed in a slightly different form and one can show that they can be written as

J i = j=p j=1 L ij ∇ -µj T + L iq ∇ 1 T J q = j=p j=1 L qj ∇ -µj T + L qq ∇ 1 T J u = L uu D (u) , (16)
where µ i is the chemical potential of the species i in the mixture and D (u) the traceless part of the deformation tensor. In a mixture of ideal gases the chemical potential by species is known (up to some additional constant) through the relation:

- µ i T = k B ln (n i ) - 3 2 ln 2πk B T m i . ( 17
)
Moreover if one assumes that the Casimir-Onsager relations are satisfied for kinetic coefficients then the following matrix:

L :=   L ij L iq 0 L qi L qq 0 0 0 L uu   (18)
must be symmetrical non positive 1 . On top of that it must be added the mass conservation equation stating that

i=p i=1 m i J i = 0 ⇒ ∀j ∈ [1, p] , i=p i=1 m i L ij = 0. ( 19
)
As a consequence the rank of the matrix L ij is at most p -1. However the TIP is not capable anymore in whole generality to say more things about L. In particular if one uses the relation [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF] and the symmetry of the matrix L (18) then whe can only infer that there are p (p + 1) /2 + 1 unknowns to search for.

Using the phenomenological point of view such fluxes are preferentially expressed as gradient of density, temperature and velocity. They read

J i = j=p j=1 D ij ∇n j + D iT ∇T, J q = j=p j=1 D qj ∇n j -D qq ∇T,
1 Following precisely the TIP, such a matrix must be actually symmetrical non negative. This discrepency is simply linked to the definition one uses for the entropy. To perfectly match the thermodynamical approach, we should have defined

H (f) := -k B i=p X i=1 Z R 3 (f i ln (f i ) -f i ) dv
where k B is the Boltzmann constant. The choice we have made for the entropy is however the most usual for mathematical reasons as it enables to consider the minimization of a convex functions rather than maximization of concave functions. As a consequence of this definition, the corresponding Onsager matrix must be symetrical non positive rather than symmetrical non negative where D ij and D iT , D qj and D qq respectively denote the Fick, Soret, Duffour and Fourier coefficients. The main reason for introducing such coefficients lies on the possibility to measure them conveniently from experimental data while kinetic coefficients cannot be obtained directly. There drawback are the lost of symmetry -Casimir-Onsager relations-on the phenomenological coefficients. Remark however that for mixture of ideal gases the one to one correspondance between those coefficients and the matrix L is straightforward. As we will see in the next section, writing the different fluxes in term of the kinetic coefficients ( 16) is the most natural choice not only for the TIP but also to express them starting from the kinetic theory.

Basic recalls on the Chapman-Enskog theory.

The Chapman-Enskog theory begins with performing a change in scale in the time and space coordinates and considers the rescaled kinetic equation:

∀i ∈ [1, p] , ∂ t f ǫ i + v • ∇ x f ǫ i = 1 ǫ k=p k=1 Q ki (f ǫ k , f ǫ i ) . ( 20
)
By letting then ǫ → 0 one searches for an asymptotical expansion of f ǫ i as [START_REF] Hamel | Kinetic model for binary gas mixtures[END_REF] remains bounded. As a consequence, at order ǫ -1 one must choose f 0 i such that the collision operator vanishes, meaning then

f 0 i + ǫf 1 i + • • • such that
∀i ∈ [1, p] , f 0 i = M i , (21)
where the hydrodynamical parameters n i , T, u defining the M i depend on space and time. Then setting f 1 i = M i g i and considering the order ǫ 0 , one must solve the equation:

∀i ∈ [1, p] , (∂ t + v • ∇ x ) [M i ] = M i L B,i (g 1 , • • • , g p ) , (22)
where

L B := (L B,1 , • • • , L B,p
) is the linearized collision operator as it as been defined in the former section. Next we aime to calculate the left-hand side of [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF] in such a way that the expression of g can be easily computed. Following the path of Levermore ([25]), the equation ( 22) may be written as follows:

Lemma 2. The linearized operator defined in ( 22) is given by the expression

L B (g) = j=p j=1 k -1 B (I -P K ) (C j ) • ∇ -µ j T + A : D (u) + B • ∇ 1 T ( 23
)
where I refers to the identity operator on L 2 (M), P K the orthogonal projection on K, while the list of tensors A, B (that belong to K ⊥ ) feature on their i th line the tensors:

(A) i = m i (v -u) ⊗ (v -u) - 1 3 v -u 2 I , (B) i = (v -u) 1 2 m i (v -u) 2 - 5 2 k B T
while the tensor D (u) reads as:

D (u) = 1 2 ∇ x u + (∇ x u) T - 1 3 (∇ x • u) I.
The proof of Lemma 2 is left to appendix. Thanks to Properties 1 and 2, such an equation has one unique solution in g ∈ K ⊥ .

3.4 Link between the Boltzmann equation and the Thermodynamic of Irreversible Processes.

While this theory has already been set in [START_REF] De Groot | Nonequilibrium Thermodynamics[END_REF] and in some way in [START_REF] Chapman | The mathematical theory of non uniform gases[END_REF] we present here a concise formalism based on the mathematical properties of the linearized operator L B and on the space C. Once the perturbative function g has been computed, one can recover the formal expressions of the different fluxes as Lemma 3. For g defined by relation 23 the fluxes J i , J u and J q are given by

J i = g, (I -P K ) (C i ) , J u = g, A , J q = g, B .
Thanks to Lemma 2, this enables to recover the expression of the kinetic coefficients. Indeed L ij , L iq , L qi , L uu and L qq can be respectively obtained from

L ij I = k -1 B L -1 B [(I -P K ) (C i )] , (I -P K ) (C j ) , (24)
L iq I = L qi I = L -1 B (B) , (I -P K ) (C i ) , (25)
L uu I ⊗ I = L -1 B (A) , A , (26)
L qq I = L -1 B (B) , B . (27)
When using the linearized Boltzmann operator, the fact that tensors finally reduce to scalar coefficients (meaning that they are spherical) is linked to the Galilean invariance of the differential cross section in the full Boltzmann collision operator. Here we see immediately that the Onsager-Casimir properties L (18) are derived from the properties of symmetry and non positiveness of the linearized operator L -1 B . It is interesting to note here, that the kinetic coefficients L ij are known as soon as we know the vectors (I -P K ) (C j ) , j ∈ [1, p]. In short this means that the Fick coefficients form a family p-1 degree of freedom. For example in the case of two species they can be computed from a single one. Remark 1. All properly defined operator in the sense of definition 2 gives a Navier-Stokes equation at the hydrodynamical limit with kinetic coefficients given by [START_REF] Kestin | Equilibrium and transport properties of the noble gases and their mixtures at low density[END_REF][START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Kosuge | Model Boltzmann equation for gas mixtures: Construction and numerical comparison[END_REF][START_REF]Morse Kinetic model equations for a gas mixture[END_REF] replacing L B by the proper linearized operator L 4 Definition of a Fick-Relaxation operator.

4.1 The idea of the relaxation: link with the BGK operator.

The main idea to construct our model is to reproduce the above connections between the thermodynamic forces (and associated kinetic repeal coefficients) and the kinetic theory in a simplest way. We have seen that every linear small perturbations around local equilibrium states are moments of the perturbation around the local Maxwellian states and kinetic coefficients are only connected to the properties of the linear interaction operator. Kinetic coefficients can be interpreted as (functions of) relaxation rates associated to those moments and any simple model should reproduce this in the hydrodynamical limit. This is the principle that was shed in light in [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF][START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF]. Then if one thinks of relaxation models the second principle is that the corresponding relaxation function(s) must satisfy an entropy minimization under constraints that are exactly the relaxation of the moments under consideration. This has permitted to recover the ES-BGK operators for monatomic ( [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF]) and polyatomic gases ( [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]). Remark that if one just takes into account the conservation laws in the second principle one falls on the classical BGK model ( [START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF]).

In the present work we restrict our approach to the Fick and Newton laws.

Fick laws must actually be obtained from the kinetic theory using moments of the distribution functions in K ⊥ (see [START_REF] Kestin | Equilibrium and transport properties of the noble gases and their mixtures at low density[END_REF]. Let w r ∈ K ⊥ and consider the moment p j=1 R 3 f j w r,j . Then replacing the Boltzmann collision operator with a relaxation operator

R (f) = ν (G -f) , (28)
the above interpretation amounts to set

ν j=p j=1 R 3 (G j -f j ) w r,j = -λ r j=p j=1 R 3 f j w r,j , (29)
where λ r is a relaxation coefficient. The choice of a basis (w r ) r=1,...,p-1 of K ⊥ is crucial as well as that of the relaxation coefficients (λ r ) r=1,...,p-1 to obtain the Fick laws at the hydrodynamical limit (µ is a free parameter that will be adjusted later on to fit the viscosity). We proceed as follows. We firstly compute the relaxation coefficients λ r and corresponding moments w r such as to obtain the above coefficients. This computation is essentially based on two arguments (theorem 1):

1. the mass fluxes must depend on the moments M j g j w r,j of a small perturbation g around a local thermodynamical equilibrium M, 2. the relaxation coefficients λ r and corresponding moments w r are only connected to the property of the Fick coefficients matrix L ij .

Secondly G is obtained through the minimization of entropy principle under the above constraints plus the conservation laws (theorem 2).

Computation of the relaxation coefficients and moments.

Before we set the main result it is important to examine the properties of the Fick coefficients matrix L ij obtained from the linearization of the true Boltzmann operator.

Lemma 4. If it is computed thanks to the Boltzmann linearized operator, the matrix L ij , (i, j) ∈ [1, p] 2 of the kinetic coefficient has a rank p -1 and for any

k ∈ [1, p -1] the matrix L ij , (i, j) ∈ [1, k]
2 has a full rank.

We postpone the proof of this lemma to appendix. An important consequence of this lemma is the following result.

Lemma 5. The symmetric non positive matrix L * ij defined by

L * ij := k B L ij ||C i || ||C j || = k B L ij m i n i k B T m j n j k B T , ∀i, j ∈ [1, p] (30)
always diagonalizes in an orthonormal basis as:

L * = W T D * W
so that the corresponding eigenvalues (d * r ) r are non null for r = 1, ..., p -1 while d * p = 0. Moreover the vectors defined by

w r = p s=1 W rs C s ||C s || , r = 1, ..., p -1 (31)
form an orthonormal basis of C while the vector

w p = s ρ s ρ C s C s = ± s=p s=1 W ps C s C s ( 32
)
is orthogonal to C (in fact belongs to K) and normalized.

Proof. Since L * is symmetrical, it may be diagonalized. The fact that it has exactly one vanishing eigenvalue is linked to its rank which is p -1. Finally, we can always permute indexes so that d * p = 0. Consider the vector ̟ whose components are √ ρ j / √ ρ. Then a direct computation shows that:

(L * ̟) i = j=p j=1 k B L ij C i C j ρ j ρ = k B C i √ ρk B T j=p j=1 L ij m j = 0
the last equality being obtained thanks to [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF]. Besides, the vector ̟ is normalized for the usual vector norm. As a consequence, when diagonalizing the matrix L * in a orthornormalized basis, then the normalized eigenvector associated to d * p = 0 is necessary equal to ±̟. Then we have

W ps = ± √ ρ s / √ ρ. At last, as C s = √ n s k B T / √
m s , a direct computation shows that the vector w p has m i (vu) / √ ρk B T on its i th line, so belongs to K. Now, since the matrix W T is orthogonal and since the family

C i / C i , i ∈ [1, p]
is othonormalized for the L 2 (M) dot product, then the family w s , s ∈ [1, p] is also orthonormalized for the L 2 (M) dot product. Finally, in the family w s , s ∈ [1, p -1], every vector is a linear combination of the C i , i ∈ [1, p] but does not belong to K so it is in C and finally constitutes an orthonormalized basis of it. Now we can set our main result which relies the eigenvalues of the matrix L * to the relaxation coefficients.

Theorem 1. Let L be the linearized operator of R (see section 2.3). Assume that [START_REF] Brull | Problem of evaporation-condensation for a two component gas in the slab[END_REF] holds and that L -1 (B) and L -1 (A) are orthogonal to C. Let d * r and w r (r ∈ [1, p]) be defined as in lemma 5, then setting

λ r = d * r -1 , λ p = 0 ( 33
)
the BGK model that satisfies relaxation equations [START_REF] Sone | Doi Kinetic theory analysis of gas flows condensing on a plane condensed phase: Case of a mixture of a vapor and noncondensable gas[END_REF] allows to recover at the hydrodynamic limit the Fick laws meaning that the corresponding mass fluxes are given by

J i = j=p j=1 L ij ∇ -µ j T .
Proof.

Step 1 Let us consider the rescaled equation

∂ t f ǫ + v • ∇ x f ǫ = 1 ǫ R (f ǫ ) (34)
and perform a Chapman-Enskog expansion Assume that (10) holds then the zeroth order equation satisfied by f 0 = M yields as usual the Euler equation. As a consequence lemma 2 remains valid and we have

L (g) = j=p j=1 k -1 B (I -P K ) (C j ) • ∇ -µ j T + A : D (u) + B • ∇ 1 T . (35)
Next assuming that L -1 (B) and L -1 (A) are orthogonal to C (see definition 1) then equation [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF] implies that mass fluxes defined by

J i = g, C i
depend only on the values of (∇µ j /T ). Our aim is to obtain at the hydrodynamic limit the Fick coefficients L ij computed from the true Boltzmann equations [START_REF] Aoki | Perthame A consistent BGK-type model for gas mixtures[END_REF], that is

J i = g, C i = p j=1 L sj ∇ -µ j T . ( 36
)
Step 2 Consider now the first order expansion of relaxation equations ( 29)

L (g) , w r = -λ r g, w r , r ∈ [1, p], (37)
If we choose the moments w r as linear combinations of the basis (C i ) i then the right hand sides of the above equations can be written in terms of the required mass fluxes (36) while the left hand sides are determined thanks to (35). So we need to find relaxation coefficients and moments such that both sides of equations (37) fit. Let us compute those terms with our choice of moments (w r ) r (lemma 5). The right hand sides read

g, w r = p s=1 W rs ||C s || J s = p s,j=1 W rs ||C s || L sj ∇ -µ j T , ( 38
)
where W is the matrix of eigenvectors defined in lemma 5. Next the left hand sides can be computed thanks to ( 23)

L (g) , w r = p j,s=1 k -1 B W rs ||C s || C j , C s ∇ -µ j T = k -1 B p j=1 W rj ||C j || ∇ -µ j T , r = 1, ..., p -1 since (I -P K ) (C j ) , w r = C j , w r ∀w r ∈ C. As concerns the case r = p it is obvious that 1 √ ρk B T j L (g) , m j C j = 0.
To summarize (37) reads

k -1 B p j=1 W rj ||C j || ∇ -µ j T = -λ r p s,j=1 W rs ||C s || L sj ∇ -µ j T , r = 1, ..., p-1,
and since those equations must be valid for all values of (∇µ j /T ) j we have

W rj = -λ r p s=1 W rs L * sj , ∀r = 1, ..., p -1, j = 1, ..., p
where L * sj is defined in equation 30. We can also complete this set of equations with 0 = p s=1 W ps L * sj , j = 1, ..., p thanks to [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF] and the symmetry of L * sj . Multiplying both side of the above equations with W tj and summing over j we obtain

δ rt = -λ r p s,j=1 W rs L * sj W tj , ∀r = 1, ..., p -1, t = 1, ..., p 0 = p s,j=1
W ps L * sj W pj , ∀t = 1, ..., p which amounts to diagonalize L * and thus define λ r = -d * r -1 , r = 1, ..., p -1 and λ p = 0 in a unique way (see lemma 5).

Definition Fick-Relaxation operator.

Theorem 2. Let φ φ φ l , l ∈ [1, p + 4] be the natural basis for K and w r , r ∈ [1, p -1] be the orthonormalized basis for C to which we associate the relaxation values λ r , r ∈ [1, p -1]. For f a list non negative non null functions with suitables integration properties, consider the set K (f) defined as

g ∈ K (f) ⇔ ∀l ∈ [1, p + 4] , i=p i=1 R 3 φ l i (g i -f i ) dv = 0, ∀r ∈ [1, p -1] , i=p i=1 R 3 w r,i g i -1 -λr ν f i dv = 0. ( 39
)
Then there exist a unique minimizer G of the entropy function H under the constraint that G should belong to K (f). Moreover, if we note by n i , u i , T the hydrodynamic parameters associated to the function f, then such a minimizer reads as:

∀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 , ( 40
)
where the velocities u i are defined by the relations

U -U = N -1 W T I - 1 ν Λ W N U -U . ( 41
)
Here U = (u 1 , ..., u p ) T , U = (u 1 , ..., u p ) T and U = (u, ..., u) T . N and Λ denote respectively the diagonal matrix with diagonal terms ( √ ρ 1 , . . . , √ ρ p ) and (λ 1 , . . . , λ p ).

Proof. The first set of constraints amounts to impose the conservation of densities for each species and the global conservation of momentum and energy. Let us now perform the computation of velocities (u i ) i . The set of constraints related to the velocities are firstly those involving moments (w r ) r=1,...,p-1 and secondly those expressing conservation of momentum (moments φ φ φ i=p+1,p+2,p+3 ). Remark that the later can be written as

i=p i=1 R 3 w p,i (g i -f i ) dv = 0
thanks to the definition of w p (32). Then the whole set of constraints on the moments (w r ) r reads

p j=1 W rj 1 ||C j || G j (v -u)dv = (1 - λ r ν ) p j=1 W rj 1 ||C j || f j (v -u)dv ⇐⇒ p j=1 W rj n j ||C j || (u j -u) = (1 - λ r ν ) p j=1 W rj n j ||C j || (u j -u), r = 1, ..., p ⇐⇒ p j=1 W rj √ ρ j (u j -u) = (1 - λ r ν ) p j=1 W rj √ ρ j (u j -u), r = 1, ..., p.
This system of equations can be written

W N (U -U) = I - 1 ν Λ W N U -U
which gives exactly the result of the theorem.

Next we must prove the existence of a (unique) minimizer of the entropy in the set K (f). Roughly speaking this is a generalization of the complex problem of minimization of the entropy under moment constraints (see [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF][START_REF] Schneider | Entropic approximation in kinetic theory[END_REF]). For example we can follow the approach Junk ( [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF]) which consists in examining the set

{(α l ) l , (β r ) r / exp l=p+4 l=1 α l φ l i + p-1 r=1 β r w r,i dv < +∞}.
Then Junk proved that if this set is open then the minimization problem has a unique solution as soon as the set of constraints is realizable in the sense that there exists a nonnegative function in K (f). Here the existence of such a function is obtained simply by exhibiting the function (40). Moreover the solution to the minimization problem must be of the exponential form and so it is G itself which ends the proof.

This result enables us to define our Fick-Relaxation operator.

Definition 3. We define the Fick-Relaxation operator R (f) by the relation

R (f) = ν (G -f) , (42)
where ν > 0 and

G = min{H(g), s.t. g ∈ K (f)}.
Remark 2. In the definition of K(f) given in equation ( 39), the first constraint leads to conservation laws whereas the second one gives the Fick law through the determination of λ r .

Remark 3. Remark that

H(G) = i n i ln n i - 3 2 ln m i 2πk B T -n
depends only on n i , m i and the temperature T of the mixture which is defined in (2). This value does not depend on each velocity u i as long as i ρ i u i = ρu. This means that it is also the minimizer of the entropy in the broader set K (f) defined by

K (f) = {g ≥ 0 a.e/ i=p i=1 R 3 φ l i (g i -f i ) dv = 0 ∀l ∈ [1, p + 4]} ( 43
)
where there is no unicity of the solution. All the minimizers are function of the form (40) with i ρ i u i = ρu. So the unicity of the solution is only the sake of specifying the set of velocity (u i ) i .

5 Properties of the Fick-Relaxation operator.

In this part we firstly prove the fundamental properties of the Fick-Relaxation operator. Next we establish that the Fick-Relaxation operator is properly defined and we compute the linearized operator. From this we easely deduce all fluxes or kinetic coefficients. The end of the section is devoted to a discussion about the indifferentiability property. More precisely Proposition 3 shows an alternative property for a BGK model between fitting Fick coefficients at the hydrodynamic limit and satisfying an indifferentiability property.

Conserved moments and H-theorem.

Here we want to establish the physical properties (7,[START_REF] Brull | Problem of evaporation-condensation for a two component gas in the slab[END_REF][START_REF] Brull | The Boltzmann equation for a two component gas in the slab for soft forces[END_REF][START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF] of the Fick relaxation operator R (f).

Theorem 3. The Fick-Relaxation operator satisfies the general properties 1. Equilibrium states:

R (f) = 0 ⇔ ∃n i , u, T s.t. ∀i ∈ [1, p] , f i = M i .
2. H-theorem:

∀f, f i ≥ 0, i=p i=1 R 3 R i (f) ln (f i ) dv ≤ 0.
3. Characterization of equilibrium:

R 3 R i (f) ln (f i ) dv = 0 ⇔ ∃n i , u, T s.t. ∀i ∈ [1, p] , f i = M i .

Conservation laws:

∀f,

f i ≥ 0, ∀Φ, i=p i=1 R 3 R i (f) Φ i dv = 0 ⇔ Φ ∈ K.
Proof. Property 1: Equilibrium states Taking into account the Maxwellian form of f = G we are let to prove that U = U =⇒ U = U. Coming back to the definition of U (41) we have

U = U =⇒ 0 = W T N -1 Λ W N U -U .
The second equation reduces to

0 = Λ W N U -U since the matrix N -1 W T is of rank p. Remark that λ r = 0 for all r ∈ [1, p -1]
so the above equation reads

0 = i W ri √ ρ i (u i -u), ∀r ∈ [1, p -1].
As concerns the case r = p recall that (32

) i W pi √ ρ i (u i -u) = 1 √ ρ i ρ i (u i -u) = 0
since u is the mean velocity of f. Gathering the two set of equations we have

W N U -U = 0 which implies that U -U = 0 since W N is of rank p.
Property 2: H theorem: The convexity of the function x → x ln x implies

i=p i=1 R 3 (G i -f i ) ln f i dv ≤ H (G) -H (f) . ( 44
)
Now since we have i ρ i u i = ρ u (which is self-evident), Remark 3 tells us that the function G f defined by

G f,i = n i (2πk B T /m i ) 3/2 exp - m i v -u i 2 2k B T , ∀i ∈ [1, p] (45)
is a minimizer of the entropy in the set K (f) (43) but also in the set K (f) defined as

g ∈ K (f) ⇔ ∀l ∈ [1, p + 4] , i=p i=1 R 3 φ l i (g i -f i ) dv = 0, ∀r ∈ [1, p] , i=p i=1 R 3 (g i -f i ) v dv = 0. Thus i=p i=1 R 3 (G i -f i ) ln f i dv ≤ H (G) -H (G f ) = 0. ( 46
)
Property 3: Characterization of Equilibrium: The nullity of entropy flux implies

0 = i=p i=1 R 3 (G i -f i ) ln f i dv ≤ H (G) -H (f) ≤ 0.
Next the second inequality is strict unless f = G f (see above) so that we have

i=p i=1 R 3 (G i -G f,i ) ln G f,i dv = 0.
A simple computation of this term gives

i=p i=1 R 3 G i (v -u i ) 2 dv - 2n i k B T m i = i=p i=1 n i (u i -u i ) 2 = 0,
where we have used the conservation of mass per species. Hence we must have U = U which implies (see property 1) U = U.

Property 4: Conservation laws: Thanks to the definition of G through a minimization problem with qualified constraints (39) the implication ⇐ is obvious. The converse may be obtained from the H-theorem. Assume that there holds:

∀f s.t.∀i, f i ≥ 0, i=p i=1 R i (f) φ i = 0
where φ is given. Then consider the particular f given as f i = exp (φ i ) ⇔ φ i = ln (f i ). Then from the characterization of the equilibrium states we get that:

∃n i , u, T s.t. ∀i ∈ [1, p] , f i = M i = exp (φ i ) so that finally φ ∈ K.
Remark 4. The stategy to obtain a relaxation operator giving exact Fourier coefficient would be to generalize the ES-BGK model ( [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF][START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF][START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF][START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]) to a mixture setting. This is however still not enough to ensure that the Fick Relaxation operator is properly defined and we need to compute the linearized operator.

5.2 Linearization of Fick-Relaxation operator.

Explicit expression.

This paragraph is devoted to the computation of the linearized operator associated to the Fick-Relaxation operator defined in (42).

Theorem 4. The linearized operator L of the Fick-Relaxation operator around a thermodynamical equilibrium writes

L = ν (P K + R • P C -I) ,
where I is the identity operator on L 2 (M), P K , P C respectively the orthogonal projection on K, C while R is the linear operator defined on C by the formula:

∀r ∈ [1, p -1] , R (w r ) = 1 -λr ν w r . Proof. Let M be a thermodynamical equilibrium then L i (g) = 1 M i lim ǫ→0 R i (M(1 + ǫg)) -R i (M) ǫ ∀i ∈ [1, p],
where g is any function in L 2 (M). Remark that thanks to theorem 3 R i (M) = 0. Next let g be fixed then there exists ǫ 0 such that for all ǫ < ǫ 0 the problem

G = min{H(h) s.t. h ∈ K (M (1 + ǫg))}
admits a solution. Indeed the set of realizable constraints K (M) is open [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF].

Let us write this solution in the following form

G i,M(1+ǫg) = exp l=p+4 l=1 α l (ǫ, g) φ l i + p-1 r=1 β r (ǫ, g) w r,i
where the pseudo-Lagrange multipliers (α l (ǫ, g)) l , (β r (ǫ, g)) r are solutions of the perturbated problem. We may compute exactly those coefficients from the form (40) given in theorem 2 and then find out that they are infinitely differentiable with respect to each variables (for ǫ < ǫ 0 ). Thus they can be expanded in ǫ

α l = α 0 l + ǫα 1 l (g) + O(ǫ 2 ) ∀l, β r = β 0 r + ǫβ 1 r (g) + O(ǫ 2 ) ∀r.
The zeroth order corresponds to the Maxwellian distribution that is

α 0 i = -µ i k B T -m i u 2 i 2k B T , α 0 p+3 = u i k B T , α 0 p+4 = - 1 2k B T α 0 l = 0 ∀l = i, p + 3, p + 4, β 0 r = 0 ∀r,
where µ i is defined in [START_REF] Gross-M | Krook Model for collision processes in gases: smallamplitude oscillations of charged two-component systems[END_REF]. Then G i,M(1+ǫg) writes

G i,M(1+ǫg) = M i 1 + ǫ p+4 l=1 α 1 l (g)φ l i + p-1 r=1 β 1 r (g) w r,i + O(ǫ 2 ) . (47)
So that

L i (g) = p+4 l=1 α 1 l (g)φ l i + p-1 r=1 β 1 r (g) w r,i -g. (48)
Now remark that the set of constraints (39) depends linearly on the moments of the function M(1 + ǫg) so that the perturbations at order 1 in ǫ of the Lagrange multipliers α 1 l (g), β 1 i (g) are also linear. We can compute them exactly by inserting the expression (47) in the constraints (39) with f = M(1 + ǫg) and considering the equalities at order one in ǫ we get

p+4 l=1 φ φ φ k , φ φ φ l α 1 l (g) = g, φ φ φ k ∀k ∈ [1, p + 4] , p-1 r=1 w s , w r β 1 r (g) = 1 - λ r ν g, w r ∀s ∈ [1, p -1] .
Here the first line corresponds exactly to the projection of g on K through the expression (48) of L i (g) while the second line amounts to project on the basis (w r ) r of C and multiply each coordinates with a factor (1 -λr ν ).

Properties of the linearized operator.

In this part we easely deduce from the explicit expression of the linearized operator the following important properties.

Proposition 1. The linearized operator associated to the Fick-Relaxation operator satisifes the following properties

i) KerL = K, ii) L is Fredholm, iii) L is self adjoint negative on K ⊥ .
Now we can finish the work by showing the last theorem.

Theorem 5. The pseudo inverse L -1 of the Fick-Relaxation operator reads as:

(∀g ∈ K ⊥ ), L -1 (g) = 1 ν (R -I C ) -1 • P C + (P C -I) (g) , (49)
where I C denotes the restriction on C of the identity operator, P C the projection on C and R has been defined in Theorem 4.

Proof. For g ∈ K ⊥ it suffices to compute L -1 • L (g) .

Hydrodynamical limit.

Here we want to compute the hydrodynamical limit at order 1 of the system of equations

∀i ∈ [1, p] , ∂ t f ǫ i + v • ∇ x f ǫ i = 1 ǫ (G ǫ i -f ǫ i ) ν . (50)
Theorem 6. Let ν be defined as

ν = nk 2 b T 2 µ
where µ = L uu is the viscosity of the mixture (see equation ( 26)). Then the hydrodynamical limit of (50) with R(f) defined in 28 gives at the order 1 the Navier-Stokes equations [START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF][START_REF] Chapman | The mathematical theory of non uniform gases[END_REF][START_REF] Ern | Multicomponent Transport Algorithms[END_REF] with the fluxes given by

J i = p j=1 L ij ∇ -µ j T , J q = L R qq ∇ 1 T , J u = ν D (u) . (51)
Here (L ij ) ij is the matrix of Fick coefficients [START_REF] Kestin | Equilibrium and transport properties of the noble gases and their mixtures at low density[END_REF] computed with the (linearized) Boltzmann equation and

L R qq = 5µk b T 2n p i=1 n i m i .
Proof. The proof is easely deduced from the fact that R(f) is properly defined according to theorem 3 and proposition 1. Hence (see remark 1) all kinetic coefficients can be computed as soon as we know the exact form of L -1 which is actually given in (49).

1. Density fluxes: We have directly

L -1 (A) = -ν -1 A and L -1 (B) = -ν -1 B ∈ C ⊥
so that theorem 1 applyes. Hence

J i = j=p j=1 L ij ∇ -µ j T ∀i ∈ [1, p]
meaning that L iq = L qi = 0 for all i ∈ [1, p] (L -1 is self-adjoint).

2. Viscosity: it is possible to recover the true viscosity µ = L uu defined by [START_REF] Kosuge | Model Boltzmann equation for gas mixtures: Construction and numerical comparison[END_REF] by adjusting the relaxation rate ν of R(f). Let L R uu be the viscosity of our model then

L R uu I ⊗ I = L -1 (A) , A = 1 ν A, A = n k 2 B T 2 ν I ⊗ I.
So obviously letting ν = n k 2 b T 2 /µ we have L R uu = µ. 3. Fourier coefficient: the "true" Fick and viscosity coefficients being recovered it is not possible to obtain the correct value of L qq [START_REF]Morse Kinetic model equations for a gas mixture[END_REF] with our model. Indeed

L R qq I = L -1 (B) , B = 1 ν B, B ⇔ L R qq = 5 µ k B T 2 n p i=1 n i m i = L qq I.
Here the density fluxes are deduced from our "natural" construction (theorem 1). That is using the concept of relaxation rates associated to ad-hoc moments of the distribution functions. To recast the hydrodynamical limit in the spirit of section 3.4 we may check the following results.

Proposition 2. The linearized operator L satisfies

∀ (i, j) ∈ [1, p] 2 , k -1 B L -1 (I -P K ) (C i ) , (I -P K ) (C j ) = L ij . (52)
Proof. First let us note the one to one correspondance between the C i and the w k as:

∀i ∈ [1, p] , C i = C i k=p k=1 W T ik w k .
From this, it is very easy to compute (I -P K ) (C i ) which is simply given as:

(I -P K ) (C i ) = C i k=p-1 k=1 W T ik w k .
Applying now the operator L -1 is now straighforward thanks to the explicit formulation and we have

L -1 (I -P K ) (C i ) = -C i k=p-1 k=1 W T ik λ k w k .
As a consequence the LHS of 3.4 reads as

-k -1 B C i C j k,l=p-1 k,l=1 1 λ k W T ik W T jl w k , w l .
Since the basis w k , k ∈ [1, p] is orthonormalized, then the sum simplifies as:

-k -1 B C i C j k=p-1 k=1 W T ik 1 λ k W kj = k -1 B C i C j k=p k=1 W T ik d * k W kj (this is because d * k = -1/λ k , k ∈ [1, p -1], while d * p = 0). That is -k -1 B C i C j k=p-1 k=1 W T ik 1 λ k W kj = C i C j k -1 B L * ij = L ij
which ends the proof.

About the indifferentiability property

Following ( [START_REF] Aoki | Perthame A consistent BGK-type model for gas mixtures[END_REF]) the indifferentiability property may be stated as follows: when in the mixture all the masses and differential cross sections are equal -which we call indifferentiable molecules-then the BGK operator for the mixture must collapse into a BGK one for a single species, that is:

∀i, j ∈ [1, p] 2 , m i = m, σ ij = σ ⇒ i=p i=1 R i (f) = R i=p i=1 f i .
This "factorization" property is an algebraic one which holds when the dilute Boltzmann operator is chosen to model the particles collisions. More precisely there is:

∀i, j ∈ [1, p] 2 , m i = m, σ ij = σ ⇒ i=p i=1 Q i (f) = Q i=p i=1 f i .
Note however that such a factorization property holds because of the bi-linearity of the dilute Boltzmann operator. But it has no particular reason to hold when the Boltzmann operator features a cubic dependancy with f which happens for instance when one considers dense situations [START_REF] Kurochkin | Enskog Kinetic Equation Modified for a Dense Soft-Sphere Gas[END_REF].

Let us now consider the following quite general form of a BGK operator for gas mixtures which we denote with

R ∀i ∈ {1, • • • , p}, R i (f) = νi ( G i -f i ). ( 53
)
Here the values of the relaxation frequencies νi may be different. We first set some quite reasonable assumptions on the above model:

1. The constraints defining the mean velocities of ( G i ) as functions of the mean velocities of (f i ) i are of the form i=p i=1

X ij R 3 G i (v -u) dv = i=p i=1 T ij R 3 f j (v -u) dv
where the matrix X and T ∈ R p × R p do not depend on v. Or in short

U -U = X -1 N -1 T N U -U (54)
(here we have kept the notations of theorem 2).

2. In the indifferentiability situation the relation u i = u holds for all i.

Let us remark that up to our knowledge the linearity of the constraints on the basis (C i ) i is the rule for all the BGK operators we have found. Those constraints are quite natural regarding to the way the true diffusion (Fick) coefficients are obtained (see [START_REF] Kestin | Equilibrium and transport properties of the noble gases and their mixtures at low density[END_REF]). Otherwise if those constraints are not linear they must be of integer degree for the sake of Galilean invariance..

The second assumption is required in most cases because the indifferentiability is obtained thanks to an "additivity" property. That is

i νi ( G i -f i ) = ν( G -f)
when all masses and cross sections are equal. For example such assumptions are easely found when G i is not a linear function of its mean velocity u i .

We are now going to draw some conclusion about the compatibility of the indifferentiability property and the Fick law for such a model. Lemma 6. Let R be a properly defined operator of the form (53) satisfying the conditions 1 and 2 above. Then when all molecules are the same the restriction to C of the linearized BGK is proportional to minus identity.

Proof. The property [START_REF] Brull | The stationary Boltzmann equation for a two component gas in the slab for different masses[END_REF] implies that the restriction of the linearized operator on the space spanned by (C i ) i actually reduces to C. Thus generalizing the calculations we have done in section 5.2.1 and using assumption 1 this restriction reads R -I C where R is a linear operator depending on T and I C is the restriction on C of the identity operator. To be consistent with the symmetry of the linearized Boltzmann operator the operator R has to be symmetrical on C while the operator R -I C must be non positive and symmetric because of the property [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]. Consequently R can be diagonalized in a proper orthogonal basis

x r , r ∈ {1, • • • , p -1} such that ∀r ∈ {1, • • • , p -1}, R (x r ) = (1 -α r )x r , α r ≥ 0.
And we can write (54) in the same form as((41)

U -U = N -1 X T (I -A) X N U -U
Here the matrix X T denotes the passage from the basis (C i ) i to the basis (x 1 , ..., x p-1 , m i C i )) and the matrix A is diagonal with a rr = α r , r ∈ {1, • • • , p-1} and a pp = 0. Next according to the second assumption the indifferentiability property requires u i = u for any i ∈ {1, • • • , p} in the limit of same masses and cross sections. As a consequence the relation (54) must be valid for each values of u i meaning that for any r ∈ {1, • • • , p}, α r = 1. Consequently R = 0.

Next we prove the following alternative. Proposition 3. Two situations are possible.

• Either there is a properly defined BGK operator of the form (53) satisfying conditions 1 and 2 above, the indifferentiability property and giving the Fick law at the hydrodynamical limit. In that case the Fick Relaxation operator R (definition 3) satisfies the indifferentiability as well. Moreover when all masses and cross sections the operator must obey a law in the indifferentiable case that is νi = ν = ν ∀i.

• Or there are no well defined BGK operator satisfying the conditions 1) and 2) and in the same time indifferentiability and Fick law.

Proof. Consider a properly defined BGK operator defined (53) satisfying Fick law and the indifferentiability condition. Its associated Fick coefficients L e R ij are given by ( 24) replacing L B with L e R (see remark 1). Thus according to lemma 6 when all masses and all differential cross sections are equal

L e R ij I = - 1 k B ν (I -P K ) (C i ) , (I -P K ) (C j ) , ν > 0. ( 55
)
If the model satisfies the Fick law as well we have (L e R ij ) i,j = (L ij ) i,j where the second matrix is defined either by [START_REF] Kestin | Equilibrium and transport properties of the noble gases and their mixtures at low density[END_REF] or by

L ij I = - 1 k B ν (R -I C ) -1 (I -P K ) (C i ) , (I -P K ) (C j ) ,
where we have used (49) and ( 52). The comparison between this equation and (55) imposes firstly that R = 0. Coming back to equation (41) this means that U = U whatever are the values of U. Hence using the definition of G i (40) the Fick-Relaxation model satisfies the indifferentiability property. But according to lemma 6, in the situation of indifferentiability, L ij is given by a formula analogous to (55) leading to ν = ν.

The last assertion is self evident.

In a situation of same masses and same differential cross sections, we do not know exactly if the kinetic coefficients obtained by the linearized Boltzmann operator can match those obtained by a BGK satisfying the indifferentiability property. But in reality we do not care because the correcteness of the linearized BGK operator is much more important, it seems, that the indifferentiability property and if they contradict, we prefer to sacrify indifferentiability.

In concrete situations masses and cross sections are always different. Therefore the indifferentiability principle is only a continuity principle. That is why we prefer to conserve the correcteness of the hydrodynamical coefficients which are more physically pertinent.

Conclusion and perspectives.

In this paper we have introduced a new relaxation operator for gas mixtures. Our construction features two ideas. The first one is that we have constructed it within the framework of moments relaxations that were introduced in [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF][START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF] as a new way to understand the ellispodal model for monospecies systems. Here we have also taken into account the relaxation on velocities species which is necessary to obtain mass diffusion, thus shedding in light the space of vectors of moments of order one orthogonal to the space of collisional invariants. The second characteristic of this work is that we have focused the construction of this relaxation operator so that its associated linear mass and momentum transports coefficients exactly match those obtained by the Chapman-Enskog expansion applied to the Boltzmann operator. While it is classical to have such concerns to recover the right Newtown viscosity and Fourrier coefficient when constructing BGK models in the case of a single specie, it is not the usual track which is followed when considering gas mixtures as we pointed it in introduction. Interestingly our relaxation operator also enjoy important properties of the Boltzmann collision terms (H theorem for instance) which are not always satisfied when pure linear models are at stake. Besides we have pointed out the remarkable link exisiting between the eigenvalues of Fick coefficient matrix and the rate of relaxation for the vanishing moments associated to functions of degree one in velocity. Up to our knowledge such a perspective has never been expressed so clearly. Finally a key aspect of the Fick Relaxation Operator we have defined lies in the very simplicity of its computation which only requires the diagonalization of the Fick matrix! It seems to us that this could be an important thing for future users if any. Unfortunately when writing these lines, we are not still able to propose a Relaxation operator that could be able to match in the linear regime all the transport coefficients. It seems to us that obtaining simultaneously the correct Fick, Newtown and Fourrier coefficients could be possible by generalizing for gas mixture the approach described in [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]. More difficult however seems the possibility to get at the same time the correct (so called "cross") kinetic coefficients (L iq ) i (18), [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]. In our model (as it is also the model of [START_REF] Aoki | Perthame A consistent BGK-type model for gas mixtures[END_REF]) such cross kinetic coefficients are zero and give uncorrect Soret and Dufour coefficients (except for Maxwellian molecules). This might be a problem for situations featuring non negligible thermodiffusion for instance. Finally the approach we have described in this paper requires the knowledge of the kinetic coefficients. Unfortunately such a knowledge is not obvious at all but combining both the numerical algorithm to compute them given in [START_REF] Ern | Multicomponent Transport Algorithms[END_REF] and the experimental work of Kestin [START_REF] Kestin | Equilibrium and transport properties of the noble gases and their mixtures at low density[END_REF] will be of precious help when actually computing the Fick Relaxation Operator. Now assume that for k ∈ [1, p] there holds i=k i=1 α i (I -P K ) (C i ) = 0 then, working on the k first line of this equation we get that:

(∀l ∈ [1, k]) i=k i=1 δ il - n i m l ρ α i = 0.
To know if this linear system on the α i has the unique trivial solution or not (and then to know if the family (I -P K ) (C i ) , i ∈ [1, k] is independant or not), one has to discuss the determinant of the k by k matrix M defined as:

(∀i, l ∈ [1, k] 2 ) M il = δ il - n i m l ρ .
Elementary algebraic calculations show that the determinant of such a matrix may be computed as

det M = i=k i=1 ρ i ρ det N, with N =        ρ ρ 1 0 • • • 0 0 ρ ρ 2 0 . . . . . . • • • . . . • • • 0 • • • 0 ρ ρ k        -        1 1 • • • 1 1 . . . • • • 1 . . . • • • • • • . . . 1 • • • • • • . . .
Now assume by contradiction that x = 0. In that case for example x 1 = 0. According to (57), for any i ∈ {1, k}, x i = 0. So (56) implies that Therefore for any j ∈ {k + 1; n}, ρ j = 0 and we get a contradition. Then it follows that N is invertible.

Next consider the case where k = p. Let x = (ρ 1 , . . . , ρ p ) = 0. In that case

ρ i ρ ρ i = k j=1 ρ j .
Then x = 0 and x ∈ KerN . So N is non invertible.

7.2 Proof of Lemma 2.

Proof. The first problem consists in computing the left-hand side of ( 22), namely any of the (∂ t + v • ∇ x ) [M i ] for i ∈ [1, p]. Now introduce the following decomposion of M i given as:

∀i ∈ [1, p] , M i = exp -µ i k B T -m i u 2 2k B T + u k B T • m i v - m i v 2 2k B T .
Then, introducing the natural basis of K, φ φ φ k , k ∈ [1, p + 4] there holds:

∀i ∈ [1, p] , (∂ t + v • ∇ x ) (M i ) = k=p+4 k=1 [(∂ t + v • ∇ x ) α k ] φ φ φ k i M i (58)

.

  Consider first the case where k < p. Let x ∈ Ker(N ). Hence by definition of N , it holds that (∀j ∈ {1; k}),

R 3 (

 3 But since w := (M 1 , • • • , M p ) satisfies the Euler equations, then we know∀l ∈ [1, p + 4] , i=p i=1 ∂ t + v • ∇ x ) (M i ) φ l i dv = 0. (59)By introducting now the left-hand side expression (58) we obtain that∀l ∈ [1, p + 4] , k=p+4 k=1 (∂ t + v • ∇ x ) (α k ) φ k , φ l = 0. (60)

If we denote by S the symmetrical matrix such that S lk = φ l , φ k then for for any vector b ∈ L 2 (w (v) dv) we have clearly

the orthogonal projection of b on K. Let us consider now

The equation (60) gives that ∀l ∈ [1, p + 4] , b 1 , φ l = -b 2 , φ l . So using the explicit expression of the orthogonal projection on K (61) we have necessary that:

Using this equality with the relation (58) enables us to see the equation [START_REF] Holway | New statistical models for kinetic theory: methods of construction[END_REF] as:

Now we need to establish the detailed expression of the right-hand side of (63). Since α k does not depend on the variable v, then k=p+4 k=1

and finally

Now multiplying the former terms by the functions φ k , summing over k and gathering the terms in

T and the others we get that the

given as the sum of the following terms:

Thanks to the symmetry of (vu) ⊗ (vu) we have

Finally, by computing the orhtogonal projection on K of the vectors having respectively

on their i th line, we get exactly the expression 22.

Proof of Lemma 3

Proof. The equation which is satisfied by the perturbation (list of) function g reads for any i ∈ [1, p] as:

where C i is the contribution on the i th line of the collision operator (Boltzmann or BGK one). Multiplying such an equation by φ l i , integrating over R 3 , summing on i and taking into account that the perturbation g is in K ⊥ gives us:

which is exactly the Euler equation which are perturbated by the flux terms v ⊗ φ φ φ l , g that need to be clarified. First let us remark that u ⊗ φ φ φ l , g = 0 so we can replace v by vu when computing v ⊗ φ φ φ l , g . To have the diffusion term associated to the density n l , we must considerate the list of function (vu) φ φ φ l = C l , l ∈ [1, p] and we get directly the following expression J l = C l , g = (I -P K ) C l , g : the fact tha t we can replace C l by (I -P K ) C l is the bracket is due to the fact that g ∈ K ⊥ Following the same idea, the diffusion of momentum is obtained when φ l i = m i v. Since g ∈ K ⊥ this may be replaced by m i (vu) so finally we get that J u = A, g . To get heat diffusion, we need to consider φ l i = 1 2 m i v 2 . This gives us an energy diffusive fluxgiven as B, g .

Proof of Lemma 4.

Proof. Note by L :,j the j th column of the matrix L ij , (i, j) ∈ [1, k] 2 and assume that j=k j=1 α j L :,j = 0. Then there is i=k i=1 j=k j=1 α i α j L ij = 0 when the matrix L ij comes from a Chapman-Enskog expansion this means:

As the operator L -1 is negative on

Now using the lemma 1, for any k < p, this implies that ∀i ∈ [1, k] , α i = 0. For k = p, one can find from lemma (1) a set of real β i , i ∈ [1, p] which are not all zero such that:


F R Petruzielo

Julien Toulouse
email: julien.toulouse@upmc.fr

C J Umrigar
email: cyrusumrigar@cornell.edu

Basis set construction for molecular electronic structure theory: Natural orbital and Gauss-Slater basis for smooth pseudopotentials

A simple yet general method for constructing basis sets for molecular electronic structure calculations is presented. These basis sets consist of atomic natural orbitals from a multi-configurational self-consistent field calculation supplemented with primitive functions, chosen such that the asymptotics are appropriate for the potential of the system. Primitives are optimized for the homonuclear diatomic molecule to produce a balanced basis set. Two general features that facilitate this basis construction are demonstrated. First, weak coupling exists between the optimal exponents of primitives with different angular momenta. Second, the optimal primitive exponents for a chosen system depend weakly on the particular level of theory employed for optimization. The explicit case considered here is a basis set appropriate for the Burkatzki-Filippi-Dolg pseudopotentials. Since these pseudopotentials are finite at nuclei and have a Coulomb tail, the recently proposed Gauss-Slater functions are the appropriate primitives. Double-and triple-zeta bases are developed for elements hydrogen through argon. These new bases offer significant gains over the corresponding Burkatzki-Filippi-Dolg bases at various levels of theory. Using a Gaussian expansion of the basis functions, these bases can be employed in any electronic structure method. Quantum Monte Carlo provides an added benefit: expansions are unnecessary since the integrals are evaluated numerically.

I. INTRODUCTION

In quantum chemistry (QC) calculations, molecular orbitals are traditionally expanded in a combination of primitive Gaussian basis functions and linear combinations of Gaussian primitives called contracted basis functions [1]. These basis sets cannot express the correct molecular orbital asymptotic behavior but are used in QC calculations to permit analytic evaluation of the twoelectron integrals [2].

Analytic integral evaluation significantly limits flexibility in basis set choice but is essential for computational efficiency in QC calculations. However, in practice, other basis function forms can be considered since an arbitrary function can be expanded in Gaussians. Of course, the fidelity of this representation is limited. An expansion in a finite number of Gaussians cannot reproduce the exponential decay of the wavefunction at large distances or the Kato cusp conditions [3] at nuclei, but it can mimic these features over a finite range.

Quantum Monte Carlo (QMC) calculations [4] offer greater freedom in choice of basis functions because matrix elements are evaluated using Monte Carlo integration. Consequently, the correct short-and long-distance asymptotics can be satisfied exactly. For systems with a divergent nuclear potential, Slater basis functions can exactly reproduce the correct electron-nucleus cusp and long-range asymptotic behavior of the orbitals. For calculations on systems with a potential that is finite at the nucleus and has a Coulomb tail, Gauss-Slater (GS) primitives [5] are the appropriate choice since they introduce no cusp at the origin and reproduce the exponential long-range asymptotic behavior of the orbitals. Despite shortcomings, traditional QC basis sets have yielded good results. The natural orbitals (NOs) from a post Hartree-Fock (HF) method are a particularly successful form of contracted function [6][7][8][9]. The simplest NO construction involves diagonalizing the one-particle density matrix from a ground state atomic calculation [6]. This construction is unbalanced due to obvious bias favoring the atom. More complicated constructions involve diagonalizing the average one-particle density matrix of several systems: atomic ground and excited states, ions, diatomic molecules, and atoms in an external electric field [7][8][9]. These constructions produce excellent results, but they are complex.

A simple but general method for constructing basis sets for molecular electronic structure calculations is proposed and tested here. The bases are combinations of the NOs obtained from diagonalizing the one-particle density matrix from an atomic multiconfigurational self-consistent field (MCSCF) calculation and primitive functions appropriate for the potential in the system. The primitives are optimized for the homonuclear dimer in coupled cluster calculations with single and double excitations (CCSD), with the intention of producing a balanced basis set. Importantly, optimal exponents for the primitive functions are shown to depend weakly on the level of theory used in the optimization. Additionally, results

show that coupling is weak between primitive functions of different angular momenta. This enables efficient determination of optimal exponents.

The utility of the above construction is demonstrated for the elements hydrogen through argon with the nondivergent pseudopotentials of Burkatzki, Filippi, and Dolg (BFD) [10]. Since these pseudopotentials are finite at the nuclei and have a Coulomb tail, the GS functions are the appropriate primitives. These pseudopotentials are chosen for demonstrated accuracy in all cases tested and because they are accompanied by a basis set. The BFD basis [10] serves as a metric for testing the new basis. The benefits of our bases extend to all electronic structure methods tested, including CCSD, HF, the Becke three-parameter hybrid density functional (B3LYP) [11], and QMC.

The main area of interest for the authors is QMC. Since QMC results depend less on basis set than traditional QC methods [5], only double-zeta (2z) and triple-zeta (3z) bases are presented.

This paper is organized as follows. Basis function form and properties are demonstrated in Sec. II. Results for calculations with the new bases are discussed in Sec. III. Concluding remarks are provided in Sec. IV. Supplementary material is provided on EPAPS [12].

II. BASIS SET

The number of basis functions for each angular momentum follows the correlation consistent polarized basis set prescription of Dunning [1]. 2z and 3z bases appropriate for the BFD pseudopotentials are generated for the elements hydrogen through argon. Since the BFD pseudopotential removes no core for hydrogen and helium, the 2z basis for these elements consists of two S functions and one P function, while the 3z basis consists of three S functions, two P functions, and one D function. Since the BFD pseudopotential removes a helium core for the first row atoms and a neon core for the second row atoms, the remaining elements lithium through argon have the same number of basis functions. In particular, the 2z basis consists of two S functions, two P functions, and one D function, while the 3z basis consists of three S functions, three P functions, two D functions, and one F function.

The bases consist of a combination of contracted and primitive functions. Since the BFD pseudopotentials are finite at the origin and have a Coulomb tail, the GS functions are the appropriate primitives. With the exception of the elements in Group 1A of the periodic table (i.e. H, Li, and Na), the basis for each element includes a single S contraction and a single P contraction combined with an appropriate number of GS primitives. Only two contractions are employed to reduce the computational cost of using this basis in QC calculations. Since elements in Group 1A of the periodic table have only one electron for the BFD pseudopotentials, a single S orbital is the ground state wavefunction, and this can be obtained exactly in HF. Thus, the basis for each element in Group 1A includes a single S contraction, no P contractions, and an appropriate number of GS primitives.

A. Contracted Functions

A contracted basis function is a linear combination of Gaussian primitives:

ϕ nlm (r, θ, φ) = i c i 2(2α i ) n+ 1 2 Γ(n + 1 2 )
r n-1 e -αir 2 Z m l (θ, φ),

(1) where r, θ, φ are the standard spherical coordinates, n is the principal quantum number, l is the azimuthal quantum number, m is the magnetic quantum number, Z m l (θ, φ) is a real spherical harmonic, c i is the i th expansion coefficient, and α i is the i th Gaussian exponent. In practice, the restriction n = l + 1 applies.

The exponents of the primitive functions that form the contracted basis functions are determined as follows. For each angular momentum for which a contraction is desired, an uncontracted basis consisting of nine eventempered primitive Gaussians is generated. For each set of uncontracted Gaussians, the minimum exponent and even-tempering coefficient are varied to minimize the CCSD energy of the atom using a Python wrapper around GAMESS [13].

An assumption of weak coupling between the different angular momenta underlies the optimization procedure. Consequently, the uncontracted basis for each angular momentum is optimized separately. This optimization is performed by calculating the CCSD energy on an initially coarse grid composed of different minimum exponents and even-tempering coefficients. Once regions of low CCSD energy are identified, a finer grid is used to obtain the final minimum exponent and even-tempering coefficient. In addition to the assumption of weak coupling, two other properties of the problem make this global optimization possible with modest computer resources; low dimensionality of search space and efficiency of atomic CCSD calculations.

Next, an atomic MCSCF calculation in a complete active space (CAS) with the optimized uncontracted basis is performed in GAMESS. For these calculations, all electrons not removed by the pseudopotential are allowed to excite. For helium, the active space consists of the orbitals from the n = 1 and n = 2 shells. For beryllium through neon, the active space includes the orbitals from the n = 2 and n = 3 shells. For magnesium through argon, the active space is composed of the orbitals from the n = 3 and n = 4 shells, with the exception of the 4D and 4F orbitals. A subset of the natural orbitals from the MCSCF calculations are used as the contracted functions of our basis.

All atomic calculations are performed in D 2h symmetry since GAMESS does not permit imposition of full rotational symmetry. Hence, different components of the same atomic subshell are not necessarily equivalent. Additionally, mixing may occur among orbitals of different angular momenta. For instance, there is mixing of S orbitals with both D 3z 2 -r 2 and D x 2 -y 2 orbitals. This anisotropy can be removed by averaging the different components of a particular subshell and zeroing out the off-diagonal blocks of the one-particle density matrix [7].

A simpler approach taken in this work is found to produce results of similar quality. For each angular momentum for which a contraction is desired, the NO with that angular momentum which has the largest occupation number is chosen. Additionally, NO elements which do not correspond to the dominant character of the orbital are zeroed out. For instance, an NO with large coefficients on the S basis functions and small coefficients on the D basis functions is considered to be dominated by S character, so the D coefficients are zeroed out. Finally, the NOs are normalized. The NOs selected in this procedure generate the contracted functions for the basis set. The expansions of the contractions are given in the supplementary material [12].

B. Gauss-Slater Primitives

GS functions [5] are defined as

ϕ ζ nlm (r, θ, φ) = N ζ n r n-1 e -(ζr) 2 1+ζr Z m l (θ, φ), (2)
where ζ is the GS exponent and N ζ n is the normalization factor. The restriction n ≥ l + 1 is imposed for GS functions. For r ≪ 1, the GS behaves like a Gaussian:

ϕ ζ nlm (r, θ, φ) ∼ = N ζ n r n-1 e -(ζr) 2 Z m l (θ, φ), (3)
and for r ≫ 1, the GS behaves like a Slater:

ϕ ζ nlm (r, θ, φ) ∼ = N ζ n r n-1 e -ζr Z m l (θ, φ). (4)
Consequently, GS functions introduce no cusp at the origin and can reproduce correct long-range asymptotic behavior of the orbitals. Unlike Gaussians and Slaters, normalization of GSs has no closed form expression. Nevertheless, normalizing an arbitrary GS is trivial with the following scaling relation between N ζ n and N 1 n :

N ζ n = ζ n+1/2 N 1 n . (5)
Values for N 1 n are given in the supplementary material [12].

Since GSs are not analytically integrable, the radial part must be expanded in Gaussians for use in QC programs that evaluate matrix elements analytically. The expansion is

ϕ ζ nlm (r, θ, φ) = i c ζ i 2(2α ζ i ) l+ 3 2 Γ(l + 3 2 ) r l e -α ζ i r 2 Z m l (θ, φ), (6)
where c ζ i is the i th expansion coefficient and α ζ i is the i th Gaussian exponent. Notice that the expansion permits the case for which n = l + 1 for the GS function. Additionally, the following scaling relations hold for the expansion coefficients and Gaussian exponents:

α ζ i = ζ 2 α 1 i (7) c ζ i = c 1 i . (8)
Once the Gaussian expansions are found for unit exponents, expansions of arbitrary GSs follow immediately from the scaling relations. For QC calculations in this paper, GSs are expanded in six Gaussians. However, if the purpose of the initial QC calculation is to generate crude starting orbitals for QMC calculations in which orbital optimization is performed, it is only necessary to expand GS primitives in a single Gaussian. In this case, the cost of QC calculations is the same for Gaussian and GS primitives. The expansions of GS functions with unit exponent in both one and six Gaussians are given in the supplementary material [12]. As mentioned above, the restriction n ≥ l + 1 is imposed for GS functions, instead of the more familiar n = l + 1 restriction imposed for Gaussian primitives. This motivates construction of two types of bases. In the first, ANO-GS, the restriction n = l+1 is enforced. In the second, ANO-GSn, for each l there can be at most a single GS primitive with a particular n. For each additional primitive with a particular l, n must be incremented.

For example, consider lithium. The 2z ANO-GS basis has one S contraction, one GS-1S function, two GS-2P functions, and one GS-3D function. On the other hand, the 2z ANO-GSn basis has one S contraction, one GS-1S function, one GS-2P function, one GS-3P function, and one GS-3D function.

A caveat to the above definition of the ANO-GSn basis is that GS-2S functions are not permitted since a single GS-2S function will introduce an undesired cusp in the wavefunction. Additionally, the 2z ANO-GS and ANO-GSn basis sets are identical for all elements except lithium and sodium. When the 2z ANO-GS and ANO-GSn basis sets are identical, the basis sets are referred to as a 2z ANO-GS/GSn basis. For both lithium and sodium, the basis sets differ because these systems have no P contractions and instead have a second P primitive for the 2z basis. This primitive is a GS-2P for the ANO-GS basis and a GS-3P for the ANO-GSn basis. Additionally, weak coupling between functions of different angular momentum causes the GS-1S and GS-3D functions in the ANO-GS bases for lithium and sodium to differ from their counterparts in the ANO-GSn bases. However, the optimal exponents differ by less than 0.01.

Optimal exponent selection for the GS primitives is discussed now. Instead of optimizing exponents for the atom as was done to generate the contractions, optimization of the GS exponents is performed for the homonuclear diatomic molecule at experimental bond length [START_REF] Huber | Constants of Diatomic Molecules[END_REF][START_REF] Gurvich | Thermodynamic Properties of Individual Substances[END_REF][START_REF]NIST Standard Reference Database Number 69[END_REF][START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF][START_REF] Grisenti | [END_REF][19][20][21][22]. This advantageously produces a balanced basis set.

Weak coupling between GS functions of different angular momenta is assumed, so the initial optimization for each angular momentum is performed separately. This assumption is validated in Figure 1, which contains plots of the CCSD energy for Si 2 while varying individual GS exponents in the 2z ANO-GS/GSn basis. Both the curve shape and exponent value which minimizes the energy vary little with fixed exponent value, signifying weak coupling between GS functions of different angular momentum. The optimization is performed at the CCSD level of theory using a Python wrapper around GAMESS. For each angular momentum, an energy landscape is defined by a grid of primitive exponents ranging from 0.1 to 6.0 with 0.1 spacing. Thorough investigation has revealed that exponents larger than 6.0 are not optimal for the systems considered. Low lying minima of this energy landscape are then handled with increasingly finer grids until energy changes are less than 0.01 mH. During this investigation of local minima, all angular momenta are handled simultaneously to account for any coupling effects. Results of this optimization are shown in Figure 2. Optimal exponents for ANO-GS and ANO-GSn bases exhibit a linear trend across each row of the periodic table. For nearly degenerate minima, the exponent following the trend in the figure is chosen as optimal, resulting in energy increase no greater than several 0.1 mH. The optimal GS exponents are given in the supplementary material [12].

In some cases, the optimal exponents for primitives with the same n and l are very close. This can lead to large equal and opposite coefficients on these basis functions when constructing molecular orbitals. Numerical problems could result, providing further motivation for the ANO-GSn basis, in which each pair of n and l is unique. However, all of our tests with the ANO-GS basis have had no numerical problems.

Finally, the optimal primitive exponents are found to depend weakly on the electronic structure method employed in the optimization, as demonstrated in Figure 3 for Si 2 with the 2z ANO-GS/GSn basis. The globally minimizing exponents are nearly equal in different methods. This exponent transferability to different levels of theory is extremely attractive for a basis set.

III. RESULTS

Section II demonstrates that the ANO-GS and ANO-GSn bases exhibit desirable properties. However, it remains to be shown that these basis sets produce accurate results. Fortunately, the basis set accompanying the BFD pseudopotential serves as a metric for testing ANO-GS and ANO-GSn basis quality. The BFD basis for elements in Groups 1A and 2A of the periodic table has recently been updated [23], but the number of functions in the new basis is inconsistent with the correlation consistent polarized basis prescription [1]. Since comparison would be difficult, their published functions are considered in this work.

Figure 4 shows the CCSD total energy gain per electron of the ANO-GS and ANO-GSn bases over the BFD bases [10] for atoms and homonuclear dimers of hydrogen through argon. Energy gains per electron tend to increase across each row of the periodic table. Both ANO-GS and ANO-GSn bases yield energy gains for most molecules and atoms. The energy gains per electron are generally larger for molecules than for atoms, and larger for the ANO-GSn basis than for the ANO-GS basis. The energy gains for the 2z bases are generally larger than for the 3z bases, as expected, since the energy left to recover becomes smaller as the basis size increases.

The ANO-GS and ANO-GSn bases also produce more accurate CCSD atomization energies than the BFD basis for the homonuclear dimers of hydrogen through argon. The 2z ANO-GS and ANO-GSn bases are identical for all elements except lithium and sodium. The GS-1S and GS-3D exponents for these elements each differ by less than 0.01 between 2z ANO-GS and ANO-GSn bases, so 2z ANO-GS and ANO-GSn are shown together as 2z ANO-GS/GSn. Exponents for GS functions of P angular momentum are not included for lithium and sodium since these elements have an extra primitive of P angular momentum.

GS-3D Exponent

FIG. 3. Change in Si2 energy for 2z ANO-GS/GSn basis shows optimal exponents depend weakly on electronic structure method (CCSD, HF, and B3LYP). Top: GS-1S exponent is varied with GS-2P and GS-3D exponents fixed at their optimal values. Middle: GS-2P exponent is varied with GS-1S and GS-3D exponents fixed at their optimal values. The large increase in energy around an exponent of 1.0 occurs since the P primitive and P contraction become nearly linearly dependent. Bottom: GS-3D exponent is varied with GS-1S and GS-2P exponents fixed at their optimal values. For Middle and Bottom, HF and B3LYP energy scale is on the right y-axis. This difference in energy scale occurs since higher angular momentum functions are less important in these effectively single-determinant theories. Atom 3z ANO-GS Atom 3z ANO-GSn Dimer 3z ANO-GS Dimer 3z ANO-GSn FIG. 4. CCSD total energy gains per electron of ANO-GS and ANO-GSn relative to the corresponding BFD basis [10] for atoms and homonuclear dimers of hydrogen through argon. Energy gains per electron tend to increase across each row of the periodic table. The 2z ANO-GS and ANO-GSn bases are identical for all elements except lithium and sodium. Differences between 2z ANO-GS and ANO-GSn results for these elements is ∼ 0.01mH, so they are shown together as 2z ANO-GS/GSn.

Figure 5 shows the fraction of experimental atomization energy recovered in CCSD for the homonuclear dimers which are not weakly bound. The 2z ANO-GS/ANO-GSn basis recovers more atomization energy than the 2z BFD basis for all dimers except those of Group 1A elements. Similarly, the 3z ANO-GSn basis recovers more atomization energy than the 3z BFD basis for the same systems, but the differences are small. The 3z ANO-GSn is on average slightly better than the 3z ANO-GS basis, the largest gains being for F 2 and Cl 2 .

For Group 1A elements, the BFD bases recover more atomization energy in CCSD than do their ANO-GS or ANO-GSn counterparts. This occurs due to inaccurate BFD energies for the atoms, as can be seen in Figure 4. However, as described above, we used the published BFD bases for these elements rather than the updated BFD bases [23] to maintain consistency.

Finally, improvements of the ANO-GS and ANO-GSn bases extend to other systems and methods. Figure 6 shows the fraction of experimental atomization energy recovered for five systems in the G2 set [START_REF] Curtiss | [END_REF] with the BFD, ANO-GS, and ANO-GSn bases in three quantum chemistry methods. For CCSD, the ANO-GS and ANO-GSn bases outperform the BFD basis for all systems. For sulfur dioxide the improvement due to the ANO-GS and ANO-GSn bases is dramatic: the 2z ANO-GS/GSn result is nearly halfway between the 2z and 3z BFD results, and the 3z ANO-GS/GSn result is nearly halfway between the 3z and 5z BFD results. ANO-GS and ANO-GSn benefits are more prominent in HF and B3LYP: for most systems, the 2z ANO-GS/GSn result is closer to the 3z BFD result than the 2z BFD result, and the 3z ANO-GS/GSn result is closer to the 5z BFD result than the 3z BFD result. Differences between results with the ANO-GS and ANO-GSn bases are small. Figure 7 shows the fraction of experimental atomization energy recovered using diffusion Monte Carlo (DMC) with the BFD, ANO-GS, and ANO-GSn bases. For each system, the DMC calculations are performed with both a single-configuration state function (single-CSF) reference (DMC-1CSF) and full-valence complete active space reference (DMC-FVCAS). However, for each of the constituent atoms in these molecules, the FVCAS and single-CSF references are equivalent. All DMC calculations are performed with a 0.01 H -1 time step and trial wavefunction obtained by optimizing Jastrow, orbital, and configuration state function (CSF) parameters (where applicable) via the linear method [28][29][30] in variational Monte Carlo. The DMC-1CSF and DMC-FVCAS calculations exhibit similar trends to the HF and B3LYP calculation for most systems: the 2z ANO-GS/GSn result is closer to the 3z BFD result than the 2z BFD result, and the 3z ANO-GS/GSn result is closer to the 5z BFD result than the 3z BFD result. Again, differences between results with the ANO-GS and ANO-GSn bases are small. There are several important points that can be made by comparing the DMC calculations of Figure 7 to the CCSD calculations of Figure 6. First, the DMC results for the atomization energies have a weaker dependence on basis size than the CCSD results. Second, for a given basis set, the most basic DMC calculations, DMC-1CSF, yield superior results compared to CCSD. In addition to yielding superior results, DMC-1CSF calculations have better computational cost scaling than CCSD calculations. Under certain assumptions, the cost of DMC-1CSF calculations scales as O(N 3 ) [31], while the cost of CCSD calculations scales as O(N 6 ) [START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF], where N is the number of electrons. However, it is important to note that the prefactor of the scaling is significantly smaller for the CCSD calculations.

Finally, our results are not the first to show that DMC calculations can produce accurate atomization energies. In particular, DMC-1CSF calculations of the entire G2 set have been performed for both pseudopotential and all-electron systems [START_REF] Grossman | [END_REF]34] and produced excellent results. Additionally, there is good agreement between the pseudopotential and all-electron results with a mean absolute deviation of about 2.0 kcal/mol over the entire G2 set [34]. Although these previous results are very good, there is room for improvement, particularly for the open shell systems. A systematic study with DMC-FVCAS calculations is currently underway in our group, which should produce results to (near) chemical accuracy for all systems in the G2 set. [START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF]24] to compare with experiment [START_REF] Huber | Constants of Diatomic Molecules[END_REF][START_REF]NIST Standard Reference Database Number 69[END_REF][START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF][START_REF] Luo | Comprehensive Handbook of Chemical Bond Energies[END_REF].

IV. CONCLUSION

A simple yet general method for constructing basis sets for molecular electronic structure theory calculations has been presented. These basis sets consist of a combination of atomic natural orbitals from an MCSCF calculation with primitive functions optimized for the corresponding homonuclear dimer. The functional form of the primitive functions is chosen to have the correct asymptotics for the nuclear potential of the system.

It was shown that optimal exponents of primitives with different angular momenta are weakly coupled. This enables efficient determination of optimal exponents. Additionally, it was demonstrated that the particular electronic structure method employed in optimization has little effect on the optimal values of the primitive exponents.

Two sets of 2z and 3z bases, ANO-GS and ANO-GSn, appropriate for the Burkatzki, Filippi, and Dolg nondivergent pseudopotentials were constructed for elements hydrogen through argon. Since these pseudopotentials do not diverge at nuclei and have a Coulomb tail, GS functions are the appropriate primitives.

It was demonstrated that both ANO-GS and ANO-GSn basis sets offer significant gains over the Burkatzki, Filippi and Dolg basis sets for CCSD, HF, B3LYP [11], and QMC calculations. Improvements were observed in both total energies and atomization energies. The latter indicates that basis sets providing a balanced description of atoms and molecules were produced by using both the atom and the dimer in the optimization. On average, the ANO-GSn basis is slightly better than the ANO-GS basis, but either is a sound choice.

In the future, these basis sets will be extended to include the transition metals, and, bases will be constructed for all-electron calculations, for which Slater functions are the appropriate primitives. with BFD, ANO-GS, and ANO-GSn bases. The 2z ANO-GS and ANO-GSn bases yield different results only for LiF. The 5z BFD* calculations do not include the G or H functions from the 5z BFD basis. Calculated atomization energies are corrected for zero point energy [START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF]24] to compare with experiment [START_REF] Huber | Constants of Diatomic Molecules[END_REF][START_REF]NIST Standard Reference Database Number 69[END_REF][START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF]27] 7. Fraction of experimental atomization energy recovered in diffusion Monte Carlo (DMC) for LiF, O2, P2, S2, and SO2 with the BFD, ANO-GS, and ANO-GSn bases. DMC calculations are performed with both a single-CSF reference (DMC-1CSF) and full-valence complete active space reference (DMC-FVCAS). The 2z ANO-GS and ANO-GSn bases yield different results only for LiF. The 5z BFD* calculations do not include the G or H functions from the 5z BFD basis. Calculated atomization energies are corrected for zero point energy [START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF]24] to compare with experiment [START_REF] Huber | Constants of Diatomic Molecules[END_REF][START_REF]NIST Standard Reference Database Number 69[END_REF][START_REF]NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101[END_REF]27]. The legend for this plot is identical to that of Figure 6.

8 FIG. 1 .

 81 FIG.1. Change in Si2 CCSD energy for 2z ANO-GS/GSn basis shows weak coupling between GS functions of different angular momenta. TOP: Energy versus GS-1S exponent for three values of the GS-2P exponent with the GS-3D exponent fixed at its optimal value. Bottom: Energy versus GS-3D exponent for three values of the GS-1S exponent with the GS-2P exponent fixed at its optimal value.

  Electron (mH) Atom 2z ANO-GS/GSn Dimer 2z ANO-GS/GSn

  Fraction of experimental atomization energy recovered in CCSD with BFD, ANO-GS, and ANO-GSn bases for the homonuclear dimers of hydrogen through argon which are not weakly bound. The 2z ANO-GS and ANO-GSn bases are identical for all elements except lithium and sodium. Differences between 2z ANO-GS and ANO-GSn atomization energies for these elements is ∼ 0.01mH, so they are shown together as 2z ANO-GS/GSn. Calculated values are corrected for zero point energy

	Fraction of Exp. Binding Energy	1.2	2z BFD 2z ANO-GS/GSn 3z BFD	3z ANO-GS 3z ANO-GSn 5z BFD
			H 2	Li 2	B 2	C 2	N 2	O 2	F 2	Na 2	Al 2	Si 2	P 2	S 2	Cl 2
	FIG. 5.

  Fraction of experimental atomization energy recovered in HF, B3LYP, and CCSD for LiF, O2, P2, S2, and SO2

	Fraction of Exp. Binding Energy	0.70 0.80 0.90 1.00 1.10	CCSD
		0.60	LiF	O 2	P 2	S 2	SO 2
	Fraction of Exp. Binding Energy	0.70 0.80 0.90 1.00 1.10	B3LYP
		0.60	LiF	O 2	P 2	S 2	SO 2
	Fraction of Exp. Binding Energy	0.20 0.40 0.60 0.80 1.00	HF			2z BFD 2z ANO-GS 2z ANO-GSn 3z BFD 3z ANO-GS 3z ANO-GSn 5z BFD*
		0.00	LiF	O 2	P 2	S 2	SO 2
		6.

  .

	Fraction of Exp. Binding Energy	0.70 0.80 0.90 1.00 1.10	DMC-FVCAS
		0.60	LiF	O 2	P 2	S 2	SO 2
	Fraction of Exp. Binding Energy	0.70 0.80 0.90 1.00 1.10	DMC-1CSF
		0.60	LiF	O 2	P 2	S 2	SO 2

FIG.

V. ACKNOWLEDGMENTS

We thank Claudia Filippi for very valuable discussions. This work was supported by the NSF (Grant Nos. DMR-0908653 and CHE-1004603). Computations were performed in part at the Computation Center for Nanotechnology Innovation at Rensselaer Polytechnic Institute.