1. Background
The optimized, cost-effective design of atmospheric-entry and hypersonic-cruise missions requires an understanding of simultaneously coexisting physical phenomena. These phenomena are related to the extreme thermal conditions of such missions, and include shock waves, excitation of internal energy modes, species interdiffusion, molecular dissociation and recombination and ionization, among others (see Anderson Reference Anderson2006).
Moreover, the instability development that ultimately leads to the transition of boundary layers from laminar to turbulent (see Mack Reference Mack1984; Morkovin Reference Morkovin1988) is strongly conditioned by these thermophysical phenomena. Laminar-to-turbulent transition is considered a potential ‘mission killer’, since the turbulent heat flux that a thermal protection system (TPS) must withstand increases by a factor of 5 with respect to the laminar one (see Wright & Zoby Reference Wright and Zoby1977). It is therefore paramount to establish physically based simplified models for the development of boundary-layer perturbations, accurately capturing their most significant features.
Shock waves can be included in the analysis by simply accounting for the compressibility of a calorically perfect gas (CPG). In order to study the implications of internal-energy-mode excitation, one must employ a nonlinear thermal law, thus assuming a thermally perfect gas (TPG). Dissociation, recombination, ionization and re-neutralization are chemical processes, and therefore require either a chemical non-equilibrium (CNE) or a local thermochemical equilibrium (LTE) assumption. Chemical non-equilibrium considers the finite rate of the reactions, whilst LTE assumes them to occur infinitely fast, reaching instantaneous equilibrium. The reactions that are and are not accounted for depend on the list of species comprising the mixture. Authors performing analyses in LTE, such as Malik & Anderson (Reference Malik and Anderson1991) or Stuckert & Reed (Reference Stuckert and Reed1994), neglected elemental demixing, observed to be negligible unless high temperatures are reached (see Rini, Vanden Abeele & Degrez Reference Rini, Vanden Abeele and Degrez2005; Rini & Vanden Abeele Reference Rini and Vanden Abeele2007). The acronym LTE is thus kept for such a traditional set of hypotheses, hereinafter referring to LTE with elemental demixing as LTEED. Finally, situations where different temperatures are necessary to describe the various groups of energy modes can be reproduced with thermal non-equilibrium (TNE) or thermochemical non-equilibrium (TCNE) flow assumptions.
Stability theories (mostly linear stability theory (LST)) have been used to investigate the early stages of perturbation development since the beginning of the twentieth century (see Mack (Reference Mack1984) for an excellent historical overview). During the last three decades, they have been successfully extended to include high-enthalpy effects. Malik (Reference Malik1989) and Malik & Anderson (Reference Malik and Anderson1991) investigated dissociation effects in a self-similar boundary layer in LTE. Second-mode instabilities were seen to shift to lower frequencies and to increase their growth rates in LTE with respect to CPG. The finite rate of reactions was, however, neglected, until Stuckert & Reed (Reference Stuckert and Reed1994) extended LST analyses to CNE. They concluded that endothermic reactions increase the region of relative supersonic flow and reduce the frequency of second-mode instabilities. Hudson, Chokani & Candler (Reference Hudson, Chokani and Candler1997) considered TCNE effects, by using a distinct temperature for the translational–rotational and the vibrational energy modes. For the cases they considered, the TCNE predictions were found to lie mid-way between the CNE and the CPG ones, suggesting that TNE stabilized the second mode. Mortensen & Zhong (Reference Mortensen and Zhong2016), and later Miró Miró & Pinna (Reference Miró Miró and Pinna2019), extended the preceding framework, featuring a non-catalytic wall boundary condition, to account for the gas–surface interaction occurring during the ablation of a TPS. Mortensen & Zhong (Reference Mortensen and Zhong2016) observed that ablation significantly affected the instability development, making the second mode more unstable, similar to what was predicted by the CPG ablation-mimicking model proposed by Miró Miró & Pinna (Reference Miró Miró and Pinna2018). Lyttle & Reed (Reference Lyttle and Reed2005), Franko, MacCormack & Lele (Reference Franko, MacCormack and Lele2010) and Miró Miró et al. (Reference Miró Miró, Beyak, Pinna and Reed2019) performed sensitivity studies on the models used for several thermophysical gas properties. The predictions of the second-mode development were seen to be strongly modified by modelling inaccuracies leading to a wrong estimation of the boundary-layer size. The transport model was seen to be the one mostly conditioning such predictions: Miró Miró et al. (Reference Miró Miró, Beyak, Pinna and Reed2019) reported that a transport model with a 10 % error with respect to the state of the art could lead to the prediction of a 38 % sooner transition onset location. Klentzman & Tumin (Reference Klentzman and Tumin2013) investigated the receptivity and stability of boundary layers with a dissociating binary oxygen mixture, assuming viscous and inviscid perturbations. They reported practically no effect of chemical-source-term perturbations on the mode shapes, whilst they did see significant differences in the growth rates. They also concluded that the second mode's characteristic time scale is much shorter than the dissociation/recombination time. Bitter & Shepherd (Reference Bitter and Shepherd2015) performed an extensive parametric study on supersonic modes in TNE conditions (yet with frozen chemistry), being later extended to TCNE conditions by Knisely & Zhong (Reference Knisely and Zhong2019a,Reference Knisely and Zhongb,Reference Knisely and Zhongc). Their main finding was that unstable supersonic modes generally appeared as a consequence of highly cooled walls. Mortensen (Reference Mortensen2018) later also observed unstable supersonic modes in configurations with strong nose bluntness. In both studies, supersonic modes were more unstable in the presence of strong dissociation. This suggests that dissociation-induced cooling of the base flow destabilizes supersonic modes, similarly to wall cooling.
Other authors have employed more physically inclusive stability models to investigate high-enthalpy effects. Malik (Reference Malik2003), Chang, Vinh & Malik (Reference Chang, Vinh and Malik1997), Johnson & Candler (Reference Johnson and Candler2005), Zanus, Miró Miró & Pinna (Reference Zanus, Miró Miró and Pinna2020) and Kline, Chang & Li (Reference Kline, Chang and Li2018) worked with linear parabolized stability equations, while Zanus & Pinna (Reference Zanus and Pinna2018) also used nonlinear parabolized stability equations, and Marxen et al. (Reference Marxen, Magin, Shaqfeh and Iaccarino2013), Ma & Zhong (Reference Ma and Zhong2004) and Mortensen & Zhong (Reference Mortensen and Zhong2016) performed direct numerical simulations. All observed that non-parallel effects are also destabilizing in chemically reacting scenarios.
The aforementioned efforts investigated cases with relatively low dissociation levels (mostly of oxygen, but rarely of nitrogen), and completely neglected ionization. Ionization has been, however, considered when investigating convective and absolute instabilities in plasma jet flows. Briggs (Reference Briggs1964) and Bers (Reference Bers1983) analysed the underlying theoretical framework, which was then applied to gas jets by Michalke (Reference Michalke1984). More recent work by Chiatto (Reference Chiatto2014) and Demange et al. (Reference Demange, Kumar, Chiatto and Pinna2018) used LST and absolute instability theory to analyse partially ionized plasma jet flows in the VKI-plasmatron facility assuming LTE conditions.
When a gas is ionized, its viscosity decreases markedly, due to the lower amount of momentum that is exchanged in collisions of ions compared to those of neutral species (see Giovangigli Reference Giovangigli1999). Given the well-known sensitivity of second-mode instabilities to viscosity, one may expect a partially ionized gas to have very different stability characteristics from a neutral gas. These differences could strongly modify the transition-onset location, and must therefore be fully understood before one can properly size the necessary TPS. Moreover, previous high-enthalpy stability analyses did not include a rigorous decoupling of the coexisting physical phenomena. Such a decoupling is necessary in order to understand the effect that each of these phenomena individually have on the instability growth, and to ultimately design solutions effectively controlling the transition dynamics. Some examples of control strategies in low-enthalpy scenarios can be found in the works of Ren, Fu & Hanifi (Reference Ren, Fu and Hanifi2016), Fransson et al. (Reference Fransson, Talamelli, Brandt and Cossu2006) and Paredes et al. (Reference Paredes, Choudhari, Li, Jewell, Kimmel, Marineau and Grossir2018).
For hypersonic flows over sharp wedges or cones at 
$0^{\circ }$ pitch and yaw, in the absence of surface excrescences, second-mode waves are predicted to be the most unstable instability mechanism (see Reed, Saric & Arnal Reference Reed, Saric and Arnal1996). A novel, but yet to be confirmed, interpretation of the second mode, introduced recently by Kuehl (Reference Kuehl2018), proposes the density gradient at the boundary-layer edge as a constraint of the acoustic wave. Traditionally, second-mode waves were commonly referred to as trapped in the region of relative supersonic flow below the relative sonic line (see Mack Reference Mack1984; Fedorov & Tumin Reference Fedorov and Tumin2011). Such an interpretation was based on the mathematical nature of the perturbation equations, but lacked the physical mechanism behind the nature of such equations. This new perspective of the second mode has subsequently been used by Sakakeeny, Batista & Kuehl (Reference Sakakeeny, Batista and Kuehl2019) to suggest why blunt-nose configurations stabilize it, and by Batista & Kuehl (Reference Batista and Kuehl2019) to explore control methodologies to damp and stabilize it.
Within this framework, this article investigates how the development of second-mode waves is affected by shock waves, internal-energy-mode excitation, species interdiffusion, molecular dissociation and ionization. Multiple boundary conditions and flow assumptions are employed, both on the laminar base-flow quantities and on the perturbation terms, in order to perform an effective decoupling of the various physical phenomena. The thermoacoustic interpretation proposed by Kuehl (Reference Kuehl2018) is also considered related to the present results. These results are consistent with this novel interpretation, but it is recognized that the work of Kuehl still requires further validation. Some LST analyses are performed within the frame of the VESTA toolkit (see Pinna Reference Pinna2013) exploiting its versatile automatic derivation and implementation capabilities, presented in Pinna & Groot (Reference Pinna and Groot2014), Pinna et al. (Reference Pinna, Miró Miró, Zanus, Padilla Montero and Demange2019) and Miró Miró (Reference Miró Miró2020). This work builds upon and extends what was presented by the authors in Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a).
2. Flow equations
Various flow assumptions on both the base-flow quantities (
$\bar {q}$) and perturbation quantities (
$q'$) are employed, and they are all summarized in table 1. The simplest base-flow assumption under consideration is CPG, where air is considered a homogeneous mixture with a constant heat capacity. The TPG assumption adds to it by relaxing the constant-heat-capacity constraint, and allowing for a nonlinear thermal law based on assuming that molecules behave like rigid rotors and harmonic oscillators (see § A.4). A TPG also assumes frozen (constant) species concentrations. Assumption CNE-5 considers a five-species mixture (N, O, NO, O
$_{\text {2}}$, N
$_{\text {2}}$), rather than a homogeneous one. Species are allowed to diffuse into one another, and react with dissociation/recombination and exchange reactions (see § A.5). Assumption LTEED-5 assumes that these reactions occur infinitely fast, such that the flow is in equilibrium, and allows for elemental demixing. Assumption LTE-5 introduces the further simplification of neglecting elemental demixing. Assumptions CNE-11, LTEED-11 and LTE-11 are equivalent to CNE-5, LTEED-5 and LTE-5, yet allowing also for ionization/re-neutralization, by assuming air to be a mixture of 11 species (N, O, NO, O
$_{\text {2}}$, N
$_{\text {2}}$, N+, O+, NO+, O
$_{\text {2}}^{+}$, N
$_{\text {2}}^{+}$, e−).
Table 1. Test matrix summarizing the combinations of base-flow and perturbation-flow assumptions under investigation.
Regarding the perturbation flow assumptions, CPG, TPG, CNE-5, LTE-5, CNE-11 and LTE-11 make the same hypotheses as their base-flow counterparts. Additionally, froz-5 and froz-11 assume chemically frozen perturbations (
$\dot {\omega }_s'=0$, where 
$\dot {\omega }_s$ is the species mass production rate), yet allow for 5- and 11-species diffusion fluxes acting on the perturbations (
$J_s'\neq 0$, where 
$J_s$ is the species mass diffusion flux). Finally, diss-11 allows for dissociation/recombination but not ionization/re-neutralization reactions acting on the perturbations (
$\dot {\omega }_{{Ion}\,s}=0$). The diss-11 assumption also accounts for the effect of diffusion fluxes on the perturbations (
$J_s'\neq 0$).
Whenever a base-flow assumption is used in conjunction with its perturbation equivalent (for example CNE-5 base-flow with CNE-5 perturbation) it is referred to with the common flow-assumption naming (CNE-5 in the mentioned example). Contrarily, if different hypotheses are employed on the base-flow and perturbation terms, then the case is designating with the base-flow assumption followed by a hyphen and the perturbation one (e.g. CNE-diss-11, or LTEED-CNE-11).
Depending on the flow assumptions, the system's equations vary slightly. These differences are overviewed in the present section.
2.1. Tensorial notation and the invariant form
The equations in this work are presented in their invariant form, making them valid with independence of the coordinate system, and featuring the metric tensor 
$g_{ij}$. It allows one to keep track not only of the variable value changes, but also of the modification of the space itself due to such basis changes. It is defined as
\begin{equation} g_{ij} = \sum_{k=1}^{3} \frac{\partial \mathscr{X}^k}{\partial x^i} \frac{\partial \mathscr{X}^k}{\partial x^j},\end{equation}
where 
$\mathscr {X}^i$ corresponds to the Cartesian coordinates and 
$x^i$ to the actual coordinate system to be employed. For a Cartesian coordinate system, 
$g_{ij}$ is simply equal to the Dirac delta function 
$\delta _{ij}$. The contravariant metric tensor (
$g^{ij}$) is simply the inverse of the covariant (
$g_{ij}$). Moreover, the expressions presented feature velocities (
$u^i$) in the non-Cartesian reference frame (
$x^i$ rather than 
$\mathscr {X}^i$). Velocities in the Cartesian reference frame (
$\mathscr {U}^i$) can be obtained through
\begin{equation} \mathscr{U}^i = \sqrt{g_{ii}}u^i,\quad \forall\ i\in[1,2,3].\end{equation}
The superscript/subscript notation refers to contravariant/covariant vectorial variables. A vectorial variable is contravariant 
$q^i$ when its components vary with the inverse transformation with respect to the basis change, i.e. they ‘contra-vary’. On the other hand, it is covariant 
$q_i$ if its components vary with the same transformation, i.e. they ‘co-vary’. A spatial covariant derivative is expressed with a comma followed by an index corresponding to the spatial direction with respect to which one is deriving: 
$q^i_{,j}$. The comma derivative notation is restricted to spatial derivatives. For the sake of notational simplicity the spatial subscripts are strictly kept to be 
$i$, 
$j$, 
$k$ and 
$l$ throughout the text. Therefore, variable subscripts containing commas followed by other symbols are not to be regarded as derivatives. The evaluation of covariant derivatives must be done also taking into consideration the curving of the space itself:
\begin{equation} q^j_{,k} = \frac{\partial q^j}{\partial x^k} + \varGamma^j_{ik} q^i,\end{equation}which features the Christoffel symbol of the second kind:
\begin{equation} \varGamma^j_{ik} = \sum_{l=1}^{3} \frac{1}{2 g^{jl}} \left(\frac{\partial g_{li}}{\partial x^k} + \frac{\partial g_{lk}}{\partial x^i} - \frac{\partial g_{ik}}{\partial x^l}\right). \end{equation}A short introduction to tensorial algebra can be found in Pinna & Groot (Reference Pinna and Groot2014) and Miró Miró (Reference Miró Miró2020). For more details, one should refer to the work of Brillouin (Reference Brillouin1964) and Aris (Reference Aris1962).
2.2. Chemical non-equilibrium
A dilute mixture of gases in CNE can be modelled with a separate mass conservation equation for each species (2.5a), three momentum equations (2.5b) and an energy equation (2.5c):
\begin{gather} \frac{\partial \rho_s}{\partial t} + (u^j \rho_s)_{,j} = - J_{s,j}^j + \dot{\omega}_s, \quad \forall\ s \in \mathcal{S}, \end{gather}
\begin{gather}\rho \frac{\partial u^i}{\partial t} + \rho u^j u^i_{,j} = -g^{ij} p_{,j} + \mathbb{T}^{ij}_{,j}, \quad \forall\ i \in [1,2,3], \end{gather}
\begin{gather}\rho \frac{\partial h}{\partial t} + \rho u^j h_{,j} = \frac{\partial p}{\partial t} + u^j p_{,j} + (\kappa^{Fr} g^{ij} T_{,i})_{,j} - \mathcal{J}^j_{,j} + g_{ik} \mathbb{T}^{kj} u^i_{,j}, \end{gather}
where 
$t$ is time, 
$\rho _s$ is the partial density of each species 
$s$, 
$\rho$ is the density of the mixture, 
$\mathcal {S}$ is the set of all species, 
$p$ is the mixture pressure, 
$\mathbb {T}^{ij}$ is the viscous stress tensor, 
$h$ is the mixture enthalpy, 
$\kappa ^{Fr}$ is the frozen thermal conductivity and 
$\mathcal {J}^j$ is the total-energy diffusion flux. The full nomenclature of this article is listed in the supplementary material available at https://doi.org/10.1017/jfm.2020.786.
Alternatively one may substitute the mass conservation equation of one of the species with the mixture continuity equation:
\begin{equation} \frac{\partial \rho}{\partial t} + (u^j \rho)_{,j} = 0.\end{equation}Additionally, those assumptions with ionizable (11-species) mixtures substitute yet another of the species mass conservation equations with the charge balance condition:
\begin{equation} \sum_{s\in\mathcal{S}} \mathcal{Z}_s X_s = 0,\end{equation}
where 
$\mathcal {Z}_s$ is the species unit charge (1 for NO+, 
$-1$ for e−, etc.) and 
$X_s$ is the species molar fraction. Equation (2.7) amounts to imposing the charge of the mixture to be locally neutral.
In order to have a well-conditioned system of equations, it is preferable to enforce (2.6) instead of the mass conservation (2.5a) of the bath species (largest mass fraction), as proposed by Stuckert (Reference Stuckert1991). Similarly, it is preferable to enforce (2.7) instead of the mass conservation of electrons, due to their small molar weight.
The viscous stress tensor is defined as
\begin{equation} \mathbb{T}^{ij} = \lambda g^{ij} u^k_{,k} + \mu \left(g^{jk} u^i_{,k} + g^{ik} u^j_{,k}\right),\end{equation}
where 
$\mu$ and 
$\lambda$ are the first and second dynamic viscosity coefficients. The full list of symbols can be found in the supplementary material. A Maxwellian reactive regime is assumed, justifying the absence of a reactive pressure term in (2.5b) (see Giovangigli Reference Giovangigli1999).
The energy diffusion flux is defined as
\begin{equation} \mathcal{J}^j = \sum_{s\in\mathcal{S}} h_s J^j_s,\end{equation}
where 
$h_s$ is the species enthalpy, and the species mass diffusion fluxes 
$J^j_s$ are defined as
\begin{equation} J_s^j = -\sum_{\ell\in\mathcal{S}} \rho_s \mathcal{D}_{s\ell} d_\ell^j, \end{equation}
where 
$\mathcal {D}_{s\ell }$ are the multicomponent diffusion coefficients of species 
$s$ in species 
$\ell$ and 
$d_\ell ^j$ are the diffusion driving forces. Barodiffusion and thermodiffusion are neglected, since they are known to have a negligible effect (see Scoggins Reference Scoggins2017). The diffusion fluxes are assumed to be ambipolar, meaning that there is no net flow of charge (see Magin & Degrez Reference Magin and Degrez2004). These assumptions, for a gas in thermal equilibrium, allow one to model the diffusion driving force as a molar concentration gradient alone:
\begin{equation} d_\ell^j = g^{ij} X_{\ell,i},\quad \forall\ \ell \in \mathcal{S}.\end{equation}
It is important to note that the action of the ambipolar electric field is implicitly contained in the definition of the multicomponent diffusion coefficients (
$\mathcal {D}_{s\ell }$); see § A.2.
2.3. Local thermodynamic equilibrium with elemental diffusion
Assuming chemical equilibrium, the molar concentration of the various species can be obtained from solving the equilibrium system. Examples of equilibrium systems for the air-5 and air-11 species mixtures can be found in Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) or Miró Miró (Reference Miró Miró2020). Such systems essentially establish a functional dependency of the species concentrations on temperature, pressure and elemental fractions: 
$X_s = X_s(p,T,\mathcal {Y}_E)$. These elemental fractions can be obtained from specific elemental mass conservation equations:
\begin{equation} \frac{\partial \rho_E}{\partial t} + (u^i \rho_E)_{,i} = - J^j_{E,j},\quad \forall\ E\in\mathcal{E},\end{equation}
where 
$\rho _E = \rho \mathcal {Y}_E$ is the partial density of each element and 
$\mathcal {E}$ is the set of all elements. Equation (2.12) has the same morphology as the species mass conservation equation (2.5a), without production terms, since there is no net production or destruction of elements. Elements may change from one species to another, but the amount of atoms of a given element cannot change. The elemental diffusion fluxes (
$J_E^j$) are obtained following the treatment in § 2.3.3 of Miró Miró (Reference Miró Miró2020).
The mixture mass, momentum and energy conservation equations remain unchanged in form with respect to their CNE version – equations (2.6), (2.5b) and (2.5c). The full system of equations to solve in LTEED is thus composed of equations (2.5b), (2.12) and (2.5c), alternatively substituting (2.12) of the bath element for (2.6).
2.4. Local thermochemical equilibrium with constant elemental composition
A very common additional hypothesis that greatly simplifies the governing equations in local thermochemical equilibrium is neglecting elemental demixing: 
$\mathcal {Y}_E = \text {cst}$. This leads to the equilibrium concentrations being exclusively a function of pressure and temperature: 
$X_s = X_s(p,T)$. If the elemental concentrations are assumed constant, their individual conservation equations (2.12) are no longer needed for the problem to have closure. Hence, only the mixture continuity equation (2.6) is required. Moreover, the energy conservation equation (2.5c) can be further simplified by neglecting both frozen and reactive barodiffusion, assuming ambipolar diffusion, and grouping the remaining terms into a reactive thermal conductivity:
\begin{equation} \rho \frac{\partial h}{\partial t} + \rho u^j h_{,j} = \frac{\partial p}{\partial t} + u^j p_{,j} + (\kappa g^{ij} T_{,i})_{,j} + g_{ik} \mathbb{T}^{kj} u^i_{,j},\end{equation}where the thermal conductivity includes both the frozen and reactive addends:
\begin{equation} \kappa = \kappa^{eq} = \kappa^{Fr} + \kappa^{Reac} = \kappa^{Fr} + \sum_{s,\ell\in\mathcal{S}} h_s\rho_s \mathcal{D}_{s\ell} \frac{\partial X_\ell}{\partial T}. \end{equation}Like in LTEED, the mixture mass and momentum conservation equations remain unchanged (2.6) and (2.5b).
2.5. Thermally perfect (frozen) or calorically perfect gas
A TPG, in general, requires the same equation set as a gas in CNE, yet neglecting the species source term (
$\dot {\omega }\approx 0$). However, when there is no surface injection of dissimilar gases, one can simplify the equation set to simply (2.6), (2.5b) and (2.13). The only difference is that the thermal conductivity appearing in the energy equation must be the frozen one 
$\kappa = \kappa ^{Fr}$.
The same applies to a CPG, the difference between a TPG and a CPG being the assumption made on the vibrational and electronic energy modes. These are neglected in CPG, rendering the internal energy and the enthalpy a linear function of temperature (see § A.4).
2.6. Closure of the equation system: thermophysical gas properties
A modelling of the various thermodynamic, transport and chemical properties is needed to provide the system closure. All flow assumptions employ the same set of thermophysical models. The transport properties are obtained from the theory of Chapman and Enskog (see § A.1). The diffusion fluxes are obtained using a self-consistent effective binary ambipolar diffusion model (§ A.2), and the necessary collisional cross-sections are approximated from polynomial-bilogarithmic fits to state-of-the-art data (see § A.3). The chemical source terms are obtained using the law of mass action, with the reaction rate constants proposed by Park, Jaffe & Partridge (Reference Park, Jaffe and Partridge2001) for the 5-species mixtures and by Park (Reference Park1993) for the 11-species ones (see § A.5). All flow assumptions feature a single temperature to describe the thermodynamic state of the gas.
It is important to bear in mind that, whilst the models presented in §§ A.1–A.5 are indistinctly applied to all flow assumptions, the values of the species concentrations needed to do so do vary. Table 2 summarizes the functional dependencies of the species (
$Y_s$) and elemental mass fractions (
$\mathcal {Y}_E$) for the various flow assumptions. The thermodynamic state variables in each case are denoted as ‘TSV’, and the functional relationships 
$f_{eq}(\cdots )$ and 
$f_{stoi}(\cdots )$ refer, respectively, to the chemical equilibrium system (see Miró Miró et al. Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a; Miró Miró Reference Miró Miró2020) and the elemental stoichiometric relations.
Table 2. Functional dependencies of the species (
$Y_s$) and elemental mass fractions (
$\mathcal {Y}_E$) for the various flow assumptions. TSV stands for ‘thermodynamic state variable’.
3. Stability analysis
In order to perform a stability analysis, all flow variables appearing in the system equations (
$q$) are decomposed into a laminar base-flow component (
$\bar {q}$) and a perturbation component (
$q'$). Equations are also simplified after substituting the corresponding ansatz, which for spatial LST leads to
\begin{equation} q(x,y,z,t) = \bar{q}(y) + \tilde{q}(y)\exp\left(\text{i}\left(\alpha x + \beta z - \omega t\right)\right) + \text{c.c.},\end{equation}
where 
$x$, 
$y$ and 
$z$ are the streamwise, wall-normal and spanwise directions, 
$t$ is time, 
$\alpha$ is a complex number of which the real part (
$\alpha _\Re$) corresponds to the streamwise wavenumber and the imaginary part (
$\alpha _\Im$) corresponds to the negative perturbation growth rate, 
$\beta$ is the spanwise wavenumber and 
$\omega$ is the perturbation frequency. Equation (3.1) is the mathematical consequence of assuming a locally parallel base flow, and periodic perturbations in the streamwise and spanwise direction and in time. The perturbation amplitude function (
$\tilde {q}$) is therefore assumed to be an exclusive function of the wall-normal direction (
$y$). Note that there is an intrinsic modelling error arising from the deployment of different assumptions on the base-flow and perturbation components. This is similar to the error associated with the evaluation of the base-flow equations with the LST base-flow ansatz (
$\bar {q}=\bar {q}(y)$), which are subtracted from the ‘base-flow plus perturbation’ equations to obtain the LST equations (see Reed et al. (Reference Reed, Saric and Arnal1996) or § 5.2 in Miró Miró (Reference Miró Miró2020)).
The study of perturbation development within a laminar flow has two clearly distinguishable steps. The first is the resolution of the steady portion of the flow, which is the unperturbed laminar component or base flow, and which provides the 
$\bar {q}$ solution. The second is the subsequent resolution of the stability equations, which spatial LST reduces to a generalized eigenvalue problem on 
$\alpha$, providing the perturbation amplitude functions (
$\tilde {q}$) as eigenvectors (see Arnal (Reference Arnal1993) for an excellent introduction to LST).
3.1. Base-flow problem
The analysis presented in this work solves the laminar base-flow problem by decoupling the inviscid and the viscous flow regions. First the oblique shock-jump relations followed by the one-dimensional Euler equations are solved in the inviscid region. Both the shock-jump relations and the Euler equations vary depending on the flow assumption employed (see § 7 in Miró Miró Reference Miró Miró2020). The inviscid wall values are subsequently imposed as a boundary condition at the boundary-layer edge, followed by the resolution of the viscous boundary-layer region. This effectively imposes a zeroth-order coupling of the inviscid and viscous solutions (see Brazier, Aupoix & Cousteix Reference Brazier, Aupoix and Cousteix1991). The viscous problem is solved after simplifying the flow equations presented in § 2 with the steady boundary-layer assumptions (see White Reference White1991):
\begin{equation} \frac{\partial }{\partial t} = 0,\quad \frac{\partial^{2} }{\partial x^{2}}, \frac{\partial }{\partial z} \approx 0,\quad \frac{\partial }{\partial x} \ll \frac{\partial }{\partial y}.\end{equation}The resulting boundary-layer equations are solved after imposing the appropriate wall boundary conditions. The no-slip and impenetrability conditions demand all velocity components to be zero, while assuming a non-catalytic wall leads to Neumann conditions on the species mass fractions:
\begin{equation} \left.\frac{\partial \bar{Y}_s}{\partial y}\right|_w = 0,\quad \forall\ s\in\mathcal{S},\end{equation}
where the 
$w$ subscript denotes the wall value. The wall is also assumed isothermal, fixing a constant value for 
$\bar {T}_w$. Equation (3.3) is obviously not necessary for the flow assumptions (CPG, TPG or LTE) where the species mass fractions are not system state variables.
The resolution of the base-flow problem is carried out with the DEKAF flow solver (see Groot et al. (Reference Groot, Miró Miró, Beyak, Moyes, Pinna and Reed2018), Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) or § 7 of Miró Miró (Reference Miró Miró2020)). It employs a pseudo-spectral collocation method in the wall-normal direction and a second-order finite-difference method in the marching direction. Moreover, it solves the equations in the 
$\eta$–
$\xi$ variables, after deploying the Illingworth transformation (see White Reference White1991):
\begin{equation} \mathrm{d}\eta = \frac{u_e}{\sqrt{2\xi}} \bar{\rho}\,\mathrm{d}y,\quad \mathrm{d}\xi = \rho_e u_e \mu_e \,\mathrm{d}x, \end{equation}which allows one to scale the domain in order to account for most of the streamwise boundary-layer growth.
3.2. Perturbation problem
The LST equations are retrieved from substituting (3.1) for all flow variables in the flow equations in § 2, and then subtracting the base-flow equation. The boundary conditions imposed on the perturbation-amplitude quantities must retain the mathematical homogeneity of the resulting LST generalized eigenvalue problem. This implies that all perturbation boundary conditions must be an exclusive function of perturbation amplitudes, without forcing terms that would otherwise change the mathematical nature of the problem.
The perturbation boundary conditions differ slightly from those on the base-flow quantities. At the wall, the no-slip and impenetrability conditions are also applicable and demand all velocity perturbations at the wall to be identically zero. Similarly, the isothermal wall assumption also makes 
$\tilde {T}_w = 0$. The pressure perturbation, however, is not zero. The LST wall-normal momentum equation ((2.5b) with 
$i=2$ after substituting (3.1) and operating) is commonly imposed at the wall as a compatibility condition (see Mack Reference Mack1984). For flow assumptions with the species concentrations as state quantities (CNE), one can either apply the LST ansatz on (3.3) or employ separate species wall-normal momentum equations for each individual species:
\begin{equation} \rho_s \frac{\partial u^i}{\partial t} + \rho_s u^j u^i_{,j} = -g^{ij} p_{s,j} + \mathbb{T}^{ij}_{,j},\quad \begin{array}{l} \forall\ s\in\mathcal{S} \\ i=2 \end{array}, \end{equation}
where 
$p_s$ is the species partial pressure. Miró Miró & Pinna (Reference Miró Miró and Pinna2017) observed that the use of the species momentum compatibility condition instead of the non-catalytic improved the condition number of the matrix system by four orders of magnitude, not modifying the stability characteristics.
Note that the contributions of species interdiffusion and momentum exchange to the momentum balance have been neglected in (3.5). This may seem like an overly restrictive simplification, but this restriction is already implicit to the usage of mixture momentum equations (2.5b) to describe the motion of all species. There is therefore no loss in generality in enforcing equation (3.5).
Regarding the free-stream perturbation boundary conditions, two different methodologies are investigated in this work. They are compared in order to identify the isolated effect of the shock wave on the instabilities.
The first one consists of ignoring the shock position, and extending the wall-normal domain far beyond it. All perturbation amplitudes must damp in the far field, and therefore one can impose Dirichlet conditions. However, in order to account for the finite size of the domain, one of the perturbation amplitudes is oftentimes liberated from the homogeneous restriction with an additional compatibility condition in the free stream. In this work, the wall-normal momentum equation is employed to account for 
$\tilde {p}$, which is liberated in the CPG, LTE and TPG solvers. Similarly, in the LTEED and CNE solvers, it is the wall-normal velocity's perturbation amplitude (
$\tilde {v}$) that is liberated in the far field, using the continuity equation (2.6) as a compatibility condition. This standard approach is hereinafter referred to as the ‘free stream Dirichlet’, and was also taken by Malik (Reference Malik1989), Hudson et al. (Reference Hudson, Chokani and Candler1997), Mortensen & Zhong (Reference Mortensen and Zhong2016) and many others.
The alternative approach is to truncate the wall-normal domain at the shock location, and to impose the Rankine–Hugoniot shock-jump relations, after substituting the LST ansatz (3.1) on both the post-shock variables and the shock height. The pre-shock region is considered to be unperturbed. Such a treatment of the far-field boundary was previously explored by Esfahanian (Reference Esfahanian1991), Stuckert & Reed (Reference Stuckert and Reed1994), Knisely & Zhong (Reference Knisely and Zhong2019b) and Pinna & Rambaud (Reference Pinna and Rambaud2013), among others. The full mathematical development with the nomenclature of this article, detailed in the supplementary material, can be found in § 6 of Miró Miró (Reference Miró Miró2020).
4. Results
The effects of ionization and dissociation on flow stability are investigated under free-stream conditions similar to three flight-envelope points in a Martian return that are presented in table 3 (see Howe Reference Howe1989). The free-stream Mach number 
${M}_\infty$ is computed based on the frozen speed of sound:
\begin{equation} a^{Fr} = \sqrt{\frac{c_p}{c_v} R T},\end{equation}
where 
$c_p$ is the frozen mixture heat capacity at constant pressure, 
$c_v$ is the frozen mixture heat capacity at constant volume and 
$R$ is the mixture-specific gas constant. The test gas is air in all cases, with a free-stream composition of 23.3 % O2 and 76.7 % N2 in mass. The considered geometry is a sharp 
$10^{\circ }$ wedge of length 
$L = 5$ m. The mean free path of oxygen and nitrogen molecules at 45 km altitude is of the order of 50 
$\mathrm {\mu }$m, thus justifying the continuum assumption. The high wall temperatures (see table 3) are clearly beyond what can be sustained by any existing TPS material. Assuming an isothermal wall results in surface heat fluxes that differ substantially from what one may encounter in realistic flow scenarios with radiative equilibrium at the surface. Moreover, hypersonic vehicles typically feature a blunt nose, where the flow is close to equilibrium, followed by a region of recombination and re-neutralization (see Malik Reference Malik2003; Miró Miró & Pinna Reference Miró Miró and Pinna2019). However, the high number of coupled phenomena in such situations largely complicates the analysis and the ultimate identification of the individual contribution of each of them. The present analysis is formulated with the goal of understanding the underlying physics in dissociating and ionizing boundary layers, rather than accurately reproducing an actual engineering problem. The intent is to identify and gain insight into the effect that individual phenomena associated with high-enthalpy gases have on the growth or stabilization of second-mode waves. Future analyses should build on these results and add more effects.
Table 3. Test matrix.
4.1. Laminar base flow
The laminar base-flow field is obtained using the DEKAF solver (see Groot et al. (Reference Groot, Miró Miró, Beyak, Moyes, Pinna and Reed2018), Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) or § 7 of Miró Miró (Reference Miró Miró2020)). As laid out in § 3.1, the inviscid flow region is characterized by solving the corresponding shock-jump relation for each flow assumption, and then employing a one-dimensional Euler solver for the CNE cases. This approach notedly neglects the shock curving induced by the dissociation of molecules in the inviscid-flow region (see § 15.3 in Anderson (Reference Anderson2006)).
The boundary-layer edge (or inviscid-wall) profiles for the first flight-envelope point (see table 3) are displayed in figure 1. Boundary-layer edge quantities are denoted with the subscript 
$e$. One can clearly see the effects of the chemical activity in the inviscid flow region, which progressively dissociates O
$_{\text {2}}$ advancing downstream (figure 1e), and consequently reduces pressure (figure 1d) and temperature (figure 1c) and accelerates the flow (figure 1b). The result is an increase in the edge Mach number (figure 1a). An interesting feature of figure 1(d) is the fact that the non-equilibrium edge pressure undershoots that of the equilibrium case. This artifact is a consequence of neglecting the curvature of the shock due to the chemical activity. In reality, this curving brings the shock closer to the surface, and increases the pressure, therefore avoiding the undershooting. However, the objective of the study is to investigate the implications of ionization as a phenomenon on instability development, rather than truthfully reproducing an engineering problem. The mentioned modelling inaccuracies are therefore accepted, since they do not conflict with such an objective. Also note that blunt geometries, like those one would employ in Martian return missions, display the opposite trend in the inviscid flow region – one has equilibrium near the tip, with molecular recombination and re-neutralization occurring downstream of it.
Figure 1. Evolution of the boundary-layer edge (or inviscid-wall) profiles for the first flight-envelope point (
${M}_\infty =45$) in table 3 and for various flow assumptions.
The other two flight-envelope points (
${M}_\infty =37.5$ and 
${M}_\infty =30$ in table 3) present practically no finite-rate chemical activity in the inviscid-flow region. The post-shock (or boundary-layer edge) values of the quantities of interest are therefore constant in the CNE base flows. The edge values for these two flight-envelope points and all flow assumptions are reported in table 4.
Table 4. Post-shock (or boundary-layer edge) values of the different base-flow quantities of interest for the second and third flight-envelope points in table 3.
It is also important to note that the test cases feature shocks that are very close to the surface. This is due to the high free-stream Mach numbers under consideration. The corresponding shock angles, relative to the surface of the 
$10^{\circ }$ wedge, are laid out in table 5.
Table 5. Relative shock angle (in degrees) with respect to the 
$10^{\circ }$ wedge surface for the flight-envelope points in table 3 and for the various flow assumptions.
The boundary-layer profiles at 
$x=1$ m, corresponding to the aforementioned edge conditions, and with the thermal wall boundary conditions in table 3 are presented in figure 2. The different levels of flow ionization are clearly visible from the electron mole-concentration profiles (figure 2, panels iv). For the first flight-envelope point, both the LTE-11 and CNE-11 assumptions predict a maximum ionization of about 10 %, which goes down to about 1 % for the second and between 0.01 % and 0.1 % for the third. In the first two flight-envelope points the most prevalent ion is N+, whilst in the third it is NO+. Aside from being ionized, the test gas is also clearly dissociated. Molecular nitrogen is completely dissociated into N close to the wall for the first and second flight-envelope points (figure 2a,b.iii), whilst approximately half of it (
$\bar {X}_\text {N} \approx 0.3$) is dissociated close to the wall for the third point (figure 2c.iii).
Figure 2. Boundary-layer profiles at 
$x=1$ m for the flight-envelope points in table 3 and for various flow assumptions.
As a consequence of the commented ionization levels for the 
${M}_\infty =45$ case, the temperature profiles for the chemically reacting flow assumptions that account for ionization (CNE-11 and LTE-11) differ significantly from those that do not (CNE5 and LTE-5). Ionization induces a cooling of the boundary layer, which is clearly visible in figure 2(b.ii). The levels of ionization in the two other flight-envelope points are not sufficient to cause major differences between the 5- and the 11-species temperature predictions (with and without ionization) – the CNE-5 and LTE-5 temperature profiles in figure 2(b,c.ii) are, respectively, coinciding with the CNE-11 and LTE-11 profiles.
Neglecting all chemical activity (as done by the CPG and TPG assumptions) leads to a clear overestimation of the boundary-layer temperature (see figure 2a–c.ii). This overestimation is more so when the excitation of internal energy modes is neglected (CPG). The cooling related to internal energy excitation, dissociation and ionization also leads to a decrease in the boundary-layer size.
4.2. Linear stability theory analyses
Stability studies are carried out using VESTA's LST solver (see Pinna Reference Pinna2013), exploiting the extended ADIT capabilities (see Pinna et al. (Reference Pinna, Miró Miró, Zanus, Padilla Montero and Demange2019) or § 8 of Miró Miró (Reference Miró Miró2020)). The wall-normal domain is discretized using the FDq-8 method with 400 points distributed with Malik's mapping (see Malik Reference Malik1990). The mapping parameter (
$y_i$) was fixed slightly above the boundary-layer height – at the wall-normal position where the transformed Illingworth variable 
$\eta$ reaches a value of 6 (see (3.4a,b)). The FDq differentiation matrices and point distributions were obtained with the library kindly provided by Dr M. Hermanns (see Hermanns & Hernández Reference Hermanns and Hernández2008; Paredes et al. Reference Paredes, Hermanns, Le Clainche and Theofilis2013).
4.2.1. Consistent base-flow and perturbation hypotheses
The first comparison features consistent hypotheses for the base-flow and the perturbation quantities, thus reducing the list of flow assumptions presented at the beginning of § 2 to CPG, TPG, CNE-5, LTE-5, CNE-11 and LTE-11. The second-mode growth rates are displayed in figure 3 for the three flight-envelope points in table 3 as a function of the perturbation frequency (
$F = \omega /2{\rm \pi}$).
All three growth-rate plots display the destabilizing nature of internal-energy-mode excitation – considered by the TPG but not by the CPG assumption. Comparing the TPG predictions with those considering finite-rate chemistry (CNE-5 and CNE-11), the trends vary slightly depending on the flight-envelope point. For the first one (figure 3a) the differences between the TPG and CNE-5 predictions suggest that dissociation is significantly destabilizing the boundary layer and increasing the range of unstable frequencies. This trend is reversed when also accounting for ionization, resulting in the differences that one can observe between the CNE-5 and CNE-11 growth-rate curves. The other two flight-envelope points (figure 3b,c), which feature lower levels of ionization (see figure 2, panels iv), do not display significant differences between the CNE-5 and CNE-11 growth rates. They are, however, distinct from the TPG predictions – displaying a slight stabilization and a shift of the most unstable mode to higher frequencies. Figure 3 also suggests that when the boundary layer reaches equilibrium conditions (LTE-5 or LTET11 assumptions), it is significantly more unstable than under frozen (TPG) or finite-rate chemistry (CNE-5 or CNE-11). This trend was also observed in the investigation of flow assumption adequacy carried out by Zanus, Miró Miró & Pinna (Reference Zanus, Miró Miró and Pinna2019). Once again, the distinct levels of ionization featured at the three flight-envelope points (see figure 2, panels iv) explain why the LTE-5 and LTE-11 predictions coincide for the third flight-envelope point (figure 3c), differ slightly for the second (figure 3b) and present major differences for the first (figure 3a).
A comparison of the growth rates for the three different flight-envelope points (figure 3) obtained with the same LTE assumption displays a non-monotonic variation. Looking at the results with Mach 45 and 37.5, it appears as if the LTE curves are converging towards the TPG. However, the Mach 30 case moves away from this trend. This is presumably linked to the variation of the wall temperature that accompanies the variation of the Mach number (see table 3).
In order to further investigate the traits observed in figure 3, figure 4 presents the temperature perturbation amplitude for the most unstable frequency (i), together with the corresponding laminar base-flow density profile (ii). The TPG density gradient is systematically stronger than the CPG one. The confinement of second-mode waves associated with a strengthened thermoacoustic impedance, within the framework proposed by Kuehl (Reference Kuehl2018), is thus a possible candidate driving the destabilization due to internal-energy-mode excitation. This claim is also supported by the stronger temperature perturbation amplitude in figure 4 (panels i) that the TPG cases display with respect to the CPG.
Figure 4. Second-mode temperature perturbation mode shapes for the most unstable frequency (i) and base-flow density (ii) at 
$x=1$ m for the flight-envelope points in table 3 and for various flow assumptions introduced at the beginning of § 2. All assumptions share the base-flow and perturbation hypotheses.
This clearly distinctive second peak in the temperature perturbation, which is visible for the CPG and TPG cases around 
$y\approx 0.01$ m, is also found, yet with a smaller amplitude, for the CNE cases in the third flight-envelope point (figure 4c.i). However, this peak in the CNE cases is further weakened for the second flight-envelope point (figure 4b.i), and reaches unrecognizable levels in the first (figure 4a.i). If the mentioned second temperature perturbation peak is a consequence of the instability confinement infringed by the thermoacoustic impedance that is a density gradient, one would expect it to be larger for the dissociating and ionizing assumptions, which display a stronger density gradient in figure 3 (panels ii). The fact that this is not the case suggests that either the thermoacoustic framework needs further improvement or there is another physical mechanism affecting the development of second-mode instabilities in dissociating and ionizing boundary layers, beyond confinement caused by the thermoacoustic impedance. This is further investigated in the following subsection (§ 4.2.2). The damping of the second temperature perturbation peak is even stronger when assuming LTE. This supports the claim that it is related to the stronger dissociation and ionization assumed in the LTE and CNE cases with respect to the TPG.
Instead of examining the stability characteristics of a single streamwise station, one can also integrate the growth rates of the various perturbation frequencies in the 
$x$ direction, thus obtaining the 
$N$-factor curves. The envelopes of these curves for frequencies between 50 and 800 kHz are displayed in figure 5 for the three points of the flight envelope. Individual frequency curves corresponding to the CPG assumption are also displayed for purely illustrative purposes. Figure 5 displays the logical consequences of the trends observed and commented for figure 3. An interesting additional behaviour that one can identify in figure 5(b) is the evolution of the CNE envelopes that coincide with the TPG one close to the leading edge, and eventually depart from it (at around 
$x=2$ m). This is a consequence of the increasing modification of the flow due to finite-rate dissociation and ionization while advancing downstream.
In order to offer a better visualization of the mentioned trends, one can quantify the contribution of the various physical phenomena to the ultimate 
$N$ factor at a given streamwise location. This is done at 
$x=4$ m, by computing the increase in 
$N$ factor when using successively more physically inclusive assumptions. The results are displayed in figure 6, where a distinction is made depending on whether chemical activity is assumed to occur at a finite rate (figure 6a), or to correspond to equilibrium conditions (figure 6b). Figure 6 displays very clearly a rather counterintuitive trend observed in figures 3 and 5: ionization is seen to stabilize the flow when considering finite-rate reactions, yet destabilize it under equilibrium conditions. The following subsections analyse various flow assumptions, featuring different hypotheses for the base-flow and perturbation quantities, in order to investigate the underlying reason for this trend reversal.
Figure 6. Variation in the maximum 
$N$ factor at 
$x=4$ m due to the different physical phenomena affecting simultaneously the base-flow and perturbation quantities. Case details are summarized in table 3. (a) Finite rate and (b) equilibrium.
4.2.2. Base-flow cooling, diffusion and chemical source term
The temperature perturbation amplitude shapes in figure 4 (panels i) suggest that there are other major actors affecting the development of instabilities. In order to shed light on what such phenomena may be and how they are interrelated, various combinations of base-flow and perturbation hypotheses are compared. Figure 7 compares the TPG against the CNE (5- or 11-species) results, as well as those corresponding to a series of intermediate flow assumptions (see the beginning of § 2).
Aside from the TPG case, all others feature chemically reacting basic-state solutions. The CNE-TPG cases assume TPG conditions for the perturbation variables 
$q'$, thus effectively neglecting the diffusion-flux and source-term perturbations (
$J_s'$ and 
$\dot {\omega }_s'$). The CNE-froz cases simply neglect the source-term perturbations, whilst the CNE-diss-11 case neglects the contribution of ionization reactions to the source-term perturbation.
One can clearly see that, for equivalent mixtures (5- or 11-species) the CNE-TPG assumption always predicts the largest range of unstable frequencies and the maximum growth-rate value. An exception to this trend is seen for the 5-species cases in the first flight-envelope point (figure 7a.i). The CNE-froz-5 case actually has a maximum growth that is larger than that of the CNE-TPG-5 case. However, the range of unstable frequencies is substantially smaller. It is therefore reasonable to expect the CNE-TPG-5 case to predict a sooner transition under the present assumptions. For the remainder of the article, the term ‘stabilizing’ or ‘destabilizing’ is used to refer to how something affects the entire frequency range, and not only the maximum growth rate. The CNE-TPG-5 case is unstable over a larger range of frequencies. Assuming that the enlargement of the unstable frequency range is consistent at all streamwise locations (as is the case), one can affirm that the CNE-TPG-5 case is more unstable than the CNE-froz-5 one under the present assumptions. Figure 7(b,c) also presents practically coinciding stability predictions with the 5- and 11-species assumptions for the second and last flight-envelope points. This indeed reinforces the idea that the finite-rate ionization occurring in them is sufficiently low so as to neglect it.
Aside from the mentioned exception, gradually accounting for diffusion flux, dissociation and ionization contributes to a reduction of the predicted range of unstable frequencies, together with the maximum growth rate. This is a clear example of the two major competing effects in dissociating and ionizing boundary layers:
(i) The destabilization due to the cooling of the base-flow field.
(ii) The stabilization due to the dissociation chemical source terms acting on the perturbations.
The present analysis suggests that ionization acting on the perturbation terms is also stabilizing, and that there is yet another actor that is actually the most important one in many cases:
(iii) The stabilization due to the diffusion fluxes acting on the perturbations.
The CNE-TPG assumptions only experience the first effect, and therefore make the most unstable of all predictions. The CNE-froz cases also include the third effect, whilst the CNE ones have all three.
Despite the framework requiring additional validation, it is of interest to explore these trends using the thermoacoustic interpretation proposed by Kuehl (Reference Kuehl2018), thus effectively placing scrutiny on the framework. To that end, figure 8 displays the temperature perturbation mode shape corresponding to the most unstable frequency (figure 8a) together with the base-flow density profile (figure 8b). The comparison is focused on the first flight-envelope point (
${M}_\infty =45$), and on 5-species mixtures but analogous results are obtained for the other flow conditions and for the 11-species mixtures. The CNE-TPG-5 assumption makes the same hypotheses on the perturbation terms as the TPG one. However, it makes significantly different base-flow hypotheses (CNE), which actually lead to predicting stronger density gradients due to dissociation (figure 8b). The CNE-TPG-5 perturbations therefore have the same characteristics as the TPG ones, yet the base flow field features a significantly stronger thermoacoustic impedance. The thermoacoustic interpretation proposed by Kuehl (Reference Kuehl2018) predicts that this must result in a destabilization of second-mode waves, which is consistent with what figure 7 displays. The temperature perturbation shapes in figure 8(a) support this interpretation – the second peak, around 
$y\approx 0.01$ m, is indeed much stronger for the CNE-TPG-5 case than it is for the TPG one. This is arguably a consequence of the commented stronger instability confinement. The effect of the diffusion fluxes is to damp this instability peak, as can be seen by a comparison of the CNE-froz and the CNE-TPG temperature perturbations in figure 8(a). This damping is likely the manifestation of the stabilizing effect of the diffusion fluxes.
It is important to note that, unlike that observed in Miró Miró et al. (Reference Miró Miró, Pinna, Beyak, Barbante and Reed2018b), the diffusion fluxes are seen to have a major impact on the stability predictions. This suggests that indeed for strongly reacting boundary layers such as those under consideration in this article, they can no longer be faultlessly disregarded.
Another interesting conclusion to draw from this comparison is that the reason why the ionizing CNE case (CNE-11) predicts a more stable boundary layer than the non-ionizing one (CNE-5) at the first flight-envelope point (
${M}_\infty = 45$) is not the modification of the mixture viscosity caused by the larger concentration of ions. If this were the case, then the CNE-TPG equivalents, which feature the same base-flow field as the CNE cases, would also experience such an ionization-related stabilization. Instead, the CNE-TPG-11 case is more unstable and over a wider range of frequencies than the CNE-TPG-5 one. The cooling of the laminar boundary layer remains the driver of the mentioned destabilization associated with ionization and distinguishing cases CNE-TPG-11 and CNE-TPG-5.
The 
$N$-factor envelopes resulting from the integration of the instabilities corresponding to the assumptions compared in figure 7 are presented in figure 9. The general trends lie within that already mentioned for figure 7, and perfectly display the three mentioned competing effects and their distinct importance for the different points within the flight envelope. By taking the different values of the 
$N$-factor envelopes at 
$x=4$ m, these effects can be very clearly quantified. The visualization of the different contributions is displayed in figure 10, where a distinction is made as to whether the mixture is assumed to be 5- or 11-species air. An interesting feature that is clearly visible from figure 10 is that for the first flight-envelope point (
${M}_\infty = 45$) diffusion is seen to have a much stronger stabilizing effect when the mixture contains ions (11 species) than when it does not (5 species). This is a plausible explanation for the overall stabilizing effect of ionization observed in figure 6(a) when assuming finite-rate chemistry.
Figure 10. Variation in the maximum 
$N$ factor at 
$x=4$ m due to the different physical phenomena. Case details are summarized in table 3.
It is also interesting to observe that, for the third flight-envelope point (
${M}_\infty =30$), the three identified competing effects happen to cancel each other out completely. The result of this is that the TPG and CNE curves in figure 9(c) coincide. This does not imply, however, that there is no chemical activity in the flow, or that it can be neglected altogether. A variation of the flow configuration could easily lead to an imbalance of these effects that circumstantially neutralize each other.
4.2.3. Equilibrium perturbations and elemental diffusion
Assuming LTE perturbations implies that they are expected to reach equilibrium instantaneously. This assumption is questionable, since, in reality, the high frequency with which second-mode waves oscillate in time may not allow for them to reach local equilibrium. Moreover, the LTE assumption (see § 2.4) neglects elemental demixing, and therefore imposes constant elemental concentration throughout the domain. This constitutes an oversimplification of the problem, since, in reality, the various species interdiffuse in a way that does not necessarily maintain a constant elemental concentration, driven by the different temperature gradients. This subsection challenges these two simplifications, by comparing LTE with LTE-CNE and LTEED-CNE assumptions. The LTE assumptions relax the hypothesis of perturbations being in equilibrium, whereas the LTE-CNE and LTEED-CNE assumptions also allow for elemental demixing in the laminar base-flow solution.
Figure 11 displays the laminar base-flow elemental mass fraction of nitrogen (i) and temperature (ii) at 
$x=1$ m for the flight-envelope points in table 3. The difference between the LTE and LTEED elemental fractions shows that elemental demixing is by no means negligible at any of the flight-envelope points. However, the fact that it has practically no influence on the temperature profiles suggests that second-mode waves may not be significantly affected by it.
The second-mode growth rates are presented in figure 12 for the LTE (figure 12a), LTE-CNE (figure 12b) and LTEED-CNE (figure 12c) flow assumptions, and for both 5- and 11-species mixtures. One can clearly see that relaxing the equilibrium perturbation assumption (figure 12b with respect to 12a) strongly modifies the predicted effect of ionization (accounted for with 11-species mixtures, but not with 5-species ones). When assuming LTE perturbations (figure 12a), ionization is predicted to destabilize the flow and to change the instability range to higher frequencies. When relaxing the LTE perturbation hypothesis (figure 12b), the instability range is still predicted to change to higher frequencies due to ionization, but the growth rates decrease, instead of increasing. This suggests that the apparent ionization-driven destabilization observable in figure 12(a) is nothing other than an artifact of assuming perturbations in equilibrium.
Allowing for elemental demixing in the laminar base flow (figure 12c with respect to 12b) does not introduce a substantial modification of the predicted instability characteristics. This suggests that second-mode waves are not strongly affected by demixing.
Only the results for the first flight-envelope point (
${M}_\infty =45$) are presented in figure 12. However, the comparisons for the other points display the same traits, yet less accentuated.
The 
$N$-factor envelopes presented in figure 13 confirm that the trends noted in figure 12 are preserved along the entire streamwise range. One can also proceed in the same way as in § 4.2.1, and quantify the differential effect of dissociation and ionization on the 
$N$ factors at 
$x=4$ m. The result is shown in figure 14, which displays the same comparison as figure 6(b), yet without assuming equilibrium perturbations. The difference between figures 14(a) and 14(b) depends on whether the base-flow solution neglected (figure 14a) or not (figure 14b) elemental demixing. The trends presented in figure 14 are indeed closer to those in figure 6(a) than they are to those in figure 6(b). This suggests that the observation made when assuming LTE (figure 6b), that ionization further destabilizes the flow in equilibrium conditions, is actually an artifact of assuming perturbations in equilibrium.
Figure 14. Variation in the maximum 
$N$ factor at 
$x=4$ m due to the different physical phenomena. Chemical activity assumed in equilibrium, with the LTE-CNE assumptions (
$\bar {J}_E=0$ yet 
$J_E'\neq 0$). Case details are summarized in table 3.
4.2.4. Supersonic modes
Second-mode instabilities have slightly distinct characteristics depending on their wave speed relative to the boundary-layer edge velocity. If the difference between the boundary-layer edge velocity and the wave speed is supersonic, one refers to them as supersonic modes, which have a series of particular characteristics (see Chuvakhov & Fedorov Reference Chuvakhov and Fedorov2016; Knisely & Zhong Reference Knisely and Zhong2019b; Bitter & Shepherd Reference Bitter and Shepherd2015). Mortensen (Reference Mortensen2018) also reported additional unstable supersonic modes that were unrelated to second-mode waves. Supersonic modes receive energy from their ‘synchronization’ with the continuous branch of the instability spectrum, similar to how the second mode does with either the fast or slow mode (following the naming convention proposed by Fedorov & Tumin (Reference Fedorov and Tumin2011)). This causes a non-exponential decay of the perturbation amplitude in the free-stream region, as well as an increase in the upper limit of the range of unstable frequencies. The latter can have a strong effect on the ultimate transition-onset predictions, since it effectively shifts the 
$N$-factor envelopes vertically (see Bitter & Shepherd Reference Bitter and Shepherd2015).
In order to better understand the various trends observed in the preceding sections, the instability growth rates of the first flight-envelope point (table 3) at 
$x=1$ m are plotted as a function of the wave speed, rather than the wave frequency, in figure 15. In order to effectively compare cases with distinct boundary-layer edge velocities and Mach numbers (see figure 1 and table 4), the horizontal axis in figure 15 is the non-dimensional wave speed relative to the sonic line based on the frozen speed of sound (4.1):
\begin{equation} \frac{c}{u_e} - \left(1-\frac{1}{{M}_e}\right) = \frac{\omega/\alpha_\Re}{u_e} - \left(1-\frac{\sqrt{\dfrac{c_{pe}}{c_{ve}}R_eT_e}}{u_e}\right). \end{equation}Therefore, instabilities will be relatively supersonic on the left of 0 in figure 15 and subsonic on the right. Figure 15 compares the same flow assumptions as in § 4.2.2, with their 5-species variant in figure 15(a) and their 11-species one in figure 15(b).
Neglecting all chemical activity in both the base-flow and perturbation quantities (TPG assumption) leads to instabilities that remain always subsonic relative to the boundary-layer edge velocity. The inclusion of dissociation and ionization on the base-flow field (CNE-TPG assumptions) increases the instability wave speed. However, considering the effect of diffusion fluxes on the perturbation terms (CNE-froz assumptions) decreases this wave speed, and makes instabilities relatively supersonic.
The appearance of unstable supersonic modes is usually associated with the boundary-layer cooling induced by either cold walls (see Mack Reference Mack1984; Bitter & Shepherd Reference Bitter and Shepherd2015; Knisely & Zhong Reference Knisely and Zhong2019b,Reference Knisely and Zhongc) or strong dissociation (see Chang et al. Reference Chang, Vinh and Malik1997; Miró Miró et al. Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a; Chang, Kline & Li Reference Chang, Kline and Li2019). However, the present results suggest that, in dissociating and ionizing scenarios with a hot wall, it is not the base-flow cooling that causes the modes to become supersonic, but rather the diffusion fluxes acting on the perturbations. It is important to stress that this claim does not conflict with the common understanding of supersonic modes, which are promoted by wall cooling. Repeating the test cases displayed in figure 15 with colder walls is expected to extend further the range of supersonic wave speeds.
The inclusion of the dissociation and ionization source-term perturbations (
$\dot {\omega }_s'$) has once again a different impact on the predictions made with 5- or 11-species mixtures. Comparing the CNE to the CNE-froz curves in figure 15, one observes that non-zero 
$\dot {\omega }_s'$ lead to a reduction of the range of unstable wave speeds. However, for 11-species assumptions (figure 15b) this actually makes the instabilities remain relatively subsonic over the entire unstable frequency range.
Similar trends are obtained when performing the same analysis for the other two flight-envelope points.
4.2.5. Transport model in CPG
A rather interesting result is obtained from a comparison of the CPG instability predictions with Sutherland's viscosity law, or with the Chapman and Enskog (CE) transport model. As discussed by Zanus et al. (Reference Zanus, Miró Miró and Pinna2019), the inaccuracy in the transport modelling associated with the use of Sutherland's law oftentimes casually counteracts the inaccuracy in the thermodynamic modelling linked to neglecting internal-energy-mode excitation and dissociation. This is once again observed for the three flight-envelope points under comparison (see table 3). Figure 16 displays the integrated 
$N$-factor envelopes as well as the individual frequency 
$N$-factor curves for the CPG case with the CE model. One can indeed see how the CPG predictions with Sutherland's law happen to lie very close to the CNE-11 ones. As pointed out by Zanus et al. (Reference Zanus, Miró Miró and Pinna2019), however, one should refrain from identifying the CPG assumption with Sutherland's viscosity law as having any modelling accuracy. It actually poses an excellent example of how one can circumstantially reach the right answer for the wrong reason.
4.2.6. Free-stream boundary condition
The test cases under consideration feature a shock that is very close to the surface (see § 4.1 and table 5). It is therefore reasonable to believe that instabilities do not have a sufficiently large wall-normal stretch so as to fully decay before reaching the shock. Also, as seen in § 4.2.4, the most unstable perturbation is a supersonic mode over a large fraction of the streamwise range. Supersonic modes are known to have a non-exponential decay in the free-stream region (see Knisely & Zhong Reference Knisely and Zhong2019a,Reference Knisely and Zhongb,Reference Knisely and Zhongc), further increasing the amplitude with which they reach the shock.
It is therefore questionable whether one should apply Dirichlet boundary conditions, as was done in the previous sections, or if it is more appropriate to use the linearized Rankine–Hugoniot relations instead. Figure 17 displays the variation in the maximum 
$N$ factor due to the imposition of the shock boundary condition, rather than the Dirichlet ones. The low values of the shock 
$N$-factor difference across all flight-envelope points support the choice of the free-stream boundary conditions made throughout this article.
Figure 17. Variation in the second-mode maximum 
$N$ factor due to the use of the shock boundary condition rather than the free-stream Dirichlet condition (see § 3.2). Test cases correspond to the flight-envelope points in table 3 and feature various flow assumptions introduced at the beginning of § 2.
One must note that the wall-normal height of the shock remains smaller than five times the boundary-layer thickness for all flow assumptions, flight-envelope points and streamwise positions. The work of Chang (Reference Chang2004) suggests that such attached shocks are expected to modify the stability characteristics. Whatever this modification may be, it is not seen to translate into meaningful differences in the 
$N$-factor envelopes for the cases considered in this work.
5. Conclusions
An extensive investigation was carried out using LST on boundary-layer, laminar, base-flow fields to determine how ionization and dissociation affect the development of second-mode instabilities within a boundary layer. Three test cases were investigated on a sharp 
$10^{\circ }$ wedge with an isothermal wall and free-stream conditions similar to three flight-envelope points in atmospheric reentry missions (Martian return). Several different flow assumptions, featuring various combinations of base-flow and perturbation hypotheses, were employed in the effective decoupling of the various physical phenomena simultaneously interacting in such scenarios.
The boundary layers of the three flight-envelope points were seen to have maximum ionization levels ranging from 0.01 % to 10 %. Successively accounting for internal-energy-mode excitation, dissociation and ionization was seen to reduce the boundary-layer size and cool it. The effect was enhanced when the chemical activity was assumed in equilibrium, rather than accounting for the finite rate of reactions. Dissociation was also seen to introduce an additional reduction in the wall density due to the decrease in the mixture molar mass. All these features contributed to increasing the thermoacoustic impedance (density gradient) near the boundary-layer edge along with a strengthening of the temperature perturbation amplitude and a destabilization of second-mode instabilities. The results analysed in this work therefore seem to be consistent with the interpretation proposed by Kuehl (Reference Kuehl2018), according to which the thermoacoustic impedance confines second-mode waves to the boundary layer. However, this interpretation still requires further validation.
The excitation of internal energy modes as well as the dissociation of molecular species were seen to destabilize the boundary layer and ultimately lead to higher 
$N$ factors. Conversely, ionization was seen to be mainly stabilizing. That is, when perturbations were assumed in non-equilibrium conditions. Restricting their behaviour by enforcing on them the equilibrium condition predicted ionization to be destabilizing. The modification of the laminar base flow due to elemental demixing was seen to not condition these findings.
The various competing effects in instability development linked to dissociation and ionization were identified and quantified for the three considered flight-envelope points:
(i) The destabilization due to the base-flow cooling induced by dissociation and ionization. This cooling also brings with it a strengthening of the thermoacoustic impedance and of the temperature perturbation amplitude near the boundary-layer edge.
(ii) The stabilization introduced by the diffusion fluxes acting on the perturbation terms. They were seen to act on the aforementioned peak in the temperature perturbation amplitude. This damping is believed to be the manifestation of the mechanism stabilizing the second mode.
(iii) A further stabilization linked to the dissociation and ionization source terms acting on the perturbation terms.
Not accounting for species diffusion, dissociation and ionization in the stability computations would imply neglecting the second and third effects, and would lead to a strong overestimation of instability growth in the considered test conditions.
Out of these three competing effects, it was the diffusion fluxes that were seen to make second-mode instabilities become supersonic in the investigated hot-wall scenarios. This contradicts the previously dominating belief that it was the dissociation-induced base-flow cooling that was actually originating supersonic modes in high-enthalpy flows. Previous studies did not decouple these two phenomena, and therefore could only speculate about the base-flow cooling being the underlying reason for it. This claim is also compatible with the common understanding of supersonic modes, which are promoted by wall cooling. It simply adds a caveat to it, in the sense that dissociation- or ionization-induced cooling of the base flow cannot be equated to wall cooling in how they affect supersonic modes. Whilst the latter promotes them, the former does not.
The modification of the predictions linked to using the linearized Rankine–Hugoniot relations on the free-stream perturbation values, rather than Dirichlet conditions, was also investigated. The ultimate modification of the 
$N$ factors was seen to be negligible. The fact that this was the case for shocks at such a small distance to the surface (1.3
$^{\circ }$–2.3
$^{\circ }$) suggests that the effect of the shock for such simple geometries may be disregarded, or that alternative, more physically accurate methodologies must be found to model the interaction between instabilities and a shock.
The results of this study provide a framework for the inclusion and assessment of other effects such as ablation and surface radiation.
Acknowledgements
The present work was supported by the Belgian FNRS through the FRIA fellowship, and by the US Air Force Office of Scientific Research (AFOSR) through the National Defense Science & Engineering Graduate (NDSEG) Fellowship Program. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of AFOSR. The authors also want to thank Dr M. Hermanns for providing the FDq library, Dr K. Groot for his contribution to the development of DEKAF's marching solver and Dr T. Magin and the rest of the mutation++ team for their valuable expert advice regarding thermophysical modelling.
Declaration of interests
The authors report no conflict of interest.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2020.786.
Appendix A. Thermophysical models
A.1. Transport model: Chapman and Enskog
The most accurate transport model consists of the first (for viscosity) and second (for the heavy-particle translational thermal conductivity 
$\kappa _H^{Trans}$) approximation to the CE molecular theory of gases using Laguerre–Sonine polynomials (see Chapman & Cowling Reference Chapman and Cowling1939; Magin & Degrez Reference Magin and Degrez2004):
\begin{equation} Q = - \left|\begin{matrix} G^Q_{s \ell} & X_s \\ X_\ell & 0 \end{matrix} \right| \bigg/ |G^Q_{s \ell}|,\quad \forall\ s,\ell \in \mathcal{H}, \end{equation}
where 
$\mathcal {H}$ is the set of all heavy species (excluding electrons) and 
$Q$ can be the viscosity 
$\mu$ or the heavy-particle translational thermal conductivity 
$\kappa ^{Trans}_H$. The elements of the matrix subsystems 
$G^Q_{s \ell }$ are detailed in Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) and Miró Miró (Reference Miró Miró2020) with the same nomenclature of this text, detailed in the supplementary material.
The other addends forming the frozen thermal conductivity are the thermal conductivity due to the internal energy modes (
$\kappa ^{Rot}$, 
$\kappa ^{Vib}$ and 
$\kappa ^{Elec}$) and, for ionized mixtures, the electron thermal conductivity (
$\kappa ^{Trans}_{{\text {e}^{{-}}}}$). Expressions for the internal-energy-mode conductivity are obtained from Eucken's relation (see Eucken Reference Eucken1913; Magin & Degrez Reference Magin and Degrez2004):
\begin{gather} \kappa^{Rot} = \sum_{s \in \mathcal{S}_{mol}} \frac{\mathscr{M}_s}{N_A} \frac{X_s c_{v s}^{Rot}}{\displaystyle\sum_{\ell \in \mathcal{H}} X_\ell/n \mathscr{D}_{s \ell}}, \end{gather}
\begin{gather}\kappa^{Vib} = \sum_{s \in \mathcal{S}_{mol}} \frac{\mathscr{M}_s}{N_A} \frac{X_s c_{v s}^{Vib}}{\displaystyle\sum_{\ell \in \mathcal{H}} X_\ell/n \mathscr{D}_{s \ell}}, \end{gather}
\begin{gather}\kappa^{Elec} = \sum_{s\in \mathcal{H}} \frac{\mathscr{M}_s}{N_A} \frac{X_s c_{v s}^{Elec}}{\displaystyle\sum_{\ell \in \mathcal{H}} X_\ell/n \mathscr{D}_{s \ell}}, \end{gather}
where 
$N_A$ is Avogadro's number, 
$\mathcal {S}_{mol}$ is the set of molecular species, 
$c_{v s}^{Mod}$ is the species heat capacity at constant volume due to each energy mode, 
$n$ is the mixture number density and 
$\mathscr {D}_{s\ell }$ is the binary diffusion coefficient of the species pair 
$s$–
$\ell$. The expressions for the product 
$n \mathscr {D}_{s\ell }$ can be found in Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) or Miró Miró (Reference Miró Miró2020) with the same nomenclature of this text, detailed in the supplementary material.
The electron thermal conductivity 
$\kappa ^{Trans}_{\text {e}^{{-}}}$ can be obtained from the second approximation to the theory of Chapman and Enskog (see Chapman & Cowling Reference Chapman and Cowling1939; Magin & Degrez Reference Magin and Degrez2004):
\begin{equation} \kappa^{Trans}_{\text{e}^{{-}}} = \frac{X_{\text{e}^{{-}}}^2}{\varLambda_{{\text{e}^{{-}}}{\text{e}^{{-}}}}^{11}}, \end{equation}
where the polynomial 
$\varLambda _{{\text {e}^{{-}}}{\text {e}^{{-}}}}^{11}$ can be found in Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) or Miró Miró (Reference Miró Miró2020) with the same nomenclature of this text, detailed in the supplementary material. The frozen thermal conductivity is thus
\begin{equation} \kappa^{Fr} = \kappa^{Trans}_H + \kappa^{Rot} + \kappa^{Vib} + \kappa^{Elec} + \kappa^{Trans}_{{\text{e}^{{-}}}}, \end{equation}
where, for non-ionized mixtures, the electron conductivity is obviously zero: 
$\kappa ^{Trans}_{{\text {e}^{{-}}}}=0$.
The Stokes hypothesis applies as long as inelastic collisions are neglected (see Ferziger & Kaper Reference Ferziger and Kaper1972; Bertolotti Reference Bertolotti1998; Giovangigli Reference Giovangigli1999):
\begin{equation} \lambda = -\tfrac{2}{3} \mu.\end{equation}A.2. Diffusion model: SCEBAD
One can accurately model the species diffusion fluxes by defining effective diffusion coefficients, similiarly to Yos (Reference Yos1963) or Fertig, Dohr & Frühauf (Reference Fertig, Dohr and Frühauf2001):
\begin{equation} \mathcal{D}_{{eff}\,s} = \frac{1-X_s}{\displaystyle\sum_{\substack{\ell \in \mathcal{S} \\ \ell \neq s}} X_\ell / \mathscr{D}_{s \ell}},\quad \forall\ s\in\mathcal{S}, \end{equation}and applying the modification proposed by Lee (Reference Lee1985) for ambipolar diffusion such as that assumed in this work:
where 
$\mathcal {S}_{ion}$ is the set of all ion species. One can enforce the self-consistency of the ambipolar diffusion fluxes with the correction proposed by Ramshaw & Chang (Reference Ramshaw and Chang1993), leading to the self-consistent effective binary ambipolar diffusion model (SCEBAD):
\begin{gather} \mathcal{D}^{ambip.}_{s \ell} = \left(\delta_{s \ell} - Y_\ell \left(1+\mathcal{Z}_\ell \frac{\mathscr{M}_{{\text{e}^{{-}}}}}{\mathscr{M}_\ell}\right)\right) \frac{1}{X_\ell}\frac{1-Y_\ell}{1-X_\ell} \left(1-\delta_{{\text{e}^{{-}}}\ell}\right) \mathcal{D}^{ambip.}_{{eff},\ell},\quad \forall\ s,\ell \in \mathcal{H}, \end{gather}
\begin{gather}\mathcal{D}^{ambip.}_{{\text{e}^{{-}}} \ell} = \left(\left(1-Y_{{\text{e}^{{-}}}}\right) \mathcal{Z}_\ell -Y_{{\text{e}^{{-}}}} \frac{\mathscr{M}_{\ell}}{\mathscr{M}_{{\text{e}^{{-}}}}} \right) \frac{1}{X_{{\text{e}^{{-}}}}}\frac{1-Y_\ell}{1-X_\ell} \left(1-\delta_{{\text{e}^{{-}}}\ell}\right) \mathcal{D}^{ambip.}_{{eff},\ell},\quad \forall\ \ell \in \mathcal{H}. \end{gather}One must bear in mind that Ramshaw's correction implicitly assumes that the diffusion fluxes of heavy species is much larger than that of electrons, allowing for a simplified expression of the electric field (Magin & Degrez Reference Magin and Degrez2004).
A.3. Collision integrals
Certain transport models require expressions for the species pairs’ collisional cross-sectional integrals. An excellent review of the physical interpretation of these integrals is given in § 7.4 of Hirschfelder, Curtiss & Bird (Reference Hirschfelder, Curtiss and Bird1954). Data for these collision integrals normally come from experiments or from computations based on the particle intermolecular force potentials. They are normally presented in tables as a function of temperature, or, for charged collisions, the reduced temperature 
$T^*$:
\begin{equation} T^* = \frac{\lambda_D k_B T}{{q}_{\text{e}^{{-}}}^2/\left(4{\rm \pi}\epsilon_0\right)},\quad \forall\ s,\ell\in\mathcal{S}, \end{equation}
where 
$k_B$ is Boltzmann's constant, 
$ {q}_{\text {e}^{{-}}}$ is the elementary charge, 
$\epsilon _0$ is permittivity of vacuum and 
$\lambda _D$ is the Debye shielding length:
\begin{equation} \lambda_D = \sqrt{\frac{\epsilon_0 k_B / {q}_{\text{e}^{{-}}}^2}{\displaystyle\sum_{s\in\mathcal{S}} \mathcal{Z}_s^2 n_s/T}}, \end{equation}
where 
$n_s$ is the species number density. There exist alternative definitions of the Debye shielding length, since it is not clear whether ions should contribute to the shielding (like in (A 10)) or not (see Scoggins Reference Scoggins2017).
When performing computational fluid dynamics simulations, it is common practice to simply interpolate between these tabled values. However, for stability analyses one needs analytical derivatives of the collision integrals with respect to temperature, prompting the creation of polynomial fits. The master expression is a higher-polynomial-order version of the one proposed by Gupta et al. (Reference Gupta, Yos, Thompson and Lee1990):
\begin{align} \ln\left(\left(\frac{T^*}{\lambda_D}\right)^2 \varOmega_{s \ell}^{(i,j)}\right) &= A_{\varOmega_{s \ell}^{(i,j)}} + B_{\varOmega_{s \ell}^{(i,j)}}\ln\left(T^*\right) + C_{\varOmega_{s \ell}^{(i,j)}}\left(\ln\left(T^*\right)\right)^2 \nonumber\\ &\quad +D_{\varOmega_{s \ell}^{(i,j)}}\left(\ln\left(T^*\right)\right)^3 + E_{\varOmega_{s \ell}^{(i,j)}}\left(\ln\left(T^*\right)\right)^4 \nonumber\\ &\quad +F_{\varOmega_{s \ell}^{(i,j)}}\left(\ln\left(T^*\right)\right)^5,\quad \forall\ s,\ell\in\mathcal{S}_{ion} \cup {\text{e}^{{-}}}, \end{align}where (
A 11a) applies to collisions between neutral species or between a neutral and a charged species, and (A 11b) applies to collisions between two charged particles. Note that the condition on the species subscripts in (
A 11a) (
) implies that either 
$s$ or 
$\ell$ must not be in the group of charged species (
$\mathcal {S}_{ion} \cup {\text {e}^{{-}}}$).
Oftentimes authors report the ratios of collision integrals 
$B^*_{s\ell }$, 
$C^*_{s\ell }$:
\begin{gather} B^*_{s\ell} =\frac{5\varOmega_{s \ell}^{(1,2)}-4\varOmega_{s \ell}^{(1,3)}}{\varOmega_{s \ell}^{(1,1)}},\quad \forall\ s,\ell \in\mathcal{S}, \end{gather}
\begin{gather}C^*_{s\ell} =\frac{\varOmega_{s \ell}^{(1,2)}}{\varOmega_{s \ell}^{(1,1)}},\quad \forall\ s,\ell \in\mathcal{S}, \end{gather}rather than the collision integrals themselves. For this reason, novel curve fits are also employed for these collision integral ratios, similarly to Gupta et al. (Reference Gupta, Yos, Thompson and Lee1990):
\begin{align} \ln\left(B^*_{s \ell}\right) &= A_{B^*_{s\ell}} + B_{B^*_{s\ell}}\ln\left(T^*\right) + C_{B^*_{s\ell}} \left(\ln\left(T^*\right)\right)^2 \nonumber\\ &\quad + D_{B^*_{s\ell}}\left(\ln\left(T^*\right)\right)^3 + E_{B^*_{s\ell}} \left(\ln\left(T^*\right)\right)^4 + F_{B^*_{s\ell}}\left(\ln\left(T^*\right)\right)^5,\quad \forall\ s,\ell \in\mathcal{S}_\text{ion} \cup {\text{e}^{{-}}}. \end{align}
Analogous ones are used for 
$C^*_{s\ell }$.
In order to obtain the coefficients 
$A$–
$F$ in (A 11) and (A 13), they are fitted to the state-of-the-art collisional data presented by:
(i) Wright et al. (Reference Wright, Bose, Palmer and Levin2005) and Wright, Hwang & Schwenke (Reference Wright, Hwang and Schwenke2007) – neutral–neutral and electron–neutral collisions;
(ii) Levin & Wright (Reference Levin and Wright2004) – ion–neutral collisions; and
(iii) Mason (Reference Mason1967) and Devoto (Reference Devoto1973) – charged–charged collisions.
A summary of the methodology used to compute the various collision integrals is presented in table 6.
Table 6. Methodology to compute the various collision integrals appearing in the transport models.
The evaluation of the curve fits outside of the tabled range of temperatures could lead to extrapolation problems. In order to avoid this, additional fictitious points are added, if necessary, to the original set of data points before performing the fitting. Depending on the dataset, the addition points are either ‘clipped off’ (zeroth-order extrapolation) or extrapolated linearly (first-order extrapolation) from the last data points. This is done to preserve eventual clear trends in the data, and to avoid spurious oscillations in the fitting due to radical changes in these trends. The curve fits of all collisions together with the original tabled data can be found in Miró Miró (Reference Miró Miró2020).
A.4. Thermal model: rigid rotor and harmonic oscillator
When a gas remains below its vibrational activation temperature it can be considered a CPG. In such instances, the enthalpy is linearly dependent on the temperature, since only the particle translational and rotational energies are considered:
\begin{equation} h = c_p T, \end{equation}
where the heat capacity (
$c_p$) is a constant that, for air, is commonly taken equal to 
$1004.5\ \textrm{J}\ \textrm{kg}^{-1}$ K
$^{-1}$.
For all other flow assumptions, one must approximate the nonlinear functional dependency of the vibrational and electronic energy modes on temperature. Such expressions are obtained from differentiating the partition functions of the different energy modes assuming molecules behave like a rigid rotor and a harmonic oscillator, and assuming species to populate the electronic energy levels according to a Boltzmann distribution (see Scoggins & Magin Reference Scoggins and Magin2014; Scoggins Reference Scoggins2017; Anderson Reference Anderson2006):
\begin{gather} e_s^{Mod} = \frac{\mathscr{R}}{\mathscr{M}_s} T^2 \frac{\partial \ln \mathcal{Q}^{Mod}_{ s}}{\partial T},\quad \forall\ s\in\mathcal{S}, \end{gather}
\begin{gather}c_{vs}^{Mod} = \frac{\mathscr{R}}{\mathscr{M}_s} T \left( 2 \frac{\partial \ln \mathcal{Q}^{Mod}_{ s}}{\partial T} + T \frac{\partial^2 \ln \mathcal{Q}^{Mod}_{ s}}{\partial T^2} \right),\quad \forall\ s\in\mathcal{S}, \end{gather}
\begin{gather} \mathcal{Q}^{TransV}_{ s} = \left( \frac{2 {\rm \pi}\mathscr{M}_s k_B T}{N_A \hbar^2} \right)^{{3}/{2}},\quad \forall\ s\in\mathcal{S}, \end{gather}
\begin{gather}\mathcal{Q}^{Rot}_{ s} = \frac{1}{\sigma_s} \left(\frac{T}{\theta_s^{Rot}}\right)^{\mathcal{L}_s/2},\quad \forall\ s \in \mathcal{S}_{mol}, \end{gather}
\begin{gather}\mathcal{Q}^{Vib}_{ s} = \sum_{m=1}^{N_{vib}} \frac{g^{Vib}_{sm}}{1-\mathrm{e}^{-\theta^{Vib}_{sm} / T}},\quad \forall\ s \in \mathcal{S}_{mol}, \end{gather}
\begin{gather}\mathcal{Q}^{Elec}_{ s} = \sum_{m=0}^{N_{elec}} g^{Elec}_{s m}\,\mathrm{e}^{-\theta^{Elec}_{s m} / T},\quad \forall\ s \in \mathcal{H}, \end{gather}
\begin{gather}\mathcal{Q}^{Int}_{{\text{e}^{{-}}}} = 2, \end{gather}
\begin{gather} h_s = \sum_{Mod} e_s^{Mod} + h_{fs}^\circ + \frac{\mathscr{R}}{\mathscr{M}_s} T,\quad \forall\ s\in\mathcal{S}, \end{gather}
\begin{gather}c_{ps} = \sum_{Mod} c_{vs}^{Mod} + \frac{\mathscr{R}}{\mathscr{M}_s},\quad \forall\ s\in\mathcal{S}, \end{gather}
where 
$\mathscr {R}$ is the universal gas constant, 
$\mathscr {M}_s$ is the species molar mass, 
$\mathcal {Q}_s^{Mod}$ is the partition function of each energy mode, 
$\sigma _s$ is the molecule steric factor (2 for symmetric and 1 for non-symmetric), 
$\mathcal {L}_s$ is the molecule linearity factor (3 for nonlinear and 2 for linear), 
$\theta _s^{Rot}$ is the rotational activation temperature, 
$\theta ^{Vib}_{sm}$ is the activation temperature of the 
$m$th vibrational mode of species 
$s$, 
$\theta ^{Elec}_{s m}$ is the activation temperature of the 
$m$th electronic energy level of species 
$s$, 
$g^{Elec}_{s m}$ is the degeneracy of the 
$m$th electronic energy level of species 
$s$, 
$g^{Vib}_{s m}$ is the degeneracy of the 
$m$th vibrational mode of species 
$s$, 
$h_{fs}^\circ$ is the species formation enthalpy at 0 K, 
$e_s^{Mod}$ is the internal energy of each species due to each energy mode and 
$c_{p s}$ is the species heat capacity at constant pressure.
Equation (A 17a) corresponds to the volumetric partition function of the translational energy, since it has been divided by the system's volume. The expression for the partition function of the electron species’ internal energy 
$\mathcal {Q}^{Int}_{{\text {e}^{{-}}}}$ (A 17e) is formulated to account for the contribution of the spin of the free electron to the species entropy (see Scoggins & Magin Reference Scoggins and Magin2014; Scoggins Reference Scoggins2017). Note that the expression of the electronic partition function (A 17d) features a summation between zero and the number of electronic modes. In reality there are infinite electronic modes, but the summation must be truncated to avoid a divergence of the degeneracies (see Scoggins Reference Scoggins2017; Bellas-Chatzigeorgis Reference Bellas-Chatzigeorgis2018).
The species properties appearing in these expressions (
$\sigma _s$, 
$\mathcal {L}_s$, 
$\theta ^{Rot}_s$, 
$g^{Vib}_{s m}$, 
$\theta ^{Vib}_{s m}$, 
$\theta ^{Elec}_{s m}$ and 
$g_{s m}$) can be found for instance in Gurvich, Veyts & Alcock (Reference Gurvich, Veyts and Alcock1989) or summarized in table 7 (in SI units) for the species considered in this work. The species formation enthalpy per mole at 298 K (
$H_{fs}^{298\ \mathrm {K}}$) is by convention zero for a series of reference species (N
$_{\text {2}}$, O
$_{\text {2}}$, e−, C(gr), He, Ne, Ar, etc.), and for all other species it is obtained from a formation reaction. The species formation enthalpies at 0 K can then be obtained using Hess’ law.
Table 7. Thermal properties of the molecular species considered in this work. Obtained from the mutation++ database (see Magin & Degrez Reference Magin and Degrez2004; Scoggins & Magin Reference Scoggins and Magin2014; Scoggins Reference Scoggins2017), which gathers values originally in Gurvich et al. (Reference Gurvich, Veyts and Alcock1989). The degeneracies and activation temperatures of the various electronic energy levels considered were also obtained from the mutation++ database, and are collected in Miró Miró et al. (Reference Miró Miró, Beyak, Mullen, Pinna and Reed2018a) and Miró Miró (Reference Miró Miró2020).
The mixture enthalpy and heat capacity are obtained by summing over the species quantities weighted with the mass fractions:
\begin{gather} h = \sum_{s \in \mathcal{S}} Y_s h_s, \end{gather}
\begin{gather}c_p = \sum_{s \in \mathcal{S}} Y_s c_{ps}. \end{gather}A.5. Chemical model
For a set of reactions 
$\mathcal {R}$ between a set of species 
$\mathcal {S}$, with reactant stoichiometries 
$\nu _{sr}'$ and product stoichiometries 
$\nu _{sr}''$,
\begin{equation} \sum_{s \in \mathcal{S}} \nu'_{sr} s \leftrightarrow \sum_{s \in \mathcal{S}} \nu''_{sr} s,\quad \forall\ r \in \mathcal{R}, \end{equation}the mass production rate of each species can be approximated by the law of mass action (see Vincenti & Kruger Reference Vincenti and Kruger1967):
\begin{equation} \dot{\omega}_s = \mathscr{M}_s \sum_{r \in \mathcal{R}} (\nu''_{sr} -\nu'_{sr})\left(k_{f r} \prod_{\ell \in \mathcal{S}} \left(\frac{\rho_\ell}{\mathscr{M}_\ell} \right)^{\nu'_{\ell r}} - k_{b r} \prod_{\ell \in \mathcal{S}} \left(\frac{\rho_\ell}{\mathscr{M}_\ell} \right)^{\nu''_{\ell r}} \right),\quad \forall\ s\in\mathcal{S}, \end{equation}
where 
$k_{fr}$ and 
$k_{br}$ are the forward and backward reaction rates, and can be defined as
\begin{gather} k_{f r} = A_r^f T^{n_{Tr}^f} \,\mathrm{e}^{-\theta_r^f/T},\quad \forall\ r\in\mathcal{R}, \end{gather}
\begin{gather}k_{b r} = \frac{k_{f r}}{K_{{eq}\,r}},\quad \forall\ r\in\mathcal{R}. \end{gather} The reaction rate constants appearing in (A 24a) are those detailed in table 8. The equilibrium constant for each reaction 
$K_{{eq}\,r}$ can be obtained, using kinetic theory, from the partition functions of the different energy modes. Assuming a rigid rotor and harmonic oscillator model for the species (see § A.4):
\begin{equation} K_{{eq}\,r} = \prod_{s \in \mathcal{S}} \left(\frac{\displaystyle\prod_{Mod} \mathcal{Q}_{ s}^{Mod}}{N_A}\right)^{\nu''_{s r} - \nu'_{s r}} \cdot \exp\left(-(\nu''_{s r} - \nu'_{s r})\frac{h_{f,s}^\circ \mathscr{M}_s}{\mathscr{R} T}\right),\quad \forall\ r\in\mathcal{R}. \end{equation}Note that in order to avoid floating-point overflow, it is desirable to evaluate these expressions in their logarithmic form.
Table 8. Forward reaction rate constants. All data are from Park (Reference Park1993), except for the neutral exchange reactions that are from Bose & Candler (Reference Bose and Candler1996, Reference Bose and Candler1997).
A.6. Thermodynamic derivatives
The resolution of the stability equations requires an expansion around zero of the perturbations of all the thermophysical properties. As a consequence of this Taylor expansion, one must also compute the derivatives of all the expressions modelling the various properties presented in this section, with respect to the thermodynamic state quantities on which they are functionally depending.
Such an endeavour is significantly error-prone, thus justifying the use of computer algebra systems that can symbolically operate and differentiate the various expressions, and implement them into executable functions. To that end, an additional module was developed within VESTA's automatic derivation and implementation tools (see § 4.6 of Miró Miró (Reference Miró Miró2020)).
