REVISITING THE EMISSION FROM RELATIVISTIC BLAST WAVES IN A DENSITY-JUMP MEDIUM

It is believed that gamma-ray bursts (GRBs) are produced by either the collapse of massive stars (Woosley 1993; MacFadyen & Woosley 1999) or the merger of compact binaries (Eichler et al. 1989), which can launch a collimated relativistic blast wave into circumburst medium. As the relativistic forward shock propagates into the surrounding medium, the circumburst electrons would be swept up and accelerated. These shocked electrons power the afterglow emission via synchrotron radiation, which can be observed for weeks at X-ray and optical frequencies. In the pre-Swift era, many afterglow light curves showed a smooth power-law decay. On the basis of the self-similar Blandford–McKee solution (BM; Blandford & McKee 1976) to the dynamics of the relativistic blast wave, a standard model of forward shock afterglow was developed (e.g., Mészáros & Rees 1997; Panaitescu et al. 1998; Sari et al. 1998; see Zhang & Mészáros 2004 and Gao et al. 2013 for recent reviews) and is generally consistent with the observations. However, some optical afterglows have shown significant temporal variabilities, with a strong deviation from the typical power-law behavior (Frontera et al. 2000). This issue has long been debated among researchers.

Recent observations by Swift/XRT (Gehrels et al. 2004) reveal that early X-ray flares are common in afterglows. These early flares are probably caused by the early activities of the central source (Burrows et al. 2005; Nousek et al. 2006). An overview of the optical afterglow samples (Liang et al. 2013) shows that many afterglows have optical bumps at a late time (usually 10⁴–10⁵ s after the trigger), while there are no corresponding significant signals at the X-ray band. Several scenarios have been proposed to explain such temporal variabilities in the GRB afterglow light curves. One proposal involves the "refreshed shocks," including the late energy injection model (Zhang & Mészáros 2002; Kong & Huang 2010; Geng et al. 2013) or the later internal shock model (Rees & Mészáros 1994; Fan & Wei 2005; Yu & Dai 2009; Wang & Cheng 2012). Notable bumps with a steep rising slope in the optical afterglow of some GRBs, like GRB 081029 (Nardini et al. 2011a) and GRB 100621A (Greiner et al. 2013), may be caused by this mechanism. The two-component jet model can also cause a late bump (Huang et al. 2006), and it has been applied to some cases (e.g., GRB 030329; Berger et al. 2003). Alternatively, a possible origin for the variabilities in the afterglows is a blast wave encountering a jump in density.

Density fluctuations near the GRB progenitor are expected because of strong winds and violent mass eruptions (Ramirez-Ruiz et al. 2001; Dai & Lu 1999) prior to the collapse of its core. Thus, many authors studied the interaction of a blast wave with density structures to check whether the density-jump (or decrease) can explain the complex variations in the afterglows. Dai & Lu (2002) calculated the emission when a relativistic blast wave interacts with a density-jump medium. In their analytical solution, a notable bump appears in the afterglow if a resultant reverse shock is relativistic. Dai & Wu (2003) applied this analytical solution to GRB 030226 to explain a rebrightening bump in its GRB afterglow. Other detailed numerical calculations (Lazzati et al. 2002; Mesler et al. 2012; Tam et al. 2005) show similar results and are applicable to some specific afterglows. However, the results from van Eerten et al. (2009, 2010) and Nakar & Granot (2007) are different—they suggest that the density-jump does not produce sharp flares in the light curves. Most recently, Gat et al. (2013) gave an analytical solution and also operated a numerical simulation to this problem, confirming the study of van Eerten et al. (2009). Thus the contradiction between the analytical result and the full hydrodynamics simulation motivated us to revisit this issue with a further semi-analytic numerical approach. In Uhm & Zhang (2014), it is concluded that the forward shock emission is not likely to reproduce significant re-brightening features (Uhm & Beloborodov 2007), and the reverse shock emission responds sensitively to the density structures in the medium (such as bumps, voids, or steps).

The dynamics of a forward shock propagating in a smooth density profile environment can be well described by the differential equations proposed by Huang et al. (1999, 2000a, 2000b). In this paper, we extend the derivation of Huang et al. (1999) to include the role of the reverse shock emerging during the interaction. Our results can help to explore whether the onset of the emission from the reverse shock can lead to the late bumps in the optical light curves. Our paper is organized as follows. In Section 2, we analyze the hydrodynamics of the blast wave expanding in the density-jump medium. In Section 3, we present the formula for the synchrotron emissivity. We present our numerical results and the comparison with previous work in Section 4. Our conclusions are summarized in Section 5. Details of our derivations for the dynamic equations are described in Appendices A and B.

After the prompt emission of GRBs, the blast wave will sweep up the ambient medium. We assume that the medium has a simple number density profile as a step function of radius R: n = n₀ for R ⩽ R₀ and n = n₁ for R ⩾ R₀ (n₁ > n₀), where R₀ is the transition radius. If the density contrast is higher than some threshold, the interaction of the forward shock with the high-density medium will lead to a reverse shock propagating into the hot shell (Sari & Piran 1995; Kobayashi et al. 1999; Kobayashi & Sari 2000). Thus the overall evolution of the blast wave should be divided into three episodes: the stage before encountering the density-jump, the period when the reverse shock is crossing the shell, and that after the reverse shock crosses the shell. We will carefully analyze these episodes in the following sections.

2.1. Before Encountering the Density Jump

Let's consider a forward shock propagating into the cold interstellar medium (ISM). The evolution of the forward shock is calculated by considering energy conservation. We assume the initial mass of the outflow is M_ej. The shock front separates the system into two regions: (1) the unshocked ISM, and (2) the shocked ISM. We treat these regions as simple homogenous shells (see Piran 1999 for an exhaustive treatment of this approach). If the outflow propagates with Lorentz factor γ₂ at time t, the thermodynamical quantities of the gas in the rest frame of region 2 are given by: $U_2^{\prime }=\psi _2 m_2 c^2=(\gamma _2-1)m_2 c^2$ (internal energy) and $P_2^{\prime }=(\hat{\gamma }-1) U_2^{\prime } =({\psi _2+2}/{3(1+\psi _2)})\psi _2m_2c^2$ (pressure), where m₂ is the total mass swept up by the shell and $\hat{\gamma }$ is the adiabatic index. Here and below, we use prime (') on variables to denote the quantities in the shock comoving frame and characters without a prime to denote quantities in the observer frame. If a fraction ₂ (radiation efficiency) of the thermal energy is radiated, then the "total" energy (not including the rest mass energy) of the gas is given by $E_2 \simeq (\gamma _2-1)(m_2+M_{\rm ej})c^2+(1-\epsilon _2)\gamma _2(U_2^{\prime }+P_2^{\prime })$. In fact, a more exact expression for the total energy is $E_2 = \lbrace \gamma _2+(1-\epsilon _2) [\hat{\gamma }\gamma _2^2-\gamma _2-(\hat{\gamma }-1)(1+\gamma _2\beta _2^2)]\rbrace m_2 c^2$ (Pe'er 2012), where β₂ is the velocity of region 2 and c is the light velocity. We use the approximation just for the simplification of the equations derived below, and this is especially helpful for the treatment of the episode in the next subsection. Because the shock-accelerated electrons carry only a fraction _e of the internal energy, the radiation efficiency of the total ejecta can be given by $\epsilon _2=\epsilon _e {t_{\rm syn}^{\prime }}^{-1}/({t_{\rm syn}^{\prime }}^{-1}+{t_{\rm ex}^{\prime }}^{-1})$ (Dai et al. 1999), where _e is the equipartition parameter for electrons, $t_{\rm syn}^{\prime }$ is the synchrotron cooling time, and $t_{\rm ex}^{\prime }=R/(\gamma _2c)$ is the comoving frame expansion time.

Following the procedure of Huang et al. (1999), we can derive the differential equation for the evolution of the Lorentz factor γ₂. When the ISM of mass dm₂ is swept up by the shell, a fraction ₂ of the gained thermal energy is assumed to be radiated, thus the lost energy is ${d}E=-\epsilon _2\gamma _2\hat{\gamma }(\gamma _2-1){d}m_2c^2.$ Substituting E with the formula above, we have

$$\begin{equation} \displaystyle {\frac{{d}\gamma _2}{{d}m_2}=-\frac{4\left(\gamma _2^2-1\right)}{8(1-\epsilon _2)\gamma _2m_2+3\epsilon _2m_2+3M_{\rm {ej}}}}. \end{equation} \tag{ 1 }$$

We have compared our Equation (1) with that of Pe'er (2012). There is almost no difference in the evolution of γ₂. So this equation can satisfactorily describe the fireball before the density-jump. Most recently, Nava et al. (2013) have proposed a realistic description for the dynamic evolution by including the adiabatic losses in a self-consistent way. However, their calculation is slightly more complicated when taking the reverse shock into consideration. We prefer to adopt the simple one here, which would not lead to significant errors on the results.

2.2. When the Reverse Shock Crosses the Shell

The interaction of the blast wave with the high-density medium is ascribed to two shocks: a reverse shock that propagates into the hot shell and a forward shock that propagates into the high-density medium. Thus the system consists of four separated regions: (1) unshocked high-density medium, (2) forward-shocked high-density medium, (3) reverse-shocked hot shell, and (4) unshocked hot shell. In this paper, the quantities (Lorentz factor, density, mass, internal energy, and pressure) of region "i" (regions 2, 3, and 4) are denoted by subscripts "i," respectively. γ_ij and β_ij are the relative Lorentz factor and dimensionless velocity of region "i" as measured in the rest frame of region "j." Similar to the derivation in Section 2.1, we have the total energy of three regions

$$\begin{equation} E_{\rm t} = E_2 + E_3 + E_4, \end{equation} \tag{ 2 }$$

$$\begin{equation} E_i = (\gamma _i-1) m_i c^2 + (1-\epsilon _i) \gamma _i \left(U_i^{\prime } + P_i^{\prime }\right), \quad i=2,3,4. \end{equation} \tag{ 3 }$$

With some derivations (see Appendix A), we can get the equation for the evolution of the Lorentz factor of regions 2 and 3 (γ₂ = γ₃):

$$\begin{equation} \displaystyle {\displaystyle\frac{{d}\gamma _2}{{d}m_2}=-\displaystyle\frac{\displaystyle\frac{4}{3} \left(\gamma _2^2-1 \right)+ f_1\displaystyle \displaystyle\frac{{d}m_3}{{d}m_2}+ (1-\epsilon _3)f_2\gamma _2\gamma _{42} \left(1-\displaystyle \displaystyle\frac{\beta _{42}}{\beta _4} \right)m_3\displaystyle \displaystyle\frac{{d}\psi _4}{{d}m_2}}{\displaystyle\frac{8}{3}(1-\epsilon _2)\gamma _2m_2+\epsilon _2m_2+(1-\epsilon _3)f_3m_3+\epsilon _3m_3}}, \end{equation} \tag{ 4 }$$

where f₁, f₂, and f₃ are functions of other variables (see Appendix A). The evolution equations for γ₄ and ψ₄ (see the definition at Appendix A) can also be obtained together.

2.3. After the Reverse Shock

After the reverse shock crosses region 4, only the forward shock is left in the system. Region 3 becomes the tail of region 2. We can get the hydrodynamical equation of γ₂ in the way similar to the section above (see Appendix B):

$$\begin{equation} \displaystyle {\frac{{d}\gamma _2}{{d}m_2}=-\frac{\displaystyle \frac{4}{3}\left(\gamma _2^2-1 \right)+(1-\epsilon _3)\displaystyle {\left[\frac{4}{3}+\frac{1}{3(1+\psi _3)^2}\right]}\gamma _2m_3\displaystyle \frac{{d}\psi _3}{{d}m_2}}{\displaystyle \frac{8}{3}(1-\epsilon _2)\gamma _2m_2+\epsilon _2m_2+(1-\epsilon _3)\displaystyle {\left[1+\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3\right]}m_3+\epsilon _3m_3}}. \end{equation} \tag{ 5 }$$

We briefly describe the radiation from the shocked electrons based on the standard model in this section. Unless noted otherwise, the quantities used in the equations in this section refer to the quantities in the shock comoving frame.

In the absence of radiation losses, the energy distribution of shock accelerated electrons behind the shock is usually assumed to be a power-law as $dN_e/d\gamma _e \propto \gamma _e^{-p}$, where γ_e is the Lorentz factor of the electron and p is the spectrum index. Since the electrons are cooled by synchrotron and inverse Compton (IC) radiation (Rybicki & Lightman 1979; Sari et al. 1998), the actual electron distribution should be given according to the following cases (Dai et al. 1999):

1. for $\gamma _{e,\rm c}\le \gamma _{e,\rm m}$,

$$\begin{equation} \frac{dN_{e}}{d\gamma _{e}}\propto \left\lbrace \begin{array}{@{}ll}\gamma _{e}^{-2},&\gamma _{e,\rm c}\le \gamma _{e}\le \gamma _{e,\rm m}\\ \gamma _{e}^{-p-1},&\gamma _{e,\rm m}<\gamma _{e}\le \gamma _{e,\rm max}\end{array}\right. \end{equation} \tag{ 6 }$$

2. for $\gamma _{e,\rm m}<\gamma _{e,\rm c}\le \gamma _{e,\rm max}$,

$$\begin{equation} \frac{dN_{e}}{d\gamma _{e}}\propto \left\lbrace \begin{array}{@{}ll}\gamma _{e}^{-p},&\gamma _{e,\rm m}\le \gamma _{e}\le \gamma _{e,\rm c}\\ \gamma _{e}^{-p-1},&\gamma _{e,\rm c}<\gamma _{e}\le \gamma _{e,\rm max}\end{array}\right. \end{equation} \tag{ 7 }$$

which are normalized by the total number of the electrons solved from the dynamic equations. The minimum, cooling, and maximum Lorentz factors of electrons are, respectively,

$$\begin{eqnarray} \gamma _{e,\rm m}&=&\epsilon _{e}\frac{p-2}{p-1}\frac{m_{\rm p}}{m_{\rm e}}(\gamma -1),\nonumber \\ \gamma _{e,\rm c}&=&\frac{6\pi m_{\rm e}c}{(1+Y)\sigma _{T}{B}^{2} (\gamma +\sqrt{\gamma ^{2}-1})(t-t_0)},\nonumber\\ \gamma _{e,\rm max}&=&\sqrt{6\pi q_{\rm e}\over \sigma _{T}B(1+Y)}, \end{eqnarray} \tag{ 8 }$$

where B is the comoving magnetic field strength, q_e is the electric charge of electrons, m_p and m_e are the mass of proton and electron, respectively, σ_T is the Thomson cross section, t₀ is the starting time of each episode in the observer frame, and Y is the Compton parameter that is defined as the ratio of the IC power to the synchrotron power. The Compton parameter of an electron with random Lorentz factor γ_e is given by $Y(\gamma _e)=(-1+\sqrt{1+4\eta _{\rm rad}\eta _{\rm KN}\epsilon _{e}/\epsilon _{B}})/2$ (He et al. 2009; Fan & Piran 2006; Wang et al. 2010), where η_rad is the fraction of the electron's energy that is radiated, η_KN is the fraction of synchrotron photons below the Klein–Nishina limit, and _e and _B are shock energy equipartition parameters for electrons and magnetic fields, respectively.

With the electron distribution determined, the synchrotron radiation flux density can be calculated following previous work (Rybicki & Lightman 1979; Huang et al. 2000a). The effect of equal-arrival-time surface (e.g., Waxman 1997; Granot et al. 1999) is considered in the calculations.

The equations above can describe the overall evolution of a blast wave encountering a sudden change in density. In order to explore the difference between our results and previous works, we first show the numerical results of dynamics. We calculate two different cases with density-jump magnitude (denoted as a) of 10 and 100 times, respectively. We further assume that the blast wave meets the density wall at γ₂ = 10 or γ₂ = 20. Below, we denote the four different calculations as γ10a10, γ10a100, γ20a10, and γ20a100 (like Gat et al. 2013). In our calculations, the typical values (e.g., Huang et al. 2000b; Freedman & Waxman 2001; Wu et al. 2003) adopted for parameters of the outflow are E_{K, iso} = 1.0 × 10⁵³ erg, θ_j = 0.1 rad, p = 2.3, _e = 0.1, _B = 0.01, Γ₀ = 300, and n₀ = 1.0 cm⁻³, where E_{K, iso} is the initial isotropic kinetic energy, and Γ₀ and θ_j are the initial Lorentz factor and half-opening angle of the jet, respectively.

Figure 1 shows the evolution of the Lorentz factors for different components during the entire time. Each panel in Figure 1 corresponds to one calculation. Before the encounter, the Lorentz factor of the forward shock satisfies the BM solution (γ∝t^−3/8) quite well. During the reverse shock crossing time, there is a short "plateau" for the evolution of γ₂. This is expected because regions 2 and 3 are being "pushed" by the faster region 4 behind them, so γ₂ will be a constant when an equilibrium between the deceleration by the medium and the acceleration from region 4 (i.e., dγ₂ ≃ 0) is established. After the reverse shock crossing, the bulk motion of regions 2 and 3 transition to a trajectory shallower than the BM solution, since they are already in the transitional period between the ultra-relativistic and the non-relativistic phases. However, in the non-relativistic phase at fairly late time, the evolution of γ₂ obeys the Sedov–Von Neumann–Taylor solution β∝R^−3/2 (Sedov 1959). Comparing the left panels with the right panels in Figure 1, it is within the expectation that the electrons in region 3 are hotter for an earlier density jump, which is more favorable to the re-brightening. The comparison between the upper panels and the lower panels shows that the larger the jump amplitude is, the stronger the reverse shock will be. However, the response in afterglow light curves depends weakly on the jump amplitude (see below).

Zoom In Zoom Out Reset image size

It is also convenient to investigate the difference of the hydrodynamics between previous results and ours. As shown in the upper left panel, the Lorentz factor ratio at the encounter time γ₄/γ₂(t = t_enc) is ∼2.2 for our result, whereas the values calculated from other analytical equations are ∼3.2 (Dai & Lu 2002) and ∼1.7 (NG: Nakar & Granot 2007), respectively. In contrast, the lower γ₂ and higher γ_{th, 3} predicted by Dai & Lu (2002) would lead to an overestimate of the flux emitted from hot electrons during this episode. The small discrepancy between NG and ours is caused by the different approximations used when dealing with the conservative equations and it does not qualitatively change our conclusions below.

Before calculating the radiation, we note the number density of hot electrons in region 3. It is crucial to take the radial structure of blast wave into account according to the discussion in van Eerten et al. (2009). In fact, van Eerten and his collaborators emphasized that the difference in the density of region 3 between the simulation and the analytical results is another significant factor leading to the contradiction of their conclusions. So it is prudent to dealing with this density in our work. Although we made the thin shell assumption when deriving the hydrodynamic equations, it is easy for us to calculate the "real" (volume-averaged) density of region 3 here in our model. Using the velocities of the forward shock and reverse shock, we can obtain the real co-moving width of region 3 and the corresponding density (see Equation (A24)), which will replace the density derived from the rough jump conditions. Likewise, the real density in our model is much lower than that predicted by jump conditions. Figure 2 depicts the comoving density profile during the encounter at the observer time t = 40, 900 s (z = 1.0) for the upper left panel in Figure 1; the real density of region 3 is ∼3 times lower than analytical one. Note that our density does not satisfy the shock jump condition, which is due to the simplified homogeneous assumption we adopted. In reality, the hydrodynamic quantities should satisfy the shock jump condition at the shock and have gradients in the shocked region. This simplified assumption will not lead to severe error.

Zoom In Zoom Out Reset image size

Afterglow light curves can be obtained by considering the radiation process described in Section 3 (assuming a redshift z = 1). Figure 3 shows the light curves in the X-ray and optical band for the four calculations in Figure 1. The emission from region 3 is found to be less than or comparable to the emission from the high-latitude early shell, no notable bumps emerge after the density-jump in these four cases. In each panel, the total light curve smoothly transits to a steep decay phase after the encounter. In fact, the total emission after the encounter is mainly determined by the curvature effect (Fenimore et al. 1996)—after the reverse shock crossing time no more fresh, hot electrons are injected into region 3. According to the scaling relation of the curvature effect F_ν∝t^{−(2 + β)} (β is the spectrum index), the temporal indices are −2.65 (slow cooling) or −3.15 (fast cooling) for p = 2.3. The slopes of the light curves in Figure 3 are just within this range. The lower component, emission from region 2, first rises steeply when γ₂ stays at the plateau phase, then decays sharply due to the decrease of γ₂ and the comoving density of region 2. Finally, it turns to a normal decay close to the BM solution. The character of the steep decay after the encounter in our model is not seen in the simulation results of Gat et al. (2013). This is the manifestation of the deviation caused by the homogenous thin shell assumption (i.e., the sharp edge of each region) in our work. With the radial structure of the blast wave considered in the full hydrodynamics simulation, the result of a shallower and smoother light curve is likely expected.

Zoom In Zoom Out Reset image size

If the Lorentz factor just before the jump was much larger (e.g., γ₂ > 20), a plateau would emerge in our model. However, this condition is not often satisfied because we were concerned with the late afterglow (t > 10⁴ s). The contribution of reverse shock component is mild and partially depends on the jump amplitude from the comparison between the upper and lower panels.

We perform the same calculations for the cases in which the blast wave travels in a stellar wind environment before the encounter and enters into a homogeneous ISM after the encounter. The visible response to the density-jump as shown in Figure 4 is smaller compared with the entire ISM case. In the wind environment, the prior temporal index of light curve in the left panel is (1 − 3p)/4 ∼ −1.5 (Zhang & Mészáros 2004) and is steeper than that of the ISM cases. Thus the smoother transition from the prior light curve to the phase after the encounter seems reasonable. Light curves for the blast wave encountering a density change of various magnitudes at a late time (γ₄ = 3) are also presented in Figure 5. There are no observable bumps for various density jumps in this case, too, as shown in Figure 5. For the blast waves traveling from a stellar wind to the ISM (right panel), the shallow behavior of optical afterglows is similar to those in Gat et al. (2013). According to the context above, the lower the encounter Lorentz factor is, the weaker the emission flux of region 3. Thus the flux decay during the encounter is naturally explained by the domination of the tail emission of region 4 in our model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A set of extended differential equations for relativistic outflows encountering a density-jump medium are derived in this paper. Using these, we obtained a more accurate dynamic (Figure 1) and calculated an appropriate density for region 3 (Figure 2) in this issue. The extended equations given here can also be applied to study the collision of two homogeneous shells after minor modifications. According to our numerical results, no high-amplitude rebrightening would emerge when the hot electrons in region 3 begin to contribute to the emission flux after encountering the density wall. van Eerten et al. (2009) and Gat et al. (2013) used adaptive mesh relativistic hydrodynamic codes to conclude that sudden transitions in circumburst density are very unlikely to be the cause of the bumps. Taking a homogeneous slab and ignoring the radial structure of the blast wave in the previous analytical result may lead to an overestimation. Although there are still some discrepancies between the simulation results and the results from our semi-analytic model—probably due to the thin shell assumption imprinted in this work—our results disfavor the density-jump origin for bumps on afterglows from another perspective.

Generally speaking, the density-jump scenario cannot explain the observed late rebrightenings in many GRBs. Other scenarios involving late refreshed shocks or late activities of the central engine may explain these rebrightenings (Margutti et al. 2011). The late internal shock model can naturally explain the X-ray flares with sharp profiles (Falcone et al. 2007). Recent works also show that the central engine (a black hole) may be re-activated after the initial burst. The fallback of material onto the central black hole after the collapse could last for a long time (Woosley 1993; MacFadyen et al. 2001; Perna et al. 2014) and lead to late central engine activities (Perna et al. 2006; Kumar et al. 2008a, 2008b). The observational evidence for this process is highlighted in Wu et al. (2013) and Yu et al. (2013). The two-component jet model scenario may also be consistent with the gradual rebrightening of some GRB afterglows (Huang et al. 2004; Peng et al. 2005; Filgas et al. 2011). Thus, it still remains to be answered which of these scenarios (refreshed shocks, late activities of central engines, two-component jets) play the key role in the observation sample. The spectral evolution during the bump may help test the models for individual GRBs (Nardini et al. 2011b). Radio observations of afterglows are also valuable for theoretical modeling (Moin et al. 2013). It is interesting to see that some bumps may be related to the central engine activities, which makes it helpful for studying the properties of progenitor stars.

We thank the anonymous referee for helpful comments and suggestions that led to an overall improvement of this study. This work was supported by the National Basic Research Program of China (973 Program, grant No. 2014CB845800, and 2013CB834900) and the National Natural Science Foundation of China (grant No. 11033002 and 11322328). X.F.W. acknowledges support by the One-Hundred-Talents Program, the Youth Innovation Promotion Association, and the Strategic Priority Research Program The Emergence of Cosmological Structures of the Chinese Academy of Sciences (grant No. XDB09000000). L.L. acknowledges support by the Swedish National Space Board, and the Erasmus Mundus Joint Doctorate Program by grant Number 2011-1640 from the EACEA of the European Commission.

We derive the differential equations for the dynamics of the system when the reverse shock exists. The total energy of regions 2, 3, and 4 are

$$\begin{equation} E_{\rm t} = E_2 + E_3 + E_4, \end{equation} \tag{ A1 }$$

$$\begin{equation} E_i = (\gamma _i-1) m_i c^2 + (1-\epsilon _i) \gamma _i (U_i^{\prime } + P_i^{\prime }), \quad i=2,3,4, \end{equation} \tag{ A2 }$$

where

$$\begin{eqnarray} U_2^{\prime }&=&\psi _2 m_2 c^2=(\gamma _2-1)m_2 c^2, \nonumber \\ U_3^{\prime }&=&\psi _3 m_3 c^2=(\gamma _{34}-1) m_3 c^2+\gamma _{34}(U_3^{\prime }+P_3^{\prime })\nonumber\\ &=&\left[\gamma _{34}\frac{4\psi _4+5}{3(1+\psi _4)}\psi _4+(\gamma _{34}-1)\right]m_3 c^2, \nonumber\\ U_4^{\prime }&=&\psi _4 m_4 c^2, \nonumber \\ P_{i}^{\prime }&=&\frac{\psi _{\rm {i}}+2}{3(1+\psi _{\rm {i}})}\psi _{\rm {i}}m_{\rm {i}}c^2. \end{eqnarray} \tag{ A3 }$$

The parameter ψ₄ is time dependent and the expression for it will be obtained below. At the time when the reverse shock emerges (R = R₀), we can get the initial value of ψ₄: $\psi _{4,0} = ({U_{4,0}^{\prime }}/{m_{4,0} c^2})=\gamma _{4,0}-1$, where γ_{4, 0} equals the value of γ₂ at R₀. Meanwhile, the total mass of regions 3 and 4 is m_{4, 0} = m₃ + m₄. Substituting Equations (A2) and (A3) into Equation (A1), the total energy of the system is

$$\begin{eqnarray} E_{\rm t}&=&\frac{4}{3}(\gamma _2^2-1)m_2 c^2-\epsilon _2\frac{4\gamma _2+1}{3}(\gamma _2-1)m_2 c^2 \nonumber \\ &&+(\gamma _3-1)m_3c^2+(1-\epsilon _3)\gamma _3\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3 m_3c^2 \nonumber \\ &&+(\gamma _4-1)m_4 c^2+\gamma _4\frac{4\psi _4+5}{3(1+\psi _4)}\psi _4 m_4 c^2, \,\,\, \end{eqnarray} \tag{ A4 }$$

where we have assumed that γ₂ = γ₃ and ₄ = 0. The unshocked portion of region 4 will not lose its energy, although a conversion of the thermal energy to the bulk kinetic energy may happen. Thus we have

$$\begin{eqnarray*} \gamma _4\left[1+\frac{4\psi _4+5}{3(1+\psi _4)}\psi _4\right] &=& \gamma _{4,0}\left[1+\frac{4\psi _{4,0}+5}{3(1+\psi _{4,0})}\psi _{4,0}\right] \nonumber\\ &=& \frac{4\gamma _{4,0}^2-1}{3}, {\rm or} \end{eqnarray*}$$

$$\begin{equation} {d}\gamma _4=-\frac{4(1+\psi _4)^2+1}{(2\psi _4+1)(2\psi _4+3)}\frac{\gamma _4}{1+\psi _4}{d}\psi _4. \end{equation} \tag{ A5 }$$

Let us have a look at the evolution of ψ (represents ψ₄) due to adiabatic expansion. We discuss it in the co-moving frame. For a system with the mass m, pressure p', and volume V', the thermal Lorentz factor is γ_th = 1 + ψ. We take the adiabatic index as $\hat{\gamma }_{{\rm th}}\simeq ({4\gamma _{{\rm th}}+1}/{3\gamma _{{\rm th}}})$, so the pressure is

$$\begin{eqnarray} p^{\prime }&=&(\hat{\gamma }_{{\rm th}}-1)e^{\prime }=(\hat{\gamma }_{{\rm th}}-1)(\gamma _{{\rm th}}-1)\rho ^{\prime }c^2 \nonumber \\ &=&\displaystyle \frac{\gamma _{{\rm th}}^2-1}{3\gamma _{{\rm th}}}\rho ^{\prime }c^2=\frac{\psi (2+\psi)}{3(1+\psi)}\rho ^{\prime }c^2. \end{eqnarray} \tag{ A6 }$$

The equation of adiabatic expansion is

$$\begin{equation} mc^2{d}\psi =-p^{\prime }{d}V^{\prime }. \end{equation} \tag{ A7 }$$

The above equation leads to the solution of

$$\begin{equation} \frac{(1+\psi)^2-1}{(1+\psi _0)^2-1}=(\frac{\rho ^{\prime }}{\rho _0^{\prime }})^{\frac{2}{3}}. \end{equation} \tag{ A8 }$$

So the evolution of ψ₄ can be written as

$$\begin{eqnarray} \frac{{d}\psi _4}{{d}m_2}&=&\frac{{d}\psi _4}{{d}\ln {\rho _4^{\prime }}} \frac{{d}\ln {\rho _4^{\prime }}}{{d}\ln {R}}\frac{{d}\ln {R}}{{d}m_2} \nonumber \\ &=&\frac{\psi _4(2+\psi _4)}{3(1+\psi _4)}\frac{{d}\ln {\rho _4^{\prime }}}{{d}\ln {R}} \frac{{d}R}{R{d}m_2}, \end{eqnarray} \tag{ A9 }$$

where

$$\begin{eqnarray} \frac{{d}\ln {\rho _4^{\prime }}}{{d}\ln {R}} &\,{=}\,&\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {R}}|_{\theta _4,\gamma _4} \,{+}\,\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\theta _4}}|_{\gamma _4,R}\frac{{d}\ln {\theta _4}}{{d}\ln {R}} \,{+}\,\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\gamma _4}}|_{R,\theta _4}\frac{{d}\ln {\gamma _4}}{{d}\ln {R}} \nonumber \\ &\,{=}\,&\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {R}}|_{\theta _4,\gamma _4} \,{-}\,\frac{\sin {\theta _4}}{1-\cos {\theta _4}}\frac{c_{s,4}}{\beta _4\gamma _4c} \,{+}\,\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\gamma _4}}|_{R,\theta _4}\frac{{d}\ln {\gamma _4}}{{d}\ln {R}}. \nonumber\\ \end{eqnarray} \tag{ A10 }$$

The second equality has included the lateral expansion of the shell (a detailed expression will be showed in Equation (A33)). Equations (A5), (A9), and (A10) lead to

$$\begin{equation} \frac{{d}\psi _4}{{d}m_2}=\frac{(1+\psi _4)\displaystyle \left[\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {R}}|_{\theta _4,\gamma _4} +\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\theta _4}}|_{\gamma _{4},R}\frac{{d}\ln {\theta _4}}{{d}\ln {R}}\right]\frac{{d}\ln {R}}{{d}m_2}}{\displaystyle {\frac{3(1+\psi _4)^2}{\psi _4(2+\psi _4)}+\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\gamma _4}}|_{R,\theta _4}}}. \end{equation} \tag{ A11 }$$

We can also get

$$\begin{eqnarray} \frac{{d}\gamma _4}{{d}R} &=&{-}\frac{\displaystyle {\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {R}}|_{\theta _4,\gamma _4} +\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\theta _4}}|_{\gamma _{4},R}\frac{{d}\ln {\theta _4}}{{d}\ln {R}}}}{\displaystyle {\frac{3(1+\psi _4)^2}{\psi _4(2+\psi _4)}\frac{(2\psi _4+1)(2\psi _4+3)}{4(1+\psi _4)^2+1}+\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\gamma _4}}|_{R,\theta _4}}} \frac{\gamma _4}{R} \nonumber \\ &=&\frac{\displaystyle {\frac{\sin {\theta _4}}{1-\cos {\theta _4}}\frac{c_{s,4}}{\beta _4\gamma _4c}-\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {R}}|_{\theta _4,\gamma _4}}}{\displaystyle {\frac{3(1+\psi _4)^2}{\psi _4(2+\psi _4)}\frac{(2\psi _4+1)(2\psi _4+3)}{4(1+\psi _4)^2+1}+\frac{\partial \ln {\rho _4^{\prime }}}{\partial \ln {\gamma _4}}|_{R,\theta _4}}} \frac{\gamma _4}{R}. \nonumber\\ \end{eqnarray} \tag{ A12 }$$

On the other hand, the radiative energy is

$$\begin{equation} {d}E_{\rm t}=-\frac{1}{3}\epsilon _2(4\gamma _2+1)(\gamma _2-1){d}m_2 c^2-\epsilon _3\gamma _3\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3{d}m_3 c^2. \end{equation} \tag{ A13 }$$

Combining Equations (A1)–(A4) and (A13), and set ₄ = 0, we get

$$\begin{eqnarray} && \frac{4}{3}(\gamma _2^2-1){d}m_2+\frac{8}{3}\gamma _2m_2{d}\gamma _2+\gamma _2\left[1+\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3\right]{d}m_3 \nonumber \\ && +m_3{d}\gamma _2 +(1-\epsilon _3)\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3m_3{d}\gamma _2+(1-\epsilon _3) \nonumber \\ && \times \left[\frac{4}{3}+\frac{1}{3(1+\psi _3)^2}\right]\gamma _2m_3{d}\psi _3 =\gamma _4\left[1+\frac{4\psi _4+5}{3(1+\psi _4)}\psi _4\right]{d}m_3 \nonumber\\ && +\epsilon _2\frac{8\gamma _2-3}{3} m_2{d}\gamma _2. \end{eqnarray} \tag{ A14 }$$

Taking γ₃₄ = (1 − β₃β₄)γ₃γ₄, γ₂ = γ₃, we have

$$\begin{eqnarray} {d}\gamma _{34} &=& \left(1-\frac{\beta _4}{\beta _2} \right)\gamma _4{d}\gamma _2 + \left(1-\frac{\beta _2}{\beta _4} \right)\gamma _2{d}\gamma _4 \nonumber\\ &=& \frac{\beta _{24}\gamma _{24}}{\beta _2\gamma _2}{d}\gamma _2+ \frac{\beta _{42}\gamma _{42}}{\beta _4\gamma _4}{d}\gamma _4. \end{eqnarray} \tag{ A15 }$$

For m₃, since

$$\begin{eqnarray} & &{d}m_2=2\pi (1-\cos {\theta _2}) R^2\rho _1 {d}R=2\pi (1-\cos {\theta _2}) R^2\rho _1\beta _2 c{d}t_b, \nonumber \\ & & {d}m_3=2\pi (1-\cos {\theta _3})R^2\rho _4^{\prime }\gamma _4(\beta _4-\beta _3)c{d}t_b, \end{eqnarray} \tag{ A16 }$$

we have

$$\begin{equation} {d}m_3=\displaystyle {\left(\frac{\beta _4}{\beta _2}-1\right)\frac{\rho _4^{\prime }}{\rho _1}\frac{1-\cos {\theta _3}}{1-\cos {\theta _2}} \gamma _4{d}m_2}, \end{equation} \tag{ A17 }$$

where ${d}t_b=\gamma _2(\gamma _2+\sqrt{\vphantom{A^A}\smash{{{\gamma _2^2-1}}}}){d}t$ is measured in the burst's frame while dt is the difference in arrival times (we neglect the (1 + z) term here). ρ₁ is the environment density at R > R₀, and

$$\begin{equation} \rho _4^{\prime }=\rho _{4,0}^{\prime }\frac{\gamma _4}{\gamma _{4,0}}\frac{R_0^3}{R^3}\frac{1-\cos {\theta _{4,0}}}{1-\cos {\theta _4}} \end{equation} \tag{ A18 }$$

is the co-moving density of region 4, and $\rho _{4,0}^{\prime }$ and θ_{4, 0} are the parameters at R₀. Inserting Equations (A3) and (A15) into (A14), we have the evolution of γ₂ during the reverse shock crossing

$$\begin{equation} \displaystyle {\displaystyle\frac{{d}\gamma _2}{{d}m_2}=-\displaystyle\frac{\displaystyle\frac{4}{3}\left(\gamma _2^2-1\right)+f_1\displaystyle \displaystyle\frac{{d}m_3}{{d}m_2}+(1-\epsilon _3)f_2\gamma _2\gamma _{42}\left(1-\displaystyle \displaystyle\frac{\beta _{42}}{\beta _4} \right)m_3 \displaystyle \displaystyle\frac{{d}\psi _4}{{d}m_2}}{\displaystyle\frac{8}{3}(1-\epsilon _2)\gamma _2m_2+\epsilon _2m_2+(1-\epsilon _3)f_3m_3+\epsilon _3m_3}}, \end{equation} \tag{ A19 }$$

where

$$\begin{equation} f_1=\gamma _2\left[1+\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3\right]-\gamma _4\left[1+\frac{4\psi _4+5}{3(1+\psi _4)}\psi _4\right], \end{equation} \tag{ A20 }$$

$$\begin{equation} f_2=\left[\frac{4}{3}+\frac{1}{3(1+\psi _3)^2}\right]\left[\frac{4}{3}+\frac{1}{3(1+\psi _4)^2}\right], \end{equation} \tag{ A21 }$$

$$\begin{eqnarray} f_3 &=& 1+\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3-\left[\!\frac{4}{3}+\frac{1}{3(1+\psi _3)^2} \!\right]\left[\! 1+\frac{4\psi _4+5}{3(1+\psi _4)}\psi _4 \!\right] \nonumber\\ && \times \frac{\beta _{42}}{\beta _2}\gamma _{42}. \end{eqnarray} \tag{ A22 }$$

Although we have made the thin shell assumption, we would calculate the co-moving width of region 3

$$\begin{equation} {d}\Delta _3^{\prime }=\frac{\gamma _3(\beta _4-\beta _3)c{d}t_b}{\frac{\gamma _3 \rho _3^{\prime }}{\gamma _4 \rho _4^{\prime }}-1} \end{equation} \tag{ A23 }$$

and its corresponding "real" (volume-averaged) density

$$\begin{equation} \rho _3^{\prime }=\displaystyle {\frac{m_3}{2\pi (1-\cos {\theta _3}) R^2 \Delta _3^{\prime }}}. \end{equation} \tag{ A24 }$$

van Eerten et al. (2009) pointed out that the radial structure of blast wave is relevant with the density-jump issue. The real density here is also a good approach to the actual case to some extent. In Equation (A23), we take $\rho _3^{\prime } = (4\gamma _{34}+3)\rho _4^{\prime }$ (jump condition) for the calculation of each small increment of $\Delta _3^{\prime }$.

We extend the jump condition to either cold (Sari & Piran 1995) or hot shell after the reverse shock emerges by

$$\begin{equation} (1+\psi _4)\left(\frac{\beta _4}{\beta _2}-1 \right)^2\frac{\rho _4^{\prime }}{\rho _1}\gamma _4^2=1, \end{equation} \tag{ A25 }$$

with the criterion for the formation of the reverse shock (i.e., β₄₂c > c_{s, 4}; sound speed of region 4 in the co-moving frame).

If the shell is hot (corresponding to the density-jump case), 1 + ψ₄ = γ₄, $\rho _4^{\prime }=4\gamma _4\rho _0$, we get

$$\begin{equation} \beta _{42}=\frac{\sqrt{\Re }}{2+\sqrt{\Re }}\beta _4, \;\;\;\; \gamma _{42}=\frac{2+\sqrt{\Re }}{\sqrt{4(1+\sqrt{\Re })\gamma _2^2+\Re }}\gamma _4, \end{equation} \tag{ A26 }$$

where ℜ ≡ (ρ₁/ρ₀), with ρ₀∝R^−k(k = 0 for ISM; k = 2 for wind) at R ⩽ R₀. In relativistic stage, γ₄ ≫ 1, $c_{s,4}=c/\sqrt{3}$, the criterion for the formation of reverse shock is $\Re >(\sqrt{3}+1)^2\simeq 7.5$. The reverse shock can be relativistic (γ₄₂ ⩾ 2) if $\sqrt{\Re }\ge \displaystyle ({6+4\sqrt{3}\beta _4}/{4\beta _4^2-3})$ or ℜ ⩾ 168 (γ₄ ≫ 2). The analysis here is consistent with the results of Dai & Lu (2002). In the non-relativistic stage, β₄ ≪ 1, $c_{s,4}=\displaystyle ({\sqrt{5}}/{3})\beta _4 c$, the criterion for the formation of reverse shock becomes $\Re >(35+15\sqrt{5})/2=34.3$.

If the shell is cold, ψ₄ = 0, we have

$$\begin{equation} \beta _{42}=\frac{\gamma _4}{\gamma _4+\sqrt{f}}\beta _4, \;\;\;\; \gamma _{42}=\frac{\gamma _4+\sqrt{f}}{\sqrt{1+f+2\gamma _4\sqrt{f}}}, \end{equation} \tag{ A27 }$$

where $f\equiv \displaystyle ({\rho _4^{\prime }}/{\rho _1})$. Since the initial shell is cold (c_{s, 4} = 0), the reverse shock can always be developed. The reverse shock can be relativistic (γ₄₂ ⩾ 2) if $\sqrt{f}\le \displaystyle ({2\sqrt{3}\beta _4-3}/{3})\gamma _4$, or $f\le \gamma _4^2/(3+2\sqrt{3})^2\simeq \gamma _4^2/42$.

However, the observed time of region 4 is shorter than that of regions 3 and 2, which is because of the difference in their bulk Lorentz factor (i.e., γ₄ > γ₃ = γ₂), while the radial increment is the same. The relation between the observed times is

$$\begin{equation} {d}t_4=\displaystyle \frac{\beta _2\gamma _2 \left(\gamma _2+\sqrt{\gamma _2^2-1}\right)}{\beta _4\gamma _4 \left(\gamma _4+\sqrt{\gamma _4^2-1} \right)}{d}t. \end{equation} \tag{ A28 }$$

We still need four more equations to complete the hydrodynamics, i.e.,

$$\begin{equation} \displaystyle \frac{{d}R}{{d}t}=\beta _2c\gamma _2 \left(\gamma _2+\sqrt{\gamma _2^2-1} \right), \end{equation} \tag{ A29 }$$

$$\begin{equation} \displaystyle \frac{{d}\theta _2}{{d}t}=\displaystyle \frac{c_{s,2} \left(\gamma _2+\sqrt{\gamma _2^2-1} \right)}{R}, \end{equation} \tag{ A30 }$$

$$\begin{equation} \displaystyle \frac{{d}\theta _3}{{d}t}=\displaystyle \frac{c_{s,3} \left(\gamma _2+\sqrt{\gamma _2^2-1}\right)}{R}, \end{equation} \tag{ A31 }$$

$$\begin{equation} \displaystyle \frac{{d}\theta _4}{{d}t}=\displaystyle \frac{c_{s,4}}{R} \displaystyle \frac{\beta _2\gamma _2 \left(\gamma _2+\sqrt{\gamma _2^2-1} \right)}{\beta _4\gamma _4}, \end{equation} \tag{ A32 }$$

or

$$\begin{equation} \frac{{d}\theta _{i}}{{d}R}=\frac{c_{s,i}}{\beta _{i}\gamma _{i}c}\frac{1}{R}, \end{equation} \tag{ A33 }$$

where the sound speeds in the co-moving frame are

$$\begin{equation} c_{s,i}^2=\hat{\gamma }_{{\rm th},i}(\hat{\gamma }_{{\rm th},i}-1)(\gamma _{{\rm th},i}-1) \displaystyle \frac{1}{1+\hat{\gamma }_{{\rm th},i}(\gamma _{{\rm th},i}-1)}c^2, \end{equation} \tag{ A34 }$$

in which γ_{th, 2} = γ₂ for region 2, γ_{th, 3} = γ₄₂(1 + ψ₄) for region 3, and γ_{th, 4} = 1 + ψ₄ for region 4.

We derive the differential equations for the dynamics of the system after the reverse shock crossing time t_across. Region 4 vanishes and the total energy of the system is now

$$\begin{equation} E_{\rm t}=E_{\rm {2}}+E_{\rm {3}}, \end{equation} \tag{ B1 }$$

with

$$\begin{equation} E_{\rm {2}}=\frac{4}{3}\left(\gamma _2^2-1 \right)m_2 c^2-\epsilon _2\frac{4\gamma _2+1}{3}(\gamma _2-1)m_2c^2, \end{equation} \tag{ B2 }$$

$$\begin{equation} E_{\rm {3}}=(\gamma _2-1)m_3 c^2+(1-\epsilon _3)\gamma _2\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3m_3c^2, \end{equation} \tag{ B3 }$$

where m₃ = m_{4, 0} = const, $\psi _{3,0}=\gamma _{34} [1+\displaystyle ({4\psi _4+5}/{3(1+\psi _4)}) \psi _4]-1|_{t=t_{{\rm across}}}$. The radiative energy is

$$\begin{equation} {d}E_{\rm t}=-\frac{1}{3}\epsilon _2(4\gamma _2+1)(\gamma _2-1){d}m_2 c^2. \end{equation} \tag{ B4 }$$

The combination of the above four equations leads to

$$\begin{eqnarray} && \frac{4}{3}\left(\gamma _2^2-1 \right){d}m_2+\frac{8}{3}\gamma _2m_2{d}\gamma _2 +m_3{d}\gamma _2 + (1-\epsilon _3)\frac{4\psi _3+5}{3(1+\psi _3)} \nonumber \\ && \quad\quad \times \psi _3m_3{d}\gamma _2+(1-\epsilon _3)\left[\frac{4}{3}+\frac{1}{3(1+\psi _3)^2}\right]\gamma _2m_3{d}\psi _3 \nonumber \\ && =\frac{1}{3}\epsilon _2(8\gamma _2-3)m_2{d}\gamma _2, \end{eqnarray} \tag{ B5 }$$

thus the hydrodynamic equation is

$$\begin{equation} \displaystyle {\displaystyle\frac{{d}\gamma _2}{{d}m_2}=-\displaystyle\frac{\displaystyle \displaystyle\frac{4}{3}\left(\gamma _2^2-1 \right)+(1-\epsilon _3)\displaystyle {\left[\displaystyle\frac{4}{3}+\displaystyle\frac{1}{3(1+\psi _3)^2}\right]}\gamma _2m_3\displaystyle \displaystyle\frac{{d}\psi _3}{{d}m_2}}{\displaystyle \displaystyle\frac{8}{3}(1-\epsilon _2)\gamma _2m_2+\epsilon _2m_2+(1-\epsilon _3)\displaystyle {\left[1+\displaystyle\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3\right]}m_3+\epsilon _3m_3}}. \end{equation} \tag{ B6 }$$

For simplicity, we set ₃ = 0, and the final result is

$$\begin{equation} \displaystyle {\frac{{d}\gamma _2}{{d}m_2}=-\frac{\displaystyle \frac{4}{3}(\gamma _2^2-1)+\displaystyle {\left[\frac{4}{3}+\frac{1}{3(1+\psi _3)^2}\right]}\gamma _2m_3\displaystyle \frac{{d}\psi _3}{{d}m_2}}{\displaystyle \frac{8}{3}(1-\epsilon _2)\gamma _2m_2+\epsilon _2m_2+\displaystyle {\left[1+\frac{4\psi _3+5}{3(1+\psi _3)}\psi _3\right]}m_3}}. \end{equation} \tag{ B7 }$$

This result is consistent with the generic model of Huang et al. (1999) if ψ₃ = 0.

The evolution of ψ₃ is model dependent. Here we give two scenarios: the first uses the work done by region 3 to region 2, and the second is considered by the adiabatic expansion of region 3. In both scenarios the decrease of the adiabatic thermal energy of region 3 can postpone the deceleration of the bulk motion of both region 3 and region 2. Here, we only discuss the adiabatic expansion scenario.

The evolution of ψ₃ also follows the adiabatic expansion of Equation (A9), and

$$\begin{equation} \rho _3^{\prime }=\rho _{3,a}^{\prime }\frac{\gamma _2}{\gamma _{2,a}}\frac{R_a^3}{R^3}\frac{1-\cos {\theta _{3,a}}}{1-\cos {\theta _3}}, \end{equation} \tag{ B8 }$$

in which $\rho _{3,a}^{\prime }\simeq 4\gamma _{42}\rho _4^{\prime }|_{t_{{\rm across}}}$, γ_{2, a}, and R_a are the parameters when the reverse shock just crosses region 4. So

$$\begin{eqnarray} \!\frac{{d}\ln {\rho _3^{\prime }}}{{d}\ln {R}} &\,{=}\,&\frac{\partial \ln {\rho _3^{\prime }}}{\partial \ln {R}}|_{\theta _3,\gamma _2} \,{+}\,\frac{\partial \ln {\rho _3^{\prime }}}{\partial \ln {\theta _3}}|_{\gamma _2,R}\frac{{d}\ln {\theta _3}}{{d}\ln {R}} \,{+}\,\frac{\partial \ln {\rho _3^{\prime }}}{\partial \ln {\gamma _2}}|_{R,\theta _3}\frac{{d}\ln {\gamma _2}}{{d}\ln {R}} \nonumber \\ &=&{-}3-\frac{\sin {\theta _3}}{1-\cos {\theta _3}}\frac{c_{s,3}}{\beta _2\gamma _2c}+\frac{{d}\ln {\gamma _2}}{{d}\ln {R}}. \end{eqnarray} \tag{ B9 }$$

We thus have

$$\begin{equation} \frac{{d}\psi _3}{{d}m_2} =\frac{\psi _3(2+\psi _3)}{3(1+\psi _3)}\frac{{d}\ln {\rho _3^{\prime }}}{{d}\ln {R}} \frac{{d}R}{R{d}m_2}, \end{equation} \tag{ B10 }$$

and

$$\begin{eqnarray} &&\frac{{d}\gamma _2}{{d}m_2}=\nonumber\\ &&-\frac{\displaystyle \frac{4}{3}\left(\gamma _2^2-1\right)-\gamma _2m_3\displaystyle \frac{\psi _3(2+\psi _3)}{3(1+\psi _3)} \left[\frac{4}{3}+\frac{1}{3(1+\psi _3)^2}\right]\left[3+\frac{\sin {\theta _3}}{1-\cos {\theta _3}}\frac{c_{s,3}}{\beta _2\gamma _2c }\right]\displaystyle \frac{{d}R}{R{d}m_2}}{\displaystyle \frac{8}{3}(1-\epsilon _2)\gamma _2m_2+\epsilon _2m_2+\displaystyle \frac{ \left(2\psi _3^2+4\psi _3+1\right)\left(8\psi _3^2+16\psi _3+9\right)}{9(1+\psi _3)^3}m_3}. \nonumber\\ \end{eqnarray} \tag{ B11 }$$