Internal Energy Relaxation Processes and Bulk Viscosities in Fluids

: Internal energy relaxation processes in ﬂuid models derived from the kinetic theory are revisited, as are related bulk viscosity coefﬁcients. The apparition of bulk viscosity coefﬁcients in relaxation regimes and the links with equilibrium one-temperature bulk viscosity coefﬁcients are discussed. First, a two-temperature model with a single internal energy mode is investigated, then a two-temperature model with two internal energy modes and ﬁnally a state-to-state model for mixtures of gases. All these models lead to a unique physical interpretation of the apparition of bulk viscosity effects when relaxation characteristic times are smaller than ﬂuid times. Monte Carlo numerical simulations of internal energy relaxation processes in model gases are then performed, and power spectrums of density ﬂuctuations are computed. When the energy relaxation time is smaller than the ﬂuid time, both the two temperature and the single-temperature model including bulk viscosity yield a satisfactory description. When the energy relaxation time is larger than the ﬂuid time, however, only the two-temperature model is in agreement with Boltzmann equation. The quantum population of a He-H 2 mixture is also simulated with detailed He-H 2 cross sections, and the resulting bulk viscosity evaluated from the Green–Kubo formula is in agreement with the theory. The impact of bulk viscosity in ﬂuid mechanics is also addressed, as well as various mathematical aspects of internal energy relaxation and Chapman–Enskog asymptotic expansion for a two-temperature ﬂuid model.

A hierarchy of thermodynamic nonequilibrium fluid models may be derived from the kinetic theory of polyatomic gas mixtures. The most general thermodynamic nonequilibrium model is the state-to-state model, in which each internal state of a molecule is independent and considered as a separate species [1][2][3][4][5][6][7][8][9][10][11][12][13]. When the internal states are lumped into energy bins, coarse-grained models are obtained, with each bin considered a separate species [14,15]. When there are partial equilibria between some of the bins or between the internal states, species internal energy temperatures may be defined, and the complexity of the model is reduced accordingly [9]. The next reduction step consists of equating some of the species' internal temperatures with the equilibrium one-temperature model ultimately obtained [9]. Each relaxation step towards a simpler and more equilibrated model then yields bulk viscosity contributions-provided the characteristic energy relaxation times are smaller than the flow times, the most complete bulk viscosity coefficient being that obtained for the equilibrium one-temperature fluid [26][27][28].
A simplified kinetic model where elastic and resonant collisions are fast but collisions exchanging translational and internal energy are slow is first considered [26]. In such a framework, the translational and internal temperatures are macroscopic quantities associated with collisional invariants of the fast collision operator. In a relaxation regime, when the characteristic time of internal energy relaxation is smaller than the flow characteristic time, the difference between the translational temperature T tr and the equilibrium temperature T becomes proportional to the divergence of the velocity field v. This leads to the apparition of a bulk viscosity coefficient κ, such that nk B (T tr − T) = −κ∇·v where n is the number density and k B the Boltzmann constant.
A more complex situation with two internal energy modes, one with a slow exchange rate and the other one with a fast exchange rate, is then investigated [26]. The translationalrapid mode temperature and the slow mode temperature are then collisional invariants of the fast collision operator. For such a model, there is a bulk viscosity due to the fast internal energy mode, as in classical one-temperature models, but part of the thermodynamic equilibrium bulk viscosity is still hidden in the slow internal mode. A detailed analysis yields that, in a relaxation regime, there are five contributions to the effective bulk viscosity, namely the fast internal mode bulk viscosity, the slow internal mode bulk viscosity, the reduced relaxation pressure and two perturbed relaxation source terms. In the thermodynamic equilibrium limit, the sum of these five contributions coincide with the one-temperature two-mode bulk viscosity. The physical interpretation of the origin of the bulk viscosity coefficient is found to be similar to that of the simplified two-temperature model [26]. This analysis may also be generalized to the situation of gas mixtures, as well as to the case of nonindependent energy modes [27].
Next, the general situation of state-to-state mixture models in which each quantum state of each species is independent is considered [28]. Relaxation equations in symmetric form are derived for the quantum state population Gibbs functions and the translational temperature. Approximate solutions of the population relaxation equations compatible with the asymptotic equilibrium limit are then obtained. At zeroth order, using a relaxation approximation, the differences between the pseudo species' chemical potentials and their equilibrium value are proportional to the divergence of the velocity field. The 'internal energy' bulk viscosity κ [01] is then recovered with the relation nk B (T tr − T) = −κ [01] ∇·v. At first order, the relaxation approximation yields a bulk viscosity that converges at thermodynamic equilibrium towards the traditional bulk viscosity.
Monte Carlo simulations of spontaneous fluctuations near thermodynamic equilibrium are then performed in order to investigate a polyatomic model gas [36][37][38]. The density fluctuation power spectrum of the model gas is evaluated by using the Boltzmann equation, as well as linearized fluid equations [7,26,[39][40][41][42]. The simplified one-temperature model, including the bulk viscosity term, then well agrees with Boltzmann equation when the internal energy relaxation time is smaller than the flow time [26,27]. When the relaxation time is larger than the flow characteristic time, however, only the two-temperature model is in agreement with Boltzmann equation.
The simplified two-temperature model is considered in Section 2, the two-temperature two-mode nonequilibrium model in Section 3, and the state-to-state model in Section 4. The Monte Carlo simulations of single gases are presented in Section 5 and that of state-to-state mixture models in Section 6. The impact in fluid mechanics and the mathematical aspects are finally addressed in Section 7.

Kinetic Framework
A single polyatomic gas is considered with the Boltzmann equation written in the form where ∂ t denotes the time derivative operator, c the particle velocity, ∇ the space derivative operator, f (t, x, c, I) the distribution function, x the spatial coordinate, I the index of the quantum energy state, J rap the rapid collision operator, J sl the slow collision operator, and the formal parameter associated with the Chapman-Enskog procedure. The complete collision operator J = J rap + J sl is in the form where in a direct collision I and J denote the indices of the quantum energy states before collision, I and J the corresponding numbers after collision, c the velocity of the colliding partner, c and c the velocities after collision, a I the degeneracy of the Ith quantum energy state, g the absolute value of the relative velocity c − c, e the unit vector in the direction of the relative velocity c − c , and σ IJI J the collision cross sections [18,22]. The dependence of f on (t, x) has been left implicit in (2) and the cross sections σ IJI J satisfy the reciprocity relations a I a J gσ IJI J dc d c de = a I a J g σ I J IJ dc d c de.
Denoting by E I the internal energy in the Ith state, the rapid collisions are either elastic without change of internal energy or resonant with ∆E = E I + E J − E I − E J = 0, E I = E I and E J = E J , whereas the slow collisions are such that Denoting by J tr−tr the operator associated with elastic collision, J int−int the operator associated with resonant collisions, and J tr−int the operator associated with collisions, such that ∆E = 0, the fast and slow collision operators are then The collisional invariants of the fast collision operator J rap are associated with particle number ψ 1 = 1, momentum ψ 1+ν = mc ν , ν ∈ {1, 2, 3} with c = (c 1 , c 2 , c 3 ) t , kinetic energy ψ 5 = ψ tr = 1 2 m|c − v| 2 and internal energy ψ 6 = ψ int = E I where m denotes the particle mass and v the fluid velocity.
The Enskog expansion is in the form f = f (0) 1 + φ + O( 2 ) and the Maxwellian distribution f (0) reads where n denotes the number density, k B the Boltzmann constant, T tr the translational temperature, T int the internal temperature, and Z int = ∑ I a I exp −E I /k B T int the partition function. There are two different temperatures, T tr and T int , in f (0) , since there are two different energy collisional invariants: ψ tr and ψ int . The scalar product ξ, ζ between two tensorial quantities ξ(t, x, c, I) and ζ(t, x, c, I) is naturally defined by where ξ·ζ is the contracted product.
The equations for conservation of mass, momentum and internal energies are then obtained by taking the scalar product of the Boltzmann equation (1) with the collisional invariants of the fast collision operator. The corresponding fluid variables are the particle number density n = ψ 1 , f = ψ 1 , f (0) or, equivalently, the mass density ρ = mn, the mass averaged velocity v with ρv = mc, f = mc, f (0) , the internal energy per unit volume of translational origin E tr = f , ψ tr = f (0) , ψ tr , and the internal energy per unit volume of internal origin E int = f , ψ int = f (0) , ψ int , or, equivalently, the translation and internal temperatures T tr and T int defined by E tr (T tr , n) = f , ψ tr and E int (T int , n) = f , ψ int . The pressure p and the internal energies E tr and E int are found in the form where E = ∑ I a I E I Z int exp −E I /k B T int is the average internal energy per particle. The corresponding translational and internal entropies and Gibbs functions are presented in [26]. Following the Chapman-Enskog procedure, the equations for conservation of mass, momentum, and internal energies are obtained in the form [9] ∂ t ρ + ∇·(ρv) = 0, where ⊗ denotes the tensor vector product, I the unit tensor, Π the viscous tensor, Q tr the translational energy heat flux, Q int the internal energy heat flux and ω int 1 the first-order energy exchange term. The transport fluxes are given by where η denotes the shear viscosity, and λ tr,tr , λ tr,int , λ int,tr , and λ int,int the thermal conductivities. The full source term ω int = ψ int , J sl = ψ int , J may be expanded as where ω int 0 is evaluated from the Maxwellian distribution f (0) and δω int 1 is the correction associated with the Navier-Stokes perturbation φ, so that the first-order source term ω int 1 is given by Finally, the pressure tensor P = pI + Π is given by and does not involve a bulk viscosity term, unlike one-temperature polyatomic gas models [16][17][18][19][20][21][22][23][24][25].

Relaxation and Bulk Viscosity
From the energy Equations (8) and (9) it is obtained at zeroth order that where the heat capacities are given by The zeroth-order source term ω int 0 is in the form Defining the nonequilibrium correction factor by ζ = , the source term ω int 0 may be rewritten in the relaxation form Subtracting the T int equation from that for T tr and using (16), the resulting equation This is a typical relaxation equation, and when τ int is smaller that the flow characteristic time, the relaxation relation T tr − T int = −τ int p∇·v/nc vl is obtained after some initial layer. The equilibrium temperature is naturally defined as the unique scalar T, such that and only this temperature T is available for the limiting one-temperature fluid model. Let- and the equilibrium limit of this coefficient is c int , which coincides with the bulk viscosity coefficient obtained independently from the Chapman-Enskog method for the equilibrium fluid [34,35]. Combining the relaxation relation with c tr (T tr − T) = c int (T − T int ) and the definition of κ, we obtain the fundamental relation The pressure tensor P is then in the form and the bulk viscosity coefficient of the one-temperature equilibrium limit fluid that only involves T has been recovered. Many authors have discussed the near thermodynamic equilibrium situation, where the internal temperature T int and the translational temperature T tr are a priori close; notably, Kohler [16], Hirschfelder Curtiss and Bird [17], Waldmann [18], Chapman and Cowling [19], Ferziger and Kapper [20], McCourt et al. [21], de Groot and Mazur [24], Keizer [25], Zhdanov [8], Nagnibeda and Kustova [9] and Brun [10].
It is also possible to establish that first-order corrections to the bulk viscosity coefficient are negligible for such a simplified two-temperature model. Discarding Burnett-type terms, the corrections indeed involve the perturbed source term δω int 1 that is in the form where φ ω is the scalar perturbed distribution function arising from the expansion [26] with φ η denoting a symmetric traceless tensor; φ λ tr and φ λ int are vectors [9,26]. However, the standard scalar basis functions φ 0010 = 3 2 − 1 2 m k B T |c − v| 2 and φ 0001 = (E − E I )/k B T int used for scalar linearized kinetic equations are both collisional invariants of the rapid collision operator, and are therefore in the nullspace of the linearized fast collision operator I rap . Therefore, the perturbed term vanishes φ ω 0 in a first approximation, so that ω int 1 ω int 0 , and there are no first-order corrections to the bulk viscosity coefficient [26].

Kinetic Framework
A single polyatomic gas with two independent internal energy modes that have different exchange rates is now considered. The first mode is assumed to have a rapid exchange rate with the translational degrees of freedom, whereas the other one is assumed to have a slow exchange rate. The internal energy in the Ith quantum state is accordingly split as where I denotes the composed index I = (I rap , I sl ), I rap the index of the quantum energy state of the rapid mode, I sl the index of the quantum energy state of the slow mode, E rap I rap the rapid mode internal energy, and E sl I sl the slow mode internal energy. It is denoted for short E rap I for E rap I rap and E sl I for E sl I sl , so that E I = E rap I + E sl I . The degeneracy of the Ith state is also denoted by a I and may be decomposed as a I = a rap I rap a sl I sl where a rap I rap and a sl I sl are the degeneracies of the fast and slow modes. The rapid collisions are all collisions, such that ∆E sl = E sl I + E sl J − E sl I − E sl J = 0, either only involving the translational and rapid mode energy or resonant with respect to the slow internal mode. Denoting by J (tr+rap)−(tr+rap) the collision operator involving solely the translational and fast internal degrees of freedom, J sl−sl the operator for resonant collision with respect to E sl , and J (tr+rap)−sl the operator for collisions, such that ∆E sl = 0, the Boltzmann equation governing the distribution f (t, x, c, I) is in the form (1) with the fast and slow collision operators given by The collisional invariants of the fast collision operator are associated with particle number ψ 1 = 1, momentum ψ 1+ν = mc ν , ν ∈ {1, 2, 3}, translational and rapid mode energy ψ 5 = ψ tr + ψ rap and slow mode energy ψ 6 = ψ sl , where ψ tr = 1 2 m|c − v| 2 , ψ rap = E rap I , and where T is the partial equilibrium temperature between the translational and fast internal degrees of freedom, T sl the temperature associated with the slow internal energy modes, T sl the internal partition function. There are two different temperatures T and T sl involved in f (0) , since there are two different energy collisional invariants ψ tr + ψ rap and ψ sl . The partition function may be decomposed as Z int = Z rap Z sl where Z rap = ∑ I rap a rap I rap exp −E rap I rap /k B T and Z sl = ∑ I sl a sl I sl exp −E sl I sl /k B T sl are the fast and slow mode partition functions.
The equations for conservation of mass, momentum and internal energies are obtained by taking scalar products of the Boltzmann equation with the collisional invariants of the fast collision operator. The extra fluid variables to consider, in addition to the particle number density n and the mass averaged velocity v, are now the energies E tr+rap = f , ψ tr + ψ rap = f (0) , ψ tr + ψ rap and E sl = f , ψ sl = f (0) , ψ sl or, equivalently, the temperatures T and T sl defined by E tr+rap (T, n) = f , ψ tr + ψ rap and E sl (T sl , n) = f , ψ sl . The pressure p and the internal energies E tr+rap and E sl are obtained in the form Z sl exp −E sl I sl /k B T sl are the average fast and slow mode internal energy per particle. The corresponding entropies and Gibbs functions are presented in [26]. The corresponding mass and momentum conservation equations are similar to (6) and (7) and are not repeated. On the other hand, the equations for conservation of internal energies are in the form where Q tr+rap denotes the translational and fast mode energy flux, Q sl the slow mode energy flux, and ω sl 1 the first-order energy exchange term. The transport fluxes are given by where p rel denotes the relaxation pressure, κ rap the fast internal energy mode bulk viscosity, η the shear viscosity and λ tr+rap,tr+rap , λ tr+rap,sl , λ sl,tr+rap , and λ sl,sl the thermal conductivities. The full source term ω sl = ψ sl , J sl = ψ sl , J may be expanded as where ω sl 0 is evaluated from the Maxwellian distribution f (0) and δω sl 1 is the correction associated with the Navier-Stokes perturbation φ, so that ω sl 1 is given by Finally, defining the pressure tensor as P = pI + Π, it is obtained that with a pressure term nk B TI, a relaxation pressure term p rel I, a bulk viscosity contribution solely associated with the fast internal modes κ rap ∇·vI and a shear viscosity term. In particular, the resulting bulk viscosity κ rap differ from that obtained at equilibrium that involves the two internal energy modes [34,35].

Relaxation and the Slow Mode Bulk Viscosity
From the energy Equations (26) and (27), it is deduced at zeroth order that where the heat capacities are given by where is the averaging operator. Defining the nonequilibrium correction factor as ζ sl = , the source term ω sl 0 may be recast in the relaxation form Subtracting the T sl equation from that of T and using (34), the resulting equation for This is a typical relaxation equation and, assuming that τ sl is smaller than the flow characteristic time, the relaxation relation at zeroth order T − T sl = −τ sl p∇·v/nc vl is obtained after an initial layer. The equilibrium temperature T is naturally defined, such that Letting The slow mode bulk viscosity is defined by κ sl = pk B c sl τ sl /c vl c vl , where c vl = c tr + c rap + c sl and κ sl may also be written Combining the relaxation relation with the later definitions yields at zeroth order that Using the state law (25) and the expression of the pressure tensor (32), yields an effective bulk viscosity in the form κ rap + κ sl that differs from the one-temperature two-mode bulk viscosity directly obtained at equilibrium [34,35] also presented in Appendix A. It is thus necessary to investigate first-order effects in order to recover the full one-temperature, two-mode equilibrium bulk viscosity.

The Effective Bulk Viscosity
First-order corrections to the relaxation relation need to be taken into account in order to recover the proper one-temperature two-mode bulk viscosity at equilibrium. From the governing Equations (26) and (27), it is deduced that at first order so that the structure of the first-order source term ω sl 1 = ω sl 0 + δω sl 1 has to be investigated. From the structure of the linearized Boltzmann equation, the perturbed distribution function φ may be expanded as where φ η is a symmetric traceless tensor, φ λ tr+rap and φ λ sl are vectors φ κ and φ ω are scalars. The coefficients φ µ , µ ∈ {η, λ tr+rap , λ sl , κ, sl}, satisfy linearized Boltzmann equations in the form I rap (φ µ ) = ψ µ with the constraints f (0) φ µ , ψ j = 0, 1 j 6, where I rap is the linearized fast collision operator and the right-hand sides ψ µ are presented in [26]. Defining W sl = ∑ J,I ,J (∆E sl ) f (0) gσ IJI J d c de the perturbed source term is then in the form δω sl 1 = f (0) φ, W sl and using the Curie principle then yields δω sl In the relaxation approximation, at first order, one may also replace ω sl 0 with its zerothorder approximation ω sl 0 ≈ −c sl p∇·v/c vl in the first-order term δω sl 1 . The perturbed source terms w κ 1 and w sl 1 then yield corrections to the temperature difference T − T sl in such a way that [26] In addition, the relaxation pressure p rel is given by where p rel is the reduced relaxation pressure so that p rel is also proportional to ∇·v and further yields a bulk viscosity contribution. Finally, using the general expression of the pressure tensor (32), the definition of p rel , as well as (41), it is found that the effective bulk viscosity in the Navier-Stokes regime is in the form [26] The detailed calculation of each term finally yields after lengthy algebra that [26] The effective bulk viscosity at thermodynamic equilibrium T = T sl = T then coincides with the one-temperature two-mode bulk viscosity derived from the Chapman-Enskog method and presented in Appendix A.
Finally, it is possible to introduce a translational temperature T tr from the relation E tr (T tr , n) = f , ψ tr . This temperature is not a collisional invariant and must be expanded in the form T tr = T tr Combining (41) and (45), it is finally obtained that T tr 1 and κ eff , defined by (43) and given by (44), are such that Therefore, the physical interpretation gained with the previous simplified model (20) is also valid in the more complex situation where there are two energy modes with different dynamics [26]. All these results may further be extended to mixtures of gases, as well as to the situation where the fast and slow energy modes are not independent [27].

Kinetic Framework
A state-to-state model for a mixture of polyatomic gases is considered [1][2][3][4][5][6][7][8][9][10][11]13]. We denote by iI the pseudo species index for the ith species in the Ith quantum state, S = {1, . . . , N s } the species indexing set, N s the number of species, Q i the ith species quantum state indexing set, Q = ∪ i∈S {i}×Q i the pseudo species indexing set, and N q the number of pseudo species. The pseudo species Boltzmann equations are written as where c iI denotes the velocity of the particle of the iI-th pseudo species, f iI (t, x, c iI ) the distribution function for the iIth pseudo species, J rap iI the rapid collision operator, J sl iI the slow collision operator, and the formal parameter associated with the Chapman-Enskog procedure. The internal energy of the iIth pseudo species is also denoted by E iI and the corresponding degeneracy by a iI . Using similar notation as in previous sections, the fast and slow collision operators are given by The collisional invariants of the fast collision operator J rap = (J rap iI ) iI∈Q are associated with the pseudo species particle numbers ψ kK = (δ ki δ IK ) iI∈Q , kK ∈ Q, momentum ψ N q +ν = (m i c iIν ) iI∈Q , ν ∈ {1, 2, 3}, where m i denotes the mass of the particles of the ith species and c iI = (c iI1, c iI2 , c iI3 ) t , and kinetic energy where v denotes the mass average mixture velocity.
The Enskog expansion reads f iI = f where T tr is the translational temperature of the mixture. Assuming that there are sufficiently inelastic collisions, the collision invariants of the slow collision operator J sl = (J sl iI ) iI∈Q , are associated with the species particle numbers ψ k = (δ ki ) iI∈Q = ∑ K∈Q k ψ kK , k ∈ S, momentum ψ N q +ν , ν ∈ {1, 2, 3}, and mixture total energy ψ en = ψ N q +4 = 1 2 m i |c iI − v| 2 + E iI iI∈Q . The scalar product ξ, ζ between two tensorial quantities ξ = (ξ iI ) iI∈Q and ζ = (ζ iI ) iI∈Q is naturally defined by where ξ iI ·ζ iI is the contracted product.
The fluid equations are obtained by taking the scalar product of Boltzmann Equation (47) with the collisional invariants of the fast collision operator. The fluid variables are the pseudo species number densities n kK = ψ kK , f = ψ kK , f (0) or, equivalently, the mass densities ρ kK = m k n kK for kK ∈ Q, the mass averaged velocity , ψ tr with n = ∑ kK∈Q n kK denoting the total number density. The pressure p, the translational energy E tr and the total internal energy E are found in the form Introducing the partition functions where h P is the Planck constant, one may then rewrite the Maxwellian distribution of the iI-th pseudo-species in the form f The entropy per particle of the iIth pseudo species is given by for iI ∈ Q and the mixture fluid entropy by S = ∑ iI∈Q n iI S iI . The Gibbs function G iI of the iIth pseudo particle is also given by G iI = k B T tr log n iI Z iI and the reduced chemical potential µ iI and µ tr by and will be used as symmetrizing variables. Following the Chapman-Enskog procedure, the conservation equations are found in the form [9] ∂ t n kK + ∇·(n kK v) + ∇·(n kK V kK ) = ω where V kK denotes the diffusion velocity of the kKth pseudo species, ω (1) kK the first-order production term of the kKth pseudo species, Π the viscous tensor and Q the heat flux.
The inelastic collisions are written for convenience as chemical reactions where M iI is the symbol of the iI pseudo species, r the reaction number, R the reaction indexing set and ν f iIr and ν b iIr the forward and backward stoichiometric coefficients of pseudo species iI in reaction r. The source term may be expanded as ω kK is evaluated with Maxwellian distributions and δω (1) kK is the first-order perturbation. The zeroth-order source terms are first obtained in the form ω kKr − ν f kKr , and the rates of progress in the where C r is an average quantity associated with chemical transition probabilities of reaction r [9,22]. Letting iI∈Q the zeroth-order source term may then be written as where µ e = (µ e iI ) iIQ denotes the equilibrium value of µ and where the property ν r , µ e = 0, r ∈ R, has been used. The perturbed source terms δω (1) kK may also be expressed in terms of the perturbed distribution φ [22,26,27].
Denoting by I rap iI the linearized fast collision operator for the iIth pseudo species and by I rap = (I rap iI ) iI∈Q the linearized operator of the mixture, the perturbed distribution and φ must satisfy the Enskog constraints f (0) φ, ψ j = 0 for 1 ≤ j ≤ N q + 4. Expanding the perturbed distribution function in terms of the gradients of the macroscopic variables, the dissipative fluxes are found in the classical form [9,22,26]. In particular, the viscous tensor reads Π = −η ∇v + (∇v) t − 2 3 (∇·v)I and there is neither a bulk viscosity term nor a relaxation pressure, since all pseudo species have a single internal state [9,22]. Defining the pressure tensor as P = pI + Π, it is obtained that so that P does not contain a bulk viscosity term, unlike one-temperature polyatomic gas mixtures [16][17][18][19][20][21][22][23][24][25].

Equilibrium Population and Bulk Viscosities
The equilibrium limit of the state-to-state model is obtained by zeroing the chemistry sources ω (0) while maintaining constant slow variables associated with the collision invariants of the full collision operator J rap + J sl , namely the species number densities n k = ∑ K∈Q k n kK , k ∈ S, momentum ρv and the total internal energy n 3 2 k B T tr + ∑ iI∈Q n iI E iI . Denoting by T the thermodynamic equilibrium temperature, the species equilibrium population is given by the Boltzmann distribution where Z int i = ∑ I∈Q i a iI exp −E iI /k B T is the equilibrium internal partition function of the ith species and where ∑ I∈Q i n iI = ∑ I∈Q i n e iI = n i , for i ∈ S. The corresponding equilibrium species average energies are E i = ∑ I∈Q i (a iI E iI /Z int i ) exp(−E iI /k B T), the internal heat capacities are given by and the equilibrium temperature is the unique scalar T, such that With the chemistry analogy, one may further introduce ψ k = (δ ik ) iI∈Q = ∑ K∈Q k ψ kK , for k ∈ S, and ν r = (ν iIr ) iI∈Q , for r ∈ R, and define Using a chemistry vocabulary, one may then say that the species are the atoms of the pseudo species. Assuming naturally that there are sufficiently energy exchanges then yields R ⊥ = A; that is, the reaction vectors ν r , r ∈ R are sufficiently numerous in order to span the maximum space A ⊥ , keeping in mind that there is no dissociation or recombination. The equilibrium pseudo species number densities µ e = (µ e iI ) iI∈Q are then the unique solution of the equilibrium conditions µ e ∈ R ⊥ or equivalently µ e ∈ A under the constraints that n k for k ∈ S and E are invariants. The reduced chemical potentials at equilibrium are indeed such that µ e iI = log(n e iI /Z e iI ) = log(n i /Z int i ) where Z e iI = Z iI (T) so that µ e = ∑ i∈S log(n i /Z int i ) ψ i and µ e , ν r = 0, r ∈ R. The usual scalar species basis functions at equilibrium are given by φ 0001k = (E i − E iI )δ ik /k B T iI∈Q and φ 0010k = ( 3 2 − 1 2 m i k B T |c − v| 2 )δ ki iI∈Q for the kth species. The basis function φ 0010k is defined for any k ∈ S, while φ 0001k are defined for any k ∈ S p where S p denotes the set of polyatomic species. The projected basis functions φ 0001k = φ 0001k − ψ en c int k /c vl , k ∈ S p , are also considered where the terms proportional to ψ en ensure that the functions φ 0001k are orthogonal to the collisional invariant ψ en of the complete collision operator in the equilibrium kinetic framework. Two bulk viscosities at equilibrium may then be defined; namely the 'internal energy' bulk viscosity κ [01] as well as the 'standard bulk viscosity' κ.
The linear system associated with the evaluation of the 'internal energy' bulk viscosity κ [01] at equilibrium is obtained by using the Galerkin variational approximation space spanned by φ 0001k , k ∈ S p . The idea behind this basis function is that the most important part of the dynamic is that associated with internal energy exchanges and not with kinetic energy [34,35]. The 'internal energy' bulk viscosity κ [01] is obtained by solving the transport is the classical bracket product [34,35]. The matrix K [01] is symmetric positive definite and letting γ = (γ i ) i∈S , the 'internal energy' bulk viscosity is given by κ [01] = − ∑ i∈S pγ i x i c int i /c vl and has been found to be accurate [34,35]. On the other hand, the standard bulk viscosity is obtained with the Galerkin variational approximation space spanned by φ 0001k and φ 0010k , k ∈ S.

Symmetric Zeroth-Order Relaxation Equations
The system of partial differential equations governing the pseudo species Gibbs functions and the translational temperature at zeroth order is written in symmetric form. Symmetrized forms are convenient for analyzing systems of partial differential equations modeling fluids, and are generally obtained by using entropic-type variable −(∂ u S) t where u denotes the vector of conservative variables. Since we are interested in the relaxation of thermodynamic variables, for convenience, we use the variables (δµ, δµ tr ) t = (µ − µ e , µ tr − µ e tr ) t where µ = (µ iI ) iI∈Q , µ e = (µ e iI ) iI∈Q , µ tr = −1/(k B T tr ), and µ e tr = −1/(k B T). Denoting for short d t = ∂ t + v·∇ the convective derivative, the governing equations at zeroth order are found in the form where a tr = ∑ iI∈Q n iI ( 3 2 k B T tr + E iI ) 2 + 3 2 pk B T tr . Denoting by N the diagonal matrix N = diag (n iI ) iI∈Q and a the vector with components a iI = n iI ( 3 2 k B T tr + E iI ), iI ∈ Q, these equations involve the symmetric positive definite matrix A = N a a t a tr . The Equations (60) and (61) imply that δµ = µ − µ e satisfies the vector partial differential Equation where N = N − a⊗a/a tr is symmetric positive definite, Λ symmetric positive semidefinite, and b = (b iI ) iI∈Q is given by b iI = n e iI (E iI − E i )/Tc vl , iI ∈ Q. Since N is positive definite and Λ is positive semidefinite, (62) is a typical vector relaxation equation and the corresponding vector relaxation relation is then in the form where Λ e is the matrix Λ at equilibrium Λ e = ∑ r∈R Λ e r ν r ⊗ν r with Λ e r = C e r exp µ e , ν f r = C e r exp µ e , ν b r and ζ e r = 1. Since the nullspace of Λ e is A, one also needs constraints to determine uniquely δµ, using ∑ I∈Q i (n iI − n e iI ) = 0, i ∈ S. In the relaxation approximation, one may linearize the constraints around equilibrium, and after some algebra, it is obtained that N e δµ,

Bulk Viscosity at Zeroth Order
Taking into account the relaxation relation (63) and the mass constraints, one is lead to decompose δµ in the vector form The term δµ R is such that Λ e δµ R = b and b is in the range R(Λ e ) since ∑ I∈Q i n e iI (E iI − E i ) = 0. An approximate solution of the constrained relaxation equations that is compatible with the equilibrium limit mixture is now sought, so that δµ R = ∑ i∈S γ i φ 0001i and δµ A = ∑ i∈S γ i ψ i . Using that δµ R ∈ (N e ) −1 A ⊥ , one first obtains from the mass constraints that In order to determine the coefficients vector γ = (γ i ) i∈S , a least square approximation of the relaxation equations is used upon writing that ∑ l∈S φ 0001i , Λ e φ 0001l γ l = φ 0001i , b for i ∈ S.
The matrix K = φ 0001i , Λ e φ 0001j i,j∈S is found to be positive definite, and a direct comparison yields that the coefficients of the both matrices K [01] and Λ e are proportional Λ e φ 0001i , φ 0001j = npK [01]ij and φ 0001i , b = −nx i c int i /c vl , so that φ 0001i , b = nβ [01]i and γ = pγ [01] . It has thus been established that the least square solution to the relaxation equations under the linearized constraints is given by where γ = (γ i ) i∈S is the solution of the transport linear system K [01] γ = nβ [01]i associated with the internal mode bulk viscosity κ [01] .
By definition of the equilibrium temperature T, one has c tr (T tr − T) + ∑ iI∈Q n iI (E iI − E i ) = 0, and linearizing the expression of the reduced potential (65) yields Multiplying by n e iI (E iI − E i ) and summing over I and summing over i ∈ S, it is finally obtained that [28] The 'internal energy mode' equilibrium bulk viscosity κ [01] has thus been recovered with the relaxation relation (67) [28]. Although the kinetic framework for state-to-state mixtures of gases is much more complex than that of previous two-temperature models, a similar physical interpretation of the bulk viscosity coefficient has been obtained with (20), (46), and (67).

Bulk Viscosity at First Order
A similar analysis may be conducted for first-order equations, but this is much more intricate analytically. The first-order relaxation equations for the pseudo species are in the form ω (0) iI + δω (1) iI = −b iI ∇·v, for iI ∈ Q, and require us to evaluate the perturbed source term δω (1) in the neighborhood of equilibrium. The resulting analytical expressions and the resulting relaxation approximation have been obtained, as well as identification of the equilibrium limit with the traditional bulk viscosity, notably using the basis vectors φ 0010k , k ∈ S, and φ 0001k , k ∈ S p , the blocks K [01] = K 0101 K 0110 , K 1010 and K 1001 of the K matrix, as well as the Schur complement K [01] = K 0101 − K 0110 (K 1010 ) −1 K 1001 .

Numerical Experiments for the Simplified Two-Temperature Model
The results derived in Section 2 are assessed against numerical experiments for a model gas. Results are obtained by solving the appropriate Boltzmann transport equation via Monte Carlo methods [7,26,37,[39][40][41][42]. The transport properties of the model system are investigated by looking at the spontaneous fluctuations at thermal equilibrium [37,[57][58][59]. Interestingly, the dynamics of spontaneous fluctuations can actually be probed by light scattering experiments [59][60][61] and molecular simulations used in order to estimate bulk viscosity coefficients [62][63][64].

Kinetic Theory of Spontaneous Fluctuations
A fluctuating gas is considered near equilibrium, and the dynamics of the fluctuations of a variable A(r, t) is investigated by using the space-time correlation function where < ... > means ensemble average and δA(r, t) = A(r, t)− < A(r, t) > is the fluctuation of the dynamic variable. For an isotropic system in thermodynamic equilibrium, the correlation function depends only on the space-time distance In particular, the quantity actually measured in light (or neutron) scattering experiments is the Laplace-Fourier transform of the correlation function of density fluctuations, the spectral density (or power spectrum) of these fluctuations [59] δn 2 (k, ω) = e i(k·r−ωt) δn 2 (r, t)drdt.
Since the equilibrium fluctuations of the fluid variables are small compared to the average values, their dynamics are governed by the same equations that govern the dynamics of the system, but linearized around the equilibrium solution. These linearized equations are then doubly Laplace-Fourier transformed to the (k, ω) space and are solved for δρ k (s = + iω). The latter is used to construct the space-time correlation function < δρ * k (0) δρ k (s) >.
Finally, this correlation function may be connected with the density fluctuation power spectrum S(k, ω), a quantity that is experimentally measurable in light-scattering experiments For thermal fluctuations in gases, the ratio of the fluctuation wavelength to the mean free path defines the flow regime (from high to low ratios: hydrodynamic, kinetic, collisionless). Different regimes are described by different values of the parameter y = (8/3 √ 2π)ρ 0 √ k B T/m/ηk, where ρ 0 is the equilibrium density, η is the shear viscosity and k is the fluctuation wavenumber. The collisionless limit corresponds to y → 0, whereas the hydrodynamic limit (k → 0) is approached for y 5. The thermal fluctuation power spectra in the hydrodynamic regime is derived as follows from the one-temperature model, as well as the two-temperature model fluid equations.

Simulation of Spontaneous Fluctuations in a Dilute Gas
The fluctuation power spectrum for a one-temperature fluid described Navier-Stokes equation is reported in the book by J.P. Boon and S. Yip [57], whereas for the simplified two-temperature model, it is reported in [26]. The general method used to derive such fluctuation in a power spectrum is to start from the fluid equation, to linearize near equilibrium, to take the Fourier transform in space and then the Laplace transform in time, and the results are typically in the form where M(k, s) and N(k, s) are polynomials in s. On the other hand, at the molecular level, thermodynamic fluctuations in gasesprovided the density is low enough that only bimolecular collisions are effective-are described by the Boltzmann equation. In the case of the Boltzmann equation, the system is described in terms of the one-particle distribution function. By linearizing the equation around the equilibrium distribution, an integro-differential equation for the space-time correlation of the fluctuations of the distribution function is obtained [65]. The density fluctuations are then readily obtained via integration over the velocity space.
For the simulation of the spontaneous fluctuations in a gas in thermodynamic equilibrium, the Direct Simulation Monte Carlo method [36][37][38] is used. DSMC is a particle simulation method that solves the nonlinear Boltzmann equation. As such, it can simulate flows in the rarefied and/or hypersonic regime that cannot be dealt with in the framework of a fluid-dynamic treatment. It can also handle situations of strong thermal nonequilibrium, where a clear hierarchy of relaxation times cannot be established and rate equation methods fail [7]. The principle of the method is the decoupling, over a small time-step, of the processes of free flight and of collisional relaxation. A number of simulated particles are moved in the simulation domain according to their velocities and prescribed boundary conditions. In the collision step, particles are made to collide inside spatially homogeneous cells. A Monte Carlo method is used to realize collision events with the appropriate frequency. The details of the molecular processes occurring in the gas system are specified by assigning the appropriate set of collision cross sections. The viscosity and diffusion coefficients of the gas can be modeled by the Variable Soft Sphere model of Koura [66].
The power spectrum of the fluctuations of the dynamic variable n(r, t) is evaluated as follows. The variable fluctuations at all sampled space-time points, δn ij = n ij − n 0 , are recorded during the simulation, n 0 being the equilibrium value. This discrete set is then Fourier transformed and squared to get the discrete power spectrum. For an isotropic medium, it is sufficient to simulate a one-dimensional spatial domain.
Concerning the simulation parameters, the typical requirements of DSMC simulations are met: • Cell width to mean free path ratio: 0.3 • Timestep to mean collision time ratio: 0.05 • Average number of simulated particles per cell: 20 As for boundary conditions, the simulation domain is kept in contact with an infinite reservoir of the gas at the specified equilibrium conditions. Care must be taken to accurately sample the statistics of the incoming particles [68].
Additionally, for obvious reasons, the number of simulated particles is much less than the number of real particles present in the physical volume. The ratio of real to simulated particle number is called the weight w of the simulated particle, and it is a constant throughout the simulation. Now, since the density fluctuations are proportional to the gas density [67], i.e., given the volume, to the number of particles, the simulated fluctuations are equal to the real fluctuations to within a factor w. Therefore, the spectrum sampled by the simulation is exactly equivalent, to within normalization factors, to the spectrum measured in light-scattering experiments. In order to reduce the statistical scatter inherent in the particle simulation method, ensemble averaging of the results is performed by averaging the results of many independent runs. This procedure also allows us to estimate the variance of the results with respect to the statistical scatter. It is worth mentioning that this procedure is amenable to implementation on a computational grid. A number of wavelengths must thus be simulated in order to increase the signal-to-noise ratio. For each spectrum, we sample 4096 time points and 64 wavelengths (corresponding to 2048 spatial points). With these parameters, we obtain a signal-to-noise ratio of 10 −3 , i.e., the high-frequency background noise lies 3 orders of magnitude below the maximum value. The GRID infrastructure allows hundreds of runs to be performed simultaneously, thus drastically reducing the global computational time. The simulations of this work, in particular, have been performed under the Compchem Virtual Organization, and more details on the computational procedures may be found in [26].

Simulations for a Model Gas
A single gas of Hard Spheres is considered with mass = 28.9641 amu, σ = 7.2 · 10 −19 m 2 . The gas has two internal energy levels with degeneracies and energies given by Molecules exchange internal energy in collision according to the simplest single-quantum, line-of-centers model where k = 1 2 µg 2 is the kinetic energy in collision. The temperature and density are chosen to be T = 285.71 K and n = 2.4 · 10 21 m −3 . The fluctuation spectra are sampled at the wavelength 2π/k = 0.02 m that gives y = 5.97, so that the probed fluctuations fall into the hydrodynamic regime.
Two situations are analyzed. In the first case, the value p 0 = 0.01 has been chosen and gives for the relaxation time τ 0 ≈ 7.0 · 10 −5 s (Z ≡ τ 0 /τ c = 49). Figure 1 shows the fluctuation power spectra for this case. The spectra are normalized to unit maximum value. In this case, the relaxation is slow enough that a relaxation approximation does not hold, and the one-temperature model fails to describe the transport properties of the system correctly. The two-temperature model, instead, gives an adequate description of the system behavior and the agreement with the DSMC simulations is satisfactory. For comparison, we also reported the spectrum predicted for the same gas when the internal energy relaxation is forbidden (frozen relaxation). Next, a situation where relaxation of internal energy is fast enough compared to the flow characteristic time (as determined by the speed of sound) is analyzed. In these conditions, we expect the one-temperature model to be accurate and that the two-temperature model will reduce to the former. The value p 0 = 0.1 has been chosen and gives for the relaxation time τ 0 = 7.0 · 10 −6 s (Z ≡ τ 0 /τ c = 4.9). Figure 2 shows the fluctuation power spectra as obtained from DSMC simulations and from the one-temperature and two-temperature models, respectively. For comparison, we also show the spectra predicted for the same gas when the bulk viscosity contribution is neglected (i.e., κ = 0). In this case, it can be seen that both models describe the DSMC results accurately. Comparison with the κ = 0 case also shows that this agreement is not trivial, since there is an important contribution of the bulk viscosity to the spectrum with κ η ≈ 1. Note, however, that the one-temperature model cannot describe the (small) change in the speed of sound that is a consequence of the finite relaxation time for internal energy.
The statistical error in the simulation results is around 4% for the slow relaxation and 12% for the fast relaxation case, respectively. Therefore, multi-temperature hydrodynamic equations, as derived from the Boltzmann transport equation, provide an adequate description of internal energy relaxation for all values of the relaxation time. Therefore, there is no need to invoke frequency-dependent transport coefficients that introduce unnecessary complications. Further, the results support the conclusion, obtained by kinetic theoretical arguments in the previous sections, that the multi-temperature model reduces to the onetemperature model when the relaxation time is small enough and that in this case, a bulk viscosity formalism is adequate. These results are also relevant in view of the renewed interest in Rayleigh-Brillouin scattering in gases made possible by the use of nonlinear optical techniques [69]. Coherent Rayleigh-Brillouin scattering is a technique capable of making localized and high signal-tonoise ratio measurements of gases from the collisionless limit to the hydrodynamic regime. CRBS data are therefore expected to become a valuable source for the study of kinetic processes in molecular gases.

Internal Energy Spectrum and Energy-Exchange Collisions
The kinetic model that has been developed for arbitrary gas mixtures S is specialized to S = {He, H 2 } for the application. The required detailed cross sections used in this work has been calculated by the quasiclassical method, with an in-house developed code, which has been tested repeatedly against accurate results from the literature [43][44][45][46][47]. The set is complete, since all the H 2 rovibrational states of the electronic ground state are considered initial and final states. Quasibound states and dissociation processes have also been considered in the trajectory calculations, even though they have not been used in the present study. Cross sections for the processes He + H 2 (v, j) → He + H 2 (w, k) with v/w initial/final vibrational states, j/k initial/final rotational states, have been calculated including both reactive (i.e., exchange) and nonreactive processes. Particular care has been taken in terms of the accuracy of trajectory calculation, with a strict checking at each step of each trajectory, in order to accurately determine the optimal time step and improve the overall computational efficiency [46]. The potential energy surface (PES) adopted in this study is the well known Muchnik-Russek [48], instead of the more recent BMP [49], used, for example, in a similar work by Kim et al. [50]. This choice is motivated by the important discrepancies found with respect to experimental data in Lee et al. [51]. Six billion trajectories, using 9.5 years of CPU, have been calculated in this way on the Muchnik-Russek potential energy surface [48], using a constant density of 50,000 trajectories per eV of collisional energy (uniformly distributed) and per Å of impact parameter, in the range 1 meV-10 eV, with stratified sampling applied. Comparisons with available accurate quantum-mechanical theoretical data [52] put in evidence a very good reliability starting from 0.1-0.5 eV, depending on the initial states, corresponding to a minimal temperature of about 2000 K for rate coefficients [53]. For lower values of translational energies, cross sections tend to be less reliable, due to problems typical of nonreactive low-energy processes treated with QCT. A specific paper on this topic is in preparation. Finally, elastic collision integrals have been taken from Ref. [54].
The theoretical results of the previous sections are here specialized to the He-H 2 mixture in thermal equilibrium conditions. The various contributions of the quantum state population as a function of temperature have been evaluated numerically. The procedure is analogous to that used in a previous study on the H-H 2 mixture [27], to which the reader is referred for further details. Since a complete set of inelastic cross sections is available for the atom-diatom collisional system only, only a simulated gas is investigated, where H 2 -H 2 collisions are elastic. Since the resulting theoretical predictions are not amenable to experimental measurements, we resort to DSMC simulation and use Green-Kubo formulas for estimating the transport coefficients of the simulated gas.
A He-H 2 mixture in thermal equilibrium conditions has been simulated with a standard DSMC code using a majorant frequency scheme [7]. For this application, a spatially homogeneous simulation is sufficient, so the code is run neglecting the particle movement, and no spatial mesh or boundary conditions need be specified. The VSS model has been used for the elastic collision cross sections, with parameters chosen to reproduce the transport coefficients in Ref. [54]. The number of simulated particles is 2000 for all cases discussed. Results have been averaged over a number of independent runs in order to reduce the statistical scatter r (typically ≈ 100). The latter is reported in the results as the simulation's error bar.

Green-Kubo Bulk Viscosity in DSMC Simulations
Linear response theory is a powerful tool for the description of the relaxation towards thermodynamic equilibrium of any system subject to small perturbations. The fluctuationdissipation theorem, in particular, connects the time correlation functions of mechanical quantities with the system transport properties. These relations are known as the Green-Kubo expressions for the transport coefficients [58]. This theory is independent of the mechanical model describing the system. It has been primarily applied to the evaluation of transport properties in Equilibrium Molecular Dynamics simulations of liquids [55] and solids (see, e.g., [56] and references therein). The Green-Kubo formulas have been applied to the evaluation of transport coefficients in DSMC simulations of dilute (i.e., ideal) gases.
For a system of volume V in equilibrium at temperature T and pressure p, the Green-Kubo formula for the bulk viscosity reads [57]: where < · · · > denotes ensemble average, whereas that of the shear viscosity reads The quantities J p xx or J p xy , the currents, are, respectively, any of the diagonal or offdiagonal components of the spatially averaged, time dependent, rank-two pressure tensor: where the sum runs over simulated particles and C = c i − v.
The equilibrium condition allows us (via the ergodic hypothesis and stationarity) to express the integral in Equation (77) as the zero-frequency limit of the current power spectrum: A p xx (ω) being the current Fourier transform.
In DSMC simulation, the current J p is sampled at discrete time points, then Fourier transformed with standard FFT algorithms and squared. Let P η v be the resulting zerofrequency value: t sim being the simulation duration, ∆ t the sampling time interval and w the weight of simulated particles (each particle representing w real particles). The factor w arises because the fluctuation amplitude is proportional to the square root of the number of simulated particles; the factor ∆ 2 t comes from the discrete Fourier transform.

Results
The DSMC simulations have been performed for an equimolar mixture of helium and hydrogen with n = 10 20 m −3 . The resulting bulk viscosity obtained from DSMC is presented in Figure 3, together with the equilibrium viscosity as a function of temperature. The equilibrium bulk viscosity has been evaluated by using the expression is the usual averaging operator and ∆E the energy jump during a collision of the simulated gas [9,26,27]. The DSMC calculations have been performed at temperatures T = 6000 K, T = 9000 K, and T = 10,000 K. This figure reveals the very good agreement between the bulk viscosity of the fluctuating quantum population and the equilibrium limit, illustrating the accuracy of the least square formulation with a variational space 'compatible with the equilibrium mixture'.

Impact in Fluid Mechanics
The impact of the viscosity has scarcely been addressed in the past literature [70][71][72][73]. Karim [70] pointed out the importance of bulk viscosity coefficients and underlined that Stokes relation is only justified for monatomic gases. The impact of bulk viscosity on hypersonic boundary layers has been investigated in depth by Emmanuel [71,72]. A review of the concept of bulk viscosity and its implications for fluids over the twentieth century has been given by Graves and Argrow [73].
More recently, Billet et al. have investigated the interaction of a shock wave with a hydrogen bubble [74]. The bulk viscosity coefficient has been found to have a major influence on both the fluid and thermal aspects. This impact originates in particular from the thickening of pressure waves by bulk viscosity and from vorticity production when pressure and density gradients are not aligned. The artificial success of the Stokes approximation has also been related to the gradient structure of the term of ∇·(κ∇·v I) = ∇(κ∇·v). This term is indeed absent in boundary layers at second order, and only has a Ma 2 influence over fluid variables. As a typical example, in a steady hydrogen-air flame, κ/η is not small; ∇·v is not small either, but the gradient structure essentially leads to replacement of the hydrodynamic perturbed pressure p by p − κ∇·v, as discussed in Billet et al. [74]. Fru et al. have further investigated the small Mach situation and established that for a long-time integration the bulk viscosity regains its major influence [75]. The impact of bulk viscosity on compressible homogeneous turbulence has been studied by Chen et al., and the compressibility of the flow found significantly reduces when bulk viscosity is involved [98]. Two-and three-dimensional simulations of the interaction of a shock wave with a shear layer have further been investigated by Boukharfane et al. [76]. Finally, the interaction of a shock wave with a droplet has recently been studied by Singh et al., and bulk viscosity has again been found to play a major role [77].
Recent investigations have also shown that large values of bulk viscosity coefficients in dilute carbon dioxide mixtures result from erroneous applications of relaxation relations out of their domain of validity [11,13].

Chapman-Enskog Method for the Simplified Two-Temperature Model
The system of partial differential Equations (6)-(9) modeling two-temperature fluids derived in Section 2 may be written in the convenient vector form where ∂ i is the partial derivation in the ith spatial direction, D the indexing set of spatial directions, and u the conservative variable given by In these Equation (82), the matrix A i is the Jacobian A i = ∂ u F i of the convective flux in the ith spatial direction F i , whereas the dissipative matrices B ij are related to the dissipative fluxes F diss i with F diss i = − ∑ j∈D B ij (u )∂ j u . The convective flux F i and dissipative flux F diss i in the ith spatial direction and source term Ω are given by 1 where and Q tr = (Q tr 1 , Q tr 2 , Q tr 3 ) t . The convective matrices A i and i ∈ D, and the the dissipation matrices B ij , which contain all the transport coefficients, are investigated in [94][95][96]. Both the diffusion terms and the internal energy relaxation time have naturally been rescaled in the form B ij and τ int = τ int where is a typical Knudsen number. Two different small parameters could more generally be introduced with d in front of the dissipation terms and for the relaxation of internal energy in (82), and an independent limit with d and may also be investigated [94][95][96]. However, it is also legitimate and convenient to investigate the simpler situation where d = [94][95][96][97]. The dependence of the solution on the parameter has been emphasized by denoting the conservative variable in the form u .
The concept of Chapman-Enskog expansion has also been extended to hyperbolic systems of partial differential equations by Liu [82] and Liu, Chen, Levermore [84] and later to hyperbolic-parabolic systems by Giovangigli and Yong [94][95][96][97]. The natural structural conditions are that there exists a mathematical entropy, taken for convenience as σ = −S/R, that is compatible with convection, diffusion as well as source terms [82,84,[94][95][96][97]. It is also required that there exists an equilibrium manifold or slow manifold E, characterized by the relation T tr = T int for the two-temperature problem, and the slow variable corresponding to (83) reads u eq = ρ, ρv, E + 1 2 ρ|v| 2 t . For each slow variable u eq , there exists a unique u eq , such that Π t eq u eq = u eq where Π eq is the projector operator on the slow manifold with matrix Π eq = [e 1 , . . . , e 4 ] with e i , 1 ≤ i ≤ 5, denoting the canonical basis vectors of R 5 , and for later use, we also introduce the projector with matrix Π rap = [e 5 ]. The slow manifold is thus parameterized by u eq , and a careful analysis shows that the equilibrium one-temperature fluid model exactly corresponds to a second-order Chapman-Enskog [94][95][96] of the partial differential Equation (82). The schematic diagram of Figure 4 is thus obtained The Chapman-Enskog expansion for the solution u to the system of partial differential Equation (82) is notably in the form where the zeroth-order term u 0 coincides with u on the slow manifold Π t eq u = Π t eq u 0 = u eq , where the Enskog constraints are in the form Π t eq u i = 0 for i ≥ 1 and where u i only depends on ∂ α u eq for |α| ≤ i. Combining (82) with (87), the equations for u i , i ≥ 0, are given by At the order −1, it is found that Ω(u 0 ) = 0, so that u 0 is an equilibrium point u 0 = u eq (u eq ) and the slow variable equations are given by The equation at zeroth-order governing u 0 are found to be a symmetrizable hyperbolic system; indeed, Euler equations for the one-temperature fluid. The perturbed term u 1 is the next solution of linearized equations with the Enskog constraints in the form and using the Euler equation in order to express ∂ t u 0 , the first-order perturbation u 1 is found in the form where M j and j ∈ D are 5x4 matrices. The reader is referred to [94][95][96][97] for more details. The first-order Navier-Stokes-type equations are then obtained as with the diffusion coefficients in the form The resulting first-order accurate governing equations for the slow variable thus involves dissipative coefficients arising from perturbed convective terms Π t eq A eq i M j , as well as inherited directly from the system out of equilibrium Π t eq B ij . Both the bulk viscosity term arising from the perturbed convective fluxes and the shear viscosity term inherited from the out-of-equilibrium viscous tensor are finally involved in the equilibrium viscous tensor. For the two-temperature model, it is indeed found that so that the bulk viscosity coefficients are exactly obtained from the perturbed out-ofequilibrium convective terms in the Chapman-Enskog method, whereas the shear viscosity and thermal conductivity contributions are directly inherited from the out-of-equilibrium model [94][95][96][97].

Multiple Time Expansions for the Simplified Two-Temperature Model
Exact mathematical convergence results have further been obtained for the simplified two-temperature system of equations in the form (82) over the full space R d for 1 ≤ d ≤ 3 [96,97]. The results and the method of proof differ depending on whether the initial data are well prepared with T tr data = T int data or ill prepared with T tr data = T int data , where (ρ data , v data , T tr data , T int data ) t denotes the variable at initial time t = 0. Both the well-prepared situation [96] and the ill-prepared situations [97] have been investigated, and only the ill-prepared situation results are summarized in this section.
The existence of solutions to the Cauchy problem with ill-prepared initial data has been established, as well as the validity of asymptotic composite expansions including initial-layer correctors. The results may be conveniently presented using the variable w = ρ, ρv, E tr + E int + 1 2 ρ|v| 2 , with the slow u and fast q components given by Multiple time expansions have been introduced in the form with an outer expansion involving the standard time t as well as initial layer correctors involving the fast time τ = t/ in the form and such that w il 0 and w il 1 are exponentially decreasing with respect to the fast time τ = t/ . The equations governing each of the asymptotic expansion coefficients u 0 , u 1 , q 0 , q 1 and q 2 , where w 0 = (u 0 , q 0 ) t , w 1 = (u 1 , q 1 ) t , w 2 = (0, q 2 ) t may be obtained. The proper initial and boundary conditions may also be written by first expanding the initial conditions where w data = (u data , q data ) t and then identifying w 0 (·, 0) + w il 0 (·, 0) = w data0 and w 1 (·, 0) + w il 1 (·, 0) = w data1 , so that in particular, u 0 (·, 0) = u data0 , q 0 (·, 0) = 0, u il 0 (·, 0) = 0, and q il 0 (·, 0) = q data0 . The mathematical structure of the hyperbolic-type equations governing u 0 and u 1 then guarantee the existence of solutions over a finite time interval [0,t ] under natural regularity assumptions [94][95][96][97]. Notably, an important role is played by the mathematical entropy σ that allows symmetrization of the corresponding systems of partial differential equations, and that is taken in the form σ = −S/R for convenience. In addition, the initial layer correctors are shown to satisfy systems of differential equations with respect to the fast time. The existence of an exponentially decreasing global solution is then obtained by using the entropy as a Lyapunov function [97]. An important tool has also been the use of an improved truncated approximation in the form w a = w 0 + w 1 + 2 w 2 + w il 0 + w il 1 , including a fast component of the second-order term w 2 = (0, q 2 ) t It is then established that the out-of-equilibrium solution exists over a finite time interval [0,t ], and the truncated approximation is second-order accurate with sup t∈[0,t ] |w (t) − w a (t)| L ∞ ≤ C 2 where C is a finite constant and |f| L ∞ = sup x∈R d |f(x)| denotes the L ∞ norm of a function f over R d .
It may also be established that the Hilbert-type expansion u 0 + u 1 and the Chapman-Enskog solution at equilibrium u eq are O( 2 ) close, so that sup t∈[0,t ] u eq (t) − u 0 (t) + u 1 (t) L ∞ ≤ C 2 .
In other words, the out-of-equilibrium slow variable u is O( 2 ) close to the Chapman-Enskog solution at equilibrium u eq , up to a first-order term u il 1 (·, t/ )) that is exponentially decreasing with respect to the fast time τ = t/ [97]. The zeroth-order initial layer corrector is also governed by the ordinary differential equation in the form ∂ τ (T tr − T int ) = − c vl c tr τ int (T tr − T int ) so that assuming for the sake of illustration that c int and τ int constant yields that T tr (τ) = T + T tr (0) − T(0) exp − c vl c tr τ τ int , and q 1 = T int τ int c tr (T tr ) 2 +c int (T int ) 2 ∇·v that yields the bulk viscosity contributions of (91). The mathematical theoretical results obtained concerning the Chapman-Enskog solution are thus in full agreement with the physical framework. Drawing a parallel with the traditional Chapman-Enskog, one may say that the slow or fluid variable is the onetemperature fluid variable u eq , the Mawellian distribution is the corresponding equilibrium two-temperature fluid variable u eq , the linearized Boltzmann equations are replaced by (88) and the resulting dissipative coefficients have two different sources with (90).

Conclusions
The relaxation of internal energy and the concept of bulk viscosity have been investigated in nonequilibrium gas models derived from the kinetic theory. Two-temperature models of single fluids and state-to-state models for mixtures of gases have been considered. When the rates for internal energy exchanges are smaller than the fluid time, a common interpretation of the apparition of the bulk viscosity coefficient has been obtained with (20), (46), and (67). The Monte Carlo simulations of spontaneous fluctuations near thermal equilibrium obtained by solving the Boltzmann equation have been shown to be in full agreement with the fluid models, provided the internal energy relaxation times are smaller than the fluid time. Numerical simulation of He-H 2 quantum populations also yields bulk viscosity coefficients in agreement with the equilibrium. Mathematical aspects of internal energy relaxation have also been discussed for the two-temperature fluid model and were found to agree with the formal expansions. In particular, Chapman-Enskog expansions from systems of partial differential equations have been discussed, as well as rigorous results in the situation of ill-prepared initial data in full agreement with the formal expansions. A practical consequence of these results is that bulk viscosity should be systematically included in the numerical simulation of compressible fluids. Future work should also consider bulk viscosity coefficients arising in coarse-grained models, as well as numerical simulations with multiple internal energy modes and, more generally, states far from thermodynamic equilibrium.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The One-Temperature Two-Mode Bulk Viscosity
In this section, we investigate the bulk viscosity associated with a one-temperature twomode polyatomic gas. The standard linear system associated with the evaluation of the twomode bulk viscosity is obtained with the Galerkin variational approximation space spanned by the orthogonal polynomials φ 0010 = 3 2 − 1 2 m k B T |c − v| 2 , φ 0001rap = (E rap − E rap I )/k B T, and φ 0001sl = (E sl − E sl I )/k B T. The general solution of the Transport Linear Systems, as well as their mathematical structure, has been investigated [34,35]. The modes are termed 'rapid' and 'slow' for notational consistency with the nonequilibrium framework of Section 3, but in the thermodynamic equilibrium framework, they are all fast. The corresponding linear system of size 3 is in the form [35] where K denotes the system matrix, K the constraint vector, α = (α 10 , α 01rap , α 01sl ) t the unknown vector, β = (β 10 , β 01rap , β 01sl ) t the right-hand side vector and the bulk viscosity is finally given by κ = α 10 β 10 + α 01rap β 01rap + α 01sl β 01sl . The matrix K is positive semidefinite with nullspace N(K) = RV, where V = (1, 1, 1) t , the constraint vector is given by K = (c tr , c rap , c sl ) t , and the right-hand side vector by β = (c rap + c sl , −c rap , −c sl ) t /c vl . We deduce from the constraint K, α = 0 that c tr α 10 + c rap α 01rap + c sl α 01sl = 0 and κ = −(c rap α 01rap + c sl α 01sl )/c tr . We may thus recast the singular linear system of size 3 into a regular linear system of size 2 involving only the unknowns α 01rap and α 01sl . Thanks to the vector relation KV = 0, the coefficients of this linear system may also be expressed solely in terms of K rap,rap , K rap,sl , and K sl,sl . We also have the relations K rap,rap After some lengthy algebra, using the reduced linear system of size 2, it is obtained that κ = 1 (c vl ) 2 (c rap ) 2 K sl,sl − 2c rap c sl K rap,sl + (c sl ) 2 K rap,rap K rap,rap K sl,sl − K rap,sl K rap,sl . (A2) Since we have to investigate the equilibrium limit of a two-temperature model in which one mode is fast and another one slow, we deduce that the coefficient K rap,rap is large, and that the cross terms K rap,sl = K sl,rap are also small. We therefore neglect the square term K sl,rap K rap,sl in the previous expression, and the limiting value of the effective nonequilibrium bulk viscosity should thus be