Dynamical Energy Dissipation of Relativistic Magnetic Bullets

The simulation box is in the external medium rest frame. For the external medium, we assume that the magnetic field and the number density are homogeneous and set as $1\ \mu$G and $n_{0}=1\ \rm{cm}^{-3}$, respectively. The temperature is also fixed as 1 MeV to stabilise the simulations, but it is sufficiently cold compared to proton rest mass energy.

From the definition of the magnetisation $\sigma$, we can write the initial magnetic field of the ejecta as

where $\sigma_{0}$ and $\Gamma_{0}=10$ are the initial magnetisation and the Lorentz factor of the ejecta, respectively. We define the deceleration radius $R_{\rm dec}$ as (Sari & Piran, 1995)

At this radius, the initial ejecta energy $E_{0}$ is expected to be converted into the energy of the shocked external medium. The inner edge of the ejecta is initially set at $R_{0}=0.11R_{\rm dec}$. We fix the initial width of the ejecta as $\Delta=R_{\rm dec}/\Gamma^{2}_{0}$. With this setup, the initial width $\Delta=0.01R_{\rm dec}$ is thicker than $R_{0}/\Gamma_{0}^{2}=0.01R_{0}$ (thick shell model), which can be realised by a significantly longer duration of the central engine activity. In our simulations, we consider only the thick shell model because the thin shell model requires a much better spatial resolution. We assume that the ejecta has constant magnetisation $\sigma_{0}$, temperature $T_{0}=100$ MeV, and Lorentz factor $\Gamma_{0}$, but the mass density follows as $\rho=\rho_{0}(r/R_{0})^{-2}$ in the initial condition. The initial energy of the ejecta is set to be $E_{0}=10^{50}$ erg, which is calculated with

Using Eqs.(1) and (2), we can write the normalisation of the mass density $\rho_{0}$ as

The deceleration time $t_{\rm dec}\simeq(R_{\rm dec}-R_{0})/c$ is $1.6\times 10^{6}$ s in our setup.

The simulation box size is $[0.1R_{\rm dec},5R_{\rm dec}]$. The outer boundary is set to be an open boundary. The inner boundary is set to be an injection boundary, where we inject a plasma with a significantly low luminosity with $\Gamma=10$ continuously. This injection can prevent the regions just behind the ejecta from becoming too dilute (avoiding a numerical floor). The number of static cells is $5\times 10^{5}$, effectively $8\times 10^{6}$ via the AMR procedure (see Appendix A). The initial width of the ejecta is resolved by almost $10^{3}$ static cells and effectively almost $10^{4}$ static cells.

Just after the start of the calculation, two rarefaction waves start to propagate from the front and the rear edges of the ejecta, respectively. As the rarefaction waves propagate, pressure gradients are formed inside the ejecta. One of them accelerates the front edge of the ejecta, while the other decelerates the rear edge of it (Granot et al., 2011). The rarefaction wave accelerates the ejecta faster for higher $\sigma_{0}$. As the forward shock propagates into the external medium, the increased thermal pressure of the shocked external medium generates the reverse shock to decelerate the ejecta.

Figure 1 shows snapshots of $\Gamma$, $\sigma$, and the total pressure around the shocks at $t=3\times 10^{5}$ s. The mass density profiles are not shown for simplicity. From right to left, we can see the forward shock (FS), where $\Gamma$ is jumping, the contact discontinuity (CD), where $\sigma$ is jumping, the reverse shock (RS), where $\Gamma$ and $\sigma$ are jumping again, the rarefaction wave accelerating the ejecta (R1), and the rarefaction wave decelerating the ejecta (R2), which forms the tail of the ejecta. Because the propagation speed of the magneto-sonic waves depends on magnetisation, the positions of R1 and RS are sensitive to $\sigma_{0}$. In Figure 1, for $\sigma_{0}=1$, the wave R1 has not yet reached the rear edge of the ejecta so that the initial flat structures of $\Gamma$ and $\sigma$ remain at 0.276–0.281 $R_{\rm dec}$, and the RS position is still close to the head of the ejecta. For $\sigma_{0}=10$, the wave R1 has already reached the rear edge of the ejecta, which eliminates the initial flat structure, and the RS position is close to the rear edge of the ejecta.

The reverse shock propagates into the ejecta. We define the time at which the reverse shock crosses the ejecta as $t_{\Delta}$. In a hydrodynamic case ($\sigma_{0}=0$), Sari & Piran (1995) provided an analytical formula of $t_{\Delta}$ as

where $\Delta$ is the shell thickness in the external medium rest frame. For arbitrary magnetised cases, Zhang & Kobayashi (2005) estimated the reverse shock crossing time $t_{\Delta}(\sigma)$ as a function of $\sigma$, which is approximated as (Giannios et al., 2008)

As shown in Figure 2, the reverse shock crossing time $t_{\Delta}$ well agrees within $\sim$ 10% with Eq. (13), where the normalisation $t_{\Delta}(0)$ is given by our simulation data of $\sigma_{0}=0$. A simulation with 10 times coarse resolution shows that the crossing time $t_{\Delta}$ becomes slightly longer. This is due to the numerical expansion of the edge of the ejecta due to the insufficient resolution (see Appendix A). In GRB afterglows, the reverse shock crossing time corresponds to the time of the peak flux from the reverse shock (Zhang & Kobayashi, 2005). From our simulations, the early peak time of the reverse shock emission for high-$\sigma$ ejecta is confirmed.

Until the deceleration radius, the forward shock can be approximated to freely expand with constant speed. Beyond $R_{\rm dec}$, the forward shock evolves adiabatically with decaying $\Gamma\beta$ approximated by the Blandford-McKee (BM) solution (Sari & Piran, 1995; Zhang & Kobayashi, 2005; Giannios et al., 2008). The BM solution is a self-similar solution for an ultra-relativistic outflow from a point-like explosion. According to this solution, the Lorentz factor just behind the shock front $\Gamma_{\rm FS}$ evolves as $\Gamma_{\rm FS}\propto t^{-3/2}$ and the downstream profile of $\Gamma$ is expressed as (Blandford & McKee, 1976)

where $R_{\rm FS}$ is the radius of the forward shock. Our purpose in this section is to check the asymptotic behaviour of our results at the late phase of expansion and discuss the difference in the case of the ejecta with a finite width.

Figure 3 shows the time evolution of $\Gamma\beta$ of the forward shock, which is estimated at the radius where the gas pressure becomes maximum behind the shock front. The black line shows the frequently used approximation for the forward shock 4-velocity, which expresses the transition from the free expansion to the BM solution at the deceleration time $t_{\rm dec}$. Because the mass increase of the shocked external medium leads to the gradual deceleration of the forward shock, a clear break due to the transition from the free expansion to the BM phase does not appear at $t=t_{\rm dec}$. However, Figure 3 shows that late breaks of $\Gamma\beta$ appear at $t_{\rm FS}\simeq 4\times 10^{6}$ s, which is almost independent of the initial magnetisation of the ejecta. As shown in the bottom panel of Figure 4, before the break time $t<t_{\rm FS}$, the profile of $\Gamma\beta$ is flatter than the BM solution just behind the forward shock because of the effect of the initial condition. This flat structure is eliminated gradually as the rarefaction wave R2 approaches the forward shock front. At the break time $t=t_{\rm FS}$, the rarefaction wave just reaches the forward shock front. After the break time $t>t_{\rm FS}$, the profile of $\Gamma\beta$ is almost the same as the BM solution Eq.(14). This break of the $\Gamma\beta$-evolution may affect lightcurves in GRB afterglows.

The time evolutions of the different energy components are shown in Figure 5. We calculate the kinetic energy component $E_{\rm kin}$, the thermal energy component $E_{\rm th}$, and the magnetic energy component $E_{\rm mag}$ from Eq. (4) as,

respectively. There are 4 characteristic times as shown in Figure 5: the deceleration time $t_{\rm dec}$ defined with Eq. (9), the reverse shock crossing time $t_{\Delta}$, the break time $t_{\rm FS}$, and the R2 crossing time at the CD $t_{\rm CD}$ defined later.

For $\sigma_{0}=0$, at the analytical deceleration time $t_{\rm dec}$, the kinetic energy and the thermal energy are roughly in equipartition. The maximum thermal energy is achieved at around $4\times 10^{6}$ s, which corresponds to the break time $t_{\rm FS}$ (the RS crossing time at the FS). After that, the thermal energy gradually decreases by adiabatic expansion.

For $\sigma_{0}=0.1$ and $\sigma_{0}=1$, the thermal and the kinetic energy evolutions are similar to the case of $\sigma_{0}=0$. The magnetic energy increases by the reverse shock compression until the reverse shock crossing time $t_{\Delta}$. At the deceleration time $t_{\rm dec}$, the thermal energy is almost equipartition with the kinetic energy for $\sigma_{0}=0.1$, while with the magnetic energy for $\sigma_{0}=1$. This means that almost half of the initial energy of the ejecta is converted to the thermal energy of the shocked external medium at $t_{\rm dec}$. For $\sigma_{0}=0.1$, the maximum thermal energy is achieved at $t_{\rm FS}$. While for $\sigma_{0}=1$, the thermal energy reaches maximum slightly before $t_{\rm FS}$. After that, the adiabatic expansion reduces the thermal energy gradually.

For $\sigma_{0}=10$, as the reverse shock is weak due to the high $\sigma$, the magnetic energy is almost constant in the initial phase. At the deceleration time $t_{\rm dec}$, the magnetic energy and the thermal energy are in equipartition. As shown in Figure 6, as the rarefaction wave R2 approaches the contact discontinuity, the flat structure of $\sigma$ in the ejecta disappears. We define the time at which the R2 crosses the CD as $t_{\rm CD}$. At $t_{\rm dec}\leqq t\leqq t_{\rm CD}$, only the magnetic energy decreases drastically, while the kinetic and the thermal energy increase. This implies that the magnetic energy is converted into the kinetic energy of the shocked external medium by the magnetic pressure, and a part of that energy is dissipated into the thermal energy by the forward shock (Granot et al., 2011; Komissarov, 2012). The maximum thermal energy is achieved at $t_{\rm CD}$. After that, the adiabatic expansion reduces the thermal energy gradually.

In GRB afterglow models, the deceleration time $t_{\rm dec}$ is assumed as the peak emission time from the forward shock (Zhang & Kobayashi, 2005). We can estimate the initial $\Gamma$ from this peak time. However, for the thick shell cases, our simulation implies that the emission has peaked at $t_{\rm FS}$ for $\sigma_{0}<1$ and at $t_{\rm CD}$ for $\sigma_{0}>1$. The difference in the magnetisation may affect the emission from not only the reverse shock but also the forward shock. Since the reverse shock dynamics might be suffered from the Rayliegh-Taylor instability at the contact discontinuity as discussed in Duffell & MacFadyen (2013) by their 2D simulation, we will further study this impact on magnetic energy dissipation by our 2D simulations in future.

The simulation box size is $[R_{0},20R_{0}]$, where $R_{0}=10^{16}\ \rm{cm}$. The inner boundary is set to be an injection boundary. The outer boundary is set to be an open boundary. The number of static cells is $2\times 10^{5}$, effectively $3.2\times 10^{6}$ via the AMR procedure (see Appendix A). For the initial condition, wind fills the simulation box with the density profile of $\rho\propto r^{-2}$.

We inject ejecta and winds alternately from the inner boundary. The Lorentz factor of the ejecta $\Gamma_{\rm jet}$ are randomly following a cutoff-Gaussian probability distribution: a Gaussian with maximum and minimum cutoffs. The typical Lorentz factor of ejecta is around 10 for Blazars (Ghisellini et al., 1993) (GRB jets have so much higher value that simulations for such ultra-relativistic outflows are difficult). Thus, we choose the mean value of $\Gamma_{\rm jet}$ is 10 and its dispersion is 5. The cutoffs are at $10\pm 5$. We define the minimum injection duration of the ejecta as

The injecting duration of the ejecta $t_{\rm jet}$ is randomly determined with a cutoff-Gaussian probability distribution with the mean value of $3t_{0}$, the dispersion of $2t_{0}$, and the cutoffs at $3t_{0}\pm 2t_{0}$. The initial magnetisation and temperature of the ejecta are fixed as $\sigma_{\rm jet}=10$ and 100 MeV, respectively. The total luminosity of the ejecta is fixed as $L_{\rm jet}=10^{45}$ erg $\rm{s}^{-1}$, which means that the density changes according to changing $\Gamma_{\rm jet}$.

For winds, the velocity, magnetisation, and temperature are fixed as $0.1c$, $\sigma_{\rm wind}=10^{-10}$, and 100 MeV, respectively. To reduce the computational cost due to AMR refinement procedures, we assume the mean injection duration of winds $t_{w}$ is significantly longer than $t_{0}$. The injecting duration of the winds $t_{\rm wind}$ is determined by a cutoff-Gaussian probability distribution with the mean value of $t_{w}=R_{0}/c$, the dispersion of $0.1t_{w}$, and cutoffs at $t_{w}\pm 0.1t_{w}$.

The total luminosity of the winds $L_{\rm wind}$ is fixed in a simulation run. We parameterise $L_{\rm wind}$ as

If the ejecta energy is totally transferred into the swept-up wind energy, we can write

For $t_{\rm jet}=3t_{0}$, $t_{\rm wind}=t_{w}$ and $\Gamma_{\rm jet}=10$, we obtain

If $\eta$ is much smaller than this critical value $\eta_{\rm crit}$, the wind luminosity may be not significant to dissipate the energy of the magnetic bullet. Conversely, we expect that most of the ejecta energy can be transferred into the energy of the shocked wind for $\eta>\eta_{\rm crit}$. To investigate the $\eta$-dependence on the outflow structure and the dissipation efficiency, we carry out four simulations listed in Table 1. We probe the difference due to different $\eta$ with models A and B, and the different $t_{\rm wind}$ with models C and D.

Snapshots of the outflows at $2\times 10^{6}$ s ($\sim 10^{3}t_{0}$) are shown in Figure 8, 9, 13, and 14 for models A, B, C, and D, respectively. The influence of the initial condition has already disappeared at this stage so the simulations are in quasi-steady states. High-temperature regions are the sites where energy dissipation occurs. Almost all of the efficient dissipation regions are concentrated around the discontinuities of the Lorentz factor. This means that efficient shock dissipation occurs through interactions between ejecta and low-magnetised plasma. Because the outflow structure is dependent on wind models, we explain the difference between models A and B as follows (see Appendix B for models C and D).

Figure 8 shows the outflow structure for model A ($\eta=10^{-4}<\eta_{\rm crit}$). The high-$\sigma$ ejecta between the low-pressure winds is gradually accelerated. This process is equivalent to the impulsive acceleration discussed in Granot et al. (2011) (see also Komissarov, 2012). The Lorentz factor of the front edge of the ejecta increases as $\Gamma\sim\Gamma_{\rm jet}(R/R_{0})^{1/3}$ until the rarefaction wave from the rear edge will catch up with it. Behind the ejecta, there exist long rarefaction tails with moderate magnetisation. The ram pressure of the winds is weak so the winds are compressed by the rarefaction tails.

The bottom left panel of Figure 8 shows the interaction between the head of an inner ejecta ($<5.87\times 10^{16}$ cm), a wind (the narrow region between high-$\sigma$ regions), and the rear tail of an outer ejecta ($>5.87\times 10^{16}$ cm). The high $\Gamma$ and temperature of the wind indicate that the magnetic energy of the inner ejecta is transferred into the kinetic and thermal energies of the wind (see Section 3 in detail). The outer jump of $\Gamma$ corresponds to the shock front propagating in the rarefaction tail of the outer ejecta. The apparent low temperature in the shocked wind region is due to thermal pressure acceleration, which transfers the thermal energy of the shocked wind into the kinetic and the thermal energy of the rarefaction tail of the outer ejecta. As shown in the bottom middle and right panels of Figure 8, the wind components are squashed at larger radii, so that they are below the resolution in our simulations.

In model A, the shocked rarefaction tails, where $\sigma<1$ and $\Gamma>10$, are regarded as kinetically dominated relativistic hot outflows, which can be efficient emission sites.

On the other hand, Figure 9 shows a very different outflow structure for model B ($\eta=10^{-2}>\eta_{\rm crit}$). It shows that the Lorentz factor of the ejecta is no larger than the value at injection. This is due to the heavy mass of the winds. The interaction between the ejecta and winds leads to a significant dissipation of the kinetic energy of the ejecta. A fraction of ejecta at outer radii has a high Lorentz factor due to impulsive acceleration. Differently from model A, the $\sigma$ values of the rarefaction tails are mostly higher than one. This is because the heavy winds prohibit rarefaction tails from expanding, which reduces their magnetisation.

The bottom left panel of Figure 9 shows a shock in the wind. Because the shock is mildly relativistic, the temperature is lower than the cases in model A. In such shocked regions, the emission efficiency may not be so high. The shock waves in the winds tend to catch up the rarefaction tail of the preceding ejecta at inner radii. Such shocks propagate in the rarefaction tail eventually. One of such shocks is shown in the bottom middle panel of Figure 9. Since the strong ram pressure of the winds prevents the preceding rarefaction tail from expanding enough to reduce its magnetisation, the shocks are weaker than those in model A, leading to less energy dissipation. However, we find succeeding acceleration of such shocks by magnetic pressure even in the rarefaction tails. Some of such shocks collide with the preceding wind again (see the bottom right panel of Figure 9), which leads to the high efficiency of the energy dissipation at outer radii. Thus, kinetically-dominated relativistic hot outflows tend to be generated at outer radii for model B. We will further discuss this radial dependency in the following section.

In the previous section, we find that interactions between magnetised ejecta and weakly magnetised plasma (winds or rarefaction tails) can produce kinetically dominated relativistic hot outflows. In this section, we estimate the magnetic dissipation efficiency. We calculate the thermal luminosity of outflows passing some radius $R$ as

Because our interest is the thermal energy of the kinetically dominated relativistic outflows, we estimate the thermal energy for only outflows with $\sigma<1$ and $\Gamma>5$ as

The result is shown in Figure 10. In all cases, the integrated thermal energy increases almost as a monotonic function of time. The maximum average dissipation efficiency $f_{\rm max}$ corresponds to the slope of the dotted lines in Figure 10 is calculated by

where $L_{\rm in}$ is the time-averaged total injection luminosity of the outflow. Without the condition of $\sigma<1$ and $\Gamma>5$, we find that the efficiency is about 18% for models A, B and C, while 14 % for model D. Irrespective of the high $\sigma$ values at injection, significantly high efficiencies are successfully obtained in our simulations.

For model A, the integrated energy increases smoothly with an efficiency of 15% independent of the radius. This is because almost all of the dissipation occurs through interactions between ejecta and rarefaction tails at any radius. On the other hand, model B shows that dissipation at $R=10R_{0}$ is more effective than that at $R=5R_{0}$. Because most outflows are mildly relativistic ($\Gamma<5$) as shown in Figure 9, the integrated energy increases step-wise at the times when a high $\Gamma$ outflow arrives. As we pointed out in Section 4.2, high $\Gamma$ shocks in the winds tend to be produced after experiencing magnetic pressure acceleration in the preceding rarefaction tails at outer radii.

For models C and D (intermediate case for $\eta$), the average dissipation efficiency becomes larger at the inner radii ($R=5R_{0}$). The wind regions are below our numerical resolution at outer radii (see Appendix B), so low magnetised outflows exist only at inner radii, where the magnetic energy can be efficiently dissipated. The maximum dissipation efficiency of model C is higher than that of model D. This is because a short wind duration suppresses the energy dissipation via the interaction with the winds.

To clarify the behaviour of the dissipation, we estimate the time-dependent dissipation efficiency at $R=5\times 10^{16}$ cm using the following definition

where $\delta t$ is set to be sufficiently long as $2\times 10^{6}$s$\sim 10^{3}t_{0}$ to produce a relatively smoothed shape of $f$. The results are plotted in Figure 11. After reaching the quasi-steady state ($\sim 10^{7}$ s), the time-dependent dissipation efficiency is almost stationary for model A (13 $\sim$ 17%) and D (7 $\sim$ 8%). On the other hand, it fluctuates significantly for model B (0 $\sim$ 15%). The behaviour in model C is intermediate (10 $\sim$ 17%). In spite of the low average dissipation efficiency for model B, the flare-like dissipation may be favourable to explain the gamma-ray flares frequently seen in blazar observations.