Energy-conserving, Relativistic Corrections to Strong Shock Propagation

TL;DR

This paper derives energy-conserving relativistic corrections to the classical Sedov-Taylor shock solution, applicable to mildly relativistic astrophysical explosions, with implications for supernovae and gamma-ray bursts.

Contribution

It introduces the first perturbative, energy-conserving relativistic corrections to the Sedov-Taylor solution for strong shocks in astrophysics.

Findings

Relativistic terms modify post-shock fluid properties and shock speed.

No relativistic correction satisfies energy conservation for certain blastwave solutions.

Corrections are derived for various adiabatic indices and ambient density profiles.

Abstract

Astrophysical explosions are accompanied by the propagation of a shock wave through an ambient medium. Depending on the mass and energy involved in the explosion, the shock velocity $V$ can be non-relativistic ($V \ll c$, where $c$ is the speed of light), ultra-relativistic ($V \simeq c$), or moderately relativistic ($V \sim few\times 0.1c$). While self-similar, energy-conserving solutions to the fluid equations that describe the shock propagation are known in the non-relativistic (the Sedov-Taylor blastwave) and ultra-relativistic (the Blandford-McKee blastwave) regimes, the finite speed of light violates scale invariance and self-similarity when the flow is only mildly relativistic. By treating relativistic terms as perturbations to the fluid equations, here we derive the $\mathcal{O}(V^2/c^2)$, energy-conserving corrections to the non-relativistic, Sedov-Taylor solution for the propagation of a strong shock. We show that relativistic terms modify the post-shock fluid velocity, density, pressure, and the shock speed itself, the latter being constrained by global energy conservation. We derive these corrections for a range of post-shock adiabatic indices $\gamma$ (which we set as a fixed number for the post-shock gas) and ambient power-law indices $n$, where the density of the ambient medium $\rho_{\rm a}$ into which the shock advances declines with spherical radius $r$ as $\rho_{\rm a} \propto r^{-n}$. For Sedov-Taylor blastwaves that terminate in a contact discontinuity with diverging density, we find that there is no relativistic correction to the Sedov-Taylor solution that simultaneously satisfies the fluid equations and conserves energy. These solutions have implications for relativistic supernovae, the transition from ultra- to sub-relativistic velocities in gamma-ray bursts, and other high-energy phenomena.

Tables

Table 1. Table 1: The eigenvalue σ 𝜎 \sigma that ensures energy conservation as a function of n 𝑛 n , which parameterizes the density profile of the ambient medium ( ρ ′ ∝ r − n proportional-to superscript 𝜌 ′ superscript 𝑟 𝑛 \rho^{\prime}\propto r^{-n} ), and the post-shock adiabatic index of the gas γ 𝛾 \gamma . Values much greater than one indicate that relativistic effects are more important for less relativistic initial shock speeds, which are achieved for small n 𝑛 n and γ 𝛾 \gamma . Cells with an ellipsis correspond to instances where the Sedov-Taylor (non-relativistic) solution possesses a finite or diverging density at the contact discontinuity, for which we find no energy-conserving solution for the relativistic corrections.

	1.05	1.10	1.15	1.20	1.25	1.30	4/3	1.40	1.45	1.50	1.55	1.60	5/3	1.70
0	$\sigma =$ 6.05	3.08	2.10	1.63	1.35	1.17	1.08	0.950	0.882	0.831	0.790	0.758	0.724	0.711
0.5	5.04	2.56	1.75	1.36	1.13	0.982	0.910	0.807	0.753	0.712	0.681	0.656	0.631	0.621
1	4.02	2.05	1.40	1.09	0.911	0.796	0.741	0.663	0.623	0.593	0.571	0.554	0.537	0.531
1.5	3.00	1.53	1.05	0.822	0.692	0.610	0.572	0.519	0.493	0.474	0.461	0.452	0.444	0.442
2	1.98	1.01	0.698	0.553	0.474	0.427	0.405	0.378	0.367	0.360	0.357	0.357	0.360	0.363
2.5	0.962	0.492	0.354	0.298	0.272	0.264	0.265	…	…	…	…	…	…	…
2.7	0.559	0.299	0.235	0.224	…	…	…	…	…	…	…	…	…	…

Equations

$$V\simeq \sqrt{\frac{E_{\mathrm{ej}}}{{\rho }_{\mathrm{a}}{r}_{\mathrm{a}}^3}}{\left(\frac{R}{r_{\mathrm{a}}}\right)}^{\frac{n-3}{2}}\simeq V_{\mathrm{ej}}{\left(\frac{R}{r_{\mathrm{a}}}\right)}^{\frac{n-3}{2}}.$$

$$R\propto t^{\frac{2}{5-n}}.$$

$$\Gamma \simeq \sqrt{\frac{E_{\mathrm{ej}}}{{\rho }_{\mathrm{a}}{r}_{\mathrm{a}}^3{c}^2}}{\left(\frac{R}{r_{\mathrm{a}}}\right)}^{(n-3)/2}\simeq \frac{V_{\mathrm{ej}}}{c}{\left(\frac{R}{r_{\mathrm{a}}}\right)}^{(n-3)/2}.$$

$$\Gamma \simeq \frac{V_{\mathrm{ej}}}{c}{\left(\frac{ct}{r_{\mathrm{a}}}\right)}^{\left(n-3\right)/2}.$$

$$\Delta t_{\mathrm{rel}}\simeq \frac{r_{\mathrm{a}}}{c}{V}_{\mathrm{ej}}^{\frac{2}{3-n}}\simeq \frac{r_{\mathrm{a}}}{c}{\left(\frac{E_{\mathrm{ej}}}{M_{\mathrm{ej}}{c}^2}\right)}^{\frac{1}{3-n}}.$$

$$V^2\simeq V_{\mathrm{ej}}^2{\left(\frac{R}{r_{\mathrm{a}}}\right)}^{n-3}\left(1+\sigma \frac{V^2}{c^2}\right),$$

$$\frac{V^2}{c^2}\simeq \frac{V_{\mathrm{ej}}^2}{c^2}{\left(\frac{R}{r_{\mathrm{a}}}\right)}^{n-3}\propto t^{\frac{2(n-3)}{5-n}},$$

$${\nabla }_{\mu }{T}^{\mu \nu }=0,$$

$$T^{\mu \nu }=\left({\rho }^{\prime }+\frac{\gamma }{\gamma -1}{p}^{\prime }\right)U^{\mu }{U}^{\nu }+p^{\prime }{g}^{\mu \nu }$$

$${\nabla }_{\mu }\left[{\rho }^{\prime }{U}^{\mu }\right]=0,$$

$$U_{\mu }{U}^{\mu }=-1.$$

$$U^{\mu }{\nabla }_{\mu }{s}^{\prime }=0,$$

$${\rho }^{\prime }{U}^{\mu }{\nabla }_{\mu }{U}^{\nu }+\frac{\gamma }{\gamma -1}{p}^{\prime }{U}^{\mu }{\nabla }_{\mu }{U}^{\nu }+{\Pi }^{\mu \nu }{\nabla }_{\mu }{p}^{\prime }=0.$$

$$U^t=\sqrt{1+U^2},$$

$$\frac{\partial }{\partial t}\left[{\rho }^{\prime }\sqrt{1+U^2}\right]+\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2{\rho }^{\prime }U\right]=0,$$

$$\sqrt{1+U^2}\frac{\partial s^{\prime }}{\partial t}+U\frac{\partial s^{\prime }}{\partial r}=0,$$

$$\sqrt{1+U^2}\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial r}+\frac{p^{\prime }U}{{\rho }^{\prime }}\left\{\sqrt{1+U^2}\frac{\partial }{\partial t}\ln \left(p^{\prime }{U}^{\frac{\gamma }{\gamma -1}}\right)+U\frac{\partial }{\partial r}\ln \left(p^{\prime }{U}^{\frac{\gamma }{\gamma -1}}\right)\right\}+\frac{1}{{\rho }^{\prime }}\frac{\partial p^{\prime }}{\partial r}=0.$$

$$\frac{\partial }{\partial t}\left[\left({\rho }^{\prime }+\frac{\gamma }{\gamma -1}{p}^{\prime }\right)\left(1+U^2\right)-p^{\prime }\right]+\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\left({\rho }^{\prime }+\frac{\gamma }{\gamma -1}{p}^{\prime }\right)U\sqrt{1+U^2}\right]=0,$$

$$\frac{\partial }{\partial t}\left[\left({\rho }^{\prime }+\frac{\gamma }{\gamma -1}{p}^{\prime }\right)U\sqrt{1+U^2}\right]+\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\left(\left({\rho }^{\prime }+\frac{\gamma }{\gamma -1}{p}^{\prime }\right)U^2+p^{\prime }\right)\right]-\frac{2p^{\prime }}{r}=0.$$

$${\rho }_2^{\prime }{U}_2^{\prime \prime }={\rho }_{\mathrm{a}}^{\prime }{U}_{\mathrm{a}}^{\prime \prime },$$

$$\left({\rho }_2^{\prime }+\frac{\gamma }{\gamma -1}{p}_2^{\prime }\right)U_2^{\prime \prime }\sqrt{1+{\left(U_2^{\prime \prime }\right)}^2}=\left({\rho }_{\mathrm{a}}^{\prime }+\frac{\gamma }{\gamma -1}{p}_{\mathrm{a}}^{\prime }\right)U_{\mathrm{a}}^{\prime \prime }\sqrt{1+{\left(U_{\mathrm{a}}^{\prime \prime }\right)}^2}$$

$$\left({\rho }_2^{\prime }+\frac{\gamma }{\gamma -1}{p}_2^{\prime }\right){\left(U_2^{\prime \prime }\right)}^2+p_2^{\prime }=\left({\rho }_{\mathrm{a}}^{\prime }+\frac{\gamma }{\gamma -1}{p}_{\mathrm{a}}^{\prime }\right){\left(U_{\mathrm{a}}^{\prime \prime }\right)}^2+p_{\mathrm{a}}^{\prime }.$$

$${\gamma }^2\left(1+{\left(U_2^{\prime \prime }\right)}^2\right)\left(U_2^{\prime \prime }-U_{\mathrm{a}}^{\prime \prime }\right)+2\gamma U_{\mathrm{a}}^{\prime \prime }\left(1+{\left(U_2^{\prime \prime }\right)}^2\right)-U_2^{\prime \prime }-U_{\mathrm{a}}^{\prime \prime }=0.$$

$$U_2=U_{\mathrm{s}}\sqrt{1+{\left(U_2^{\prime \prime }\right)}^2}+U_2^{\prime \prime }\sqrt{1+U_{\mathrm{s}}^2}.$$

$${\rho }_2^{\prime }=-\frac{U_{\mathrm{s}}}{U_2^{\prime \prime }}{\rho }_{\mathrm{a}}^{\prime },$$

$$p_2^{\prime }=-\frac{\gamma -1}{\gamma }\frac{U_{\mathrm{s}}}{U_2^{\prime \prime }}\left(\frac{\sqrt{1+U_{\mathrm{s}}^2}}{\sqrt{1+{\left(U_2^{\prime \prime }\right)}^2}}-1\right){\rho }_{\mathrm{a}}^{\prime }.$$

$$U_2(R)=\frac{2}{\gamma +1}\left(1-\frac{\gamma -1}{2}\frac{{\gamma }^2+1}{{\left(\gamma +1\right)}^3}{U}_{\mathrm{s}}^2\right)U_{\mathrm{s}},$$

$${\rho }_2^{\prime }(R)=\frac{\gamma +1}{\gamma -1}\left(1+\frac{2\gamma }{{\left(\gamma +1\right)}^3}{U}_{\mathrm{s}}^2\right){\rho }_{\mathrm{a}}^{\prime }(R),$$

$$p_2^{\prime }(R)=\frac{2}{\gamma +1}\left(1-\frac{{\gamma }^2\left(\gamma -1\right)}{{\left(\gamma +1\right)}^3}{U}_{\mathrm{s}}^2\right){\rho }_{\mathrm{a}}^{\prime }(R)U_{\mathrm{s}}^2,$$

$$E_{\mathrm{tot}}=4\pi {\int }_{R_{\mathrm{c}}}^R{T}^{00}{r}^{2}dr=4\pi {\int }_{R_{\mathrm{c}}}^R\left(\left({\rho }^{\prime }+\frac{\gamma }{\gamma -1}{p}^{\prime }\right)\left(1+U^2\right)-p^{\prime }\right)r^{2}dr.$$

Full text

gamma-ray burst: general — hydrodynamics — methods: analytical — relativistic processes — shock waves — supernovae: general

1 Introduction

A core-collapse supernova is, in the now-classic picture, initiated by the “bounce” of the overpressured, protoneutron star that forms from the collapse of the iron core of a massive star (Colgate & White, 1966). From the protoneutron star bounce is launched a shock wave, which propagates through and unbinds the stellar envelope to yield the supernova. If it is sufficiently energetic, the shock promptly plows through the star following the bounce and liberates the gas; alternatively, if it is not energetic enough to overcome the ram pressure of the infalling material and the dissociation of heavy nuclei in the core of the star (Arnett, 1982), the shock “stalls” at small radii but can be revived by some means (e.g., neutrino heating, convective instabilities behind the shock, a dynamic instability of the standing shock, or the magnetorotational amplification of magnetic fields; respectively, e.g., Bethe & Wilson 1985; Burrows et al. 1995; Blondin et al. 2003; Mösta et al. 2015). If the shock fails to be revived, with sufficient angular momentum an accretion disc forms around the natal black hole that can, through the combination of bipolar outflows and liberated accretion energy, unbind the remaining stellar envelope in the collapsar picture of a long gamma-ray burst (Woosley, 1993; MacFadyen & Woosley, 1999; Woosley & Bloom, 2006). Finally, even in the absence of sufficient angular momentum, a failed supernova generates a secondary, weak shock in the outer layers of the star from the mass lost to neutrinos during the de-leptonization of the core (Nadezhin, 1980; Lovegrove & Woosley, 2013; Piro, 2013; Coughlin et al., 2018a; Fernández et al., 2018; Coughlin et al., 2018b).

All of these explosion scenarios involve the formation and expansion of a shock wave into its surroundings, and this shock leaves in its wake a “sea” of post-shock fluid (as, of course, do explosion scenarios not initiated by the collapse of a massive star, such as compact object mergers; e.g., Li & Paczyński 1998; Levinson et al. 2002; Nakar & Piran 2011; Abbott et al. 2017). One of the most useful techniques for describing the spatial and temporal evolution of the post-shock gas and of the shock itself is self-similarity. This mathematical technique exploits the scale invariance of the fluid equations and, in the absence of any temporal or spatial scales of the ambient medium, the necessary scale invariance of the solutions to those equations (e.g., Ostriker & McKee 1988).

Among the best-known examples of a self-similar solution to the fluid equations is the Sedov-Taylor (ST) blastwave (Sedov 1959; Taylor 1950). The ST blastwave describes the propagation of an energy-conserving, strong (Mach number much greater than one) shock into an ambient medium that possesses a power-law density profile. The conservation of energy implies that there is a unique shock speed $V$ that can be directly related to the initial energy of the explosion, the impulsive injection of which initiated the explosion in the first place. The ST solution is also non-relativistic, in that the shock speed is assumed to be much less than the speed of light and the energy is Newtonian, and hence terms of order $V^{2}/c^{2}$ that enter into the relativistic fluid equations are ignored ($c$ being the speed of light). The ST blastwave can be used to describe terrestrial explosions and can also constrain the age of supernova remnants (Chevalier, 1976).

In the other, extreme limit of an ultra-relativistic explosion – where the shock velocity is nearly equal to the speed of light – Blandford & McKee (1976) derived a distinct, energy-conserving, self-similar solution to the relativistic fluid equations. The Blandford-McKee (BMK) solution is the ultra-relativistic analog of the Sedov-Taylor blastwave, in that the conservation of energy implies that there is a unique shock Lorentz factor $\Gamma=(1-V^{2}/c^{2})^{-1/2}$ that is relatable to the explosion energy. The BMK blastwave is ultra-relativistic in the sense that the solution only accounts for terms in the fluid equations to order $\mathcal{O}(1/\Gamma^{2})$. While there are at present no (known) terrestrial applications of this solution, gamma-ray bursts should exhibit some phase of shock propagation appropriate to ultra-relativistic speeds, and simulations have found evidence of the transition to a Blandford-McKee-type phase of explosion (Kobayashi & Sari, 2000; Duffell & MacFadyen, 2013; Xie et al., 2018).

In this paper we are interested in the behavior of an energy-conserving explosion that is between the limits of Newtonian and ultra-relativistic. The specific question we ask is: when the flow is mildly relativistic, such that the shock speed is $V\sim few\times 0.1c$, how do relativistic effects modify the Sedov-Taylor solution? There are particularly energetic supernovae and lower-luminosity gamma-ray bursts that can produce marginally relativistic shock speeds (Soderberg et al., 2006; Drout et al., 2011; Corsi et al., 2014, 2016; Whitesides et al., 2017). For these modestly-relativistic flows, how do relativistic corrections modify the shock velocity scaling predicted from the ST blastwave, and how are the post-shock velocity, density, and pressure profiles altered from the self-similar solutions provided by the ST solution?

In this paper, we answer these questions by considering special relativistic terms as perturbations to the non-relativistic fluid equations and, correspondingly, the solutions to those equations. In Section 2, we first give some basic considerations of the problem, and we derive order of magnitude estimates for the corrections to the shock velocity that are induced by relativistic motion. In Section 3 we present a rigorous perturbation analysis of the shock jump conditions and the fluid equations to leading relativistic order, and from those equations we derive relativistic, non-self-similar corrections to the post-shock velocity, density, and pressure that result from the velocity scale established by the finite speed of light. We also demonstrate that there is a unique, relativistic correction to the shock velocity that results from the requirement that the energy – which includes relativistic terms – be exactly conserved behind the shock front. We summarize and conclude in Section 5.

2 General Considerations and Order of Magnitude Estimates

We characterize an explosion by the ejection of an amount of mass $M_{\rm ej}$ with a corresponding energy $E_{\rm ej}$, which can be combined to yield a characteristic velocity $V_{\rm ej}=\sqrt{E_{\rm ej}/M_{\rm ej}}$. As this material encounters an ambient medium medium, a forward shock is generated that initially propagates at the same characteristic speed $V_{\rm ej}$. Once the shock entrains sufficient inertia from its surroundings that the initial mass of the explosion is forgotten, the forward shock propagation settles into a self-similar state such that the boundary conditions at the shock (namely the jump conditions; see Section 3.3) govern the entire post-shock evolution of the flow, while the characteristic ambient density $\rho_{\rm a}$ and length scale $r_{\rm a}$ and the shock energy $E_{\rm ej}$ dictate the propagation of the shock itself. If the density of the ambient medium falls off as a power-law with radial power-law index $n$, then, during this self-similar phase, energy conservation in the non-relativistic limit implies that the shock velocity $V$ is related to the shock position $R$ via

[TABLE]

The last line in this expression follows from the fact that the self-similar state is only reached once the mass entrained from the ambient medium is comparable to the initial mass of the explosion, i.e., $\rho_{\rm a}r_{\rm a}^{3}\simeq M_{\rm ej}$. This expression can be rearranged and integrated to solve for the shock position as a function of time, which yields

[TABLE]

This temporal scaling of the shock position is the well-known, Newtonian-energy-conserving result derived independently by Sedov (1959) and Taylor (1950)

We see from Equation (1) that when the mass involved in the explosion is small or the energy imparted to that mass is large, both of which are possible in astrophysical contexts, the characteristic ejecta velocity can exceed the speed of light. In this limit the Newtonian approach breaks down, and instead of being constrained by a characteristic velocity, the shock can be parameterized by its Lorentz factor $\Gamma=(1-V^{2}/c^{2})^{-1/2}$. In the ultra-relativistic limit, the post-shock inertia is dominated by the internal energy of the gas, which scales as $\rho_{\rm a}\Gamma^{2}$, and mass conservation implies that the post-shock gas is swept into a thin shell of width $\sim R/\Gamma^{2}$ (Blandford & McKee, 1976). To leading order in the shock Lorentz factor, relativistic energy conservation then dictates that the Lorentz factor of the shock satisfies

[TABLE]

In this expression $R\simeq c\,t$, and it therefore follows that the Lorentz factor of the relativistic, self-similar shock evolves temporally as

[TABLE]

In between these two limits – Sedov-Taylor when $\sqrt{E_{\rm ej}/M_{\rm ej}}\ll c$ and Blandford-McKee when $\sqrt{E_{\rm ej}/M_{\rm ej}}\gg c$ – how does the flow behave? In addition to the examples of hyperenergetic supernovae considered in Section 1, it can also be seen from Equation (4) that – even if the blastwave is extremely energetic and begins in the Blandford-McKee regime – the deceleration of the shock Lorentz factor (assuming $n<3$, above which the Sedov-Taylor solution does not exist; we will always adopt $n<3$ in this paper, and return to the case of $n>3$ briefly in the conclusions) implies that the ultra-relativistic approximation only holds for a finite time $\Delta t_{\rm rel}$. Setting $\Gamma=1$ in Equation (4) and rearranging shows that this timescale is

[TABLE]

For fiducial values of $E_{\rm ej}\simeq 10^{52}$ erg, $M_{\rm ej}\simeq 0.1M_{\odot}$, $n=0$, and setting $r_{\rm a}\simeq N\times R_{\odot}$, we find $\Delta t_{\rm rel}\simeq 1\times N$ sec. Therefore, even for large values of $r_{\rm a}$ (or correspondingly small values of the ambient density), the ultra-relativistic phase of self-similar shock propagation can be very short lived. Relativistic terms that modify the solution from the non-relativistic, Sedov-Taylor phase – to which the flow eventually asymptotes – will then be present during the transition from ultra- to non-relativistic shock expansion.

At the order of magnitude level, the importance of relativistic corrections to the Sedov-Taylor blastwave can be understood by noting that the leading-order, relativistic modifications to the fluid equations appear as $\propto\mathcal{O}(v^{2}/c^{2})$, where $v$ is the three-velocity of the fluid; we derive these corrections explicitly in Section 3, but such a scaling is reasonable from the observation that the four-velocity, which differs from the three velocity by a factor of $\Gamma\simeq 1+\mathcal{O}(v^{2}/c^{2})$, transforms in a covariant (i.e., tensor-like) sense and therefore enters manifestly into the relativistic fluid equations. There will therefore be corrections of this same order to the jump conditions at the shock front, and hence the lowest-order, relativistic corrections to the post-shock velocity, density, and pressure will be of the order $V^{2}/c^{2}$. There will also be modifications to the conserved energy that enter as $v^{2}/c^{2}$ (again, we show this explicitly below, but this feature follows naturally from the covariant nature of the four-velocity over the three-velocity). Therefore, in order to satisfy energy conservation, there must also be relativistic corrections of the order $V^{2}/c^{2}$ that modify the propagation of the shock itself. We thus expect that during the marginally-relativistic phase when $V_{\rm ej}/c\lesssim 1$, the shock speed will be characterized by

[TABLE]

where $\sigma$ is an unknown but otherwise pure number. Since the shock is assumed to be only mildly relativistic, it follows that this relativistic term can be approximated as

[TABLE]

where the last line follows from the expression for $R(t)$ from the Sedov-Taylor solution (Equation 2). Thus, for $n=0$, the relativistic corrections fall off as $\propto t^{-6/5}$, and become – as one would anticipate – less important as the shock decelerates owing to the entrainment of mass from the ambient medium. Note, however, that this power-law decline in the importance of relativistic effects is shallower than that predicted from the Blandford-McKee solution alone, being $\Gamma\propto t^{-3/2}$ for $n=0$ (Equation 4); relativistic corrections to non-relativistic shock propagation can therefore be longer lived than might be anticipated by extrapolating the Blandford-McKee solution to the limit of $\Gamma=1$. Furthermore, as the density of the ambient medium falls off more steeply, relativistic effects remain important for longer periods of time (indeed, in the limit that $n\rightarrow 3$, the relativistic corrections are permanent modifications to the shock speed).

While these order of magnitude estimates provide useful diagnostics for probing the importance of relativistic effects and the rate at which they should depreciate in time, they cannot be used, for example, to determine the constant $\sigma$ that enters into the correction for the shock velocity in Equation (6). Furthermore, while these simple estimates indicate that the corrections to the post-shock fluid variables (the velocity, density, and pressure) enter at the level $V^{2}/c^{2}$, they tell us nothing about the spatial dependence of these corrections. To understand these aspects of the problem, we now turn to a quantitative analysis of the relativistic fluid equations and treat the leading-order, relativistic corrections in the perturbative limit.

3 Fluid Equations and Perturbative Approach

3.1 General continuity equations

The equations of hydrodynamics in covariant and differential form are given by

[TABLE]

where $\nabla_{\mu}$ is the covariant derivative and

[TABLE]

is the energy-momentum tensor of the fluid; here primes denote quantities measured in the comoving frame (i.e., in the frame where the fluid is instantaneously at rest), $\rho^{\prime}$ is the fluid density, $p^{\prime}$ is the gas pressure, $U^{\mu}$ is the fluid four-velocity, and $g^{\mu\nu}$ is the inverse of the metric, and we also adopted an adiabatic equation of state such that the internal energy $e^{\prime}$ is related to the pressure $p^{\prime}$ via $e^{\prime}=p^{\prime}/(\gamma-1)$. In this equation and for the remainder of this paper we adopt the Einstein summation convention, meaning that repeated upper and lower indices imply summation over the repeated index, and we will also use units where the speed of light is equal to one. The continuity equation, which ensures the conservation of mass, is written in covariant form as

[TABLE]

and to maintain the invariance of the line element the contraction of the four velocity with itself must be conserved:

[TABLE]

Equations (8) – (11) constitute the inviscid conservation laws that govern the evolution of the components of the four velocity, the mass density, and the pressure of the fluid in any arbitrary geometry.

While one can deal directly with the individual components of Equation (8), two additional, useful equations are obtained by contracting this expression with $U_{\nu}$ and $\Pi^{\beta}_{\,\,\nu}=U^{\beta}U_{\nu}+g^{\beta}_{\,\,\nu}$, which respectively select out the time-like and space-like components of the equations (note that $U_{\beta}\Pi^{\beta}_{\,\,\nu}=0$). Doing so and performing some simple manipulations gives the entropy equation

[TABLE]

where $s^{\prime}=\ln\left[p^{\prime}/(\rho^{\prime})^{\gamma}\right]$ is the entropy, and the momentum equations

[TABLE]

3.2 Equations in spherical symmetry

Here we restrict solutions to the fluid equations to be spherically symmetric and irrotational, such that the metric is given by $g_{\mu\nu}=\text{diag}\left\{-1,1,r^{2},r^{2}\sin^{2}\theta\right\}$, where $r$ is radial distance from the origin and $\theta$ is the standard polar angle in spherical coordinates, and the only non-vanishing components of the four-velocity are $U^{r}\equiv U$ and $U^{t}$. From Equation (10) we can relate the time component of the four-velocity to the radial component via

[TABLE]

and using this expression the continuity, entropy, and radial momentum equations respectively become

[TABLE]

[TABLE]

[TABLE]

Two other equations that will be useful in the following sections are the total energy and total radial momentum equations, which are given by the $\nu=t$ and $\nu=r$ components of Equation (8). These respectively read

[TABLE]

[TABLE]

3.3 Strong shock jump conditions

Equations (15) – (17) govern the evolution of the post-shock fluid. Assuming that the shock generates neither mass, energy, nor momentum, the fluxes of these quantities must be conserved across the shock in the comoving frame of the shock itself. From Equations (15), (18), and (19), conserving these fluxes yields the following three jump conditions:

[TABLE]

[TABLE]

[TABLE]

Here $U^{\prime\prime}$ denotes the radial component of the four-velocity in the comoving frame of the shock, quantities with a subscript 2 are post-shock fluid quantities, and those with a subscript “a” pertain to the ambient medium. We will assume that the shock is sufficiently supersonic that the ambient pressure can be ignored, which is equivalent to stating that the energy behind the blastwave is much greater than the internal energy of the ambient medium. In this case, these three equations can be combined into the following cubic to be solved for $U^{\prime\prime}_{2}$ in terms of $U^{\prime\prime}_{\rm a}$:

[TABLE]

Of the three algebraic solutions to this equation, only one is purely real for all $U_{\rm a}^{\prime\prime}$ and hence physical. If we denote the radial component of the lab-frame four-velocity of the shock by $U_{\rm s}$, then from Lorentz transformations it follows that $U^{\prime\prime}_{\rm a}=-U_{\rm s}$, while the radial component of the lab-frame post-shock velocity $U_{\rm 2}$ is given by

[TABLE]

With the solution for the comoving, post-shock fluid velocity in terms of the lab-frame velocity of the shock from Equation (23), Equation (24) yields the lab-frame post-shock fluid velocity for arbitrary shock speeds. In the nonrelativistic limit where $U_{\rm s}\ll 1$, the full expression reduces to $U_{2}/U_{\rm s}=2/(\gamma+1)$, while in the ultra-relativistic limit where $1+U_{\rm s}^{2}\simeq U_{\rm s}^{2}$ it becomes $U_{2}/U_{\rm s}=\sqrt{2/\gamma-1}$. Figure 1 shows the exact solution for the ratio $U_{2}/U_{\rm s}$ (purple, solid) when $\gamma=4/3$, the non-relativistic limit of 6/7 (black, dashed), and the ultra-relativistic limit of $1/\sqrt{2}$ (black, dotted).

The post-shock comoving density is simply found by inverting Equation (20)

[TABLE]

while the pressure is obtained by rearranging Equation (22) and using Equation (25) to remove the dependence on $\rho^{\prime}_{2}$:

[TABLE]

As for the post-shock velocity, these jump conditions reduce to the known non-relativistic and ultra-relativistic limits when $U_{\rm s}\ll 1$ and $U_{\rm s}\gg 1$. The solid, green curve in the left panel of Figure 1 shows the full solution for the post-shock pressure normalized by $U_{\rm s}^{2}$, and the horizontal, dot-dashed line gives the ultra-relativistic limit of $2-\gamma=2/3$ when $\gamma=4/3$ (the dashed, horizontal line shows the non-relativistic limit of $p^{\prime}/U_{\rm s}^{2}=2/(\gamma+1)=6/7$). The right panel of Figure 1 illustrates the full solution for the comoving density normalized by the Lorentz factor of the shock (solid, blue curve) when $\gamma=4/3$, the dashed, horizontal line shows the non-relativistic limit of $(\gamma+1)/(\gamma-1)=7$, and the dotted line depicts the ultra-relativistic limit of $\sqrt{2\gamma-\gamma^{2}}/(\gamma-1)=2\sqrt{2}$.

In this paper we are interested in the consequences of mildly relativistic velocities on the propagation of a strong shock and the post-shock fluid, i.e., where $U_{\rm s}/c\lesssim few$ (recall that $U_{\rm s}$ is the radial component of the four-velocity, not the three-velocity). In this case, we can Taylor expand our expressions for the post-shock fluid quantities to the lowest, non-zero order in $U_{\rm s}$ beyond the non-relativistic limit. Doing so, we find that Equations (24) – (26) become

[TABLE]

[TABLE]

[TABLE]

where for clarity we included in these relations the fact that they hold at the location of the shock position, $R$. In agreement with our arguments in Section 2, these expressions contain modifications to the non-relativistic jump conditions at the order $U_{\rm s}^{2}/c^{2}$.

3.4 Relativistic, Conserved energy

To initiate the outward motion of the fluid and the formation of the strong shock, there is an assumed-impulsive injection of a large amount of energy $E_{\rm ej}$ (large relative to the ambient internal energy and any gravitational potential energy) into the ambient medium, which for typical supernovae is on the order of $E\simeq 10^{50-52}$ erg (though for failed and very weak supernovae, the shock energy is $E_{\rm ej}\simeq 10^{47-48}$ erg or less). Owing to the fact that the ambient medium is assumed to be pressureless, there is no source of thermal energy as the shock moves out. Therefore, under the assumption that there are no sources or sinks of energy interior to the flow (or that, if there is a compact object such as a black hole or neutron star accreting matter, the binding energy drained by the compact object is small compared to the initial energy), this initial energy must be conserved as the shock propagates outward.

It is tempting to relate this conserved energy to the integral of the energy density over the volume enclosed by the shock, i.e.,

[TABLE]

(The lower bound on this integral, $R_{\rm c}$, can be non-zero for certain combinations of the ambient density profile and the adiabatic index of the post-shock flow, where the solution terminates in a contact discontinuity; see Section 4 below.) However, this energy is not conserved, because the shock sweeps up rest mass energy from the ambient medium as it moves outward. Indeed, this additional source of energy can be seen directly by integrating Equation (18) over the volume enclosed by the shock and rearranging:

[TABLE]

The right-hand side is the change in the total energy budget due to the addition of rest mass energy from the ambient medium. Thus, the total energy $E_{\rm tot}$ is not equivalent to the initial, injected energy associated with the explosion. However, we see from integrating the continuity equation from $R_{\rm c}$ to $R(t)$ that

[TABLE]

where

[TABLE]

Using this expression for the right-hand side of Equation (31) and rearranging then gives

[TABLE]

where

[TABLE]

This energy is manifestly conserved as the shock advances into the ambient medium, and is the exact, relativistic analog of the energy behind the blast that is relatable to the initial energy of the explosion.

Including lowest-order, relativistic corrections to this integral for the conserved energy gives

[TABLE]

The first term is the familiar, Newtonian expression for the total (kinetic plus internal) energy behind the blastwave, and the Sedov-Taylor solution exactly conserves this quantity. The second term is the lowest-order relativistic correction, and scales as $\propto U^{4}$ (recall that the pressure behind the shock, by virtue of the shock jump conditions, is comparable to the ram pressure of the shock; specifically see Equation 29). These relativistic corrections to the energy imply that, if the energy is to be absolutely conserved, there must also be corrections to the shock velocity that account for these additional terms.

3.5 Self-similar equations and perturbed solutions to lowest relativistic order

When $U\ll 1$, Equations (15) – (17) reduce to the well-known Euler equations in spherical symmetry. Here, however, we are interested in the corrections to these equations that are induced by mildly relativistic flow, and hence we need to account explicitly for these corrections. The Sedov-Taylor blastwave also adopts the change of variables

[TABLE]

which removes the time dependence of the boundary conditions at the shock front and yields self-similar solutions of the form

[TABLE]

It is also useful to introduce the time-like variable

[TABLE]

so the temporal and spatial derivatives then transform as

[TABLE]

where we used the fact that the shock three-velocity is related to the four-velocity via $dR/dt=U_{\rm s}/\sqrt{1+U_{\rm s}^{2}}$. The continuity, entropy, and radial momentum equations then become

[TABLE]

[TABLE]

[TABLE]

We will now expand these equations to lowest, post-Newtonian order. Before doing so, however, we note that there are exact, self-similar solutions (being the Sedov-Taylor solutions) in the Newtonian case, and hence we expect that the velocity, density, and pressure of the fluid should be approximately given by these solutions but with perturbations that are introduced from the velocity scale set by the speed of light. Moreover, investigating the shock jump conditions (27) – (29), we see that relativistic terms modify the post-shock fluid quantities at the order $U_{\rm s}^{2}$. We therefore seek solutions for the velocity, comoving density, and pressure that are of the form

[TABLE]

[TABLE]

[TABLE]

From the shock jump conditions, the functions satisfy the following boundary conditions at the shock:

[TABLE]

[TABLE]

In addition to satisfying these boundary conditions, energy must also be globally conserved. Returning to Equation (36) and using these forms for the velocity, density, and pressure, the conserved energy is given by

[TABLE]

If the second term in this expression were absent, then the energy would be conserved if

[TABLE]

where $E_{*}$ is a constant that scales with the total energy, and this is just the familiar velocity-radius relation for the Sedov-Taylor blastwave. However, because the second, relativistic term modifies the energy, this scaling cannot hold exactly, as otherwise the relativistic corrections would violate energy conservation. There must therefore be relativistic corrections to the shock velocity, and energy conservation dictates that these corrections must be of the form

[TABLE]

where $\sigma$ is a dimensionless number. Inserting this expression into the above equation for the total energy and requiring that the relativistic terms cancel exactly yields

[TABLE]

For a given $\sigma$, the functions $f_{1}$, $g_{1}$, and $h_{1}$ are completely specified from the fluid equations. This fourth boundary condition, which enforces global energy conservation, will therefore only be satisfied for a certain value of $\sigma$ (for given $n$ and $\gamma$), making $\sigma$ an “eigenvalue” from the standpoint that it is uniquely determined by this additional, integral constraint.

The governing equations for the self-similar functions can now be obtained by inserting Equations (44) – (46) and (51) into Equations (41) – (43) and keeping the lowest-order, surviving terms in $U_{\rm s}$. Doing so and performing some algebra yields the following three equations for the non-relativistic quantities:

[TABLE]

[TABLE]

[TABLE]

These three equations can be integrated three times exactly, with the pressure being analytically related to the velocity and density as

[TABLE]

The equations for the perturbed functions are given by

[TABLE]

[TABLE]

[TABLE]

Equation (57) can be integrated, and using the boundary conditions at the shock front gives the following expression for the relativistic correction to the density:

[TABLE]

The integral constraint that determines the eigenvalue, Equation (52), can also be written as a fourth boundary condition on the functions $f_{1}$, $g_{1}$, and $h_{1}$: subtracting the continuity from the energy equation and integrating from $R_{c}(t)$ to $R(t)$ gives

[TABLE]

where we used Equation (60) to remove the dependence on $g_{1}$ and all of the functions are evaluated at $\xi_{\rm c}$. By virtue of Equation (56), this expression only contains a relativistic correction, and for energy-conserving solutions we therefore require that the term in braces (multiplied by $\xi_{\rm c}^{2}$) be equal to zero. Either this fourth boundary condition or the integral constraint (52) can be used to determine $\sigma$.

3.6 Shock position and unperturbed coordinates

Equation (51) can be rearranged and integrated numerically to yield the shock position as a function of time. However, we can also use the assumed-smallness of terms of order $U_{\rm s}^{2}/c^{2}$ to decompose the shock position and velocity into their non-relativistic and relativistically-corrected parts; denoting the non-relativistic shock position and three-velocity as $R_{0}$ and $V_{0}=dR_{0}/dt$ and their relativistically-perturbed counterparts as $R_{1}$ and $V_{1}$, we find

[TABLE]

where we introduced factors of $c^{2}$ for clarity. This expression demonstrates, as expected, that the unperturbed shock velocity and position are related via the standard, energy-conserving prescription for the Sedov-Taylor blastwave:

[TABLE]

where we defined $V_{\rm i}^{2}=E_{*}/(4\pi\rho_{\rm a}r_{\rm a}^{3})$ as the unperturbed velocity of the shock when the shock position coincides with $R=r_{\rm a}$, while the relativistic corrections to the shock position and velocity satisfy

[TABLE]

This equation can be integrated to yield, if $n\neq 1$,

[TABLE]

[TABLE]

while if $n=1$

[TABLE]

[TABLE]

The second term in parentheses in Equations (65) and (66) is a consequence of initial conditions, and arises from the fact that scale invariance allows us to define the relativistic corrections to the shock position to be zero at $t=0$. Interestingly, if $\sigma\equiv 1$, then the relativistic corrections to the shock velocity and position are exactly zero. This effect arises from a competition between the increase in the four-velocity generated by positive $\sigma$, and time dilation that reduces the three-velocity from the four-velocity – when $\sigma=1$ these two effects exactly balance to yield no relativistic correction to the shock velocity.

Following Coughlin et al. (2019), we wrote our solutions for the relativistic corrections to the velocity, density, and pressure of the post-shock fluid in terms of the true shock position and velocity and the total self-similar variable $\xi=r/R(t)$. While formally correct to order $U_{\rm s}^{2}/c^{2}$, these expressions (specifically Equations 44 – 46) also contain terms that are of a higher order than $U_{\rm s}^{2}/c^{2}$. We can remove these additional terms by rewriting the solutions in terms of the “unperturbed” self-similar variable $\xi_{0}=r/R_{0}(t)$ and using Equations (65) and (66) to write the corrections to the shock position and velocity in terms of their non-relativistic counterparts; the resulting expressions, which are identical to Equations (44) – (46) to order $U_{\rm s}^{2}/c^{2}$, are

[TABLE]

[TABLE]

[TABLE]

Finally, the three-velocity of the fluid is

[TABLE]

and the lab-frame density is given by

[TABLE]

4 Solutions

Here we present the numerical solutions for the functions $f_{1}$, $g_{1}$, $h_{1}$, and the eigenvalue $\sigma$ that satisfy the differential Equations (58) and (59), with Equation (60) relating $g_{1}$ to $f_{1}$, and the fourth, energy-conserving boundary condition, given either by the integral constraint (52) or Equation (61). Since the functions $f_{1}$, $g_{1}$, and $h_{1}$ satisfy the boundary conditions at the shock front, given by Equation (48), we can numerically integrate Equations (58) and (59) from $\xi=1$ for a given, initial guess for $\sigma$. We then perturb the guess for $\sigma$ and calculate the change in the energy residual, i.e., we determine how well the new value of $\sigma$ satisfies the fourth boundary condition given by Equation (52) or (61), which motivates the choice for the next $\sigma$ that will better satisfy the fourth boundary condition. In this way, we iteratively determine the eigenvalue $\sigma$ that globally conserves the energy behind the blastwave.

The left panel of Figure 2 shows the self-similar velocity, density, and pressure for the Sedov-Taylor blastwave when $\gamma=4/3$ and $n=0$, corresponding to a radiation-pressure dominated, post-shock fluid and a constant ambient density. The linear velocity, constant pressure, and approximately zero density near the origin are familiar features of the Sedov-Taylor blastwave. The right panel of this figure gives the relativistic corrections to the velocity, density, and pressure for this combination of $n$ and $\gamma$. The eigenvalue that results in the exact conservation of energy, including the relativistic terms, is $\sigma\simeq 1.078$. The right panel of this figure also shows that, very near the shock front, the perturbation to the density is positive, and this is just a consequence of the jump conditions at the shock. However, the perturbation to the density drops very steeply from the shock front inward, and becomes negative and reaches a relative minimum near $\xi\simeq 0.95$. This qualitative behavior implies that the material behind the shock becomes increasingly confined to a region very near the shock front, and that relativistic effects cause the material to be swept into an even thinner shell than the one predicted by the Sedov-Taylor blastwave alone. We also see that the perturbation to the pressure, while it does asymptote to a constant near the origin, shows much more variability than the unperturbed solution for $\xi\gtrsim 0.1$. This feature demonstrates that the total pressure behind the blast wave shows more spatial variation when relativistic effects are included, which is a familiar property of the ultra-relativistic, Blandford-McKee blastwave.

The left panel of Figure 3 gives the correction to the shock position (solid, purple) and the correction to the shock speed (dashed, red), each normalized by its non-relativistic counterpart, as functions of time in units of $r_{\rm a}/c$. Here we set $n=0$, corresponding to a constant ambient density, the post-shock adiabatic index to $\gamma=4/3$, and the initial shock velocity to $V_{\rm i}/c=0.5$. We see that the relativistic correction to the shock velocity is initially positive, corresponding to an increase in the shock position over the non-relativistic value, while at late times the correction to the shock speed changes sign; this behavior is due to the competition between the effects of positive-$\sigma$, which increases the four-velocity (see Equation 51), and time dilation, which reduces the three-velocity over the four-velocity. In the asymptotic limit of $t\rightarrow\infty$, both of these corrections decay to zero and the flow settles into the non-relativistic regime.

The right panel of Figure 3 shows the post-shock fluid three-velocity, normalized by the initial shock speed, as a function of normalized radial position behind the shock front (see Equation 72). As for the left panel of this figure, here the ambient density is constant ($n=0$), the adiabatic index is $\gamma=4/3$, and the initial shock speed is $V_{\rm i}/c=0.5$. The different, colored curves correspond to the times shown in the legend, and the black, dashed curve is the Sedov-Taylor solution for this combination of $n$ and $\gamma$ (and is identical to the blue curve in the left panel of Figure 2). We see that near the shock front the post-shock velocity is slightly increased, which is reasonable given the slight increase in the shock velocity itself, as demonstrated in the left panel of this figure. However, at small radii the velocity falls significantly below the non-relativistic prediction, and also displays more nonlinear behavior near the origin.

The left panel of Figure 4 illustrates the normalized, lab-frame, post-shock density as a function of spherical radius $r$ normalized by the shock radius, while the right panel of this figure gives the post-shock pressure normalized by the shock velocity (see Equations 73 and 71 respectively). As for the right panel of Figure 3, the different colored curves correspond to the times in the legend, the black, dashed curves are the Sedov-Taylor prediction, and we set $n=0$ (constant ambient density), $\gamma=4/3$, and $V_{\rm i}/c=0.5$ (note that the $x$-axis in the left panel is compressed to highlight the behavior near the shock; for $\xi_{0}\lesssim 0.7$, all of the functions rapidly approach zero and show little deviation from one another). We see that both the post-shock density and pressure increase above the Sedov-Taylor solution near the shock front, but, as is also true for the post-shock velocity, each of these quantities declines and falls below the non-relativistic prediction at a radius not far behind the shock front. The post-shock pressure also shows significant deviation from the nearly-constant value expected from the Sedov-Taylor blastwave. These findings confirm that relativistic effects tend to compress the fluid into a more confined region immediately behind the shock.

The left Panel of Figure 5 gives the post-shock, relativistic correction to the velocity profile of the fluid, the right panel of this figure shows the relativistic correction to the post-shock, comoving density, and the left panel of Figure 6 illustrates the post-shock correction to the pressure, and all of these panels set $n=0$ (constant density ambient medium). The different curves in each of these figures correspond to the adiabatic indices shown in the legend. We see that, while all of these curves show the same qualitative trends, relativistic effects become amplified as the adiabatic index of the gas decreases: the magnitude of the velocity reduction is more pronounced; the material behind the shock becomes increasingly compressed to the shock itself; and the pressure has increased variation near the shock, possesses more nonlinear behavior, and the decrease near the origin is enhanced. These findings – that relativistic effects become more important for smaller adiabatic indices – are consistent with the fact that the eigenvalue $\sigma$ increases as $\gamma$ decreases (see Table 1).

The right panel of Figure 6 shows the relativistic correction to the shock velocity, normalized by the non-relativistic shock speed, as a function of time in units of $r_{\rm a}/c$. The different curves are appropriate to the adiabatic indices shown in the legend, the ambient density profile is constant ($n=0$), and we set the initial, non-relativistic shock speed to $V_{\rm i}/c=0.5$. We see that, when $\gamma$ is small, the initial velocity increases compared to the non-relativistic one, but becomes negative after a time of $ct/r_{\rm a}\simeq 2.57$ (the exact time at which this occurs can be derived analytically from Equation 66). However, for values of $\gamma\gtrsim 1.4$, this behavior inverts, with the initial correction to the velocity being negative at early times and transitioning to positive after $ct/r_{\rm a}\simeq 2.57$. This inversion occurs mathematically because $\sigma$ crosses $\sigma=1$ at this location, which is where the relativistic correction to the three-velocity equals zero (this feature can be seen by inserting $\sigma=1$ into Equation 51), and – as we noted above – arises physically because the three-velocity of the fluid is affected by time dilation.

Figure 7 illustrates the self-similar correction to the post-shock four-velocity for $\gamma=4/3$ and the range of $n$ shown in the legend (the density profile of the ambient medium falls off with spherical radius $r$ as $\rho^{\prime}\propto r^{-n}$). We see that, as $n$ increases, the relativistic correction to the velocity decreases in magnitude, and the solution with $n=2.5$ ends at a contact discontinuity – coinciding with a point in the flow where the non-relativistic velocity profile satisfies $f(\xi_{\rm c})=\xi_{\rm c}$, such that the fluid elements at these locations are stationary with respect to the shock. The right panel of this figure shows the correction to the pressure profile behind the shock, again with $\gamma=4/3$. As $n$ increases from [math] to $2$, the pressure profile exhibits progressively less deviation behind the shock, and the overall magnitude of the correction is reduced. When $n=2.5$, the post-shock correction to the pressure equals zero at the contact discontinuity, which is also where the non-relativistic pressure equals zero.

The left panel of Figure 8 gives the correction to the self-similar, comoving density behind the shock for $\gamma=4/3$ and the power-law indices of the ambient medium, $n$, in the legend. As for the pressure, the relativistic contributions tend to be less pronounced as $n$ increases from [math] to 2, with less of the variation being confined to the immediate vicinity of the shock and the magnitude of the variation reduced. However, when $n=2.5$, we see that the magnitude of the correction increases again, and the function $g_{1}$ actually diverges weakly at the location of the contact discontinuity. The right panel of this figure shows, for reference, the Sedov-Taylor solution for a subset of $n$, with the vertical, dashed lines coinciding with the location of the contact discontinuity. We see that, when $n=2.5$, the non-relativistic density at the contact discontinuity equals zero, while for $n=2.7$ the function $g_{0}$ diverges at the contact discontinuity.

The relativistic correction to the lab-frame, three-velocity of the shock $V_{1}$, plotted relative to the non-relativistic velocity $V_{0}$ as a function of time, is shown in the left panel of Figure 9; here we set $V_{\rm i}/c=0.5$, where $V_{\rm i}$ is the Newtonian shock speed at $t=0$. Different curves correspond to different radial power-law indices of the density of the ambient medium, $n$, such that $\rho^{\prime}_{\rm a}\propto r^{-n}$. The correction is negative for larger values of $n$, which signifies that time dilation actually reduces the shock speed below the Sedov-Taylor prediction in spite of the fact that the velocity is increased in the comoving frame of the non-relativistic shock. Relativistic effects remain important for longer times as $n$ increases, and this occurs because the Sedov-Taylor velocity – to the square of which, as shown in Equation (66), the relativistic correction is proportional – declines less slowly as the density profile of the ambient medium steepens (which arises physically from the reduced momentum flux across the shock for larger $n$).

The eigenvalue $\sigma$, which conserves the relativistic energy behind the shock and relates the shock speed to the explosion energy, is shown as a function of $\gamma$ in the right panel of Figure 9. Each curve corresponds to the power-law index of the ambient medium, $n$, shown in the legend, and the horizontal, dashed line simply indicates where $\sigma=1$ for clarity. The fact that $\sigma$ is always positive implies that observers moving with the non-relativistic, Sedov-Taylor shock speed see a relativistically-boosted shock velocity in that frame. However, it is only for $\sigma>1$ that observers in the lab frame (i.e., the frame in which the ambient medium is at rest) also measure a positive, relativistic increase to the shock three-velocity; for $\sigma<1$ time dilation reduces the velocity below the Sedov-Taylor solution, as can be seen directly in the left panel of this figure. The values of $\sigma(n,\gamma)$ which were used to make this figure are given in Table 1.

It can be shown from the integrals of Equations (53) – (55) that the self-similar, Sedov-Taylor density scales as

[TABLE]

which implies that when the Sedov-Taylor solution ends in a contact discontinuity, the density diverges at that location if $n>6/(\gamma+1)\simeq 2.57$ for $\gamma=4/3$ (see also Goodman 1990). We do not find any solution for $\sigma$ that satisfies global energy conservation for power-law indices that are steeper than this value. Mathematically, solutions in this regime do not exist because we can combine the integral constraint on the energy (52) and Equation (60) to show that, if the relativistic correction to the energy is to remain finite, then we must have111Note that this is not required when the density equals zero at the contact discontinuity, as in this situation the density scales as $g_{1}\propto(1-\xi_{\rm c})^{-\alpha}$ near $\xi_{\rm c}$ with $\alpha<1$; thus, while the correction to the comoving density diverges at $\xi_{\rm c}$, it does so in a way that yields a finite, relativistic correction to the energy. $f_{1}(\xi_{\rm c})=\xi_{\rm c}(\xi_{\rm c}^{2}-1)/2$. However, from Equation (61) and the fact that $h_{0}(\xi_{\rm c})=0$, we see that we must also have $h_{1}(\xi_{\rm c})=0$ to maintain energy conservation. The system is therefore over-constrained when the density diverges at the contact discontinuity, and these two boundary conditions at $\xi_{\rm c}$ will not, in general, be satisfied simultaneously for a single $\sigma$. In particular, we find that the solution for which $f_{1}(\xi_{\rm c})=\xi_{\rm c}(\xi_{\rm c}^{2}-1)$ possesses a finite, but non-zero pressure at the contact discontinuity, meaning that these solutions cannot simultaneously satisfy both the integral constraint (52) and Equation (61).

It is also not surprising from a physical standpoint that these diverging-density solutions are problematic. For one, it was shown by Goodman (1990) that such solutions are unstable to aspherical perturbations, as the decelerating nature of the fluid and the density inversion renders the contact discontinuity susceptible to the Rayleigh-Taylor instability. These Sedov-Taylor blastwaves therefore cannot be manifested in any physical (i.e., one with permissible angular deviations from spherical symmetry) scenario. The diverging-density solutions also violate the self-similar hypothesis that the flow is predominantly characterized by the physical conditions at a single point within the flow: for small $n$ and $\gamma$, the vast majority of the mass is contained very near the shock front. However, when the density profile becomes inverted, most of the mass is concentrated near the contact discontinuity, and the causal connectedness of the shocked fluid implies that the shock “knows” this property of the inner flow. It is therefore likely that, for the diverging-density solutions, the physical conditions both at the shock front and near the contact discontinuity remain important for establishing the long-term behavior of the shock, and the solution may be fundamentally non-self-similar.

It is interesting to note that the ultra-relativistic, Blandford-McKee blastwave avoids this issue, as the comoving pressure behind the shock transforms as a higher power of the shock Lorentz factor than the comoving density. The contribution of the kinetic energy to the total energy behind the blastwave is therefore dropped in their solution, and only the internal energy must be conserved. Thus, while the Blandford-McKee solution formally exists in this regime, it is likely that the integrated kinetic energy remains important for the dynamics (i.e., by virtue of the boundary conditions at the shock, it is true in the ultra-relativistic limit that the kinetic energy is sub-dominant to the internal energy, but that may not be true deeper within the flow; the non-existence of relativistic corrections to the Sedov-Taylor solutions suggests that this may be the case for, at least, smaller values of the Lorentz factor).

Finally, one could argue that our prescription for the behavior of the relativistic corrections was too restrictive, and that while solutions of the form given by Equations (44) – (46) do not exist in this regime, there may be other, more general solutions that do satisfy energy conservation and the modified boundary conditions at the shock. Instead of writing the expressions for the perturbations as in (44) – (46), we could have parameterized them as

[TABLE]

and similarly for the other variables, and therefore maintained the derivatives with respect to $\chi$ in Equations (57) – (59) (recall that $\chi=\ln R$ is a time-like variable). However, it is difficult to see how relativistic effects (to order $V^{2}/c^{2}$) could modify the solution in a way other than one that scales as $V^{2}/c^{2}$, as this is the only physical smallness parameter introduced in the problem. Therefore, any additional time dependence in the expression $f_{1}(\xi,\chi)$ should, instead, be regarded as a higher-order correction to the self-similar solution; we also note that this is precisely the motivation for expanding the self-similar functions as $f(\xi,\chi)\simeq f_{0}(\xi)+U_{\rm s}^{2}f_{1}(\xi)$ – the second term accounts for the time dependence that is induced by relativistic effects, and that time dependence can only physically be of the form $U_{\rm s}^{2}$. While we acknowledge that this is not a rigorous proof of the non-existence of more general solutions, we find it suggestive that no such solutions exist.

5 Summary and Conclusions

When the shockwave from an astrophysical explosion is strong (Mach number much greater than one) and non-relativistic, and therefore characterized by a shock speed $V$ much less than the speed of light $c$, the Sedov-Taylor, energy-conserving blastwave provides an analytic, self-similar solution for the temporal evolution of the shock itself, and the time and space-dependent evolution of the post-shock velocity, density, and pressure. In the other, extreme-relativistic limit where $V\simeq c$, the Blandford-McKee blastwave gives the energy-conserving evolution of the shock Lorentz factor and the post-shock fluid quantities. In between these two extremes, the finite speed of light introduces an additional velocity scale into the problem, which destroys the pure self-similarity of the solutions.

In this paper, we analyzed the leading-order, relativistic corrections to the fluid equations – which enter as $\mathcal{O}(V^{2}/c^{2})$ – to understand the effects that such relativistic terms have on the non-relativistic, Sedov-Taylor solution for strong shock propagation. By treating such terms as perturbations (i.e., ignoring nonlinear terms that enter as higher powers of $V^{2}/c^{2}$), we showed that there are relativistic corrections to the Sedov-Taylor solutions for the post-shock fluid quantities that vary (consistent with expectations) as $V^{2}/c^{2}$. In particular, we demonstrated that the radial component of the post-shock, fluid four-velocity can be written as $U=U_{\rm s}\left\{f_{0}(\xi)+U_{\rm s}^{2}/c^{2}\times f_{1}(\xi)\right\}$, where $U_{\rm s}$ is the shock four-velocity, $f_{0}$ is the Sedov-Taylor solution, $\xi=r/R(t)$ with $R(t)$ the shock position, and $f_{1}$ is a function that is self-consistently determined from the fluid equations and the relativistic jump conditions at the shock. The function $f_{1}$, and the analogous functions $g_{1}$ and $h_{1}$ for the density and pressure, respectively, induce more nonlinear behavior to the post-shock velocity, further compress the post-shock material to the immediate vicinity of the shock itself, and generate greater variation in the post-shock pressure as compared to the Sedov-Taylor, non-relativistic limit (see Figures 2 – 4). These are all features of the ultra-relativistic, Blandford-McKee blastwave, where the pressure declines rapidly behind the shock and all of the material is swept into a shell of width $\Delta R=R/U_{\rm s}^{2}$.

In addition to the post-shock fluid quantities, we also determined the relativistic correction to the velocity of the shock itself. We denoted this additional correction by an “eigenvalue” $\sigma$, such that the relativistically-corrected shock four-velocity is written implicitly as $U_{\rm s}^{2}R^{3-n}=E\left(1+\sigma U_{\rm s}^{2}/c^{2}\right)$. When $\sigma\equiv 0$, one recovers the familiar relationship between the shock velocity and position that guarantees the conservation of the blast energy, $E$. However, owing to the existence of relativistic corrections to the energy, $\sigma$ cannot be exactly zero, and there must therefore be corrections to the shock velocity that maintain total (i.e., including relativistic terms to order $U_{\rm s}^{2}/c^{2}$) energy conservation. For all of our solutions, the value of $\sigma$ was found to be positive, implying an increase to the shock four-velocity from relativistic effects; equivalently, observers moving at the non-relativistic shock speed measure a small, slight increase to the shock velocity, and hence the true shock position leads the non-relativistic one in the comoving frame of the non-relativistic shock. Nevertheless, the lab-frame three-velocity is reduced from the four-velocity by time dilation, which counterbalances the relativistically-boosted effect of positive $\sigma$, and if $\sigma\equiv 1$ these two effects exactly offset to yield a shock three-velocity that is identical to the Newtonian, Sedov-Taylor value. Interestingly, for a radiation-pressure dominated, $\gamma=4/3$ post-shock equation of state and a constant-density ambient medium, the shock three-velocity is slightly increased over the non-relativistic value, while for all declining density profiles with $n>0.5$ – where the ambient density $\rho^{\prime}_{\rm a}$ falls off with spherical radius $r$ as $\rho^{\prime}_{\rm a}\propto r^{-n}$ – the lab-frame velocity is reduced below the non-relativistic value (see Figure 9 and Table 1 for the values of $\sigma$ for a range of $n$ and $\gamma$).

Once the power-law index of the ambient density profile equals or exceeds the critical value $n_{\rm cr}(\gamma)=6/(\gamma+1)$, the density at the contact discontinuity present in the Sedov-Taylor solution diverges. For Sedov-Taylor blastwaves with $n>n_{\rm cr}$, we do not find any solution for $\sigma$ that maintains a finite and conserved relativistic correction to the energy, which mathematically follows directly from the nature of the solutions for $f_{1}$, $g_{1}$, and $h_{1}$ and the integral that maintains the conservation of energy (see Equations 52 and 61). It is possible that more general, time-dependent solutions (i.e., those that do not assume the form given by Equations 44 – 46) could be found in this regime that do conserve energy. However, based on the argument that additional time dependence from the self-similar, Sedov-Taylor solutions is itself seeded by relativistic effects, it is difficult to see how any extra time dependence would not be in the form of higher-order (than $U_{\rm s}^{2}/c^{2}$) terms. We therefore find it unlikely that such generalized solutions exist at the leading relativistic order.

Here we described the leading-order, relativistic corrections to the fluid flow, which enter into the fluid equations and the boundary conditions as $U_{\rm s}^{2}/c^{2}$. One can regard our solutions as the first in a series expansion of the fluid equations in $U_{\rm s}^{2}$, and the next-order solution for (for example) the four-velocity would be $U=U_{\rm s}\left\{f_{0}(\xi)+U_{\rm s}^{2}f_{1}(\xi)+U_{\rm s}^{4}f_{2}(\xi)\right\}$; one could then, by expanding the fluid equations to the next order, derive self-consistent equations for $f_{2}$, $g_{2}$, and $h_{2}$, and the boundary conditions at the shock could be found by expanding the general jump conditions to the next-highest order. In principle, one should also be able to construct the next-order solution for the Blandford-McKee solution, and approach the problem from the other, ultrarelativistic direction by finding the next-highest-order correction in $1/\Gamma$.

In this paper we focused on ambient density profiles less steep than $\rho^{\prime}_{\rm a}\propto r^{-3}$. For steeper density profiles, the shock enters an accelerating regime, and the self-similar solutions for the post-shock fluid quantities are provided by Waxman & Shvarts (1993) (see also Koo & McKee 1990). In this case, the self-similar flow is constrained to lie between a sonic point within the interior of the flow and the shock front, and energy and mass are drained into a non-self-similar, inner region that is causally disconnected from the fluid near the shock. One could apply all of the formalism developed in this paper to derive the leading-order, relativistic corrections to such accelerating, self-similar flows, and the resulting equations for the functions $f_{1}$, $g_{1}$, and $h_{1}$ would, in fact, appear almost identical to Equations (57) – (59), with the exception that various factors of $n-3$ – which result from energy conservation – would be replaced by a numerical factor that is determined from the self-similar, accelerating solutions. When the shock accelerates, such corrections actually come to dominate the self-similar solution at late times, owing to the increasing nature of $U_{\rm s}^{2}/c^{2}$, and one could interpret this result by saying that such self-similar solutions are “unstable” to relativistic corrections.

These solutions for the relativistically-corrected shock speed and post-shock fluid quantities could be used for generating more accurate models for the late-time lightcurves of long gamma-ray bursts as well as the lightcurves of energetic supernovae. In particular, the relativistic beaming induced by the marginally-relativistic velocity, and the prediction for the time over which the shock speed declines to sub-relativistic speeds, would yield correspondingly different break timescales for the lightcurve of the event and peaks in the synchrotron spectrum (e.g., Sari et al. 1998; De Colle et al. 2012). The event AT2018cow (Prentice et al., 2018; Rivera Sandoval et al., 2018; Ho et al., 2019; Kuin et al., 2019; Margutti et al., 2019; Perley et al., 2019), tentatively an extreme example of the class of fast-rising transients (Drout et al., 2014), also provided evidence of a moderately relativistic outflow with speed $\sim 0.1c$; the model presented here could be combined with current multiwavelength data to further constrain properties of the progenitor and surrounding medium.

This work was supported by NASA through the Einstein Fellowship Program, Grant PF6-170170. I thank Brian Metzger, Eliot Quataert, and Jonathan Zrake for useful discussions.