Weak Shock Propagation with Accretion. II. Stability of Self-similar Solutions to Radial Perturbations

The generation and propagation of a shock wave is central to many otherwise-distinct, astrophysical phenomena, including supernovae, star formation, and galaxy evolution. One of the most powerful methods for describing the time- and space-dependent evolution of a fluid in the presence of a shock is that of self-similarity (see Ostriker & McKee 1988 for a review). This mathematical technique yields exact solutions to the fluid equations—a set of coupled, nonlinear, partial differential equations—under a set of basic and nonrestrictive assumptions. Moreover, self-similar solutions reproduce the properties of shocks to a high degree of accuracy in a wide range of contexts: the Sedov–Taylor blast wave (Taylor 1950; Sedov 1959), perhaps one of the best-known examples of a self-similar solution to the fluid equations, can be used to predict the energy released from an atomic bomb and a supernova explosion alike.

Recently, Coughlin et al. (2018b, hereafter Paper I) described a new mode of self-similar shock propagation, which, dissimilar from the Sedov–Taylor blast wave, is valid when the shock Mach number is only marginally greater than 1 and the fluid is situated in the gravitational field of a point mass. These solutions predict not only the outward motion of gas immediately behind the shock but also the accretion of matter onto the point mass. The energy budget of the shocked fluid is also modified by the binding energy of the ambient medium and the accretion of binding energy by the black hole, neither of which is present in the Sedov–Taylor blast wave. These solutions were also shown to accurately reproduce the propagation of a shock in the hydrogen envelope of a massive star following a failed supernova, in which a weak shock is still generated from the mass radiated in neutrinos during core collapse (Nadezhin 1980; Lovegrove & Woosley 2013; Coughlin et al. 2018a; Fernández et al. 2018).

While the utility of self-similar solutions has been demonstrated both mathematically and phenomenologically, the question of their global stability is of paramount importance for determining the asymptotic, temporal behavior of shock waves. In particular, while a self-similar solution provides an exact solution to the fluid equations and can accurately reproduce the propagation of a shock in certain contexts, the evolution of deviations from the self-similar solution must also be understood to determine whether or not the self-similar solution characterizes the state of a shock at late times. If small deviations in the fluid variables from the self-similar prescription are found to grow asymptotically, then the self-similar solution is unstable and cannot describe the long-term evolution of a shock.

In this paper we analyze the stability of self-similar shocks to radial perturbations, focusing primarily on the solutions presented in Coughlin et al. (2018a) (though we also analyze the stability of the Sedov–Taylor blast wave in the Appendix). For ease of notation, throughout the remainder of this paper we refer to the self-similar solutions presented in Paper I as the Coughlin, Quataert, & Ro (CQR) solutions. In Section 2 we present the basic equations and develop the formalism for analyzing the evolution of perturbations on top of the CQR solution, and in Section 3 we write the solutions to these equations in the form of eigenmodes. In Section 4 we present the eigenmodes for specific density profiles of the ambient medium, and we show that the CQR solutions are extremely weakly unstable, with perturbations to the shock position growing with time t as t^α with α ≲ 0.12. We compare the results of the eigenmode analysis to simulations in Section 5, showing very good agreement between the two (the simulations will be described in more detail in Paper III of this series). In Section 6, we show that there is a second self-similar solution to the fluid equations when the shock Mach number is exactly 1; this solution corresponds to a rarefaction wave (RW) that travels at the sound speed and likely represents the asymptotic state of shocks with Mach number less than the CQR value. We discuss the physical origin of the instability in Section 7, and we summarize and conclude in Section 8.

In this section, we derive the self-similar solutions that describe a shocked fluid as it expands in the gravitational field of a compact object, and we present the equations that govern the evolution of perturbations on top of those self-similar solutions. Our notation will be very similar to that of Paper I, but we will introduce some small differences that are useful for the perturbation analysis of this paper. For ease of reference, we explicitly point out these differences in Table 1.

Table 1.  A Summary of the Variables Defined in Paper I and Those Introduced Here

Variable	Paper I	This Paper
Time-dependent shock position, velocity	${r}_{\mathrm{sh}},{v}_{\mathrm{sh}}$	R, V
Ambient length scale, density, pressure	${r}_{0},{\rho }_{0},{p}_{0}$	${r}_{{\rm{a}}},{\rho }_{{\rm{a}}},{p}_{{\rm{a}}}$
Self-similar shock velocity	$V\sqrt{{GM}/{r}_{\mathrm{sh}}}$	${V}_{{\rm{c}}}\sqrt{{GM}/R}$
Self-similar post-shock velocity	$v=V\sqrt{{GM}/r}f(\xi )$	$v=V(t){f}_{0}(\xi )$
Self-similar post-shock density	$\rho ={\rho }_{0}{\left(r/{r}_{0}\right)}^{-n}g(\xi )$	$\rho ={\rho }_{{\rm{a}}}{\left(R(t)/{r}_{{\rm{a}}}\right)}^{-n}{g}_{0}(\xi )$
Self-similar post-shock pressure	$p={\rho }_{0}{GM}/r{\left(r/{r}_{0}\right)}^{-n}h(\xi )$	$p={\rho }_{{\rm{a}}}V{\left(t\right)}^{2}{\left(R(t)/{r}_{{\rm{a}}}\right)}^{-n}{h}_{0}(\xi )$

Download table as:  ASCIITypeset image

We assume that the gas is spherically symmetric, in hydrostatic equilibrium, and the gravitational field is dominated by a compact object of mass M. The continuity, radial momentum, and entropy equations that govern the evolution of the fluid velocity, v, density, ρ, and pressure, p, are then

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}\rho v\right]=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial v}{\partial t}+v\displaystyle \frac{\partial v}{\partial r}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial r}=-\displaystyle \frac{{GM}}{{r}^{2}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\mathrm{ln}(s)+v\displaystyle \frac{\partial }{\partial r}\mathrm{ln}(s)=0,\end{eqnarray} \tag{ 3 }$$

where s = p/ρ^γ is the specific entropy of the gas with adiabatic index γ. Further assume that there is a shock propagating through the medium with time-dependent position R(t) and velocity $V(t)={dR}/{dt}$. If the pre-shock gas obeys an adiabatic equation of state and we are sufficiently far from any surface, then the density and pressure of the pre-shock gas, ρ₁ and p₁, respectively, satisfy

$$\begin{eqnarray}&&{\rho }_{1}={\rho }_{{\rm{a}}}{\left(\displaystyle \frac{r}{{r}_{{\rm{a}}}}\right)}^{-n},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{p}_{1}=\displaystyle \frac{1}{n+1}\displaystyle \frac{{GM}}{r}{\rho }_{{\rm{a}}}{\left(\displaystyle \frac{r}{{r}_{{\rm{a}}}}\right)}^{-n},\end{eqnarray} \tag{ 5 }$$

where ρ_a is the density of the ambient medium at radius r_a and n is a free parameter (note that n has often been denoted by ω in past works on the Sedov–Taylor blast wave). The shock jump conditions, which guarantee the continuity of energy, momentum, and mass across the shock, further yield a set of boundary conditions that must be satisfied by the fluid at the shock front. If we denote the pre- and post-shock adiabatic indices by γ₁ and γ₂, respectively, then these boundary conditions are (Ostriker & McKee 1988; Paper I)

$$\begin{eqnarray}\begin{array}{rcl}v(R) & = & \left\{1-\displaystyle \frac{{\gamma }_{2}}{{\gamma }_{2}+1}\left[1+\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right]\right.\\ & & +\,\left.\displaystyle \frac{{\gamma }_{2}-1}{{\gamma }_{2}+1}\sqrt{{\left(\displaystyle \frac{{\gamma }_{2}}{{\gamma }_{2}-1}\right)}^{2}{\left(1+\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right)}^{2}-\left(\displaystyle \frac{{\gamma }_{2}+1}{{\gamma }_{2}-1}\right)\left(1+\displaystyle \frac{2{\gamma }_{1}}{{\gamma }_{1}-1}\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right)}\right\}\times V,\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{rcl}\rho (R) & = & \displaystyle \frac{{\gamma }_{2}+1}{{\gamma }_{2}-1}\left\{\displaystyle \frac{{\gamma }_{2}}{{\gamma }_{2}-1}\left[1+\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right]\right.\\ & & -\,{\left.\sqrt{{\left(\displaystyle \frac{{\gamma }_{2}}{{\gamma }_{2}-1}\right)}^{2}{\left(1+\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right)}^{2}-\left(\displaystyle \frac{{\gamma }_{2}+1}{{\gamma }_{2}-1}\right)\left(1+\displaystyle \frac{2{\gamma }_{1}}{{\gamma }_{1}-1}\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right)}\right\}}^{-1}\times {\rho }_{{\rm{a}}}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n},\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}p(R) & = & \displaystyle \frac{1}{{\gamma }_{2}+1}\left\{1+\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right.\\ & & +\,\left.\left({\gamma }_{2}-1\right)\sqrt{{\left(\displaystyle \frac{{\gamma }_{2}}{{\gamma }_{2}-1}\right)}^{2}{\left(1+\displaystyle \frac{{GM}}{\left(n+1\right){{RV}}^{2}}\right)}^{2}-\left(\displaystyle \frac{{\gamma }_{2}+1}{{\gamma }_{2}-1}\right)\left(1+\displaystyle \frac{2{\gamma }_{1}}{{\gamma }_{1}-1}\displaystyle \frac{{GM}}{{{RV}}^{2}}\right)}\right\}\times {\rho }_{{\rm{a}}}{V}^{2}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}.\end{array}\end{eqnarray} \tag{ 8 }$$

These expressions are for arbitrary n, γ₁, and γ₂, though they simplify considerably when γ₁ = γ₂ and when the shock Mach number is much greater than 1. Here, however, we will make neither of those assumptions.

Since the shock position changes with time, these boundary conditions take place at a temporally evolving radius. This observation suggests that we make the change of variables

$$\begin{eqnarray}&&r\to \xi ,\quad \xi =\displaystyle \frac{r}{R(t)},\end{eqnarray} \tag{ 9 }$$

as the boundary conditions are then satisfied at a single value of ξ for all t. Also, while it is not immediately obvious that it is beneficial to do so, we will make the additional change of variables

$$\begin{eqnarray}&&t\to \tau ,\quad \tau =\mathrm{ln}\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right).\end{eqnarray} \tag{ 10 }$$

We emphasize that the R(t) appearing in these definitions is the true shock position, and we will return to the ramifications of this choice (compared to using the unperturbed, or self-similar, shock position) in later sections. In terms of these variables, Equations (1)–(3) become

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial \tau }-\xi \displaystyle \frac{\partial \rho }{\partial \xi }+\displaystyle \frac{1}{V}\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left({\xi }^{2}\rho v\right)=0,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial v}{\partial \tau }-\xi \displaystyle \frac{\partial v}{\partial \xi }+\displaystyle \frac{1}{V}v\displaystyle \frac{\partial v}{\partial \xi }+\displaystyle \frac{1}{V}\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial \xi }=-\displaystyle \frac{{GM}}{{VR}}\displaystyle \frac{1}{{\xi }^{2}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{ln}s}{\partial \tau }-\xi \displaystyle \frac{\partial \mathrm{ln}s}{\partial \xi }+\displaystyle \frac{1}{V}v\displaystyle \frac{\partial \mathrm{ln}s}{\partial \xi }=0,\end{eqnarray} \tag{ 13 }$$

where temporal derivatives here are with respect to constant ξ. Note that we did not yet assume anything about the temporal dependence of the shock position—we only used the fact that $V={dR}/{dt}$ to derive the above three equations.

We now write the fluid quantities in the following forms:

$$\begin{eqnarray}&&\begin{array}{r}v\,=\,V\left\{{f}_{0}(\xi )+{f}_{1}(\xi ,\tau )\right\},\,\rho \,=\,{\rho }_{{\rm{a}}}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}\left\{{g}_{0}(\xi )+{g}_{1}(\xi ,\tau )\right\},\,p\,=\,{\rho }_{{\rm{a}}}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}{V}^{2}\left\{{h}_{0}(\xi )+{h}_{1}(\xi ,\tau )\right\},\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray*}&&\begin{array}{l}\Rightarrow \,\mathrm{ln}s\simeq \mathrm{ln}{V}^{2}{R}^{-n+n{\gamma }_{2}}+\mathrm{ln}\left(\displaystyle \frac{{h}_{0}}{{g}_{0}^{{\gamma }_{2}}}\right)+\displaystyle \frac{{h}_{1}}{{h}_{0}}-{\gamma }_{2}\displaystyle \frac{{g}_{1}}{{g}_{0}}\equiv \mathrm{ln}{V}^{2}{R}^{-n+n{\gamma }_{2}}+{s}_{0}(\xi )+{s}_{1}(\xi ,\tau ),\end{array}\end{eqnarray*}$$

where in the last line we kept only first-order corrections in the assumed-small ratios ${g}_{1}/{g}_{0}$ and ${h}_{1}/{h}_{0}$. There is also an additional constant that appears in the definition of the entropy that is irrelevant to the equations. With these parameterizations, the fluid equations become

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {g}_{1}}{\partial \tau }-n\left({g}_{0}+{g}_{1}\right)-\xi \displaystyle \frac{\partial }{\partial \xi }\left({g}_{0}+{g}_{1}\right)+\,\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left[{\xi }^{2}\left({g}_{0}+{g}_{1}\right)\left({f}_{0}+{f}_{1}\right)\right]=0,\end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{1}}{\partial \tau }+\left(\displaystyle \frac{1}{2}\displaystyle \frac{d}{d\tau }\mathrm{ln}{V}^{2}\right)\left({f}_{0}+{f}_{1}\right)-\xi \displaystyle \frac{\partial }{\partial \xi }\left({f}_{0}+{f}_{1}\right)+\left({f}_{0}+{f}_{1}\right)\displaystyle \frac{\partial }{\partial \xi }\left({f}_{0}+{f}_{1}\right)+\displaystyle \frac{1}{{g}_{0}}{\left(1+\displaystyle \frac{{g}_{1}}{{g}_{0}}\right)}^{-1}\displaystyle \frac{\partial }{\partial \xi }\left({h}_{0}+{h}_{1}\right)=-\displaystyle \frac{{GM}}{{{RV}}^{2}{\xi }^{2}},\end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {s}_{1}}{\partial \tau }+\displaystyle \frac{d}{d\tau }\mathrm{ln}{V}^{2}-n+n{\gamma }_{2}-\xi \displaystyle \frac{\partial }{\partial \xi }\left({s}_{0}+{s}_{1}\right)+\left({f}_{0}+{f}_{1}\right)\displaystyle \frac{\partial }{\partial \xi }\left({s}_{0}+{s}_{1}\right)=0.\end{array}\end{eqnarray} \tag{ 17 }$$

These equations can be satisfied exactly with all subscript-1 quantities identically zero if

$$\begin{eqnarray}&&\displaystyle \frac{{GM}}{{{RV}}^{2}}=\displaystyle \frac{1}{{V}_{{\rm{c}}}^{2}}\quad \iff \quad V={V}_{{\rm{c}}}\sqrt{\displaystyle \frac{{GM}}{R}},\end{eqnarray} \tag{ 18 }$$

where V_c is an unspecified parameter. If the shock velocity satisfies precisely this scaling, then the solutions to the fluid equations (and the boundary conditions) are self-similar, i.e., the solutions can be written as separable in t and ξ, and we recover the following three ordinary differential equations for the unperturbed quantities:

$$\begin{eqnarray}&&-{{ng}}_{0}-\xi {g}_{0}^{{\prime} }+\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{d}{d\xi }\left({\xi }^{2}{g}_{0}{f}_{0}\right)=0,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&-\displaystyle \frac{1}{2}{f}_{0}-\xi {f}_{0}^{{\prime} }+{f}_{0}{f}_{0}^{{\prime} }+\displaystyle \frac{1}{{g}_{0}}{h}_{0}^{{\prime} }=-\displaystyle \frac{1}{{V}_{{\rm{c}}}^{2}{\xi }^{2}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&n{\gamma }_{2}-n-1+\left({f}_{0}-\xi \right){s}_{0}^{{\prime} }=0.\end{eqnarray} \tag{ 21 }$$

These three equations are the same as those derived in Appendix A of Paper I, and the boundary conditions at ξ = 1 can be read off from Equations (6)–(8). The parameter V_c, as shown in Paper I, is determined by the Requirement that the solutions smoothly pass through a critical point and is of the order 1. These shocks are therefore only mildly supersonic, and they result in both outward motion immediately behind the shock and accretion onto the point mass at the origin. We refer the reader to Paper I for a much more detailed discussion of the properties of the self-similar solutions.

Importantly, this analysis demonstrates that deviations from the self-similar velocity prescription induce finite perturbations to the fluid quantities. In other words, if the shock velocity varies as anything other than Equation (18), then perturbations to the fluid quantities must exist in order to satisfy the fluid equations. We therefore parameterize the shock velocity by

$$\begin{eqnarray}&&\displaystyle \frac{{GM}}{{{RV}}^{2}}=\displaystyle \frac{1}{{V}_{{\rm{c}}}^{2}}\left[1-2\zeta \left(\tau \right)\right],\end{eqnarray} \tag{ 22 }$$

where ζ(τ) is assumed to be a small correction that encodes the deviation of the shock velocity from the self-similar one (the sign convention and the factor of 2 will be explained in the next section). Keeping only first-order terms in ζ in Equations (15)–(17) then gives the following three equations for the perturbations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {g}_{1}}{\partial \tau }-{{ng}}_{1}-\xi \displaystyle \frac{\partial {g}_{1}}{\partial \xi }+\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left({\xi }^{2}{g}_{0}{f}_{1}+{\xi }^{2}{f}_{0}{g}_{1}\right)=0,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{1}}{\partial \tau }-\displaystyle \frac{1}{2}{f}_{1}-\xi \displaystyle \frac{\partial {f}_{1}}{\partial \xi }+\displaystyle \frac{\partial }{\partial \xi }\left({f}_{0}{f}_{1}\right)+\displaystyle \frac{1}{{g}_{0}}\displaystyle \frac{\partial }{\partial \xi }\left({h}_{0}{s}_{1}+{\gamma }_{2}{h}_{0}\displaystyle \frac{{g}_{1}}{{g}_{0}}\right)-\displaystyle \frac{{g}_{1}}{{g}_{0}^{2}}\displaystyle \frac{\partial {h}_{0}}{\partial \xi }=-{f}_{0}\displaystyle \frac{d\zeta }{d\tau }+\displaystyle \frac{2}{{V}_{{\rm{c}}}^{2}{\xi }^{2}}\zeta ,\end{array}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {s}_{1}}{\partial \tau }+\left({f}_{0}-\xi \right)\displaystyle \frac{\partial {s}_{1}}{\partial \xi }+{f}_{1}{s}_{0}^{{\prime} }=-2\displaystyle \frac{d\zeta }{d\tau }.\end{eqnarray} \tag{ 25 }$$

We can cast these equations in a slightly simplified form if we further define

$$\begin{eqnarray}&&{F}_{1}\equiv {g}_{0}{\xi }^{2}{f}_{1},\quad {G}_{1}\equiv \displaystyle \frac{{g}_{1}}{{g}_{0}},\end{eqnarray} \tag{ 26 }$$

which are the perturbation to the mass flux and the relative perturbation to the density, respectively. Written in matrix form, our final form for the perturbation equations is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial \tau }\left(\begin{array}{c}{s}_{1}\\ {F}_{1}\\ {G}_{1}\end{array}\right)\,+\,\left(\begin{array}{ccc}{f}_{0}-\xi & 0 & 0\\ {\xi }^{2}{h}_{0} & {f}_{0}-\xi & {\gamma }_{2}{\xi }^{2}{h}_{0}\\ 0 & {\xi }^{-2}{g}_{0}^{-1} & {f}_{0}-\xi \end{array}\right)\displaystyle \frac{\partial }{\partial \xi }\left(\begin{array}{c}{s}_{1}\\ {F}_{1}\\ {G}_{1}\end{array}\right)+\left[\begin{array}{ccc}0 & {\xi }^{-2}{g}_{0}^{-1}{s}_{0}^{{\prime} } & 0\\ {\xi }^{2}{h}_{0}^{{\prime} } & \displaystyle \frac{3}{2}-n+2{f}_{0}^{{\prime} } & {\xi }^{2}({\gamma }_{2}-1){h}_{0}^{{\prime} }\\ 0 & 0 & 0\end{array}\right]\,\left(\begin{array}{c}{s}_{1}\\ {F}_{1}\\ {G}_{1}\end{array}\right)\,=\,\left(\begin{array}{c}-2\dot{\zeta }\\ \displaystyle \frac{2{g}_{0}}{{V}_{{\rm{c}}}^{2}}\zeta -{f}_{0}{g}_{0}{\xi }^{2}\dot{\zeta }\\ 0\end{array}\right),\end{array}\end{eqnarray} \tag{ 27 }$$

where primes denote differentiation with respect to ξ and dots with respect to τ. The boundary conditions for s₁, F₁, and G₁ are determined by the shock jump conditions, specifically Equations (6)–(8), by inserting our expression for the shock velocity (Equation (22)) into the right-hand sides, Taylor-expanding about ζ = 0, and keeping only first-order terms. The resulting expressions are lengthy in general, but they are of the form

$$\begin{eqnarray}\begin{array}{rcll}{s}_{1}(1,\tau ) & = & {a}_{{\rm{s}}}\zeta (\tau ),\ {F}_{1}(1,\tau )={a}_{{\rm{F}}}\zeta (\tau ),\ {G}_{1}(1,\tau )\,= & {a}_{{\rm{G}}}\zeta (\tau ),\end{array}\end{eqnarray} \tag{ 28 }$$

where the a are constants that depend on γ₁, γ₂, and n (and V_c, though this is also a function of γ₁, γ₂, and n). For the specific case where γ₁ = γ₂ = 1 + 1/n, which is a good approximation when the shock is propagating through the hydrogen envelope of a massive star (see the discussion at the end of Section 2 of Paper I), these constants reduce to

$$\begin{eqnarray}\begin{array}{rclcl}{a}_{{\rm{s}}} & = & -\displaystyle \frac{2\left(1+2n\right)\left(-2+3{{nV}}_{{\rm{c}}}^{2}+2{n}^{2}{V}_{{\rm{c}}}^{2}\right)}{n\left(2+{V}_{{\rm{c}}}^{2}\right)\left(-1+2{{nV}}_{{\rm{c}}}^{2}+2{n}^{2}{V}_{{\rm{c}}}^{2}\right)},\,{a}_{{\rm{F}}} & = & {a}_{{\rm{G}}}=\displaystyle \frac{4}{2+{V}_{{\rm{c}}}^{2}}.\end{array}\end{eqnarray} \tag{ 29 }$$

Equation (27) appears underconstrained, as the number of equations (three) is fewer than the number of unknowns (four). The resolution to this apparent issue is that the determinant of the matrix multiplying the derivatives in Equation (27) is

$$\begin{eqnarray}\begin{array}{l}\det \left(\begin{array}{ccc}{f}_{0}-\xi & 0 & 0\\ {\xi }^{2}{h}_{0} & {f}_{0}-\xi & {\gamma }_{2}{\xi }^{2}{h}_{0}\\ 0 & {\xi }^{-2}{g}_{0}^{-1} & {f}_{0}-\xi \end{array}\right)=\left({f}_{0}-\xi \right)\left[{\left({f}_{0}-\xi \right)}^{2}-\displaystyle \frac{{\gamma }_{2}{h}_{0}}{{g}_{0}}\right].\end{array}\end{eqnarray} \tag{ 30 }$$

As shown in Paper I, the second factor in this expression equals zero at a critical point ξ_c and characterizes the location in the unperturbed, self-similar solution where the gas goes from sub- to supersonic infall (note that, because of the background flow and the expansion of the shock, it is the relative velocity f₀ − ξ that determines the location of the sonic point). We see that this position also marks a special location for the perturbations, and if the perturbations are to remain finite at the critical point,⁴ there is a fourth boundary condition that must be satisfied by the functions at this location within the flow. This fourth boundary condition serves as an additional constraint that closes the system of equations.

Before moving on to the solutions to the above equations, we would like to stress the fact that, in all of the formalisms that we introduced to analyze radial perturbations to self-similar solutions, we never invoked or explicitly made reference to the unperturbed shock position or velocity. In particular, Equation (22)—which parameterizes deviations from pure self-similarity—only introduces the dimensionless parameter ζ, which itself is written as a function of τ (which is just the log of the true shock position). The consistency of the perturbation approach is only dependent on the smallness of ζ. We will return to the relevance of this point and its advantages (compared to working with the perturbations to the shock position and velocity individually) when discussing the previous literature in Section 3.5 and when making comparison to simulations in Section 5.

Equation (27) describes the general space- and time-dependent evolution of the perturbations imposed on top of the self-similar profile. We seek solutions to this equation of the form

$$\begin{eqnarray}&&\zeta ={\zeta }_{\sigma }{e}^{\sigma \tau },\end{eqnarray} \tag{ 31 }$$

where ζ_σ is a constant and σ is an undetermined parameter. We can satisfy the differential equations if we let

$$\begin{eqnarray}&&{F}_{1}={\zeta }_{\sigma }{F}_{\sigma }(\xi ){e}^{\sigma \tau },\quad {G}_{1}={\zeta }_{\sigma }{G}_{\sigma }(\xi ){e}^{\sigma \tau },\quad {s}_{1}={\zeta }_{\sigma }{s}_{\sigma }(\xi ){e}^{\sigma \tau },\end{eqnarray} \tag{ 32 }$$

where the functions F_σ, G_σ, and s_σ satisfy the boundary conditions at the shock (Equation (28)). Notice, however, that the amplitude ζ_σ scales out of the problem, because the boundary conditions and the right-hand side of Equation (27) are all proportional to this parameter. The value of σ is then determined solely by the fourth boundary condition, making σ an "eigenvalue" from the standpoint that it is uniquely specified by the continuity of the variables through the sonic point. The eigenvalue equation reads

$$\begin{eqnarray}&&\begin{array}{l}\sigma \left(\begin{array}{c}{s}_{\sigma }\\ {F}_{\sigma }\\ {G}_{\sigma }\end{array}\right)\,+\,\left(\begin{array}{ccc}{f}_{0}-\xi & 0 & 0\\ {\xi }^{2}{h}_{0} & {f}_{0}-\xi & {\gamma }_{2}{\xi }^{2}{h}_{0}\\ 0 & {\xi }^{-2}{g}_{0}^{-1} & {f}_{0}-\xi \end{array}\right)\displaystyle \frac{\partial }{\partial \xi }\left(\begin{array}{c}{s}_{\sigma }\\ {F}_{\sigma }\\ {G}_{\sigma }\end{array}\right)+\left[\begin{array}{ccc}0 & {\xi }^{-2}{g}_{0}^{-1}{s}_{0}^{{\prime} } & 0\\ {\xi }^{2}{h}_{0}^{{\prime} } & \displaystyle \frac{3}{2}-n+2{f}_{0}^{{\prime} } & {\xi }^{2}({\gamma }_{2}-1){h}_{0}^{{\prime} }\\ 0 & 0 & 0\end{array}\right]\,\left(\begin{array}{c}{s}_{\sigma }\\ {F}_{\sigma }\\ {G}_{\sigma }\end{array}\right)\,=\,\left(\begin{array}{c}-2\sigma \\ 2\displaystyle \frac{{g}_{0}}{{V}_{{\rm{c}}}^{2}}-{f}_{0}{g}_{0}{\xi }^{2}\sigma \\ 0\end{array}\right).\end{array}\end{eqnarray} \tag{ 33 }$$

Equation (31) by itself cannot characterize the general deviation of the shock position from the self-similar one, because at any given time we can measure the shock velocity, acceleration, and all higher-order derivatives, and a single value of ζ_σ will not be able to accommodate all of these constraints. However, if there are infinitely many eigenvalues, then we can write

$$\begin{eqnarray}&&\zeta =\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{e}^{\sigma \tau },\end{eqnarray} \tag{ 34 }$$

and we can relate the set {ζ_σ} to the properties of the shock at a given time.⁵ The total solution for F₁ is then

$$\begin{eqnarray}&&{F}_{1}=\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{F}_{\sigma }{e}^{\sigma \tau },\end{eqnarray} \tag{ 35 }$$

and similarly for the other variables. It is straightforward to show that if σ is a solution to the eigenvalue equation, then the complex conjugate σ* is also a solution. Therefore, to ensure that the shock position and post-shock fluid quantities are purely real, we require ζ_σ = (ζ_σ*)*. This can be achieved, with both the shock position and the post-shock fluid variables having no imaginary components, by letting ζ_σ itself be purely real.

3.1. Shock Position and Velocity

The defining equation for ζ, given by Equation (22), can be rearranged to yield the differential equation for the position of the shock as a function of time:

$$\begin{eqnarray}&&{\left(1-2\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{e}^{\sigma \tau }\right)}^{1/2}\displaystyle \frac{d}{{dt}}\left({e}^{\tfrac{3\tau }{2}}\right)=\displaystyle \frac{3}{2}\displaystyle \frac{{V}_{{\rm{c}}}\sqrt{{GM}}}{{r}_{{\rm{a}}}^{3/2}}.\end{eqnarray} \tag{ 36 }$$

We reiterate that the ζ_σ in this expression are determined by the deviations from the shock velocity, acceleration, and all higher-order derivatives at τ = 0. Therefore, once these quantities are specified, this differential equation can simply be integrated numerically to solve for the position of the shock given some initial condition.

While the numerical integration of the above equation for R(t) is straightforward, since ζ(τ) is assumed to be small, we can also write R = R₀ + R₁, V = V₀ + V₁, with R₀ and V₀ the solutions to the zeroth-order terms in Equation (36) and R₁ and ${V}_{1}={{dR}}_{1}/{dt}$ satisfying the first-order corrections induced by finite ζ. Taylor-expanding about ζ and using our solution for ζ in terms of the eigenmodes in Equation (22) then gives

$$\begin{eqnarray}&&{V}_{0}={V}_{{\rm{c}}}\sqrt{\displaystyle \frac{{GM}}{{R}_{0}}},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{V}_{1}}{{V}_{0}}+\displaystyle \frac{1}{2}\displaystyle \frac{{R}_{1}}{{R}_{0}}=\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{e}^{\sigma \tau }.\end{eqnarray} \tag{ 38 }$$

Due to the arbitrariness of the temporal origin, we can, without loss of generality, define the perturbation to the position of the shock to be zero at some τ₀, which we can also define to be zero owing to the scale invariance of the problem (recall that $\tau =\mathrm{ln}(R/{r}_{{\rm{a}}})$). The unperturbed shock position is therefore

$$\begin{eqnarray}&&{R}_{0}(t)={r}_{{\rm{a}}}{\left(1+\displaystyle \frac{3}{2}\displaystyle \frac{{V}_{{\rm{c}}}\sqrt{{GM}}}{{r}_{{\rm{a}}}^{3/2}}t\right)}^{2/3},\end{eqnarray} \tag{ 39 }$$

and the solutions for the perturbations to the shock position and velocity are

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{1}}{{R}_{0}}=\displaystyle \sum _{\sigma }\displaystyle \frac{{\zeta }_{\sigma }}{\tfrac{3}{2}+\sigma }\left({e}^{\sigma \tau }-{e}^{-\tfrac{3}{2}\tau }\right),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{V}_{1}}{{V}_{0}}=\displaystyle \sum _{\sigma }\displaystyle \frac{{\zeta }_{\sigma }}{\tfrac{3}{2}+\sigma }\left[\left(1+\sigma \right){e}^{\sigma \tau }+\displaystyle \frac{1}{2}{e}^{-\tfrac{3}{2}\tau }\right].\end{eqnarray} \tag{ 41 }$$

Note that, in Equation (40), the solution for the perturbation to the shock position appears not as a sum of modes, but as a sum of differences between modes with eigenvalues σ and −3/2. However, the latter is not a true eigenvalue, as it is merely a consequence of the scale invariance of the problem—we can define the perturbation to the shock position to be zero at τ = 0, and the second factor in parentheses in Equation (40) simply enforces this initial condition. The insignificance of the −3/2 "mode" is even more apparent if we simply numerically integrate the differential equation for R(t) (Equation (36)), as in this case there is no direct appearance of this temporal dependence.

3.2. Post-shock Fluid Quantities in Unperturbed Variables

The total solutions for the fluid velocity, density, and pressure behind the shock generated by perturbations to the shock position are

$$\begin{eqnarray}\begin{array}{rcl}v & = & V(t)\left[{f}_{0}(\xi )+{f}_{1}(\xi ,\tau )\right]\,=\,V(t)\left[{f}_{0}(\xi )+\displaystyle \frac{1}{{g}_{0}{\xi }^{2}}\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{e}^{\sigma \tau }{F}_{\sigma }{e}^{\sigma \tau }\right],\end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}\begin{array}{rcl}\rho & = & R{\left(t\right)}^{-n}\left[{g}_{0}(\xi )+{g}_{1}(\xi ,\tau )\right]\,=\,R{\left(t\right)}^{-n}\left[{g}_{0}(\xi )+{g}_{0}(\xi )\displaystyle \sum {\zeta }_{\sigma }{e}^{\sigma \tau }{G}_{\sigma }\right],\end{array}\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}\begin{array}{rcl}p & = & R{\left(t\right)}^{-n}V{\left(t\right)}^{2}\left\{{h}_{0}(\xi )+{h}_{1}(\xi ,\tau )\right\}=R{\left(t\right)}^{-n}V{\left(t\right)}^{2}\,\left[{h}_{0}(\xi )+{h}_{0}(\xi )\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{e}^{\sigma \tau }\left({s}_{\sigma }+\gamma {G}_{\sigma }\right)\right],\end{array}\end{eqnarray} \tag{ 44 }$$

where V(t) is the total shock velocity, ξ = r/R is in terms of the total shock position, and $\tau =\mathrm{ln}(R/{r}_{{\rm{a}}})$ is the log of the true shock position. Again, once ζ_σ is specified, we can use our derived expressions for the total shock position and velocity to completely constrain the deviations from self-similarity using the eigenmodes F_σ and write the result in terms of the physical coordinates r and t.

As we saw in the previous subsection, we can also decompose the shock velocity and shock position into perturbed and unperturbed parts, $V={V}_{0}+{V}_{1}$ and $R={R}_{0}+{R}_{1}$. We can therefore use these first-order corrections to derive expressions for the velocity, density, and pressure profiles in terms of the unperturbed shock velocity, unperturbed shock position, and unperturbed self-similar variable ${\xi }_{0}=r/{R}_{0}(t)$, which are equivalent to using the true shock velocity, true shock position, and total self-similar variable to first order in ζ. The result is

$$\begin{eqnarray}&&v(r,t)={V}_{0}\left[{f}_{0}({\xi }_{0})+\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{f}_{\sigma }({\xi }_{0},\tau )\right],\end{eqnarray} \tag{ 45 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\sigma }({\xi }_{0},\tau ) & = & \left[\displaystyle \frac{\left(1+\sigma \right){f}_{0}-{\xi }_{0}{f}_{0}^{{\prime} }}{\tfrac{3}{2}+\sigma }+\displaystyle \frac{{F}_{\sigma }}{{g}_{0}{\xi }_{0}^{2}}\right]{e}^{\sigma \tau }+\displaystyle \frac{1}{\tfrac{3}{2}+\sigma }\left(\displaystyle \frac{1}{2}{f}_{0}+{f}_{0}^{{\prime} }{\xi }_{0}\right){e}^{-\tfrac{3}{2}\tau }.\end{array}\end{eqnarray} \tag{ 46 }$$

Similarly, the density and pressure are

$$\begin{eqnarray}&&\rho (r,t)={\rho }_{{\rm{a}}}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{-n}\left[{g}_{0}({\xi }_{0})+\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{g}_{\sigma }\right],\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}\begin{array}{rcl}{g}_{\sigma }({\xi }_{0},\tau ) & = & \left(-\displaystyle \frac{{{ng}}_{0}+{\xi }_{0}{g}_{0}^{{\prime} }}{\tfrac{3}{2}+\sigma }+{g}_{0}{G}_{\sigma }\right){e}^{\sigma \tau }+\displaystyle \frac{{{ng}}_{0}+{\xi }_{0}{g}_{0}^{{\prime} }}{\tfrac{3}{2}+\sigma }{e}^{-\tfrac{3}{2}\tau },\end{array}\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&p(r,t)={\rho }_{0}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{-n}{V}_{0}^{2}\left[{h}_{0}({\xi }_{0})+\displaystyle \sum _{\sigma }{\zeta }_{\sigma }{h}_{\sigma }({\xi }_{0},\tau )\right],\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}\begin{array}{rcl}{h}_{\sigma }({\xi }_{0},\tau ) & = & \left[-\displaystyle \frac{{{nh}}_{0}+{\xi }_{0}{h}_{0}^{{\prime} }-2\left(1+\sigma \right){h}_{0}}{\tfrac{3}{2}+\sigma }\right.+\left.\,\gamma {h}_{0}{G}_{\sigma }+{h}_{0}{s}_{\sigma }\right]{e}^{\sigma \tau }+\displaystyle \frac{\left(n+1\right){h}_{0}+\xi {h}_{0}^{{\prime} }}{\tfrac{3}{2}+\sigma }{e}^{-\tfrac{3}{2}\tau }.\end{array}\end{eqnarray} \tag{ 50 }$$

Equations (45)–(50) are the self-consistent, first-order expressions for the post-shock velocity, density, and pressure in terms of the physical radial coordinate r. To lowest order in ζ, these expressions are identical to Equations (42)–(44), which use the true shock position, velocity, and self-similar variable.

3.3. Total Energy and Energy Flux

The energy equation takes the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\mathscr{E}}}{\partial t}+\displaystyle \frac{\partial {{\mathscr{F}}}_{}}{\partial r}=0,\end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray}&&{\mathscr{E}}=\left(\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{p}{\rho }-\displaystyle \frac{{GM}}{r}\right){r}^{2}\rho \end{eqnarray} \tag{ 52 }$$

and

$$\begin{eqnarray}&&{{\mathscr{F}}}_{}=\left(\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{p}{\rho }-\displaystyle \frac{{GM}}{r}\right){r}^{2}\rho v\end{eqnarray} \tag{ 53 }$$

are the energy density and energy flux. Integrating from r = 0 to R +  and taking the limit as $\epsilon \to 0$ gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{sh}}}{\partial t}=4\pi V{{\mathscr{E}}}_{{\rm{a}}}(R)+4\pi {{\mathscr{F}}}_{}(0),\end{eqnarray} \tag{ 54 }$$

where

$$\begin{eqnarray}&&{E}_{\mathrm{sh}}=4\pi {\int }_{0}^{R}{\mathscr{E}}{dr}\end{eqnarray} \tag{ 55 }$$

is the total energy behind the shock and

$$\begin{eqnarray}&&{{\mathscr{E}}}_{{\rm{a}}}(R)={{GMr}}_{{\rm{a}}}{\rho }_{{\rm{a}}}\displaystyle \frac{1-\gamma }{\gamma }{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{1-n}\end{eqnarray} \tag{ 56 }$$

is the energy density of the ambient medium at the location of the shock. After expanding out to first order, this expression becomes

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {E}_{\mathrm{sh}}}{\partial t}=\displaystyle \frac{1-\gamma }{\gamma }4\pi {{GMr}}_{{\rm{a}}}{\rho }_{{\rm{a}}}{V}_{0}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{1-n}\,\left[1+\displaystyle \frac{{V}_{1}}{{V}_{0}}+\left(1-n\right)\displaystyle \frac{{R}_{1}}{{R}_{0}}\right]+4\pi {{\mathscr{F}}}_{0}(0)+4\pi {{\mathscr{F}}}_{1}(0),\end{array}\end{eqnarray} \tag{ 57 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{\mathscr{F}}}_{0} & = & {{GMr}}_{{\rm{a}}}{\rho }_{{\rm{a}}}{V}_{0}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{1-n}{f}_{0}{g}_{0}{\xi }^{2}\,\left(\displaystyle \frac{1}{2}{V}_{{\rm{c}}}^{2}{f}_{0}^{2}+\displaystyle \frac{\gamma }{\gamma -1}{V}_{{\rm{c}}}^{2}\displaystyle \frac{{h}_{0}}{{g}_{0}}-\displaystyle \frac{1}{\xi }\right)\end{array}\end{eqnarray} \tag{ 58 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{{\mathscr{F}}}_{1} & = & {{GMr}}_{{\rm{a}}}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{1-n}{V}_{0}{\xi }^{2}{g}_{0}{f}_{0}\left\{\left(\displaystyle \frac{1}{2}{V}_{{\rm{c}}}^{2}{f}_{0}^{2}+\displaystyle \frac{\gamma }{\gamma -1}{V}_{{\rm{c}}}^{2}\displaystyle \frac{{h}_{0}}{{g}_{0}}-\displaystyle \frac{1}{\xi }\right)\right.\,\left[\displaystyle \frac{{g}_{1}}{{g}_{0}}+\displaystyle \frac{{f}_{1}}{{f}_{0}}+\left(2-n\right)\displaystyle \frac{{R}_{1}}{{R}_{0}}+\displaystyle \frac{{V}_{1}}{{V}_{0}}\right]\\ & & +{V}_{{\rm{c}}}^{2}{f}_{0}{f}_{1}+{V}_{{\rm{c}}}^{2}{f}_{0}^{2}\displaystyle \frac{{V}_{1}}{{V}_{0}}+\left.\,\displaystyle \frac{\gamma }{\gamma -1}{V}_{{\rm{c}}}^{2}\displaystyle \frac{{h}_{0}}{{g}_{0}}\left(2\displaystyle \frac{{V}_{1}}{{V}_{0}}+\displaystyle \frac{{h}_{1}}{{h}_{0}}-\displaystyle \frac{{g}_{1}}{{g}_{0}}\right)+\displaystyle \frac{1}{\xi }\displaystyle \frac{{R}_{1}}{{R}_{0}}\right\}.\end{array}\end{eqnarray} \tag{ 59 }$$

The total correction to the energy therefore has two terms, one arising from the additional sweeping up of binding energy from the ambient medium, and another from the change in the flux of energy onto the black hole.

3.4. Total Mass and Accretion Rate

Multiplying the continuity Equation (1) by r² and integrating from zero to R +  gives, after taking the limit as $\epsilon \to 0$,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {M}_{\mathrm{sh}}}{\partial t}=4\pi V{\rho }_{{\rm{a}}}(R){R}^{2}+\dot{M}(0),\end{eqnarray} \tag{ 60 }$$

where

$$\begin{eqnarray}&&{M}_{\mathrm{sh}}={\int }_{0}^{R}4\pi {r}^{2}\rho {dr}\end{eqnarray} \tag{ 61 }$$

is the total mass contained behind the shock and

$$\begin{eqnarray}&&\dot{M}=4\pi {r}^{2}\rho v\end{eqnarray} \tag{ 62 }$$

is the mass flux. Expanding to first order, we find

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {M}_{\mathrm{sh}}}{\partial t} & = & {\rho }_{{\rm{a}}}{r}_{0}^{2}{V}_{0}{\left(\displaystyle \frac{{R}_{0}}{{r}_{0}}\right)}^{2-n}\left[1+\displaystyle \frac{{V}_{1}}{{V}_{0}}+\left(2-n\right)\displaystyle \frac{{R}_{1}}{{R}_{0}}\right]+{\dot{M}}_{0}(0)+{\dot{M}}_{1}(0),\end{array}\end{eqnarray} \tag{ 63 }$$

where

$$\begin{eqnarray}&&{\dot{M}}_{0}=4\pi {\rho }_{{\rm{a}}}{r}_{{\rm{a}}}^{2}{V}_{0}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{2-n}{f}_{0}{g}_{0}{\xi }^{2}\end{eqnarray} \tag{ 64 }$$

and

$$\begin{eqnarray}&&{\dot{M}}_{1}={\dot{M}}_{0}\left[\displaystyle \frac{{f}_{1}}{{f}_{0}}+\displaystyle \frac{{g}_{1}}{{g}_{0}}+\displaystyle \frac{{V}_{1}}{{V}_{0}}+\left(2-n\right)\displaystyle \frac{{R}_{1}}{{R}_{0}}\right].\end{eqnarray} \tag{ 65 }$$

As was true for the energy, the total mass behind the shock is modified by both the change in the rate at which material is swept up and the change in the removal of material at the origin. This expression also shows that the perturbation to the accretion rate onto the black hole is

$$\begin{eqnarray}&&{\dot{M}}_{\bullet ,1}=-{\dot{M}}_{0}\,\mathop{\mathrm{lim}}\limits_{\xi \to 0}\left[\displaystyle \frac{{f}_{1}}{{f}_{0}}+\displaystyle \frac{{g}_{1}}{{g}_{0}}+\displaystyle \frac{{V}_{1}}{{V}_{0}}+\left(2-n\right)\displaystyle \frac{{R}_{1}}{{R}_{0}}\right].\end{eqnarray} \tag{ 66 }$$

The minus sign appears here because the rate at which matter leaves the post-shock fluid equals the rate at which matter is accreted by the black hole.

3.5. Comparison to Previous Approaches

The previous section set out a general formalism for analyzing the temporal and spatial evolution of perturbations to the self-similar shock position and the fluid behind the shock. The general solutions were written in terms of eigenmodes, with the eigenfrequencies determined by the continuity of the fluid variables through a sonic point. In this section we use this formalism to derive the eigenfrequencies and the eigenfunctions for different values of n, γ₁, and γ₂ to determine the stability of the self-similar solutions.

4.1. A Specific Example: n = 2.5, ${\gamma }_{1}={\gamma }_{2}=1.4$

Here we focus on the case where the ambient density falls off as ρ ∝ r^−2.5, and the pre- and post-shock adiabatic indices are equal to their polytropic indices of γ = 1 + 1/n = 1.4. These parameters were used to compare to weak shock propagation in the hydrogen envelope of a yellow supergiant in Paper I and therefore have physical relevance. The critical shock velocity is V_c ≃ 1.202, and the top left panel of Figure 1 shows the solutions for the self-similar velocity, density, and pressure that pass through the critical point ξ_c ≃ 0.629 and result in accretion at the origin. The asymptotic behavior of these functions near the origin is ${f}_{0}(\xi )\propto {\xi }^{-1/2}$, so the velocity approaches freefall; ${g}_{0}(\xi )\propto {\xi }^{-3/2}$, so the density approaches its time-independent, freefall value; and ${h}_{0}{(\xi )\propto {g}_{0}(\xi )}^{\gamma }$, so the flow is isentropic (the latter condition actually holds over the entirety of the post-shock gas, and the flow is exactly isentropic everywhere; this is not true if ${\gamma }_{2}\ne 1+1/n$). These scalings are depicted in the top right panel of Figure 1.

Zoom In Zoom Out Reset image size

4.1.1. Growing Eigenmode

To determine the eigenfrequencies, we integrate Equation (33) numerically from the shock front (ξ = 1) inward. An arbitrarily chosen σ will result in the divergence of the solutions near the sonic point, so we iteratively perturb the value of the eigenfrequency until the fluid quantities smoothly pass through the critical point. We find that there is one purely real and positive eigenfrequency that both satisfies the boundary conditions at the shock and results in the continuity of the functions through the critical point when n = 2.5 and ${\gamma }_{1}={\gamma }_{2}=1.4$, and its value is

$$\begin{eqnarray}&&{\sigma }_{}\simeq 0.175311.\end{eqnarray} \tag{ 67 }$$

We do not find any other purely real or growing modes (though there are complex and decaying modes; see below). The bottom left panel of Figure 1 shows the eigenmodes for the velocity, f_σ, the density, ${g}_{0}{G}_{\sigma }$, and the pressure, ${h}_{0}({s}_{\sigma }+\gamma {G}_{\sigma })$, for this growing mode, which satisfy Equation (33) and smoothly pass through the critical point at ${\xi }_{{\rm{c}}}\simeq 0.629$. The fact that these functions approach finite numbers near the origin implies that the perturbations to the velocity, density, and pressure scale identically to their unperturbed quantities near the black hole.

Because this eigenvalue is positive and real, the amplitudes of the perturbations grow with time, and at late times this mode dominates over all of the others. Returning to Equations (40) and (41), good approximations for the late-time evolution of the shock position and velocity are therefore

$$\begin{eqnarray}&&R\simeq {R}_{0}\left(1+\zeta \displaystyle \frac{1}{\tfrac{3}{2}+\sigma }{e}^{\sigma \tau }\right),\quad {V}_{1}\simeq {V}_{0}\left(1+\zeta \displaystyle \frac{1+\sigma }{\tfrac{3}{2}+\sigma }{e}^{\sigma \tau }\right).\end{eqnarray} \tag{ 68 }$$

If we further recall that $\tau =\mathrm{ln}({R}_{0}/{r}_{{\rm{a}}})$ and that ${R}_{0}\simeq {t}^{2/3}$, then using σ ≃ 0.175 allows us to rewrite the above as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}R}{{R}_{0}}\propto \displaystyle \frac{{\rm{\Delta }}V}{{V}_{0}}\propto {t}^{0.117},\end{eqnarray} \tag{ 69 }$$

where ΔR and ΔV are the differences between the full and unperturbed shock position and velocity, respectively. These expressions demonstrate that, consistent with the unstable nature of the mode, perturbations to the shock on top of the self-similar solutions grow at late times and cause the solutions to asymptotically deviate from their self-similar prescriptions. However, the rate at which these deviations grow is extremely slow, both because the growth is in the form of a power law, and not an exponential, and because the power-law index itself is very close to zero; it would take, for example, roughly 10 orders of magnitude in time for a perturbation to grow by an additional factor of 10. Therefore, for perturbations that are not so large that nonlinear effects start to become important, the shock properties should be very well reproduced by Equation (68) even for very late times; we validate this notion in Section 5, where we compare these predictions to numerical simulations.

4.1.2. Higher-order Modes

As we noted above, the eigenmode approach is only consistent if there is an infinite number of modes, as this allows us to completely describe the shock position and all higher-order derivatives at a given time. Finding the higher-order modes numerically is challenging, as the divergence of the functions near the singular point becomes very weak over certain ranges of the eigenvalue. To minimize the numerical uncertainty related to the divergence of the functions near the critical point, we first reduce the set of eigenvalue Equations (33) analytically to a single, second-order equation for ${F}_{\sigma }^{{\prime\prime} }$ (the equation is only second order when γ = 1 + 1/n and the unperturbed entropy is exactly constant). We then use L'Hôpital's rule to solve for the relationship between F_σ and ${F}_{\sigma }^{\prime}$ at the critical point, and we integrate inward from the shock front and seek solutions that satisfy this constraint. However, near the critical point we replace the full form of the right-hand side of the equation, which diverges for solutions that do not satisfy the critical point condition, with the reduced right-hand side that we obtain from L'Hôpital's rule under the assumption that the function is continuous through the critical point. Therefore, the equation always integrates continuously through the critical point, but its value at the critical point will only be correct if it satisfies the correct relation deduced from L'Hôpital's rule. In this way we assign a meaningful, finite value of the function at the critical point and minimize the uncertainty related to spurious numerical divergences that arise from maintaining the full right-hand side.

Following this approach, we find a series of eigenmodes that can be ordered by the magnitude of the increasing value of the (absolute value of the) imaginary part of the eigenvalue; doing so, the next eigenvalue that we find is σ ≃ −2.53 ± 0.602i, which decays and oscillates. The next-highest mode we find is σ ≃ −3.86 ± 1.43i; beyond this, the real part of each mode becomes slightly more negative with increasing order, and the imaginary part increases by roughly a factor of one for each higher mode. Defining the eigenfrequency as $\sigma ={\sigma }_{{\rm{r}}}+i{\sigma }_{{\rm{i}}}$, the left panel of Figure 2 shows the real part of the frequency as a function of mode number for the first 25 modes, while the right panel shows the imaginary part (there are also solutions with −σ_i).

Zoom In Zoom Out Reset image size

The eigenfunctions will also have real and imaginary parts; defining ${F}_{\sigma }\equiv {F}_{\sigma ,{\rm{r}}}+{{iF}}_{\sigma ,{\rm{i}}}$, the left (right) panel of Figure 3 shows the real (imaginary) part of the mass flux F_σ for the first five eigenmodes. The fact that these functions all equal a constant near the black hole demonstrates that the perturbations all approach freefall near the origin. Although we have not shown them, the higher-order modes for the relative density G_σ and the entropy s_σ all show similar behavior, each exhibiting somewhat oscillatory behavior and approaching a constant near the origin.

Zoom In Zoom Out Reset image size

Unlike standard eigenmodes, each higher-order mode for the mass flux does not possess an additional zero crossing. The reason for this is that these modes do not form an orthogonal basis; equivalently, the operator equation for ${F}_{\sigma }^{{\prime\prime} }$, while second order, is not in Sturm–Liouville form.

Because these higher-order modes are complex, one expects the shock position to exhibit oscillations in time. Returning to our general expression for the shock position (Equation (40)) and expanding σ in terms of its real and imaginary parts, we find

$$\begin{eqnarray}\begin{array}{rcl}R & = & {R}_{0}\left\{1+\displaystyle \sum _{\sigma }\displaystyle \frac{{\zeta }_{\sigma }}{{\left(\tfrac{3}{2}+{\sigma }_{{\rm{r}}}\right)}^{2}+{\sigma }_{{\rm{i}}}^{2}}\,\left[{e}^{{\sigma }_{{\rm{r}}}\tau }\left(\left(\displaystyle \frac{3}{2}+{\sigma }_{{\rm{r}}}\right)\cos \left({\sigma }_{{\rm{i}}}\tau \right)+{\sigma }_{{\rm{i}}}\sin \left({\sigma }_{{\rm{i}}}\tau \right)\right)\right.-\left.\left(\displaystyle \frac{3}{2}+{\sigma }_{{\rm{r}}}\right){e}^{-\tfrac{3}{2}\tau }\right]\right\},\end{array}\end{eqnarray} \tag{ 70 }$$

where we let ${\zeta }_{\sigma }\to {\zeta }_{\sigma }/2$ to maintain consistency with the case where ${\sigma }_{{\rm{i}}}=0$. In line with our expectations, a nonzero imaginary part of the eigenfrequency introduces temporal oscillations into the position of the shock. However, unlike the case where the background flow is static, these oscillations proceed logarithmically in time (i.e., $\tau =\mathrm{ln}{R}_{0}$, so $\cos \left({\sigma }_{{\rm{i}}}\tau \right)\simeq \cos \left({\sigma }_{{\rm{i}}}\mathrm{ln}t\right)$, and their periods grow exponentially).

Figure 4 shows the relative change in the shock position (left) and the shock velocity (right) as we increase the relative importance of the second-order mode. As we increase the ratio ζ₂/ζ₁ (we set ζ₁ = 0.1 for definiteness, but this only changes the scale of the y-axis in the perturbative limit), we see that oscillations start to become more pronounced for both the shock position and the velocity. Because the real component of the second-order mode is large, the oscillations are damped fairly quickly, and only the growing mode dominates at times t/T_a ≳ 10.

Zoom In Zoom Out Reset image size

4.2. Growing Modes for Other $n,\gamma$

The previous section demonstrated that, when n = 2.5 and γ₁ = γ₂ = 1 + 1/n = 1.4, the self-similar solutions are weakly unstable, with perturbations growing as ∝ t^0.117 at late times. Here we investigate the nature of the instability for other n and γ₁ = γ₂ = γ, and we derive the power-law index for the growing modes as a function of these parameters. We focus exclusively on the case when γ₁ = γ₂ because, when this condition is violated, the solutions behave qualitatively differently and, as can be seen from the boundary conditions at the shock front (Equations (6)–(8)), the trivial solution when the shock Mach number equals 1 stops existing (since γ₁ $\ne$ γ₂ is inconsistent with a ${\mathscr{M}}=1$ solution in which nothing changes at the "shock").

Figure 5 shows the unstable eigenvalue σ, multiplied by 2/3 to show the asymptotic temporal scaling of the perturbations (recall that the perturbations were parameterized as $\propto {R}^{\sigma }\propto {t}^{2\sigma /3}$), as a function of n for the γ shown in the legend. This plot demonstrates that the growth rate of the perturbations exhibits some dependence on the power-law index of the ambient medium, equaling zero when n = 2 and when n = 3.5; as noted in Paper I, these values of n are special because they are identical to the Sedov–Taylor blast wave (n = 2) and the shock reduces to an RW (i.e., the shock moves at the local sound speed and is an infinitesimal perturbation on the ambient gas; n = 3.5). In between these two extremes, the growth rate of the perturbations reaches a relative maximum of 2σ/3 ≲ 0.12 for n ≃ 2.5.

Zoom In Zoom Out Reset image size

Interestingly, we see that the value of the growth rate is almost completely independent of the adiabatic index, with the most noticeable difference being a slight trend with the peak growth rate and the value of n at which that peak is obtained. Specifically, smaller γ result in a maximum growth rate at a smaller value of n. However, the differences between the curves in this figure are extremely small. This weak dependence on the adiabatic index is not too surprising, as the critical velocity V_c is not heavily modified by γ, and the functions f₀, g₀, and h₀ are also not too affected by this value. The black dashed curve in this figure is the fit

$$\begin{eqnarray}&&\displaystyle \frac{2\sigma }{3}=1.58{\left(n-2\right)}^{0.87}-1.5\left(n-2\right),\end{eqnarray} \tag{ 71 }$$

which, aside from slightly overestimating the growth rate for small and large n, does a very good job at reproducing the numerically obtained curves. Table 2 gives the specific value of the growth rate for a range of n and different values of γ.

Table 2.  The Left Column gives the Power-law Index of the Ambient Medium, and Each Cell Gives the Critical Shock Speed V_c of the Unperturbed Self-similar Solution and the Growth Rate of the Unstable Mode (the Temporal Scaling of the Perturbations Is ${t}^{2\sigma /3}$)

n ($\rho \propto {r}^{-n}$)	${\gamma }_{1}={\gamma }_{2}=1+1/n$	${\gamma }_{1}={\gamma }_{2}=4/3$	γ1 = γ2 = 3/2	γ1 = γ2 = 5/3
n = 2.01	$\{{V}_{{\rm{c}}},\sigma \}=\{8.72,0.0104\}$	{3.92, 0.0246}	{8.78, 0.0103}	{11.5, 0.00967}
2.1	{2.93, 0.0828}	{2.29, 0.103}	{3.03, 0.0813}	{3.59, 0.0759}
2.2	{2.05, 0.131}	{1.76, 0.143}	{2.15, 0.128}	{2.48, 0.0122}
2.3	{1.64, 0.158}	{1.47, 0.165}	{1.74, 0.156}	{1.97,0.150}
2.4	{1.38, 0.172}	{1.28, 0.175}	{1.48, 0.170}	{1.66,0.166}
2.5	{1.20, 0.175}	{1.14, 0.177}	{1.30, 0.174}	{1.44,0.173}
2.6	{1.07, 0.172}	{1.02, 0.172}	{1.16, 0.172}	{1.28, 0.172}
2.7	{0.960, 0.163}	{0.934, 0.162}	{1.04, 0.164}	{1.14, 0.165}
2.8	{0.874, 0.149}	{0.860, 0.149}	{0.955, 0.151}	{1.04, 0.153}
2.9	{0.801, 0.133}	{0.796, 0.133}	{0.877, 0.135}	{0.953, 0.137}
3.0	{0.740, 0.114}	{0.740, 0.114}	{0.810, 0.116}	{0.876, 0.119}
3.1	{0.688, 0.0938}	{0.692, 0.0939}	{0.752, 0.0955}	{0.809, 0.0973}
3.2	{0.642, 0.0718}	{0.649, .0720}	{0.701, 0.0730}	{0.780, 0.0742}
3.3	{0.602, 0.0487}	{0.610, 0.0488}	{0.655, 0.0493}	{0.698, 0.0500}
3.4	{0.566, 0.0247}	{0.576, 0.0247}	{0.614, 0.0249}	{0.651, 0.0251}
3.49	{0.538, 0.00249}	{0.547, 0.00248}	{0.581, 0.00248}	{0.613, 0.00249}

Note. Different columns are for different values of the adiabatic index of the gas ${\gamma }_{1}={\gamma }_{2}=\gamma$. Near n = 2, V_c grows to asymptotically large values, and the growth rate approaches zero. Near n = 3.5, V_c approaches the sound speed of $\sqrt{\gamma /(n+1)}$, and the growth rate again nears zero.

Download table as:  ASCIITypeset image

The previous subsections outlined a general framework for analyzing the perturbations to the self-similar, weak shock solutions presented in Paper I. We applied this formalism to the specific case when the density of the ambient medium declines as ρ ∝ r^−2.5 and γ₁ = γ₂ = 1.4 and demonstrated that the lowest-order eigenmode—which characterizes the long-term evolution of perturbations on top of the self-similar solutions—grows as a very weak power law in time. A more general analysis showed that this weak instability should be present for other n and γ as well. In this section, we investigate the validity of the perturbation theory by comparing to the results of idealized, hydrodynamical simulations.

5.1. Simulation Setup

We used the hydrodynamics code flash (v 4.3; Fryxell et al. 2000) to analyze the propagation of a weak shock into an ambient medium that is in hydrostatic equilibrium with a central mass, is non-self-gravitating, and has a density profile of ρ ∝ r⁻ⁿ and adiabatic index γ₁ = 1 + 1/n. The adiabatic index of the post-shock gas was fixed to γ₂ = 1 + 1/n, and hence the setup is identical to the scenario in which we expect the CQR, self-similar solutions to arise.

We generate an outward-propagating shock by initializing a pressure wave on top of the hydrostatic fluid. The pressure is set to $p(t=0)=(1+\delta p)\times {p}_{1}({r}_{{\rm{a}}})$, between the inner boundary, r_a, and wave front ${r}_{{\rm{f}}}=1.5\times {r}_{{\rm{a}}}$. The density and velocity are unchanged. Outflow boundary conditions are used at the inner boundary, and ghost cells are populated by linearly extrapolating from the innermost grid cells on the domain. This construction enables material below r_f to begin falling out of the domain, due to the lack of pressure support, while an outward-propagating shock forms at the overpressurized wave front.

The simulation domain is on a uniform grid that spans between r_a = 1 and r_max = 10³; therefore, the initial pressure wave is spatially small. The simulations present here have a minimum grid resolution of δr ≃ 10⁻² except for the case where δp = 0.1, which has δr ≃ 0.25 × 10⁻². Table 3 gives the parameters for each simulation analyzed here; the first column gives the initial pressure perturbation δp, the second column gives the difference between the shock Mach number and 1, $\delta {\mathscr{M}}$, when the shock reaches a distance of 100 r_a, and the third column is the resolution $\delta r/{r}_{{\rm{a}}}$ (grid spacing) used in the simulation. In the following, we label the simulation by the second column in this table, $\delta {\mathscr{M}}$. In a forthcoming paper (Paper III) we will present more results of these and other simulations.

Table 3.  The Parameters for Each Simulation Done with flash

δp	$\delta {\mathscr{M}}$	δr/ra
0.1	0.78	6.25 × 10−4
0.25	0.99	0.01
0.275	1.0	0.01
0.6	1.7	0.01
1.0	2.6	2.5 × 10−3

Note. The initial pressure perturbation δp, the difference between the Mach number and 1 when the shock reaches a distance of 100r_a, δM, and the grid spacing adopted for each simulation, δr/r_a.

Download table as:  ASCIITypeset image

5.2. Results

The left panel of Figure 6 shows ζ as a function of time in units of ${T}_{{\rm{a}}}={r}_{{\rm{a}}}^{3/2}/\sqrt{{GM}}$ measured from the simulation with $\delta {\mathscr{M}}=0.78$ (orange, solid curve) and the predicted scaling from the growing mode (black dashed curve). The orange curve is obtained by using the numerically obtained shock position and velocity, R and V, to construct

$$\begin{eqnarray}&&\zeta =\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{{{GMV}}_{{\rm{c}}}^{2}}{{{RV}}^{2}}\right).\end{eqnarray} \tag{ 72 }$$

The dashed curve in this figure is given by

$$\begin{eqnarray}&&\zeta ={\zeta }_{\sigma }{e}^{\sigma \mathrm{ln}R},\end{eqnarray} \tag{ 73 }$$

where σ ≃ 0.175 is the growing eigenvalue; ζ_σ is chosen to match the asymptotic scaling approached by the numerical solution, and in this case it is given by ζ_σ = −0.0301. The negative sign of ζ_σ arises here because the initial perturbation yields a shock with a Mach number that falls below the CQR value. The right panel of this figure gives the numerically obtained shock position (orange solid), the self-similar shock position (black dotted), and the perturbed solution obtained by including the growing mode (black dashed). For the latter case, which agrees extremely well with the numerical solution, we used ζ_σ = −0.0301 in Equation (36) and integrated the equation numerically to obtain the shock position.

Zoom In Zoom Out Reset image size

In addition to the shock position, the perturbation analysis of Section 3 predicts that there should also be corrections to the velocity, density, and pressure of the fluid behind the shock. Moreover, these corrections are determined by the amplitude of the perturbation to the growing mode, ζ, meaning that there is one parameter that should account for the asymptotic deviation of the fluid velocity, density, and pressure from their self-similar forms. The left panel of Figure 7 gives the numerically obtained velocity profile as a function of r/r_a (solid, lightly colored curves), the self-similar solution (dotted curves), and the perturbed solutions with the one-parameter correction (dashed curves) at the times shown in the legend for $\delta {\mathscr{M}}=0.78$. The self-similar solution is given by

$$\begin{eqnarray}&&{v}_{\mathrm{ss}}=V(t){f}_{0}(\xi ),\end{eqnarray} \tag{ 74 }$$

where V(t) is the numerically obtained shock velocity and $\xi =r/R(t)$ is the total self-similar variable, while the perturbed velocity profile is

$$\begin{eqnarray}&&{v}_{\mathrm{pert}}=V(t)\left[{f}_{0}(\xi )+{\zeta }_{\sigma }{e}^{\sigma \tau }{f}_{1}(\xi )\right],\end{eqnarray} \tag{ 75 }$$

where ζ_σ = −0.0301, σ = 0.175311, and $\tau =\mathrm{ln}(R(t)/{r}_{{\rm{a}}})$. From this figure, we see that there are small differences between the numerical and self-similar velocity profiles, and these differences are accounted for well by the contribution from the unstable mode.

Zoom In Zoom Out Reset image size

The right panel of Figure 7 shows the difference between the numerical and self-similar velocity profiles normalized by the shock velocity (solid, lightly colored curves) and the prediction from the unstable mode (dashed curves). From this figure, we see that the discrepancy between the self-similar solution and the numerical one is very accurately reproduced by the growing mode. However, there are small differences that become more pronounced at later times, and these differences likely arise from nonlinear terms that are not accounted for in the linear perturbation analysis.

The left panel of Figure 8 shows the numerically obtained density profile (solid, lightly colored), the self-similar density profile (dotted), and the one-parameter correction at the times shown in the legend, and the left panel of Figure 9 gives the analogous curves for the post-shock pressure; the definitions of the self-similar functions are

$$\begin{eqnarray}&&{\rho }_{\mathrm{ss}}={R}^{-n}{g}_{0}(\xi ),\quad {p}_{\mathrm{ss}}={V}^{2}{R}^{-n}{h}_{0}(\xi ),\end{eqnarray} \tag{ 76 }$$

while the perturbed functions are

$$\begin{eqnarray}\begin{array}{rclcl}{\rho }_{\mathrm{pert}} & = & {R}^{-n}\left[{g}_{0}(\xi )+{\zeta }_{\sigma }{e}^{\sigma \tau }{g}_{1}(\xi )\right],\,{p}_{\mathrm{pert}} & = & {V}^{2}{R}^{-n}\left[{h}_{0}(\xi )+{\zeta }_{\sigma }{e}^{\sigma \tau }{h}_{1}(\xi )\right].\end{array}\end{eqnarray} \tag{ 77 }$$

For these figures ζ_σ = −0.0301 is set by the perturbations to the shock position, and σ ≃ 0.175 is the eigenvalue of the growing mode. As was true for the velocity profile, we see that there are small but noticeable differences between the predictions of the self-similar solutions and those obtained from the simulations. The right panels of Figures 8 and 9 show the differences between the numerical solutions and the self-similar solutions (lightly colored, solid) and the prediction from the growing mode (dashed) and demonstrate that the differences are accurately reproduced by the eigenfunctions of the growing mode.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left panel of Figure 10 illustrates ζ as a function of time obtained numerically (solid curves) and the one-parameter correction (dashed curves) for the range of $\delta {\mathscr{M}}$ shown in the legend. In this case, the constant ζ_σ that determines ζ(0) increases as $\delta {\mathscr{M}}$ increases. In particular, denoting the amplitude of the perturbation particular to simulation $\delta {\mathscr{M}}$ by ${\zeta }_{\delta {\mathscr{M}}}$, for this plot we set ζ_0.78 = −0.0301, ζ_0.99 = 0.0206, ζ_1.0 = 0.0291, ζ_1.7 = 0.117, and ζ_2.6 = 0.176. These values provide reasonable fits to the numerical position of the shock. At late times, we see that simulations with relatively large $\delta {\mathscr{M}}$ produce ζ that differ noticeably from the prediction of the perturbation analysis, and the deviations start to become appreciable at a time where ζ ≳ few × 0.1. This breakdown is expected, as in these cases the nonlinear terms are no longer ignorable.

Zoom In Zoom Out Reset image size

The right panel of Figure 10 shows the numerical shock position (solid), the perturbed shock position (dashed), and the self-similar solution (dotted), where the perturbed solution is found by integrating Equation (36) with the growing mode included. As expected from the left panel of this figure, we see that the numerical solutions match the perturbed solutions well when $\delta {\mathscr{M}}$ is small, but there are significant deviations from the predictions of the perturbation theory when $\delta {\mathscr{M}}$ becomes large.

While the amplitude of the growing mode, ζ_σ, is in principle determined from the eigenmode expansion of all of the derivatives of the shock position at a given time, we also expect ζ_σ to scale monotonically with the difference between the Mach number of the shock and the one enforced by the CQR solution. Taking the Mach number at a given instant in time or an average over some temporal range also should not greatly modify the scaling, as the Mach number only depends very weakly on time (and varies as ${t}^{2\sigma /3}$) when the solution is close to the CQR one. Performing a quadratic fit to $\zeta ({\rm{\Delta }}{\mathscr{M}})$, where ${\rm{\Delta }}{\mathscr{M}}$ is the difference between the shock Mach number and the CQR value when the shock reaches R = 100, we find

$$\begin{eqnarray}&&\zeta =0.19{\rm{\Delta }}M-0.053{\rm{\Delta }}{M}^{2}.\end{eqnarray} \tag{ 78 }$$

There is only a very small difference (at the 1% level) if we fit the average ${\rm{\Delta }}{\mathscr{M}}$ over the duration of the simulation.

In the previous section we demonstrated that, when the shock Mach number is larger than the one imposed by the CQR solution, the instability causes the shock position and velocity to asymptotically grow. In this case the shock position and the post-shock quantities should asymptotically approach the energy-conserving, Sedov–Taylor blast wave (or, for sufficiently large n, the Waxman–Shvarts solution; Waxman & Shvarts 1993), and we plan to investigate this transition in detail in a follow-up paper.

However, we also showed that when the shock Mach number is below the CQR prediction, then the value of ζ multiplying the growing mode is negative. Therefore, in this regime the difference between the shock Mach number and the CQR value becomes increasingly negative and, if the linear theory remained valid into the nonlinear regime, would predict that the shock becomes subsonic at late times. However, the Mach number cannot become less than unity, resulting in the dissolution of the shock and the formation of a pressure wave of finite thickness, as the conservation of wave luminosity implies that such a wave should steepen back into a shock (Dewar 1970; Ro & Matzner 2017; Coughlin et al. 2018a). This finding then raises the question as to the asymptotic state of these solutions as they become ever weaker.

One can show that there is, in fact, a second, exact solution to the self-similar Equations (19)–(21)) that have ${V}_{{\rm{c}}}=\sqrt{\gamma /(n+1)}$, or a Mach number of exactly 1. In this case, the critical point is at the location of the shock itself, and one of the solutions is the trivial solution ${f}_{0}=0$, ${g}_{0}={\xi }^{-n}$, ${h}_{0}={\xi }^{-n-1}/\gamma$ (i.e., the ambient medium retains perfect hydrostatic balance everywhere). However, owing to the nature of the critical point, there is a second solution that satisfies f₀(1) = 0, g₀(1) = 1, h₀(1) = 1/γ, meaning that the "post-shock" fluid is still in hydrostatic balance precisely at the location of the sound-crossing horizon, but with derivatives that differ from those of the trivial solution.

These solutions adopt the same self-similar variable as the CQR blast wave, and therefore the self-similar equations are identical:

$$\begin{eqnarray}&&-{{ng}}_{{\rm{r}}}-\xi {g}_{{\rm{r}}}^{\prime} +\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{d}{d\xi }\left({\xi }^{2}{g}_{{\rm{r}}}{f}_{{\rm{r}}}\right)=0,\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&-\displaystyle \frac{1}{2}{f}_{{\rm{r}}}-\xi {f}_{{\rm{r}}}^{\prime} +{f}_{{\rm{r}}}{f}_{{\rm{r}}}^{\prime} +\displaystyle \frac{1}{{g}_{{\rm{r}}}}{h}_{{\rm{r}}}^{\prime} =-\displaystyle \frac{1}{{V}_{{\rm{r}}}^{2}{\xi }^{2}},\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&n\gamma -n-1+\left({f}_{{\rm{r}}}-\xi \right)\displaystyle \frac{d}{{dr}}\mathrm{ln}\left(\displaystyle \frac{{h}_{{\rm{r}}}}{{g}_{{\rm{r}}}^{\gamma }}\right)=0.\end{eqnarray} \tag{ 81 }$$

Here we are letting γ₁ = γ₂ = γ, which is necessary for recovering the trivial solution (i.e., the trivial solution mandates that the wave does nothing, and setting γ₁ $\ne$ γ₂ by definition implies that the wave has changed something). The definitions of the self-similar functions f_r, g_r, and h_r in terms of the velocity, density, and pressure are also the same,

$$\begin{eqnarray}\begin{array}{rclcl}v & = & V(t){f}_{{\rm{r}}}(\xi ),\quad \rho ={\rho }_{{\rm{a}}}{\left[\displaystyle \frac{R(t)}{{r}_{{\rm{a}}}}\right]}^{-n}{g}_{{\rm{r}}}(\xi ),\,p & = & {\rho }_{{\rm{a}}}V{\left(t\right)}^{2}{\left[\displaystyle \frac{R(t)}{{r}_{{\rm{a}}}}\right]}^{-n}{h}_{{\rm{r}}}(\xi ),\end{array}\end{eqnarray} \tag{ 82 }$$

and the shock velocity and position are likewise related in the same manner,

$$\begin{eqnarray}&&{V}^{2}={V}_{{\rm{r}}}^{2}\displaystyle \frac{{GM}}{R},\end{eqnarray} \tag{ 83 }$$

where ${V}_{{\rm{r}}}=\sqrt{\gamma /(n+1)}$, i.e., the "shock" moves out at the sound speed. The jump conditions at the shock front, Equations (6)–(8), also give the boundary conditions obeyed by the functions f_r, g_r, and h_r, and when ${V}_{{\rm{r}}}=\sqrt{\gamma /(n+1)}$ they become

$$\begin{eqnarray}&&{f}_{{\rm{r}}}(1)=0,\quad {g}_{{\rm{r}}}(1)=1,\quad {h}_{{\rm{r}}}=1/\gamma .\end{eqnarray} \tag{ 84 }$$

Using L'Hôpital's rule in Equations (79)–(81) shows that the derivatives of the nontrivial solutions are

$$\begin{eqnarray}\begin{array}{rclcl}{f}_{{\rm{r}}}^{{\prime} }(1) & = & \displaystyle \frac{5+2n}{2+2\gamma },\ {g}_{{\rm{r}}}^{{\prime} }(1)=\displaystyle \frac{5-2n\gamma }{2+2\gamma },\,{h}_{{\rm{r}}}^{{\prime} }(1) & = & -\displaystyle \frac{2+2n-3\gamma }{2\gamma +2{\gamma }^{2}}.\end{array}\end{eqnarray} \tag{ 85 }$$

These RW solutions are thus in hydrostatic equilibrium at radii exterior to the sound-crossing horizon and undergo collapse onto the point mass at the origin for radii interior to that radius. These solutions were also found by Kazhdan & Murzina (1994), who focused on the case of n ≤ 2.

The left panel of Figure 11 shows the solution for the self-similar velocity, multiplied by ξ^1/2, for the range of n shown in the legend, and for this figure we let γ = 1 + 1/n. We see that, near the position of the RW, the velocity is zero and the gas is in hydrostatic equilibrium. Immediately inside of ξ = 1, however, the velocity adopts a nontrivial value, signifying the fact that the sound wave conveys to the overlying envelope that there is accretion at the origin. The fact that these solutions all approach ∼ξ^−1/2 near the origin confirms the fact that they do, indeed, accrete. These solutions also exist for n < 2 and n > 3.5, unlike the CQR blast wave.

Zoom In Zoom Out Reset image size

The right panel of Figure 11 gives the velocity, normalized by the circular velocity of the point mass, as a function of radius normalized by the radius to which the RW has propagated. The solid curves give the RW solutions for a subset of n in the left panel, while the dashed curves are the CQR solutions for the same n. It is apparent that while the functions f_r differ fairly substantially, the velocity profile resulting from the different RWs are all quite similar; this finding implies that the functions f_r are almost the same function aside from a scaling factor. We also see that the RW solutions and the CQR solutions predict different behavior near the location of the shock front. However, all of the curves equal $v/\sqrt{{GM}/r}=-\sqrt{2}$ at the origin, which demonstrates that the velocity always approaches freefall onto the black hole independent of the solution near the shock front.

The left panel of Figure 12 gives the normalized density as a function of dimensionless radius for the RW (solid) and CQR (solutions) and the set of n shown in the legend. In all cases the density approaches the freefall scaling of ρ ∝ r^−3/2 near the origin, but the proportionality factor changes depending on the type of solution (CQR or RW) and n. We also see that, as n nears 2, the density of the CQR solution falls below that of the RW solution. The origin of this behavior is that, as $n\to 2$, the CQR solution converges to the Sedov–Taylor blast wave, where the kinetic energy is infinite relative to the gravitational binding energy. In this limit, then, the point mass effectively vanishes from the solution, and there is no finite radius at which the solution approaches freefall. The black dot-dashed and red dashed curves in this figure show, respectively, the Sedov–Taylor solution for n = 2 and γ = 1.5 and the CQR solution when n = 2.001 (and γ = 1 + 1/2.001), which illustrates directly the notion that the freefalling nature of the solutions drops out as $n\to 2$.

Zoom In Zoom Out Reset image size

The right panel of Figure 12 gives the Bernoulli parameter of the fluid, given by

$$\begin{eqnarray}&&{\mathscr{B}}=\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{p}{\rho }-\displaystyle \frac{{GM}}{r},\end{eqnarray} \tag{ 86 }$$

normalized by the square of the shock velocity. We see that, for the CQR solution, the Bernoulli parameter is characterized by positive and negative regions, the former being confined to the fluid near the shock front, the latter being near the point mass. The RW solutions, on the other hand, have a negative Bernoulli parameter everywhere and only approach exactly zero near the location of the shock front.

Because the RW and CQR solutions adopt the same definitions for the self-similar velocity, density, and pressure and the equations describing the evolution of the fluid quantities are the same, we can employ exactly the same formalism as outlined in Sections 2 and 3 to assess the stability of the RW solutions. In particular, the matrix equation for the eigenmodes describing perturbations to the RW is identical to Equation (33) if we define

$$\begin{eqnarray}&&\begin{array}{r}v(r,t)\,=\,V\left[{f}_{0}(\xi )+{e}^{\sigma \tau }{f}_{\sigma }(\xi )\right],\,\rho (r,t)\,=\,{\rho }_{{\rm{a}}}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}\left[{g}_{0}(\xi )+{e}^{\sigma \tau }{g}_{\sigma }(\xi )\right],\,p\,=\,{V}^{2}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}\left[{h}_{0}+{e}^{\sigma \tau }{h}_{\sigma }(\xi )\right],\,\displaystyle \frac{{GM}}{{{RV}}^{2}}\,=\,\displaystyle \frac{1}{{V}_{{\rm{c}}}^{2}}\left(1-2{\zeta }_{\sigma }{e}^{\sigma \tau }\right),\end{array}\end{eqnarray} \tag{ 87 }$$

with the functions f₀, g₀, and h₀ given by the RW solutions and ${V}_{{\rm{c}}}=\sqrt{\gamma /(n+1)}$. The boundary conditions on the functions f_σ, g_σ, and h_σ are given by

$$\begin{eqnarray}&&{f}_{\sigma }(1)={g}_{\sigma }(1)=\displaystyle \frac{4}{\gamma +1}{\zeta }_{\sigma },\quad {h}_{\sigma }(1)=\displaystyle \frac{2(\gamma -1)}{\gamma (\gamma +1)}{\zeta }_{\sigma },\end{eqnarray} \tag{ 88 }$$

which result from the shock jump conditions.

In this case, the fourth boundary condition that determines the eigenvalue—being that the solutions be continuous through the sonic point—also takes place at the location of the RW. Therefore, at the sonic point we know the values of the perturbations, the values of the unperturbed functions, and, from Equation (85), the values of the derivatives of the unperturbed functions, which implies that we can determine the eigenvalue analytically. Doing so yields

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{1}{4}\left(n-\displaystyle \frac{7}{2}\right).\end{eqnarray} \tag{ 89 }$$

When n < 3.5, the eigenvalue is negative, which implies that the solutions are stable; when n = 3.5, the solutions are marginally stable, which agrees with the predictions of the CQR solution for this value of n (see Figure 5); and when n > 3.5, the solutions are unstable, meaning that any small perturbation to the amplitude of the RW will cause the Mach number to grow. This expression is also completely independent of the adiabatic index of the gas.

The results of Section 4 demonstrated that the lowest-order eigenmode of the CQR blast wave, which characterizes the asymptotic deviation of the shock position and the post-shock velocity, density, and pressure from their self-similar expressions, is weakly unstable, meaning that perturbations imposed on top of the self-similar solution grow as very shallow power laws in time. The numerical simulations performed with flash also validated this finding, and differences between the numerically obtained shock position and the self-similar solution were shown to grow approximately as ${\rm{\Delta }}R/{R}_{0}\propto {t}^{0.117}$ when the density profile of the ambient medium falls off as ρ ∝ r^−2.5 (and the pre- and post-shock adiabatic indices are equal to γ₁ = γ₂ = 1 + 1/2.5 = 1.4). The post-shock fluid quantities were also very well matched by the solutions resulting from the perturbation analysis.

Paper I made the following heuristic argument to suggest that the CQR self-similar solutions are likely stable: the Mach number of the CQR blast wave is uniquely specified by the self-similar solutions and is on the order of unity when the power-law index of the density profile of the ambient medium is not too close to 2. If the initial Mach number is greater than the one imposed by the self-similar solutions, then the energy of the fluid immediately behind the shock is slightly greater than the self-similar value (by virtue of the jump conditions). However, the rate at which binding energy is added to the fluid—being $\dot{E}={E}_{{\rm{a}}}V$, where E_a is the binding energy of the ambient medium—is slightly larger. Therefore, it seemed plausible that the shock would self-regulate by sweeping up binding energy at a greater rate, thereby reducing the energy budget of the post-shock fluid to the required, self-similar value. An analogous argument holds when the Mach number is slightly below the self-similar one.

With the solutions for the perturbations to the CQR blast wave, we can now directly assess the validity of this argument: from the terms in parentheses on the right-hand side of Equation (57), we see that there are two contributions to the perturbation to the binding energy swept up by the shock. The first is proportional to V₁/V₀ and corresponds to the term that, in Paper I, was argued to contribute to the self-regulation of the shock—increases in the velocity correspond to a greater influx of the binding energy of the ambient medium. The second term, proportional to $(1-n){R}_{1}/{R}_{0}$, was not considered in the heuristic assessment of Paper I; this perturbation arises from the fact that the binding energy of the ambient medium has a radial dependence and falls off as ${E}_{{\rm{a}}}\propto {GM}/r\times \rho {r}^{2}\propto {r}^{1-n}$. Thus, if we perturb the position of the shock to larger values, then the material at the new location of the shock is less bound than it was at the original location, which reduces the rate at which binding energy is entrained by the flow. This term, therefore, has a destabilizing effect on the propagation of the shock.

The ability of the shock to self-regulate is therefore determined by the ratio of the second term, $(n-1){R}_{1}/{R}_{0}$, to the first term, ${V}_{1}/{V}_{0}$. In situations where the magnitude of this ratio is greater than (less than) 1, we expect that the shock is unstable (stable) to perturbations to the energy flux. For a given mode, this ratio is given by (from Equations (40) and (41))

$$\begin{eqnarray}&&\displaystyle \frac{(n-1){R}_{1}{V}_{0}}{{V}_{1}{R}_{0}}=\displaystyle \frac{n-1}{1+\sigma }.\end{eqnarray} \tag{ 90 }$$

We thus see that, for small values of σ, the magnitude of this ratio can be greater than 1, and for all of the unstable modes we find that this number is, indeed, slightly greater than unity (e.g., for n = 2.5, the above ratio is ≃1.28). Ironically, this situation arises because the instability is so weak, meaning that the perturbation to the velocity is too small to significantly reduce the post-shock energy budget. Indeed, for larger values of σ, which induce correspondingly larger variations in the velocity, the magnitude of this ratio is less than 1 (e.g., for the second mode for n = 2.5, which has σ ≃ −2.5 + 0.6i, the modulus of this ratio is ≃0.92).

For the unstable solutions, when ζ is positive, the total energy entrained by the shock as it moves out increases asymptotically compared to the self-similar solution, which can be shown directly by integrating the differential equation for $\dot{E}$ (Equation (57)). In this case the shock position and the post-shock quantities should asymptotically approach the energy-conserving, Sedov–Taylor blast wave. In the other limit, when ζ < 1, the relative energy decreases asymptotically to ever more negative values, and in this regime the shock likely transitions to the RW solutions presented in Section 6. We plan to investigate these transitions in detail in a follow-up paper.

For the range of n that satisfy 2 < n < 3.5, there are three exact, self-similar solutions that describe the propagation of shocks in the vicinity of a compact object, and only one—the CQR solution—has a Mach number that is both greater than 1 and less than infinity. Figure 13 illustrates this situation schematically. The green curve in this figure uses the solution for the CQR Mach number when γ₁ = γ₂ = 1 + 1/n (though, as shown by Figure 5, there is not much dependence on the adiabatic index). There is a tempting analogy to draw between the existence of these three solutions and the three types of orbits that are possible in a Keplerian potential, being unbound (positive energy), bound (negative energy), and marginally bound (zero energy). The last type of orbit is unstable, as infinitesimal perturbations to the energy generate either bound or unbound orbits. Furthermore, the Sedov–Taylor blast wave has positive energy and a positive Bernoulli parameter everywhere, the RW has a negative energy and a negative Bernoulli parameter everywhere, and the energy of the CQR solution approaches exactly zero in the limit of large time and has both a positive and negative Bernoulli parameter throughout the post-shock fluid. The CQR solution therefore appears to represent the fluid analog of a marginally bound orbit in a Keplerian potential.

Zoom In Zoom Out Reset image size

This paper is the second of a series that investigates the propagation of weak shocks in the gravitational field of a compact object. Here we have analyzed the stability of new, self-similar solutions presented in Paper I. These solutions, which for brevity we refer to as the "CQR solutions," differ from the Sedov–Taylor solutions in a number of ways; for one, the Mach number of the CQR shock is not infinite but is a number—of order one—that is set by the smooth passage of the fluid quantities through a sonic point within the flow. In addition to outward motion immediately behind the shock, the CQR self-similar solutions contain a stagnation point within the flow, interior to which the fluid falls back onto the point mass at the origin. These solutions therefore simultaneously describe the outward motion of fluid behind a shock and the accretion of matter onto a compact object, and the accretion rate is a function of the properties of the ambient gas. We refer the reader to Paper I for more details of the self-similar solutions.

In this paper we first developed a general formalism for describing the radial perturbations to the shock position and post-shock fluid quantities (the velocity, density, and pressure) and showed that the deviations of the shock properties from self-similarity can be characterized in terms of "eigenmodes." While our approach is similar to that adopted by previous authors, it differs in a number of notable ways (see Section 3.5) and more fully exploits the scale invariance of the problem (which simplifies the stability analysis, not just for the self-similar solution of Paper I but also for the Sedov–Taylor blast wave; see the Appendix). Using this methodology, we showed that the CQR solutions are weakly unstable, with perturbations to the shock position and velocity growing as power laws in time with power-law index ≲0.1 (i.e., perturbations grow with time t as t^α, with α ≲ 0.1).

To investigate the validity of our analytical formalism and results, we ran a suite of one-dimensional, high-resolution simulations with flash (Fryxell et al. 2000), in which we followed the evolution of a shock propagating through a non-self-gravitating atmosphere with density profile ρ ∝ r^−2.5 in equilibrium with a point mass. These simulations confirmed the general predictions of the perturbation theory: when the shock Mach number was very close to that predicted by the CQR solution, the shock evolution closely tracked the predicted ∝ t^2/3 shock scaling, and the post-shock fluid velocity, density, and pressure were described well by the CQR self-similar solutions. However, small deviations between the numerically obtained shock position and the CQR prediction grow with time at a rate that agrees very well with the linear perturbation theory prediction.

When the Mach number of the shock at late times is greater than the CQR value, the perturbation theory predicts that the Mach number should asymptotically grow at a very slow rate, and the solution likely eventually approaches the energy-conserving, Sedov–Taylor blast wave. On the other hand, when the Mach number is less than the CQR value, the linear theory predicts that the Mach number should eventually fall below unity. To understand the long-term behavior of such solutions, we showed that there is a third solution to the self-similar equations with a "shock" Mach number of exactly 1, which corresponds to an RW that travels outward at the local sound speed and informs the hydrostatic medium of the infall at the center. It is likely that these RW solutions characterize the nonlinear, asymptotic evolution of weak shocks with Mach numbers below the CQR value.

Paper I argued that the CQR solutions are likely stable, as larger shock velocities (compared to the CQR value) correspond to greater post-shock energy budgets, but the rate at which binding energy is swept up by the ambient medium is also amplified; it therefore seemed reasonable that these shocks would be capable of self-regulating energetically. However, in addition to this stabilizing term, there is an additional, destabilizing term in the energy equation that is related to the fact that the ambient medium is less bound at larger radii (the ambient energy density falls off as ∝ r^n − 1 when the density scales as ρ ∝ r⁻ⁿ). Thus, if one perturbs the shock position to a larger value, the binding energy at that new radius is slightly reduced, which inhibits the ability of the shock to self-regulate. We find that the magnitude of this second, destabilizing term can slightly exceed that of the stabilizing, velocity-dependent term. It is this slight difference that renders the CQR solutions unstable.

In Paper I, it was shown that the CQR self-similar solution very accurately reproduced the simulation results of Fernández et al. (2018), who numerically studied the propagation of a weak shock through the envelope of a yellow supergiant (evolved from the zero-age main sequence with the stellar evolution code mesa; Paxton et al. 2013, 2015, 2018) following its failure to explode in a traditional supernova. In their simulation, the shock was generated by the sudden decrease of the gravitational field that follows the formation of the proto–neutron star during core collapse. The failed explosion of a supergiant therefore represents a generic, physical scenario in which one expects such weak shocks to arise. In addition to the fact that the density and pressure profiles of the hydrogen envelope of the yellow supergiant were not exactly power laws, the numerical simulation of Fernández et al. (2018) included the effects of a nonadiabatic equation of state and self-gravity, which can all be considered perturbations on top of the conditions necessary to maintain the CQR solution. How, then, did the CQR self-similar solution so accurately describe the propagation of the shock, as well as the time- and space-dependent evolution of all the post-shock fluid variables, over three decades in radius (and four decades in time) when, as we showed here, it is linearly unstable?

The resolution of this apparent contradiction is that the perturbations grow extremely slowly; therefore, as long as the deviations from self-similarity are initially small (e.g., the difference between the shock Mach number and the CQR value is not too large), then the time for the perturbations to grow into the nonlinear regime can be much longer than the time for the shock to propagate through the hydrogen envelope of a supergiant. Thus, the CQR solutions are effectively stable, provided that the initial energy involved in the explosion is not too much greater or less than the binding energy of the hydrogen envelope of the star. Indeed, for the YSG analyzed in Fernández et al. (2018), the neutrino-induced mass loss generates a predicted shock energy of E ≃ 5 × 10⁴⁷ erg, while the total binding energy of the hydrogen envelope is E_bind ∼ 8 × 10⁴⁷ erg (see the left panel of Figure 14 of Paper I).

In this paper, we restricted our perturbation analysis to purely radial modes, to which the CQR solutions—in contrast to the Sedov–Taylor solutions (see the Appendix)—are unstable. One can generalize the perturbation approach to include nonspherical modes by decomposing the shock front into spherical harmonics, which in turn induce aspherical, post-shock perturbations on the fluid quantities (see, e.g., Vishniac 1983). However, the physical origin of the instability found in this paper for radial modes is related to the fact that, as the shock moves outward, it sweeps up progressively less binding energy from the ambient medium. This effect is maximized geometrically when the perturbations are spherical. For example, an ℓ = 1 mode generates a "sloshing" motion of the shock, meaning that roughly half of the shock receives a positive perturbation in radial velocity, while the other half receives a negative radial velocity perturbation. Therefore, the net reduction in the rate at which binding energy is supplied to the post-shock fluid is smaller in this case, which will likewise result in a smaller growth rate of any potential instability. We therefore expect that, while nonradial modes may be unstable, radial perturbations are likely the most unstable from the energy standpoint considered in this paper (which neglects, e.g., buoyancy, which can lead to distinct, unstable modes). We do, however, plan to analyze the stability of the CQR solution to nonspherical perturbations in a future investigation.

E.R.C. acknowledges support from NASA through the Einstein Fellowship Program, grant PF6-170170. This work was supported in part by a Simons Investigator award from the Simons Foundation (EQ) and the Gordon and Betty Moore Foundation through grant GBMF5076.

Software: FLASH (Fryxell et al. 2000).

In this appendix, we use the formalism for analyzing the stability of the CQR solution to assess the stability of the Sedov–Taylor blast wave to radial perturbations, and we compare our results to those of Ryu & Vishniac (1987). The Sedov–Taylor blast wave adopts the same functional form for the self-similar variables as the CQR solution, namely, that the velocity, v, the density, ρ, and the pressure, p, are written self-similarly as

$$\begin{eqnarray}\begin{array}{rclcl}v & = & V\left({f}_{0}+{f}_{1}\right),\quad \rho ={\rho }_{{\rm{a}}}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}\left({g}_{0}+{g}_{1}\right),\,p & = & {\rho }_{{\rm{a}}}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{-n}{V}^{2}\left({h}_{0}+{h}_{1}\right),\end{array}\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray*}\begin{array}{rcl}\Rightarrow \quad s & = & \displaystyle \frac{p}{{\rho }^{{\gamma }_{2}}}=\mathrm{ln}\left({V}^{2}{R}^{-n+n{\gamma }_{2}}\right)+\mathrm{ln}\left(\displaystyle \frac{{h}_{0}}{{g}_{0}^{{\gamma }_{2}}}\right)+\displaystyle \frac{{h}_{1}}{{h}_{0}}-\displaystyle \frac{{\gamma }_{2}{g}_{1}}{{g}_{0}}\equiv \mathrm{ln}\left({V}^{2}{R}^{-n+n{\gamma }_{2}}\right)+{s}_{0}+{s}_{1},\end{array}\end{eqnarray*}$$

where R is the true shock position and $V={dR}/{dt}$ is the shock velocity. Since the change of variables is the same, the governing equations are identical, except that we ignore the gravitational term on the right-hand side of the momentum equation:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {g}_{1}}{\partial \tau }-n\left({g}_{0}+{g}_{1}\right)-\xi \left({g}_{0}^{{\prime} }+\displaystyle \frac{\partial {g}_{1}}{\partial \xi }\right)+\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left[{\xi }^{2}\left({g}_{0}+{g}_{1}\right)\left({f}_{0}+{f}_{1}\right)\right]=0,\end{array}\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{1}}{\partial \tau }+\left(\displaystyle \frac{1}{2}\displaystyle \frac{d}{d\tau }\mathrm{ln}{V}^{2}\right)\left({f}_{0}+{f}_{1}\right)-\xi \left({f}_{0}^{{\prime} }+\displaystyle \frac{\partial {f}_{1}}{\partial \xi }\right)+\left({f}_{0}+{f}_{1}\right)\displaystyle \frac{\partial }{\partial \xi }\left({f}_{0}+{f}_{1}\right)+\displaystyle \frac{1}{{g}_{0}}{\left(1+\displaystyle \frac{{g}_{1}}{{g}_{0}}\right)}^{-1}\,\displaystyle \frac{\partial }{\partial \xi }\left({h}_{0}+{h}_{1}\right)=0,\end{array}\end{eqnarray} \tag{ 93 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {s}_{1}}{\partial \tau }+\displaystyle \frac{d}{d\tau }\mathrm{ln}{V}^{2}-n+n\gamma -\xi \displaystyle \frac{\partial }{\partial \xi }\left({s}_{0}+{s}_{1}\right)+\left({f}_{0}+{f}_{1}\right)\displaystyle \frac{\partial }{\partial \xi }\left({s}_{0}+{s}_{1}\right)=0.\end{array}\end{eqnarray} \tag{ 94 }$$

The boundary conditions are still governed by the jump conditions at the shock front, namely, Equations (6)–(8), but here we adopt the strong-shock limit, which makes them extremely simple:

$$\begin{eqnarray}&&{f}_{0}(\tau ,1)+{f}_{1}(\tau ,1)=\displaystyle \frac{2}{\gamma +1},\end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&{g}_{0}(\tau ,1)+{g}_{1}(\tau ,1)=\displaystyle \frac{\gamma +1}{\gamma -1},\end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray}&&{h}_{0}(\tau ,1)+{h}_{1}(\tau ,1)=\displaystyle \frac{2}{\gamma +1}.\end{eqnarray} \tag{ 97 }$$

Notice that these boundary conditions are this simple only because we are defining the self-similar variable in terms of the total shock position and velocity; if one instead opts to initially break up the shock velocity and position into its unperturbed and perturbed quantities, these become more complicated (see Equation (5.17) of Ryu & Vishniac 1987).

Investigating Equations (92)–(97), we see that we can satisfy both the differential equations and the boundary conditions with all subscript-1 quantities identically equal to zero if

$$\begin{eqnarray}&&{V}^{2}\propto {R}^{\alpha },\end{eqnarray} \tag{ 98 }$$

with α any arbitrary constant. However, the energy behind the shock must be conserved for any realistic system (under this set of assumptions, where we ignore gravity and the Mach number is infinite), and we can write this energy as

$$\begin{eqnarray}\begin{array}{rcl}E & = & 4\pi {\displaystyle \int }_{0}^{R(t)}\left(\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{p}{\rho }\right)\rho {r}^{2}{dr}\,=\,4\pi {\rho }_{{\rm{a}}}{r}_{{\rm{a}}}^{3}{V}^{2}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{3-n}{\displaystyle \int }_{0}^{1}\left({g}_{0}+{g}_{1}\right)\,\left[{\left({f}_{0}+{f}_{1}\right)}^{2}+\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{{h}_{0}+{h}_{1}}{{g}_{0}+{g}_{1}}\right]{\xi }^{2}d\xi .\end{array}\end{eqnarray} \tag{ 99 }$$

There are combinations of n and γ for which the Sedov–Taylor blast wave only extends to a finite inner radius that serves as a contact discontinuity (see Goodman 1990), and in this case one replaces the lower bound in the integral with some ξ₀; for simplicity, however, we focus here on the case when the entire post-shock fluid can be described by the Sedov–Taylor self-similar solutions. We therefore see that if we let

$$\begin{eqnarray}&&4\pi {\rho }_{{\rm{a}}}{r}_{{\rm{a}}}^{3}{V}^{2}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{3-n}={E}_{0},\end{eqnarray} \tag{ 100 }$$

with E₀ a constant, then the energy is conserved and all of the perturbations can be set to zero. Thus, for the Sedov–Taylor blast wave we write

$$\begin{eqnarray}&&4\pi {\rho }_{{\rm{a}}}{r}_{{\rm{a}}}^{3}{V}_{}^{2}{\left(\displaystyle \frac{R}{{r}_{{\rm{a}}}}\right)}^{3-n}={E}_{0}\left[1+2{\zeta }_{}(\tau )\right]\end{eqnarray} \tag{ 101 }$$

$$\begin{eqnarray*}&&\Rightarrow \,\mathrm{ln}{V}^{2}\simeq \left(n-3\right)\tau +2{\zeta }_{}(\tau ),\end{eqnarray*}$$

where ζ(τ) is a small quantity that encodes the perturbations to the position and velocity of the shock, and the second line follows in the perturbative limit (i.e., we dropped all nonlinear terms in ζ). The equations for the unperturbed quantities are then

$$\begin{eqnarray}&&-{{ng}}_{0}-\xi {g}_{0}^{{\prime} }+\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{d}{d\xi }\left({\xi }^{2}{g}_{0}{f}_{0}\right)=0,\end{eqnarray} \tag{ 102 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\left(n-3\right){f}_{0}+\left({f}_{0}-\xi \right){f}_{0}^{{\prime} }+\displaystyle \frac{{h}_{0}^{{\prime} }}{{g}_{0}}=0,\end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray}&&-3+n\gamma +({f}_{0}-\xi )\displaystyle \frac{d}{d\xi }\mathrm{ln}\left(\displaystyle \frac{{h}_{0}}{{g}_{0}^{\gamma }}\right)=0,\end{eqnarray} \tag{ 104 }$$

while those for the perturbations read

$$\begin{eqnarray}&&\displaystyle \frac{\partial {g}_{1}}{\partial \tau }-{{ng}}_{1}-\xi \displaystyle \frac{\partial {g}_{1}}{\partial \xi }+\displaystyle \frac{1}{{\xi }^{2}}\displaystyle \frac{\partial }{\partial \xi }\left({\xi }^{2}{g}_{0}{f}_{1}+{\xi }^{2}{f}_{0}{g}_{1}\right)=0,\end{eqnarray} \tag{ 105 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{1}}{\partial \tau }+\displaystyle \frac{1}{2}(n-3){f}_{1}-\xi \displaystyle \frac{\partial {f}_{1}}{\partial \xi }+\displaystyle \frac{\partial }{\partial \xi }\left({f}_{0}{f}_{1}\right)+\displaystyle \frac{1}{{g}_{0}}\displaystyle \frac{\partial }{\partial \xi }\left({h}_{0}{s}_{1}+{\gamma }_{2}{h}_{0}\displaystyle \frac{{g}_{1}}{{g}_{0}}\right)-\displaystyle \frac{{g}_{1}}{{g}_{0}^{2}}{h}_{0}^{{\prime} }=-\dot{\zeta }{f}_{0},\end{array}\end{eqnarray} \tag{ 106 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {s}_{1}}{\partial \tau }+\left({f}_{0}-\xi \right)\displaystyle \frac{\partial {s}_{1}}{\partial \xi }+{f}_{1}{s}_{0}^{{\prime} }=-2\dot{\zeta }.\end{eqnarray} \tag{ 107 }$$

The boundary conditions for the unperturbed quantities are the usual strong-shock jump conditions,

$$\begin{eqnarray}&&{f}_{0}(1)={h}_{0}(1)=\displaystyle \frac{2}{\gamma +1},\quad {g}_{0}(1)=\displaystyle \frac{\gamma +1}{\gamma -1},\end{eqnarray} \tag{ 108 }$$

while the perturbations satisfy homogeneous boundary conditions,

$$\begin{eqnarray}&&{f}_{1}(1,\tau )={g}_{1}(1,\tau )={h}_{1}(1,\tau )=0.\end{eqnarray} \tag{ 109 }$$

Equations (105)–(107) can be written in matrix form as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial \tau }\left(\begin{array}{c}{s}_{1}\\ {F}_{1}\\ {G}_{1}\end{array}\right)\,+\,\left(\begin{array}{ccc}{f}_{0}-\xi & 0 & 0\\ {\xi }^{2}{h}_{0} & {f}_{0}-\xi & \gamma {\xi }^{2}{h}_{0}\\ 0 & {\xi }^{-2}{g}_{0}^{-1} & {f}_{0}-\xi \end{array}\right)\displaystyle \frac{\partial }{\partial \xi }\left(\begin{array}{c}{s}_{1}\\ {F}_{1}\\ {G}_{1}\end{array}\right)+\left[\begin{array}{ccc}0 & {\xi }^{-2}{g}_{0}^{-1}{s}_{0}^{{\prime} } & 0\\ {\xi }^{2}{h}_{0}^{{\prime} } & \displaystyle \frac{1}{2}-\displaystyle \frac{n}{2}+2{f}_{0}^{{\prime} } & {\xi }^{2}(\gamma -1){h}_{0}^{{\prime} }\\ 0 & 0 & 0\end{array}\right]\,\left(\begin{array}{c}{s}_{1}\\ {F}_{1}\\ {G}_{1}\end{array}\right)\,=\,\left(\begin{array}{c}-2\dot{\zeta }\\ -{f}_{0}{g}_{0}{\xi }^{2}\dot{\zeta }\\ 0\end{array}\right).\end{array}\end{eqnarray} \tag{ 110 }$$

Notice that this is exactly the same as Equation (27) if we make the substitutions $\zeta \to -\zeta$ and $n\to -1+2n$ and the gravitational term on the right-hand side is missing.

As was true for the perturbations to the CQR self-similar solutions, Equation (110) appears overconstrained because of the introduction of the variable ζ. There is, however, an additional constraint that these variables must satisfy: taking the energy equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\mathscr{E}}}{\partial t}+\displaystyle \frac{\partial {\mathscr{F}}}{\partial r}=0,\end{eqnarray} \tag{ 111 }$$

where

$$\begin{eqnarray}&&{\mathscr{E}}=\left(\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{p}{\rho }\right)\rho {r}^{2}\end{eqnarray} \tag{ 112 }$$

is the energy density and

$$\begin{eqnarray}&&{\mathscr{F}}=\left(\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{p}{\rho }\right)\rho {r}^{2}v\end{eqnarray} \tag{ 113 }$$

is the energy flux, and integrating from zero to $R(t)+\epsilon$ and taking the limit as $\epsilon \to 0$ gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}={\mathscr{F}}(0),\end{eqnarray} \tag{ 114 }$$

where

$$\begin{eqnarray}&&E={\int }_{0}^{R}\left(\displaystyle \frac{1}{2}{v}^{2}+\displaystyle \frac{1}{\gamma -1}\displaystyle \frac{p}{\rho }\right)\rho {r}^{2}{dr}\end{eqnarray} \tag{ 115 }$$

is the total—including the perturbations—energy behind the shock. Since we are not introducing any sources or sinks of energy, this quantity must be exactly conserved, meaning that the perturbations to this quantity must also be identically zero. Thus, if we expand the flux into its perturbed and unperturbed components, we recover the boundary condition

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\xi \to 0}\left(\displaystyle \frac{{F}_{\sigma }}{{g}_{0}}+{f}_{0}{\xi }^{2}{s}_{\sigma }+\gamma {f}_{0}{\xi }^{2}{G}_{\sigma }\right)=\mathop{\mathrm{lim}}\limits_{\xi \to 0}\left({\xi }^{2}{f}_{\sigma }+{f}_{0}{\xi }^{2}{h}_{\sigma }\right)=0,\end{eqnarray} \tag{ 116 }$$

where in the last line we wrote this condition in terms of the self-similar velocity and pressure (Equation (91)). Notice that, if the first term in parentheses and the geometric factors are ignored, this boundary condition is in agreement with that of Ryu & Vishniac (1987), being that the correction to the pressure vanishes at the origin. However, their boundary condition is not in fact correct since ${f}_{0}{\xi }^{2}$ goes to zero at the origin for the Sedov–Taylor blast wave; Equation (116) is more general, and any, weakly diverging correction to the pressure (one that does not diverge faster than ∝1/ξ³, since the unperturbed, Sedov–Taylor velocity satisfies f₀∝ ξ) will satisfy the conservation of energy provided that the velocity does not diverge as 1/ξ² or faster.

We now search for eigenmodes of the form

$$\begin{eqnarray}\begin{array}{rclcl}\zeta & = & {\zeta }_{\sigma }{e}^{\sigma \tau },\quad {F}_{1}={\zeta }_{\sigma }{e}^{\sigma \tau }{F}_{\sigma },\,{G}_{1} & = & {\zeta }_{\sigma }{e}^{\sigma \tau }{G}_{\sigma },\quad {s}_{1}={\zeta }_{\sigma }{e}^{\sigma \tau }{s}_{\sigma }.\end{array}\end{eqnarray} \tag{ 117 }$$

As was true for the CQR solution, ζ_σ scales out of the problem, and the value of σ is determined by the requirement that the functions satisfy the fourth boundary condition at the origin. The eigenvalue equation is

$$\begin{eqnarray}&&\begin{array}{l}\sigma \left(\begin{array}{c}{s}_{\sigma }\\ {F}_{\sigma }\\ {G}_{\sigma }\end{array}\right)\,+\,\left(\begin{array}{ccc}{f}_{0}-\xi & 0 & 0\\ {\xi }^{2}{h}_{0} & {f}_{0}-\xi & {\gamma }_{2}{\xi }^{2}{h}_{0}\\ 0 & {\xi }^{-2}{g}_{0}^{-1} & {f}_{0}-\xi \end{array}\right)\displaystyle \frac{\partial }{\partial \xi }\left(\begin{array}{c}{s}_{\sigma }\\ {F}_{\sigma }\\ {G}_{\sigma }\end{array}\right)+\left(\begin{array}{ccc}0 & {\xi }^{-2}{g}_{0}^{-1}{s}_{0}^{{\prime} } & 0\\ {\xi }^{2}{h}_{0}^{{\prime} } & \displaystyle \frac{1}{2}-\displaystyle \frac{n}{2}+2{f}_{0}^{{\prime} } & {\xi }^{2}({\gamma }_{2}-1){h}_{0}^{{\prime} }\\ 0 & 0 & 0\end{array}\right)\,\left(\begin{array}{c}{s}_{\sigma }\\ {F}_{\sigma }\\ {G}_{\sigma }\end{array}\right)\,=\,\left(\begin{array}{c}-2\sigma \\ -{f}_{0}{g}_{0}{\xi }^{2}\sigma \\ 0\end{array}\right).\end{array}\end{eqnarray} \tag{ 118 }$$

Notice that, if σ = 0, the right-hand side of this matrix equation is zero. Moreover, if we let F_σ = G_σ = s_σ = 0, then the differential equations and, since they are homogeneous, the boundary conditions are satisfied. This is therefore an exact, neutrally stable solution for the perturbations to the Sedov–Taylor blast wave. From our definition of ζ (Equation (101)), we find that the unperturbed and perturbed shock positions satisfy

$$\begin{eqnarray}&&\begin{array}{l}4\pi {\rho }_{{\rm{a}}}{r}_{{\rm{a}}}^{3}{V}_{0}^{2}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{3-n}\left(1+2\displaystyle \frac{{V}_{1}}{{V}_{0}}\right)\left[1+\left(3-n\right)\displaystyle \frac{{R}_{1}}{{R}_{0}}\right]={E}_{0}\left(1+2\zeta \right),\end{array}\end{eqnarray} \tag{ 119 }$$

$$\begin{eqnarray}&&4\pi {\rho }_{{\rm{a}}}{r}_{{\rm{a}}}^{3}{V}_{0}^{2}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{3-n}={E}_{0},\end{eqnarray} \tag{ 120 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{1}}{{R}_{0}}=\displaystyle \frac{2\zeta }{5-n}\left(1-{e}^{-\tfrac{5-n}{2}\tau }\right),\end{eqnarray} \tag{ 121 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{V}_{1}}{{V}_{0}}=\displaystyle \frac{2\zeta }{5-n}\left(1+\displaystyle \frac{3-n}{2}{e}^{-\tfrac{5-n}{2}\tau }\right).\end{eqnarray} \tag{ 122 }$$

From these expressions, we see that, similar to the CQR solution, the perturbations to the shock position and velocity are written as the difference between the eigenmode (which here has an eigenvalue of σ = 0) and a second "mode" that has an eigenvalue of $-(5-n)/2$. However, this second mode is just a consequence of the initial conditions and is not a true eigenvalue from the standpoint that it does not characterize the temporal dependence of ζ. If we further use the fact that the unperturbed shock position for the Sedov–Taylor blast wave scales as ${R}_{0}(t)\propto {t}^{2/(5-n)}$ (this can be verified directly by integrating Equation (120)), then this second mode decays as t⁻¹, which is identical to the temporal scaling of the trivial, decaying mode of the CQR solution.

While the perturbations to the post-shock fluid quantities may appear to be trivial (i.e., ${F}_{\sigma }={G}_{\sigma }={s}_{\sigma }\equiv 0$), this is in fact not the case because of our definitions of the velocity, density, and pressure in terms of the total shock velocity and position. Specifically, the solutions to first order in the unperturbed self-similar variable, or physical radius r, are

$$\begin{eqnarray}\begin{array}{rcl}v({\xi }_{0},\tau ) & = & {V}_{0}\left[{f}_{0}+\displaystyle \frac{2\zeta }{5-n}\left(1+\displaystyle \frac{3-n}{2}{e}^{-\tfrac{5-n}{2}\tau }\right)\right.\,\left.{f}_{0}-\displaystyle \frac{2\zeta }{5-n}\left(1-{e}^{-\tfrac{5-n}{2}\tau }\right){\xi }_{0}{f}_{0}^{{\prime} }\right],\end{array}\end{eqnarray} \tag{ 123 }$$

$$\begin{eqnarray}&&\begin{array}{l}\rho ({\xi }_{0},\tau )={\rho }_{{\rm{a}}}\,{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{-n}\left[{g}_{0}-\displaystyle \frac{2\zeta }{5-n}\left(1-{e}^{-\tfrac{5-n}{2}\tau }\right)\left({{ng}}_{0}+{\xi }_{0}{g}_{0}^{{\prime} }\right)\right],\end{array}\end{eqnarray} \tag{ 124 }$$

$$\begin{eqnarray}\begin{array}{rcl}p({\xi }_{0},\tau ) & = & {\rho }_{{\rm{a}}}{V}_{0}^{2}{\left(\displaystyle \frac{{R}_{0}}{{r}_{{\rm{a}}}}\right)}^{-n}\left[{h}_{0}+\displaystyle \frac{4\zeta }{5-n}\left(1+\displaystyle \frac{3-n}{2}{e}^{-\tfrac{5-n}{2}\tau }\right){h}_{0}\right.-\left.\displaystyle \frac{2\zeta }{5-n}\left(1-{e}^{-\tfrac{5-n}{2}\tau }\right)\left({{nh}}_{0}+{\xi }_{0}{h}_{0}^{{\prime} }\right)\right].\end{array}\end{eqnarray} \tag{ 125 }$$

These solutions can be shown to satisfy Equations 4.18(a)–4.18(d) of Ryu & Vishniac (1987) and the boundary conditions at the shock in terms of the unperturbed self-similar variable (their Equations 5.17(a)–5.17(d)) when the perturbations are spherically symmetric (ℓ = 0).

Notice that our solution for the perturbation to the pressure does not go to zero at the origin, which is the fourth boundary condition imposed by Ryu & Vishniac (1987). However, our boundary condition—which demands that the flux of energy at the origin be zero and hence that the total energy be conserved—is satisfied for any nondiverging pressure and velocity (specifically Equation (116)). Our solutions, being written in terms of the well-behaved, unperturbed solutions, are clearly convergent near the origin, and Figure 14 shows this convergence explicitly for n = 0 and γ = 1.5. On the other hand, the solutions that satisfy the vanishing pressure boundary condition of Ryu & Vishniac (1987), which are shown in the right panel of this figure (and have a decaying eigenmode of σ ≃ −0.3495), possess a perturbation to the velocity that diverges as ∝1/ξ² near the origin. Therefore, the solutions found by Ryu & Vishniac (1987) do not globally conserve energy owing to their finite energy flux at the origin.

Zoom In Zoom Out Reset image size