Consistent lattice Boltzmann model for multicomponent mixtures

A new lattice Boltzmann model for multicomponent ideal gas mixtures is presented. The model development consists of two parts. First, a new kinetic model for Stefan- Maxwell diffusion amongst the species is proposed and realized as a lattice Boltzmann equation on the standard discrete velocity set. Second, a compressible lattice Boltzmann model for the momentum and energy of the mixture is established. Both parts are consistently coupled through mixture composition, momentum, pressure, energy and enthalpy whereby a passive scalar advection-diffusion coupling is obviated, unlike in previous approaches. The proposed model is realized on the standard three-dimensional lattices and is validated with a set of benchmarks highlighting various physical aspects of compressible mixtures. Stefan-Maxwell diffusion is tested against experiment and theory of uphill diffusion of argon and methane in a ternary mixture with hydrogen. The speed of sound is measured in various binary and ternary compositions. We further validate the Stefan-Maxwell diffusion coupling with hydrodynamics by simulating diffusion in opposed jets and the three-dimensional Kelvin-Helmholtz instability of shear layers in a two-component mixture. Apart from the multicomponent compressible mixture, the proposed lattice Boltzmann model also provides an extension of the lattice Boltzmann equation to the compressible flow regime on the standard three-dimensional lattice.

A new lattice Boltzmann model for multicomponent ideal gas mixtures is presented. The model development consists of two parts. First, a new kinetic model for Stefan-Maxwell diffusion amongst the species is proposed and realized as a lattice Boltzmann equation on the standard discrete velocity set. Second, a compressible lattice Boltzmann model for the momentum and energy of the mixture is established. Both parts are consistently coupled through mixture composition, momentum, pressure, energy and enthalpy whereby a passive scalar advection-diffusion coupling is obviated, unlike in previous approaches. The proposed model is realized on the standard three-dimensional lattices and is validated with a set of benchmarks highlighting various physical aspects of compressible mixtures. Stefan-Maxwell diffusion is tested against experiment and theory of uphill diffusion of argon and methane in a ternary mixture with hydrogen. The speed of sound is measured in various binary and ternary compositions. We further validate the Stefan-Maxwell diffusion coupling with hydrodynamics by simulating diffusion in opposed jets and the three-dimensional Kelvin-Helmholtz instability of shear layers in a two-component mixture. Apart from the multicomponent compressible mixture, the proposed lattice Boltzmann model also provides an extension of the lattice Boltzmann equation to the compressible flow regime on the standard three-dimensional lattice.

Introduction
The lattice Boltzmann method (LBM) is a recast of fluid dynamics into a fully discrete kinetic system of designer particles with the discrete velocities c i , i = 0, . . . , Q − 1, fitting into a regular space-filling lattice, with the kinetic equation for the populations f i (x, t) following a simple algorithm of "stream along links c i and collide at the nodes x in discrete time t". Since its inception Higuera & Jiménez 1989), LBM has evolved into a versatile tool for the simulation of complex flows including transitional flows (Dorschner et al. 2017), flows in complex moving geometries (Dorschner et al. 2016), compressible flows (Frapolli et al. 2016b;Dorschner et al. 2018), multiphase flows (Mazloomi et al. 2015(Mazloomi et al. , 2017Wöhrwag et al. 2018) and rarefied gas (Shan et al. 2006), to mention a few recent instances; see Succi (2018); Krüger et al. (2017) for a discussion of LBM and its application areas.
In view of extensive development, it seems surprising that the multicomponent gas mixtures so far resisted a significant advancement in the LBM context. Importance of compressible mixtures is hard to overestimate because they are a prerequisite for combustion applications (Williams 1985). However, incorporating even the basic mechanism of multicomponent diffusion in gases, the Stefan-Maxwell diffusion, remains an essentially unsolved problem in the LBM context, in spite of that the Stefan-Maxwell diffusion is itself a derivative of Boltzmann's kinetic theory (Chapman & Cowling 1990). It is worth reminding that the Stefan-Maxwell diffusion mechanism is well recognized as a fundamental feature of gas mixtures, supported by experiment (Toor 1957;Duncan & Toor 1962;Arnold & Toor 1967) and molecular dynamics simulations (Wheeler & Newman 2004;Krishna & van Baten 2005). As highlighted by Krishna & Wesselingh (1997), the Stefan-Maxwell diffusion is more subtle than the conventional Fick's model. The latter implies that any component in a mixture moves from higher to lower concentration regions. The Stefan-Maxwell model, on the other hand, accounts for binary interaction between each of the species pairs through pairwise diffusion coefficients and can lead to counter-intuitive effects such as uphill diffusion when a component in a ternary mixture moves from lower to the higher concentration region (Toor 1957;Duncan & Toor 1962;Arnold & Toor 1967). Among applications of Stefan-Maxwell diffusion in the conventional computational fluid dynamics, we mention recent studies of diffusion in fuel cells (Hsing & Futerko 2000;Stockie et al. 2003;Suwanwarangkul et al. 2003).
However, the majority of existing LBM models for the multicomponent mixtures Chiavazzo et al. 2009;Feng et al. 2019;Huang et al. 2019;Hosseini et al. 2018;Lin et al. 2017) are bound to use the Fick diffusion model rather than the Stefan-Maxwell. To the best of our knowledge, the only LBM realization of the Stefan-Maxwell diffusion was reported in a very recent work of Chai et al. (2019); However, the two-dimensional LBM of Chai et al. (2019) is restricted by the isothermal and isobaric assumptions and can thus not provide a basis for the development of a compressible mixture LBM. Another obstacle arises at the coupling of diffusion to the transport of momentum and energy. The simplest way of tackling multicomponent mixtures with the LBM is by representing the dynamics of the species by a advectiondiffusion equation (see, e. g. Chiavazzo et al. (2009) and references therein). In this approach, the species are treated as passive scalars, advected with the fluid velocity and the species do not influence the fluid or other species. The passive scalar viewpoint on LBM for mixtures was adopted and extended in a number of recent proposals (Feng et al. 2019;Huang et al. 2019;Hosseini et al. 2018;Lin et al. 2017). However, apart from the inability of incorporating the Stefan-Maxwell diffusion, the passive scalar approach suffers from a more fundamental shortcoming, namely thermodynamic inconsistency. For example, the aforementioned models do not readily recover the correct heat flux in the multicomponent system (Chapman & Cowling 1990;Williams 1985;Bird et al. 2006) and miss the enthalpy flux due to diffusion.
In this paper, we revisit the LBM construction for a compressible multicomponent mixture, focusing on a thermodynamically consistent coupling between the Stefan-Maxwell diffusion and momentum and energy transfer in the system. We begin in Sec. 2 with setting up a kinetic system for the species in the M -component mixture. The construction follows the path of so-called quasi-equilibrium relaxation models (Gorban & Karlin 1994;Ansumali et al. 2007); see Arcidiacono et al. (2007) in the context of isothermal mixtures. Here, we significantly extend the quasi-equilibrium kinetic model for the species to a generic ideal gas equation of state and, unlike in the earlier approach of Arcidiacono et al. (2007), enabling the Stefan-Maxwell constitutive relation. After a short summary of nomenclature in Sec. 2.1, the species kinetic equations are introduced in Sec. 2.2, in the continuous time-space setting. We show in Sec. 2.3 that the proposed kinetic equations recover the Stefan-Maxwell diffusion together with the barodiffusion in the hydrodynamic limit. The species kinetic equations are realized on the standard set of discrete velocities in Sec. 2.4. In Sec. 2.5, we derive the lattice Boltzmann scheme for the species kinetic equations following the technique of integration along characteristics introduced by He et al. (1998). This concludes the first part of the model development.
We continue in Sec. 3 with a mean-field lattice Boltzmann formulation of the mixture momentum and energy. After a summary on the mixture energy and enthalpy in Sec. 3.1, we present a generic two-population lattice Boltzmann equation for the mixture. We note that the mean-field approach requires only two lattice Boltzmann equations, one for the mixture density and momentum and another one for the energy. While the two-population LBM is established approach for a single-component compressible fluid (Frapolli et al. 2015;Saadat et al. 2019), the application of the two-population techniques to the mixture requires a modification of the nonequilibrium fluxes discussed in Sec. 3.2. The mixture density, momentum and energy equations are presented in Sec. 3.3 while details of their derivation with the Chapman-Enskog analysis (Chapman & Cowling 1990) are summarized in the Appendix A. The two-population mixture LBM is realized on the standard lattice in Sec. 3.4 where we extend the two-dimensional compressible LBM of Saadat et al. (2019) to three-dimensional mixtures. Finally, in Sec. 3.5 we discuss the coupling between the M LBM equations for the species and the double-population mean field mixture LBM. The resulting LBM provides a reduced description of the Mcomponent mixture with M + 2 tightly coupled lattice Boltzmann equations, unlike a standard kinetic approach which would require 2 × M kinetic equations.
In Sec. 4, the LBM model is validated with a number of select benchmarks. After a summary of general aspects of numerical implementation in Sec. 4.1, we present a simulation of diffusion of argon and methane ternary mixture with hydrogen in the Loschmidt tube apparatus in Sec. 4.2, along with the classical experiment of Arnold & Toor (1967) and theoretical discussion of Krishna & Wesselingh (1997). We show that the LBM simulations reproduce in a quantitative fashion the experimentally observed features of the Stefan-Maxwell diffusion such as uphill and osmotic diffusion and the diffusion barrier (Toor 1957;Duncan & Toor 1962;Arnold & Toor 1967;Krishna & Wesselingh 1997). The coupling between hydrodynamics and diffusion is validated in a counterflow diffusion in opposed jets in Sec. 4.3 and the speed of sound measurements are presented in Sec. 4.4 for probing the compressible flow aspect of the model. Finally, a simulation of the three-dimensional Kelvin-Helmholtz instablity in a binary mixture is reported in Sec. 4.5 as a test for the performance of the proposed LBM in a complex flow. Conclusions are drawn in Sec. 5.

Composition and equation of state of ideal gas mixture
We begin with introducing some nomenclature and notation. Let us consider a mixture composed of M ideal gases. The composition is described by the species densities ρ a , a = 1, . . . , M , while the mixture density is (2.1) Equivalently, the mixture composition is defined by the mixture density ρ and the M − 1 independent mass fractions Y a , With the molar mass of the component m a , the mean molar mass m depends on the composition, The equation of state (EoS) provides a relation between the pressure P , the temperature T and the composition, Here, R is the specific gas constant that contains the information about the composition of the gas by way of the mean molar mass m, where R U ≈ 8.314 kJ/K · kmol is the universal gas constant. Thus, for a mixture of ideal gases, the specific gas constant R is a function of local composition and changes in space and time. The pressure of an individual component P a is related to the pressure of the mixture P through Dalton's law of partial pressures as follows: where the mole fraction of a component X a is related to its mass fraction Y a as A consequence of Dalton's law of partial pressure is that M a=1 P a = P . Combined with the equation of state (2.4), the partial pressure P a takes the form, where R a is the specific gas constant of the component, (2.9) With these thermodynamic relations in mind, we proceed to setting up the kinetic equations that recover the Stefan-Maxwell diffusion in the macroscopic limit.

Kinetic equation for the species
In this section, we set up kinetic equations which recover the Stefan-Maxwell diffusion model for a M -component ideal gas mixture. Each component is described by a set of populations f ai , a = 1, . . . , M , corresponding to the discrete velocities c i , i = 0, . . . , Q−1. The proposed kinetic equation for each of the species a is written as, (2.10) Here, ρ a is the density of the component a, which is defined as the zeroth moment of populations f ai , Furthermore, a symmetric set of relaxation parameters θ ab = θ ba shall be related to the binary diffusion coefficients below. The equilibrium f eq ai and the quasi-equilibrium f * ai populations will be fully defined in section 2.4. Here, we only need to specify the conditions for the low-order moments thereof. To that end, let us introduce the partial momenta ρ a u a as first moments of the species' populations, (2.12) The quasi-equilibrium populations f * ai must satisfy the following set of constraints, (2.14) The momenta of the components sum up to the mixture momentum, M a=1 ρ a u a = ρu. (2.15) The equilibrium populations f eq ai have to verify the following set of constraints: In Eq. (2.18), the partial pressure P a (2.8) depends on the temperature T , which is obtained from the mixture kinetic equations of Sec. 3. Finally, the quasi-equilibrium distribution must match the equilibrium if the species velocity equals the velocity of the mixture, u a = u: Some comments are in order: (i) The M -component kinetic system satisfies M + D conservation laws, where D is the space dimension: The densities ρ 1 , . . . , ρ M and the vector of fluid momentum ρu are locally conserved fields.
(ii) Thanks to the matching condition (2.19), the relaxation term on the right hand side of Eq. (2.10) vanishes only at the equilibrium. We now proceed with the identification of the relaxation parameters θ ab in terms of the Stefan-Maxwell binary diffusion coefficients.

Hydrodynamic limit of kinetic equations for the species
Evaluation of the zeroth and of the first moments of the kinetic equation (2.10) results in the balance equations for the species densities and species velocities, (2.20) (2.21) Here, P a is the partial pressure tensor, Upon summation over the components in (2.20), we arrive at the continuity equation for the mixture density, while the summation over components in (2.21) results in the mixture momentum balance, where P is the mixture pressure tensor, (2.25) The low-order closure relation for the species balance equation (2.20) is established by considering a perturbation around the equilibrium, where the perturbation δu a satisfies the consistency condition, To first order, upon substitution into (2.21), we get the constitutive relation for the diffusion velocity u a , (2.28) Upon summation over the species, and by taking into account Dalton's law in the equilibrium pressure tensor (2.18), the compressible Euler equation for the flow velocity is established, (2.29) By elimination of the time derivative in (2.28), we get the constitutive relation as follows: (2.30) The constitutive relation (2.30) becomes the Stefan-Maxwell diffusion equation once the relaxation parameters θ ab are identified in terms of the binary mass diffusion coefficients D ab as follows: (2.31) Summarizing, kinetic equations for the species (2.10) recover the Stefan-Maxwell law of diffusion in the hydrodynamic limit, with both the diffusion due to non-uniformity of the species concentration and the barodiffusion taken into account. The present model does not include thermodiffusion as it should be expected by the simplicity of the relaxation term. We comment that the above derivation of (2.30) assumes the validity of the equation of state. The latter, in turn, depends on the temperature derived from the mixture energy equation, and which shall be introduced in Sec. 3. We now proceed with finalizing the continuous time-space kinetic equations by identifying the equilibrium and the quasiequilibrium populations.

Realization on the standard lattice
The above kinetic model is realized on the standard three-dimensional D3Q27 lattice, where D = 3 stands for three dimensions and Q = 27 is the number of discrete velocities: (2.32) Following Karlin & Asinari (2010), we define a triplet of functions in two variables, ξ and ζ > 0, For a vector ξ = (ξ x , ξ y , ξ z ), we consider a product-form associated with the discrete velocities c i (2.32), (2.34) The equilibrium f eq ai and the quasi-equilibrium f * ai are represented with the product-form (2.34) by choosing ξ = u or ξ = u a , respectively, and by assigning ζ = R a T in both cases: (2.36) One can readily verify that, the equilibrium (2.35) and the quasi-equilibrium (2.36) satisfy all the constraints put forward in sec. 2.2. We now proceed with the lattice Boltzmann discretization of the kinetic equations (2.10).

Kinetic equations in the relaxation form
With the mass diffusivity (2.31), the kinetic equation (2.10) is written as It can readily be seen that, in its present form, Eq. (2.37) is not well suited for numerical implementation. Indeed, in an actual problem, the density of some species can be small or even vanishing if a particular gas component is absent at some location. This is an inconvenience rather than a failure since vanishing of the density is compensated by the simultaneously vanishing molar fraction in the product X a X b . Hence, we first transform equation (2.37) in order to eliminae this artifact. Substituting the equation of state (2.4) into Eq. (2.37), we get (2.38) Equation (2.38) is equivalent to equation (2.37) and does not suffer from a spurious division by a vanishing density. Furthermore, it proves convenient to recast Eq. (2.38) in a relaxation form. To that end, let us define characteristic times τ ab = τ ba , and let us introduce the relaxation times τ a , Finally, let us introduce a shorthand notation, With these definitions, the kinetic equation (2.38) can be rearranged as follows: The species kinetic equations (2.42) are now cast into the relaxation form, familiar from previous lattice Boltzmann models: The right hand side comprises the conventional relaxation term and a source term (2.41). The latter depends on the populations only through the local densities, momenta and the temperature, as prescribed by the local equilibrium and quasi-equilibrium populations (2.35) and (2.36). Hence, kinetic equation (2.42) is amenable to a lattice Boltzmann discretization in time and space.

Derivation of the lattice Boltzmann equation
Following a procedure first introduced by He et al. (1998), we integrate (2.42) along the characteristics and apply the trapezoidal rule on the right hand side to obtain, Next, we introduce transformed populations k ai , Let us evaluate the pertinent moments of the transform (2.44). Summation over the discrete velocities gives, where we have specified that the density ρ a (f ) on the left hand side is defined using the original populations f ai while the density ρ a (k) is defined by the zeroth moment of the k-populations, Thus, the species densities do not alter under the populations transformation, and the specification can be dropped: ρ a = ρ a (f ) = ρ a (k). On the other hand, evaluating the first moment of (2.44) gives, where the species velocity u a (k) is defined by the k-populations in a conventional way, Summation over the species in (2.47) shows that the momentum ρu is also an invariant of the transform (2.44): Since the term F ai vanishes at equilibrium, and also thanks to the invariance of the local conservation (2.45) and (2.49), the equilibrium is the fixed point of the map (2.44): (2.52) The last term in Eq. (2.52) is spelled out as follows: The quasi-equilibrium population f * ai in the expression F ai (2.41) depends on the species velocity u a (f ). The latter, unlike the mixture velocity u(f ) (2.49), is not invariant under the transform to the k-populations. Rather, u a (f ) has to be evaluated using the linear relation (2.47) in terms of u b (k) by solving a M × M linear system for each of the spatial components.

Summary of the lattice Boltzmann equation for the Stefan-Maxwell diffusion
For convenience, we summarize the lattice Boltzmann equation for the species. We return to a more conventional notation and rename k ai to f ai in (2.52), The last term is written where we have introduced the diffusion velocity of the components, V a , a = 1, . . . , M . The latter are determined by (2.47) which can be recast as follows: (2.55) The field ρ a u a in the right hand side of (2.55) is defined by the moment relation as before, The M +D independent conservation laws of the M -component system (2.53) correspond to the mass of each component and the momentum of the mixture. The corresponding locally conserved fields are the species densities ρ a and the momentum flux ρu, The conservation law of the fluid mass is the implication of the mass conservation of each component. The density of the mixture ρ is defined by the densities of the components, while the flow velocity u is, (2.60) Note that, the velocity of the component is defined as the sum of the diffusion velocity and the flow velocity, V a + u, in the quasi-equilibrium populations contributing in (2.54), rather than as the first moment of the populations (2.56). Finally, we remind that the temperature dependence in the equilibrium and the quasi-equilibrium populations has to be supplied by the energy equation to be discussed in the next section.
3. Lattice Boltzmann model of mixture momentum and energy 3.1. First law of thermodynamics for ideal gas mixture For further references and notation, we open this section with a summary of the first law of thermodynamics for ideal gas mixtures. The caloric equation of state of a singlecomponent ideal gas provides for the specific mole-based internal energy of species a: ( 3.1) Here,C a,v is the specific heat at constant volume. Thus, the specific enthalpy reads, whereC a,p is the specific heat at constant pressure defined by Mayer's relation, Proceeding from the mole-basis onto the mass-basis, the specific heats are defined relative to the molar mass, while the mass-based specific internal energy and enthalpy are, Finally, the Mayer relation in the mass-basis reads, where the gas constant R a is defined by (2.9). Switching to the case of a M -component mixture, the mixture internal energy ρU is defined on the mass-basis as follows: (3.9) The specific mixture internal energy U can be rewritten, where the specific heat at constant volume is the mass-averaged value over the composition, Similarly, the mixture enthalpy ρH is defined as, (3.12) while the specific mixture enthalpy H can be transformed in the manner of Eq. (3.10), The specific heat at constant pressure reads, while both the specific heats satisfy the Mayer relation, with the mixture gas constant R defined by Eq. (2.5). Below, we shall formulate the lattice Boltzmann equation for the mixture density, momentum and energy for a generic case of temperature-dependent specific heats of the components.

Two-population lattice Boltzmann equation for the mixture
Point of departure is a lattice Boltzmann model for a single-component ideal gas with variable Prandtl number and adiabatic exponent. To that end, several suitable singlecomponent lattice Boltzmann models exist in the literature; here we mention compressible LBM by Frapolli et al. (2015Frapolli et al. ( , 2016a and by Saadat et al. (2019). The common feature of these single-component models is the use of the double-population construction, the idea first introduced in the context of incompressible thermal convective LBM in the classical paper by He et al. (1998) and further expanded in Guo et al. (2007); Karlin et al. (2013);. One set of populations, commonly quoted as f -populations, represents the density and momentum as the locally conserved fields of the corresponding f -LBM equation while another set, the g-populations represent the energy as the local conservation of the g-LBM kinetics. The coupling between the f -and g-LBM equations is also well understood and enables the realization of an adjustable Prandtl number and adiabatic exponent. Various realizations differ by the choice of the discrete velocities of the f -and g-sets; in particular, the compressible LBM of Frapolli et al. (2015) employs higher-order lattices with higher isotropy while the two-dimensional model developed in Saadat et al. (2019) uses the standard lattice with correction terms to compensate for insufficient isotropy.
Whichever single-component double-population model is taken as the starting point for representing a multi-component mixture, the central question is how to modify it. Note that, this question would not arise if one would follow the conventional approach by extending the already available M species LBM equations of sec. 2 to represent the energy equation of the mixture. However, with the double-population approach, this would lead to 2 × M lattices since the lattice for each component would need to be doubled to represent the energy of that component. On the contrary, the mean-field approach pursuit here avoids the kinetic representation of partial energies, instead it addresses only the total energy of the mixture by a single g-set. This requires only M + 2 lattices, M for the species and two for the mixture momentum and energy.
Below, we refer to the f -populations as the momentum lattice, and the g-populations as the energy lattice. For the momentum lattice, the locally conserved fields are the density and the momentum of the mixture, For the energy lattice, the locally conserved field is the total energy of the mixture, Here, the total energy ρE is the sum of the mixture internal energy ρU (3.10) and the kinetic energy ρu 2 /2, Since the mixture internal energy (3.10) depends on the composition, it is the first instance where the species kinetic equations become coupled with the kinetic equations for the mixture. Conversely, the temperature of the mixture is computed from the integral equation, (3.20) In the simplest case, the specific heats of all components can be approximated by constants and the temperature becomes (3.21) The temperature evaluated from solving (3.20) is used as the input in the definition of the pressure P in the equilibrium and quasi-equilibrium constraints of the species lattice Boltzmann system. This furnishes the input from the energy lattice into the species lattices.
We comment that, in the present section, the mixture density (3.16) and momentum (3.17) are defined self-consistently in the sense of f -populations of the momentum lattice. On the other hand, quantities carrying the same physical meaning were independently and also self-consistently defined earlier using the species populations, Eqs. (2.59) and (2.58), respectively. Doubling of the conservation of the total mass and momentum is the feature of the intermediate steps of the construction during which the species subsystem and the mixture subsystem are set up independently from one another. At the end of the construction, the doubling of the conservation shall be resolved through a coupling of both the species and the mixture subsystems in Sec. 3.5.
The lattice Boltzmann equations for the momentum and for the energy lattice are patterned from the single-component developments, where relaxation parameters ω and ω 1 shall be related to the viscosity and thermal conductivity below. We now proceed with specifying the constraints on the equilibrium populations f eq i and g eq i , and the quasi-equilibrium g * i in order that the system (3.22) and (3.23) recovers the momentum and energy equations of the mixture.
First, the equilibrium populations must satisfy the D + 2 conservation laws, Second, the equilibrium pressure tensor P eq and the tensor of equilibrium third-order moments Q eq of the momentum lattice must verify the Maxwell-Boltzmann relations in order to recover the compressible flow momentum equation, where overline denotes symmetrization. Similarly, the equilibrium mixture energy flux q eq and the second-order moment tensor R eq pertinent to the energy lattice are, The mixture equation of state P (2.4), the mixture gas constant R (2.5) and the specific enthalpy of the mixture H (3.13) entering the constraints (3.27), (3.28), (3.29), (3.29) and (3.30) depend linearly on the composition through the local mass fractions Y a .
To that end, the constraints on the equilibrium populations of the mixture momentum and energy lattices is a straightforward extension of those of the single-component doublepopulation LBM for compressible flow, where the ideal gas equation of state, the internal energy and the enthalpy are merely replaced by their mixture-averaged counterparts. A major difference comes next with the constraints for the quasi-equilibrium. The zeroth, the first-and the second-order moments of the quasi-equilibrium populations g * i , or the quasi-equilibrium energy ρE * , the energy flux q * and the flux of the energy flux R * , respectively, have to satisfy the following relations: (3.33) The first and the third of these quasi-equilibrium constraints, Eqs. (3.31) and (3.33), as well as the first and the second terms in the quasi-equilibrium energy flux (3.32) are again the direct extension of the single-component LBM. Specifically, the two first terms in (3.32), comprising the energy flux q and the pressure tensor P , are needed to decouple the viscosity from thermal conductivity, and to maintain a variable Prandtl number, in both the single-and multcomponent cases. The remaining two terms in the quasi-equilibrium energy flux (3.32), q diff and q corr are specific to the multicomponent case and appear due to the mean-field approach to the energy representation. The interdiffusion energy flux q diff reads as follows: where the diffusion velocities V a are defined according to Eq. (2.55). The interdiffusion energy flux contributes the enthalpy transport due to diffusion and hence it vanishes in the single-component case. The effect of the inter-diffusion energy flux is typically significant at the initial stages of the diffusion process and cannot be neglected. Finally, the correction flux q corr reads, and is explained by the following consideration: The thermal flux is the mixture-average of the component thermal fluxes, q th = M a=1 Y a q th a , where q th a = −τ P C a,p ∇T is the Fourier law for the component, τ is a parameter of no importance to the current consideration. On the other hand, in the single-component LBM, the thermal flux is q th sc = −τ P ∇H sc , and with the single-component enthalpy H sc it returns the Fourier law in this case. However, the extension of the single-to the multicomponent case so far invokes only the replacement of the single-component enthalpy with the "lumped" mixture enthalpy and without any correction one gets, Thus, apart from the mixture-averaged Fourier law q th , the thermal flux also contains a spurious term. The spurious term is eliminated by the correction flux q corr (3.37), where the prefactor is chosen by considering the hydrodynamic limit; see Appendix A. The correction flux vanishes if all components are thermodynamically indistinguishable, that is, if all species have the same specific heat. In many cases, the correction flux contributes negligibly, for example, for air at moderate temperatures where the standardair assumptions for diatomic molecules holds to a good approximation.

Hydrodynamic limit of the two-population lattice Boltzmann model for mixtures
Constraints on the pertinent equilibrium and quasi-equilibrium moments (3.27), (3.28), (3.29), (3.30), (3.31), (3.32), (3.36), (3.37) and (3.33) are sufficient to study the hydrodynamic limit of the two-population lattice Boltzmann system (3.22) and (3.23) without a complete specification of the equilibrium and the quasi-equilibrium populations. The analysis follows the route of expanding the propagation to second order in the time step δt and evaluating the moments of the resulting expansion. Details of the derivation are included in Appendix A, here we present the final result: The continuity equation: (3.38) The momentum equation: Here, the pressure tensor π reads, where S is the strain rate, The dynamic viscosity µ and the bulk viscosity ς are related to the relaxation parameter ω, The energy equation: Here, the heat flux q reads, (3.45) The first term is the Fourier law of thermal conduction, with thermal conductivity λ related to the relaxation parameter ω 1 , The second term in (3.45) is the interdiffusion energy flux. Some comments are in order: • The continuity, the momentum and the energy equations are the standard equations for a multicomponent compressible mixtures (Williams 1985;Bird et al. 2006).
• The bulk viscosity vanishes if all components are monatomic,C a,v = DR U /2.
• Introducing the thermal diffusivity α = λ/ρC p and the kinematic viscosity ν = µ/ρ, the Prandtl number becomes, (3.47) • Using the equation of state of the mixture (2.4) in (3.42) and (3.46), relaxation parameters ω and ω 1 are expressed in terms of dynamic viscosity and thermal conductivity, . (3.49) Finally, the dynamic viscosity µ and the thermal conductivity λ of the mixture at any point is evaluated as a function of the local composition by using the methods described in Wilke (1950) and Mathur et al. (1967), respectively: where µ a are the dynamic viscosity of the components while the dimensionless factor φ ab is given by the equation, (3.51) The thermal conductivity of the mixture λ is calculated from the thermal conductivity of the components λ a , (3.52)

Equilibrium and quasi-equilibrium
In order to finalize the construction of the lattice Boltzmann equations for the mixture, one needs to specify the choice of the momentum and the energy lattices, and to provide the corresponding equilibrium and quasi-equilibrium populations. To that end, the singlecomponent lattice Boltzmann models satisfying the moment constraints of Sec. 3.2 are known in the literature. These employ higher-order lattices with a relatively large number of discrete velocities such as D2Q49 (Q = 49 in two dimensions, Frapolli et al. (2015)) or D3Q39 (Q = 39 in three dimensions, Frapolli et al. (2020)).
In this paper, we develop the standard D3Q27 lattice realization as in the above case of the species LBM of Sec. 2.4. We thus consider a two-dimensional compressible single-  For the evaluation of the equilibrium g eq i and of the quasi-equilibrium g * i of the energy lattice, we proceed with the following ansatz G i , parameterized with a scalar θ ∈ [0, 1], a scalar M 0 , a vector m and a second-order tensor M , (3.54) Here, the weights w i are calculated in the product form as, based on the fundamental triplet, (3.56) Furthermore in (3.53), Z is a vector with the components Here, M αα is the diagonal component of the second-order tensor M , while the components of vectors B i are defined as follows: (3.58) Note that, when the parameter θ is set to the lattice reference temperature θ = 1/3, the term in (3.57) vanishes and the remaining term (3.54) becomes the familiar secondorder Grad's approximation (Grad 1949). By construction, the form (3.53) satisfies the moment relations, for any θ: (3.59) The equilibrium and the quasi-equilibrium populations g eq i and g * i are defined with the help of the form (3.53) by specifying the parameters θ = RT , M 0 = ρE and M = R eq (3.30) in both cases, and m = q eq (3.29) or m = q * (3.32) for the equilibrium or the quasi-equilibrium, respectively, following Saadat et al. (2019): We shall now proceed with identifying the equilibrium of the momentum lattice and the modification of the lattice Boltzmann equation necessary for the D3Q27 model.

Augmented lattice Boltzmann equation for the momentum lattice
Equilibrium populations of the momentum lattice f eq i are evaluated in the conventional way with the help of the product-form (2.34) and using ξ = u and ζ = RT , It is well known that the diagonal element of the equilibrium third order moment of the momentum lattice Q eq ααα cannot satisfy the required moment relation (3.28). This happens due to the lattice constraint, c 3 iα = c iα , which makes the diagonal thirdorder moments linearly dependent on the momentum, cf., e. g. (Karlin & Asinari 2010). Following (Saadat et al. 2019), we consider the augmented lattice Boltzmann equation on the momentum lattice as follows, where X is the vector with the components α = x, y, z, and where the components of vectors A i are defined as, This completes the realization of the mixture momentum and energy lattice Boltzmann equations on the standard D3Q27 lattice. In the next section we shall specify the coupling between the species and the mixture lattice Boltzmann subsystems.
3.5. Matching of mixture density and momentum: Weak and strong coupling

Weak coupling
Summarizing, the lattice Boltzmann model for a compressible M -component mixture of ideal gas on the standard D3Q27 lattice consists of M species lattices where the lattice Boltzmann equation is given by Eq. (2.53), and the momentum and energy lattice Boltzmann equations (3.63) and (3.23). In total, the M + 2 lattice Boltzmann equations are tightly coupled, as has been already specified above: The temperature from the energy lattice is provided to the species lattices through species equilibrium (2.35) and quasiequilibrium (2.36), but also in the Stefan-Maxwell temperature-dependent relaxation rates (2.31). On the other hand, the mass fractions from the species lattices are used to compute the mixture energy and enthalpy in the equilibrium and the quasi-equilibrium of the momentum and energy lattices. Another coupling is the input of species diffusion velocities into the quasi-equilibrium population of the energy lattice via the interdiffusion flux (3.36). Finally, the momentum and the energy lattices are coupled in the standard way already present in the single-component setting. This entire set of interconnections between the species, and the momentum and energy lattices shall be termed the weak coupling.

Strong coupling
With the two subsystems, of the species and the mixture, first constructed independently from each other and after that being coupled in the way described above, we are left with two independent definitions of the mixture density and the mixture momentum: On the one hand, the mixture density ρ (3.16) and the mixture momentum ρu (3.17) are defined as the moments of the populations f i . On the other hand, the same quantities are defined with the species populations as the sum of partial densities and partial momenta. The number of the conservation laws for the species subsystem is M + D, while for the mixture subsystem it is D + 2. The total number of the conservation laws in the weakly coupled combined system is M + 2D + 2. Thus, the weakly coupled system is in excess of D + 1 conservation laws as compared to the M + D + 1 conservation laws of the mixture.
We consider the mixture density and mixture momentum defined by the momentum lattice as primary. The over-determination is then resolved by replacing the density of a selected component (here, the M th ) ρ M by the deficit of density once the M − 1 other components are taken into account. Similarly, the momentum of the M th component ρ M u M accounts for the deficit of the mixture momentum once the momenta of the other species are counted: In other words, the density and momentum of the M th component is no more an independent field but is slaved by the corresponding mixture quantities and the rest of the mixture composition. This means that the lattice Boltzmann equation for the M th component becomes purely relaxational, stripped of its conservation law. At the same time, the total momentum conservation is slaved by the momentum conservation of the momentum lattice. The number of independent conservation laws in this strongly coupled system is thus, (3.68) and corresponds to the locally conserved fields, ρ 1 , . . . , ρ M −1 , ρ, ρu and ρE. While the assignment of the "slaved" component M is not unique, it is advisable to select the component which carries the majority of mass in the mixture. It should be emphasised that the strong coupling is only necessary to avoid an overdetermined system. In practice, performing simulations under the weak coupling does not result in a deviation of mass or momentum except in cases where the boundary conditions introduce such a deviation. The matching of the density and the momentum of the mixture under the weak coupling condition thus highlights the intrinsic consistency of the proposed mean-field model.

Overview of numerical implementation
In order to validate various physical aspects of the proposed lattice Boltzmann model, we consider four benchmarks as follows: • Diffusion in a ternary gas mixture. This test case validates the Stefan-Maxwell diffusion and exhibits effects such as a diffusion barrier and uphill diffusion, which cannot be captured by Fick's diffusion assumption.
• Diffusion in opposed jets. Here, we verify the coupling between the hydrodynamics and the diffusion model.
• Speed of sound measurement in a mixture. This test case further validates the compressible model.
• Three-dimensional Kelvin-Helmholtz instability. Finally, this canonical benchmark demonstrates the extension of the fully coupled model to three dimensions and therefore shows viability for complex flows including turbulence. In all cases, we consider the temperature-independent specific heats for each component and set the reference temperature T 0 = 0 in (3.6), so that the internal energy is U a = C a,v T and the mixture internal energy is, The speed of sound c s is defined as, where both the adiabatic exponent γ = C p /C v and the specific gas constant R depend on the mixture composition. In what follows, we use the acoustic scaling: The speed of sound (4.1) at a specified reference composition (typically, at the equilibrium) and specified temperature is used to make velocity non-dimensional, unless otherwise stated. The characteristic length is given in the respective setup. Transport coefficients including dynamic viscosity, thermal conductivity and diffusivity were derived from the GRI-Mech 3.0 mechanism (Smith et al. 1999); acoustic scaling was used to convert these to lattice units. Finally, the second-order accurate isotropic lattice operators proposed by Thampi et al. (2013) were used for the evaluation of spatial derivatives in the correction to the heat flux (3.37) as well as in the isotropy correction (3.64).

Diffusion in a ternary gas mixture
Classical experiment on diffusion in ternary mixture of hydrogen H 2 , argon Ar and methane CH 4 in a Loschmidt tube apparatus was performed by Arnold & Toor (1967); Results of the experiment of Arnold & Toor (1967) were later analysed in depth by Krishna & Wesselingh (1997). Experiment of Arnold & Toor (1967) highlighted a number of, in part counter-intuitive, features of the Stefan-Maxwell diffusion. It is therefore natural to test our mixture model against the experiment of Arnold & Toor (1967).
The strongly coupled version of the three-dimensional lattice Boltzmann model for Stefan-Maxwell diffusion was realized on a quasi-one-dimensional domain with 864×1×1 grid points. In order to represent a closed tube, the bounce-back boundary condition was used for all populations at each end of the tube, while periodic boundary conditions were applied in the other two directions. The same initial composition of the mixture as in the experiment of Arnold & Toor (1967) was used; the setup was initialized with a uniform atmospheric pressure and temperature T = 300 K. It should be stressed that the temperature was not stipulated to be fixed during the simulation. Rather, the fully coupled thermo-hydrodynamic system maintained the isobaric and isothermal conditions by the initial and boundary conditions. As in the experiment, the evolution of the composition was represented by the average mole fraction of each species in the left and in the right halves of the tube. The non-dimensional time was used to represent the data, t N D = t/t s , where t s is the time required for the sound wave to traverse the domain in a reference equilibrium composition.
In our first numerical experiment, the mole fractions of the species in the left and in the right halfs of the tube were initiated as in the case 1T of Arnold & Toor (1967); see Fig. 1: Left X H2 = 0.491, X Ar = 0.509, X CH4 = 0.000 Right X H2 = 0.000, X Ar = 0.485, X CH4 = 0.515 (4.2) Time evolution of hydrogen H 2 and of methane CH 4 follows Fick's diffusion law: the species from the higher concentration side reduce in mole fractions as they move towards the low concentration side. Thus, CH 4 can be seen moving from right to left and H 2 in the opposite direction, both species eventually attaining a uniform concentration.
However, the behaviour of argon Ar cannot be explained by Fick's law. Although Ar has a negligible concentration gradient due to the initial conditions (4.2), it does start diffusing, see Fig. 2. This phenomenon was termed osmotic diffusion by Toor (1957). Osmotic diffusion is said to occur when the rate of diffusion of a component is not zero even though its concentration gradient is negligible; this would correspond to an infinite Fick's diffusivity. The concentration of Ar keeps on growing in the left section even though its concentration is higher in the left section itself, see Fig. 2. The effect was termed uphill diffusion (or reverse diffusion) in Toor (1957) because the component diffuses in the direction of increase of its concentration; in Fick's picture this would amount to negative diffusivity. The reverse diffusion is seen to proceed for some time and then flattens at t N D ≈ 400. At this point in time, an appreciable concentration gradient is built up but the diffusion is negligible, see Fig. 3. The effect was termed a diffusion barrier in (Krishna & Wesselingh 1997), the point at which the diffusion rate of a component vanishes even though its concentration gradient does not. This would mean zero Fick's diffusivity. After the diffusion barrier, the ordinary Fick's diffusion sets in and proceeds downhill of the concentration gradient until the uniform steady state is reached, see Figs. 4 and 5. Fig. 6 shows the evolution of the average mole fractions of the species in the left and in the right halves of the tube. The effects just mentioned were observed in the experiment of Arnold & Toor (1967) and in an earlier similar experiment by Duncan & Toor (1962). Krishna & Wesselingh (1997) provided explanation by drawing an analogy between the frictional drag and binary diffusion coefficients of pairs of species. According to Krishna & Wesselingh (1997), the Stefan-Maxwell diffusivity plays a role of an inverse drag coefficient. The binary diffusion coefficient between Ar and H 2 is 8.14543×10 −5 m 2 /s while that between Ar and CH 4 is 2.17321 × 10 −5 m 2 /s. This means that the frictional drag exerted on Ar by CH 4 is much greater than that exerted on Ar by H 2 . Thus, CH 4 drags Ar along with it during the initial period when the flux of CH 4 from right to left is large, causing the uphill diffusion of Ar. The transport of CH 4 from right to left eventually reduces because the driving force causing it reduces due to the reduction in concentration gradient of CH 4 . At the same time, the increasing concentration of Ar creates a driving force for Ar to diffuse downhill. A balance is reached at the point of diffusion barrier after which the drag force caused by CH 4 is overcome and Ar starts diffusing downhill of its concentration in a Fick's fashion.
It is apparent from Fig. 6 that the lattice Boltzmann simulation was able to correctly capture the experimentally observed phenomena. For a more quantitative assessment, we compare simulation results with the linearized theory of multicomponent mass transfer  (Arnold & Toor 1967), Eq. (4.2). Averaged mole fractions of hydrogen H 2 , argon Ar and methane CH 4 in the left half and the right half sections of the tube. The figure shows reverse diffusion of argon. Symbol: present simulation; Line: theory (Arnold & Toor 1967).
proposed in Arnold & Toor (1967). The theory relies upon a semi-analytical solution of the one-dimensional diffusion equations for the average mole fractions using linearized Stefan-Maxwell relation for the diffusion fluxes, and was shown to match the experiment in a quantitative fashion (Arnold & Toor 1967). For the purpose of this study, we numerically solved the equations of the linearized theory using Python. As is evident from Fig. 6, the results of the simulation agree well with the linearized theory, both in terms of magnitudes of the mole fractions as well as the time.
Continuing along the lines of the experimental study, the simulation was repeated with a different set of initial conditions, corresponding to the case 2T of Arnold & Toor (1967): Left X H2 = 0.512, X Ar = 0.000, X CH4 = 0.448 Right X H2 = 0.000, X Ar = 0.485, X CH4 = 0.515 (4.3) According to the theory of the inverse relation between mass diffusivity and drag (Krishna & Wesselingh 1997), methane should now exhibit uphill diffusion due to the flux of argon.  (Arnold & Toor 1967).
This is indeed what was seen in the experiments of Arnold & Toor (1967) as well as in our simulations, see Fig. 7. Simulations are in good agreement with the linearized theory.
A final but equally important situation is the one marked as case 3T in the experiment (Arnold & Toor 1967),

Left
X H2 = 0.512, X Ar = 0.000, X CH4 = 0.448 Right X H2 = 0.491, X Ar = 0.509, X CH4 = 0.000 (4.4) The binary diffusivity between Ar and H 2 is 8.14543×10 −5 m 2 /s while that between CH 4 and H 2 is 7.37433 × 10 −5 m 2 /s. The diffusivities are comparable and thus the interaction of H 2 with Ar is very similar to the interaction of H 2 with CH 4 . In Fig. 8, the results from the lattice Boltzmann simulation as well as the linearized theory show a nearly-Fickian diffusion of H 2 . Hydrogen however does show a small but nevertheless clear tendency to accumulate in the right half of the tube. This is possibly due to a slightly greater drag exerted on H 2 by CH 4 thanks to a somewhat smaller diffusivity between the pair.  (Arnold & Toor 1967), Eq. (4.4). Averaged mole fractions of hydrogen H 2 , argon Ar and methane CH 4 in the left half and the right half sections of the tube.. The figure shows near-Fickian diffusion of hydrogen. Symbol: present simulation; Line: theory (Arnold & Toor 1967).
For an additional validation, we compare the composition map of the simulation with the composition map of the experiment, since the composition map for the Stefan-Maxwell diffusion is independent of time (Duncan & Toor 1962). Fig. 9 verifies that the composition paths of the simulations agree well with that of all the three experiments of Arnold & Toor (1967). The composition path that would be followed by a purely Fickian diffusion is also marked in Fig. 9 as 'Path w/o reverse diffusion' for the purpose of contrast. The setups eventually attain a homogeneous composition at the 'Equilibrium' points on the composition map located midway of the Fickian lines, at the intersection with the Stefan-Maxwell trajectory. It can be again seen in the composition map that even for the case 3T , the diffusion path of hydrogen is almost yet not purely Fickian.

Diffusion in opposed jets
In order to assess the coupling between the diffusion and the hydrodynamic systems we consider the case of planar opposed jets.  Arnold & Toor, 1967: 1T Arnold & Toor, 1967: 2T Arnold & Toor, 1967: Composition path of the averaged mole fractions of H 2 , Ar and CH 4 in the left half and the right half section of the tube. The composition path shows three different cases, each with a reverse diffusion of Ar, CH 4 and a nearly Fickian diffusion of hydrogen. Lines: present simulation; Symbol: experiment of Arnold & Toor (1967).
that studied in Arcidiacono et al. (2007). It consists of two facing jets of fluid with equal momentum and different compositions. As shown in Fig. 10, the simulation is performed on a grid of size L x × L y × L z = 200 × 400 × 1 points, with the distance between the nozzles L x = 200. For the inlets, the incoming populations are replaced by the equilibrium distributions while the the outlets are modelled by making the derivative normal to the boundary zero. For the solid vertical boundaries at y < 0.1L y and y > 0.9L y , a free-slip boundary condition is used. The compositions of the jet streams are, Left X H2 = 0.1, X N2 = 0.85, X O2 = 0.0, X H2O = 0.05 Right X H2 = 0.0, X N2 = 0.90, X O2 = 0.1, X H2O = 0.00 (4.5) The solution of the present model at steady state is compared to the solution produced by the 'CounterflowDiffusionFlame' function of the open source package Cantera (Goodwin et al. 2018). This function computes a steady state solution to counter-flow diffusion flame using a reduced one-dimensional similarity solution, as derived in section 6.2 of Kee et al. (2005). In order to get a solution comparable with the Stefan-Maxwell formulation of the lattice Boltzmann model, the reactions are turned off in Cantera and the transport model is set to 'Multi', which accounts for the pairwise diffusion between the species. As can be seen in Fig. 11, the LBM solution for the mole fractions of all the components as well as the scaled velocity agree well with solution produced by Cantera. It should be noted that this test case is regarded severe in Arcidiacono et al. (2007) since the diffusion of the components proceeds against the velocity of the bulk flow. For example, hydrogen and water from the left nozzle diffuse against the bulk flow on the right side, upstream towards the right nozzle. The good agreement of the results indicates that the coupling of the Navier-Stokes and the Stefan-Maxwell models is consistent and correct.

Speed of sound
As a standard test for a compressible flow LBM, we verify that the model correctly reproduces the speed of sound c s (4.1). The speed of sound was measured for the following four compositions S1-S4: (S1) R = 0.046897 X H2 = 0.491, X Ar = 0.509, X CH4 = 0.000 (S2) R = 0.0333655 X H2 = 0.200, X Ar = 0.700, X CH4 = 0.100 (S3) R = 0.026625 X H2 = 0.000, X Ar = 0.900, X CH4 = 0.100 (S4) R = 0.0283563 X H2 = 0.200, X Ar = 0.100, X C3H8 = 0.700 (4.6) Composition S1 in (4.6) is chosen from one of the cases in section (4.2) whereas the case S3 is chosen arbitrarily in order to test for the sound speed in a composition with a considerable difference in mole fractions. The cases S2 and S4 are chosen to verify the speed of sound in a ternary mixture and in the presence of heavier gases such as propane, respectively. The test is performed by tracking a small perturbation in pressure ∆P = 10 −5 at a specified temperature T . Fig. 12 compares the measured speed of the propagation of the perturbation with the theoretical speed of sound prediction (4.1). The lattice Boltzmann model correctly recovered the sound speed over a tested range of temperatures from T min = 0.025 to T max = 0.8, in lattice units. The tested temperature range is characterized by the ratio of the temperatures, T max /T min = 32 which is sufficiently large for many applications. Temperature between T = 0.2 to T = 0.5 in lattice units was used for most of the simulations presented in this paper.

Kelvin-Helmholtz instability
Without a pretence of an in-depth study of shear layers instabilities in this paper, the final example presents a three-dimensional simulation of the classical Kelvin-Helmholtz instability (KHI) in order to validate the proposed model towards its possible use for high Reynolds number simulations. Similar to the setup in San & Maulik (2018), we simulate the Kelvin-Helmholtz instability in a periodic domain of the size L x ×L y ×L z = 800 × 800 × 200 lattice grid points. The domain is split into three sections in the flownormal direction, where the initial conditions for a two-component mixture of nitrogen N 2 and water vapor H 2 O are as follows: u x = 0.1M c , X H2O = 0.9, X N2 = 0.1, for 0 y < 0.25L y , u x = −0.1M c , X H2O = 0.1, X N2 = 0.9, for 0.25L y y < 0.75L y , u x = 0.1M c , X H2O = 0.9, X N2 = 0.1, for 0.75L y y L y . (4.7) The velocity in the normal direction u y and in the span-wise direction u z is given by, respectively, u y = 2|u x |0.01 sin(2πx/L x ), (4.8) The initial composition of the binary mixture (4.7) is so chosen as to equilibrate at the 50/50 equilibrium composition X N2 = X H2O = 0.5 in the absence of any flow. The initial condition (4.8) and (4.9) further introduces a small perturbation in both the normal and the span-wise directions with a wavelength equal to the length of the domain and a magnitude of one percent of the relative shear velocity. As defined by Leep et al. (1993), the convective Mach number M c is the Mach number relative to the frame of reference of the simulation. The relative Mach number M r based on the relative velocity across the shear layers is M r = 0.2 according to the initial conditions (4.7). The Reynolds number with respect to the viscosity of the bottom-most layer, with M r as the velocity scale and L y as the length scale is Re = 11963.46. The mole fractions are chosen 0.1 and 0.9 to make the test more severe. We define an eddy turnover time t e = L x /U r , with the initial relative velocity U r = 2|u x |.  We present contours and iso-surfaces of the equilibrium mole fraction of nitrogen X N2 = 0.5 at different times in Figs. 13 to 16. The iso-surfaces are colored by the x-component of the convective Mach number M c to provide a visual indication of the direction of motion.
Soon after the initial condition, the normal perturbation breaks the symmetry of the flow and the shear layer begins to curl up into a vortex without a significant span-wise deformation. This is evident in Fig. 13, where the three-dimensionality of the flow is visible only in the nonuniform span-wise velocity of the iso-surface. As the simulation proceeds, the flow in Fig. 14 develops anti-symmetric vortices which are also visibly deformed in the span-wise direction. This process continues and the vortices stretch and deform over time, which is visible in both the contours of the mole fraction of N 2 as well as the iso-surfaces in Fig. 15. The flow eventually becomes more chaotic, forming smaller-scale structures in Fig. 16. It should not be forgotten that the components  are also undergoing diffusion during this mixing, as evident from the smearing of the contour values over time. While the above observations are inline with what is typically observed in the literature, it is also important to verify the energy distribution across the scales of the flow. To that end, we measure the turbulent kinetic energy spectrum, which shows the expected −5/3 Kolmogorov scaling in the inertial subrange in Fig. 17. This additionally validates our model and outlines a path towards complex multicomponent flow simulations.

Conclusion
Let us consider a "good" lattice Boltzmann setting for a single-component gas. From the past experience, this would imply a two-population LBM because the second pop-ulation would be ultimately needed for an adequate description of the energy, leaving aside a special case of monatomic gas. It is then possible to envision, following the rule of the conventional kinetic theory (Chapman & Cowling 1990), a valid lattice Boltzmann model for a mixture of such gases, with the total number of kinetic equations equal to 2 × M for the M -component mixture since each component needs to be represented by its individual two-population LBM.
On the contrary, in this paper we proposed a lattice Boltzmann framework for multicomponent mixtures of ideal gases with a more realistic number of coupled lattice Boltzmann equations. We addressed two equally important aspects. First, we proposed a new LBM system for the Stefan-Maxwell diffusion and barodiffusion comprising M lattice Boltzmann equations. Second, we proposed a reduced, mean-field description of the mixture momentum and energy using the two-population setting. The resulting framework consists of M + 2 lattice Boltzmann equations rather than 2 × M as it would be if the detailed and not the mean-field approach to the energy of the mixture would have been pursuit.
Special attention was devoted to the consistent thermodynamic coupling of the above two sub-systems in such a manner that the hydrodynamic limit is not compromised. The proposed framework was realized on the standard three-dimensional lattice using an extension of the compressible model of Saadat et al. (2019). Specific to the multicomponent problem, the interdiffusion energy flux was added in a natural way to the heat flux to recover the correct energy equation while a counter-flux was introduced to remove the spurious contribution to the Fourier law inevitably arising with the mean-field approach to the energy description. While we focused on the multicomponent case, the proposed realization is also an extension of the augmented, compressible LB model proposed by Saadat et al. (2019) to a general form in three dimensions, which can of course also be used for single-component flows.
The simulation of the diffusion in a ternary mixture demonstrated that the proposed LBM correctly accounts for binary interaction between species. The coupling of diffusion to hydrodynamics was assessed by computing diffusion in opposed jets and the basic compressibility features were demonstrated through the speed of sound simulation at various compositions. Finally, the simulation of the three-dimensional shear layer instability in a binary mixture with a high composition contrast indicates that the proposed method can be useful for direct numerical simulations of complex flows.
All of the above gives us grounds to believe that the proposed multicomponent framework fills the gap in the development of the lattice Boltzmann method and is a first step towards reactive flow applications which will be the focus of our future studies.