{'fraction_non_alphanumeric': 0.03578531713194837, 'fraction_numerical': 0.011033201225266031, 'mean_word_length': 4.431207127651186, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How to build coarse-grain transport models consistent from the kinetic to fluid regimes

In this paper, we examine how to build coarse-grain transport models consistently from the kinetic to fluid regimes. The internal energy of the gas particles is described through a state-to-state approach. A kinetic equation allows us to study transport phenomena in phase space for a non-homogeneous gas mixture. Internal energy excitation is modeled using a binary collision operator, whereas the gas chemical processes rely on a reactive collision operator. We obtain an asymptotic fluid model by means of a Chapman-Enskog perturbative solution to the Boltzmann equation in the Maxwellian reaction regime. The macroscopic conservation equations of species mass, mixture momentum, and energy are given, as well as expressions of the transport properties. Reversibility relations for elementary processes are formulated in the coarse-grain model at the kinetic level and are enforced in the collision algorithm of the direct simulation Monte Carlo method used to solve the kinetic equation. Furthermore, respecting these reversibility relations is key to deriving a fluid model that is well-posed and compatible with the second law of thermodynamics. Consistency between the kinetic and fluid simulations is assessed for the simulation of a shock wave in a nitrogen gas using the uniform rovibrational collisional coarse-grain model. The kinetic and fluid simulations show consistency for the macroscopic properties and transport fluxes between both regimes.

(Dated: 12 November 2020) In this paper, we examine how to build coarse-grain transport models consistently from the kinetic to fluid regimes. The internal energy of the gas particles is described through a state-to-state approach. A kinetic equation allows us to study transport phenomena in phase space for a non-homogeneous gas mixture. Internal energy excitation is modeled using a binary collision operator, whereas the gas chemical processes rely on a reactive collision operator. We obtain an asymptotic fluid model  [11][12][13] and, at most, one-dimensional flows behind inviscid normal shocks 14 . Electronic-specific state-to-state CFD models have been developed to simulate electronic excitation and partial ionization in gases such as argon 15 , although in this case the number of discrete internal energy levels was much smaller than for the aforementioned molecular systems. Equivalent DSMC studies are even less common. Bruno et al 16 were the first to incorporate QCT-derived vibrationalspecific N 2 -N cross sections into a DSMC solver to study internal energy exchange and dissociation of nitrogen across a normal shock. To date, the only state-to-state DSMC simulations using a complete set of rovibrational-specific reaction cross sections for the N 2 -N system have been carried out by Kim and Boyd 17 .
One appealing way to reduce the computational cost of state-to-state nonequilibrium reacting flow calculations has been to develop coarse-grain models. The details of the reduction vary, but can broadly be classified into vibrational-specific, 18,19 energy bin, [20][21][22] hybrids of both, 8,23 or more recently adaptive grouping of rovibrational states. [24][25][26] The basic concept is always to approximate the behavior of the full kinetic database with a much smaller set of cross sections/rate coefficients by grouping together many individual processes. In addition to air chemistry, the approach has also been applied to electronic-specific simulations of argon plasma 27 . On the DSMC side, coarse-grain models have been investigated as well. 28,29 In every case, the lumping-together of internal energy levels leads to a reduction in the number of associated state-to-state reaction rate coefficients / cross sections and greatly reduces the cost of simulations.
But with the size reduction of the kinetic databases also comes a loss of fidelity in the thermodynamic and chemical-kinetic description, especially if the binning strategy chosen is inadequate. The problem was recognized early on with the uniform rovibrational collisional (URVC) bin model 20 . Its main underlying assumption of "freezing" the populations of rovibrational levels within a given bin to constant values (short recap in Sec. II) was shown to produce a set of mass balance equations at the coarse-grain level, which would be incompatible with micro-reversibility relations linking forward and backward rates between rovibrational states. This meant that the fluid equations could not satisfy the second law of thermodynamics, and the formulation of an associated entropy equation with unambiguous non-negative entropy production terms was not possible. Thus, one would not be guaranteed to retrieve the correct thermodynamic equilibrium with the URVC model. Furthermore, it was shown 30 that for the original formulation of the URVC bin model a large number of bins (up to 50) was required to approach the chemical dynamics of the full state-to-state system.
For the inviscid fluid limit, this problem was first remedied by allowing the level populations within each bin to assume a Boltzmann population at the gas temperature T (so-called Boltzmann binning 21,22 ), while the bin populations themselves could still be relaxing toward equilibrium. With this change, the gas would be guaranteed to eventually assume the correct equilibrium populations regardless of the number of bins used. Building on this approach, socalled Maximum Entropy grouping 23-26 adopted the use of bin-specific, or "group-internal" temperatures to allow the level populations within each bin even more degrees of freedom for local adjustment. This made it possible to almost exactly match the thermodynamic equilib-rium and reaction rates of full state-to-state calculations with as few as two, or three overall bins. As a consequence, most research has so far concentrated on refining the coarse-grain models to best approximate the full chemical kinetics in the inviscid limit.
In the few cases where viscous phenomena have been taken into account, 22,31 the transport properties were assumed to be independent of the molecules' internal energy states and were computed based on the current state-of-the-art collision integrals 32 . It has however been theorized [33][34][35][36] that transport properties should exhibit a dependence on the internal energy distributions. Indeed, in the state-to-state framework such a dependence appears naturally when deriving the Navier-Stokes equations as asymptotic solutions to the Boltzmann equation. Therefore, it must be taken into account for a well-posed fluid model. None of the coarse-grain models proposed so far have really addressed this issue. The aforementioned Boltzmann bin and Maximum entropy reductions also imply that the partition function of each energy bin must be temperature-dependent. This may work well within the context of a CFD fluid model, where detailed balance relations involve ratios of temperature-dependent forward and backward rate coefficients, but breaks down in DSMC 37 , where these same relations have to be expressed in terms of collision energy-dependent cross sections. Such coarse-grain models effectively require individual molecules in the gas to "be aware" of the surrounding temperature, which is not compatible with the kinetic-scale description.
In this work we formulate a coarse-grain model, which works within the context of DSMC and simultaneously allows for the derivation of viscous fluid equations with consistent transport terms. We propose a small, but important change 38 to the early URVC bin model 20 .
Instead of attempting to enforce micro-reversibility relations for all rovibrational levels, we will postulate that such relations must hold only for the bin populations. We will therefore sacrifice some of the fine-grain detail of the original system in exchange for a simpler coarse-grain model. Our model still assumes constant populations for all rovibrational levels lumped into a given bin. Its main usefulness lies not so much in the ability to reproduce the full chemical kinetics with the smallest number of bins, but in the rather simple manner with which reversibility relations can be expressed in terms of coarse-grain cross sections.
Furthermore, it allows us to postulate a Boltzmann equation at the coarse-grain level, which forms the starting point for a straightforward application of the Chapman-Enskog method to derive the fluid equations. As part of this, we obtain expressions for the transport properties directly based on the same coarse-grain cross sections appearing in the kinetic equation. As a consequence, the resulting transport properties are fully consistent with the corresponding coarse-grain DSMC collision model 29  (3) Verify the consistency of the hydrodynamic (Euler, or Navier-Stokes eqs.) and kinetic (Boltzmann eq.) coarse-grain models by simulation of normal shocks in nitrogen with CFD and DSMC methods. Assess the degree to which continuum breakdown across the shock causes the flow fields in the hydrodynamic and kinetic solutions to depart from one another.
This paper is organized as follows. In Sec. II, we introduce the coarse-grain model for inelastic processes in molecular gas mixtures and recall its main features. In Sec. III, we discuss the governing kinetic equation and detail its constituting terms, including reversibility relations. In Sec. IV, we apply the Chapman-Enskog method to derive the corresponding fluid equations, along with expressions for all necessary transport and chemical source terms.
In addition, we show that the entropy production terms due to transport and chemistry are always non-negative and thus the coarse-grain fluid equations satisfy the second law of thermodynamics. In Sec. V, we apply the coarse-grain model to reveal the structure of normal shock waves in a reacting gas mixture using three distinct simulation techniques. We first obtain the flow field in the inviscid limit by solving the system of master equations coupled to total momentum and energy balances behind the shock front. Then, we solve the full fluid equations across the shock with added viscous terms (Navier-Stokes) by means of the Finally, in Sec. VI, we state the conclusions and discuss possible future work.

II. COARSE-GRAIN MODEL FOR N3 SYSTEM
Throughout the remainder of this paper, we will consider as an example a mixture of molecular and atomic nitrogen, with both species in their ground electronic states, and use a set of cross sections derived from QCT calculations on an ab initio PES for the N 2 (v, J) + N system, originally compiled at NASA Ames Research Center 6,41 . The 9390 rovibrational levels of the N 2 molecule in its ground electronic state have been grouped together into a much smaller number of discrete internal energy bins according to the uniform rovibrational collisional (URVC) energy bin model 20,29,38 . As a result, our mixture is composed of energy bins (labeled with indices k = 1, 2, . . . , N bins ) that encompass all bound-and pre-dissociated levels i ∈ I N 2 , plus atomic nitrogen in the ground electronic state.
The number of molecules per unit volume populating bin k is the sum over all level populations belonging to it, i.e.: n k = i∈I k {n i }. The core assumptions in the URVC reduction are that the level populations within a bin are fixed by the relation n i = n k a i /a k and that each bin possesses an internal energy defined as the weighted average over the energies of its constituting rovibrational levels, i.e.: E k = i∈I k {a i E i }/a k . Here, the overall degeneracy of each bin is the sum over degeneracies of all rovibrational levels belonging to it: a k = i∈I k {a i }.
The set K N 2 contains the indices pointing to every one of the bins. For simplicity, it is assumed that atomic nitrogen only occupies a single internal energy state. The full set of N s = 1 + N bins (pseudo)-species in the mixture then becomes S = {N, N 2 (k) ∀ (k ∈ K N 2 )}.
This reduction effectively replaces the highly resolved representation of molecular nitrogen's thermodynamic state provided by the full set of level populations n i , (i ∈ I N 2 ) with a similar, but lower-resolution one that only relies on the bin populations n k , (k ∈ K N 2 ). By applying the URVC binning approach, the level-specific reaction rate/cross section data from the Ames database are condensed into bin-resolved rate coefficients/cross sections, first for inelastic collisions between molecular and atomic nitrogen: Here we have labeled the forward rate coefficient for the transition between molecules populating bins N 2 (k) and N 2 (l) as k Ef k→l (i.e. when Eq. (1) is read from left to right), whereas the backward rate coefficient in the opposite sense is labeled as k Eb k→l . Second, we have dissociation/recombination of an N 2 (k)-molecule by collision with an N-atom: where we have labeled the rate coefficients for dissociation k Df k and k Db k recombination respectively. Third, N 2 (k) + N-collisions in which no transition to another bin occurs, are referred to as "intra-bin scattering". Such processes are the equivalent of elastic collisions in our framework, since no internal energy is exchanged. We make the rather strong as- However, this does not constitute a problem for the purposes of our comparisons, as long as the simplification is done in a consistent manner when evaluating the collision terms of the Boltzmann equation and when calculating the dissipative fluxes and chemical reaction rates in the Navier-Stokes equations.

COARSE-GRAIN MODEL
At the kinetic scale, the evolution of the gas mixture is governed by a system of Boltzmann equations: Here, f i = f i (x, c i , t) are the velocity distributions of the N-atoms and the N 2 (k)molecules populating each one of the discrete internal states k ∈ K N 2 . The distributions depend on position x in physical space, particle velocity c i and time t. The term The terms on the right hand side are the collision operators. Together they account for any changes in the velocity distributions due to collisions between the mixture species.
The precise mathematical form of these operators depends on the collision types considered. In the present work we take into account the processes listed in Table I. There the collision types have been sub-divided into so-called fast and slow processes, based on their relative time scales. The fast scattering processes are responsible for driving the mixture toward a Maxwell-Boltzmann distribution at a common kinetic temperature T (i.e. thermalization) and for diffusive transport phenomena, whereas the slow processes can be either excitation/deexcitation reactions (responsible for relaxation of internal energy) and molecular dissociation-recombination reactions. The slow processes typically involve some energy threshold and any individual collision is far less likely to produce a significant change in the colliding particles' states than the fast collision types. This is reflected in the relative sizes of the associated cross sections. The fast processes possess differential cross sections typically orders of magnitude greater than the slow ones, i.e. σ slow ≪ σ fast . The sub-division into fast and slow processes is of little concern when Eq. (3) is solved directly, e.g. by means of the DSMC method. However, as discussed in Sec. IV B, the associated difference in time scales is exploited to derive the corresponding governing equations at the hydrodynamic scale.

A. Fast collision operators
The fast collision operator in Eq.
The partial collision operators: all possess the same structure for the fast processes listed in Table I. The integral in Eq. (4) is short notation for a three-fold integral over velocity space, plus a surface integral over the unit sphere. We take the dependence of f i on x, c i and t to be implicit. The variables Slow collision processes represent the velocity distributions of pseudo-species i and j evaluated at the "postcollision" particle velocities c ′ i and c ′ j respectively (i.e. the right-hand side of the collision as written in Table I). Conversely, the unprimed f i represent the distribution evaluated at the "pre-collision" particle velocities (left-hand side) in the same table.
The collision operator in Eq. (4) is made up of two competing terms: one involving the product f i f j , which accounts for depletion (negative sign) of f i due to collisions in the forward sense, and another one involving f ′ i f ′ j accounts for simultaneous replenishment (positive sign) by inverse collisions. The term in parentheses is multiplied in Eq. (4) by the magnitude of the pre-collision relative velocity g = |c i − c j | and the differential scattering cross section σ ij = σ ij (g, ω).
The differential cross section may in general depend both on g and on the orientation For the fast processes reversibility is enforced at the coarse-grain level and we postulate that the cross sections at both "ends" of the collision must verify the relation: where This is what allows us to combine the contributions of depleting and replenishing collisions in Eq. (4) into a single integral. Strictly speaking, Eq. (5) will only hold for elastic collisions, i.e. those where g = g ′ and no change in internal energy states occurs. This is the case for the N-N collisions at the top of Table I and also for the other two fast processes we have defined. Since N 2 (k)-N and N 2 (k)-N 2 (l) intra-bin scattering comprises all possible transitions between rovibrational levels within a given bin, they are not true elastic collisions. However, this energy is only exchanged within the bin, so we effectively treat them as if they were elastic collisions and in our coarse-grain model Eq. (5) is assumed to hold true for all fast collision types.
Although the differential cross section σ ij may in general depend on both g and ω, for the calculations discussed in Sec. V, we will neglect their dependence on the latter. This allows us to replace the differential cross sections in Eq. (4) with their integral counterparts 42 σ I ij (g) = S 2 σ ij (g, ω) dω and employ the variable hard sphere (VHS) model 43 for isotropic scattering in our DSMC calculations. As discussed in App. C, the choice of scattering model has a direct effect on the transport properties of the corresponding Navier-Stokes calculations. We should note that employing transport coefficients based on the VHS model in CFD calculations of viscous flows is rather unusual, since much more accurate methods are available 44,45 . In fact, several researchers have gone the opposite route 39,40,46 and "calibrated" the VHS, or similar cross sections in their DSMC codes with the state-of-the art transport collision integrals. In the present work, we base our transport properties on collision integrals derived from the VHS model (see App. C) to ensure consistency with our DSMC calculations, thus making the comparisons in Sec. V D more straightforward.

B. Slow collision operators for N and N 2 (k)
The slow collision operators account for all types of reactive collisions in the broader sense of our coarse-grain state-to-state description. The general mathematical form of reactive collision terms has been derived in Sec. 4.2.5 of Giovangigli 34 and here we merely write down the particular cases applicable to the slow processes listed in Table I. Operator C k (f ) appears in all rows of Eq. (3) involving the pseudo-species N 2 (k). It is itself composed of two separate terms, C k (f ) = C E k (f ) + C D k (f ). The first one accounts for the effect of excitation/deexcitation on f k 47 : Here, σ l,N k,N = σ l,N k,N (g, ω) is the differential cross section for the transition of an N 2 (k)+N pair into an N 2 (l)+N collision pair. The ratio of degeneracies a l /a k corresponding to postand pre-collision internal energy states N 2 (k) and N 2 (l) appears multiplying the post-collision distributions to account for detailed balance between forward (i.e. excitation) and the backward (i.e. deexcitation) reactions. For the excitation-deexcitation reaction, the detailed balance relation expressed in terms of forward and backward differential cross sections takes on the form: Again, we will assume isotropic scattering for all such collisions and replace the differential cross section with their counterparts integrated over all deflection angles: Here, σ Ef k→l (g) = S 2 σ l,N k,N (g, ω) dω is the integrated excitation cross section evaluated at the "pre-collision" relative speed g = |c k − c N | and σ Eb k→l (g ′ ) represents the integrated cross section for deexcitation from bin N 2 (l) to bin N 2 (k) evaluated at the "postcollision" relative speed g ′ = |c ′ l − c ′ N |. Energy conservation implies that the relation g ′ = g 2 + 2(E k − E l )/µ N 2 ,N must hold between pre-and post-collision pairs. Here, is the reduced mass for the N 2 -N collision pair. Notice also that the summation in Eq. (6) excludes the term (k = l), because this corresponds N 2 (k)-N intra-bin scattering, which we consider belonging to the fast processes.
The second term contributing to C k (f ) is due to dissociation-recombination reactions: This expression is more complex than Eqs. (4) and (6), because it involves a three-body interaction (the three N atoms after dissociation). This is reflected in the triple product of "post-collision" distribution functions f N appearing as part of the replenishing term in Eq. (9). Notice that instead of being "primed", these three f N are each identified by a unique overbar to distinguish them from one another. Equation (9) now involves a 12-fold integral in velocity space. The factor W 3N k,N is referred to 34,48 as the "reaction probability" for the dissociation-recombination reaction (in the forward sense), even though it has dimensions of time 8 × length −6 . Unlike in Eqs. (4) and (6), it is not straightforward to write Eq. (9) in terms of a differential, or integrated cross section.
The factors β k and β N , which appear in Eq. (9) multiplying the replenishing term are "statistical weights" of the colliding species: where h P is Planck's constant, m N 2 , m N are the molecular masses (m N 2 = 4.65 × 10 −26 kg for all N 2 (k) and m N = 1 2 m N 2 for atomic nitrogen) and a k , a N again the degeneracies of pseudospecies N 2 (k) and of N respectively. The ratio of statistical weights appears in Eq. (12) to account for detailed balance between the forward (i.e. dissociation) and backward (i.e. recombination) reactions. Analogous to the case for excitation-deexcitation just discussed, the terms in Eqs. (9) and (12) accounting for dissociation-recombination have been written exclusively in terms of the forward probability W 3N N,k , i.e. in the left-to-right sense as written in Table I). This is possible, because we have postulated the existence of a reversibility relation for this three-body interaction: The statistical weight of atomic nitrogen appears on both sides of Eq. (11) with an exponent equal to its stoichiometric coefficient right and left of Eq. (2), but simplifies once substituted into Eq. (9). Notice also that Eq. (11) implies that the dimensions of W k,N 3N are now just time 5 .
Finally, when considering the Boltzmann equation for atomic nitrogen, C N (f ) accounts for the effect of N + N 2 (k) dissociation-recombination on f N and assumes the form: Every element of the sum in Eq. (12) is composed of two integrals. Both share the same structure as the one in Eq. (9). The first one is focused on atomic nitrogen on the left of Eq.
(2) and accounts for depletion of this species due to dissociation and its simultaneous replenishment due to recombination. The second integral does the same, but is focused on one of the three N-atoms on the right of Eq. (2). It accounts for depletion of any of the three N-atoms due to recombination and their simultaneous replenishment due to dissociation, hence the minus sign multiplying the integral. The factor 3 appears, because one must account cumulatively for the loss of the three nitrogen atoms on the right-hand side of Eq. (2).
Writing down Eq. (3) and the associated collision terms is a useful framework for deriving the macroscopic equations in Sec. IV. However, in this work we only solve the Boltzmann equation indirectly, by means of the particle-based DSMC method. In this approach the behavior of the collision terms has to be translated into a collision algorithm, which has been detailed previously in Ref. 29 .

C. Macroscopic flow variables in terms of velocity distributions
The set of kinetic equations represented by Eq. The mass density of every species is given by: with individual species number densities following from n i = ρ i /m i . Mixture number and mass densities are calculated as n = i∈S {n i } and ρ = i∈S {ρ i } respectively. The hydrodynamic velocity of the gas is given by: Diffusion velocities of each species are given by: where C i = c i − u represent the peculiar velocities of particles belonging to species i ∈ S.
By definition, the diffusion velocities always verify the constraint The kinetic stress tensor is obtained as: The pressure tensor can be split into an isotropic and a remaining anisotropic contribution P = p I − τ , where p is the hydrostatic pressure, I stands for the unit tensor and τ is the viscous stress tensor. The hydrostatic pressure is calculated as 1/3 of the trace of P, e.g.
in Cartesian coordinates p = 1 3 (P xx + P yy + P zz ). Next, we invoke the perfect gas law to introduce the translation temperature as T = p/(nk B ). Since we are dealing with a dilute gas mixture, we may express the composition in terms of partial pressures p i = x i p, where x i = n i /n are the species mole fractions. Alternatively, the mixture composition can be expressed in terms of mass fractions It is an implicit function of the number densities n k , as explained Appendix C of Ref. 29 .
The total energy per unit volume in terms of the distribution is given by: where the E i represent the internal energies of each species i ∈ S. In our coarse-grained state-to-state description, they correspond to the bin-averaged energies E k for each internal state N 2 (k), ∀ k ∈ K N 2 and E N to the 0-K energy of formation of atomic nitrogen. For consistency with our prior definitions 29,31,38 , we set E N = D 0 /2, where D 0 = 9.75 eV is the heat of dissociation per N 2 -molecule from the ground rovibrational level as given by the NASA Ames N3 diatomic potential 49 . Notice that the kinetic temperature T and Eq. (17) are related to one another through ρE = 1 Finally, the mixture heat flux is the flux of kinetic and internal energy transported with every particle along each Cartesian direction: For the exact expressions used to evaluate Eqs. (13)-(18) in our DSMC calculations, refer to App. A.

IV. HYDRODYNAMIC DESCRIPTION FOR COARSE-GRAIN MODEL
In this section we discuss the macroscopic balance equations used to model the flow at the hydrodynamic scale. They are derived from Eq. (3) by applying the Chapman-Enskog method. 34,50,51 Here we give a quick overview of this procedure for our particular application.

A. Chapman-Enskog method for coarse-grain model
We introduce suitable reference quantities at the kinetic and macroscopic level to perform a dimensional order-of-magnitude analysis 52 of Eq. (3). This allows us to re-write it in its non-dimensional form:D where Kn = λ 0 /L 0 is a pseudo-Knudsen number based on reference mean free path λ 0 and macroscopic length scale L 0 . The scaling for arriving at the compressible Navier-Stokes equations is to select Kn ∼ ε ≪ 1. The fast and slow processes in Eq. (19) are assumed to occur at time scales different enough to require separate reference cross sections and the Maxwellian reaction regime 34 is obtained assuming that σ slow ∼ ε 2 σ fast . Applying this scaling is a choice, which ultimately determines the structure of the resulting hydrodynamic equations. Expressed in terms of the small parameter ε and reverting back to dimensional variables for convenience, we will thus seek solutions to Eq. (3) in the continuum limit of the form: Performing an Enskog expansion around the zero-order velocity distributions f 0 i in terms of the small parameter ε: i are perturbation functions, and substituting back into Eq. (20) yields: where Solving Eq. (21) at order ε −1 (corresponding to the fastest time scale) yields the equilibrium, or Maxwell-Boltzmann distribution. Here we have defined the macroscopic moments: species mass density The equality f 0′ i f 0′ j = f 0 i f 0 j allows us now to write the linearized fast operator as: Averaging Eq. (21) at order ε 0 over pseudo-species mass, momentum and energy leads to the Euler equations for the non-reacting gas mixture: Notice that due to the choice of constraints, the definitions of macroscopic moments in the Chapman-Enskog solution slightly differ from those introduced in Sec. III C at the kinetic scale. However, out of convenience here we will use the same symbols for both definitions.
Definitions for ρ N 2 , etc. and corresponding number densities follow the same pattern as in Sec. III C. Note also that, given the scaling in Eq. (20), the slow collision operators do not contribute to the solution at order ε 0 , and thus no chemical source terms appear on the right hand side of Eq. (23).
B. Macroscopic balance (Navier-Stokes) equations for coarse-grain system including viscous and chemical source terms With f 0 i known, we go back to solving Eq. (21) at order ε 0 for the first-order perturbations φ = (φ i ) i∈S : Uniqueness of the solution is ensured through the constraint that the perturbations do not contribute to the macroscopic moments, i.e.: (22), we express all resulting time derivatives of macroscopic flow variables in terms of spatial gradients by re-arranging Eqs. (23)- (25). The result is a linear combination of the transport forces 34 , i.e. gradients in flow velocity, species partial pressure and temperature: where we have defined the driving forces for diffusion of species j as: d j = (∇ x p j )/p. The individual contributions to Ψ i are given by: It can be shown that the unique solution to Eq. (26) is a linear combination of the transport fluxes: The tensorial functions φ η = (φ η i ) i∈S and vectorial functions In the continuum, or hydrodynamic limit the complete governing equations are finally obtained by averaging Eq. (21) at order ε 1 over pseudo-species mass, total momentum and energy: Here, Eq. (33) represents the set of continuity equations for every pseudo-species i ∈ S.
The structure of the chemical source terms on the right hand side is discussed in more detail in Sec. IV D. The transport fluxes for pseudo-species mass, momentum and energy appearing in Eqs. (33)- (35) are given in the Chapman-Enskog approximation by: respectively. We discuss the manner in which these fluxes are evaluated in Sec. IV C.

C. Transport fluxes
Solving the kinetic equations (26) subject to the constraint j∈S {j j } = 0 to ensure mass conservation. The Stefan-Maxwell matrix entries in Eq. (39) are: In the absence of external force fields, all remaining driving forces for diffusion appear on the right hand side of Eq. (39). The previously introduced linearly dependent driving forces can be decomposed as In order to evaluate the entries of the Stefan-Maxwell matrix, one must supply the bi- The expression for χ i = χ i (p j ∀(j ∈ S), T ) is given in Chapter 5 of Ref. 34 and in App. B.
Following the structure of the matrix for the thermal conductivity transport system, it can be shown that their sign is not defined, but that i∈S {χ i } = 0 must hold 53 . Alternatively, the diffusion fluxes can be expressed in terms of multi-component diffusion coefficients , and is the pseudo-inverse of the Stefan-Maxwell matrix appearing in Eq. (39).
The viscous stress tensor τ appearing in Eqs. (34) and (35) takes on the form: where η is the mixture shear viscosity and S = 1 The first term on the right hand side is the contribution due to heat conduction q cond . It is the product of the mixture thermal conductivity λ and the temperature gradient. The second term q diff accounts for heat transfer by diffusion of enthalpy of each mixture component, The expression for the thermal conductivity λ = λ (p i ∀(i ∈ S), T ) is given in Chapter 5 of Ref. 34 and in App. B. An alternative formulation for the heat flux is to use the partial thermal conductivity λ and the thermal diffusion coefficients θ i , (i ∈ S).
Both formulations are equivalent, but the one chosen here is advantageous to study the entropy production in Sec. IV E. Following the structure of the matrix for the thermal conductivity transport system, it can be shown that λ > 0, provided that some conditions on the collision integral data are met 53 . Note that within the state-to-state formalism there is no need to consider Eucken's correction to the thermal conductivity 50 , because transfer of internal energy is implicitly taken into account through the chemical production rates ω i , i ∈ S. The third term formally accounts for heat transfer induced by concentration gradients (Dufour effect). It is the complement to the Soret effect appearing in Eq. (39).
However, note that it is also being neglected in the Navier-Stokes calculations presented in The necessary routines for the solution of the transport systems have been implemented in the Mutation++ 54 thermodynamic and transport library, which is tightly coupled to the Navier-Stokes flow solver used to generate the results of Sec. V B and V D.

D. Chemistry source terms
The terms on the right hand side of Eq. (33) represent the mass production terms for atomic nitrogen and every N 2 (k) respectively. For the latter, both excitation-deexcitation and dissociation-recombination reactions contribute to the source term: ω k = ω E k + ω D k . These two contributions are obtained by averaging Eqs. (6) and (9) Normalized with the respective molecular masses, these terms take on the following form: for excitation-deexcitation and: for dissociation-recombination. For atomic nitrogen only the dissociation-recombination reactions contribute to the source term. Taking the moments of Eq. (12) in analogous manner yields ω N = m N C N (f 0 ) dc N , and can be simplified to: In previous work 38 the coarse-grain reaction cross sections σ Ef k→l (g) and σ Df k (g) were fitted to an analytical form consistent with Arrhenius-type expressions for the corresponding rate coefficients k Ef k→l (T ) and k Df k (T ) appearing in Eqs. (44)- (46). Special care was taken to ensure consistency between the kinetic and hydrodynamic description. This meant that the reversibility relations postulated to exist between forward and backward cross sections/probabilities as discussed in Sec. III whereas the recombination rate coefficient appearing in Eqs. (45) and (46) is obtained as: Here, the partition function per unit volume of each pseudo-species i has the form:

E. Entropy equation and sign of the chemical entropy production term
A macroscopic balance equation for the entropy per unit volume based on thermodynamic considerations is derived in Chapter 2.6 of Giovangigli 34 . In the form applicable to our case it reads: where terms on the left-hand side represent the (1) local time rate of change of entropy, On the right hand side of Eq. (49) the volumetric entropy production rate can be split up into Υ = Υ tran + Υ chem , i.e. entropy production due to (a) transport phenomena and (b) chemical reactions 55 . General expressions for both terms have been derived by Giovangigli,34 and here we recall only the terms relevant for our fluid model. The first production term can be written as: Given the structure of the first two terms on the right hand side of Eq. (50) and the fact that η, λ > 0 it can be easily seen that they must always be non-negative. The third term contains as factors D ij the components of the multi-component diffusion matrix. Its properties guarantee that the associated entropy production term will always remain nonnegative. Thus, Υ tran ≥ 0 must hold for any physically realizable flow.
Now, for the particular set of reactions given by Eqs. (1) and (2), it is worthwhile to have a closer look at the entropy production due to chemical reactions: It is a function on the Gibbs free energies per unit mass and the chemical source terms appearing on the right hand side of Eq. (33).
Following the general procedure outlined by Giovangigli 34 , it is possible to show that Υ chem ≥ 0 for all cases, in accordance with the second law of thermodynamics. The key to demonstrating this lies in re-writing the chemical production terms in the symmetric form (Sec. 4.6.6 of Giovangigli 34 ), where the rate coefficients for the excitation-deexcitation and dissociation-recombination reaction become k s Consistency between these production rates in symmetric form and the original notation of Eqs.  1) and (2) verifying detailed balance. This, in turn, implies that the backward rate coefficients for excitation-deexcitation and dissociation-recombination must be computed according to Eqs. (47) and (48) respectively. After some algebraic manipulation, one arrives at the final form: for the entropy production due to chemical reactions. Here, we have defined the re-

V. INTERNAL ENERGY EXCITATION AND DISSOCIATION ACROSS NORMAL SHOCK WAVE
In this section we present simulation results for a steady, normal shock wave. We apply three distinct numerical approaches and compare them in terms of their degree of physical fidelity. In order to formulate a discretized version of the macroscopic balance equations amenable to numerical solution, we re-write Eqs. (33)- (35) for the unsteady, one-dimensional case in the form: where U = (ρ i (i ∈ S), ρu x , ρE) T is the vector of conserved variables, F = (ρ i u x (i ∈ S), ρu 2 x + p, ρu x (E + is the inviscid flux vector and

A. Hydrodynamic inviscid solution based on Finite difference ODE method
To obtain a first estimate of the thermo-chemical non-equilibrium region, we simulate the normal shock following a steady-state, one-dimensional inviscid approach. When such conditions are assumed, the time derivatives ∂U/∂t and the diffusive transport fluxes F d in Eq. (52) all vanish. This makes it possible to re-cast the original set of equations into an ordinary differential equation (ODE) system: where the solution vector is now given by P Two different supersonic free-stream conditions are considered. For the high-speed case, we impose a free-stream velocity of u 1 = 10 km · s −1 , while for the low-speed case we use u 1 = 7 km · s −1 . All other parameters, such as free-stream temperature, pressure and composition are the same for both cases. The higher-speed conditions are listed in Table II, where they are labeled as (1) pre-shock. In the ODE approach the shock is not captured by the numerical method. It is instead replaced by a sudden jump in flow conditions at x = 0, which only affects the translational mode. Therefore, the analytical Rankine-Hugoniot jump relations with specific heat ratio γ = 5/3 are used to predict the non-equilibrium post-shock state (state (1a) in Table II). While the kinetic temperature reaches T = 62550 K behind the discontinuity, the internal temperature and composition remain frozen at the free-stream values. Thus, the initial bin mass fractions y k , (k ∈ K N 2 ) in Eq. (53) are made to follow a Boltzmann distribution at T int = 300 K. The ODE algorithm then marches along x starting from state (1a). Notice that the free stream contains a non-zero amount of atomic nitrogen, even though the gas in equilibrium at 300 K should only consist of N 2 -molecules.
We add a small amount of N to the free-stream gas to trigger internal energy exchange and dissociation processes, since only reactions induced by N-N 2 collisions are taken into account by the chemical source terms of Eqs. (44)- (46). The pre-and post-shock conditions for the low-speed case are listed in Table III. Due to the lower post-shock temperature, the gas does not dissociate to the same degree as at the high-speed conditions and about 1/3 of the post-shock gas remains in the form of molecular nitrogen. We carry out four separate simulations at the high-and low-speed conditions respectively. The first simulations provide  Tables II and III,  We obtain such close agreement with the full model, because we are using energy bins of variable, instead of constant width. In previous work 29 we were able to show that switching to variably-sized bins allows us to closely match the thermodynamic properties of the full model with a much smaller number of URVC bins. In particular, using more bins of smaller width to group together the lowest-energy rovibrational levels is advantageous to accurately capture the internal energy content of the cold free stream.
Mass density and temperature profiles for the high-and low-speed cases are plotted in   Tables II and III,  In Fig. 1a we plot profiles of mixture density ρ (continuous lines) and molecular nitrogen ρ N 2 (dotted lines) for the high-speed case. The behavior in all four cases is very similar and the main differences are confined to the region immediately behind the shock front.
Each one of the four ρ N 2 -profiles reaches its maximum several millimeter downstream of the discontinuity, before dissociation begins to consume the remaining molecular nitrogen.
The reference solution for the full N3 system exhibits the quickest response to the shock, whereas the coarse-grained systems lag behind. The response becomes slower with decreasing number of bins. In Fig. 1b we examine the corresponding temperature profiles. The kinetic temperature T quickly decreases from its initial value of 62550 K to about 30000 K in the first 2 − 3 mm behind the discontinuity. Simultaneously, the internal temperature rises from its free-stream value of 300 K to a maximum of about 25000 K in the same distance, In Fig. 2a we now show the density profiles for the low-speed case. Again, all four systems follow the same general behavior. Whereas in the high-speed case practically all molecular nitrogen eventually dissociated behind of the shock front, at these lower-speed conditions the N 2 -profiles remain fairly flat further downstream. However, the trend is now reversed, in the sense that the 10-bin system is the quickest to react to the shock, whereas the response becomes slower as the number of bins is increased all the way up to the full system. Figure 2b shows the corresponding temperatures for the low-speed case. With a length of approximately 5 m, the post-shock non-equilibrium region is now almost two orders of magnitude greater than in Fig. 1b. A closer look suggests that at these lower- in the opposite sense relative to the high-speed case.
By studying these two flow conditions with the inviscid ODE method we found that the relaxation region for the high-speed case extends for about 10 cm and for the lowspeed case roughly 5 m from the discontinuity. This helps us size the domain and to adjust the computational parameters for the Navier-Stokes and DSMC calculations discussed in Secs.V B and V C. Furthermore, we see that the coarse-grained model has an influence on the evolution of the gas state in the post-shock region and these profiles diverge to some degree from the reference solution. As would be expected, the closest agreement with the full system is observed for the cases with the largest number of bins (837), while the biggest differences are observed for the 10-bin cases. However, these deviations become less severe further downstream of the initial discontinuity. discontinuity in flow parameters. The portion of the flow field left of the discontinuity is initialized to the pre-shock equilibrium state, whereas to its right the post-shock equilibrium state is imposed (recall Tables II and III for

Navier-Stokes (black lines) and inviscid ODE flow field of Sec. V A (red lines). All profiles
shown are for the high-speed condition using the 10-bin coarse-grained system. Density profiles are shown first in Fig. 3. The origin of the x-axis lies at the location of the initial discontinuity for the Euler cases. Due to numerical diffusion in the FV approach this discontinuity is captured over an extent of 2-3 cells (see close-up in Fig. 3(b)). However, the grid has been carefully refined in the vicinity to ensure that this adverse numerical effect remains minimal. This is confirmed by the excellent agreement of the FV-Euler and inviscid ODE density profiles over the remainder of Fig. 3 Fig. 4(b)). The jump in kinetic temperature is captured well by the FV method, as is its peak value for the inviscid case. Again, viscous effects act to smooth out these flow features and diffuse the shock front upstream. In the Navier-Stokes profile the gas temperature begins to depart from its pre-shock value about 0.003 m ahead of the initial discontinuity and reaches a lower maximum (T max ≈ 51800 K for Navier-Stokes vs.

C. Normal shock solution with DSMC
In this section we describe how the normal shock for both the high-and low-speed conditions was simulated using the DSMC method 3  To ensure that the shock front builds up at a well-defined location within the domain, we generate initial particles corresponding to the pre-shock equilibrium state (1) in the region left of the initial discontinuity and particles corresponding to post-shock equilibrium state (2) to the right of this location. This becomes the point where the supersonic free stream is "tripped" into transitioning to the post-shock equilibrium state and marks the initial location of the standing shock. As the simulation progresses, this discontinuity is smoothed out by particle transport. Once this phase is complete, the steady-state flow parameters are gathered from the DSMC particles and further refined through time-averaging. The location of the initial discontinuity is somewhat arbitrary, but if it is placed too close to either boundary, random walk may push the shock front out of the domain before steadystate macro-parameters can be extracted. Given that our primary goal is to observe as much of the relaxation region behind the shock, we place it as close as is reasonable to the left boundary. By setting L u = 3 cm (see Table IV), we make sure to leave ample space (i.e. 6 000 cells) between the inlet and the location of the initial discontinuity. Notice that for the high-speed condition only parameters for the 10-bin and 100-bin systems are listed in Table IV. Due to the greater computational cost of the DSMC method compared to the ODE approach of Sec. V A and the Navier-Stokes calculations of Sec. V B, no DSMC simulations for the higher-resolution 837-bin case and the full database were carried out.
In both high-speed simulations, the DSMC particle weight is set to ensure that at least 20 particles are present in every upstream cell. Due to the rise in density across the shock, there are ≈ 540 particles per cell in the downstream region. For the 100-bin case the domain length is reduced to L u = 3 cm and L d = 10 cm respectively. This reduction is justified, as we are still able to capture the full relaxation region, while significantly reducing the computational expense.
Two complementary measures are taken to reduce the statistical noise inherent in DSMC flow fields. For the two high-speed cases in Table IV we  While for the high-speed case we could comfortably contain the entire shock within 40 000 collision cells, this was not feasible for the low-speed case. Fortunately, for our purposes it is not necessary to simulate the entire post-shock relaxation region with DSMC. As was seen for the high-speed case, most of the diffusive effects are only appreciable within a narrow region surrounding the shock front. By concentrating on this portion we managed to significantly reduce the domain size. To accomplish this, we resort to the approach described by Strand and Goldstein 62 , where the normal shock is treated as inherently unsteady. The supersonic free stream is fed into the domain on the left boundary, while a specular wall reflects all particles on the boundary to the right. This stagnates the incoming flow and generates a shock wave moving from right to left into the undisturbed gas upstream. Unlike in the previous set-up, the reference frame is now attached to the post-shock equilibrium gas, implying that u ′ 2 = 0. Therefore, in order to obtain the desired post-shock thermodynamic conditions of Table III in our simulation, we adjust the inflow velocity to u ′ 1 = u 1 − u 2 . Once the shock front has left the near-wall region, it begins to take on its steady-state structure and travels upstream at approximately u shock = −u 2 . At this point macroscopic flow parameters can be sampled and individual samples time-averaged to reduce statistical noise. Since the shock is continuously moving upstream, these instantaneous samples have to be displaced to a common origin before time-averaging. Again, we resort to the procedure described in Ref. 62 to define a common reference location for all profiles.
At the low-speed condition the higher post-shock density and lower temperature (see Table III) impose more stringent constraints on the collision cell-and time step size. Thus, in the rightmost column of Table IV, several simulation parameters were adjusted accordingly.
Just as for the high-speed condition, ensemble-averaging over 64 independent simulations is used to reduce the statistical scatter in the instantaneous samples.

D. Comparison Navier Stokes vs. DSMC
We now examine the flow fields obtained through the methods described in Secs. V B and V C. First, in Figs. 5 and 6 we compare DSMC profiles obtained using the 100-bin (blue dot on blue line) and 10-bin (red line) systems to Navier-Stokes profiles with the 10-bin system (black square on black line) at the high-speed conditions.
We start with the gas density profiles in Fig. 5a. The DSMC and Navier-Stokes curves have been translated on the x-axis, such that the initial rise in density occurs at the same location for all three profiles. The location of the origin is arbitrary, but the same convention is used consistently in all flow parameter plots in Figs. 5 and 6. Focusing on the 10-bin system, both the DSMC and Navier-Stokes density profiles show close agreement, except for a weak increase of the density slope in the DSMC result at x ≈ 0, which is absent from the Navier-Stokes curve. The Navier-Stokes density profile exhibits a quicker and more uniform initial rise, before intersecting the DSMC profile at x ≈ 0.003 m.
Next, in Fig. 5b we compare the corresponding kinetic and internal temperatures. Here, the differences between both methods are more apparent. The maximum T -value obtained with DSMC (10 bins) is T max ≈ 58800 K, which lies roughly 7000 K above the corresponding peak for Navier-Stokes. Incidentally, both maxima lie very close to one another, at x ≈ −0.002 m. By contrast, the maximum T int -values for all three curves are much closer to one another, with the Navier-Stokes curve slightly leading the DSMC profiles. The most noticeable difference is that both T -curves for DSMC begin to rise farther upstream and more gradually than the Navier-Stokes profile. Back in Fig. 5a, we also plot the partial density of N 2 using dotted lines. As was observed in Fig. 1a for the inviscid case, there is an initial rise in ρ N 2 across the shock, before dissociation kicks in and gradually consumes the molecular nitrogen further downstream. At these high temperatures, the post-shock gas is almost entirely made up of atoms. Here, Navier-Stokes predicts dissociation occurring slightly ahead of the corresponding DSMC (10 bins) curve. This is consistent with the lower kinetic temperature observed for Navier-Stokes in Fig. 5b. We now move on to Fig. 6 and the comparison of flow parameters associated with diffusive transport at the high-speed condition. In Fig. 6a we first show the mass diffusion flux of N 2 along the x-direction. For the two DSMC curves and the single Navier-Stokes result j x,N 2 is calculated as the mass-weighted average over all internal energy bins, i.e.: j x, The corresponding mass diffusion flux of atomic nitrogen: j x,N = ρ N u d N (not shown) is equal in magnitude, but opposite in sign. The peak of j x,N 2 captured by the DSMC and Navier-Stokes methods with 10 bins agrees to within less than 5%, although in the Navier-Stokes profile this maximum appears slightly ahead of the DSMC curve. For the 100-bin DSMC case, the peak diffusion flux lies about 10% below the corresponding 10-bin DSMC value, but at almost exactly the same x-location. Next, in Fig. 6b we plot the three normal components of the viscous stress tensor. For our one-dimensional flow configuration only the velocity derivative ∂u x /∂x becomes non-zero across the shock. As a consequence, the only components of τ in Eq. (34), which take on non-zero values turn out to be τ xx = 4 3 η (∂u x /∂x) and τ yy = τ zz = − 2 3 η (∂u x /∂x). Both DSMC and the Navier-Stokes profiles reach their maxima at essentially the same x-location. The DSMC stress profiles are slightly more spread out than their Navier-Stokes counterparts. The ratio τ xx,max /τ yy,max yields exactly −2 for the Navier-Stokes profiles, in accordance with the analytical expressions for τ xx and τ yy . The same ratio of −2 is maintained for the DSMC profiles, although the peak viscous stresses obtained with Navier-Stokes lie about 34% above the corresponding DSMC values. As can be seen by comparing the two DSMC profiles, the number of bins has practically no effect on the shape of the viscous stress profiles. Finally, in Fig. 6c we compare q x , i.e. the heat flux component along the flow direction. Both DSMC and the Navier-Stokes profiles exhibit their peak negative values (due to heat being transferred upstream across the shock front) at roughly the same x-location. However, the maximum flux for DSMC is nearly −22.1 MW/m 2 , while for the Navier-Stokes result it only reaches −16.7 MW/m 2 . As was the case for the kinetic temperature in Fig. 5b, the DSMC heat flux profiles are noticeably more diffuse and begin to deviate from zero much sooner upstream than their Navier-Stokes counterpart. A second smaller, but positive peak appears in all three q x -profiles further downstream. Thus, some amount of heat is also being transferred from the shock front in the downstream direction.
It is interesting to note that the location of this second, positive peak in q x nearly coincides with the maximum in j x,N 2 reported in Fig. 6a for all three calculations. One might thus assume that "diffusion of enthalpy" plays a significant role in shaping the heat flux profile in this region. In order to answer this question we have decomposed the Navier-Stokes (solid black lines) result into q cond x , i.e. its contributions due to heat conduction (dash-dotted line) and q diff x , i.e. its contribution due to diffusion of enthalpy (dotted line). It turns out that the second peak observed in the q The overall close agreement between the DSMC and Navier-Stokes profiles in Figs. 5 and 6 is somewhat surprising. Given the strong deceleration, the molecular velocity distributions across the shock obtained with DSMC will deviate significantly from the Chapman-Enskog distribution, on which the Navier-Stokes solution is based. Thus, one might have expected a greater difference between both results. Another noteworthy aspect is that, apart from minor differences in the mixture and partial density profiles, the 10-bin and 100-bin DSMC flow fields exhibit almost the same behavior. This is in contrast with what was observed in system. Recall from Sec. V A that at 7 km · s −1 the post-shock chemical nonequilibrium region extends much farther downstream than at 10 km · s −1 . However, here we focus on the region immediately surrounding the shock front, where the strongest thermo-chemical nonequilibrium is observed. Thus, density, temperature and in particular mixture composition do not fully reach their post-shock equilibrium values in the x-range plotted. However, the moments associated with viscous and diffusive phenomena adjust much more quickly and are fully contained within the region shown.
In Fig. 7a we begin by plotting density profiles. As was done for the high-speed case, the DSMC and Navier-Stokes profiles have been aligned such that the initial rise in density occurs at a common x-location. For both the DSMC and Navier-Stokes calculations the overall gas density ρ is represented by solid lines, whereas ρ N 2 is shown using dotted lines.
One can see two distinct "bumps" in both ρ-profiles, with the first one appearing at the same x-location with both methods. Near the second bump further downstream, the two ρcurves begin to diverge, and beyond this point the DSMC profile remains slightly above the corresponding Navier-Stokes curve. Up until the second bump in the ρ-profiles dissociation plays only a minor role. But past this point the amount of atomic nitrogen begins to rapidly increase, while ρ N 2 remains almost constant. In Fig. 7b we plot the corresponding temperature profiles. As was seen for the high-speed case in Fig. 5b, the peaks in kinetic temperature T appear at almost the same x-location for both DSMC and Navier-Stokes.
Of course, given the significantly lower total enthalpy of the flow, the peak T -values are much lower than for the high-speed case. At T max ≈ 31100 K, DSMC predicts a somewhat higher peak value than Navier-Stokes, where a maximum of ≈ 28200 K is reached. Similar to the high-speed case, the kinetic temperature profile from DSMC in Fig. 7b is more diffuse and exhibits a more gradual initial rise than the Navier-Stokes curve. The location of the T int -maximum appears almost exactly at the same x-location and both values differ by less than 2% (DSMC: T int ≈ 15300 K vs. Navier-Stokes: T int ≈ 15100 K). Slightly different behavior is seen downstream of this point, with the common DSMC temperature decreasing somewhat faster than in the Navier-Stokes profile. It is worth noting that both methods predict the highest kinetic temperature about 0.005 m upstream of the point where significant amounts of N-atoms begin to be produced. In fact, for both methods the location in Fig. 7a where the ρ and ρ N 2 profiles begin to diverge coincides with the peak in T int observed in Fig. 7b, and beyond which the translational and internal temperatures reach a common value. This suggests that at these lower-speed conditions a noticeable "incubation length" for dissociation exists and that dissociation primarily occurs under near-equilibrium conditions downstream of the shock front. Overall, DSMC predicts slightly quicker dissociation of N 2 than the Navier-Stokes calculation. This can be seen by comparing the density profiles in Fig. 7a. The behavior of the temperature profiles in Fig. 7b is consistent with this fact.
Since in the DSMC calculation a slightly larger number of endothermic dissociation reactions remove a greater amount of energy from the translational and internal modes, the DSMC temperature stays below the Navier-Stokes profile past the initial shock front.
Next, in Fig. 8 we compare the flow parameters associated with diffusive transport for the low-speed case. First, in Fig. 8a we examine the diffusion fluxes of N 2 along the xdirection. Here, slightly different behavior between DSMC and the Navier-Stokes profiles are apparent. The diffusion flux for N 2 obtained with DSMC exhibits two distinct peaks, one at x ≈ −0.001 m and another closer to x = 0.0075 m. This behavior is exactly mirrored for N, although with opposite sign (not shown). By contrast, in the Navier-Stokes solution the first peak does not appear at all. Furthermore, the maxima in predicted j x,N 2 lie at about 0.0065 kg · m/s for DSMC vs. 0.005 kg · m/s for Navier-Stokes.
In Fig. 8b we plot the three normal components of the viscous stress tensor for the low-speed case. The magnitudes of these stresses are approximately half of those for the high-speed case, but follow the same general behavior. Both for DSMC and Navier-Stokes we retrieve precisely τ xx,max /τ yy,max = −2, but the ratio between the peak values is now [τ xx,max ] NS /[τ xx,max ] DSMC = 1.23. In a slight departure from the high-speed case, the normal stresses do not immediately return to zero downstream of their peaks. Instead, a small plateau forms in both the DSMC and Navier-Stokes profiles.
Finally, in Fig. 8c we compare the heat flux profiles for the low-speed shock. The peak heat flux for DSMC was observed to be −8.06 MW/m 2 , whereas it was −5.70 MW/m 2 in the Navier-Stokes result. This amounts to a ratio [q max ] NS /[q max ] DSMC = 0.708, as opposed to 0.756 for the high-speed case. As was the case for the high-speed case, the DSMC and Navier-Stokes profiles agree in general shape, but differ somewhat in the location and magnitude of their maxima. As had been observed for the high-speed case, the initial departure from zero begins further upstream and is more gradual in DSMC than in the Navier-Stokes profile.
Past the initial negative peak in q x , both profiles exhibit a second, slightly positive overshoot downstream of the shock front. This peak, or plateau is much less pronounced and more spread out than in the high-speed case. Again, in Fig. 8c we have split up the Navier-Stokes profile into contributions due to heat conduction (dash-dotted line) and diffusion of enthalpy  fluxes are less significant when compared to the high-speed case, but the general effect is still present at this condition.

VI. CONCLUSIONS
We have presented the procedure to build a coarse-grain fluid model incorporating internal energy exchange and nonequilibrium chemistry fully consistent with the gas-kinetic description. The resulting hydrodynamic equations are equipped with dissipative transport and chemical source terms that are rigorously derived from the collision operators of the underlying kinetic equation.
We have used a state-to-state approach, which allows for detailed description of inelastic processes in a gas mixture. A set of coarse-grain cross sections and corresponding rate coefficients derived from the NASA Ames ab initio database for the N 2 (v, J)-N system was Our calculations reveal rather close agreement between the Navier-Stokes and DMSC predictions. This is somewhat surprising given the free-stream Mach numbers we studied (Ma ∞ ≈ 28 and 20 respectively). At such extreme conditions one could have expected that the inability of the Navier-Stokes solutions to fully reproduce the strong translational nonequilibrium effects across the shock (i.e. bi-modal velocity distributions) would cause them to more noticeably deviate from the DSMC results. However, even though the macroscopic flow properties predicted with DSMC are more diffuse than in the Navier-Stokes calculations, all major features appear in both solutions at nearly the same x-location and are comparable in magnitude. With regard to resolving the shock structure with Navier-Stokes, our study reinforces the notion that employing transport properties consistent with the corresponding scattering cross sections is fundamental to obtaining close agreement with kinetic-scale solutions. Furthermore, our DMSC calculations suggest that the sensitivity of the flow field to the number of bins used in our coarse-grain model is greatly attenuated when viscous and diffusive transport effects are included.
It should be recalled that the URVC bin model we employ assumes a constant average energy for all rovibrational levels and freezes their relative populations within a bin. This is rather restrictive and clearly not the ideal reduction strategy. On the other hand, these constraints make deriving the associated fluid equations rather simple from a mathematical viewpoint. At this stage it is not clear whether equivalent asymptotic solutions can be derived in the same manner for other existing coarse-grain models. In particular, if each bin is assumed to have an associated temperature (e.g. Boltzmann bins), it is not straightforward to translate the model into the Chapman-Enskog framework, because of the need for compatibility of this ansatz with associated scaling of the collision operators. A rigorous treatment of the internal energy for multi-temperature gases is still an open problem for the Chapman-Enskog method. Other types of closures for transport phenomena, such as as the Maximum Entropy closure, 63,64 may be more natural.

ACKNOWLEDGMENTS
The authors would like to thank Dr. Federico Bariselli for his contributions to the improvement of the URVC bin model and to Dr. Alessandro Munafò for the codes used to generate the hydrodynamic solutions. We would also like to thank Dr. R. L. Jaffe and Dr. D. W. Schwenke from NASA Ames Research Center for access to the kinetic database used in this work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request is obtained through a classical discretization in phase space, such as by Muafò et al. 65 , a discretized version of the distribution function for each mixture component is obtained at every x-location. This can then be numerically integrated over velocity space to yield the macroscopic moments. By contrast, in the DSMC method the distribution function is implicitly represented by a finite number of simulated particles taking on random velocities, while physical space is discretized into an array of contiguous cells. Therefore, one directly estimates the local macroscopic moments by averaging over the ensemble of DSMC particles in each cell. In this section we give the equivalent expressions for macroscopic flow field variables resulting from the DSMC simulations discussed in Secs. V C and V D.
The mass densities defined by Eq. (13) in a given DSMC cell are calculated as: where W p = N real /N sim is the "particle weight" relating real gas molecules to DSMC simu-  Appendix B: Transport linear systems Linear systems of size N s must be solved to compute the mixture transport coefficients.
For shear viscosity one must solve j∈S {G η ij α η j } = x i , ∀ (i ∈ S), with the viscosity matrix

Appendix C: Collision integrals for viscous transport properties
Our goal is to ensure consistency between the transport phenomena modeled at the hydrodynamic scale of Navier-Stokes and the kinetic-scale DSMC simulations. To this end one must compute the relevant transport properties using collision integrals consistent with the set of cross sections and scattering laws used in DSMC.
For sake of simplicity, in this work we assume that all (pseudo-) species involved in the fast processes of Table I scatter isotropically. We use the variable hard sphere (VHS) model of Bird 3,43 for N-N elastic scattering and N 2 (k)-N 2 (l) intra-bin scattering. For the VHS model the differential cross section takes on the form: Stephani et al. 39 and listed in Table V. Furthermore, in Eq. (C1) k B is Boltzmann's constant, µ ij = m i m j /(m i + m j ) is the reduced mass for the species pair ij and Γ (.) is the gamma function. Since the VHS model assumes isotropic scattering, its differential cross section is actually independent of the post-collision deflection angle χ. The corresponding integrated cross section is obtained as σ I ij = 2π π 0 σ ij (g, χ) sin χ dχ: The integrated cross section is required at the moment of computing the collision probability [σ I ij (g) · g] pair /[σ · g] max for a given collision pair in the No Time Counter scheme 66 of DSMC. Note that in using a single set of parameters for all N 2 (k)-N 2 (l) collision pairs, we have implicitly assumed that the cross sections for all molecule-molecule intra-bin collisions possess the same value, regardless of the pre-collision internal states N 2 (k) and N 2 (l). − ω ij f ij (T ) (C7) with the common factor f ij (T ) = πd 2 ref,ij (T /T ref ) 1/2−ω ij . For the fast N 2 (k)-N intra-bin collisions, we also assume isotropic scattering. However, instead of defining the cross sections in terms of VHS parameters, we determine them directly based on the coarse-grained cross section database of Ref. 38 . An analytical expression for the integrated cross section σ I N 2 (k),N (g) = σ E k→k (g) was proposed in that reference and can be written as: where the notation σ E k→k is shorthand for the intra-bin scattering cross section of collision the pair N 2 (k) + N and A E k→k , b E k→k are parameters derived in Ref. 38 by post-processing the NASA Ames N3 database. Given the assumption of isotropic scattering, the corresponding