NO FLARES FROM GAMMA-RAY BURST AFTERGLOW BLAST WAVES ENCOUNTERING SUDDEN CIRCUMBURST DENSITY CHANGE

The initial conditions for an adiabatic blast wave ("jet") formed by a GRB are described by the BM solution in spherical coordinates. For this paper, we use a conical section truncated at a certain opening angle for the initial setup instead of the full spherical solution, for at early times, spreading has not yet occurred. Using the initial conditions for a jet flowing in the radial direction with an opening angle, θ₀, and a jet energy, E_jet, we obtain the isotropic equivalent energy, E_iso:

$$\begin{equation} E_{{\rm jet}} = E_{{\rm iso}}(1 - \cos \theta _0)\approx E_{{\rm iso}}\theta _0^2/2. \end{equation} \tag{ 1 }$$

The initial radius, time, and dynamics of the interaction of the jet with the circumburst medium are calculated using the BM self-similar solution. For this paper, we represent the circumburst medium interior to the density change using a power-law density profile which is described as:

$$\begin{equation} \rho _{{\rm ext}} = \rho _{{\rm ref}}\left(\frac{r}{R_{{\rm ref}}}\right)^{-k}, \end{equation} \tag{ 2 }$$

where ρ_ref is the circumburst density at a reference radius, R_ref, and k = 2  or  0 to represent the stellar wind or ISM environment respectively. An ISM is a region of constant density, and throughout this paper we use "ISM" and "homogeneous medium" interchangeably. The isotropic energy of the system is given by (BM):

$$\begin{equation} E_{{\rm iso}} = \frac{8\pi \rho _{{\rm ext}}\gamma ^2 t^3 c^5}{17 - 4k}, \end{equation} \tag{ 3 }$$

where γ is the fluid Lorentz factor just behind the shock, and $\Gamma = \sqrt{2}\gamma$ is the Lorentz factor of the shock itself. Substituting Equation (2) into Equation (3), and using the approximation that at ultra-relativistic speeds, R₀ ≈ cT₀, the initial radius of the shock is:

$$\begin{equation} R_0 = \left(\frac{(17-4k)E_{{\rm iso}}}{8\pi \rho _0 R_{{\rm ref}}^k c^2 \gamma _0^2}\right)^{\frac{1}{3-k}}. \end{equation} \tag{ 4 }$$

Using these expressions, we can solve for the initial time of the simulation:

$$\begin{equation} T_0 = \frac{R_0}{c}\left(1 + \frac{1}{2(4-k)\Gamma _0^2}\right). \end{equation} \tag{ 5 }$$

As the jet propagates into the circumburst medium described by Equation (2), the Lorentz factor and the radius of the shock over time are well known—the radius will follow a similar solution to Equation (5). The subscript of "FS" is used to denote the front of the shock, or the forward shock:

$$\begin{equation} R_{{\rm FS}} = cT\left(1- \frac{1}{2(4-k)\Gamma ^2}\right). \end{equation} \tag{ 6 }$$

The Lorentz factor evolves as

$$\begin{equation} \Gamma = \Gamma _0\left(\frac{T}{T_0}\right)^{-\frac{3-k}{2}}. \end{equation} \tag{ 7 }$$

The initial conditions at T₀ of our numerical simulations use the full BM profile, but for simplicity in our analytical solution, we assume the blast wave is a homogeneous shell of constant density within the shell, and the back of the shell being at a radius of R_back:

$$\begin{equation} R_{{\rm back}} = R_{{\rm FS}}\left(1 - \frac{1}{2(3-k)\Gamma _{{\rm FS}}^2}\right). \end{equation} \tag{ 8 }$$

Equation (8) is derived to the leading order of 1/γ² using the approximation that the total swept up mass is contained within the shell.

To simulate the encounter with a change in external medium, we model the external density profile as a piecewise function:

$$\begin{equation} \rho _{{\rm ext}} = \left\lbrace \begin{array}{@{}rl}\rho _{{\rm wind}}&\mbox{ $\gamma >\gamma _{{\rm enc}}$} \\ \rho _{{\rm ISM}}&\mbox{ $\gamma \le \gamma _{{\rm enc}}$}\end{array},\right. \end{equation} \tag{ 9 }$$

where γ_enc is the fluid Lorentz factor at which the encounter with the new environment occurs. The dynamics before the encounter are described by the equations above, but the dynamics during the encounter and after the encounter are very different. The simulation does not necessarily return to a BM self-similar solution for a relativistic blast wave in a homogeneous medium, for the BM solution is a slow attractor and a perturbation in the simulation, such as a density change, can cause large deviations from the BM solution (Gruzinov 2000). Also, after the encounter, the blast wave may no longer be relativistic, which prevents the blast wave from evolving toward the BM solution in a homogeneous medium, and a solution similar to ST is expected instead.

To analytically model the dynamics of a blast wave in a stellar wind environment encountering a density change followed by a homogeneous medium, we extend the derivations performed in PW and NG. With our extensions, the analytical solution is applicable for a wide range of density jumps as well as density drops.

For simplicity, we model the shocks of the blast wave as simple homogeneous shells. During the encounter, the shock breaks up into three different regions—regions 2, 3, and 4 of Figure 1. To get the full analytical solution of the simulation, the pressures, densities, Lorentz factors, and radii of the three regions (regions 2, 3, and 4, which are henceforth referenced by subscripts in this paper) are needed.

Zoom In Zoom Out Reset image size

Across the contact discontinuity, the pressure and fluid Lorentz factors are equal, meaning p₂ = p₃, and γ₂ = γ₃ = Γ_CD. Before the reverse shock passes the "back" of the shock, the fluid in region 4 has no knowledge of the encounter, meaning the density, pressure, and Lorentz factor are calculated as if no encounter has occurred. The rest of the fluid has experienced the encounter, and to calculate the values analytically, a convenient factor to define is that of NG Equation (3), ψ, which can be derived from conservation of mass:

$$\begin{equation} \psi ^2 \equiv \left[\frac{\gamma _{4}(T_{{\rm enc}})}{\gamma _2(T_{{\rm enc}})}\right]^2 = \frac{3a - 4}{\sqrt{12/a}(a-1) - 1} \end{equation} \tag{ 10 }$$

$$\begin{equation} a = \frac{\rho _{{\rm ISM}}}{\rho _{{\rm wind}}(R = R_{{\rm enc}})}. \end{equation} \tag{ 11 }$$

Using this factor, the Lorentz factor of the fluid behind the forward shock is described by a modified version of Equation (6) from NG:

$$\begin{equation} \gamma _2 = \gamma _{{\rm enc}}\left[\psi ^2 + \left(\frac{17 - 4k}{17} \right)a\left(\left(\frac{T}{T_{{\rm enc}}}\right)^3 - 1\right) \right]^{-1/2}, \end{equation} \tag{ 12 }$$

where Equation (12) differs from NG Equation (6) in that γ₂∝γ_enc as opposed to γ₂∝γ₄ of NG Equation (6).

To calculate the radius of the forward shock to the leading order of $1/\gamma _2^2$, we integrate the velocity of the forward shock:

$$\begin{eqnarray} R_{{\rm FS}}(T) &= &R_{{\rm enc}} + c\int ^T_{T_{{\rm enc}}}\beta _{{\rm FS}} dt \nonumber\\ &\cong & R_{{\rm enc}} + cT_{{\rm enc}} a\frac{17 - 4k}{68\gamma _{{\rm enc}}^2} \left(\frac{T}{T_{{\rm enc}}} - \frac{3}{4} - \frac{T^4}{4T_{{\rm enc}}^4}\right) \nonumber\\ &&+ c T_{{\rm enc}}\left(\frac{T}{T_{{\rm enc}}} - 1\right)\left(1 - \frac{\psi ^2}{4\gamma _{{\rm enc}}^2} \right). \end{eqnarray} \tag{ 13 }$$

The radius of the contact discontinuity is calculated in the same fashion as Equation (13), except with β_FS replaced by β_CD.

The pressure in regions 2 and 3 are identical, and are found using the strong shock jump conditions:

$$\begin{equation} p_2 = p_3 = \frac{2}{3}\rho _{{\rm ISM}}\Gamma _{{\rm FS}}^2c^2. \end{equation} \tag{ 14 }$$

One important aspect to note is that the reverse shock is not necessarily strong or relativistic in the frame of the fluid in region 4: if the forward shock is relativistic, the reverse shock is not (see PW for details). However, if the blast wave encounters a strong density jump causing the forward shock to become non-relativistic, the reverse shock can form with relativistic speeds. Since we cannot assume a relativistic reverse shock, we use the full shock jump conditions yielding the following equation for the density in region 3:

$$\begin{equation} \rho _3 = \frac{4p_3\bar{\gamma }_3 \rho _4}{3p_4+ p_2}, \end{equation} \tag{ 15 }$$

where $\bar{\gamma }_3$ is the Lorentz factor of the fluid in region 3 from the frame of the unshocked fluid in region 4:

$$\begin{equation} \bar{\gamma }_3 = \frac{1}{2}\left[ \frac{\gamma _4}{\gamma _3} + \frac{\gamma _3}{\gamma _4}\right]. \end{equation} \tag{ 16 }$$

Lastly, we need to calculate the velocity and radius of the reverse shock. We use conservation of mass to obtain the velocity of the reverse shock:

$$\begin{equation} \beta _{{\rm RS}} = \frac{\rho _3\gamma _3\beta _3 - \rho _4\gamma _4\beta _4}{\rho _3\gamma _3 -\rho _4\gamma _4}. \end{equation} \tag{ 17 }$$

To calculate the radius of the reverse shock, we use conservation of mass as well. The total mass in regions 3 and 4 is constant after the encounter, meaning that the integral over these regions equals the swept-up mass prior to the encounter:

$$\begin{equation} \int _{R_{{\rm back}}}^{R_{{\rm RS}}}\!r^2\rho ^{\prime }_4\, {d}r + \int _{R_{{\rm RS}}}^{R_{{\rm CD}}}\!r^2\rho ^{\prime }_3\, {d}r = \int _0^{R_{{\rm enc}}}\! r^2 \rho _{{\rm wind}}(r)\, {d}r. \end{equation} \tag{ 18 }$$

Equation (18) is important because we know the densities in regions 3 and 4, we have calculated R_back from Equation (8), we know R_CD from Equation (13) but with β_FS replaced by β_CD, and we know the total mass swept up prior to the encounter (the right hand side of Equation (18)). This leaves us with an equation for the radius of the reverse shock:

$$\begin{equation} R_{{\rm RS}}^3 (\rho ^{\prime }_4 - \rho ^{\prime }_3) = \frac{3}{3-k}R_{{\rm enc}}^3\rho _{{\rm wind}}(R_{{\rm enc}}) \,{+}\, R_{{\rm back}}^3\rho ^{\prime }_4 - R_{{\rm CD}}^3\rho ^{\prime }_3.\quad \end{equation} \tag{ 19 }$$

Equations (18) and (19) contain ρ', which is the density of the fluid in the lab frame, ρ' = γρ, whereas ρ, with no prime, is the co-moving density.

Solving for R_RS in Equation (19) is only applicable when the reverse shock has not yet passed the back of the shock. Once the reverse shock passes the back of the shock, Equation (18) no longer applies, and the new integral to solve is:

$$\begin{eqnarray} \int _{R_{{\rm back}}}^{R_{{\rm CD}}}\!r^2\rho ^{\prime }_3\, {d}r &=& \int _0^{R_{{\rm enc}}}\! r^2 \rho _{{\rm wind}}(r)\, {d}r \nonumber\\ R_{{\rm back}}^3 &=& R_{{\rm CD}}^3 - \frac{3}{3-k}R_{{\rm enc}}^3\frac{\rho _{{\rm wind}}(R_{{\rm enc}})}{\rho ^{\prime }_3}, \end{eqnarray} \tag{ 20 }$$

where we now solve for R_back instead of R_RS, and R_back is no longer calculated from Equation (8) because region 4 has been shocked by the reverse shock and has knowledge of the encounter. The blast wave continues on with these two regions of Equation (20)—one containing the newly swept-up mass in region 2, and one containing the mass swept up prior to the encounter—resulting in the blast wave never completely returning to the BM solution of a blast wave in an ISM.

In the following section, Section 2.3, we compare our analytical solution with the numerical results for various types of encounters.

Our numerical simulations demonstrate the resulting dynamics of an afterglow blast wave in a stellar wind environment encountering a sudden change in density and an ISM environment. We use the ram parallel RHD code (Zhang & MacFadyen 2006, 2009) with a second-order weighted scheme. We simulate two different size jumps ("walls"), of factors 4 and 100, and one drop, of factor $\frac{1}{100}$. We use these specific initial conditions to numerically simulate a wide range of encounters thought to potentially cause light curve flares. However, as we explain in Section 3, none of these scenarios result in flares on the time scale of 0.1–10 days. Before delving into the resulting light curves, we first explain the dynamics and numerical results.

We set up six simulations, all of which have a starting time and radius equivalent to the fluid Lorentz factor behind the front of the shock, γ = 15. This starting Lorentz factor was chosen to ensure $\gamma >\frac{1}{\theta _0}$, where θ₀ = 0.1 rad (typical for afterglow jets; Frail et al. 2001) for our simulations and that it is the half opening angle of the jet.

• 1.  
  Simulation no encounter:

  $$\begin{equation*} \rho _{{\rm ext}} = \rho _0\left(\frac{r}{R_{{\rm ref}}}\right)^{-2}. \end{equation*}$$
• 2.  
  Simulation γ5a4:

  $$\begin{equation*} \rho _{{\rm ext}} = \left\lbrace \begin{array}{@{}lr}\rho _0\left(\frac{r}{R_{{\rm ref}}}\right)^{-2} & \gamma >5 \\ 4\rho _0\left(\frac{R_{\gamma =5}}{R_{{\rm ref}}}\right)^{-2} &\gamma \le 5. \end{array}\right. \end{equation*}$$
• 3.  
  Simulation γ5a100:

  $$\begin{equation*} \rho _{{\rm ext}} = \left\lbrace \begin{array}{@{}lr}\rho _0\left(\frac{r}{R_{{\rm ref}}}\right)^{-2} &\gamma >5 \\ 100\rho _0\left(\frac{R_{\gamma =5}}{R_{{\rm ref}}}\right)^{-2} &\gamma \le 5. \end{array}\right. \end{equation*}$$
• 4.  
  Simulation γ10a4:

  $$\begin{equation*} \rho _{{\rm ext}} = \left\lbrace \begin{array}{@{}ll}\rho _0\left(\frac{r}{R_{{\rm ref}}}\right)^{-2} &\gamma >10 \\ 4\rho _0\left(\frac{R_{\gamma =10}}{R_{{\rm ref}}}\right)^{-2} &\gamma \le 10. \end{array}\right. \end{equation*}$$
• 5.  
  Simulation γ10a100:

  $$\begin{equation*} \rho _{{\rm ext}} = \left\lbrace \begin{array}{@{}ll}\rho _0\left(\frac{r}{R_{{\rm ref}}}\right)^{-2} &\gamma >10 \\ 100\rho _0\left(\frac{R_{\gamma =10}}{R_{{\rm ref}}}\right)^{-2} &\gamma \le 10. \end{array}\right. \end{equation*}$$
• 6.  
  Simulation $\gamma 5a\frac{1}{100}$:

  $$\begin{equation*} \rho _{{\rm ext}} = \left\lbrace \begin{array}{@{}ll}\rho _0\left(\frac{r}{R_{{\rm ref}}}\right)^{-2} &\gamma >5 \\ 0.01\rho _0\left(\frac{R_{\gamma =5}}{R_{{\rm ref}}}\right)^{-2} &\gamma \le 5. \end{array}\right. \end{equation*}$$

We first ran all six simulations in 1D, and then ran three simulations in 2D, choosing to run those which would be most informative to answering the questions outlined in the Introduction.

Figure 2 shows the analytical solution of the density of simulation γ5a4 in 1D plotted against the numerical result of simulation γ5a4 in 1D, showing the correspondence of our analytical model with our numerical results of the shock formation at the encounter. Our analytical model assumes homogeneous shells with constant density between the shocks, which is why the analytical solution in Figure 2 has a clearly marked "back" of the shock, whereas the numerical simulation uses the full BM profile and has no exact "back" of the shock. The small discrepancy between the analytical and numerical R_FS shown in Figure 2 is roughly ΔR_FS = 9.6 × 10¹⁵ cm. Our analytical solution is accurate to the order of 1/γ² = 0.0855 at the time shown in Figure 2 and (ΔR_FS/cT) = 3.893 × 10⁻⁴ ≪ 1/γ², meaning that the analytical model conforms to the numerical results to the required accuracy. Also, it is apparent from Figure 2 that the densities within each region are not constant, but actually have a positive slope. This is a known discrepancy between the assumed homogeneous shells and the full BM solution, and leads to a slight overprediction of the flux during the encounter by the analytical solution (van Eerten et al. 2009).

Zoom In Zoom Out Reset image size

Figures 3 and 4 show our analytical solutions plotted with the simulation results of the Lorentz factors behind the forward shocks as they evolve. The small deviations toward the beginning of the simulations (at smaller radii) in these figures are from lack of resolution. The resolution required to resolve the blast wave at these early times is very high. However, for the simulations with encounters at γ₂ = 5, the discrepancy disappears well before the encounter. For the simulations with the encounter at γ₂ = 10, the lack of resolution does affect the convergence of the numerical and analytical results, but qualitatively, the results are the same. The resolutions of our simulations are discussed further in Appendix A. The deviations at large radii (at very late times) in Figures 3 and 4 are from the Lorentz factor of the forward shock becoming non-relativistic, and our analytical model having accuracy of order 1/γ².

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the examples shown in Figures 2–4, it is apparent that our analytical solution accurately predicts the dynamics of a 1D blast wave in a stellar wind encountering a change in density, be it a jump or drop, followed by an ISM environment. Also, it is important to note that when the post-encounter forward shock remains relativistic, the fluid Lorentz factor at which the encounter occurs does not qualitatively change the post-encounter dynamics in 1D, as seen in Figure 3. The post-encounter dynamics are simply scaled to reflect the larger-encounter Lorentz factor. This is also clear from the analytical model, since neither the change in Lorentz factor at the encounter nor the pre- and post-encounter slopes depend on the Lorentz factor.

The 1D simulations show that the dynamics follow the analytical predictions that the Lorentz factor of the forward shock will drop immediately as it encounters a higher density region, and that the amount by which it drops is proportional to the change in density at the encounter. This is shown in Figure 5 for simulations γ5a4, γ5a(1/100), and γ10a100. Figure 5 depicts the Lorentz factor in the color coding, with the complete radial fluid profile along the central axis of the blast wave (plotted in the horizontal) over time (plotted in the vertical). This figure shows the evolution of the forward shock over time before and after the encounter as well as the formation of the reverse shock at the encounter and its post-encounter evolution. Figure 5 clearly illustrates the varying behaviors of the forward and reverse shocks for the various simulations. Simulations γ5a100 and γ10a4 show similar behavior. The left column in Figure 5 shows the Lorentz factor over time and radius for simulations γ5a4 and γ10a100 in 1D. If the blast wave encounters a drop in density, the Lorentz factor of the blast wave will increase (i.e., the blast wave will speed up), again in proportion to the size of the drop, which is seen in the top right panel in Figure 5.

Zoom In Zoom Out Reset image size

The effects of the magnitude of the density change on the dynamics of the reverse shock are shown in Figure 5. The top left panel in Figure 5 depicts the Lorentz factor for the simulation with a small density increase of factor 4. Shortly after the encounter, the reverse shock slowly travels away from the forward shock toward the origin, but this is not an immediate process like that of a larger density jump shown in the bottom left panel in Figure 5. For a large density jump, almost immediately after the encounter, the reverse shock has a large enough Lorentz factor, compared to the Lorentz factor of the forward shock, to break away from the forward moving fluid and travel backward toward the origin. This difference in reverse shock Lorentz factors of simulation γ5a4 and simulation γ10a100 is from the forward shock becoming immediately non-relativistic after the encounter with a large density jump (of factor 100) versus the much smaller decrease in the forward shock Lorentz factor with a smaller density jump (of factor 4). As explained in Section 2.2, the amount by which the forward shock Lorentz factor drops is proportional to the Lorentz factor with which the reverse shock forms. For simulation γ5a(1/100), which has a density decrease of factor 100, the reverse shock takes much longer to overcome the forward moving fluid and start traveling backward toward the origin. This is seen in the top right panel in Figure 5. The Lorentz factor of the forward shock increases, causing the Lorentz factor of the reverse shock shortly after the encounter to be very small in comparison with the forward shock Lorentz factor.

A closer look at the bottom left panel in Figure 5 reveals the reflecting-like behavior of the inner boundary as the reverse shock travels toward the origin and is seemingly reflected back. This behavior is independent of boundary conditions or the inner boundary position and results from the increase in pressure caused by the reverse shock traveling toward the origin. This increase in pressure creates conditions similar to a fireball, causing a new blast wave to be formed. This is not seen with simulations γ5a4, γ10a4, or γ5a(1/100) which have smaller density jump factors/a density drop, but is seen with simulation γ5a100 (wind termination shock when γ₂ = 5 with a density jump of 100) yielding the conclusion that even at slower speeds, the strong, large density jump still slows the forward shock enough to cause a strong reverse shock to form and in turn, causes conditions similar to a fireball. Henceforth in this paper, we will refer to this newly formed blast wave caused by the strong reverse shock as the "secondary blast wave" for simplicity and clarity.

The formation of the secondary blast wave is also seen in our 2D model of simulation γ10a100. However, the fluid instabilities that follow the encounter with the wind termination shock in 2D slow the secondary blast wave causing it never to encounter the wind termination shock, and in turn, the continuous reflecting back and forth is never realized. This is shown in the bottom right panel in Figure 5, which displays the fluid Lorentz factor over time and the radius of a 1D slice of the 2D simulation γ10a100 along the axis of the jet.

In order to study the effect of the density change on the post-encounter blast wave collimation, we ran simulations γ5a4, γ10a100, and γ5a(1/100) in 2D. Figure 6 shows the results of our 2D simulations. We find that when the blast wave encounters a sudden density increase that is not very large (of factor 4, like that of simulations γ5a4 and γ10a4), the fastest, most energetic fluid will puncture the wall and continue through to the new medium (top left panel in Figure 6). This behavior results in the blast wave recollimating after the encounter as only the fluid with the highest pressure can pass through to the new region and open a pathway for the rest of the blast wave to continue through to the new environment. This is most clearly depicted in the top right panel in Figure 6, which shows the percentage of energy enclosed at various angles over time for simulation γ5a4 in 2D. Before the encounter, the blast wave spreads out, which is seen as the angles encompassing certain percentages of the energy increase over time at times less then T_enc in the top right panel in Figure 6. At the time of the encounter, the most energetic fluid passes through to the new environment and recollimates, which compresses a large percentage of the blast wave energy into a smaller angle. This accounts for the decrease in angle encompassing roughly half of the total energy after the encounter in the top right panel in Figure 6. Also, if there were no encounter, the angle encompassing roughly 80% of the total energy would continue to increase as it did prior to the encounter in the top right panel in Figure 6, but instead, the angle stays relatively constant. This is because the fluid of higher energy (depicted with the colors yellow through red in the top left panel in Figure 6) no longer spreads sideways but tries to travel through the wall to the new medium. The only part of the blast wave that continues to spread sideways is the lowest energy fluid (shown in green and light blue in the top left panel in Figure 6), which results in the continual increase of the angle encompassing 100% of the total energy of the blast wave for simulation γ5a4 in 2D (shown in the top right panel in Figure 6) after the encounter.

Zoom In Zoom Out Reset image size

The dynamics of a blast wave encountering a large density increase (on the order of 100, like that of simulations γ5a100 and γ10a100) are different from that of a smaller density increase. The most energetic fluid punctures the wall and continues on to the next medium, similar to the most energetic fluid encountering a smaller density increase, but instead of opening up a pathway for the rest of the blast wave to move through, a strong reverse shock is formed. This strong reverse shock shocks and spreads the fluid that has not yet passed through to the new medium, forming vortices that are trapped within the stellar wind environment and never pass through to the ISM. The middle left panel in Figure 6 and Figure 7 illustrate these dynamics showing the specific energy of simulation γ10a100 in 2D at two times after the encounter. The middle left panel in Figure 6 depicts the blast wave a short time after the encounter showing that the most energetic fluid has passed through to the ISM, and the reverse shock is traveling back through the fluid in the stellar wind, shocking and spreading the fluid. At the much later time of Figure 7, the highest energy fluid that passed through to the ISM has began to expand as it assimilates to the new medium. Also, the reverse shock has now passed through the fluid behind the density jump, caused the formation of vortices, and created the condition for the secondary blast wave to form. The remains of secondary blast wave can be seen in Figure 7 by the vortices near the horizontal axis at radii around 2 × 10¹⁸ cm.

Zoom In Zoom Out Reset image size

The middle right panel in Figure 6 shows the percentages of energy encompassed in various angles, depicting that the highest concentration of energy is contained within the angle holding the fluid that passes through the density jump, and a very minimal amount of energy is contained within the vortices that have been formed in the stellar wind region. The chaotic behavior of the region in purple and pink in the middle right panel in Figure 6 is from those vortices as well.

Lastly, the dynamics of the simulation with the sudden density drop do not show a recollimation with the increase in the Lorentz factor. The fluid speeds up and as it does, it spreads at a rate slightly faster than in the stellar wind environment. This is shown in the bottom row of Figure 6.

With the understanding of the dynamics of the blast wave encountering a sudden change in density, the next section discusses the resulting light curves.

Figure 8 depicts that there are no observable flares in the flux emitted from a blast wave encountering a small jump, large jump, or large drop. This figure shows the light curves calculated from the 1D simulation and analytical model as well as the light curves calculated from the 2D numerical simulations for on- and off-axis observer angles. For all cases studied in this paper, there were no observable flares.

Zoom In Zoom Out Reset image size

The simulations demonstrate that a blast wave encountering a wind termination shock of a small increase in density (of factor 4) does not result in any disruption in the flux observed at early times. In fact, for a jump of factor 4 (simulations γ5a4 and γ10a4), the light curve gradually transitions to the slope expected in the new environment, confirming the results from the spherical cases of NG and van Eerten et al. (2009). This is seen with the on-axis light curves of the 1D and 2D simulations (top left panel in Figure 8) and with the light curves for a wide range of observer angles of the 2D simulation γ5a4 (top right panel in Figure 8). The encounter is first observable at roughly 46 days for the simulations with γ_enc = 5, and for our simulations of γ_enc = 10, the first time the encounter is observed, T_{obs, enc} ≈ 3 days. These times are much later than any flares observed from Swift and follow from our choice for γ_enc, which was selected for numerical reasons. However, the time at which the encounter is observed is inversely proportional to $\Gamma _{{\rm enc}}^4$:

$$\begin{equation} T_{{\rm obs, enc}} = \frac{T_0 \Gamma _0^2}{4\Gamma _{{\rm enc}}^4}. \end{equation} \tag{ 23 }$$

This scaling is calculated using Equations (6) and (7), and noting that at (and before) the encounter k = 2. It follows that modeling encounters at larger Lorentz factors, of say γ_enc = 25, corresponding to T_{obs, enc} ≈ 0.07 days, would result in the observed encounter being in the range of times that flares are observed (see, e.g., Zhang et al. 2006). Nevertheless, from the top right panel in Figure 9, it is apparent that the light curves for blast waves encountering changes in external density profiles at varying Lorentz factors scale for 1D dynamics. Thus the light curve behavior at the time of the encounter for our numerical simulations of γ_enc = 5  and  10 are indicative of the behavior of encounters at larger Lorentz factors. To ensure that this is accurate, we analytically model encounters at γ_enc = 25, and discuss the light curves from those dynamics throughout this section. The light curves calculated from our analytical model use the same linear radiative transfer algorithm as the light curves calculated from the numerical simulations.

Zoom In Zoom Out Reset image size

For simulation γ5a(1/100), the flux observed follows the same qualitative results as simulation γ5a4 at early times—the observed flux transitions smoothly to the slope of the flux observed solely in the ISM environment after the encounter. Even though the blast wave's velocity increases after the encounter (Figure 4 and top right panel in Figure 5), a flare during the encounter is not observed. Instead, the flux drops more steeply as it evolves toward its lower post-encounter base level. This is shown in the two figures in the bottom row in Figure 8 where the figure on the bottom left depicts the flux observed on axis with the blast wave using the 1D and 2D simulations along with the analytical solution, and the latter displays the flux observed at various observer angles for the 2D simulation.

It is apparent from the light curves of simulations γ5a4 and γ5a(1/100), that as the magnitude of the flux of a blast wave in the ISM environment lowers, there is no rebrightening. Next, we discuss the light curves of simulation γ10a100, which are of a blast wave encountering an ISM environment of much higher density.

With the magnitude of the flux of a blast wave traveling solely in the ISM environment being much higher than the magnitude of the flux from a blast wave in the wind (because the magnitude is a function of the density in that region), one might first assume that a flare must be observed as the flux from the blast wave at the encounter evolves to the new slope and magnitude. However, this is not correct—there is still no rebrightening at the time of the encounter. This is shown in the middle row of Figure 8. The middle right panel in Figure 8 shows the light curves for the 2D simulation γ10a100 plotted against the light curves for the 2D simulation no encounter for various observer angles. This figure illustrates the deviation of the light curve from the blast wave encountering the density jump from the light curve of no encounter at all. The magnitude of the flux decreases as the blast wave travels into the new medium resulting in no flaring activity. The off axis light curve calculated at θ_obs = 0.2 rad in the right panel in the middle row of Figure 8 shows a bump at around 100 days. Although the flux does not increase at this point, it does deviate from the expected slope. This bump in the light curve is observable at a wide range of frequencies, from the optical to the X-ray, and exhibits similar behavior at all frequencies at which it is observed. The bump could be considered a rebrightening, however, it is only seen at observer angles larger than the jet opening angle, and it is unlikely that the prompt emission will be observed at this angle due to the extreme beaming needed for the gamma radiation.

Figure 9 depicts the transition of the light curve as the blast wave changes environments as well as the overall scaling of the light curves. The top left panel in Figure 9 shows analytical light curves for a blast wave traveling in a stellar wind environment with γ_enc = 25 for two scenarios: a density jump of factor 100 (the solid blue line) and a density drop of factor 10,000 (the dashed blue line). These light curves are plotted against the light curve for a blast wave traveling solely in the stellar wind environment (green line), a blast wave traveling solely in the homogeneous medium corresponding to the density jump of factor 100 (red solid line), and a blast wave traveling solely in the homogeneous medium corresponding to the density drop of factor 10,000 (dashed red line). As seen in this figure, the observed flux after the encounter slowly transitions to the slope of the flux observed from the ISM environment, resulting in no sudden increase in flux. The flux does not immediately jump up to the new magnitude because the blast wave has been slowed by the encounter and is not energetic enough to cause a rebrightening. For the case of a drop, the flux observed from a blast wave in the new ISM environment is at a lower magnitude than the flux observed from a blast wave in the wind environment. This causes the flux observed from the blast wave encountering the new environment to decrease, instead of increase and, after the encounter, to evolve in the new environment.

In addition, this figure shows an interesting feature, that can also be seen in Figure 8, of a shallowing of the light curve. The light curve transitions from the slope of a blast wave in a stellar wind to the less steep slope of a blast wave in an ISM. This late-time shallowing of the light curve is not the turnover of a steep decay into a plateau of the canonical light curve, but a different transition at later times. This is an observed feature of GRB afterglows (e.g., Evans et al. 2009; Li et al. 2012; Zaninoni et al. 2013) and a change in circumburst medium is one of the few ways for the light curve of the blast wave to transition to a new slope.

From the light curves of simulations γ5a4, γ10a100, and γ5a(1/100), as well as the light curves from the simulations not shown here (γ5a100 and γ10a4), we conclude that there is no rebrightening at times on the order of minutes to days for a blast wave encountering a change in circumburst medium. The simulations show that when the magnitude of the flux observed in the ISM is lower than the flux observed in the stellar wind environment at the time of the encounter (top left panel in Figure 9) the blast wave conforms to the new slope and there is no rebrightening. When the magnitude of the flux observed from a blast wave in an ISM environment is higher than the magnitude of the flux observed from a blast wave in the stellar wind environment, the blast wave still smoothly transitions to the new environment. In all cases, though, the flux observed from the blast wave after the encounter is lower than the flux observed from a blast wave that has only traveled in the ISM.

To test whether the light curves for our simulations are an accurate depiction of a blast wave encountering a circumburst density change, or if they only represent the light curves for encounters at fluid Lorentz factors of γ ⩽ 10, with changes in density of factors ⩽100 and with θ₀ = 0.1, we analytically solve for the light curves of various blast wave encounters for various encounter Lorentz factors, γ_enc, density change factors, and jet half opening angles, θ₀. The top right panel in Figure 9 shows light curves for 1D analytical models of blast waves traveling in a stellar wind environment that encounter a density change of factor 10 followed by a homogeneous medium with the encounters happening at varying Lorentz factors. When the factor by which the density changes is kept fixed, the shape of the light curve remains independent of the actual time of the encounter and light curves for different encounter times can simply be shifted to match, as can be seen from the figure.

The bottom left panel in Figure 9 shows the light curves for blast waves traveling in the stellar wind environment that encounter a density change of varying magnitudes at γ_enc = 25. We show here density changes ranging from 10⁻⁸ to roughly correspond with the density drops used in Mesler et al. (2012), to 10³. For larger density jumps, the blast wave quickly becomes non-relativistic, causing our analytical model to become less accurate. Lastly, the bottom right panel in Figure 9 shows light curves for blast waves traveling in a stellar wind environment that encounter a density jump of factor 10 at γ_enc = 25 for various jet half opening angles, θ₀. The only light curve in this figure that has some resemblance of a flare is that of the blue line, which has an opening angle very close to γ_enc ∼ 1/θ₀, meaning that for any flare to possibly occur around the time of the encounter, γ_enc = 1/θ₀. This is the setup of our simulations γ10a4 and γ10a100: γ_enc = 10, and θ₀ = 0.1, and no flare is observed in 1D or 2D.

From this analysis, we conclude that for a blast wave traveling in a stellar wind environment that encounters a sudden change in circumburst environment, there is no observable rebrightening regardless of the size of the jump, drop, or fluid Lorentz factor at the time of the encounter that we have studied. There is a bump in the light curve at the observation angle of θ_obs = 0.2 rad for the light curve of simulation γ10a100 in 2D, but θ_obs > θ₀ and thus this is not a likely indication of an observable flare. We do not see a flare from the forward shock, nor do we see a flare caused by the reverse shock at times corresponding to the afterglow flares seen by Swift.

We study a single spherically symmetric change in the circumburst environment. However, the actual environment may have a more complicated structure with instabilities (e.g., van Marle & Keppens 2012) or thin circumburst shells (Mesler et al. 2012). A study of asymmetric structures (i.e., density clumps) around the progenitor would require three-dimensional simulations and is beyond the scope of this paper. We do look at the interaction with circumburst shells, and our results are discussed further in Appendix B. In this case also, we have been unable to reproduce the flares reported by Mesler et al. (2012).

The simulations of a density jump did show flares at times long after the encounter occurred (at times much later than plotted in Figure 8). At such late times (on the order of 10⁴ days), these flares are beyond X-ray observations, but may be of interest at radio frequencies. There are two contributors to these flares—the fireball conditions created from the strong reverse shock that causes a secondary blast wave to form, shown in the bottom row of Figure 5, and the encounter of the receding jet with the wind termination shock. We discuss the contribution of the secondary blast wave first, followed by a discussion of the receding jet.

Figure 10 shows the very late time flare observed from the 1D and 2D simulations of simulation γ10a100. Comparing the time at which the flare begins for simulation γ10a100 in 1D and 2D, shown in Figure 10, with the time and radius of the secondary blast wave shown in the bottom left panel in Figure 5, it is apparent that the secondary blast wave is a cause of the rebrightening. This is confirmed by the left panel in Figure 10, which shows the light curve calculated for simulation γ10a100 in 1D using the entire blast wave (blue line), against the light curve observed from the parts of the blast wave at radii larger than the encounter radius (green line). The left panel in Figure 10 clearly illustrates that the flare is caused by the fluid inside the encounter radius at late times, and the most energetic fluid in the 1D simulation that is at radius smaller than R_enc at late times is the secondary blast wave.

Zoom In Zoom Out Reset image size

The left panel in Figure 10, although explaining the flare for the 1D case of simulation γ10a100, does not explain the small flares seen in the 1D conical outflow for simulation γ5a4 nor the flares for the 2D simulations γ5a4 and γ10a100.

To fully understand the cause of this late time flare, we must emphasize that it is seen at a wide range of observer angles and frequencies (right panel in Figure 10) at the same observer time, and is observable at the same time as the flare from the secondary blast wave. The secondary blast wave cannot be the cause of the flare for the 2D simulation because it is quickly slowed by the surrounding material (seen in the bottom right panel in Figure 5), and is also not a characteristic of simulation γ5a4 in 1D or 2D. Moreover, the flare observed in simulation γ5a4 is observable in the 2D light curves and the 1D light curve that assumed conical outflow, but not in the 1D light curve assuming spherical outflow.

The cause of the late time flare is the receding jet that encounters the sudden density jump (but not from a density drop). The reason this event seemingly occurs at the same time as the secondary blast wave is because the sum of the distance traveled by the reverse shock and the flux from the secondary blast wave is the same as the distance traveled by the flux observed from the receding jet encountering the density jump. In short, the density jump results in a sudden increase in flux traveling in the opposite direction of the blast wave itself, meaning a spike in the light curve will not be observed at the time of the encounter because the light from the encounter is traveling away from the observer. This is confirmed by the green lines in the right panel of Figure 10, which shows that the light curve observed only from the forward jet does not result in a very late time flare whereas the light curve observed from both the forward and the receding jet does have a flare. We obtain the same results when analyzing simulation γ10a100 at different frequencies, and by analyzing the very late time flare of simulation γ5a4 in 2D (and 1D conical outflow). The reason for the differing behavior of the green line and the red and blue lines of the right panel of Figure 10 is that the radio frequency corresponds to a different spectral regime than the visible and X-ray frequencies at early times. The radio frequency at early times is below the synchrotron peak critical frequency, ν_m, meaning the light curve scales with t⁰ in the wind scenario (e.g., Granot & Sari 2002). This early time radio light curve behavior is showing a different spectral regime with ν_m passing through the observed band at the turnover.

The receding jet, or counter jet, has been known to be observable at very late times (Granot & Loeb 2003; Li & Song 2004; Zhang & MacFadyen 2009; Wang et al. 2009), with its magnitude being a function of the density structure into which the forward and receding jets are traveling (Wang et al. 2009). What differentiates the flares from the receding jet shown in this paper and the visibility of the receding jet in previous work is that this flare is from the encounter of the receding jet with the sudden higher density region, as opposed to the counter jet being observable as it becomes non-relativistic and isotropic (Li & Song 2004) or as it travels into a differing circumburst medium than for the forward jet, or as both the forward and receding jet travel into a highly dense medium (Wang et al. 2009). In addition, the flare in this paper is observable at a wide range of observer angles, in contrast to the prediction that the receding jet is only observed at optimal angles (Granot & Loeb 2003). The very late time flare observed in this work is from the sudden decrease in the Lorentz factor of the jets at the encounter, causing a flare to be observed as opposed to observing the receding jet itself.

From our simulations, the only flare calculated is that at much later times and much lower magnitudes in flux than those observed by Swift. The only possible flare on the time scale of observed flares that we are able to observe with a sudden change in the circumburst medium is the small bump in the light curve seen at around 30 days at θ_obs = 0.2 rad in the middle right panel in Figure 8, but this occurs only for θ_obs > θ₀ and thus is unlikely to be observed.