Magnetized Reverse Shock: Density-fluctuation-induced Field Distortion, Polarization Degree Reduction, and Application to GRBs

TL;DR

This paper investigates how density fluctuations in the environment of gamma-ray bursts can distort magnetic fields in reverse shocks, reducing polarization, and proposes a method to estimate the outflow's magnetization parameter from polarization data.

Contribution

It introduces relativistic MHD simulations of magnetized jets in clumpy media to quantify polarization reduction and relates it to the outflow's magnetization parameter, σ.

Findings

Density fluctuations cause a power-law reduction in polarization degree.

A relationship between polarization reduction and σ allows estimation of outflow magnetization.

Simulations demonstrate the impact of environmental inhomogeneities on magnetic field order.

Abstract

The early optical afterglow emission of several gamma-ray bursts (GRBs) shows a high linear polarization degree (PD) of tens of percent, suggesting an ordered magnetic field in the emission region. The light curves are consistent with being of a reverse shock (RS) origin. However, the magnetization parameter, $\sigma$, of the outflow is unknown. If $\sigma$ is too small, an ordered field in the RS may be quickly randomized due to turbulence driven by various perturbations so that the PD may not be as high as observed. Here we use the "Athena++" relativistic MHD code to simulate a relativistic jet with an ordered magnetic field propagating into a clumpy ambient medium, with a focus on how density fluctuations may distort the ordered magnetic field and reduce PD in the RS emission for different $\sigma$ values. For a given density fluctuation, we discover a clear power-law relationship between the relative PD reduction and the $\sigma$ value of the outflow. Such a relation may be applied to estimate $\sigma$ of the GRB outflows using the polarization data of early afterglows.

Tables

Table 1. Table 1 The initial parameters of the simulation models

Group 1: ${\Gamma }_4=10$
Models	${\sigma }_4$	$B_{y,4}$	${\Gamma }_4$	$n_4$	$n_1$	$P_4$	$P_1$
0	0	0	10	100	1	1	0.01
1	0.1	31.62	10	100	1	1	0.01
2	${10}^{-2}$	10	10	100	1	1	0.01
3	${10}^{-3}$	3.162	10	100	1	1	0.01
4	${10}^{-4}$	1	10	100	1	1	0.01
5	${10}^{-5}$	0.3162	10	100	1	1	0.01
Group 2: ${\Gamma }_4=20$
Models	${\sigma }_4$	$B_{y,4}$	${\Gamma }_4$	$n_4$	$n_1$	$P_4$	$P_1$
6	0.1	63.24	20	100	1	1	0.01
7	${10}^{-2}$	20	20	100	1	1	0.01
8	${10}^{-3}$	6.324	20	100	1	1	0.01
9	${10}^{-4}$	2	20	100	1	1	0.01
10	${10}^{-5}$	0.6324	20	100	1	1	0.01

Equations

$$\frac{n_4}{n_1}=\frac{{\Gamma }_{21}^2-1}{{\Gamma }_{34}^2-1},$$

Full text

1. Introduction

The jet composition and energy dissipation, particle acceleration and radiation mechanisms of GRBs have not been identified and thus are subject to debate (Kumar & Zhang, 2015, for a recent review). Leading models include the matter-dominated fireball shock models (e.g. Goodman, 1986; Paczynski, 1986; Rees & Meszaros, 1994; Mészáros & Rees, 2000) and the Poynting-flux-dominated magnetic dissipation models (e.g. Lyutikov & Blandford, 2003; Zhang & Yan, 2011; McKinney & Uzdensky, 2012; Deng et al., 2015, 2016). Growing evidence suggests that the GRB jets carry a dynamically important magnetic field, so that magnetic reconnection may play an important role in producing observed GRB prompt emission and afterglow emission.

One piece of such evidence comes from the observations of the early optical afterglow of GRBs, which is dominated by a rapidly decaying ($\propto t^{-2}$) segment early on from the external reverse shock (RS) (Akerlof et al., 1999; Mészáros & Rees, 1999; Sari & Piran, 1999; Gomboc et al., 2008; Gao et al., 2015). Modeling such a component suggests that the RS is more magnetized than the forward shock (FS), so that the GRB jet is more magnetized than the ambient medium (Fan et al., 2002; Kumar & Panaitescu, 2003; Zhang et al., 2003; Zhang & Kobayashi, 2005). This is consistent with the scenario that prompt emission is powered by magnetic reconnection and turbulence with a moderate magnetization parameter $\sigma$ in the prompt emission region (Zhang & Yan, 2011). On the other hand, when the jet is decelerated by an ambient medium (when the RS propagates through the ejecta), the $\sigma$ value of the ejecta cannot be very high. Otherwise, the RS would be weak or suppressed (Zhang & Kobayashi, 2005; Mimica et al., 2009, 2010; Mizuno et al., 2009), contrary to the observations.

Direct proof of such a scenario comes from the polarization observations of the early optical emission of some GRBs. Steele et al. (2009) discovered an $\sim 10\%$ polarization degree (PD) of the early RS optical emission from GRB 090102. Mundell et al. (2013) for the first time measured the “polarization light curve” of GRB 120308A, which shows an $\sim 28\%$ PD at 4 minutes after the GRB, which is reduced to $\sim 16\%$ PD over the subsequent 10 minutes as the FS emission contribution becomes stronger. These observations suggest that an ordered magnetic field in the RS region exists (Lan et al., 2016). However, the $\sigma$ value of the RS emission cannot be easily diagnosed based on the PD and light-curve data.

For a relativistic outflow, strong turbulence is likely developed in the flow due to various perturbation mechanisms (e.g. density clumps in the ambient medium or intrinsic irregularity within the flow). If $\sigma\ll 1$, even if an ordered magnetic field configuration may exist initially, turbulence may quickly randomize the field configurations, which would greatly reduce the ordered field component and the PD of synchrotron radiation. As a result, based on the observed PDs of GRBs 090102 and 120308A, one may be able to place a lower limit on the $\sigma$ values of the outflows. To do this, detailed numerical simulations are needed.

The main goal of this work is to investigate the resilience of the ordered magnetic fields in a relativistic jet against the density-fluctuation perturbations from the ambient medium as well as the resultant reduction of the PD of synchrotron radiation due to such perturbations. To achieve this goal, we have chosen to solve a simplified problem, i.e. to simulate how an ambient density fluctuation alters the initially parallel magnetic fields in a relativistic jet. A more realistic problem would invoke a partially ordered magnetic field configuration (after the GRB prompt emission phase, which resulted from dissipation of magnetic fields) in a conical jet (in which a Rayleigh-Taylor (RT) instability would play a role to further destroy ordered magnetic fields). In any case, our results present a “lower limit” of the relative PD reduction for a certain $\sigma$ value, which may be potentially used to constrain the “lower limit” of $\sigma$ in a GRB outflow once a PD is measured. The simulations are presented in Section 2. The simulation results, especially a relationship between the relative PD reduction and the initial $\sigma$ value of the ejecta are presented in Section 3. The results are summarized in Section 4 along with a discussion on their applications to the GRB problem.

2. Problem setup

We use the relativistic MHD code “Athena++” (White et al., 2016) to simulate the dynamical evolution of the outflow, and the “3DPol” multi-zone polarization-dependent ray-tracing radiation code (Zhang et al., 2014) to model the PD of synchrotron radiation in the RS region.

We perform our simulations in the “lab” frame. We set up a jet that propagates in the $z$ direction with a relativistic speed into an ambient medium at rest (Fig. 1). In order to isolate the problem in testing the resilience of ordered magnetic fields against density perturbation, we adopt several simplified assumptions. First, the jet has a planar geometry, i.e., the cross-section of the jet does not expand with time (as is the case of a conical jet). This would reduce other factors, such as the RT instability (Duffell & MacFadyen, 2013), to deform ordered magnetic fields. Second, we assume a parallel magnetic field configuration within the simulation domain in the unshocked region. In principle, after prompt emission, the magnetic field configuration is expected to be somewhat distorted. Adopting a simple magnetic field configuration, on the other hand, can remove other effects and allow us to concentrate on the density perturbation effect on destroying the ordered fields. The global magnetic field configuration of a GRB ejecta may not be approximated as the simple parallel configuration we use. However, a GRB observer only observes emission from a $1/\Gamma$ cone, in which a toroidal magnetic field configuration (e.g. Lyutikov et al., 2003; Deng et al., 2015) can be approximated as being approximately parallel (e.g. Lazzati, 2006).

For the planar geometry problem, it is convenient to simulate the jet propagation problem using a 2D rather than a 3D box111Even though the simulation box is essentially 2D, all the vectors are solved in 3D.. The effective resolution of the simulations is $\Delta y\approx\Delta z\approx 0.0012L_{0}$, where $L_{0}$ is the length normalization factor between code units and physical units. The entire simulation frame has $20L_{0}$ in the $z$ direction and $0.156L_{0}$ in the $y$ direction. Since the simulation box may be regarded as a patch inside the $1/\Gamma$ cone of the entire thin shell in both of the radial and the transverse directions, we set up a periodic boundary condition in the $y$ direction and an inflow boundary condition in the lower boundary of the $z$ direction.

As the jet runs into the medium, four regions are developed: (1) unshocked medium, (2) shocked medium, (3) shocked jet, and (4) unshocked jet. Our initial conditions define the parameters in regions 1 and 4. We use the subscripts “4” and “1” in the following parameters to denote the unshocked jet and the unshocked medium, respectively, $\Gamma_{4}$ is the initial jet Lorentz factor, $n$ is the number density, and $P$ is the thermal pressure. The parameters are all in the code units (Table 1). One can convert them to the physical units by assigning $L_{0}$ and the typical magnetic field strength $B_{0}$ to certain values. For example, if one takes a typical physical parameters $L_{0}=3\times 10^{13}$ cm, one then has the time normalization parameter $t_{0}=L_{0}/c=10^{3}$ s in the frame, which corresponds to $t_{\rm obs}\sim t_{0}/\Gamma_{4}^{2}=10$ s in the observer frame. Our simulation lasts for $19t_{0}$, which corresponds to $\sim 190$ s in the observer frame, a time scale shorter than the deceleration time (or RS crossing time) for $\Gamma_{4}=10$. This is also consistent with the estimate that the simulation scale $20L_{0}=6\times 10^{14}$ cm is shorter than the thickness of the shell $R_{\rm dec}/\Gamma^{2}=(10^{15}~{}{\rm cm})R_{\rm dec,17}\Gamma_{4,1}^{-2}$, so that the simulation box is only part of the global jet.

For a conical jet at the RS crossing time (which is also the deceleration time), the shock jump conditions at both FS and RS gives a relation (Sari & Piran, 1995; Zhang & Kobayashi, 2005)222To derive this formula, the adiabatic index $\hat{\gamma}$ has been assumed to be $\Gamma$-dependent, i.e. $\hat{\gamma}=(4\Gamma+1)/3\Gamma$ (e.g. Kumar & Granot, 2003; Uhm, 2011).

[TABLE]

where $\Gamma_{21}$ and $\Gamma_{34}$ are the relative Lorentz factors between upstream and downstream at the FS and RS, respectively. For a thin shell ($\Gamma_{34}\gtrsim 1$), which is relevant for most GRB problems, one gets a rough relation $n_{4}/n_{3}\sim\Gamma_{4}^{2}$. Since our simulation does not introduce a conical jet (for which $n_{4}/n_{1}$ decreases with time), we set up $n_{4}/n_{1}=\Gamma_{4}^{2}$ by hand to mimic the physical condition at the deceleration time.

The initial setup of the number density distribution is shown in Fig. 1a. The high density ejecta ($n_{4}=100$) and the low density ambient ($n_{1}=1$) are separated at the position $z=1$. Some density clumps are introduced in the ambient medium to mimic density fluctuation. The density profile of each clump is delineated in the form $n=n_{1}+C\cdot e^{-r}$, where r is the distance away from the center of the clump, and the displacement between the centers of two neighboring clumps are $y_{disp}=0.04L_{0}$ and $z_{disp}=2L_{0}$, and $C=100$ is the amplitude of the fluctuation333The degree of distortion depends on the distribution and amplitude of the clumps. A denser distribution of the clumps tends to distort the RS more significantly, leading to a larger PD reduction for a certain initial $\sigma$ value. We choose a relatively sparse distribution of the clumps to set a more conservative “lower limit” of the relative PD reduction.. The initial magnetic field configuration is assumed to be uniform in the jet, with $B_{y,4}$ only in the co-moving frame. In the lab frame, an electric field component in the $x$ direction ($E_{x,4}$) with a similar amplitude as $B_{y,4}$ (self-consistently determined by $B_{y,4}$ and $\Gamma_{4}$) is introduced so that a Poynting flux is coupled with the jet in the lab frame. Different values of $B_{y,4}$ are explored, which correspond to different initial values of the magnetization factor $\sigma_{4}$ in the ejecta. Since the initial thermal energy is much smaller than the kinetic energy, $\sigma_{4}$ can be calculated as $\sigma_{4}\approx B_{y,4}^{2}/4\pi\Gamma_{4}^{2}nm_{p}c^{2}$, where $m_{p}$ is the proton mass and $c$ is the speed of light. The typical Lorentz factor ($\Gamma_{4}$) of GRBs at the deceleration phase is around several hundred (e.g. Liang et al. 2010). Such a large Lorentz factor is not easy to achieve for stable simulations. We set up two sets of simulations with $\Gamma_{4}=10,20$, respectively444The maximum $\Gamma$ that the code can handle is higher and varies depending on the values of other parameters, e.g. $\sigma$. (Fig. 1b). Although these values are smaller than the true Lorentz factor of GRBs at the deceleration phase, the simulation can nonetheless catch the key physics explored in this paper. The two different values $\Gamma_{4}$ adopted in our simulations can reflect the dependence of results on the initial Lorentz factor.

The parameters for all the simulation models are listed in Table 1. Model 0 is the reference model without density fluctuation and magnetic field in the jet, which is used to test the basic evolution process of the blastwave. The simulation results are fully consistent with analytical solutions, suggesting the validity of relativistic Athena++. Models 1-5 have the same density fluctuations but different $\sigma_{4}$ (and hence $B_{y,4}$) in the ejecta. This allows us to explore how magnetic field distortion and the resulting PD vary with $\sigma_{4}$. Models 6-10 have a different initial Lorentz factor ($\Gamma_{4}$=20) corresponding to Models 1-5 ($\Gamma_{4}$=10), which is a preliminary test of the effect of $\Gamma_{4}$.

The $\sigma_{4}$ values are set to be below unity based on two considerations. (1) Significant magnetic dissipation has likely occurred during the prompt emission phase (Zhang & Yan, 2011; Deng et al., 2015), so that the $\sigma$ of the outflow has significantly reduced beyond the prompt emission radius. There is a continuous decrease of $\sigma$ as the jet expands and at the deceleration radius, $\sigma_{4}$ is likely below unity (e.g. Gao & Zhang, 2015). (2) The RS is weak or does not exist if $\sigma_{4}$ is above unity (Zhang & Kobayashi, 2005; Mimica et al., 2009; Mizuno et al., 2009). Since observationally one does observe the RS emission, $\sigma_{4}$ cannot be too high.

3. results

3.1. The representative case: Model 3

We first present the detailed evolution of the representative case of Model 3 ($\sigma_{4}=10^{-3}$). The zoom-in contour plots of the number density (Fig. 2b) and the two components of the magnetic field, $B_{y}$ (Fig. 2c) and $B_{z}$ (Fig. 2d) are presented at $t=8t_{0}$, where $t_{0}=L_{0}/c$ is the time normalization factor between the code and physical units. As a comparison, we also show the density contour of Model 0 in Fig. 2a. One can see that the surface of the contact discontinuity (CD, the separation between the yellow and the red regions) is highly distorted in Model 3 due to the perturbation of the density clumps. Also, in the RS region, the initially pure ordered magnetic field $B_{y}$ is significantly distorted and a significant $B_{z}$ component is generated. This significant topological change of the magnetic field would inevitably lead to the reduction of the synchrotron radiation PD in the RS emission.

Figure 3 shows the synchrotron radiation PD evolution of Model 3. To calculate the polarization evolution, we use the same method used in our previous work (Deng et al., 2016; Zhang et al., 2016). We first define the region between the CD and the RS in the surface density contour as the region where RS emission comes from. In order to handle the highly distorted CD surface, we develop a numerical approach to define the complex CD profile based on a sharp density discontinuity in the numerical cells. We then uniformly inject particles with a non-evolving power-law energy distribution ($N_{e}=N_{e,0}\gamma_{e}^{-2.5}$, where $\gamma_{e}\in(10,10^{6})$ and $N_{e,0}$ is an arbitrary normalization factor555We do not simulate the light curve and spectral evolution in our radiation calculation. The simulated PD and PA evolutions do not depend on the particle number normalization factor ($N_{e,0}$), which can be taken as an arbitrary number in our polarization calculations.) in the entire RS region. Next, we use the “3DPol” multi-zone polarization-dependent ray-tracing radiation code (Zhang et al., 2014) to calculate the evolution of the PD from the RS region based on our MHD simulation results. The polarization calculations are carried out by post-processing the MHD simulation results with particles injected every two time steps in the MHD simulation. From Fig. 3, we show that the PD quickly (in $\sim 2t_{0}$) drops from the initial high value to a lower value. This is caused by the distortion of magnetic fields due to density fluctuations. In long term, the PD gradually recovers due to the resilience of the magnetic field lines. However, we believe that such a recovering effect may be artificial. By considering a more realistic global magnetic configuration (Deng et al., 2015) and the Rayleigh-Taylor instability in the conical shell (Duffell & MacFadyen, 2013), one would expect that the distorted $B$ configuration may not be easily recovered. In the following discussion, we focus on the minimum value in the PD curve and use it to define the reduction of PD from the original value.

3.2. Relationship between PD reduction and $\sigma$

We run other Models with different initial $\sigma_{4}$ for each group of $\Gamma_{4}$, and display the initial PD and final PD in Fig. 4a. For the uniform magnetic field distribution and a planar geometry, the maximum PD allowed for synchrotron radiation, i.e. $\Pi=(p+1)/(p+7/3)\simeq 72\%$ for $p=2.5$ (Rybicki & Lightman, 1979), can be achieved. After density perturbation, one can see that the minimum PD in the RS region clearly scales with $\sigma_{4}$ for each group of $\Gamma_{4}$. The smaller the $\sigma_{4}$ value, the lower the minimum PD. In reality, the field configuration after the prompt emission phase must have been distorted due to magnetic dissipation processes such as internal-collision-induced magnetic reconnection and turbulence (Zhang & Yan, 2011; Deng et al., 2015), the initial $\Pi$ in the RS crossing phase must be lower than 72%. As a result, the relative PD reduction, defined as $\xi=(\Pi_{0}-\Pi_{\rm min})/\Pi_{0}$, may be more relevant. This is plotted in Fig. 4b, which clearly shows a decrease of the relative reduction with increasing $\sigma_{4}$ for each group of $\Gamma_{4}$. Fits to this power-law regime ($\sigma_{4}=0.1$ removed) give empirical relations $\xi=6.0\sigma_{4}^{-0.21}$ for the group of $\Gamma_{4}=10$, and $\xi=1.0\sigma_{4}^{-0.37}$ for the group of $\Gamma_{4}=20$.

The behavior may be understood in terms of pressure and energy. One may consider the problem in the rest-frame of the RS region, in which the RS region is a target with growing thickness (as the RS front propagates into the jet) while the density clumps are bullets penetrating into the target. In our simulations, the ram pressure ($P_{r}$) and kinetic energy (${\cal E}_{K}$) of the density clumps remain the same for different models. The co-moving magnetic pressure ($P_{B}={B^{\prime}}^{2}/8\pi$) varies with $\sigma$ in different models, and the co-moving magnetic energy (${\cal E}_{B}$) varies with $\sigma$ and also with time as the volume of the RS region expands. In all our models, $P_{r}$ is designated to be greater than $P_{B}$, so that the magnetic field configuration can be distorted. When $\sigma$ is small, the magnetic field is dynamically unimportant. The magnetic amplification factor across the RS front depends on $\Gamma$ only when the initial density and pressure ratio are fixed, so that the magnetic pressure in the RS region would scale linearly with $\sigma$ in logarithmic space for each group of $\Gamma_{4}$. This is indeed seen from the numerical data in Fig. 4c. Since the magnetic pressure is the main stabilizer of the magnetic configuration against density-fluctuation perturbation, and the relative PD reduction therefore follows a power-law relation with $\sigma$ in the $\sigma\ll 1$ regime. We also find that the power-law indices are different between two groups with the different initial Lorentz factor, which means that the indices may be a function of the initial Lorentz factor, the initial density and pressure ratio, and the amplitude of the ambient density perturbation, which need additional investigation in the future.

When $\sigma_{4}$ is higher than 0.01, the magnetic field starts to affect the RS dynamics, which would decrease the RS strength and the compression ratio and in the meantime increase the reverse shock speed in the CD frame. Although the magnetic pressure increases more slowly (Fig. 4b) due to the weaker RS, the thickness of the RS region and hence the magnetic energy increase more rapidly due to the faster RS speed (Zhang & Kobayashi, 2005; Mimica et al., 2009; Mizuno et al., 2009). Even though the field configuration close to CD is distorted more, the field configuration from the far end of CD (close to the RS front) is less perturbed due to the larger volume of the RS region so that, on average, the field configuration in the entire RS region is less perturbed. This explains the roll-over in the $P_{B}-\sigma$ relation in Fig. 4c, and the relatively small PD reduction at higher $\sigma$ in Fig. 4b.

4. conclusions and discussion

In this paper, we explore how resilient a magnetic field is in a relativistic jet against density perturbation from the ambient medium in view of the recent observations of moderately high linear PD in some GRB early optical afterglows (Steele et al., 2009; Mundell et al., 2013). The motivation is to set a possible lower limit on the $\sigma$ value of the outflow based on the data. In view of the complication of the GRB RS problem, we limit ourselves to a cylindrical-like jet with a planar geometry and an initially uniform ordered magnetic field configuration. By introducing a perturbation due to density fluctuations from the ambient medium, we investigate how the distortion of magnetic configuration and reduction of PD depend on the $\sigma$ parameter of the jet. As expected, our numerical simulations using the relativistic Athena++ code show more significant reduction of PD with decreasing $\sigma$. The numerical results show an interesting power-law behavior between the relative PD reduction and $\sigma$, which may be applied to place constraints on the jet composition of GRBs. For example, if at the end of GRB prompt emission the PD is of $\sim(10\%-40\%)$ (e.g. Deng et al., 2016), the detection of $\sim 10\%$ PD in GRB 090102 (Steele et al., 2009) and $\geq 28\%$ in GRB 120308A (Mundell et al., 2013) at the early stage of the RS dominant phase666In this paper, we do not simulate the PD and PA data of GRB 120308A (Mundell et al., 2013), which requires a more detailed simulation of a conical jet interacting with an ambient medium with the emission from both the RS and FS represented. Such a simulation is beyond the scope of this paper. suggest that the PD is not reduced significantly in the RS region. Since turbulence is ubiquitous in a GRB environment (Zhang & Yan, 2011), the preservation of such a moderately high degree of PD suggests a relatively strong $B$ field in the RS region. The low-$\sigma$ cases (e.g. $\sigma_{4}<10^{-3}$) may be ruled out. This conclusion is consistent with light curve modeling of GRB early afterglows that shows a relatively large magnetic microphysics parameter ratios between FS and RS (Fan et al., 2002; Kumar & Panaitescu, 2003; Zhang et al., 2003; Mimica et al., 2010; Gao et al., 2015).

Other effects not included in the current simulations include the possible RT instability in a conical jet geometry, the initial more complicated magnetic configuration directly adopted from our previous simulations of prompt emission (Deng et al., 2015), and the 3D effect that in principle can lead to more tangled fields. We plan to investigate these effects in future more detailed studies. In any case, all of these parameters tend to further destroy the ordered magnetic field configurations and reduce the PD in the RS emission. The PD values obtained in this paper are conservative upper limits. The $\sigma$ constraints using the results obtained from this paper can be regarded as conservative lower limits.

We are grateful to Christopher White, Xuening Bai, Gustavo Rocha da Silva, Zhaohuan Zhu, Jiming Shi, Yanfei Jiang and Chang-Goo Kim for the important technical support during the usage of the “Athena++” MHD code, Haocheng Zhang for the helpful instruction regarding the “3DPol” polarization calculation code, and Fan Guo, Haocheng Zhang, Xiangrong Fu, and Shangfei Liu for helpful discussions and suggestions. An anonymous referee is greatly appreciated for very helpful comments. This work is supported by NASA through grants NNX15AK85G and NNX14AF85G, and by the LANL/LDRD program and Institutional Computing Programs at LANL and by DOE/Office of Fusion Energy Science through CMSO.