Wave-driven Shocks in Stellar Outbursts: Dynamics, Envelope Heating, and Nascent Blast Waves

In general, unlike the simple problem, the cross-sectional area, fluid density, and sound speed may all vary with r. A further modification is thus required to adapt the equal-area method to more complicated environments like stellar envelopes, with their spherical geometries and strongly stratified density profiles. For this, Lighthill (Waves in Fluids, Lighthill 1978, ch. 2) proposed to generalize Whitham's method by replacing J_± with an approximate invariant (a linearized adiabatic invariant, to be precise) of acoustical motions in the more complicated problem. Lighthill chose the approximate invariant $\delta p/\sqrt{Z}$, where δ p is the pressure perturbation and Z(r) = ρ₀ c₀/A is the acoustic impedance. However we prefer to use the "wave-front amplitude,"

$$\begin{eqnarray}&&{{ \mathcal R }}_{\pm }=\sqrt{{\rho }_{0}{c}_{0}A}\,\,\delta {J}_{\pm }\end{eqnarray} \tag{ 2 }$$

(Matzner & Ro 2020, hereafter MR20). Here, δ J_± is the acoustical perturbation of J_± from the background state. Our quantity ${{ \mathcal R }}_{\pm }$ is also an approximate invariant. Its advantage is that it provides a more direct connection to J_± and also to the acoustical wave luminosity, which for waves of linear amplitude is

$$\begin{eqnarray}&&{L}_{w}=\displaystyle \frac{1}{4}({{ \mathcal R }}_{+}^{2}-{{ \mathcal R }}_{-}^{2}).\end{eqnarray} \tag{ 3 }$$

Conservation of ${{ \mathcal R }}_{+}$ along characteristics in a purely outward-traveling wave (${{ \mathcal R }}_{-}=0$) is therefore equivalent to the conservation of L_w in a wave that neither reflects nor dissipates.

In the generalized method, the waveform shears at a changing rate that depends on the local conditions. This can be accommodated by defining a new evolution coordinate y(r) to replace t, as we describe in more detail in Section 3. Before we do, we wish to review a couple of results from MR20 that help to clarify the validity of the generalized method.

First, it is important to note that the quantity ${{ \mathcal R }}_{\pm }$ is not strictly conserved along characteristics, but rather varies in linear theory according to

$$\begin{eqnarray}&&\displaystyle \frac{d}{{{dr}}^{\pm }}{{ \mathcal R }}_{\pm }={{ \mathcal R }}_{\mp }\displaystyle \frac{\partial }{\partial r}\mathrm{ln}\sqrt{Z},\end{eqnarray} \tag{ 4 }$$

where d/dr^± is the derivative along C^± with respect to r. The right-hand side causes reflection whenever the impedance Z varies significantly on scales of order the wavelength, but is completely negligible when it does not (the high-frequency limit). Therefore, Lighthill's generalized equal-area method requires that the wave not only be linear in amplitude but high in frequency relative to its surroundings. This is not a strong constraint in the process of shock formation, because the dynamical timescales are much shorter in the core than in the envelope of an evolved star, and because a wave must be low amplitude to travel many wavelengths before it forms a shock. However, if the wave were to meet a rapid change in Z (at a composition boundary, for instance), then one would need to evaluate what fraction is transmitted and apply our technique to the transmitted wave. Reflection may also affect shocks that accelerate near the stellar surface, as we discuss in Section 6.3.

Second, Equation (4) shows that ${{ \mathcal R }}_{\pm }$ is an adiabatic invariant of the linearized problem, but shock formation is a nonlinear process. In MR20, we show that nonlinear adiabatic invariants, when these can be found, almost always tend toward ${{ \mathcal R }}_{\pm }$ in the linear limit. More generally, we find that the nonlinear error one incurs by assuming that ${{ \mathcal R }}_{\pm }$ is constant is second order in the wave amplitude and also proportional to the logarithmic gradient of a combination of environmental quantities. Note also that the "conserved area" $\int {{ \mathcal R }}_{+}({t}_{i}){{dt}}_{i}$ has units of the square root of the action, which should be conserved in adiabatic evolution. These results give us confidence that Lighthill's proposed generalization is sound.

A shock forms where ∂f/∂t becomes infinite; this occurs first on those characteristics, labeled t_i,s , that are local maxima of the compression rate and therefore local maxima of $f^{\prime} ({t}_{i})$. (This is a necessary but not always sufficient criterion: it is possible for such characteristics to collide with another shock before forming their own.) The shock formation radius R_s is found from the relation

$$\begin{eqnarray}&&{y}_{s}\equiv y({R}_{{\rm{s}}})=\displaystyle \frac{1}{f^{\prime} ({t}_{i,s})}.\end{eqnarray} \tag{ 10 }$$

Nearby characteristics will collide with the shock: these correspond to the sections of the waveform excised in the equal-area method.

The shock formation criterion, Equation (10), was derived by Lighthill (1978, his Equation (254)), and rederived by Lin & Szeri (2001), Tyagi & Sujith (2003, 2005, who noted agreement with Lighthill's criterion) and us (Ro & Matzner 2017). As we discuss in MR20, it assumes the wave's slope is not altered by reflection at some location where Z changes rapidly; this is more of a restriction for internal shocks than head shocks.

To address the trajectory R(t) of a shock after it forms, we follow the characteristic-tracking procedure of Friedrichs (1948) and Whitham (1974). The velocity V = dR/dt of a weak shock is just the average of the characteristic velocities immediately to each side of it:

$$\begin{eqnarray}&&V=\displaystyle \frac{{(v+c)}_{u}+{(v+c)}_{d}}{2},\end{eqnarray} \tag{ 11 }$$

to first order in v/c₀ (Courant & Friedrichs 1948), where u and d refer to the states immediately upstream and downstream of the shock, respectively. This condition enforces the equal-area rule we discussed in Section 2.1, but in differential form.

In order to evaluate V, we require the rate at which characteristics arrive at the shock front from either side:

$$\begin{eqnarray*}&&{\left.\displaystyle \frac{{{dt}}_{i}}{{dt}}\right|}_{s}={\partial }_{t}{t}_{i}+V{\partial }_{r}{t}_{i}=\displaystyle \frac{v+c-V}{v+c}{\partial }_{t}{t}_{i};\end{eqnarray*}$$

the second equality uses ∂_t t_i  = − (v + c)∂_r t_i . Substituting Equation (8),

$$\begin{eqnarray*}&&\left(1-f^{\prime} ({t}_{i})y\right){\left.\displaystyle \frac{{{dt}}_{i}}{{dt}}\right|}_{s}\simeq \displaystyle \frac{v+c-V}{v+c},\end{eqnarray*}$$

which can be evaluated on either side of the shock. Combining this with Equation (11) and using $(q+1)\sqrt{{L}_{p}/{L}_{\max }}{\left.d/{dt}\right|}_{s}\,={\left.d/{dy}\right|}_{s}$ gives coupled differential equations for the initial times of characteristics as they arrive at the shock:

$$\begin{eqnarray}&&2\left(1-f^{\prime} ({t}_{{id}})y\right)\displaystyle \frac{{{dt}}_{{id}}}{{dy}}\simeq f({t}_{{id}})-f({t}_{{iu}}),\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&2\left(1-f^{\prime} ({t}_{{iu}})y\right)\displaystyle \frac{{{dt}}_{{iu}}}{{dy}}\simeq f({t}_{{iu}})-f({t}_{{id}}).\end{eqnarray} \tag{ 12b }$$

The shock trajectory, described by t_iu (y) and t_id (y), is found by integrating these equations from the point of shock formation, t_iu  = t_id  = t_i,s . When it is desirable to work with dimensionless solutions scaled to the formation of a given shock, we note that Equations 12(a) and 12(b) remain valid and unchanged with the replacements ${t}_{i}\to {\tilde{t}}_{i}$ and $y\to \tilde{y}$, where

$$\begin{eqnarray*}&&{\tilde{t}}_{i}\equiv \displaystyle \frac{{t}_{i}-{t}_{i,s}}{{y}_{s}},\,\,\,\,\,\,\tilde{y}(r)\equiv \displaystyle \frac{y(r)}{{y}_{s}}=\displaystyle \frac{Y(r)}{{Y}_{s}},\end{eqnarray*}$$

and Y_s means Y(R_s ).

To specify the deposition of wave energy, we define a normalized energy variable $\tilde{E}=E/({L}_{p}{y}_{s})$. Equation (1) becomes

$$\begin{eqnarray}&&\displaystyle \frac{d{\tilde{E}}_{\mathrm{diss}}}{{dr}}=\displaystyle \frac{{\left({\rm{\Delta }}f\right)}^{3}}{6}\displaystyle \frac{d\tilde{y}}{{dr}}=\displaystyle \frac{d{\tilde{E}}_{\mathrm{diss}}}{d\tilde{y}}\displaystyle \frac{d\tilde{y}}{{dr}}.\end{eqnarray} \tag{ 13 }$$

We see from this that shock dissipation can be broken into two parts. The first, $d{\tilde{E}}_{\mathrm{diss}}/d\tilde{y}={({\rm{\Delta }}f)}^{3}/6$, depends only on the form of the initial pulse and the evolution coordinate $\tilde{y}$. The second, $d\tilde{y}/{dr}$, depends only on the wave amplitude and the structure of the background envelope.

Consider a sinusoidal wave of N cycles, $f=\sin \omega {t}_{i}$ for 0 < ω t_i  < 2N π. We assume that it begins and ends in a state of compression, so that head and tail shocks form in addition to its internal shocks. All of these shocks share the common value y_s  = 1/ω, so $\tilde{y}=\omega y$ and one can define ${\tilde{t}}_{i}=\omega ({t}_{i}-{t}_{i,s})$ for any shock labeled by t_i,s . The example in Figure 1 consists of two cycles.

4.1.1. Internal Shocks

Internal shocks evolve according to Equation (15) with b = 1 and f(t_i,s ) = 0, so $\tilde{y}({\tilde{t}}_{i})={\tilde{t}}_{i}/\sin ({\tilde{t}}_{i})$ and

$$\begin{eqnarray}&&{\tilde{t}}_{i}={\mathrm{sinc}}^{-1}\displaystyle \frac{1}{\tilde{y}(r)}.\end{eqnarray} \tag{ 16 }$$

This coincides with what one would obtain by construction for the internal shock in Figure 1.

Symmetry around the internal shocks implies u_u  = −u_d and so Δf = 2f_d . Evaluating Equation (13),

$$\begin{eqnarray}&&\displaystyle \frac{d{\tilde{E}}_{\mathrm{diss}}}{d\tilde{y}}=\displaystyle \frac{4}{3}{\left[\sin \left({\mathrm{sinc}}^{-1}\displaystyle \frac{1}{\tilde{y}}\right)\right]}^{3},\end{eqnarray} \tag{ 17 }$$

and the cumulative dissipated energy is

$$\begin{eqnarray}&&{\tilde{E}}_{\mathrm{diss}}={\tilde{t}}_{i}-\sin {\tilde{t}}_{i}\cos {\tilde{t}}_{i}-\displaystyle \frac{2}{3}{\tilde{t}}_{i}{\left(\sin {\tilde{t}}_{i}\right)}^{2},\end{eqnarray} \tag{ 18 }$$

which can be converted to ${\tilde{E}}_{\mathrm{diss}}(\tilde{y})$ using Equation (16). We obtained this by integrating the wave energy in the equal-area portions of the initial waveform that will be removed by the shock (see Figure 1); we then checked for numerical consistency with Equation (17).

The total energy deposited by an internal shock is ${\tilde{E}}_{\mathrm{diss}}(\infty )={\int }_{1}^{\infty }(d{\tilde{E}}_{\mathrm{diss}}/d\tilde{y})d\tilde{y}=\pi$, or in dimensional terms, E_diss( ∞ ) = π L_p /ω. This matches the initial wave energy carried by characteristics that will intersect the shock. We will also be interested in what fraction of the wave energy survives at late times. At large values of $\tilde{y}$, $d{\tilde{E}}_{\mathrm{diss}}/d\tilde{y}\to (4{\pi }^{3}/3){\tilde{y}}^{-3}$, so the residual wave energy that will ultimately interact with the shock is

$$\begin{eqnarray}&&{\tilde{E}}_{w}(\tilde{y})\to \displaystyle \frac{2}{3}{\pi }^{3}{\tilde{y}}^{-2}.\end{eqnarray} \tag{ 19 }$$

4.1.2. Head and Tail Shocks

Head and tail shocks evolve in opposite directions but otherwise in exactly the same way (in linear theory), so we focus on the case of a head shock. Using b = 2 and f(t_i,s ) = 0, Equation (15) yields $\tilde{y}=2/(1+\cos {\tilde{t}}_{i})$, so

$$\begin{eqnarray}&&{\tilde{t}}_{i}={\cos }^{-1}\left(\displaystyle \frac{2}{\tilde{y}}-1\right).\end{eqnarray} \tag{ 20 }$$

In this case Δf = f(t_i );  using $\sin [{\cos }^{-1}(x)]=\sqrt{1-{x}^{2}}$,

$$\begin{eqnarray*}&&{\rm{\Delta }}f=2\displaystyle \frac{\sqrt{\tilde{y}(r)-1}}{\tilde{y}(r)};\end{eqnarray*}$$

and by Equation (13)),

$$\begin{eqnarray}&&\displaystyle \frac{d{\tilde{E}}_{\mathrm{diss}}}{d\tilde{y}}=\displaystyle \frac{4}{3}\displaystyle \frac{{\left(\tilde{y}-1\right)}^{3/2}}{{\tilde{y}}^{3}}.\end{eqnarray} \tag{ 21 }$$

The total energy deposited by a head or tail shock, ${\tilde{E}}_{\mathrm{diss}}(\infty )=\pi /2$, is only half of what is deposited at an internal shock, because external shocks interact with characteristics from only half a wave period. However, external shocks expend this energy much more slowly. For large $\tilde{y}$, external shocks obey $d{\tilde{E}}_{\mathrm{diss}}/d\tilde{y}\to (4/3){\tilde{y}}^{-3/2}$, so the residual wave energy is

$$\begin{eqnarray}&&{\tilde{E}}_{w}(\tilde{y})\to \displaystyle \frac{8}{3}{\tilde{y}}^{-1/2},\end{eqnarray} \tag{ 22 }$$

in sharp contrast to internal shocks (Equation (19)). For this reason, internal shocks deliver more energy to radii of order the shock formation radius, while leading and trailing shocks have a much greater potential to create strong shocks in stellar atmospheres. Figure 2 compares the normalized residual energy ${\tilde{E}}_{w}(\tilde{y})$ between internal and head/tail shocks.

Zoom In Zoom Out Reset image size

Although we concentrate on purely sinusoidal wave trains, it is important to realize that the asymptotic properties of shock decay are not sensitive to the details of the waveform. For head and tail shocks, for instance, the limit $f({t}_{{id}})\to 2{\tilde{y}}^{-1/2}$ is a special case of the result

$$\begin{eqnarray}&&f{\left(y\right)}^{2}=\displaystyle \frac{2}{y}\int f({t}_{i})\,{{dt}}_{i},\end{eqnarray} \tag{ 23 }$$

which can be obtained from the equal-area construction (see, for instance, Whitham), where the range of integration includes all those characteristics that will interact with the shock. This construction is valid so long as the driving pulse is finite and significantly shorter than a preceding quiet period. (For it to hold at y, the quiet period must last a duration yf(y).)

For our light, polytropic, power-law atmosphere, we consider ρ₀ ∝ r⁻ⁿ and adopt the Keplerian gravity of a central point mass, requiring ${p}_{0}\propto {r}^{-(n+1)}\propto {\rho }_{0}^{1+1/n}$, and c₀ ∝ r^−1/2. We adopt n = 5/2. Our grid spans a factor of 2 in radius (5.6 in density), which we resolve with 1.4 × 10⁵ cells. Into this we launch a train of sinusoidal waves ($v({r}_{a})={v}_{p}\sin \omega t$ for t > 0) with ω = 490r_a /c₀(r_a ), so the initial wavelength is r_a /78, or 1840 cells. The initial amplitude is v_p /c₀ = 1/100. At this amplitude and wavelength, the linear high-frequency limit holds very well, and the numerical resolution is sufficient. As we show in Figure 3, the simulation agrees very precisely with our analytical prediction for the radius of shock formation (R_s  = 1.15r_a ), as well as the development of the head shock and all subsequent internal shocks.

Zoom In Zoom Out Reset image size

The comparison is equally successful when we initialize a wave of only one sinusoidal cycle, which develops head and tail shocks but no internal shocks; see Figure 4. (We have explored other waveforms, such as sawtooth and reverse-sawtooth waves. These also show excellent agreement, in the linear regime, between analytical and numerical solutions.)

Zoom In Zoom Out Reset image size

In the bottom panels of Figures 3 and 5, we show the local maximum value of L_w in units of the input peak value L_p , demonstrating that wave luminosity is conserved in the acoustic regime and decays as expected during shock evolution.

Zoom In Zoom Out Reset image size

For a polytropic stellar model, we adopt the self-gravitating Lane-Emden solution of radius R_*, in which ${p}_{0}\propto {\rho }_{0}^{1+1/n};$ again, we choose the polytropic index n = 5/2. We arrange our grid to span radii 0.1R_*−0.999R_* with 1.5 × 10⁵ equally spaced cells. Over this range, ρ₀ varies by a factor of 2.6 × 10⁸ and ${L}_{\max }$ varies by a factor of 2.9 × 10¹¹. At the inner boundary, we introduce sinusoidal waves of low amplitude (v/c₀ = 10⁻²) and high frequency (R_*/200 in wavelength) so that each wave is initially resolved by approximately 800 cells. This simulation is displayed in Figure 5. Again, a head shock and multiple internal shocks form at the predicted location (R_s  = 0.147R_*) and evolve as predicted in the linear regime.

In contrast to the power-law atmosphere, however, the sharp drop in ${L}_{\max }$ near the stellar surface causes the head shock to become nonlinear. The departure from linear theory is visible by 0.98R_*, where flow downstream of the head shock reaches v/c₀ ∼ 0.5, and grows from there. Linear theory remains remarkably accurate for flow Mach numbers of order unity. New features, such as the presence of significant postshock velocity, become visible as the head shock becomes strong.

We note that the internal shocks remain weak as the head shock strengthens. Indeed, the internal shocks continue to decelerate (because V ≈ c₀ for a weak shock) as the head shock begins to accelerate; for this reason, a gap opens behind the advancing head shock. (Tail shocks differ from head shocks in the nonlinear regime: they cannot strengthen and advance to the same degree.) Although our simulation volume does not capture the region outside R_*, we can be confident that the head shock develops sufficient strength to eject matter from the stellar surface.

The most striking result of our analysis is the sharp difference in the decay properties of head and tail shocks, compared to internal shocks, in the simple model of a sinusoidal wave train composed of a finite number of cycles. Consider the decay phase in which $f({t}_{i,d})\to [\pi {\tilde{y}}^{-1},2{\tilde{y}}^{-1/2}]$ for [internal, head/tail] shocks, and let us define a characteristic local propagation time

$$\begin{eqnarray*}&&{t}_{\mathrm{char}}(r)\equiv Y(r){L}_{\max }^{1/2}\end{eqnarray*}$$

to eliminate Y. Note that t_char need not be monotonic in r. With this and our other definitions, the asymptotic shock strength is given by

$$\begin{eqnarray}\displaystyle \frac{{v}_{d}^{2}}{{c}_{0}^{2}}=\displaystyle \frac{{L}_{w}({t}_{i,d})}{{L}_{\max }}\to \left\{\begin{array}{ll}\tfrac{{\pi }^{2}}{{\omega }^{2}{t}_{\mathrm{char}}^{2}} & \text{internal shock}\\ \tfrac{4}{\omega {t}_{\mathrm{char}}}\tfrac{{L}_{p}^{1/2}}{{L}_{\max }{\left(r\right)}^{1/2}} & \text{head/tail shock}\end{array}\right..\end{eqnarray} \tag{ 24 }$$

Internal shocks, remarkably, do not depend in this limit on the initial wave amplitude; this implies a universal asymptotic shock-heating behavior that depends only on ω and local conditions, which we shall put to use in Section 6.2. Internal shocks dissipate wave energy relatively rapidly, but nevertheless strengthen in dimensionless terms in regions where t_char(r) declines outward.

Head and tail shocks, in contrast, retain some sensitivity to the input luminosity and strengthen once the combination ${t}_{\mathrm{char}}{(r){L}_{\max }(r)}^{1/2}$ begins to decline. Head and tail shocks evolve differently in the nonlinear regime; head shocks accelerate more readily to become nascent blast waves.

Modeling the input wave as a pure sinusoidal wave train is highly idealized, so we consider the internal shocks and head or tail shocks as limiting cases of what would be observed in a real star. We anticipate that a continuous but incoherent wave source will tend to resemble the internal shocks studied here, perhaps with a rare head- or tail-like shock arising by chance. Transient events, on the other hand, may tend to resemble the head and tail shocks but also produce internal shocks in their ring-down phases. We defer a full analysis to later work.

Fuller & Ro (2018) present a useful differential equation for the dissipation of wave energy at internal shocks, building on results from Landau & Lifshitz (1959), Ulmschneider (1970), and Mihalas & Mihalas (1984). We check this formula by examining the asymptotic behavior of internal shocks in sinusoidal waves: $\tilde{y}\to \pi /f$ for large $\tilde{y}$, from Equation (16). Differentiating this and using definitions for $\tilde{y}$ and f to express the result in terms of L_w , and then using A ρ₀ dr = dm (defining the mass coordinate m), we find

$$\begin{eqnarray}&&{{dL}}_{w}=-\displaystyle \frac{\gamma +1}{\pi }\omega {c}_{0}^{2}{\left(\displaystyle \frac{{L}_{w}}{{L}_{\max }}\right)}^{3/2}{dm}.\end{eqnarray} \tag{ 25 }$$

Note, L_w here means L_w (t_i,d )—that is, it refers to the conditions at the shock. The phase-averaged luminosity ${\bar{L}}_{w}$ is the more relevant quantity. Given that L_w  ∝ f², ${\bar{L}}_{w}={L}_{w}({t}_{i,d})/3$ in the asymptotic "N-wave" state; therefore,

$$\begin{eqnarray}&&d{\bar{L}}_{w}=-\displaystyle \frac{\gamma +1}{\pi }\sqrt{3}\omega {c}_{0}^{2}{\left(\displaystyle \frac{{\bar{L}}_{w}}{{L}_{\max }}\right)}^{3/2}{dm}.\end{eqnarray} \tag{ 26 }$$

The dissipation rate presented by Fuller & Ro (2018) agrees with this result, up to an inconsistency in their Equation (17): on its left-hand side, the symbol L_w refers to ${\bar{L}}_{w}$, while on its right-hand side L_w means ${L}_{w}({t}_{{id}})={\bar{L}}_{w}/3$. Note also that Equation (26) is accurate with our definition of ${L}_{\max }$, which is double the definition used by Fuller & Ro.

Expression (26) is convenient in form, as it can be augmented to account for nonshock dissipation, as Fuller & Ro have done. However, it has a couple of drawbacks. It applies only to the phase of asymptotic decay, so it cannot accurately describe the bulk of wave heating. It can only be applied in the region $\tilde{y}\gt 1$ where shocks exist; but even there, the initial value of ${\bar{L}}_{w}$ depends on what survives the early phase of shock evolution.

For these reasons, we recommend instead using the exact form of ${\bar{L}}_{w}(\tilde{y})$ that results from recognizing that ${\bar{L}}_{w}$ is proportional to the residual wave energy ${\tilde{E}}_{w}$ as shown in Figure 2:

$$\begin{eqnarray}&&\displaystyle \frac{{\bar{L}}_{w}(\tilde{y})}{{\bar{L}}_{w}(0)}=1-\displaystyle \frac{{\tilde{E}}_{\mathrm{diss}}({\tilde{t}}_{i}(\tilde{y}))}{\pi },\end{eqnarray} \tag{ 27 }$$

for $\tilde{y}\gt 1$, where ${\tilde{E}}_{\mathrm{diss}}({\tilde{t}}_{i})$ is given by Equation (18) and ${\tilde{t}}_{i}(\tilde{y})={\mathrm{sinc}}^{-1}(1/\tilde{y})$.

Or, for an even more convenient alternative, we suggest an approximation based on the common asymptotic evolution of internal shocks in Equation (24), using ${\bar{L}}_{w}={L}_{w}({t}_{i,d})/3$ :

$$\begin{eqnarray}&&{\bar{L}}_{w}(r)\simeq {\left[{\bar{L}}_{w}{\left(0\right)}^{-k}+{\left(\displaystyle \frac{{\pi }^{2}{L}_{\max }(r)}{3{\omega }^{2}{t}_{\mathrm{char}}{\left(r\right)}^{2}}\right)}^{-k}\right]}^{-1/k},\end{eqnarray} \tag{ 28 }$$

which respects both limits and remains within 10% for k = 0.64.

Both the exact and approximate expressions (Equations (27) and 28) assume the input wave is a pure uninterrupted sinusoid that remains in the linear, high-frequency regime, so its shock dissipation can be predicted by the results of Section 4.1.1. For a more sophisticated approach, one must start from a more realistic input waveform. These expressions also ignore other types of dissipation, such as radiation diffusion; we (Ro & Matzner 2017) have found this to be valid when ${\bar{L}}_{w}$ exceeds the radiation diffusion luminosity of the background state.

We wish to make a rough estimate for the mass ejected by a head shock that strengthens in the outer stellar envelope, in the weak-wave limit for which this occurs very close to the stellar surface (z ≡ (R_* − r₀)/R_* ≪ 1, where R_* is the stellar radius). We begin with ejection by a weak acoustic pulse, in which shock formation also occurs in the outer envelope, and then consider the case of a central point explosion in the next subsection.

Our analysis here is only asymptotic, however, as we ignore radiation diffusion, which would truncate the acceleration toward v_esc *, the stellar escape velocity, for sufficiently low wave energies (Linial et al. 2020). We also require that the wave be sufficiently radial that lateral pressure gradients do not alter the character of the flow (Matzner et al. 2013; Salbi et al. 2014; Linial & Sari 2019).

The outer polytropic scalings are $A=4\pi {R}_{* }^{2}$, ρ₀ = ρ_0h zⁿ , and c₀ = c_0h z^1/2, where the "h" subscript indicates the outer normalization. From these, we derive the external mass m_h zⁿ⁺¹ (where ${m}_{h}=4\pi {R}_{* }^{3}{\rho }_{0h}/(n+1)$), maximum luminosity ${L}_{\max }\,={L}_{\max ,h}{z}^{n+3/2}$ (where ${L}_{\max ,h}=4\pi {R}_{* }^{2}{\rho }_{0,h}{c}_{0,h}^{3}$), and characteristic time t_char = t_char,h z^1/2 (where t_char,h  = (γ + 1)(2n + 1)R_*/8c_0,h )).

For a reasonable approximation to the ejected mass, we first identify the value of z = z_S at which the shock becomes strong, in the sense that v_d (z_S ) = c₀(z_S ), using Equation (24) to extrapolate the head shock out of the weak-shock regime. We determine the corresponding shock velocity ${v}_{s}({z}_{S})\,=\left(\gamma +1+\sqrt{17+\gamma (2+\gamma })\right){c}_{0}({z}_{S})/4$. We then apply self-similar scalings for a strong planar shock, ${v}_{s}\propto {\rho }_{0}^{-\beta }$ (Sakurai 1960) until the ejection condition f_sph Cv_s  = v_esc * is met. Here, β ≃ 0.19 is the shock acceleration index, and C ∼ 2 is the postshock acceleration factor for purely planar, nonrelativistic flow (Sakurai 1960; Ro & Matzner 2013). The factor f_sph accounts for the difference between spherical and planar flow, which matters for finite but small E_in because of the slow nature of postshock acceleration (Litvinova & Nadezhin 1990); for nongravitating flow, f_sph ≃ (1 − {0.51, 0.34}z^1/3) for $n=\left\{\tfrac{3}{2},3\right\}$ (Matzner & McKee 1999).

We must also allow for the possibility of wave reflection, in the case where the high-frequency limit does not hold throughout its transition to a strong, accelerating shock. (Reflection is not an issue in the strong regime, as Sakurai's solution accounts for it self-consistently.) The high-frequency limit requires ω zR_*/c₀ ≫ 1;  evaluating this at the weak-to-strong transition,

$$\begin{eqnarray}&&\omega \displaystyle \frac{{z}_{S}{R}_{* }}{{c}_{0}({z}_{S})}={\left[\displaystyle \frac{{L}_{p}}{{L}_{\max ,h}}{\left(\displaystyle \frac{{t}_{\mathrm{char},h}{c}_{0h}}{4{R}_{* }}\right)}^{2}{\left(\displaystyle \frac{{c}_{0h}}{\omega {R}_{* }}\right)}^{3+2n}\right]}^{\tfrac{1}{5+2n}}.\end{eqnarray} \tag{ 29 }$$

As an analysis of wave reflection is beyond the scope of this paper, we capture it using a velocity transmission prefactor f_tm ≤ 1, where f_tm = 1 in the absence of reflection. Then,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{m}_{\mathrm{ej}}}{{m}_{h}} & = & {\left[{f}_{\mathrm{tm}}{f}_{\mathrm{sph}}K\displaystyle \frac{{c}_{0h}}{{v}_{\mathrm{esc}* }}{\left(\displaystyle \frac{{L}_{p}}{{L}_{\max ,h}}\displaystyle \frac{16}{{\omega }^{2}{t}_{\mathrm{char},{\rm{h}}}^{2}}\right)}^{\displaystyle \frac{1+2\beta n}{5+2n}}\right]}^{{\gamma }_{p}/\beta }\\ & \longrightarrow & {\left(\displaystyle \frac{3.47{f}_{\mathrm{tm}}{f}_{\mathrm{sph}}{c}_{0h}}{{v}_{\mathrm{esc}* }}\right)}^{7.18}{\left(\displaystyle \frac{{L}_{p}}{{L}_{\max ,h}}\displaystyle \frac{1}{{\omega }^{2}{t}_{\mathrm{char},h}^{2}}\right)}^{1.38}\end{array}\end{eqnarray} \tag{ 30 }$$

The coefficient K ≃ Cv_s (z_S )/c₀(z_S ) depends on our approximate matching between weak and strong limits. The evaluation on the second line of (30) is for radiation-dominated flow in an n = 3 polytrope (γ = γ_p  = 4/3).

Our results can also be applied to the similar problem of mass ejection caused by an explosion of small but finite strength in the stellar center, which has been studied by Wyman et al. (2004) and Linial et al. (2020) for polytropes, and by Dessart et al. (2010), Owocki et al. (2019), and Kuriyama & Shigeyama (2020) for realistic stellar models. We consider polytropes for simplicity. The explosion is strong until it decelerates to about c₀; for a complete polytrope, this occurs on a time ${t}_{\mathrm{dec}}\simeq {[2/{(5{c}_{0c})]}^{5/3}({E}_{\mathrm{in}}/{\rho }_{0c})}^{1/3}$ for an explosion of energy E_in that is small compared to the binding energy; the subscript c indicates conditions at the center, and "dec" refers to the conditions at the deceleration radius.

While there is no shock-free phase in this problem, we can, for a good approximation, adopt a head shock model in which the shock forms at the deceleration radius (r_s  ≃ r_dec); the energy in the wave pulse is  ∼ E_in/2 (accounting for heat deposited within r_dec); and the pulse duration, or wave half-period, is π/ω ≃ t_dec. Correspondingly, L_p  ≃ E_in/t_dec.

We require the normalized coordinate $\tilde{y}(r)=Y(r)/{Y}_{s}\,=Y(r)/Y({r}_{\mathrm{dec}})$. For the denominator, we note that Y(r) is almost constant in a uniform medium: $Y({r}_{\mathrm{dec}})={\rm{\Lambda }}(\gamma +1){(4\pi {c}_{0c}{\rho }_{0c}^{5})}^{1/2}/2;$ we drop Λ, which is a logarithm of r_dec. For the numerator, $Y{(z)=2[(\gamma +1)/(2n+1)]{{zR}}_{* }/[{c}_{0}(z){L}_{\max }(z)}^{1/2}]$ in the subsurface region (z ≪ 1).

We now apply the asymptotic head shock formula ${v}_{d}^{2}/{c}_{0}^{2}=(4/\tilde{y}){L}_{p}/{L}_{\max }$ to identify z_S as the depth at which the shock regains its strength (v_d  = c₀) and proceed to estimate m_ej as before. The full expression is too cumbersome to write, but we note that the ejected mass scales $\propto {({f}_{\mathrm{tm}}{f}_{\mathrm{sph}})}^{{\gamma }_{p}/\beta }{E}_{\mathrm{in}}^{\alpha }$, where

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{4{\gamma }_{p}(1+2\beta n)}{3\beta (5+2n)}.\end{eqnarray} \tag{ 31 }$$

The coefficient can be determined from the polytropic structure, which fixes ${({\rho }_{0h}/{\rho }_{0c})}^{1/n}={({c}_{0h}/{c}_{0c})}^{2}={(\xi \theta ^{\prime} )}_{{\xi }_{n}},$ where θ(ξ) is the Lane-Emden function with outer coordinate ξ_n . For the case n = 3/2 and γ = 5/3, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{\mathrm{ej}}}{{M}_{* }}\simeq {\left(2.23{f}_{\mathrm{tm}}{f}_{\mathrm{sph}}\right)}^{7.52}{\left(\displaystyle \frac{{E}_{\mathrm{in}}}{{E}_{\mathrm{bind}}}\right)}^{2.09},\end{eqnarray} \tag{ 32 }$$

where E_bind is the negative total energy of the polytrope. This is consistent with the results of Linial et al. (2020) if Kf_tm f_sph increases slowly with E_in/E_bind.

We suspect that the threshold energy observed by Wyman et al. (2004) is a consequence of finite numerical resolution in that work. However, radiation diffusion sets a true lower limit for m_ej. For an analysis of this point, see Linial et al. (2020).

A "failed supernova" is a collapsing massive star in which the stalled protoneutron star bounce shock fails to revive and eject the envelope. Instead, the protoneutron star is expected to collapse into a black hole because the total accreted mass will exceed the Tolman–Oppenheimer–Volkoff limit. The nearly instantaneous liberation of ∼0.3 M_⊙ in neutrinos from the formation of the protoneutron star causes the overpressurized envelope to accelerate outward at all radii, at least initially. The expansion launches an outward-propagating spherical acoustic pulse that steepens into a weak shock. In the wake of the shock is an expanding rarefaction wave that informs the exterior envelope of a collapse in the stellar core.

Red and yellow supergiant stars are well approximated by the light polytropic envelope structure described in Section 5.1, with index n = 2.5, over several decades of r. For all positive values of n, there are at least two self-similar solutions to the flow structure: one in which the weak shock becomes strong and is described by the Sedov–Taylor solutions, and another where its strength vanishes and it merges with the leading edge of the rarefaction wave. For the range 2 < n < 3.5, Coughlin et al. (2018) discovered a third set of intermediate self-similar solutions involving weak shocks of a finite strength. However, Coughlin et al. (2019) show by linear analysis that these self-similar weak-shock solutions are linearly unstable to perturbations in the shock strength, and Ro et al. (2019) show that these shocks asymptotically approach the rarefaction and Sedov–Taylor self-similar solutions.

Why do these intermediate solutions only exist for 2 < n < 3.5? The lower limit presumably derives from the fact that, for n < 2, the Sedov–Taylor blast wave weakens over time, as it implies V/c₀ ∝ t^{(n−2)/(5−n)}. For n ≥ 3.5, Ro et al. (2019) found that all shocks strengthen, evolving toward the Sedov–Taylor limit. Using our solutions from Section 4, we find that a head shock obeys ${v}_{d}/{c}_{0}\propto r\sqrt{A{\rho }_{0}{c}_{0}}\propto {r}^{(n-7/2)/4}$: it strengthens over time when n > 3.5. In these conditions, a weak shock tends to outrun outward characteristics, causally disconnecting itself from the collapsing interior.