# Acoustic stability of nonadiabatic high-energy-density shocks

César Huete, Francisco Cobos-Campos, Ernazar Abdikamalov, Serge Bouquet

► To cite this version:

César Huete, Francisco Cobos-Campos, Ernazar Abdikamalov, Serge Bouquet. Acoustic stability of nonadiabatic high-energy-density shocks. Physical Review Fluids, 2020, 5 (11), 10.1103/PhysRevFluids.5.113403. hal-03209540

HAL Id: hal-03209540

https://hal.science/hal-03209540v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
This is a postprint version of the following published document:

Huete, C., Cobos-Campos, F., Abdikamalov, E., & Bouquet, S. (2020). Acoustic stability of nonadiabatic high-energy-density shocks. *In Physical Review Fluids, 5*(11)

DOI:10.1103/physrevfluids.5.113403

©2020 American Physical Society
Acoustic stability of non-adiabatic high-energy-density shocks

Csar Huete

Fluid Mechanics Group, Escuela Politécnica Superior,
Universidad Carlos III de Madrid, SP*

Francisco Cobos-Campos

Instituto de Investigaciones Energéticas y Aplicaciones Industriales,
Universidad de Castilla-La Mancha, 13071 Ciudad Real, SP

Ernazar Abdikamalov

Department of Physics, Nazarbayev University, Nur-Sultan 010000, KA

Serge Bouquet

CEA/DAM/DIF, Bruyres-le-Chtel, F91297 Arpajon Cedex, FR

(Dated: June 17, 2020)
Abstract

As firstly predicted by D’yakov and Kontorovich (DK), an initially disturbed shock front may exhibit different asymptotic behaviours depending on the slope of the Rankine-Hugoniot curve. Adiabatic and isolated planar shocks traveling steadily through ideal gases are stable, in the sense that any perturbation on the shock front decays in time with the power $t^{-3/2}$ (or $t^{-1/2}$ in the strong-shock limit). While some gases whose equation of state cannot be modelled as a perfect gas, as those governed by van der Waals forces, may induce constant-amplitude oscillations to the shock front in the long-time regime, fully unstable behaviours are seldom to occur due to the unlikely conditions that the equation of state must meet. In this work, it has been found that unstable conditions are might be found when the gas undergoes an endothermic or exothermic transformation behind the shock. In particular, it is reported that constant-amplitude oscillations can occur when the amount of energy release is positively-correlated to the shock strength and, if this correlation is sufficiently strong, the shock turns to be fully unstable. The opposite highly-damped oscillating regime may occur in negatively-correlated configurations. The mathematical description then adds two independent parameter to the regular adiabatic index $\gamma$ and shock Mach number $M_1$, namely: the total energy added/removed and its sensitivity with the shock strength. The formulation in terms of endothermic or exothermic effects is extended, but not restricted, to include effects associated to ionization, dissociation, thermal radiation, and thermonuclear transformations, so long as the time associated to these effects is much shorter time than the acoustic time.

<footnote>* chuete@ing.uc3m.es</footnote>
# I. INTRODUCTION

Motivated by the dominant role of shock waves in different areas, as may be high-speed aerodynamics and inertial confinement fusion, there exists an exhaustive theoretical work accumulated on shock dynamics [1–19] since the pioneering works of D’yakov [20] and Kontorovich [21]. The literature accumulated on the Richtmyer-Meshkov instability [22–30] deserves special recognition. By way of contrast, non-adiabatic waves have received considerably less attention, with regular detonations standing among any other type of supersonic reactive front due to its direct application in propulsion engines and safety issues [31–38]. Further examples of non-adiabatic supersonic fronts are found in shocks that can induce a phase change on their path [39–42], ionization [6, 7], dissociation [6, 8] and radiation [43].

In the astrophysical context, non-adiabatic shocks can develop in a wide range of phenomena, including thermonuclear detonations formed in thermonuclear supernovae [44–48]. Further examples can be found in the so-called core-collapse supernovae (CCSNe), a type of supernovae powered by stellar core collapse. The shock wave evolves through a phase of accretion shock, in which it becomes endothermic due to the nuclear dissociation process taking place across the shock [49–52]. This phase is pivotal in determining the parameters of the supernova explosion [53–57]. A distinctive common feature of these two fronts is that energy variations (either exothermic and endothermic) may likely depend on thermodynamical absolute values, and thus on the shock intensity. A clear example is found in the energy lost by radiation in the so-called radiative shocks [43, 58–64].

In most of the cases, based on the different scales associated to the shock thickness and the post-shock transformation zone, the phenomena occurring in the supersonic front can be split into the precursor shock wave and the attached reacting layer. For situations where the transformation right behind the shock does not modify the mass and momentum fluxes, the Rayleigh-Michelson line remains the same regardless the non-adiabatic transformation occurring behind the precursor shock. Then, the final equilibrium state must lay on some point along this line. For the case addressed here, attention will be paid in the global structure only, where net changes (shock plus non-adiabatic transformation) are included in the “non-adiabatic shock”, thereby being considered a pure discontinuity. Classical D’yakov-Kontorovich (DK) theory states that, in certain conditions related to the Rankine-Hugoniot (RH) slope, the shock oscillates permanently with the corresponding generation of entropic
and rotational perturbations downstream, as well as the emission of constant-amplitude sonic waves, the latter being commonly referred as spontaneous acoustic emission (SAE). Beyond certain limit the shock may behave in a fully unstable manner, with any perturbation growing exponentially with time. This long-time dynamics, which can be investigated via normal-mode analysis, will be studied here by means of Laplace transform method as it provides information of the transient evolution. The function that governs the shock dynamics is the DK parameter, here written as

$$ \Gamma_s = \frac{p_2 - p_1}{V_1 - V_2} \left. \frac{dV}{dp} \right|_2 = \frac{\tan \beta_{RM}}{\tan \beta_{RH}}, \quad (1) $$

directly relates with the slopes of the RH curve in the final state, $\tan \beta_{RH}$, and the slope of Rayleigh-Michelson line, $\tan \beta_{RM}$ (see sketch in Fig.1). The subscripts 2 and 1 in (1) identify downstream and upstream properties, respectively, and the functions $p$ and $V$ define the dimensional pressure and the specific volume, respectively.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Sketch of the Rankine-Hugoniot curve and the Rayleigh-Michelson line for a non-adiabatic shock. Different processes will place the final post-shock state (orange point) in a different position.</figcaption></figure>

When a shock moves steadily and it is isolated from piston-like coupling effects, linear theory dictates that any perturbation to the isolated shock surface decays on a long time scale in accordance with the power law $t^{-3/2}$ [2]. However, shock moving through gases that obey non-ideal equation of state can oscillate with constant amplitude, as shown by Bates [9] and Wouchuk [12], who predicted the possibility of DK-instability in gases governed by van
der Waals forces. Further examples are associated to ionization [6–8] and relativistic effects [18]. In a recent work [38] it was shown that permanent oscillations and SAE can also occur in strength-sensitive detonations while highly-damped oscillations have been also associated to non-ideal gases [11] and shocks that induce endothermic transformation in the gas state [52]. In this work, as well as in the references mentioned in this paragraph, the base-flow is taken uniform and the shock isolated. Consequently, the instability thresholds should be adequately modified in more realistic conditions resulting from unsteadiness [65, 66] and/or acoustic coupling with the supporter mechanism needed in non-exothermic scenarios [16].

The objective of this work is to extend, with a unified formulation, the analysis of the stability of non-adiabatic shocks traveling in a perfect gas when endothermic or exothermic transformations take place. The formulation allows writing other effects, as might be molecule structure changes, in terms of equivalent non-adiabatic effects. The resolution, presented is a self-contained form, is based on the Laplace transform to facilitate the analysis for both long-time and transient regimes. The methodology employs linear perturbation analysis in the thin-shock limit, where the perturbation wavelength $2\pi/k$ is much larger than the shock thickness $\ell$, which comprises the precursor adiabatic shock and the following non-adiabatic transformation layer. Since the DK instability parameters ultimately depend on the RH curve slope in the final state, the analysis can be simplified by accounting the shock strength (propagating Mach number), adiabatic index, energy added/subtracted per unit mass (either positive or negative) and the variation of energy release/taken with the shock strength. That is, results do not depend on the inner structure profiles but on the final equilibrium conditions. Hereafter the term non-adiabatic will comprise any type of energy gain or loss across the shock, e.g., reactive, dissociating, radiating or ionizing.

For the reader's convenience, the outline of the paper is stated in the following. First, the formulation of the problem is presented in section II in terms of one-dimensional steady properties. The effect of energy variations across the shock and the limits of interest are discussed in detail. It follows the mathematical details of the linear model in terms of the Laplace transform. Both long-time asymptotic expressions and the transient formulae are given for the shock front evolution. Theory is later applied, in section III, to cases of interest in high-energy-density physics. They include ionization, dissociation, radiative and nuclear-transformation shocks. A short summary of outcomes is offered in section IV. Appendices A-D provide the derivation of the non-adiabatic phenomena used in section III.
## II. LINEAR THEORY FOR THE SHOCK FRONT EVOLUTION

### A. Base flow equations

It is defined an isolated shock whose relative speed with respect to the upstream and downstream flow is $u_1$ and $u_2$, respectively. The subscripts 1 and 2 will refer to upstream and downstream flow variables, respectively, which are identified with a prime in dimensional form. They include: velocity $u$, density $\rho$, pressure $p$, and enthalpy $h$. The conservation equations read

$$ \rho_1 u_1 = \rho_2 u_2, \qquad (2) $$

$$ p_1 + \rho_1 u_1^2 = p_2 + \rho_2 u_2^2, \qquad (3) $$

$$ h_1 + \frac{u_1^2}{2} + q = h_2 + \frac{u_2^2}{2}, \qquad (4) $$

where the factor $q$ defines the energy per unit mass delivered to the fluid particles in case of considering exothermic effects ($q > 0$) or the energy per unit mass subtracted to the fluid particles when endothermic effects are considered ($q < 0$). In general, there exists the possibility of having three different functions for the first adiabatic exponent $\Gamma_1 = (\partial \ln p)/(\partial \ln \rho)_{ad}$, the ratio of the specific heats $\gamma = c_p/c_v$, and the enthalpy coefficient $g$, as occurs in partially dissociated or ionized gases [67, 68]. The definition of the latter $g = h/e = h/(h - p/\rho)$, where $h$ and $e$ are the specific enthalpy and internal energy, respectively. It allows the energy equation be written

$$ \frac{g_1}{g_1 - 1} \frac{p_1}{\rho_1} + \frac{u_1^2}{2} + q = \frac{g_2}{g_2 - 1} \frac{p_2}{\rho_2} + \frac{u_2^2}{2}, \qquad (5) $$

which renders

$$ \left( \frac{\rho_2}{\rho_1} \right)^{-1} = \frac{u_2}{u_1} = \frac{1 + g_1^{-1} M_1^{-2}}{1 + g_2^{-1}} \left( 1 - \sqrt{1 - \kappa} \right) \qquad (6) $$

as the specific volume variation of the fluid particles, where

$$ \kappa = \frac{(1 - g_2^{-2}) [1 + 2 (g_1 - 1)^{-1} M_1^{-2}]}{(1 + g_1^{-1} M_1^{-2})^2} + \frac{2(1 - g_2^{-2}) M_1^{-2}}{(1 + g_1^{-1} M_1^{-2})^2} \frac{q}{g_1 p_1 / \rho_1}, \qquad (7) $$

is an auxiliary function and where $M_1 = u_1/(g_1 p_1 / \rho_1)^{1/2}$ is the scaled propagation speed that reduces to the Mach number for $\Gamma_1 = g_1 = \gamma_1$.

From (6), it is immediate to see that $\kappa = 1$ is the maximum energy gained for a given propagation speed, akin to Chapman-Jouget condition in detonation waves, and $\kappa = 0$
is the maximum energy taken from the bulk for a given propagation speed that leaves a quasi-static flow (infinite density) downstream. While the former is an achievable limit, the latter is just a theoretical bound that could be approached in extreme over-densification. Regular RH jump conditions are retrieved for $q = 0$ and $g_1 = g_2 = \Gamma_1 = \gamma$, which yields $\kappa = 1 - (1 + \gamma M_1^2)^{-2}$. Assuming known the mass compression ratio $R_s = \rho_2/\rho_1$ given in (6), the pressure ratio across the shock is determined by the Rayleigh-Michelson line

$$P_s = \frac{p_2}{p_1} = 1 + g_1 M_1^2 (1 - R_s^{-1}) \qquad (8)$$

that can be used to compute the effective post-shock Mach number

$$M_2 = \frac{u_2}{\sqrt{g_2 p_2 / \rho_2}} = \frac{M_1}{\sqrt{P_s R_s}} \sqrt{\frac{g_1}{g_2}} \qquad (9)$$

The temperature jump across the shock, written as a function of the density and pressure ratios, involves the corresponding equations of state, yielding

$$T_s = \frac{T_2}{T_1} = \frac{P_s}{R_s} \frac{1}{1 + \alpha}, \qquad (10)$$

where $\alpha$ refers to the relative variation of number of particles by ionization or dissociation phenomena (see next section for details) through $p_2 = (1 + \alpha)R_g T_2 \rho_2$, where $R_g$ is the gas constant in the upstream flow. Clearly, these equations call for complementary constitutive relations to couple the variation of the post-shock variables with the functions $q$, $g_2$ and $\alpha$.

Admitting that it may have some cost in the accuracy, the following simplification will allow writing the problem formulation in terms of much fewer parameters, with the corresponding gain in generality. For example, the effect associated to the change in $g$ and the corresponding endothermic or exothermic phenomena taking place across the shock are conveniently collected as

$$Q = \frac{\rho_1 (g_1^2 - 1)}{2g_1 p_1} q + \frac{(g_1 + 1)}{2g_1 (g_2 + 1)} \frac{P_s}{R_s} \Delta g, \qquad (11)$$

to compute the “effective” order-of-unity energy change. The term $\Delta g = g_2 - g_1$ is equivalent to an endothermic (exothermic) process when $g_2 < g_1$ ($g_2 > g_1$). This makes possible to get an effective adiabatic conditions ($Q = 0$) when existing exothermic transformations $q > 0$ but $\Delta g < 0$, and vice-versa. The above expression can be further simplified when the fluid can be modelled as a calorically perfect gas, $\Gamma_1 = \gamma = g$, since enthalpy reduces to $h(\gamma - 1)\rho = \gamma p$. The variable $\gamma$ is the adiabatic index associated to the slope of the adiabatic
relationship $p \sim \rho^{\gamma}$, that is used to define the speed of sound $a^2 = \gamma p/\rho$. When the upstream gas conditions do not differ significantly from standard conditions, such assumption is fully justified. Then, the deviation of the regular RH curve can be computed through the balance term $Q$ in the energy equation, which yield

$$ \mathcal{R}_s = \frac{\rho_2}{\rho_1} = \frac{u_1}{u_2} = \frac{(\gamma + 1) \mathcal{M}_1^2}{\gamma \mathcal{M}_1^2 + 1 - \sqrt{(\mathcal{M}_1^2 - 1)^2 - 4Q\mathcal{M}_1^2}} \quad (12) $$

for the mass compression ratio, from which the rest of the variables can be determined.

When the energy change across the shock does not depend on the shock properties, the above expressions can be used to set the limits of the shock conditions in both exothermic and endothermic scenarios. For instance, it is known that in absence of external influences the minimum propagation velocity of a exothermic shock (detonation) is the so-called Chapman-Jouget condition that occurs when $\mathcal{M}_2 = 1$, namely

$$ \mathcal{M}_1 = \mathcal{M}_{cj} = \sqrt{1 + \mathcal{Q}} + \sqrt{\mathcal{Q}}. \quad (13) $$

Alternatively, it can be written in terms of maximum energy release for a given propagation Mach number. Another limit appears in the opposite endothermic case, which is associated to an infinite mass compression ratio and a zero velocity $u_2 = 0$. Then the theoretical limits for the energy release are given by $Q_{\min} \le Q \le Q_{\max}$, where

$$ Q_{\min} = -\frac{\gamma + 1}{4} [(\gamma - 1)\mathcal{M}_1^2 + 2] \quad (14) $$

refers to minimum value of $Q$ (maximum energy subtracted) and

$$ Q_{\max} = \frac{(\mathcal{M}_1^2 - 1)^2}{4\mathcal{M}_1^2} \quad (15) $$

to the maximum value of $Q$ released to the fluid particles. These cases correspond to $\kappa = 1$ and $\kappa = 0$, respectively. Notice that $Q_{\min} < -(\gamma + 1)^2/4$ for $\mathcal{M}_1 > 1$.

When the energy variation across the shock depends on the shock intensity, the amount of energy delivered or taken from the fluid should be modeled according to the particular phenomenon taking place $Q = \text{fun}(\rho, T)$. It is found, for example, that accretion shocks formed in CCSNe are able to break heavy nuclei and the energy employed in this nuclear dissociation scales with the upstream energy flux [49, 50]. Further examples of the same sort can be found in supersonic fronts that induce phase change [39, 40]. Likewise, detonations moving through reactive mixtures made of air and sprays exhibit a shock-dependent behavior
as the amount of fuel actually available would depend on the shock evaporation capacity [41]. Thermonuclear detonations also depend on the shock intensity, as the amount of nuclei fused depends on temperature and density [44, 46, 47]. It is clear that $Q$ is a shock-dependent function when $\Delta g \neq 0$, even if the process is strictly adiabatic $Q = 0$.

In general, how energy variations scale with the flow properties will depend on the particular phenomenon, but the stability thresholds are given by the RH slope in the post-shock state. This reduces the parametric domain to three main gas properties (four if accounting for the adiabatic index $\gamma$): the shock Mach number $M_1$, the dimensionless energy release $Q$, and the local variation of the energy with the shock strength $\epsilon = (dQ)/(dM_1)$, defined as energy sensitivity parameter. Therefore, when the factor $q$ depends on local properties, as might be density and/or temperature (and pressure through the equation of state), post-shock values cannot be generally given in terms of independent expressions as a function of the shock strength, but rather as non-linear combination of them. For example, let's take temperature as the dominant parameter affecting the variation of energy in the form $q \sim Be^{-\phi/T_s}$, where $\phi$ stands for the dimensionless reference temperature below which the energy changes are negligible and $B$ is the potential of reference. With this law in mind, which is a reasonable assumption in many processes, as discussed later, final post-shock conditions are governed by an implicit relationship that cannot be reduced to any polynomial function. Once $T_s(M_1)$ is determined as a function of the shock Mach number, so is the function $q(M_1)$, and are the density ratio, pressure change and post-shock Mach number through (8)-(12). More complicated functions for $q$ are necessary to describe, as seen later, dissociation and ionization processes. Irrespective of the functional dependence that describes the variation of energy in the gas, it can ultimately be related (either numerically or analytically) to shock Mach number $M_1$, as well as the DK parameter $\Gamma_s$, conveniently rewritten as

$$ \Gamma_s = \frac{\gamma M_1^2}{R_s^2} \left(\frac{\partial P_s}{\partial R_s}\right)^{-1} = \frac{M_1^2(1+2Q) - 1 + \sqrt{(M_1^2-1)^2 - 4M_1^2Q - \epsilon M_1^3}}{M_1^2(M_1^2-1-2Q) + M_1^2\sqrt{(M_1^2-1)^2 - 4M_1^2Q - \epsilon M_1^3}}, \quad (16) $$

where the factor $\epsilon$ accounts for the sensitivity of the energy release with the shock strength. Notice that this function does not depend on $\gamma$ explicitly, as it is absorbed in the definition of $Q$. For a given shock with strength $M_1$ moving through a gas defined by the adiabatic index $\gamma$, the slope of Rayleigh-Michelson line is explicitly given by $\tan \beta_{RM} = \gamma M_1^2$, which
is always greater than unity. It is the energy changes what determine the slope of the RH equation, $\tan \beta_{\text{RH}} = \gamma M_1^2 / \Gamma_s$ and therefore, what settles the instability limits.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Left: DK parameter $\Gamma_s$ as a function of the dimensionless energy $Q$ for a strong shock $M_1 = 10$ and $\gamma = 1.4$, and for different values of energy sensitivities $\varepsilon$. Right: Qualitative representation of the RH angle $\beta_{\text{RH}}$.</figcaption></figure>

Figure 2 shows the DK parameter $\Gamma_s$ as a function of the dimensionless energy $Q$ for a strong shock ($M_1 = 10$) in air ($\gamma = 1.4$), and for different values of energy-change sensitivities $\varepsilon$. The angle $\beta_{\text{RH}}$ is qualitatively represented on the right hand side. It is observed that, ranging $\varepsilon$ between -10 and 10 suffices to find the most distinguished scenarios. To begin with the well-known adiabatic shock, the DK parameter for $Q = \varepsilon = 0$, $M_1 = 10$ and $\gamma = 1.4$ is $\Gamma_s = 0.01$, a small positive value that yields $\beta_{\text{RH}} < \pi/2$. This is readily seen as the crossing point with $\Gamma_s = 0$ (orange circles) occurs before. In the adiabatic strong-shock limit, $M_1 \gg 1$, the value of $\Gamma_s$ tends to zero and the slope approaches $\beta_{\text{RH}} \sim \pi/2$.

The effect of constant energy release can be evaluated from the $\varepsilon = 0$ curve. Exothermicity increases the value $\Gamma_s$ and decreases $\beta_{\text{RH}}$. On the other hand, exothermic shocks may exhibit a reversion in the RH curve, as the value of $\Gamma_s$ may become negative. This effect is amplified when $\varepsilon > 0$ and reduced for $\varepsilon < 0$. More over, the possibility of finding
a fully horizontal RH slope would occur for large values of $q$ and $\epsilon$. Endothermic shocks can also display a positive RH curve slope, as found in ionization and dissociation shocks [6–8, 10]. When taking $\epsilon$ as a purely independent parameter it is necessary to explore the values related to distinguished transformations behind the shock. For example, the value of $\beta_{RH}$ approaches $\pi/2$ when the numerator in (16) gets closer to zero, and $\beta_{RH} \sim \pi$ when the denominator gets closer to zero. For the Chapman-Jouget condition, or $Q = Q_{max}$, where the corresponding energy sensitivities are the same and equal to $(M_1 - M_1^{-3})/2$. For such conditions, both numerator and denominator approach zero with the same trend, thereby providing a finite value for the RH slope, which is widely known to be the slope of the Rayleigh-Michelson line $\beta_{RH} = \beta_{RM}$. This value critical value, which corresponds to $\sim M_1/2$ in the strong shock limit, is represented with a cyan circle in Fig.2.

## B. Shock transient evolution and stability limits

### 1. Shock transient evolution by the Laplace transform technique

To study the evolution of the planar shock, it is assumed an initial ripple of the form $\psi_s(y, t = 0) = \psi_{s0} \cos(ky)$, with $\psi_{s0}$ being the initial amplitude, $k$ its perturbation wavenumber, and $y$ the coordinate transverse to the shock direction, namely $x$. The wavenumber is employed to scale the spatial and temporal variables, $\hat{x} = kx$, $\hat{y} = ky$, and $\tau = a_2kt$. The functions of interest are linearized by expanding them around the small factor $\psi_{s0}k$ to yield

$$ \hat{p}(\hat{x}, \hat{y}, \tau) = \frac{1}{\psi_{s0}k} \frac{p(\hat{x}, \hat{y}, \tau) - p_2}{\rho_2 a_2^2}, \qquad (17) $$

$$ \hat{\rho}(\hat{x}, \hat{y}, \tau) = \frac{1}{\psi_{s0}k} \frac{\rho(\hat{x}, \hat{y}, \tau) - \rho_2}{\rho_2}, \qquad (18) $$

$$ \hat{u}(\hat{x}, \hat{y}, \tau) = \frac{1}{\psi_{s0}k} \frac{u(\hat{x}, \hat{y}, \tau) - (u_1 - u_2)}{a_2}, \qquad (19) $$

$$ \hat{v}(\hat{x}, \hat{y}, \tau) = \frac{1}{\psi_{s0}k} \frac{v(\hat{x}, \hat{y}, \tau)}{a_2} \qquad (20) $$

for pressure, density, longitudinal velocity and transverse velocity, respectively. The functions identified with the hat symbol are order-of-unity variables. Similarly, the dimensionless order-of-unity shock ripple is defined as $\xi_s = \psi_s/\psi_{s0}$.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Sketch of the perturbed shock wave and the perturbation variables in the shock reference frame, with $\ell \ll \psi_{s0} \ll 2\pi/k$.</figcaption></figure>

The dimensionless order-of-unity functions are used to write

$$ \frac{\partial \hat{\rho}}{\partial \tau} + \frac{\partial \hat{u}}{\partial x} + \hat{v} = 0, \qquad (21) $$

$$ \frac{\partial \hat{u}}{\partial \tau} + \frac{\partial \hat{p}}{\partial x} = 0, \qquad (22) $$

$$ \frac{\partial \hat{v}}{\partial \tau} - \hat{p} = 0, \qquad (23) $$

$$ \frac{\partial \hat{p}}{\partial \tau} - \frac{\partial \hat{\rho}}{\partial \tau} = 0 \qquad (24) $$

for the dimensionless conservation equations for mass, x-momentum, y-momentum, and energy, provided that $\hat{p}$ and $\hat{v}$ are always proportional to $\cos(y)$ and $\sin(y)$, respectively. One boundary condition is given by the values at the shock $\hat{\phi}(x = M_2\tau) = \hat{\phi}_s$, which are determined by linearized RH equations

$$ \frac{d\xi_s(\tau)}{d\tau} = \frac{\mathcal{R}_s}{\mathcal{R}_s - 1} \frac{1 - \Gamma_s}{2M_2} \hat{p}_s(\tau), \qquad (25) $$

$$ \hat{u}_s(\tau) = \frac{1 + \Gamma_s}{2M_2} \hat{p}_s(\tau), \qquad (26) $$

$$ \hat{\rho}_s(\tau) = \frac{\Gamma_s}{M_2^2} \hat{p}_s(\tau), \qquad (27) $$

$$ \hat{v}_s(\tau) = M_2(\mathcal{R}_s - 1) \xi_s(\tau), \qquad (28) $$

with use made of $\xi_s \sim \cos(y)$ in (28). The function $\Gamma_s$, associated to the slope of the RH curve, has been previously given in (16).

The other boundary condition is provided by the isolated-shock assumption, which translates into not considering the effect of the acoustic waves reaching the shock front from
behind. For this condition to be true, the shock must be sufficiently far from driving conditions, i.e., $\tau \gg M_2^{-1} - 1$, a condition met when $M_1/M_{cj} \sim 1$. Besides, the linear theory and the thin-shock assumptions set the following limits of validity: $\ell \ll \psi_{s0} \ll k^{-1}$.

Following the analysis shown in [38], omitted here for the sake of conciseness, the most relevant information relative to the asymptotic behavior can be inferred from the Laplace transform of the shock ripple amplitude, namely

$$ \Xi_s(s) = \int_0^\infty \xi_s(r) e^{-sr} dr = \frac{\sqrt{s^2+1} + \sigma_b s}{s\sqrt{s^2+1} + \sigma_b s^2 + \sigma_c} \quad (29) $$

applied over the variable $r = \tau\sqrt{1 - M_2^2}$, with the auxiliary factors defined as

$$ \sigma_b = \frac{1 + \Gamma_s}{2M_2} \quad \text{and} \quad \sigma_c = \frac{\mathcal{R}_s M_2}{1 - M_2^2} \frac{1 - \Gamma_s}{2}. \quad (30) $$

It is found that $\sigma_b < -1$ is associated to an exponential unstable growth of the shock ripple amplitude. For $-1 < \sigma_b < \sigma_c$ the shock oscillates with constant amplitude, with the associated emission of constant-amplitude acoustic waves (SAE). For values of $\sigma_b = \sigma_c$ the characteristic asymptotic decay rate of the shock oscillations amplitude is proportional to $\tau^{-1/2}$. For $\sigma_c < \sigma_b < \sigma_c + 1/(4\sigma_c)$ the shock front oscillates towards the asymptotic planar solution with an amplitude that decays in time with the power law $\tau^{-3/2}$, and for $\sigma_c + 1/(4\sigma_c) < \sigma_b < -1$ the approach towards the asymptotic decay rate occurs faster than that occurring in regular conditions, although with the same power law $\tau^{-3/2}$.

A simple way to evaluate the approach towards the asymptotic regime of the disturbed shock evolution is the use of Bessel functions, written in the form

$$ \Xi_s(s) = \sum_{\nu=0}^{\infty} D_{\nu} \frac{e^{-\nu \sinh^{-1} s}}{\sqrt{s^2 + 1}}, \quad (31) $$

in the Laplace variable $s$. The coefficients $D_\nu$, obtained with the aid of (29), are given by the following recurrence relation

$$ D_{2\nu} = D_{2\nu-4} \frac{1 - \sigma_b}{1 + \sigma_b} + D_{2\nu-2} \frac{2\sigma_b - 4\sigma_c}{1 + \sigma_b}, \quad (32) $$

which calls for the first values

$$ D_0 = 1, \quad D_2 = \frac{2}{1 + \sigma_b} + 2\frac{\sigma_b - 2\sigma_c}{1 + \sigma_b}, \quad \text{and} \quad D_4 = \frac{2(1 - \sigma_b)}{1 + \sigma_b} + 2\frac{\sigma_b - 2\sigma_c}{1 + \sigma_b} D_2 \quad (33) $$

to be initiated, with the odd values for $D_{2\nu+1}$ being found to be null.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Asymptotic temporal evolution of the shock ripple amplitude $\xi_s(\tau)$ for $\gamma = 1.4$, $M_1 = 10$, $Q = 15$, and different values of the energy sensitivity parameter $\epsilon$. Dashed lines correspond to asymptotic regression trends.</figcaption></figure>

The temporal evolution, given by

$$ \xi_s(\tau) = \sum_{\nu=0}^{\infty} D_{\nu} J_{\nu}(r = \tau \sqrt{1 - M_2^2}), \qquad (34) $$

is computed in Fig.4 for $M_1 = 10$, $\gamma = 1.4$, and $Q = 15$, and some distinguished values of $\epsilon$. Computations correspond to the triangle markers placed along the gray-dashed line in next Fig.5 and the color-code employed to identify the regimes is the same as in Fig.5. The energy sensitivity parameter ranges from -10 to 10, passing through all distinguished regimes. The points correspond to the instant value of the shock ripple amplitude (chosen at arbitrary maxima for the oscillating cases) and the lines refer to the asymptotic trends. For stable cases, including the constant-oscillating case, it is observed how the temporal evolution of the shock front approaches the long-time asymptotic trends in 10 times the characteristic time $k/a_2$. The unstable case, on the other hand, rapidly follows the asymptotic exponential trend, as seen in the log-scale inset.
When the shock ripple grows exponentially, i.e., for $\sigma_b < -1$, the long-time expression for the shock-ripple amplitude is

$$ \xi_s(\tau \gg 1) = \xi_s^0 \exp \left( \sqrt{1 - M_2^2} \sigma \tau \right), \qquad (35) $$

where the pre-exponential factor and the asymptotic growth rate are determined by

$$ \xi_s^0 = \frac{m-1}{4m} + \frac{\sigma_c(\sigma_c - \sigma_b)}{m\sqrt{2[1 + 2\sigma_c(\sigma_c - \sigma_b) + m]}}, \qquad (36) $$

and

$$ \sigma = \sqrt{\frac{1 - 2\sigma_c\sigma_b + m}{2(\sigma_b^2 - 1)}} \qquad (37) $$

respectively, with $m = \sqrt{1 + 4\sigma_c(\sigma_c - \sigma_b)}$. For the case considered in Fig.4 ($\gamma = 1.4$, $M_1 = 10$, $Q = 15$) the growth rate renders $\sigma\sqrt{1 - M_2^2} = 8.744$, which agrees fairly well with the regression line computed in Fig.4.

On the other hand, when the shock evolves with constant-amplitude oscillations in the long-time regime ($-1 < \sigma_b < \sigma_c$), the asymptotic behaviour is directly given by the function

$$ \xi_s(\tau \gg 1) = \xi_s^\infty \cos \left( \sqrt{1 - M_2^2} \omega \tau \right), \qquad (38) $$

where the amplitude and the frequency are given by

$$ \xi_s^\infty = \frac{m-1}{2m} + \frac{\sigma_c}{m} \sqrt{\frac{m-1-2\sigma_b(\sigma_c-\sigma_b)}{2(1-\sigma_b^2)}} \qquad (39) $$

and

$$ \omega = \sqrt{\frac{2\sigma_b\sigma_c - 1}{2(1 - \sigma_b^2)} \left( \frac{m}{2\sigma_b\sigma_c - 1} - 1 \right)}, \qquad (40) $$

respectively. The value of the asymptotic amplitude is $\xi_s^\infty = 0.207$ for the case considered in Fig.4, which presents a perfect agreement with the long-time amplitudes given by the equation describing the transient evolution (34). From (39) and (40) it is seen the transition towards the unstable case previously shown, as they become singular for $\sigma_b = -1$.

## 2. Stability limits

The stability limits can be also derived by the corresponding normal-mode analysis, as done in pioneering references [20] and [21]. In particular, the unstable, sonic radiating and highly-damped regimes are delimited by the condition $\Gamma_s = \varphi_{uns}$, $\varphi_{rad}$ or $\varphi_{dam}$, where
$$ \varphi_{\text{uns}} = -(1 + 2\mathcal{M}_2^2), \qquad (41) $$

$$ \varphi_{\text{rad}} = \frac{1 - \mathcal{M}_2^2(\mathcal{R}_s + 1)}{1 + \mathcal{M}_2^2(\mathcal{R}_s - 1)}, \qquad (42) $$

$$ \varphi_{\text{dam}} = \frac{\mathcal{R}_s \mathcal{M}_2^2 - (1 - \mathcal{M}_2^2)^{3/2} \sqrt{1 - \mathcal{R}_s^{-1}}}{\mathcal{R}_s \mathcal{M}_2^2 + 1 - \mathcal{M}_2^2}, \qquad (43) $$

respectively. It can be seen that $\varphi_{\text{uns}}$, $\varphi_{\text{rad}}$ and $\varphi_{\text{dam}}$ depend on the global jump conditions, that is, they do not depend on the local slope of the Rankine Hugoniot curve. For the sake of clarity, a sketch is depicted in Fig.5 along with the iso-curves for $\varepsilon$ shown previously in Fig.2. On the right side there is a bar indicating the stability limits in terms of the traditional nomenclature ($\varphi$) and that naturally given by the Laplace transform ($\sigma$).

<figure><img src="image_6.png" /><figcaption>FIG. 5. Left: stability regions as a function of the energy release $q$ and DK parameter $\Gamma_s$. Iso-curves of $\varepsilon$ as shown in Fig.2. Triangles correspond to the computations shown in the appendix. Right: sketch of the distinguished regimes for both $\Gamma_s$ and $\sigma_b$ parameters.</figcaption></figure>

It is readily seen that, for the unstable region to be reached (zone in red), the value of $\Gamma_s$ must be negative (reverse RH curve), as widely known result. In terms of the energy sensitivity, for this to happen the value of $\varepsilon$ must be positive and above $(\mathcal{M}_1 - \mathcal{M}_1^{-3})/2$. That is, $\varepsilon$ must be of the same order as $\mathcal{M}_1$. On the other hand, the radiating condition (zone in orange), conditions for the energy sensitivity are much less restrictive, with small
values being sufficient to enter into the non-decaying region. Regular shocks, here defined by $Q = \varepsilon = 0$ lie in the green zone. For this case, as the computation in Fig.5 is made for $M_1 = 10$, a strong-shock case, the corresponding value of the DK-parameter $\Gamma_s \gtrsim \varphi_{\text{rad}}$. It is for negative values of $q$ or $\varepsilon$ when the shock may enter into the highly-damped regime (zone in blue). It is observed that curves $\Gamma_s = \varphi_{\text{rad}}$ and $\Gamma_s = \varphi_{\text{dam}}$ collapse towards unity in the CJ-regime, given at $Q = Q_{\text{max}}$. It highlights the fact the sonic condition downstream is a distinctive regime, where transverse acoustic waves do not have effective time move along the shock surface.

Once defined the distinctive limits, the associated critical values of $\varepsilon$ can be given, provided that $\Gamma_s$ is defined in (16). For example, the unstable threshold $\Gamma_s = \varphi_{\text{uns}}$ gives

$$ \varepsilon_{\text{uns}} = - \frac{\varphi_{\text{uns}} \frac{\mathcal{R}_s^2}{\gamma \mathcal{M}_1^2} \frac{\partial \mathcal{P}_s}{\partial \mathcal{M}_1} - \frac{\partial \mathcal{R}_s}{\partial \mathcal{M}_1}}{\varphi_{\text{uns}} \frac{\mathcal{R}_s^2}{\gamma \mathcal{M}_1^2} \frac{\partial \mathcal{P}_s}{\partial \mathcal{Q}} - \frac{\partial \mathcal{R}_s}{\partial \mathcal{Q}}} \quad (44) $$

that admits a simple expression in the strong-shock limit, namely

$$ \varepsilon_{\text{uns}} (M_1 \gg M_{\text{cj}}) = \frac{2 + \sqrt{2\gamma(\gamma - 1)}}{\gamma + 1} M_1 + O(M_1^{-1}). \quad (45) $$

It dictates that $\varepsilon_{\text{uns}}$ must be of the order of $M_1$ for the strong shock to be fully unstable. Similarly, when value of the energy variations are dominant, either for $Q = Q_{\text{min}}$ or $Q = Q_{\text{max}}$, the value of $\varphi_{\text{uns}}$ turns $-1$ and $-3$, respectively, which renders

$$ \varepsilon_{\text{uns}} (Q_{\text{min}}) = \frac{(\gamma^2 - 2\gamma - 3 - 4Q_{\text{min}}) \sqrt{\gamma + 1}}{2\sqrt{2(\gamma - 1)(-2Q_{\text{min}} - \gamma - 1)}} = (\gamma + 1) \frac{\mathcal{M}_1^2 + 1}{2\mathcal{M}_1} \quad (46) $$

and

$$ \varepsilon_{\text{uns}} (Q_{\text{max}}) = 2\sqrt{Q_{\text{max}}} \left[ 1 + \sqrt{Q_{\text{max}}} \left( 1 - \sqrt{Q_{\text{max}} + 1} \right) \right] = \frac{\mathcal{M}_1^4 - 1}{2\mathcal{M}_1^3}, \quad (47) $$

which also places the critical energy sensitivity of the order of $M_1$.

The function $\varepsilon_{\text{uns}}$ is computed in Fig.6 as a function of $Q$ for $\gamma = 1.4$ (left) and $\gamma = 1.01$ (right) and for different values of the shock Mach number $M_1$. The corresponding limits for $\varepsilon_{\text{uns}}(Q_{\text{min}})$ and $\varepsilon_{\text{uns}}(Q_{\text{max}})$ are also plotted (dashed lines). While the latter does not change with the adiabatic index, the former is highly affected, as evidenced from (46), which shows a pole for $Q_{\text{min}} = -(\gamma + 1)/2$, but this condition is never achieved since $Q_{\text{min}} < -(\gamma + 1)^2/4$ for $M_1 > 1$. As predicted through the asymptotic values shown in (46) and (47), the fully unstable regime calls for very high values of $\varepsilon$.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Function $\epsilon_{uns}$ as a function of $Q$ for $\gamma = 1.4$ (left) and $\gamma = 1.01$ (right) and for different values of the shock Mach number $M_1$</figcaption></figure>

Similarly, the critical value of energy sensitivity for the radiation condition is

$$ \epsilon_{\text{rad}} = - \frac{\varphi_{\text{rad}} \frac{\mathcal{R}_s^2}{\gamma \mathcal{M}_1^2} \frac{\partial \mathcal{P}_s}{\partial \mathcal{M}_1} - \frac{\partial \mathcal{R}_s}{\partial \mathcal{M}_1}}{\varphi_{\text{rad}} \frac{\mathcal{R}_s^2}{\gamma \mathcal{M}_1^2} \frac{\partial \mathcal{P}_s}{\partial Q} - \frac{\partial \mathcal{R}_s}{\partial Q}} = \frac{\mathcal{M}_1^2 - 1 + \sqrt{(\mathcal{M}_1^2 - 1)^2 + 4Q\mathcal{M}_1^2}}{2\mathcal{M}_1^3}, \quad (48) $$

which does not depend on the adiabatic index. It is also instructive to evaluate $\epsilon_{rad}$ in the strong-shock limit, i.e.,

$$ \epsilon_{\text{rad}} (M_1 \gg M_{\text{cj}}) = \frac{1}{M_1} - \frac{1+Q}{M_1^3} + O(M_1^{-4}), \quad (49) $$

where it is found that radiation condition is rather feasible for strong shocks, it requires very small values of $\epsilon$ as it scales with $M_1^{-1}$ in this limit. This trend does not change in the limits $Q = Q_{max}$ and $Q = Q_{min}$, since

$$ \epsilon_{\text{rad}} (Q_{\text{max}}) = \frac{\sqrt{Q_{\text{max}}}}{(\sqrt{Q_{\text{max}}} + \sqrt{1+Q_{\text{max}}})^2} = \frac{\mathcal{M}_1^2 - 1}{2\mathcal{M}_1^3}, \quad (50) $$

$$ \epsilon_{\text{rad}} (Q_{\text{min}}) = \frac{(\gamma + 1)\sqrt{\gamma^2 - 1}}{2\sqrt{2}\sqrt{-2Q_{\text{min}} - \gamma - 1}} = \frac{\gamma + 1}{2M_1}. \quad (51) $$

To complete the asymptotic analysis, the critical value of energy sensitivity for the highly-
<figure><img src="image_8.png" /><figcaption>FIG. 7. Critical values of energy change sensitivity, $\varepsilon_{\text{rad}}$ and $\varepsilon_{\text{dam}}$, as a function of $Q$ for different propagation Mach numbers $M_1$ and for $\gamma = 1.4$ (left) and $\gamma = 1.01$ (right).</figcaption></figure>

damped condition is found to be

$$ \varepsilon_{\text{dam}} = - \frac{\varphi_{\text{dam}} \frac{\mathcal{R}_s^2}{\gamma \mathcal{M}_1^2} \frac{\partial \mathcal{P}_s}{\partial \mathcal{M}_1} - \frac{\partial \mathcal{R}_s}{\partial \mathcal{M}_1}}{\varphi_{\text{dam}} \frac{\mathcal{R}_s^2}{\gamma \mathcal{M}_1^2} \frac{\partial \mathcal{P}_s}{\partial Q} - \frac{\partial \mathcal{R}_s}{\partial Q}}, \quad (52) $$

which does not admit a simple expression, except for the strong-limit

$$ \varepsilon_{\text{dam}} (M_1 \gg M_{\text{cj}}) = 2 \frac{\sqrt{\gamma} - 1}{\sqrt{\gamma} + 1} M_1 + O(M_1^{-1}) \quad (53) $$

and the limiting cases for $Q = Q_{\text{max}}$

$$ \varepsilon_{\text{dam}} (Q_{\text{max}}) = -2 \sqrt{Q_{\text{max}}} \left( \sqrt{1 + Q_{\text{max}}} - \sqrt{Q_{\text{max}}} \right) = - \frac{(M_1^2 - 1)^2}{2M_1^3} \quad (54) $$

and $Q = Q_{\text{min}}$, the latter being the same as the radiating condition $\varepsilon_{\text{dam}} (Q_{\text{min}}) = \varepsilon_{\text{rad}} (Q_{\text{min}})$.

The values of $\varepsilon_{\text{rad}}$ and $\varepsilon_{\text{dam}}$ are computed in Fig.7 as a function of the energy release $Q$ for different propagation Mach numbers $M_1$ and for $\gamma = 1.4$ (left) and $\gamma = 1.01$ (right), where the longitudinal domain covers $Q_{\text{min}} < Q < Q_{\text{max}}$. In solid lines are displayed the curves associated to the DK limit, $\varepsilon_{\text{rad}}$, while the values of $\varepsilon_{\text{dam}}$ are shown in dashed lines. It is readily seen that the region in between the solid and dashed lines refer to regular decay
of the shock oscillations. This region collapses for $Q = Q_{\text{min}}$, since $\varepsilon_{\text{rad}} = \varepsilon_{\text{dam}}$ in this limit. The solid-line curves associated to $\Gamma_s = \varphi_{\text{rad}}$ do not depend on the adiabatic index, as seen in (48), so the only difference between left and right panels is the stretching of the longitudinal domain given by the dependence of $Q_{\text{min}}$ with $\gamma$, since $Q_{\text{max}}$ only depends on the shock strength. For each value of $M_1$, the region above the limit $\Gamma_s = \varphi_{\text{rad}}$ corresponds to the permanent-oscillations regime. The dashed-line curves associated to $\Gamma_s = \varphi_{\text{dam}}$ do depend on the adiabatic index and the region below them for each Mach number is assigned to highly-damped oscillations. The end point on the right of each curve corresponds to the associated CJ value, $Q = Q_{\text{max}} = \sqrt{M_1^2 - 1} + \sqrt{M_1^2}$.

Qualitatively, it is observed that the preferred region to get permanent oscillations in the first quadrant, for $Q > 0$ and $\varepsilon > 0$. Contrarily, the most likely configuration to get oscillations that decay very fast is when $Q < 0$ and $\varepsilon < 0$, although these criteria are not necessarily true, as depicted in Fig.7. As for the fully unstable conditions shown in Fig.6, they are restricted to very high positive energy sensitivities with the shock strength.

### III. APPLICATION TO HIGH-ENERGY-DENSITY SHOCKS

The analysis presented in the preceding section is applicable under the following assumptions: the non-adiabatic transformations of the fluid particles occur in a very thin layer right behind the shock for the whole supersonic front to be considered a fluid discontinuity. This transformation is allowed to depend on the local properties and, therefore, on the shock Mach number. Besides, pressure, density and temperature are related through the perfect gas equation of state in the perturbation variables in the equilibrium states, but *ad hoc* non-ideal thermodynamic changes can be incorporated across the RH equations, on condition that they become negligible in the equilibrium state downstream.

One example of this type is a gaseous detonation. In this case, although heat release rate is mainly affected by local properties behind the precursor shock wave, the overall heat released also depends on the thermodynamic conditions once equilibrium conditions are reached. Sufficiently strong changes in pressure or temperature could promote the formation of different radicals, in which case the global thermodynamic equilibrium involving a combination of different species could affect the overall amount of heat released. When the detonation is driven by reactions that depend on local thermodynamic functions, the
heat released depends on the shock strength. Certainly, any particular case would call for a
specific model to compute the variations in the global heat released with the shock intensity.
However, because the pressure and temperature scale with the shock strength squared for
sufficiently strong shocks, such as those expected in detonations, variations of the global heat
released with the thermodynamic variables can be collected as term $Q = Q_0(1 + \delta Q M_1^2)$,
where $\delta q$ accounts for the dimensionless variation of heat release with two independent
thermodynamic functions.

The aim of this section is to find and analyze examples where the above-mentioned
hypothesis are fulfilled with sufficient accuracy. The goal is to characterize flow transfor-
mations that may appear in strong-shock conditions in terms of the governing parameters:
upstream adiabatic index $\gamma$, shock Mach number $M_1$, dimensionless energy change $q$ (that
accounts for changes in the gas molecular structure) and energy sensitivity with the shock
strength $\epsilon$. They include: excitation of the vibrational degrees of freedom, molecular dis-
sociation, ionization, electro-magnetic radiation. A discussion of nuclear transformations is
also included. For the sake of clearness, the dominant effects that distinguish any of these
phenomena are studied independently, although some of them may likely coexist in real con-
ditions. For the sake of generality and in order to place the focus on the energy variations,
a perfect-gas equation of state is assumed, admitting this may not be accurate for some
high-energy-density scenarios [69].

A. Energy function for ionization, dissociation, radiation, and nuclear reactions

1. *Ionization of monotomic gases*

When a shock is sufficiently strong to ionize the gas, assumed monatomic, part of the
electron shell is removed from the atom electronic cloud. This is an endothermic effect that
requires some energy to be accomplished. Therefore, due to the non-monotonic behavior of
the shock adiabatic, it sensibly to wonder if ionizing shocks could generate SAE. In fact,
Glass and Liu detected density and electrons' concentration perturbations far behind strong
shocks traveling in inert argon [70] and krypton [71]. In their shock-tube experiments, they
observed that such instabilities appeared if the shocks were sufficiently strong. Nevertheless,
these observed instabilities are not predicted using the classical DK criteria. It was not until
the 90's, when Mond and Rutkevich theoretically shed some light on this phenomenon. In their works [6, 7], they demonstrated that ionization by electron-atom collisions and thermal non-equilibrium between the electrons and heavy particles are indispensable to get SAE. Like in this work, they reduce the gas-dynamical shock together with the relaxation zone to an effective surface of discontinuity, called ionizing shock. Certainly, such approximation, robust for gas-dynamical shocks, is not appropriate when the characteristic length of the perturbations is of the same order (or smaller) than the ionizing shock thickness. For this particular case, the width of the relaxation zone, where the quantities change continuously due the inelastic processes (excitation of the internal molecular energy levels, dissociation and ionization of the molecules or atoms, etc), might not be negligible. Therefore, in order to consider the ionizing shock as a discontinuous zero-length front, it must be assumed that the wavelengths of the emitted acoustic waves are much larger than the width of the relaxation zone. Therefore, the acoustic perturbations should be small enough such that $kl \ll 1$, where k is the distinctive perturbation wavelength and $\ell$ is the characteristic length for the variation of zero-order flow variables within the relaxation zone.

For this case $Q$ is restricted to the energy per unit mass subtracted to the fluid particles because of single ionization in monatomic gases. The excitation of electronic levels are neglected since they barely affect the bulk-flow properties nor the stability limits [7]. The parameter $I_i^1$ defines the minimum amount of energy required to remove the most loosely bound electron of an isolated neutral gaseous atom. As it depends on the number of electron removed, a general form to express the n$^{th}$ ionization is $X^{(n-1)+} + I_i^n \rightarrow X^{n+} + e^-$ where $X^{(n-1)}$ is any atom capable of ionization, $X^{n+}$ is that atom with the n$^{th}$ electron removed, and $e^-$ is the removed electron. As detailed in appendix A, the energy employed to ionize the gas can be written in terms of the following dimensionless variables,

$$Q = -\alpha_i \frac{\gamma^2 - 1}{2\gamma} \phi_i. \qquad (55)$$

which is found to be proportional to the ionization degree $\alpha_i$ and the dimensionless energy parameter $\phi_i = I_i^1/(k_B T_1)$, where $k_B$ is the Boltzmann constant. The problem formulation is closed, for a given Mach number $M_1$ and adiabatic index $\gamma$, with the aid of (12), (8), (10) (for $\alpha = \alpha_i$), (55), and the corresponding Saha equation for the dissociation degree (A2). They form a system of five equations for $R_s, P_s, T_s, Q$ and $\alpha_i$. The equation of state depends on the ionization degree as a result of the contribution of the partial electron pressure to the
total pressure of the ionized plasma [72], on condition that thermal equilibrium is reached.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Energy removed by ionization $Q$, associated energy sensitivity $\epsilon$ and degree of ionization $\alpha_i$ as a function of the Mach number $M_1$ for argon. Triangle: value of $\alpha_i$ according to Ref. [70]</figcaption></figure>

To establish a connection with the shock stability theory presented before, it is necessary to evaluate the value of $q(M_1)$ along with the value of the energy change sensitivity $\epsilon = (dQ)/(dM_1)$. Figure 8 shows the energy removed by ionization $Q$, associated energy sensitivity parameter $\epsilon$, and degree of ionization $\alpha_i$ as a function of the Mach number $M_1$ for argon. As expected, the function $\alpha_i$ approaches unity for high Mach numbers, and the parameter $Q$ follows the same trend but multiplied by the factor $\sim 325$ given by (55). Complete ionization is found to happen for $M_1 \sim 37$ as shown in Fig. 12 and highlighted in the inset. However, this is found to occur for very strong shocks, for which the shocked gas may undergo subsequent ionization and photoionization processes [73]. Such a multiple ionization might result in the appearance of a sequence of maxima in the function $R_s(P_s)$ shown in Fig. 12. The parameter $\epsilon$ reaches a nearly plateau region at intermediate shock intensities after which it rapidly drops as the $Q$ becomes constant. For this case, ionization and thermal equilibrium are established at each time during the perturbation. However, thermal equilibrium between electrons and the heavy particles may not be achieved during one period of the perturbation when the oscillation frequency is sufficiently high.
## 2. Dissociation of diatomic gases

In standard conditions, diatomic gases are characterized by $\gamma = 7/5$ as two rotational plus three translational degrees of freedom enter into play. For very low temperatures the gas behaves as it were monotonic as molecules stop rotating, thereby yielding $\gamma = 5/3$. By way of contrast, if temperature increases (~3000 K for H₂), the two-atom molecule behaves as a pure oscillator and the idealized vibrational mode should be included by adding two degrees of freedom. They correspond to potential and kinetic energy components. In total there exist 7 degrees of freedom that render $\gamma = 9/7$ for high-temperature conditions. At some point, inter-atom vibrations are sufficiently strong to split apart the two atoms of the molecule, a process that is called dissociation. If all molecules have been dissociated the gas becomes monatomic and $\gamma = 5/3$ is recovered.

For this case the dimensionless energy per unit of mass subtracted to the fluid particles

$$Q = -\frac{12}{35} \left[ \alpha_d \left( \frac{T_s}{2} + \phi_d \right) + (1 - \alpha_d) \frac{\phi_v}{e^{\phi_v/T_s} - 1} \right] \quad (56)$$

is not simply proportional to the degree of dissociation $\alpha_d$ as postshock properties enter into play through the factor $T_s$. The dimensionless factors $\phi_d$ and $\phi_v$ refer to the corresponding dissociation and vibrational potentials, as shown in appendix B. In obtaining (56) it has been made use of $\gamma_1 = 7/5$ (upstream), $\gamma_2 = 7/5$ (downstream non-dissociated gas) and $\gamma_3 = 5/3$ (downstream dissociated gas). The assumption of $\gamma_2 = 7/5$ stems from the fact that vibrational degrees of freedom are neglected. The problem formulation is closed, for a given Mach number and the adiabatic indices, with the aid of (12), (8), (10) (for $\alpha = \alpha_d$), (56), (B2), which comprise a system of five equations for $R_s$, $P_s$, $T_s$, $Q$ and $\alpha_d$.

Figure 9 shows the energy removed by dissociation $Q$, associated energy sensitivity $\epsilon$ and degree of dissociation $\alpha_d$ as a function of the Mach number $M_1$ for N₂ [10]. Unlike previous ionization process, the endothermic process does not cease when gas runs out of diatomic molecules $\alpha_d = 1$ since there exists a term proportional to $T_s$ in (56) as a result of the change in the adiabatic index. This effect would be modified if the change in the adiabatic index as a result of potential vibrations is accounted, i.e., $\gamma_2 = 9/7$ or, more accurately, a temperature-dependent function is employed to describe $\gamma_2(T_s)$.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Energy removed by dissociation $Q$, associated energy sensitivity $\epsilon$ and degree of dissociation $\alpha_d$ as a function of the Mach number $M_1$ for $N_2$ [10].</figcaption></figure>

### 3. Electromagnetic radiation

Several astrophysical phenomena involve radiative shocks (RS's) with various physical characteristics regarding their spatial structure or their dynamics [43, 58–64]. As an example, they occur in supernovae, either in core-collapse or thermonuclear [45]. For gravitational or core-collapse supernova explosion (CCSNe), the RS is nearly adiabatic at the early time of the propagation of the shock after bounce [74] and once it is recognized that the adiabatic coefficient $\gamma$ of a fluid with matter and photons is $\gamma = 4/3$, the material compression, $R_s \equiv \rho_2/\rho_1$, at the shock front is close to the value $R_s = (\gamma + 1)/(\gamma - 1) = 7$ [43, 74]. In contrast, thermonuclear SN experience a thermonuclear detonation (see next section) with energy release behind and at the shock front making the adiabatic approximation not anymore valid for the modeling of the propagation of the RS. Nevertheless, in both cases radiation plays a significant role and the energy and pressure of the photons should be taken into account in the momentum and energy balances.

RS's are also of great importance in the evolution of the interstellar medium. In that case, the radiative energy (resp. radiation pressure) is however negligible compared to the thermal energy (resp. to the thermal pressure). Indeed, since the ISM is mostly optically
thin, photons leave the downstream flow but, at the same time, they produce a cooling in the bulk of the fluid. As a consequence, a compression above $R_s \sim 100$ can be achieved in the shocked material although the shock front is adiabatic with a compression $R_s \sim 4$ (the ISM has $\gamma = 5/3$) [75]. This cooling can be described by an energy loss per unit mass, $Q$, at the shock discontinuity. Barenblatt and his colleagues [76, 77] have accounted such an effect in a phenomenological way by choosing a linear dependence $q(T_2) \propto T_2$ and by taking an effective adiabatic coefficient, $\gamma_{\text{eff}}$, at the discontinuity smaller than $\gamma$. With such a free parameter, the authors evidence analytically the high densification of the material at the shock front due to cooling and they show that for $\gamma_{\text{eff}} \to 1$, one gets $R_s \gg 1$. This linear approximation for $Q$ is however not relevant for an actual RS where the losses at the density jump should be described by the radiation flux $F_r = -f \sigma T_2^4$, where $\sigma$ is the Stefan-Boltzmann constant and where $f$ is a positive numerical factor (the minus sign accounts for the energy leaving the downstream flow). Furthermore, for the sake of simplicity the radiation energy and the radiation pressure will be neglected compared to the thermal energy and the thermal pressure of the gas, respectively.

For an optically thick material, the photon distribution obeys the black body law and $f = 1$, and for an optically thin fluid, one has is of the order of $L/\lambda \ll 1$ where $L$ is the longitudinal extension of the downstream flow and where $\lambda$ is the photon mean free path [59, 63]. As derived in the appendix C, the dimensionless energy loss becomes

$$Q = -\phi \frac{\gamma^2 - 1}{4M_1} T_s^4 \qquad (57)$$

where $\phi = 2f\sigma T_1^4/(\rho_1 a_1^3)$. In the strong-shock limit, the radiative flux scale as $F_r(M_1 \gg 1) = -\rho_1 a_1^3 M_1^3/2$. Replacing the Mach number by its expression gives $F_r(M_1) = \rho_1 u_1^3/2$, which is equivalent as saying that all kinetic energy of the incoming material is radiated away from the downstream flow. The high value of the temperature behind the shock front then results from both the strong compression of the material and the thermal energy of the upstream gas. On the other side, the radiative dimensionless energy $Q(M_1 \gg 1) = -(\gamma^2 - 1)M_1^2/4$, shows an interesting behavior, as it does not depend on $\phi$ and it equals $Q_{\text{min}}$ in (14) in the limit $M_1 \gg 1$. That is, the energy radiated is the maximum energy the fluid particles can lose. In this condition, pressure remains finite (by Rayleigh line) $p_2 = p_1(1 + \gamma M_1^2)$ but velocity remains stagnant in the shock reference frame as a result of the density unbounded rise. Consequently, the energy sensitivity is $\varepsilon(M_1 \gg 1) = -(\gamma^2 - 1)M_1/2$, which is always
negative, thereby anticipating a decaying regime for the shock oscillations.

<figure><img src="image_11.png" /><figcaption>FIG. 10. Energy loss by radiation $Q$ and the associated energy sensitivity $\epsilon$ as a function of the Mach number $M_1$ for $\gamma = 5/3$.</figcaption></figure>

Figure 10 shows the energy loss by radiation $Q$ and the associated energy sensitivity $\epsilon$. It is clearly identified to distinguished regimes. The relatively low Mach number regime for which energy release is very small and consequently, post-shock flow variables remains effectively unaltered, and the latter associated to the maximum compression ratio determined by $q(M_1 \gg 1) = Q_{min}$. The transition occurs for shock intensities of the order of $M_1 \sim \phi^{-1/5}$.

### 4. Thermonuclear transformations

At the end of their lives, massive stars develop central cores made of iron nuclei and generate electrons. The iron core becomes unstable when it accumulates more than ~1.4 solar masses, a limit set by the maximum pressure that the degenerate electrons can provide. A rapid collapse to a proto-neutron star ensues, launching a strong shock into the stellar envelope. As it ventures outwards, the shock has to break up iron nuclei, exhausting its energy and quickly turning into a stagnant accretion shock. In order to produce a supernova, known as core-collapse supernova (CCSNe), the shock has to revive in less than a second and then expel the stellar envelope. Despite the decades of effort, the exact details of how to
revive the shock is not fully understood. Heating by neutrinos emitted by the proto-neutron star and the interaction of the shock with multi-dimensional hydrodynamic phenomena, such as turbulence and standing accretion shock instability, are expected to play a major role in its revival (see, e.g., [55–57] for recent reviews). As a first approach, the breaking of the heavy nuclei can be incorporated into the non-adiabatic shock model akin to the molecular dissociation presented before. For this reason, this effect is not computed in this section, although further details can be found in appendix D.

Unlike CCSNe, thermonuclear supernovae develop as a the nuclear ignition temperature is achieved inside white dwarf stars composed primarily of C and O nuclei and detached degenerate electrons [45, 47]. At high temperatures, thermonuclear sub-products may endothermically dissociate with the corresponding energy subtractions. Therefore, the net energy balance across the thermonuclear SNe must account for both the coupled endothermic and exothermic processes coexisting behind the shock. Besides, in realistic situations the nuclear fuel does not undergo full burning at the shock, an effect that can be modelled by means of an ad-hoc contribution ~exp(-φ<sub>eff</sub>/T<sub>s</sub>), which converges to unity with temperature as 1 - φ<sub>eff</sub>/T<sub>s</sub>. This is roughly similar to solving the associated Saha equation. Low values of φ<sub>eff</sub> are associated with efficient burning rates even at low temperature, while the contrary applies for high values of φ<sub>eff</sub>, which demands very high temperatures to get order-of-unity values of the fuel mass fraction consumption. A rough estimate for this factor in the strong-shock limit is φ<sub>eff</sub> ~ M<sup>2</sup><sub>1</sub>, which leads to partial burning in the whole domain but takes into account the temperature-dependent character of the process. This is used (see appendix D for details) to write the net dimensionless energy release as

$$Q = \frac{\gamma^2 - 1}{2\gamma} \exp\left(-\frac{\phi_{eff}}{T_s}\right) (\phi_{tn} - X_{\alpha}\phi_{nd}), \quad (58)$$

where φ<sub>tn</sub> is the dimensionless contribution of the exothermic thermonuclear reactions and φ<sub>nd</sub> measures the energy subtracted by dissociation processes, on condition that X<sub>α</sub> ≪ 1. There exists a positive correlation between the amount of heat release and shock intensity, on condition that X<sub>α</sub> is small. As found in gaseous detonations [38], this may lead to constant-amplitude oscillations of the perturbed shock.
## B. Stability limits for ionization, dissociation, radiation, and nuclear reactions

Figure 11 shows the DK parameters $\Gamma_s$, $\varphi_{\text{rad}}$ and $\varphi_{\text{dam}}$ as a function of the shock Mach number for adiabatic shocks, gaseous detonations as in [38], ionization shocks, dissociation shocks, radiative shocks and thermonuclear detonations. The parameters chosen are specified in previous sections. For the six cases computed, it can be seen that fully unstable regimes are never achieved.

For adiabatic shocks, Fig.11 (a), the function $\Gamma_s$ is always in the green region, except for $M_1 \gg 1$ where it approaches the function $\varphi_{\text{rad}}$. This is consistent with the regular oscillation decay $t^{-3/2}$ ($t^{-1/2}$ in the strong-shock limit) of the shock amplitude oscillations [2]. Panel (b) shows the case of gaseous detonations, whose heat release is modelled with $Q = Q_0(1 + \delta Q M_1^2)$ for as small positive value of the correlation parameter $\delta Q$. It shows the transition to the spontaneous acoustic emission, as predicted in [38]. Negative values of $\delta q$ will place the curve $\Gamma_s$ above $\varphi_{\text{dam}}$, in the highly-damped regime. Results for $\delta Q = 0$ [32, 35, 36]

For ionization shocks, the parameter $\Gamma_s$ is plotted in Fig.11 (c) as a function of the Mach number $M_1$ for argon at $T_1 = 300$ K and $p_1 = 5$ Torr., along with the functions $\varphi_{\text{rad}}$ and $\varphi_{\text{dam}}$. It can be observed that the shock is always stable for any value of shock strength, in agreement with [7] when thermal equilibrium is imposed. Moreover, the shock enters into the highly-damped oscillations regime. In their works [6, 7], it was shown that ionization by electron-atom collisions and thermal nonequilibrium between the electrons and heavy particles are indispensable to get SAE.

Dissociation shocks, according to the computations made for nitrogen in preshock standard conditions, exhibit a similar behaviour, as seen in Fig.11 (d). The shock decays in a regular way (green zone) except in a medium region ($10 \gtrsim M_1 \gtrsim 30$, for this particular case), where its oscillations decay faster ($\Gamma_s > \varphi_{\text{rad}}$). Although not shown, beyond $M_1 = 30$, the curve for $\Gamma_s$ asymptotically approaches the curve $\varphi_{\text{rad}}$. This result agrees with [10], where it was shown that solely dissociation process is not enough to produce spontaneous acoustic emission, and that ionization and non-thermal equilibrium are necessary for the shock to exhibit that behavior. The qualitative difference between Figs.11 (c) and (d) stems from two facts: on one side the post-shock Mach number (9) is affected by having $g_2 \neq \gamma$, on the other side ionization-induced energy variations tend to a constant when the shock Mach
<figure><img src="image_12.png" /><figcaption>FIG. 11. DK parameter $\Gamma_s$ and limits $\phi_{rad}$ and $\phi_{dam}$ as a function of the Mach number $M_1$ for adiabatic shocks (a), gaseous detonations as in [38] (b), ionization shocks in argon (c), dissociation shocks in molecular nitrogen (d), radiative shocks for $\gamma = 5/3$ (e), and thermonuclear shocks for $\gamma = 4/3$ (f).</figcaption></figure>
number increases, while the energy variation grows unbounded as a result of the dissociation
process for the change in the adiabatic index. That is, $\epsilon = 0$ and $\epsilon < 0$ for ionization and
dissociation shocks, respectively, when $M_1 \gg 1$. In connection with Fig.5, which was com-
puted for $\gamma = 7/5$, a strong shock in thermal equilibrium will be placed in the green-or-blue
region placed on the left side ($Q < 0$) and below the axis $\Gamma_s = 0$.

Radiation effects are displayed in Fig.11 (e) for different radiation constants, $\phi = 10^{-1}$ and $10^{-5}$. None of these cases, and none in any other within the range $0 \le \phi \le 1$, will render spontaneous acoustic emission. It is rather the opposite, the endothermic effect $q < 0$ and its correlation with the shock intensity $\epsilon < 0$ makes the shock oscillations decay faster than the adiabatic case towards its final planar form for sufficiently strong shocks. This is in part due to the isolated effect of the radiation losses, where, unlike previous cases, no change in the equation of state has been considered. Moreover, the effect of pressure radiation could also enter into play, but the model should be rewritten to account not only changes in the energy equation but also the momentum equation.

Thermonuclear exothermic shocks, displayed in Fig.11 (f), are found to render sponta-
neous acoustic emission when the factor $\phi_{eff} \sim 10$. This effect is similar to that found in
[38] for gaseous detonations, where the detonation may oscillate permanently when heat
release correlation with the shock intensity was sufficiently strong. There it was found a
positive feedback that heat release exert on the detonation dynamics, which is not unlike
the Rayleigh criterion for the thermo-acoustic instability. Nuclear dissociation alone, not
computed for the sake of conciseness, was previously found to lie onto regular and highly-
damped oscillating regime [52]. The formulation presented here supports these results.

IV. CONCLUSIONS

This study is an extension of the previous work [38] that considered the D'yakov-
Kontorovich instability in gaseous detonations. The present work includes endothermic
and exothermic effects and the general formulation does not restrict itself to any particular
form for the dependence of energy $Q$ with the shock Mach number $M_1$. The theory does
restrict to low frequency oscillations whose characteristic length is much larger than the
non-adiabatic shock thickness and the perfect gas equation of state. Besides, constant base
flow and isolated boundary condition are assumed, which reduce the limits of validity of
the results [15, 16]. The effect of supporting boundary conditions, along with the effect of
radiation pressure and energy, is left for a future work.

In terms of general properties, it is found that the fully unstable regime, that predicting an
exponential growth of the shock front perturbations, calls for a very high energy sensitivity
$\epsilon = (dq)/(dM_1)$. That is, the energy variations associated to the change in the shock
intensity must be of the order of the shock Mach number. The acoustic radiating condition,
or spontaneous acoustic emission, is less demanding in regard to the energy sensitivity.
Small-to-moderate values of $\epsilon$, of the order of $M_1^{-1}$, may suffice to put the shock in constant
oscillation, even in the flow undergoes an endothermic transformation across the shock.
The opposite limit associated to a strong damping is easy to achieve when, although not
restricted to, endothermic configurations. Analytical expressions for $\epsilon$ in the distinguished
limits are provided, along with computations of the transient evolution of the shock front.

The general theory is also applied to different contexts of interest in high-energy-density
physics, they are ionization, dissociation, radiative and thermonuclear shocks. The model
allows the inclusion of any phenomenon that may alter the energy balance across the shock,
as it is the change in the adiabatic index or the excitation of vibrational degrees of freedom
in the molecules. The global energy change $Q$ is then analyzed in terms of its sensitivity
$\epsilon$ to study the stability limits. It is observed that endothermic processes generally place
the shock into regular and highly-damped oscillations, on condition that there exits ther-
mal equilibrium downstream. It applies for ionization, dissociation, which are in agreement
with [6, 10], but also radiation phenomena. Exothermic effects are more likely to develop
an unstable behaviour when energy release is positively correlated with the shock inten-
sity, as found in thermonuclear detonations, which aligns with previous results on gaseous
detonations [38].

ACKNOWLEDGMENTS

Work produced with the support of a 2019 Leonardo Grant for Researchers and Cul-
tural Creators, BBVA Foundation. FCC has received support from MINECO (Grant No.
ENE2016-75703-R) JCCM (Grant No. SB-PLY/17/180501/000264). EA thanks Rodrigo
Fernandez for productive discussions. EA acknowledges support from grants BR05236454
and AP05135753 of the Ministry of Education and Science of Republic of Kazakhstan.
Appendix A: Derivation of the monatomic-gas ionization energy equation

This section will consider the first ionization only. The value of $I_i^1$ is tabulated and can be found for many substances. For example, $I_i^1 = 15.7596$ eV for the first ionization energy of argon. This effect can be included through the one-dimensional energy conservation equation, namely

$$ \frac{\gamma}{\gamma - 1} \frac{p_1}{\rho_1} + \frac{u_1^2}{2} = \frac{\gamma}{\gamma - 1} \frac{p_2}{\rho_2} + \frac{u_2^2}{2} + \alpha_i \frac{I_i^1 N_A}{W} \quad (A1) $$

provided that thermodynamic equilibrium is achieved downstream and that the adiabatic index for the monatomic gas $\gamma = 5/3$ remains invariant since electrons and ions also have three degrees of freedom. The parameters $W$ and $N_A$ are the molecular weight of the species and is the Avogadro number, respectively. The parameter $\alpha_i = n_e/(n_a + n_i)$ defines the degree of ionization, where the terms $n_e$, $n_a$ and $n_i$ are the concentrations of free electrons, neutral atoms and positive ions, respectively. Neutral conditions for single ionization are met when $n_e = n_i$. The influence of radiation, can be arguably ignored in the range of parameters of interest since the characteristic length of the radiation-due changes is much longer than the width of the relaxation zone.

Assuming that temperature is not high enough for the fluid particles to undergo photoionization, the process can be fully described by accounting for ionization resulting from ions, atoms and electron collisions. Then, the degree of ionization can be written, in terms of the thermodynamic parameters, using the Saha equation [6]

$$ \frac{\alpha_i^2}{1 - \alpha_i} = \frac{G_i}{n_i + n_a} \left( \frac{m_e k_B T_2}{2\pi \hbar^2} \right)^{3/2} \exp \left( -\frac{I_i^1}{k_B T_2} \right), \quad (A2) $$

where $k_B$ is the Boltzmann constant, $\hbar$ is the reduced Planck constant, and $m_e$ is the mass of the electron. The function $G_i$ depends on the partition functions of ions and atoms. Since the electronic excitation is neglected, a constant value could be assumed. In fact, we consider $G_i = g_{0e}g_{0i}/g_{0a}$, where $g_{0e} = 2$, $g_{0i}$ and $g_{0a}$ are the statistical weights of the ground states of electrons, single-charge ions and atoms, respectively.

In terms of dimensionless variables, the degree of ionization is readily obtained after simple manipulation of (A2), reading

$$ \alpha_i = 2\sqrt{\Psi_i} \left( \sqrt{\Psi_i + 1} - \sqrt{\Psi_i} \right), \quad (A3) $$
where the function

$$ \psi_i(\mathcal{T}_s, \mathcal{R}_s) = B_i \frac{\mathcal{T}_s^{3/2}}{\mathcal{R}_s} \exp\left(-\frac{\phi_i}{\mathcal{T}_s}\right), \qquad (A4) $$

includes the dependence with the temperature and density jumps across the shock, $\mathcal{T}_s$ and $\mathcal{R}_s$, respectively. The constants

$$ \phi_i = \frac{I_i^1}{k_B T_1} \quad \text{and} \quad B_i = \frac{G_i}{4} \left( \frac{m_e}{2\pi\hbar^2} \right)^{3/2} \frac{(k_B T_1)^{5/2}}{p_1} \qquad (A5) $$

are functions of the upstream flow conditions. Then, along with the parameters describing the inelastic processes, the Hugoniot curve shape and the stability criteria are functions of the upstream flow properties through the parameters $\phi_i$ and $B_i$, which render ~ 610 and ~ 214, respectively, for argon at $T_1 = 300$ K and $p_1 = 5$ Torr [6]. The latter somewhat computes how effective ionization is by inter-atom and atom-electron collisions, while the former defines the incipient temperature condition for ionization to become non-negligible. The implicit algebraic system of equations is closed, provided that $\mathcal{R}_s$, $\mathcal{T}_s$ are defined in (12), (10), and the dimensionless energy change is given in (55).

<figure><img src="image_13.png" /><figcaption>FIG. 12. Density jump $\mathcal{R}_s$ as a function of $\mathcal{P}_s$ (left) and $\mathcal{T}_s$ (right) for $\gamma = 5/3$, $M_1 = (1, 48)$. Conditions chosen to be similar to [6].</figcaption></figure>

In Fig.12 it is computed the RH jump conditions for a broad range of incident Mach numbers. The effect of ionization is clearly distinguished, as it enhances the gas compression
in the initial stage. As expected, the degree of compression reaches by ionizing shocks exceed the classical limit of $R_s = 4$ [78]. The value and position of the compression ratio peak agrees with the result experimentally obtained by Glass and Liu [70] and numerically by Refs. [6, 79]. For sufficiently high pressures, the density jump drops towards its adiabatic classical limit, which translates into a reversion in the RH curve, or $\beta_{RH} > \pi/2$ as in Fig.2. The inflection point placed in the high-pressure region separates two distinguished regimes: on the left hand side ionization dominates while on the right hand side ionization becomes negligible as it can be taken accomplished.

## Appendix B: Derivation of the diatomic-gas dissociation energy equation

The interest is placed in shock-driven transformations. Stream-wise integration of the energy conservation equation across the dissociation shock gives

$$ \frac{\gamma_1}{\gamma_1 - 1} \frac{p_1}{\rho_1} + \frac{u_1^2}{2} = \frac{p_2}{\rho_2} + \frac{u_2^2}{2} + (1 - \alpha_d) \left( \frac{1}{\gamma_2 - 1} \frac{p_2}{\rho_2} + \frac{T_v N_A k_B}{W} \frac{1}{e^{T_v/T_2} - 1} \right) \\ + \alpha_d \left( \frac{2}{\gamma_3 - 1} \frac{p_2}{\rho_2} + \frac{I_d N_A}{W} \right) \qquad (B1) $$

where it has been considered the situation where the excitation of the vibration and dissociation of diatomic molecules dominate the inelastic process. The factor $\alpha_d$ refers to the dissociation degree, $I_d$ is the Bond-dissociation energy, and $T_v = (2\pi\hbar\nu)/k_B$ is the threshold temperature for the vibrational excitation mode, where $\nu$ is the characteristic oscillation frequency of the molecule. Ideally, the adiabatic indices refer to $\gamma_1 = 7/5$ and $\gamma_3 = 5/3$, but $\gamma_2$ is a function of temperature with the two limiting cases for $\gamma_2$ being 7/5 and 9/7, depending on the degree of excitation of the vibrational modes. This equation can be easily extended to include ionization by adding the corresponding term proportional to $\alpha_d\alpha_i$, provided that ionization is a process that usually initiates after dissociation. This effect is omitted here for the sake of conciseness.

Assuming that the characteristic time of dissociation is much smaller than the period of the perturbations, the dissociation degree is determined by the Saha equation [8, 10]

$$ \frac{\alpha_d^2}{1 - \alpha_d} = 4B_d \frac{\sqrt{T_s}}{R_s} \left[ 1 - \exp\left(-\frac{\phi_v}{T_s}\right) \right] \exp\left(-\frac{\phi_d}{T_s}\right) \qquad (B2) $$

where

$$ B_d = \frac{G_d \phi_r}{4} m_A^{5/2} \frac{(\pi k_B)^{3/2}}{\hbar^3} \frac{T_1^{3/2}}{\rho_1}, \quad \phi_r = \frac{T_r}{T_1}, \quad \phi_v = \frac{T_v}{T_1} \quad \text{and} \quad \phi_d = \frac{I_d}{k_B T_1} \qquad (B3) $$
are the independent dimensionless parameters that characterize the diatomic gas. The pre-exponential factor $G_d = g_{0a}^2/g_{0m}$ where $g_{0a}$ and $g_{0m}$ are the statistical weights of the ground states of atoms and diatomic molecules, respectively, the mass of the atom $m_A$ and the characteristic temperature for rotational excitation $T_r$ through $\phi_r$. The factors $\phi_v$ and $\phi_d$ are the dimensionless characteristic vibrational and dissociation temperatures, respectively. It is found $B_d \sim 3.5 \ 10^6$, $\phi_r \sim 0.01$, $\phi_v \sim 11.3$ and $\phi_d \sim 380$ for N$_2$ at $T_1 = 300$ K at $p_1 = 1$ atm. As a matter of example, $T_r = 87.5$ K, $T_v = 5986$ K, and $T_d = 52438$ K for H$_2$ [8].

The degree of dissociation can be easily reduced to dimensionless variables to yield $\alpha_d = 2\sqrt{\psi_d}(\sqrt{\psi_d} + 1 - \sqrt{\psi_d})$, where the function

$$ \psi_d(\mathcal{T}_s, \mathcal{R}_s) = B_d \frac{\sqrt{\mathcal{T}_s}}{\mathcal{R}_s} \left[ 1 - \exp\left(-\frac{\phi_v}{\mathcal{T}_s}\right) \right] \exp\left(-\frac{\phi_d}{\mathcal{T}_s}\right) \quad (B4) $$

includes dependence with the temperature and density jumps across the shock, $\mathcal{T}_s$ and $\mathcal{R}_s$, respectively. In contrast to ionization shocks, the factor $\mathcal{T}_s$ enters with the power 1/2 and two exponential factors appear. The one associated to dissociation kills the function $\psi_d$ at low temperatures, while the one associated to vibration does likewise at very high temperatures. Once the system of equations is solved, and then the degree of dissociation is given by a relation with the shock strength $\alpha_d(M_1)$, the dimensionless energy per unit of mass subtracted to the fluid particles can be computed, see (56).

## Appendix C: Derivation of the electromagnetic radiation energy equation

With negligible radiation energy and pressure, the first two RH equations (2) and (3) are unchanged but (4) becomes

$$ \rho_1 u_1 \left( \frac{\gamma}{\gamma - 1} \frac{p_1}{\rho_1} + \frac{u_1^2}{2} \right) + F_r = \rho_2 u_2 \left( \frac{\gamma}{\gamma - 1} \frac{p_2}{\rho_2} + \frac{u_2^2}{2} \right), \quad (C1) $$

where the minus sign in $F_r$ indicates that fluid particles loss energy, since $\rho_1 u_1 = \rho_2 u_2 > 0$. According to the definition in (11), the radiative dimensionless energy

$$ Q = -f \frac{\gamma^2 - 1}{2\rho_1 a_1^3} \frac{1}{M_1} \sigma T_2^4 \quad (C2) $$

is a definite negative function. In addition to the explicit dependence with respect to $M_1$ in $Q$, the Mach number will also appear implicitly through $T_2$.

Equation (C1) is nonlinear and it provides the dependence of the dimensionless temperature $\mathcal{T}_s = T_2/T_1$ in terms of the Mach number $M_1$. This dependence can be obtained by
the proper manipulation of the conservation equations. Continuity equation (2) is inserted in (3) and in (C1) to give

$$ \mathcal{T}_s - \mathcal{R}_s^{-1} [1 + \gamma(1 - \mathcal{R}_s^{-1})\mathcal{M}_1^2] = 0, \quad (C3) $$

$$ (\gamma - 1)\phi \mathcal{T}_s^4 + 2\mathcal{M}_1 \mathcal{T}_s + [(\gamma - 1)(\mathcal{R}_s^{-2} - 1)\mathcal{M}_1^2 - 2] \mathcal{M}_1 = 0, \quad (C4) $$

as the system of equations for the two unknowns $\mathcal{R}_s^{-1}$ and $\mathcal{T}_s$, where the positive dimensionless constant $\phi$ reads

$$ \phi = 2f \frac{\sigma \mathcal{T}_1^4}{\rho_1 a_1^3}. \quad (C5) $$

Simple manipulation renders a polynomial of degree nine for $\mathcal{R}_s^{-1}$. However, if the linear term in (C4) is replaced by (C3), a simpler expression for temperature, namely

$$ \mathcal{T}_s^4 = \frac{1}{\phi} \frac{\gamma + 1}{\gamma - 1} (1 - \mathcal{R}_s^{-1}) (\mathcal{R}_{s,ad}^{-1} - \mathcal{R}_s^{-1}) \mathcal{M}_1^3 \quad (C6) $$

is obtained, where the factor $\mathcal{R}_{s,ad}$, which corresponds to the adiabatic condition given in (12) for $Q = 0$, is conveniently introduced. Equation (C6) is found to be very instructive because it allows the derivation of the major contribution in the dependence on $\mathcal{M}_1$ for the scaling of $\mathcal{T}_s$ in terms of the Mach number, irrespectively of the value of $\mathcal{R}_s^{-1}$. Indeed, when $\mathcal{M}_1$ varies, $\mathcal{R}_s^{-1}$ remains in the very small range $0 < \mathcal{R}_s^{-1} < \mathcal{R}_{s,ad}^{-1}$, and, therefore, its net contribution in (C6) is almost insensitive to the value of $\mathcal{M}_1$. It is readily seen that radiation loses produce an over-densification at the shock front, thereby rendering $1 < \mathcal{R}_{s,ad} < \mathcal{R}_s$ to ensure a positive value to the right hand side of (C6). Once $\mathcal{T}_s(\mathcal{M}_1)$ is obtained from (C6) and inserted in (C3), the ratio $\mathcal{R}_s(\mathcal{M}_1) \sim [1 + \gamma(1 - \mathcal{R}_s^{-1})\mathcal{M}_1^2]/\mathcal{T}_s(\mathcal{M}_1)$ provides the dependence of the density jump with $\mathcal{M}_1$, provided that variations of $1 - \mathcal{R}_s^{-1}$ are small.

To further simplify the analysis, the strong shock approximation, $\mathcal{M}_1 \gg 1$, can be justifiably applied, thereby allowing the writing of the downstream quantities in terms of power laws of $\mathcal{M}_1$. Anticipating that $\mathcal{R}_s \gg \mathcal{R}_{s,ad}$ and replacing the latter by its strong-shock asymptotic value $(\gamma + 1)/(\gamma - 1)$, equation (C6) gives immediately $\mathcal{T}_s^4 \sim \mathcal{M}_1^3/\phi$ yielding $\mathcal{T}_s(\mathcal{M}_1) \sim \phi^{-1/4} \mathcal{M}_1^{3/4}$, with the power $3/4 < 1$ being the result of the net contribution of the adiabatic shock $\mathcal{T}_s \propto \mathcal{M}_1^2$ and the counter-case scenario for fully RS, $\mathcal{T}_s \propto \mathcal{M}_1^{1/2}$, where both radiation energy and the radiation pressure are taken into account [43]. It is not surprising that the RS with the radiative flux only leads to a temperature growth with $\mathcal{M}_1$ slower than for the adiabatic case but faster than for the fully radiative regime.
Using the strong-shock limit derived for post-shock temperature in (C3) renders $R_s \sim \gamma \phi^{1/4}(1 - R_s^{-1}) M_1^{5/4}$, which can be further reduced to $R_s(M_1) \sim \gamma \phi^{1/4} M_1^{5/4}$ for $R_s \gg 1$. Unlike adiabatic shocks, the mass-compression ratio grows unbounded for $M_1 \gg 1$, thereby rendering $1 - R_s^{-1} \sim 1$ and $R_{s,ad}^{-1} - R_s^{-1} \sim (\gamma - 1)/(\gamma + 1)$, as previously assumed. As a direct consequence, the downstream velocity goes to zero and matter accumulates just behind the shock front, resulting in a diverging mass density. On the other hand, the pressure ratio is $P_s = R_s T_s \sim \gamma M_1^2$, where it is noticed that $\phi$ is missing, thereby being the same as that for an adiabatic shock.

<figure><img src="image_14.png" /><figcaption>FIG. 13. Density $R_s$ (left) temperature $T_s$ (right) jumps as a function of $M_1$.</figcaption></figure>

The functions $R_s$ and $T_s$ are plotted in Fig.13 as a function of the shock strength $M_1$ for both whole-range (solid) and asymptotic strong-shock (dashed) conditions. By simple inspection in the corresponding strong-shock limits, it is observed that $M_1 \sim \phi^{-1/5}$ marks the region where functions $T_s$ and $R_s$ detach from the curve ($\phi = 0$) towards their corresponding asymptotic trends derived in the previous paragraph.

## Appendix D: Derivation of endothermic and exothermic nuclear reactions

When considering nuclear dissociation in CCSN shocks as an isolated phenomenon, and assuming nuclear statistical equilibrium, the mass fraction of each of the free neutrons, free
protons, and alpha particles is given by the nuclear Saha equation. For simplicity, it is further assumed only alpha particles are present in the cold upstream flow, while in the post-shock flow, due to high temperatures, a fraction of alpha particles are dissociated into free protons and neutrons. This is indeed a simplified picture, although the qualitative aspects of the dissociation energy are captured sufficiently well. A simple justification is given by the fact that the binding energy of nucleons in alpha particles is $E_{\text{bin}} \sim 7.075$ MeV per baryon, which is within $\sim 20\%$ of other nuclei present in flow such as oxygen or iron.

The mass fraction of free nucleons can be obtained from the Saha equation [80]. Using the fitting function of [81], the approximate solution

$$X_{\alpha} = B_{nd1} \left( \frac{T_s^3}{R_s^2} \right)^{3/8} \exp \left( -\frac{\phi_{nd1}}{T_s} \right) \quad (\text{D1})$$

can be used to compute the mass fraction of alpha particles on condition that $X_{\alpha} < 1$, since (D1) grows unbounded for $T_s \gg 1$. The coefficients are

$$\phi_{nd1} = \frac{7.075}{T_{1,\text{MeV}}}, \quad B_{nd1} \sim 2.8 \left( \frac{T_{1,\text{MeV}}^3}{\rho_{1,9}^2} \right)^{3/8}, \quad (\text{D2})$$

where $\rho_{1,9}$ is the pre-shock density in units of $10^9$ g/cm$^3$ and $T_{1,\text{MeV}}$ is the temperature in the pre-shock region in units of MeV. For accretion shocks in CCSNe, $\rho_{1,9}$ and $T_{1,\text{MeV}}$ are $\sim 1$. The dimensionless dissociation energy is then $Q = -X_{\alpha}\phi_{nd1}(\gamma^2 - 1)/(2\gamma)$. The exponential dependence indicates that nuclear dissociation does occur effectively for temperatures sufficiently large, i.e., when $\phi_{nd1}/T_s \sim 1$ or smaller. It is then straightforward to relate it with the energy release sensitivity relative to the shock strength $\epsilon$. Notice that energy subtracted from by nuclear dissociation is a complex phenomena that depends non-linearly on the changing upstream conditions. This process is typically simplified to be proportional to the hydrodynamic incoming energy through the Bernoulli parameter [49–51]. In such cases, the planar shock is found to be always stable [52].

Calculating the energy release in thermonuclear SNe is complicated for a number of reasons, including the large number of reactions involved and coupling to hydrodynamics [45]. A rough estimate can be obtained in the following way. For simplicity, it is assumed that the nuclear fuel burns into $^{56}$Ni at the shock. For white dwarfs that mainly consist of $^{12}$C and $^{16}$O before explosion, this process releases 0.79 MeV per baryon (or 44.4 MeV per $^{56}$Ni nucleus), which is just the difference in the binding energies of $^{56}$Ni and two $^{12}$C and $^{16}$O nuclei. At high temperatures, $^{56}$Ni may endothermically dissociate into smaller nuclei,
a process that can take up to 1.6 MeV per baryon, when $^{56}$Ni split into alpha particles. Heavier sub-products will lead to smaller energy subtractions. Therefore, the net energy balance across the thermonuclear SNe must account for both the coupled endothermic and exothermic processes coexisting behind the shock. For $X_{\alpha} \ll 1$, the corresponding Saha equation for the dissociation process can approximated as

$$X_{\alpha} \sim B_{nd2} \left( \frac{T_s^3}{R_s^2} \right)^{13/28} \exp \left( -\frac{\phi_{nd2}}{T_s} \right), \qquad (D3)$$

where

$$\phi_{nd2} = \frac{1.6}{T_{1,\text{MeV}}}, \quad B_{nd2} \sim 3876 \left( \frac{T_{1,\text{MeV}}^3}{\rho_{1,9}^2} \right)^{13/28}, \qquad (D4)$$

and $\rho_{1,9}$ is the pre-shock density in units of $10^9$ g/cm$^3$. In the context of thermonuclear supernovae, $\rho_{1,9}$ and $T_{1,\text{MeV}}$ are roughly of the order of 1, so $\phi_{nd2}$ is the order of unity while $B_{nd2} \sim 10^3 - 10^4$. The corresponding energy releases is offered in (58).

[1] N. Freeman, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences **228**, 341 (1955).

[2] P. Zaidel', Journal of Applied Mathematics and Mechanics **24**, 316 (1960).

[3] G. Fowles and G. Swan, Physical Review Letters **30**, 1023 (1973).

[4] G. Fowles, The Physics of Fluids **24**, 220 (1981).

[5] G. Fraley, The Physics of Fluids **29**, 376 (1986).

[6] M. Mond and I. Rutkevich, Journal of Fluid Mechanics **275**, 121 (1994).

[7] M. Mond, I. Rutkevich, and E. Toffin, Physical Review E **56**, 5968 (1997).

[8] J. W. Bates and D. C. Montgomery, Physics of Fluids **11**, 462 (1999).

[9] J. W. Bates and D. C. Montgomery, Physical Review Letters **84**, 1180 (2000).

[10] M. Mond and I. Rutkevich, Physics of Fluids **14**, 1468 (2002).

[11] J. W. Bates, Physical Review E **69**, 056313 (2004).

[12] J. Wouchuk and J. L. Cavada, Physical Review E **70**, 046303 (2004).

[13] J. W. Bates, Physics of Fluids **19**, 094102 (2007).

[14] A. V. Konyukhov, A. P. Likhachev, V. E. Fortov, K. V. Khishchenko, S. I. Anisimov, A. M. Oparin, and I. V. Lomonosov, JETP letters **90**, 18 (2009).

[15] J. Bates, Journal of fluid mechanics **691**, 146 (2012).
[16] J. Bates, Physical Review E **91**, 013014 (2015).

[17] N. Alferez and E. Touber, Journal of Fluid Mechanics **814**, 185 (2017).

[18] N. Wetta, J.-C. Pain, and O. Heuzé, Physical Review E **98**, 033205 (2018).

[19] E. Touber and N. Alferez, Journal of Fluid Mechanics **864**, 807 (2019).

[20] S. D'yakov, Zh. Eksp. Teor. Fiz **27**, 288 (1954).

[21] V. Kontorovich, Zh. Eksp. Teor. Fiz **33**, 1525 (1957).

[22] R. D. Richtmyer, Communications on pure and applied mathematics **13**, 297 (1960).

[23] E. Meshkov, Fluid Dynamics **4**, 101 (1969).

[24] A. L. Velikovich, Physics of Fluids **8**, 1666 (1996).

[25] M. Brouillette, Annual Review of Fluid Mechanics **34**, 445 (2002).

[26] K. Nishihara, J. Wouchuk, C. Matsuoka, R. Ishizaki, and V. Zhakhovsky, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences **368**, 1769 (2010).

[27] J. Wouchuk, Physical Review E **63**, 056303 (2001).

[28] J. Wouchuk, Physics of Plasmas **8**, 2890 (2001).

[29] F. Cobos-Campos and J. Wouchuk, Physical Review E **93**, 053111 (2016).

[30] F. Cobos-Campos and J. Wouchuk, Physical Review E **96**, 013102 (2017).

[31] J. J. Erpenbeck, The Physics of Fluids **7**, 684 (1964).

[32] G. Joulin and P. Vidal, Collection Alea Saclay, Monographs and texts in statistical physics **1**, 493 (1998).

[33] M. Short and D. S. Stewart, Journal of Fluid Mechanics **382**, 109 (1999).

[34] R. Daou and P. Clavin, Journal of Fluid Mechanics **482**, 181 (2003).

[35] P. Clavin and F. A. Williams, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences **370**, 597 (2012).

[36] P. Clavin and G. Searby, *Combustion Waves and Fronts in Flows: Flames, Shocks, Detona-tions, Ablation Fronts and Explosion of Stars* (Cambridge University Press, 2016).

[37] C. Huete, T. Jin, D. Martínez-Ruiz, and K. Luo, Physical Review E **96**, 053104 (2017).

[38] C. Huete and M. Vera, Journal of Fluid Mechanics **879**, 54 (2019).

[39] H. Smolders and M. Van Dongen, Shock Waves **2**, 255 (1992).

[40] A. Guha, Physics of fluids **6**, 1893 (1994).

[41] N. Smirnov, V. Betelin, A. Kushnirenko, V. Nikitin, V. Dushin, and V. Nerchenko, Acta
Astronautica **87**, 14 (2013).

[42] D. Vimercati, G. Gori, and A. Guardone, Journal of Fluid Mechanics **847**, 266 (2018).

[43] S. Bouquet, R. Teyssier, and J. P. Chieze, The Astrophysical Journal Supplement Series **127**, 245 (2000).

[44] A. Khokhlov, Monthly Notices of the Royal Astronomical Society **239**, 785 (1989).

[45] W. D. Arnett, *Supernovae and nucleosynthesis: An investigation of the history of matter, from the Big Bang to the present*, Vol. 7 (Princeton University Press, 1996).

[46] A. Khokhlov, E. S. Oran, and J. C. Wheeler, The Astrophysical Journal **478**, 678 (1997).

[47] V. N. Gamezo, A. M. Khokhlov, and E. S. Oran, Physical Review Letters **92**, 211102 (2004).

[48] F. K. Röpke and S. A. Sim, Space Science Reviews **214**, 72 (2018).

[49] R. Fernández and C. Thompson, The Astrophysical Journal **697**, 1827 (2009).

[50] R. Fernández and C. Thompson, The Astrophysical Journal **703**, 1464 (2009).

[51] C. Huete, E. Abdikamalov, and D. Radice, Monthly Notices of the Royal Astronomical Society **475**, 3305 (2018).

[52] C. Huete and E. Abdikamalov, Physica Scripta **94**, 094002 (2019).

[53] J. M. Blondin, A. Mezzacappa, and C. DeMarino, The Astrophysical Journal **584**, 971 (2003).

[54] T. Foglizzo, P. Galletti, L. Scheck, and H.-T. Janka, The Astrophysical Journal **654**, 1006 (2007).

[55] H.-T. Janka, Annual Review of Nuclear and Particle Science **62**, 407 (2012).

[56] A. Burrows, Reviews of Modern Physics **85**, 245 (2013).

[57] B. Müller, M. Viallet, A. Heger, and H.-T. Janka, The Astrophysical Journal **833**, 124 (2016).

[58] R. P. Drake, in *High Energy Density Laboratory Astrophysics* (Springer, 2005) pp. 49–59.

[59] R. P. Drake, *High-energy-density physics* (Springer-Verlag, 2006).

[60] R. P. Drake, IEEE Transactions on plasma science **35**, 171 (2007).

[61] R. P. Drake, Physics of Plasmas **14**, 043301 (2007).

[62] A. Reighard and R. Drake, Astrophysics and Space Science **307**, 121 (2007).

[63] C. Michaut, E. Falize, C. Cavet, S. Bouquet, M. Koenig, T. Vinci, A. Reighard, and R. Drake, Astrophysics and Space Science **322**, 77 (2009).

[64] R. G. McClarren, R. P. Drake, J. Morel, and J. P. Holloway, Physics of Plasmas **17**, 093301 (2010).

[65] D. Ryu and E. Vishniac, The Astrophysical journal **313**, 820 (1987).
[66] E. Vishniac and D. Ryu, The Astrophysical journal **337**, 917 (1989).

[67] B. Ahlborn and M. Salvat, Zeitschrift für Naturforschung A **22**, 260 (1967).

[68] J. Strachan and B. Ahlborn, Australian Journal of Physics **28**, 395 (1975).

[69] H. Nieuwenhuijzen, C. De Jager, M. Cuntz, A. Lobel, and L. Achmad, Astronomy and Astro-physics **280**, 195 (1993).

[70] I. Glass and W. Liu, Journal of Fluid Mechanics **84**, 55 (1978).

[71] I. Glass, W. Liu, and F. Tang, Canadian Journal of Physics **55**, 1269 (1977).

[72] Y. B. Zel'dovich and Y. P. Raizer, *Physics of Shock Waves and High-Temperature Hydrody-namic Phenomena* (Dover Publications, 2002).

[73] N. M. Kuznetsov, Sov. Phys. Usp. **32**, 993 (1989).

[74] H. A. Bethe, The Astrophysical Journal **490**, 765 (1997).

[75] B. T. Draine and C. F. McKee, Annual review of astronomy and astrophysics **31**, 373 (1993).

[76] G. Barenblatt, R. Guirguis, M. Kamel, A. Kuhl, A. Oppenheim, and Y. B. Zel'Dovich, Journal of Fluid Mechanics **99**, 841 (1980).

[77] G. Barenblatt, *Scaling, self-similarity, and intermediate asymptotics: dimensional analysis and intermediate asymptotics*, Vol. 14 (Cambridge University Press, 1996).

[78] M. A. Liberman and A. L. Velikovich, *Physics of shock waves in gases and plasmas* (Springer Verlag Springer Series on Electrophysics, 1986).

[79] I. Rutkevich and M. Mond, Journal of Propulsion and Power **10**, 906 (1994).

[80] L. S. Shapiro and S. A. Teukolsky, *Black Holes, White Dwarfs and Neutron Stars* (John Wiley & Sons, New York U. S. A., 1983).

[81] R. Fernández, A. Tchekhovskoy, E. Quataert, F. Foucart, and D. Kasen, Mon. Not. Roy. Astron. Soc. **482**, 3373 (2019).