text:
# Stability of an Ultra-Relativistic Blast Wave in an External Medium with a Steep Power-Law Density Profile

## I Introduction

The self-similar solutions of relativistic blast waves are of much interest because of their recent applications to the study of Gamma-Ray Bursts (GRBs). A sudden release of a large amount of energy within a small volume results in a strong explosion that drives a relativistic shock into the surrounding medium. At late times the blast wave approaches a self-similar phase whereby its speed and the distribution of the pressure, density, and velocity of the gas behind the shock front do not depend on the length and time scales of the initial explosion, but only on the explosion energy and the properties of the unshocked external medium. The self-similar solutions describing this phase have been first studied by Blandford and McKee  (hereafter, BMK). We list their central results, which are relevant to this paper, in section §II.A. Note that in the BMK solution the total energy released in the explosion $`E`$ is the only relevant parameter.

The BMK solution is only valid for $`k<4`$, where $`k`$ is the power-law index of the radial density profile of the external medium, i.e., $`\rho _1r^k`$. When $`k4`$, the similarity variable defined by Blandford and McKee  is no longer appropriate. Even in the range $`3k<4`$, the validity of the BMK solution is not justified, because the mass contained behind the shock front diverges if the density profile of the shocked fluid is described by the BMK solution.

The self-similar solutions for steep density profiles with a power-law index $`k4`$ were derived recently by Best & Sari The more extreme case of an exponential density profile has been discussed by Perna & Vietri .. The derivation of these solutions is similar to that in the non-relativistic regime. The self-similar solutions to a non-relativistic blast wave were discovered independently by Sedov  , Von Neumann , and Taylor . The so-called ”Sedov-Von Neumann-Taylor blast wave” solutions exist only for $`k<5`$, but Waxman & Shvarts  showed that in the range $`3k<5`$ these solutions fail to describe the asymptotic flow because the energy diverges; instead they found second-type self-similar solutions for $`3<k<5`$ as well as for $`k5`$. The new class of non-relativistic, self-similar solutions describe the flow in a limited spatial region $`D(t)rR(t)`$, where $`R(t)`$ is the shock radius and $`D(t)`$ coincides with a $`C_+`$ characteristic so that the flow inside the region $`r<D(t)`$ does not affect the flow in the outer self-similar region. The self-similar solution has to cross the sonic line into the region where the $`C_+`$ characteristic can not catch-up with the shock front. The solution describes a shock accelerating with a temporal dependence whose power-law index is uniquely determined by requiring that the self-similar solution cross the sonic line at a singular point. Note that in these second-type self-similar solutions the total energy released in the explosion $`E`$ is not a relevant parameter. Although the energy in the self-similar part of the flow approaches a constant as time diverges, the fraction of the explosion energy $`E`$ carried by the self-similar component depends on the details of the initial conditions. Thus, contrary to the BMK case, dimensional arguments can not be used to determine the power-law index of the temporal dependence. Instead, the singular point determines the temporal power-law index. In the ultra-relativistic regime, the second-type self-similar solutions for $`k4`$ can similarly be obtained by requiring that they cross the sonic line at a singular point. Best & Sari  found that these self-similar solutions exist for $`k>5\sqrt{3/4}`$ and describe accelerating shock waves. However, the properties of the flow in the self-similar region, such as the energy and mass contained in the region, were not discussed.

In this paper, we rederive the self-similar solutions of ultra-relativistic blast waves for $`k>4`$ using a different self-similar variable and discuss the properties of the flow in the self-similar regime. Our main goal is to study the stability of these self-similar solutions. The stability of the Waxman-Shvarts self-similar solutions in the non-relativistic regime was studied by Sari, Waxman & Shvarts . They found that shocks accelerating at a rate larger than a critical value and corresponding to solutions that diverge in finite time, are unstable for small and intermediate wavenumbers. Shocks that accelerate at a rate smaller than the critical rate are stable for most wavenumbers. The acceleration rate can be quantified by the measure $`\delta =R\ddot{R}/\dot{R}^2`$, where the dots denote time derivatives and $`R(t)`$ is the radius of the shock front. This measure provides the fractional change of the velocity over a characteristic time scale for evolution ($`R/\dot{R}`$). Solutions that diverge in finite time have $`\delta >1`$ while others have $`\delta <1`$. Thus, when shocks accelerate sufficiently fast they become unstable.

In the following sections we study the stability of the self-similar solutions of ultra-relativistic blast waves for steep density profiles with a power-law index $`k>4`$. The self-similar solutions are described in §II. We list the BMK solutions for $`k<4`$ in §II.A and derive the self-similar solutions for $`k>4`$ in §II.B. In §III, we discuss the properties of the self-similar flow and calculate the energy and mass contained in the self-similar regime. In §IV and §V, we study the stability of the self-similar solutions. Finally, we summarize our main results in §VI.

## II Self-Similar Solutions

### A BMK solutions for $`k<4`$

For pedagogical reasons, we first briefly outline the derivation of the self-similar solutions of relativistic blast waves for $`k<4`$ by Blandford & McKee. For a complete derivation, the reader is referred to the original paper .

Assuming an ultra-relativistic equation of state, $`p=(1/3)e`$, where $`p`$ and $`e`$ are the pressure and energy density measured in the fluid frame, the equations describing a relativistic, spherically-symmetric, perfect fluid can be written as,

$$\frac{d}{dt}(p\gamma ^4)=\gamma ^2\frac{p}{t},$$
(1)
$$\frac{d}{dt}\mathrm{ln}(p^3\gamma ^4)=\frac{4}{r^2}\frac{}{r}(r^2\beta ),$$
(2)
$$\frac{n^{}}{t}+\frac{1}{r^2}\frac{}{r}(r^2n^{}\beta )=0,$$
(3)
where $`n^{}`$ is the density as measured in the laboratory frame, $`\gamma `$ and $`\beta `$ are the Lorentz factor and velocity of the fluid, and

$$\frac{d}{dt}\frac{}{t}+\beta \frac{}{r}$$
(4)
is the convective derivative. Throughout this paper we set the speed of light $`c`$ to unity. Assuming that the blast wave is ultra-relativistic so that the Lorentz factor of the shock front $`\mathrm{\Gamma }`$ and the shocked fluid $`\gamma `$ are much larger than unity, we only search for solutions accurate to the lowest order in $`\gamma ^2`$ and $`\mathrm{\Gamma }^2`$.

The effective thickness of the blast wave is approximately $`R/\mathrm{\Gamma }^2`$, where $`R`$ is the radius of the shock front. Thus an appropriate choice of similarity variable is

$$\xi =\left(1\frac{r}{R}\right)\mathrm{\Gamma }^20.$$
(5)
Next we assume that the external medium has a scale-free, power-law density profile $`\rho _1r^k`$. Ignoring radiative losses, the total energy contained in the shocked fluid remains constant and so the Lorentz factor of the shock front evolves adiabatically as a power law,

$$\mathrm{\Gamma }^2t^m,m>1.$$
(6)
Keeping only terms up to order $`O(\mathrm{\Gamma }^2t)`$, the shock radius is then given by

$$R=t\left[1\frac{1}{2(m+1)\mathrm{\Gamma }^2}\right].$$
(7)
A more convenient similarity variable can be defined as

$$\chi =1+2(m+1)\xi =[1+2(m+1)\mathrm{\Gamma }^2]\left(1\frac{r}{t}\right).$$
(8)
In terms of $`\chi `$, the pressure, velocity, and density in the shocked fluid can be written as

$$p=\frac{2}{3}w_1\mathrm{\Gamma }^2f(\chi ),$$
(9)
$$\gamma ^2=\frac{1}{2}\mathrm{\Gamma }^2g(\chi ),$$
(10)
$$n^{}=2n_1\mathrm{\Gamma }^2h(\chi ),$$
(11)
where $`\chi 1`$, $`w_1`$ and $`n_1`$ are the enthalpy and number density of the unshocked external medium. We assume that the unskocked external medium is cold, so that $`w_1`$ equals the energy density $`\rho _1`$. The jump conditions for a strong ultra-relativistic shock are satisfied by the boundary conditions

$$f(1)=g(1)=h(1)=1.$$
(12)
For an adiabatic impulsive blast wave, equations (1)–(3) admit a simple analytical solution, first derived by BMK

$$f=\chi ^{(174k)/(123k)},$$
(13)
$$g=\chi ^1,$$
(14)
$$h=\chi ^{(72k)/(4k)},$$
(15)
for

$$m=3k>1.$$
(16)
### B Self-similar solutions for $`k>4`$

In searching for self-similar solutions for $`k>4`$, we assume that the Lorentz factor of the shock front still obeys a power law, $`\mathrm{\Gamma }^2t^m`$ with $`m<1`$. When $`m<1`$, the similarity variable $`\chi `$ defined in equation (8) (and used by Best & Sari ) could be negative. For convenience we will use $`\xi `$, defined in equation (5), instead as our similarity variable. If at an initial time $`t_0`$, the shock radius is $`R_0`$ and the Lorentz factor of the shock front is $`\mathrm{\Gamma }_0`$, then at a later time $`t`$, to $`O(\mathrm{\Gamma }^2t)`$ the shock radius is given by

$$R=R_0+t\left[1\frac{1}{2(m+1)\mathrm{\Gamma }^2}\right]t_0\left[1\frac{1}{2(m+1)\mathrm{\Gamma }_0^2}\right].$$
(17)
We can rewrite this equation as

$$R=t\left[1\frac{1}{2(m+1)\mathrm{\Gamma }^2}\right]+a,$$
(18)
where $`a`$ is a constant dictated by the initial conditions. This equation for $`R`$ with $`m<1`$ differs from equation (7) by a constant $`a`$. However, we can choose the initial time $`t_0`$ such that $`a`$ is equal to zero. This is appropriate because of two reasons. First, the self-similar solutions are valid at much later times $`t>>t_0`$, thus the effect of the special choice of $`t_0`$ can be ignored. Second, what matters in the derivation of the self-similar solutions is the derivative of $`R`$, instead of $`R`$ itself. When $`a=0`$, the similarity variable becomes

$$\xi =\left(1\frac{r}{R}\right)\mathrm{\Gamma }^2=\mathrm{\Gamma }^2\frac{r}{t}\left[\mathrm{\Gamma }^2+\frac{1}{2(m+1)}\right].$$
(19)
Note that we have ignored higher order terms in $`\mathrm{\Gamma }^2`$ in the above expression.

Similarly to equations (9)–(11), we write the pressure, velocity, and density in the shocked fluid as

$$p=\frac{2}{3}w_1\mathrm{\Gamma }^2f(\xi ),$$
(20)
$$\gamma ^2=\frac{1}{2}\mathrm{\Gamma }^2g(\xi ),$$
(21)
$$n^{}=2n_1\mathrm{\Gamma }^2h(\xi ),$$
(22)
where $`\xi 0`$ and the boundary conditions,

$$f(0)=g(0)=h(0)=1,$$
(23)
correspond to the jump conditions for a strong ultra-relativistic shock.

We can now treat $`\mathrm{\Gamma }^2`$ and $`\xi `$ as two new independent variables in place of $`r`$ and $`t`$, and get

$$t\frac{}{t}=m\mathrm{\Gamma }^2\frac{}{\mathrm{\Gamma }^2}+\left[\mathrm{\Gamma }^2\frac{1}{2}\frac{m}{m+1}(m+1)\xi \right]\frac{}{\xi },$$
(24)
$$t\frac{}{r}=\left[\mathrm{\Gamma }^2+\frac{1}{2(m+1)}\right]\frac{}{\xi },$$
(25)
$$t\frac{d}{dt}=m\mathrm{\Gamma }^2\frac{}{\mathrm{\Gamma }^2}\left[\frac{1}{2}+(m+1)\xi \frac{1}{g}\right]\frac{}{\xi }.$$
(26)
In deriving the above equations, we have assumed that the blast wave is ultra-relativistic so that $`\mathrm{\Gamma }>>1`$ and $`\gamma >>1`$. Thus we only keep terms of the lowest contribution order in $`\mathrm{\Gamma }^2`$ and $`\gamma ^2`$.

Substituting equations (24)–(26) into equations (1)–(3), we obtain the following differential equations for $`f`$, $`g`$, and $`h`$:

$`2(3m+k)g+[g+2(m+1)g\xi +2]{\displaystyle \frac{1}{f}}{\displaystyle \frac{df}{d\xi }}`$ (27)
$`+2[g+2(m+1)g\xi 2]{\displaystyle \frac{1}{g}}{\displaystyle \frac{dg}{d\xi }}=0,`$ (28)
$`2(5m+3k8)g+3[g+2(m+1)g\xi 2]{\displaystyle \frac{1}{f}}{\displaystyle \frac{df}{d\xi }}`$ (29)
$`+2[g+2(m+1)g\xi +2]{\displaystyle \frac{1}{g}}{\displaystyle \frac{dg}{d\xi }}=0,`$ (30)
$`2(m+k2)g+[g+2(m+1)g\xi 2]{\displaystyle \frac{1}{h}}{\displaystyle \frac{dh}{d\xi }}`$ (31)
$`+2{\displaystyle \frac{1}{g}}{\displaystyle \frac{dg}{d\xi }}=0.`$ (32)
Using a new variable,

$$y=[1+2(m+1)\xi ]g,$$
(33)
we rewrite equations (28)–(32) as follows

$$\frac{1}{g}\left(\frac{1}{f}\frac{df}{d\xi }\right)=\frac{2[4(2m+k2)(m+k4)y]}{y^28y+4},$$
(34)
$$\frac{1}{g}\left(\frac{1}{g}\frac{dg}{d\xi }\right)=\frac{2[(7m+3k4)(m+2)y]}{y^28y+4}.$$
(35)
$`{\displaystyle \frac{1}{g}}`$ $`\left({\displaystyle \frac{1}{h}}{\displaystyle \frac{dh}{d\xi }}\right)=2[(y^28y+4)(y2)]^1`$ (36)
$`\times [(m+k2)y^2(10m+8k12)y`$ (37)
$`+(18m+10k16)].`$ (38)
One solution to the above equations is obtained for $`y=1`$ and $`m=3k`$, and corresponds to the BMK solution. From the definition of $`y`$, $`y=[1+2(m+1)\xi ]g`$, and the requirement that $`g`$ be positive it follows that this solution is only valid for $`m>1`$, i.e., $`k<4`$.

In our search for possible solutions with $`k>4`$, we start by analyzing equations (34)–(38). The right-hand-side of these equations diverges to infinity if $`y^28y+4=0`$. This corresponds to two singular points, $`y_1=42\sqrt{3}=0.536`$ and $`y_2=4+2\sqrt{3}=7.464`$. In addition, equation (38) has another singular point at $`y_3=2`$. The solution to equation (34) can bypass the singular points $`y_1`$ and $`y_2`$ if the numerator on the right-hand-side of the equation vanishes at $`y_1`$ or $`y_2`$. This gives

$$m_{1,2}=\frac{84k+(k4)y_{1,2}}{8y_{1,2}}.$$
(39)
It is easy to prove that when $`m=m_1`$, equations (35) and (38) will also bypass the singular point $`y_1`$ (the numerator in the right hand side of each equation is equal to zero at $`y_1`$). The same is true for $`m=m_2`$ and the singular point $`y_2`$. We will show below that when $`m=m_1`$ and $`k`$ is bigger than a critical value $`k_c`$, we have $`y1`$, thus equations (34)–(38) are able to bypass the singular point $`y_1`$ and never reach $`y_2`$ and $`y_3`$. The critical value $`k_c`$ can be calculated by setting $`m_1`$ equal to $`3k`$, the $`m`$ value corresponding to a BMK solution. For $`y_1=42\sqrt{3}`$, we get

$$m_1=12\sqrt{3}20+(32\sqrt{3})k.$$
(40)
Thus $`m_1=3k`$ gives us

$$k_c=5\frac{\sqrt{3}}{2}=4.134$$
(41)
The value of $`m_1`$ corresponding to $`k_c`$ is $`m_{1c}=2+\sqrt{3}/2=1.134`$. Thus when $`k>k_c`$, we have $`m_1<m_{1c}`$.

We now prove that when $`m=m_1`$ and $`k>k_c`$, we always have $`y1`$. Using equation (35) and the definition of $`y`$ in equation (33) we obtain

$`{\displaystyle \frac{dy}{d\xi }}`$ $`=`$ $`2(m+1)g(\xi )+{\displaystyle \frac{y}{g(\xi )}}{\displaystyle \frac{dg}{d\xi }}`$ (42)
$`=`$ $`{\displaystyle \frac{2g(\xi )[y^2+(m3k+12)y4(m+1)]}{y^28y+4}}.`$ (43)
When $`m=m_1`$, the above equation can be rewritten as

$$\frac{dy}{d\xi }=\frac{2g(yb)}{yy_2},$$
(44)
where

$$b=410\sqrt{3}+2\sqrt{3}k=1+2\sqrt{3}(kk_c).$$
(45)
When $`k>k_c`$, we have $`b>1`$. We also have $`g>0`$ and $`y_2=4+2\sqrt{3}>1`$, and so the right hand side of equation (44) is negative when $`y1`$. Since the boundary condition is $`y(\xi =0)=1`$, $`y(\xi )`$ must be a monotonically decreasing function of $`\xi `$ with $`y(\xi )1`$. The asymptotic behavior of $`y(\xi )`$ can be derived as follows. When $`\xi `$ is large, $`y`$ is negative and $`|y|`$ is large so that equation (44) can be approximated as

$$\frac{dy}{d\xi }2g\frac{1}{m_1+1}\frac{y}{\xi },$$
(46)
where we have used the approximation $`y2(m_1+1)\xi g`$ for large $`\xi `$. Equation (46) yields

$$y\xi ^{1/(m_1+1)}.$$
(47)
Note that when $`k>k_c`$, we have $`m_1<m_{1c}<1`$. Thus the exponent in the above power law is always positive.

It can be proven that when $`kk_c`$, equations (34)–(38) can not bypass all the singular points with either $`m_1`$ or $`m_2`$, and so the BMK solution is the only possible solution. When $`k>k_c`$, the equations can not bypass all the singular points with $`m_2`$. But this by itself is not sufficient for justifying that $`m_1`$ is the only viable choice. What if the solutions cut-off at some radius before reaching any singular points? We know that in order not to run into divergences of the energy or mass of the system, the solutions must be truncated at some radius $`r_+`$ (or $`\xi _+`$ in terms of the similarity variable), which should coincide with a $`C_+`$ characteristic line. The $`C_+`$ characteristic guarantees that the flow in the inner region $`r<r_+`$ will not influence the flow in the self-similar region $`r_+rR`$. This $`C_+`$ characteristic should not overtake the shock front in finite time, otherwise the self-similar region will eventually disappear. This argument has been applied in the non-relativistic case . We will prove below that in order to get to the regime where the $`C_+`$ characteristics can not catch the shock front, the solutions have to pass the singular point $`y_1`$, making $`m_1`$ the only viable choice.

First, let us derive the equation for a $`C_+`$ characteristic. We use $`v`$ to denote the fluid velocity in the laboratory frame. The sound speed in the fluid frame is $`u_s^{}=1/\sqrt{3}`$. Thus the sound speed in the laboratory frame $`u_s`$ is given by

$`u_s`$ $`=`$ $`{\displaystyle \frac{u_s^{}+v}{1+vu_s^{}}}=1{\displaystyle \frac{\sqrt{3}1}{\sqrt{3}+1}}{\displaystyle \frac{1}{2\gamma ^2}}`$ (48)
$`=`$ $`1{\displaystyle \frac{\sqrt{3}1}{\sqrt{3}+1}}{\displaystyle \frac{1}{\mathrm{\Gamma }^2g(\xi )}},`$ (49)
where we only keep the first-order term in $`\gamma ^2`$. Thus, a $`C_+`$ characteristic is described by

$$\frac{dr_+}{dt}=1\frac{\sqrt{3}1}{\sqrt{3}+1}\frac{1}{\mathrm{\Gamma }^2g(\xi _+)}.$$
(50)
We can rewrite this equation in terms of the similarity variable $`\xi _+`$. Using the definition of $`\xi `$ in equation (5), and the relations, $`dR/dt=11/(2\mathrm{\Gamma }^2)`$, $`d\mathrm{\Gamma }^2/dt=m\mathrm{\Gamma }^2/t`$, we obtain

$$\frac{dr}{dt}=1\frac{1}{2\mathrm{\Gamma }^2}\frac{(m+1)}{\mathrm{\Gamma }^2}\xi \frac{t}{\mathrm{\Gamma }^2}\frac{d\xi }{dt},$$
(51)
where we only keep the first-order term in $`\mathrm{\Gamma }^2`$. Substituting equation (51) into equation (50), we get

$$t\frac{d\xi _+}{dt}=\frac{\sqrt{3}1}{\sqrt{3}+1}\frac{1}{g(\xi _+)}(m+1)\xi _+\frac{1}{2},$$
(52)
which describes the evolution of a $`C_+`$ characteristic. We can further rewrite this equation as

$$t\frac{d\xi _+}{dt}=\frac{1}{2g(\xi _+)}(y_+y_1),$$
(53)
where

$$y_+=[1+2(m+1)\xi _+]g(\xi _+).$$
(54)
Equation (53) implies that when $`y_+>y_1`$, the right-hand-side of the equation is negative and so $`\xi _+`$ will decrease with time and the $`C_+`$ characteristic will approach the shock front. Only when $`y_+<y_1`$, $`\xi _+`$ will increase with time and the $`C_+`$ characteristic will not overtake the shock front. Also notice that the self-similar solution has the boundary condition $`y(\xi =0)=1>y_1`$. We thus proved that in order to get to the regime where $`C_+`$ characteristics can not overtake the shock front, the self-similar solution must pass through the singular point $`y_1`$, and therefore $`m_1`$ is the only viable choice.

We can now attempt to obtain the self-similar solutions for equations (34)–(38). For $`m=m_1`$, these equations become

$$\frac{1}{g}\left(\frac{1}{f}\frac{df}{d\xi }\right)=\frac{2(m_1+k4)}{yy_2},$$
(55)
$$\frac{1}{g}\left(\frac{1}{g}\frac{dg}{d\xi }\right)=\frac{2(m_1+2)}{yy_2},$$
(56)
$$\frac{1}{g}\left(\frac{1}{h}\frac{dh}{d\xi }\right)=\frac{2(m_1+k2)(yd)}{(yy_2)(y2)},$$
(57)
where

$$d=4+\sqrt{3}+\frac{3+2\sqrt{3}}{k+\sqrt{3}4}.$$
(58)
Treating $`y`$ as the independent variable instead of $`\xi `$ and making use of equation (44), equations (55)–(57) can be rewritten as

$$\frac{1}{f}\frac{df}{dy}=\frac{m_1+k4}{yb},$$
(59)
$$\frac{1}{g}\frac{dg}{dy}=\frac{m_1+2}{yb},$$
(60)
$$\frac{1}{h}\frac{dh}{dy}=\frac{(m_1+k2)(yd)}{(y2)(yb)}.$$
(61)
The boundary conditions are $`f(y=1)=g(y=1)=h(y=1)=1`$. Equations (59) and (60) have the solutions

$$f=\left(\frac{by}{b1}\right)^{m_1+k4},$$
(62)
$$g=\left(\frac{by}{b1}\right)^{m_1+2}.$$
(63)
A special case is obtained at $`k=6`$ ($`m_1=2`$) for which $`f(y)=g(y)=1`$. When $`b2`$ ($`k(15\sqrt{3})/3`$), equation (61) has the solution

$$h=\left[\left(\frac{1}{2y}\right)^{d2}\left(\frac{by}{b1}\right)^{db}\right]^{(m_1+k2)/(2b)}.$$
(64)
When $`b=2`$ ($`k=(15\sqrt{3})/3`$), equation (61) has the solution

$`h`$ $`=`$ $`(2y)^{m_1+k2}`$ (66)
$`\times \mathrm{exp}\left[(m_1+k2)(d2)\left({\displaystyle \frac{1}{2y}}1\right)\right].`$
In general, the functions $`f(\xi )`$, $`g(\xi )`$ and $`h(\xi )`$ do not admit simple analytical forms. Their values can be derived numerically from equations (62)–(66). For example, $`g(\xi )`$ satisfies the implicit algebraic equation

$$g(\xi )=\left(\frac{b[1+2(m_1+1)\xi ]g(\xi )}{b1}\right)^{m_1+2}.$$
(67)
But generally, we can derive the analytical forms for the asymptotic behaviors of $`f(\xi )`$, $`g(\xi )`$ and $`h(\xi )`$ in the limit of large $`\xi `$. In this limit, $`y`$ is negative and $`|y|`$ is large, and so equation (63) yields

$$g\left(\frac{y}{b1}\right)^{m_1+2}\left[\frac{2(m_1+1)\xi g}{b1}\right]^{m_1+2}.$$
(68)
We can solve $`g(\xi )`$ from the above equation and get

$$g(\xi )\alpha _g\xi ^{(m_1+2)/(m_1+1)},$$
(69)
where

$$\alpha _g=\left[\frac{2(m_1+1)}{b1}\right]^{(m_1+2)/(m_1+1)}.$$
(70)
Using equation (69) we obtain the asymptotic form of $`y(\xi )`$ for large $`\xi `$

$$y(\xi )\alpha _y\xi ^{1/(m_1+1)},$$
(71)
where

$$\alpha _y=\left[\frac{2(m_1+1)}{(b1)^{m_1+2}}\right]^{1/(m_1+1)}.$$
(72)
Using equation (71), we can derive the asymptotic form of $`f(\xi )`$ for large $`\xi `$ from equation (62) and get

$$f(\xi )\alpha _f\xi ^{(m_1+k4)/(m_1+1)},$$
(73)
where

$$\alpha _f=\left[\frac{2(m_1+1)}{b1}\right]^{(m_1+k4)/(m_1+1)}.$$
(74)
Similarly, the asymptotic form of $`h(\xi )`$ for large $`\xi `$ can be derived from equations (64) and (66). We obtain

$$h(\xi )\alpha _h\xi ^{(m_1+k2)/(m_1+1)},$$
(75)
where

$$\alpha _h=\left[\frac{2(m_1+1)}{1(b1)^{(m_1+1)(d2)/(2b)}}\right]^{(m_1+k2)/(m_1+1)}$$
(76)
if $`b2`$, and

$`\alpha _h`$ $`=`$ $`\left[{\displaystyle \frac{2(m_1+1)}{(b1)^{m_1+2}}}\right]^{(m_1+k2)/(m_1+1)}`$ (78)
$`\times \mathrm{exp}[(m_1+k2)(d2)]`$
if $`b=2`$.

Equation (53), which describes the evolution of the $`C_+`$ characteristic, can also be rewritten using the variable $`y`$. Using equation (44), we obtain

$$t\frac{dy_+}{dt}=\frac{(y_+y_1)(y_+b)}{y_+y_2}.$$
(79)
We have proven earlier that when $`k>k_c`$ one finds $`b>1`$ and $`y1`$. Thus the sign of the right-hand-side of the above equation is decided by the term $`(y_+y_1)`$. If a $`C_+`$ characteristic emerges from the region $`y<y_1`$, i.e., $`y_+(t=t_0)<y_1`$, we know that $`y_+(t)`$ will decrease with time or equivalently $`\xi _+(t)`$ will increase with time and the $`C_+`$ characteristic will not catch-up with the shock front.

Equations (55)–(57) can also be solved numerically for different values of $`k`$. We plot the results for $`k=5.5`$ ($`m_1=1.77`$) and $`k=6.5`$ ($`m_1=2.23`$) in Figure 1. When $`k_c<k<6`$, the function $`f(\xi )`$ decreases with increasing $`\xi `$ while the function $`g(\xi )`$ increases with increasing $`\xi `$. This implies that when moving inwards away from the shock front, the pressure decreases while the Lorentz factor increases. For $`k>6`$ the situation is reversed. For all $`k>k_c`$, the function $`h(\xi )`$ increases with $`\xi `$ implying that the density always increases when moving inwards away from the shock front.

## III Properties of the flow in the self-similar region

We now examine the properties of the flow in the self-similar region bounded by $`\xi =0`$ and $`\xi =\xi _+(t)`$, where $`\xi _+(t)`$ coincides with a $`C_+`$ characteristic that emerges from the region $`y<y_1`$. The evolution of the $`C_+`$ characteristic is described by equation (79) with the initial condition $`\xi _+(t=t_0)=\xi _0`$, and correspondingly $`y_+(t=t_0)=y_0<y_1`$. By solving equation (79), we get the following equation for $`y_+(t)`$,

$$\left(\frac{y_1y_+}{y_1y_0}\right)^{(y_2y_1)/(by_1)}\left(\frac{by_+}{by_0}\right)^{(by_2)/(by_1)}=\frac{t}{t_0}.$$
(80)
We can now derive the asymptotic behavior of the $`C_+`$ characteristic at $`t\mathrm{}`$. In this limit $`|y_+(t)|`$ is large, and equation (80) yields

$$y_+(t)\alpha _0t,$$
(81)
where

$$\alpha _0=(y_1y_0)^{(y_2y_1)/(by_1)}(by_0)^{(by_2)/(by_1)}/t_0.$$
(82)
Substituting equation (71) into equation (81), we get

$$\xi _+(t)\alpha _+t^{(m_1+1)},$$
(83)
where

$$\alpha _+=\left(\frac{\alpha _0}{\alpha _y}\right)^{(m_1+1)}.$$
(84)
The initial condition for the $`C_+`$ characteristic, namely the value of the $`C_+`$ characteristic with which $`\xi _+(t)`$ coincides, contains the information about the initial explosion.

The shock front is accelerating with a power-law temporal dependence of its Lorentz factor, $`\mathrm{\Gamma }^2t^{m_1}`$. How does the $`C_+`$ characteristic propagate? From equations (69) and (83) we find that when $`t\mathrm{}`$,

$$\gamma ^2(\xi _+)=\frac{1}{2}\mathrm{\Gamma }^2g(\xi _+)t^2.$$
(85)
We see that irrespective of the value of $`k`$, the $`C_+`$ characteristic always accelerates as $`\gamma ^2t^2`$. We can also calculate the thickness of the self-similar region. Using equations (5) and (83) we obtain when $`t\mathrm{}`$,

$$Rr_+=\frac{R\xi _+}{\mathrm{\Gamma }^2}\mathrm{const}.$$
(86)
Note that when $`2<m_1<m_{1c}`$, the $`C_+`$ characteristic accelerates faster than the shock front, but because the Lorentz factors of both surfaces are accelerating as power-laws of time, the $`C_+`$ characteristic can never catch-up with the shock front. Instead, the distance between the two surfaces approaches a constant value at late times.

We can now examine the energy and mass contained in the self-similar region. The energy contained in the spherical shell between $`\xi =0`$ and $`\xi =\xi _+(t)`$ is given by

$$E=_{r_+}^R16\pi p\gamma ^2r^2𝑑r=\frac{16\pi }{3}w_1\mathrm{\Gamma }^2R^3_0^{\xi _+}f(\xi )g(\xi )𝑑\xi .$$
(87)
Using equations (73) and (69), we can calculate the above integral for large values of $`\xi _+`$. This gives,

$`E`$ $``$ $`{\displaystyle \frac{16\pi }{3}}\left({\displaystyle \frac{m_1+1}{m_1+k3}}\right)\alpha _f\alpha _gw_1\mathrm{\Gamma }^2R^3`$ (89)
$`\times \xi _+^{(m_1+k3)/(m_1+1)}.`$
In deriving the above result we have used the fact that when $`k>k_c`$, we have $`m_1+k3>0`$ and $`m_1+1<0`$, so that the exponent of $`\xi _+`$ in the above equation $`(m_1+k3)/(m_1+1)`$ is positive. When $`t\mathrm{}`$, $`\xi _+`$ is given by equation (83). In addition $`\mathrm{\Gamma }^2t^{m_1}`$, $`w_1=\rho _1R^k`$. Thus when $`t\mathrm{}`$,

$$E\mathrm{const}.$$
(90)
The total number of particles contained between $`\xi =0`$ and $`\xi =\xi _+(t)`$ is given by

$$N=_{r_+}^Rn^{}4\pi r^2𝑑r=8\pi n_1R^3_0^{\xi _+}h(\xi )𝑑\xi .$$
(91)
Using equations (75) and (83), we obtain that for $`t\mathrm{}`$,

$`N`$ $`8\pi \left({\displaystyle \frac{m_1+1}{k3}}\right)\alpha _hn_1R^3\xi _+^{(k3)/(m_1+1)}`$ (93)
$`\mathrm{const}.`$
We have thus proven that both the energy and mass contained between the $`C_+`$ characteristic and the shock front will approach constant values as $`t\mathrm{}`$. The situation is similar to the non-relativistic case .

## IV Approximate (Analytic) Stability Analysis

In order to analyze the stability of the self-similar solutions obtained in §II.B, we first follow an analytic approach (in §IV) based on the assumptions of variable separation and a fixed $`\mathrm{\Gamma }_0`$, where $`\mathrm{\Gamma }_0`$ is the unperturbed Lorentz factor of the shock front. As we will explain later, these assumptions limit the generality of the results. We then use numerical simulations (in §V) to directly solve the evolution of the perturbations without those assumptions. The numerical simulations demonstrate that the results obtained using the analytic approach are qualitatively valid.

### A Derivation of linear perturbation equations

For the analytic approach to the stability analysis of the self-similar solutions, we use linear perturbation analysis similar to that used in the non-relativistic case . We start from the equations of motion for an ideal relativistic fluid:

$$\frac{}{t}(n^{})+(n^{}𝐯)=0,$$
(94)
$$\gamma ^2(e+p)\left[\frac{𝐯}{t}+(𝐯)𝐯\right]+\left(p+𝐯\frac{p}{t}\right)=0,$$
(95)
$$\frac{}{t}\left(\frac{\gamma ^{4/3}p}{n^{4/3}}\right)+(𝐯)\left(\frac{\gamma ^{4/3}p}{n^{4/3}}\right)=0,$$
(96)
where $`𝐯`$ and $`\gamma `$ are the fluid velocity and Lorentz factor in the perturbed solution measured in the laboratory frame, $`e`$ and $`p`$ are the energy density and pressure in the perturbed solution measured in the fluid frame, and $`n^{}`$ is the fluid density in the perturbed solution measured in the laboratory frame. We use the Eulerian perturbation approach, i.e., the perturbed quantities are defined as the difference between the perturbed solution and the unperturbed one in the same spatial point. Therefore we define the perturbed hydrodynamic quantities as

$$\delta p(r,\theta ,\varphi ,t)=p(r,\theta ,\varphi ,t)p_0(r,t),$$
(97)
$$\delta 𝐯(r,\theta ,\varphi ,t)=𝐯(r,\theta ,\varphi ,t)v_0(r,t)\widehat{𝐫}=\delta v_r\widehat{𝐫}+\delta 𝐯_T,$$
(98)
$$\delta n^{}(r,\theta ,\varphi ,t)=n^{}(r,\theta ,\varphi ,t)n_0^{}(r,t),$$
(99)
where the quantities with subscript “0” are the unperturbed values. Substituting the above quantities into equations (94)–(96), we obtain the following linear perturbation equations,

$`{\displaystyle \frac{}{t}}\delta n^{}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\{r^2[{\displaystyle \frac{1}{2\gamma _0^4}}(1+{\displaystyle \frac{1}{2\gamma _0^2}})n_o^{}\delta \gamma ^2`$ (100)
$`+(1{\displaystyle \frac{1}{2\gamma _0^2}})\delta n^{}]\}+n^{}_0_T\delta 𝐯_T=0,`$ (101)
$`{\displaystyle \frac{2}{\gamma _0^4}}\left[\left(1+{\displaystyle \frac{1}{2\gamma _0^2}}\right){\displaystyle \frac{\gamma _0^2}{t}}+{\displaystyle \frac{\gamma _0^2}{r}}\right](\gamma _0^2\delta p+p_0\delta \gamma ^2)`$ (102)
$`+4\gamma _0^2p_0\{{\displaystyle \frac{}{t}}\left[{\displaystyle \frac{1}{2\gamma _0^4}}(1+{\displaystyle \frac{1}{2\gamma _0^2}})\delta \gamma ^2\right]`$ (103)
$`+\left(1{\displaystyle \frac{1}{2\gamma _0^2}}\right){\displaystyle \frac{}{r}}\left[{\displaystyle \frac{1}{2\gamma _0^4}}\left(1+{\displaystyle \frac{1}{2\gamma _0^2}}\right)\delta \gamma ^2\right]`$ (104)
$`+{\displaystyle \frac{1}{4\gamma _0^8}}(1+{\displaystyle \frac{1}{\gamma _0^2}})\left({\displaystyle \frac{\gamma _0^2}{r}}\right)\delta \gamma ^2\}`$ (105)
$`+{\displaystyle \frac{}{r}}\delta p+\left(1{\displaystyle \frac{1}{2\gamma _0^2}}\right){\displaystyle \frac{}{t}}\delta p`$ (106)
$`+{\displaystyle \frac{1}{2\gamma _0^4}}\left(1+{\displaystyle \frac{1}{2\gamma _0^2}}\right)\left({\displaystyle \frac{p_0}{t}}\right)\delta \gamma ^2=0,`$ (107)
$`4p_0\gamma _0^2[{\displaystyle \frac{}{t}}\delta 𝐯_T+(1{\displaystyle \frac{1}{2\gamma _0^2}}){\displaystyle \frac{}{r}}\delta 𝐯_T`$ (108)
$`+(1{\displaystyle \frac{1}{2\gamma _0^2}}){\displaystyle \frac{1}{r}}\delta 𝐯_T]+_T\delta p+({\displaystyle \frac{p_0}{t}})\delta 𝐯_T=0,`$ (109)
$`{\displaystyle \frac{}{t}}\left[{\displaystyle \frac{\gamma _0^{4/3}p_0}{n_0^{4/3}}}\left({\displaystyle \frac{2}{3}}{\displaystyle \frac{\delta \gamma ^2}{\gamma _0^2}}+{\displaystyle \frac{\delta p}{p_0}}{\displaystyle \frac{4}{3}}{\displaystyle \frac{\delta n^{}}{n_0^{}}}\right)\right]`$ (110)
$`+\left(1{\displaystyle \frac{1}{2\gamma _0^2}}\right){\displaystyle \frac{}{r}}\left[{\displaystyle \frac{\gamma _0^{4/3}p_0}{n_0^{4/3}}}\left({\displaystyle \frac{2}{3}}{\displaystyle \frac{\delta \gamma ^2}{\gamma _0^2}}+{\displaystyle \frac{\delta p}{p_0}}{\displaystyle \frac{4}{3}}{\displaystyle \frac{\delta n^{}}{n_0^{}}}\right)\right]`$ (111)
$`+{\displaystyle \frac{1}{2\gamma _0^4}}\left(1+{\displaystyle \frac{1}{2\gamma _0^2}}\right){\displaystyle \frac{}{r}}\left({\displaystyle \frac{\gamma _0^{4/3}p_0}{n_0^{4/3}}}\right)\delta \gamma ^2=0,`$ (112)
where we have used the relations

$$\frac{v_0}{t}=\frac{1}{2\gamma _0^4}\left(1+\frac{1}{2\gamma _0^2}\right)\frac{\gamma _0^2}{t},$$
(113)
$$\frac{v_0}{r}=\frac{1}{2\gamma _0^4}\left(1+\frac{1}{2\gamma _0^2}\right)\frac{\gamma _0^2}{r},$$
(114)
$$\delta v_r=\frac{1}{2\gamma _0^4}\left(1+\frac{1}{2\gamma _0^2}\right)\delta \gamma ^2,$$
(115)
and the operator $`_T(\widehat{\theta }/r)(/\theta )+(\widehat{\varphi }/r)(/\varphi )`$ acts as follows on a scalar $`\mathrm{\Psi }`$ and a vector $`𝐟`$:

$$_T\mathrm{\Psi }=\frac{1}{r}\frac{\mathrm{\Psi }}{\theta }\widehat{\theta }+\frac{1}{r\mathrm{sin}\theta }\frac{\mathrm{\Psi }}{\varphi }\widehat{\varphi },$$
(116)
$$_T𝐟=\frac{1}{r\mathrm{sin}\theta }\frac{}{\theta }(\mathrm{sin}\theta f_\theta )+\frac{1}{r\mathrm{sin}\theta }\frac{f_\varphi }{\varphi }.$$
(117)
Since the unperturbed quantities satisfy equations (20)–(22), we write

$$p_0=\frac{2}{3}\rho _1\mathrm{\Gamma }_0^2f(\xi ),$$
(118)
$$\gamma _0^2=\frac{1}{2}\mathrm{\Gamma }_0^2g(\xi ),$$
(119)
$$n_0^{}=2n_1\mathrm{\Gamma }_0^2h(\xi ),$$
(120)
where $`\mathrm{\Gamma }_0`$ is the unperturbed Lorentz factor of the shock front, and $`\xi `$ is the similarity variable defined as

$$\xi =\left(1\frac{r}{R_0}\right)\mathrm{\Gamma }_0^2,$$
(121)
where $`R_0`$ is the unperturbed radius of the shock front. We further define the perturbation variables as

$$\delta \gamma ^2(r,\theta ,\varphi ,t)=\frac{1}{2}\mathrm{\Gamma }_0^2\delta g(\xi )Y_{lm}(\theta ,\varphi )X(t),$$
(122)
$$\delta 𝐯_T(r,\theta ,\varphi ,t)=\frac{1}{\mathrm{\Gamma }_0^2}\delta v_T(\xi )\stackrel{~}{}_TY_{lm}(\theta ,\varphi )X(t),$$
(123)
$$\delta p(r,\theta ,\varphi ,t)=\frac{2}{3}\rho _1\mathrm{\Gamma }_0^2\delta f(\xi )Y_{lm}(\theta ,\varphi )X(t),$$
(124)
$$\delta n^{}(r,\theta ,\varphi ,t)=2n_1\mathrm{\Gamma }_0^2\delta h(\xi )Y_{lm}(\theta ,\varphi )X(t),$$
(125)
where the operator $`\stackrel{~}{}_T\widehat{\theta }(/\theta )+\widehat{\varphi }(1/\mathrm{sin}\theta )(/\varphi )`$. Note that the variables $`\xi `$ and $`t`$ are separated in above definitions of the perturbations, and so we consider only “global” perturbations . The function $`X(t)`$ measures the amplitude of the perturbation relative to the unperturbed values.

Substituting equations (118)–(120) and (122)–(125) into equations (101)–(112), we obtain

$`\left[q+2km_1+{\displaystyle \frac{2(m_1+2)}{(yy_2)}}\right]{\displaystyle \frac{\delta h}{h}}{\displaystyle \frac{(y2)}{2}}\left({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{h}}{\displaystyle \frac{d\delta h}{d\xi }}\right)`$ (126)
$`+\left[{\displaystyle \frac{4(m_1+2)}{(yy_2)}}+{\displaystyle \frac{2(m_1+k2)(yd)}{(yy_2)(y2)}}\right]{\displaystyle \frac{\delta g}{g}}`$ (127)
$`\left({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{g}}{\displaystyle \frac{d\delta g}{d\xi }}\right)+{\displaystyle \frac{l(l+1)}{\mathrm{\Gamma }_0^2}}\delta v_T=0,`$ (128)
$`\left[k3m_1+q+{\displaystyle \frac{2(m_1+2)(y2)}{(yy_2)}}\right]{\displaystyle \frac{\delta f}{f}}`$ (129)
$`{\displaystyle \frac{(y+2)}{2}}\left({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{f}}{\displaystyle \frac{d\delta f}{d\xi }}\right)`$ (130)
$`+\left[2q{\displaystyle \frac{2(m_1+2)(y4)}{(yy_2)}}{\displaystyle \frac{2(m_1+k4)}{(yy_2)}}\right]{\displaystyle \frac{\delta g}{g}}`$ (131)
$`(y2)\left({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{g}}{\displaystyle \frac{d\delta g}{d\xi }}\right)=0,`$ (132)
$`\left[2(m_1+q+1){\displaystyle \frac{2(m_1+k4)}{(yy_2)}}\right]g\delta v_T`$ (133)
$`(y2){\displaystyle \frac{d\delta v_T}{d\xi }}{\displaystyle \frac{\delta f}{f}}=0,`$ (134)
$`[{\displaystyle \frac{2}{3}}q{\displaystyle \frac{2}{3}}{\displaystyle \frac{(m_1+2)(y4)}{(yy_2)}}+{\displaystyle \frac{2(m_1+k4)}{(yy_2)}}`$ (135)
$`{\displaystyle \frac{8}{3}}{\displaystyle \frac{(m_1+k2)(yd)}{(yy_2)(y2)}}]{\displaystyle \frac{\delta g}{g}}{\displaystyle \frac{(y2)}{3}}({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{g}}{\displaystyle \frac{d\delta g}{d\xi }})`$ (136)
$`+\left[q{\displaystyle \frac{(m_1+k4)(y2)}{(yy_2)}}\right]{\displaystyle \frac{\delta f}{f}}{\displaystyle \frac{(y2)}{2}}\left({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{f}}{\displaystyle \frac{d\delta f}{d\xi }}\right)`$ (137)
$`{\displaystyle \frac{4}{3}}\left[q{\displaystyle \frac{(m_1+k2)(yd)}{(yy_2)}}\right]{\displaystyle \frac{\delta h}{h}}`$ (138)
$`+{\displaystyle \frac{2}{3}}(y2)\left({\displaystyle \frac{1}{g}}{\displaystyle \frac{1}{h}}{\displaystyle \frac{d\delta h}{d\xi }}\right)=0,`$ (139)
where we have used equations (55)–(57) and the following relations

$$\frac{d\mathrm{\Gamma }_0^2}{dt}=m_1\frac{\mathrm{\Gamma }_0^2}{t},$$
(140)
$$\frac{\xi }{r}=\frac{1}{t}\left[\mathrm{\Gamma }_0^2+\frac{1}{2(m_1+1)}\right],$$
(141)
$$\frac{\xi }{t}=\frac{1}{t}\left[\mathrm{\Gamma }_0^2\frac{m_1}{2(m_1+1)}\xi (m_1+1)\right],$$
(142)
$$\frac{d\rho _1}{dt}=k\frac{\rho _1}{t}\left[1\frac{m_1}{2(m_1+1)\mathrm{\Gamma }_0^2}\right],$$
(143)
$$\frac{dn_1}{dt}=k\frac{n_1}{t}\left[1\frac{m_1}{2(m_1+1)\mathrm{\Gamma }_0^2}\right],$$
(144)
$$r=t\left[1\frac{1}{2(m_1+1)\mathrm{\Gamma }_0^2}\frac{\xi }{\mathrm{\Gamma }_0^2}\right].$$
(145)
Note that in deriving equations (128)–(139) we have assumed that there is no perturbation in the external medium. Moreover, in order to separate variables, $`X(t)`$ has to be a power law in time, $`X(t)t^q`$, where $`q`$ defines the temporal evolution of the perturbation amplitude. If the real part of $`q`$ is positive then the perturbation grows, while if the real part of $`q`$ is negative then the perturbation decays.

In equation (128), the term $`l(l+1)/\mathrm{\Gamma }_0^2`$ is associated with causality, namely the fact that a perturbation can only propagate at a speed $`c/\mathrm{\Gamma }_0`$ in the transverse direction and hence expand across a maximum opening angle of $`1/\mathrm{\Gamma }_0`$. Since $`\mathrm{\Gamma }_0`$ is a function of time, it is not possible to achieve a complete separation of variables for this equation in contrast with the non-relativistic case. However, for any constant value of $`\mathrm{\Gamma }_0`$ we can still calculate the power-law index for the growth of the perturbation, $`q`$. These results are meaningful if we find $`q>|m_1|`$, so that perturbations grow on a time scale shorter than the time scale for changes in $`\mathrm{\Gamma }_0`$. Therefore, the assumptions of variable separation and fixed $`\mathrm{\Gamma }_0`$ limit the generality of the results. However, even if we find $`q<|m_1|`$, we should still be able to gain an insight into some qualitative properties of the perturbation amplitude evolution.

Equations (128)–(139) are a complete set of first-order differential equations for $`\delta f`$, $`\delta g`$, $`\delta h`$, $`\delta v_T`$. After some algebraic manipulations, one may write the equations for the first order terms $`d\delta f/d\xi `$, $`d\delta g/d\xi `$, $`d\delta h/d\xi `$ and $`d\delta v_T/d\xi `$ in the following matrix form

$`(y^28y+4){\displaystyle \frac{d}{d\xi }}\left(\begin{array}{c}\delta f\\ \delta g\\ \delta h\\ \delta v_T\end{array}\right)`$ (150)
$`=𝐀(q,k,l(l+1)/\mathrm{\Gamma }_0^2,\xi )\left(\begin{array}{c}\delta f\\ \delta g\\ \delta h\\ \delta v_T\end{array}\right),`$ (155)
where $`𝐀`$ is a $`4\times 4`$ matrix. Note that $`(y^28y+4)=(yy_1)(yy_2)`$. Thus, the solutions for the perturbation variables must pass the same singular point (or the sonic line), $`y_1=42\sqrt{3}`$, as the unperturbed variables. Therefore, the value of $`q`$ can be found by requiring that the solutions pass through the singular point $`y_1`$. This is very similar to the non-relativistic case .

In order to numerically integrate the differential equations (155) and derive $`q`$ we need to specify the boundary conditions at the shock front when the shock is perturbed. Since the relativistic jump conditions across the shock front must be satisfied, we have

$$p=\frac{2}{3}\mathrm{\Gamma }^2\rho _1,$$
(156)
$$n^{}=2\mathrm{\Gamma }^2n_1,$$
(157)
$$\gamma ^2=\frac{1}{2}\mathrm{\Gamma }^2,$$
(158)
where $`\mathrm{\Gamma }`$ is the Lorentz factor of the perturbed shock front. By linearizing these boundary conditions with respect to the perturbed quantities, we find

$$\delta p+\left(\frac{p_0}{r}\right)\delta R=\frac{4}{3}\mathrm{\Gamma }_0^4\rho _1\frac{d}{dt}\delta R\frac{2k}{3}\mathrm{\Gamma }_0^2\rho _1\frac{\delta R}{R_0},$$
(159)
$$\delta n^{}+\left(\frac{n_0^{}}{r}\right)\delta R=4\mathrm{\Gamma }_0^4n_1\frac{d}{dt}\delta R2k\mathrm{\Gamma }_0^2n_1\frac{\delta R}{R_0},$$
(160)
$$\delta \gamma ^2+\left(\frac{\gamma _0^2}{r}\right)\delta R=\mathrm{\Gamma }_0^4\frac{d}{dt}\delta R,$$
(161)
where $`\delta R`$ is the deviation of the perturbed shock radius $`R`$ from the unperturbed shock radius $`R_0`$. In deriving equations (159)–(161), we used the relations $`\rho _1R^k`$ and

$$\delta \mathrm{\Gamma }^2=2\mathrm{\Gamma }_0^4\frac{d}{dt}\delta R,$$
(162)
where $`\delta \mathrm{\Gamma }^2`$ is the deviation of the square of the perturbed shock Lorentz factor $`\mathrm{\Gamma }^2`$ from the square of the unperturbed shock Lorentz factor $`\mathrm{\Gamma }_0^2`$.

We now define

$$\delta R(\theta ,\varphi ,t)=\eta \frac{1}{\mathrm{\Gamma }_0^2}R_0Y_{lm}(\theta ,\varphi )X(t),$$
(163)
where $`\eta `$ is a scale factor that can have an arbitrary value; for convenience we set $`\eta =1`$. Substituting equations (163), (118)–(120), (122), (124) and (125) into equations (159)–(161), we find

$$\delta f=\frac{df}{d\xi }+2(m+q+1),$$
(164)
$$\delta h=\frac{dh}{d\xi }+2(m+q+1),$$
(165)
$$\delta g=\frac{dg}{d\xi }+2(m+q+1).$$
(166)
Another boundary condition results from the requirement that the tangential velocities must be continuous across the shock front, yielding

$$\delta 𝐯_T=\frac{1}{R_0}\stackrel{~}{}_T\delta R.$$
(167)
Substituting equations (123) and (163) into this equation, we get

$$\delta v_T=1.$$
(168)
Equations (164)–(166) and (168) are the four boundary conditions necessary to solve the perturbation equations.

### B Numerical results

Based on the derivations presented in the previous subsection, we may now examine the stability of different modes for different values of $`k`$. As a particular example, we consider the case of $`k=5.5`$. We derive $`q`$ for different values of the mode wavenumber $`l`$ by integrating the differential equations (155) from the shock front to its interior and requiring that the solutions pass through the singular point. The results are shown in Figure 2. The top panel shows the real component of $`q`$ ($`\mathrm{Re}[q]`$) as a function of $`\sqrt{l(l+1)}/\mathrm{\Gamma }_0`$, where $`\mathrm{\Gamma }_0`$ is treated as a scaling factor. As mentioned before, $`\mathrm{Re}[q]`$ determines the growth rate of the perturbation. In the bottom panel, we plot the imaginary component of $`q`$ ($`\mathrm{Im}[q]`$), which provides the oscillation frequency of the perturbation. Figure 2 separates the behavior of $`q`$ into three different regimes:

* In the regime of small $`l/\mathrm{\Gamma }_0`$ ($`0<l/\mathrm{\Gamma }_0<0.87`$), $`q`$ is a real number and $`\mathrm{Re}[q]`$ is positive, implying that the perturbation grows monotonically in time. The value of $`\mathrm{Re}[q]`$ increases as $`l`$ increases. Note that $`q`$ vanishes in the limit of $`l0`$. This result can be derived analytically by comparing two unperturbed spherical solutions with different parameters.
* In the regime of intermediate $`l/\mathrm{\Gamma }_0`$ ($`0.87<l/\mathrm{\Gamma }_0<17`$), $`q`$ is a complex number and $`\mathrm{Re}[q]`$ is positive, implying that the perturbation grows while oscillating. As $`l`$ increases the real part $`\mathrm{Re}[q]`$ decreases while the imaginary part $`\mathrm{Im}[q]`$ increases. Note that the transition between the real and imaginary solutions for $`q`$ occurs at $`l/\mathrm{\Gamma }_00.87`$. This result follows from causality. When the wavelength of the perturbation ($`1/l`$) is smaller than $`1/\mathrm{\Gamma }_0`$, the maximum angular separation of two regions that can interact with each other, the perturbation can oscillate.
* Finally in the regime of large $`l/\mathrm{\Gamma }_0`$ ($`l/\mathrm{\Gamma }_0>17`$), $`q`$ is a complex number and $`\mathrm{Re}[q]`$ is negative, implying that the perturbation decays while oscillating. The value of $`\mathrm{Re}[q]`$ decreases (so the absolute value of $`\mathrm{Re}[q]`$ increases) as $`l`$ increases while the value of $`\mathrm{Im}[q]`$ increases as $`l`$ increases.

The actual evolution of a perturbation is shaped by the fact that $`\mathrm{\Gamma }_0`$ increases with time as $`\mathrm{\Gamma }_0^2t^{m_1}`$. If initially the wavenumber of the perturbation is sufficiently large so that it is in the regime of large $`l/\mathrm{\Gamma }_0`$, the perturbation will start to decay while oscillating. As time progresses, $`l/\mathrm{\Gamma }_0`$ decreases and so both $`\mathrm{Re}[q]`$ and $`\mathrm{Im}[q]`$ decrease, the perturbation decays with slower speed and oscillates on longer timescales. As soon as the perturbation enters the regime of intermediate $`l/\mathrm{\Gamma }_0`$, it starts to grow slowly over time and oscillate on even longer timescales. The growth rate slowly increases over time, but is always limited by the rather small upper bound, $`\mathrm{Re}[q]0.38`$. Eventually, the perturbation enters the regime of small $`l/\mathrm{\Gamma }_0`$ and grows slowly without oscillating. As $`t`$ increases, the growth rate approaches zero, and so the perturbation saturates. Therefore, perturbations with large wavenumbers (short wavelengthes) grow when $`1l/\mathrm{\Gamma }_010`$ only by a modest factor. In the case of intermediate wavenumbers, the perturbation goes through the two regimes of intermediate and small $`l/\mathrm{\Gamma }_0`$. Therefore it grows slowly with some initial oscillations, but soon afterwards it stops oscillating and saturates. Perturbations with small wavenumbers stay in the regime of small $`l/\mathrm{\Gamma }_0`$. The perturbation grows slowly without oscillating at the beginning but soon saturates.

The above results remain qualitatively the same for all values of $`k>4.134`$.

## V Numerical Simulations

We have verified the above behavior by a direct integration of the partial different equations which determine the evolution of the perturbation variables, without assuming separability of the solutions with respect to $`\xi `$ and $`t`$. Instead of equations (122) – (125), we redefined the perturbation variables as

$$\delta \gamma ^2(r,\theta ,\varphi ,t)=\frac{1}{2}\mathrm{\Gamma }_0^2\delta g(\xi ,t)Y_{lm}(\theta ,\varphi ),$$
(169)
$$\delta 𝐯_T(r,\theta ,\varphi ,t)=\frac{1}{\mathrm{\Gamma }_0^2}\delta v_T(\xi ,t)\stackrel{~}{}_TY_{lm}(\theta ,\varphi ),$$
(170)
$$\delta p(r,\theta ,\varphi ,t)=\frac{2}{3}\rho _1\mathrm{\Gamma }_0^2\delta f(\xi ,t)Y_{lm}(\theta ,\varphi ),$$
(171)
$$\delta n^{}(r,\theta ,\varphi ,t)=2n_1\mathrm{\Gamma }_0^2\delta h(\xi ,t)Y_{lm}(\theta ,\varphi ).$$
(172)
Equations (128) – (139) were then replaced by four partial differential equations (PDEs) for the perturbation variables $`\delta f(\xi ,t)`$, $`\delta g(\xi ,t)`$ $`\delta h(\xi ,t)`$ and $`\delta v_T(\xi ,t)`$. We then solved for the evolution of these perturbation variables by numerically integrating the PDEs with appropriate initial values. In our numerical simulations, the outer boundary ($`\xi =0`$) is the shock front where the shock jump conditions are assumed to be satisfied. We can still define $`\delta R`$ as in equation (163). Then at the outer boundary the perturbation variables satisfy

$$\delta f(0,t)=\frac{df}{d\xi }X(t)+2(m_1+1)X(t)+2t\frac{dX(t)}{dt},$$
(173)
$$\delta h(0,t)=\frac{dh}{d\xi }X(t)+2(m_1+1)X(t)+2t\frac{dX(t)}{dt},$$
(174)
$$\delta g(0,t)=\frac{dg}{d\xi }X(t)+2(m_1+1)X(t)+2t\frac{dX(t)}{dt},$$
(175)
$$\delta v_T(0,t)=X(t).$$
(176)
The inner boundary is chosen to be sufficiently large so as to cover the entire similarity region which is bounded by a inner $`C_+`$ characteristic. This way, the values of the perturbation variables at the inner boundary can not affect the shock front.

The numerical simulations gave us the same behavior for the perturbations as the previously mentioned analytical results for the growth rate $`q`$. In Figure 3, we show the evolution of $`X(t)`$ in a numerical simulation with $`k=5.5`$ and $`ł/\mathrm{\Gamma }_0=75`$. We plot $`X(t)`$ in the range $`t[1,10]`$ in the top panel and $`t[1,100]`$ in the bottom panel. Note that $`X(t)`$ describes the relative displacement of the perturbed shock radius from the unperturbed value and we set its initial value to be $`X(t=1)=1.0`$. From Figure 3 we see that $`X`$ oscillates over time over increasingly longer timescales and stops oscillating at late times. Its amplitude first decreases, then slowly increases and finally saturates. Overall it grows by a factor of $`10`$. These results are consistent with the previous discussion on the three regimes for the evolution of the perturbations.

## VI Conclusions

We have derived the self-similar solutions for an ultra-relativistic blast wave in an external medium with a density profile $`\rho _1r^k`$ and $`k>4`$. The solutions exist for $`k`$ larger than a critical value $`k_c=4.134`$. They describe the flow in the self-similar region bounded by the shock front and a $`C_+`$ characteristic. The shock front accelerates with Lorentz factor $`\mathrm{\Gamma }^2t^{m_1}`$ and $`m_1<1.134`$, while the $`C_+`$ characteristic accelerates with Lorentz factor $`\gamma ^2t^2`$. The energy and mass contained inside the self-similar region approach constant values as time diverges.

We have found that at large wavenumbers the perturbations first decay, then grow slowly over time and eventually saturate. The initial decay and the intermediate growth are accompanied by temporal oscillations. These small wavelength perturbations grow when $`1l/\mathrm{\Gamma }_010`$ with an overall factor of $`10`$. At intermediate wavenumbers, the perturbations first grow slowly and then saturate. The initial growth is also accompanied by temporal oscillations. At small wavenumbers the perturbations grow monotonically in time but soon saturate. Our results also apply to expanding relativistic jets as long as the opening angle of the jet is larger than the inverse of its Lorentz factor.

In the collapsar model of gamma-ray bursts, a collimated relativistic outflow is generated due to the collapse of the core of a massive star. The outflow approaches the stellar envelope at a modest semi-relativistic speed but is expected to accelerate significantly across the sharp density gradient at the surface of the star. Our results indicate that in the breakout phase perturbations are close to being stable in spherical symmetry. It is still possible, however, that the lateral expansion of the jet at breakout would be accompanied by instabilities. These instabilities may produce variations in the Lorentz factor of the jet needed in the internal shock model. They may also be responsible for the complex light curves observed in most GRBs. Current numerical simulations lack adequate resolution at the stellar surface to follow the shock breakout and confirm the instabilities. We leave a detailed study of the instabilities associated with the lateral expansion of the jet for future work.

#### Acknowledgments.

This work was supported in part by grants from the Israel-US BSF (BSF-9800343) and NSF (AST-0071019, AST-0204514).