In this paper, we propose a construction of a new BGK model generalizing the Ellipsoidal Statistical Model ([2], [19]) to the context of gas mixtures. The derivation of the model is based on the introduction of relaxation coefficients associated to some moments and the resolution of a minimization problem as in ([7], [8], [10]). We obtain in this work, an ESBGK model for gas mixtures satisfying the fundamental properties of the Boltzmann collision operator (conservation laws, H theorem, equilibrium states,. . .) and that is able to give a range of Prandtl numbers including the indifferentiability situation. subject classifications 35Q20, 35Q35 1. Introduction The complexity of the nonlinear Boltzmann operator suggests to introduce simpler kinetic models. Hence, the BGK model ([17]) which is a well known simplified kinetic relaxation model has been introduced. Its features are to replace the complicated integral of collisions by a relaxation model while keeping some important physical and mathematical properties of the interaction term (con-servation laws, H theorem, equilibrium states,. . .). The interest of these models is that they are easier to handle numerically and less costly at a computational point of view. However the extension of the BGK model to multi-component gases meets fundamental difficulties. For example, the hydrodynamic limit is much more complicated since phenomena such as diffusion or thermal diffusion must be considered. Many BGK models have been proposed up to now, but they are still unsatisfactory and have been constructed on the same manner. That is they must reproduce the rates of exchanges of impulsion and energy of the Boltzmann operator for Maxwell molecules between the different species. For example one polynomial model has been proposed in [13]. But its relaxation function being a polynom, the nonnegativity of the distribution function is not satisfied. However this model satisfies the indiffer-entiability principle stated in [13] inherited from the Boltzmann operator. That is when all the molecules have the same masses and their cross sections are equal, then the system of equations reduces to a single one by adding the distribution functions. A BGK model enjoying good mathematical properties has been derived in [1]. The main idea is to introduce only one BGK operator per species whose macrosopic parameters reproduce the interaction between each species with the others. But this model leads to uncorrect transport coefficients at the hydrodynamic limit ([21]). In [21], the authors constructed a BGK model which coincides with the Grad moments of the linearized Boltzmann operator. However the relaxation function of the model being polynomial, the nonnegativity of the distribution function is not sure. Hence the aim of this paper is to construct a BGK model that is able to recover correct transport coefficients and enjoying good mathematical properties. The question is important because in general the authors compute the hydrodynamic limit from their own model and eventually compare with the correct fluid model ([18]). In the situation of single and monospecies a BGK model leading to the correct Prandtl number has been proposed in ([19]). This model is called Ellipsoidal Statistical Model (ESBGK). But the proof of the H theorem has been shown later on ([2]). Recently,Read less <