Weibel-Mediated Shocks Propagating into Inhomogeneous Electron-Positron Plasmas

TL;DR

This study uses PIC simulations to show that inhomogeneous upstream plasmas cause collisionless shocks to sustain stronger magnetic fields downstream, explaining observed GRB afterglow features.

Contribution

It demonstrates that upstream density inhomogeneities significantly influence magnetic field strength in collisionless shocks, a novel insight into shock physics.

Findings

Magnetic fields decay more slowly in inhomogeneous media.

Sound and entropy waves are generated by shock interactions.

Temperature anisotropy can account for observed magnetic field strengths.

Abstract

The external forward shock emitting the gamma-ray burst (GRB) afterglow is collisionless, and it is mediated by the Weibel instability which generates the magnetic field. The GRB afterglow shows that the magnetic field in the large downstream region is much stronger than the shock-compressed pre-shock field. However, particle-in-cell (PIC) simulations of relativistic shocks propagating into homogeneous media show that the Weibel generated field decays near the shock front. Some GRB observations and theoretical studies suggest that the preshock medium is inhomogeneous. We perform the PIC simulation of a relativistic shock propagating into inhomogeneous plasma. It is found that the post-shock magnetic field decays slowly compared with the homogeneous case. Sound waves and entropy waves are also generated by the shock-wave interaction, and temperature anisotropy is produced by the sound wave in the downstream region. The free energy of the temperature anisotropy is large enough to explain the observed field strength. Our results show that the upstream density fluctuation has a significant effect in the downstream region of collisionless shocks even if the wavelength of the upstream inhomogeneity is much larger than the kinetic scale.

Figures

Equations

$$\frac{c}{{\omega }_{\mathrm{pe}}}=\sqrt{\frac{\Gamma m_{\mathrm{e}}{c}^2}{4\pi n_0{e}^2}},$$

$$n(x,y)=n_0[1+0.5\sin (2\pi x/{\lambda }_x)].$$

$$\kappa \sim \frac{1}{3}{r}_{\mathrm{g}}c=\frac{1}{3\sqrt{8}}\frac{c^2}{{\omega }_{\mathrm{pe}}}{ϵ}_{\mathrm{B}}^{-1/2},$$

$$t_{\mathrm{decay}}$$

$$T_x$$

$$T_y$$

$${\gamma }_{\max }$$

$$t_{\mathrm{grow}}$$

$${ϵ}_{\mathrm{sw}}$$

$${\lambda }^{\mathrm{u}}<R_{\mathrm{dec}}^{\mathrm{u}}=1.2\times 10^{17}{\left(\frac{E_{\mathrm{iso},53}}{n_0^{\mathrm{u}}}\right)}^{1/3}{\Gamma }_{\mathrm{sh},2}^{\mathrm{u}-2/3}\mathrm{cm},$$

$$t_{\mathrm{lifetime}}=\frac{{\lambda }_{\mathrm{down}}}{⟨v⟩}+t_{\mathrm{dis}}.$$

$${\lambda }^{\mathrm{u}}>$$

Full text

††slugcomment: Submitted to ApJ: July 19, 2019

1. Introduction

The energy dissipation in relativistic shocks is important in the study of the production of relativistic particles (cosmic rays) and the photon emissions in high-energy astronomical transients. Such phenomena occur in dilute environments, so that the relativistic shocks arise through the kinetic processes in collisionless plasmas. The Weibel instability, which occurs in the collisonless plasmas with temperature anisotropy (Weibel, 1959), has been actively investigated in this field. It generates magnetic fields and would be important in the particle acceleration in the relativistic collisionless shocks (Medvedev & Loeb, 1999; Spitkovsky, 2008b; Sironi et al., 2013; Matsumoto et al., 2017). The Weibel instability also works for the excitation of non-relativistic shocks and especially is crucial role in the particle acceleration (Kato & Takabe, 2010; Matsumoto et al., 2015). Recently, the generation of the Weibel mediated shock has been also studied using high-power lasers (Fox et al., 2013; Huntington et al., 2015; Ross et al., 2017).

Observations of gamma-ray burst (GRB) afterglows have shown that they arise from the synchrotron radiation from energetic electrons in post-shock regions (Sari et al., 1998). Various theoretical modelings have suggested that the fraction of the magnetic-field energy density converted from kinetic energy density, $\epsilon_{\rm B}$, widely ranges from $\sim 10^{-8}$ to $\sim 10^{-2}$ with the mean of about $10^{-5}$, with an assumption that the density of $1\ {\rm cm}^{-3}$ (Santana et al., 2014). The lower limit is still larger than the shock-compressed value of interstellar magnetic field ($\sim 10^{-9}$), so that the field amplification is necessary. Moreover, since afterglows continue to be bright over a period of days to years, the field should remain strong in flowing downstream (Gruzinov & Waxman, 1999; Piran, 2005; Katz et al., 2007). In principle, the value of $\epsilon_{\rm B}$ may be constrained from the observed ratio of the peaks in the photon energy spectrum of the synchrotron and inverse Compton emissions. However, it still has large uncertainty because of the small number of events detected in GeV-TeV energy bands. Therefore, it is important to theoretically estimate the value of $\epsilon_{\rm B}$ through the micro-processes at the shock for understanding of the radiation mechanism of the GRB afterglows.

So far, the nonlinear evolution of the magnetic field generated by the Weibel instability has been investigated for shocks propagating into homogeneous plasmas (Silva et al., 2003; Kato, 2005; Chang et al., 2008; Keshet et al., 2009; Medvedev et al., 2011; Lemoine, 2015; Takamoto et al., 2018; Ruyer et al., 2018). However, it is currently unknown whether the Weibel generated field can survive long enough to explain the GRB afterglow emission. According to the two dimensional pair plasma shock simulation in Chang et al. (2008), the generated field rapidly decays, and the decay rate is $\epsilon_{\rm B}\approx(t\omega_{\rm pe})^{-1}$, where $\omega_{\rm pe}$ is the plasma frequency of electrons measured in the downstream rest frame. This means that the Weibel field can not persist during the GRB afterglow emission lasting for $t\omega_{\rm pp}\sim 10^{7}$, where $\omega_{\rm pp}$ is the plasma frequency of protons measured in the downstream rest frame (Sironi & Goodman, 2007). In the long-term simulation ($t\omega_{\rm pe}\sim 10^{4}$) performed by Keshet et al. (2009), the field is moderate in decaying because a wavelength of the magnetic field becomes large owing to the high-energy particles accelerated at the shock. The decay rate also depends on the time evolution of the spectrum of the field fluctuation (Chang et al., 2008; Lemoine, 2015). Chang et al. (2008) gave the analytic form of the decay rate as, $\epsilon_{\rm B}\propto(t\omega_{\rm pe})^{-2(p+1)/3}$, where p is the long-wavelength spectral index, and $p\simeq 1/2$ at late time ($t\omega_{\rm pe}\sim 10^{3}$).

Previous local PIC simulations with periodic boundary conditions have shown that the current filaments more easily merge in the three-dimensional case than in the two-dimensional case. The wavelength becomes longer and the decay rate decreases in both cases of the electron-positron (Silva et al., 2003) and electron-ion plasmas (Takamoto et al., 2018). Medvedev et al. (2011) suggested that the decay exponent is $-2$ in the pair plasma, although the simulation time ($t\omega_{\rm pe}\simeq 50$) is much shorter than the timescale of GRB afterglows. Recent three-dimensional simulation with larger simulation box in the flow direction has shown that the current filament in electron-ion plasma is broken by the drift-kink mode before being large scale (Ruyer et al., 2018). These simulations were, however, performed with the periodic boundary condition along the flow direction, which is different from the non-periodic shock structure, and no effect of particle acceleration is considered.

Meanwhile, the observed large dispersion of $\epsilon_{\rm B}$ appears to be a key ingredient in understanding the origin of the amplified magnetic field. If the field strength and plasma temperature in upstream region are similar to those in Galactic ISM, then both Alfvén and sonic Mach numbers are much larger than unity at early afterglows, where we assume the Lorentz factor of the flow $\sim 100$. Then, the downstream shock structure does not depend on these quantities, that is, they are not important parameters in calculating the GRB afterglow emission (Sari et al., 1998). More important is the upstream number density to determine the downstream magnetic field. Santana et al. (2014) showed that the amplification of the magnetic field has the weak dependence of the upstream density; ${\rm AF}\propto n_{\rm 0}^{0.21}$, where AF is an amplification factor characterizing how mach the seed magnetic field was amplified beyond the shock compression. Hence, one can see that the value of $\epsilon_{\rm B}$ is not determined by the upstream density alone, and there should be another piece controlling the downstream state.

The upstream medium of the shock wave is generally inhomogeneous. For example, a progenitor of the long GRB is thought to be a massive star. Then, the environment around them has some density fluctuations produced by the stellar wind or some inhomogeneity of interstellar medium (ISM). The characteristic length scale of the density fluctuations is uncertain; for example, it is on the order of 100 pc for ISM turbulence (Armstrong et al., 1995), or smaller for circumstellar matter (CSM) (Chugai & Danziger, 1994; Smith et al., 2009; Yalinewich & Zwart, 2019). There are some events suggesting the upstream density fluctuations (Wang & Loeb, 2000; Lazatti et al., 2002; Dai & Lu, 2002; Heyl & Perna, 2003; Schaefer et al., 2003; Nakar et al., 2003).

Then, the large dispersion of $\epsilon_{\rm B}$ may come from the inhomogeneity of the surrounding environment, that is, the wavelength of the upstream fluctuation potentially becomes an important parameter in describing the downstream magnetic field. In fluid scale, the field amplification process depends on the environment around a progenitor star. Magnetohydrodynamics (MHD) simulations have shown that the magnetic field is amplified by turbulent dynamo when the shock is propagating into the inhomogeneous media (Sironi & Goodman, 2007; Inoue et al., 2011; Mizuno et al., 2011). In this case, the decay timescale of the field is the eddy turn over time.

However, when the turbulent kinetic energy is dissipated, particle diffusion or kinetic plasma instability is excited in collisionless plasma. In addition, diffusive motion of high-energy particles may potentially impact on the shock structure. With fluid approximation, a linear analysis of an interaction between a shock and an incident entropy wave shows that entropy and subsonic modes are excited by the mode conversion in front of and behind the shock (Mckenzie & Westphal, 1968; Lemoine et al., 2016). Shu et al. (1989) and Arshed & Hoffman (2013) studied a shock tube problem in the case of the interaction between the shock and some waves. They claimed that an upstream incident entropy wave is compressed at the shock front, being advected stably in the downstream region. However, it is not clear whether the incident waves passing through the shock front can persist in kinetic scale when the particle diffusion is non-negligible. In this case, the density fluctuations may rapidly decay before the magnetic field is amplified by the turbulence in large scale. Moreover, the energy dissipation by the particle diffusion is expected to lead to the plasma heating and particle acceleration. The kinetic plasma instability may convert the turbulent kinetic energy to the magnetic energy. In this paper, using first principles approach, we investigate possible effects of the upstream density fluctuation on the evolution of downstream magnetic field.

In our previous study (Tomita & Ohira, 2016), two-dimensional local PIC simulation was performed to see the nonlinear evolution of the Weibel instability in electron-positron plasma with an anisotropic density structure, which mimics the downstream region of a relativistic shock propagating into inhomogeneous media. We showed that particles escape from an anisotropic high-density region, leading a temperature anisotropy large enough for the Weibel instability to make the magnetic field fluctuation.

This paper is organized as follows. In Section 2, we conduct two-dimensional PIC simulation of relativistic shocks propagating into the inhomogeneous media for the first time. Our simulation shows that the Weibel magnetic field has larger scale and higher magnitude than that for the case of the homogeneous density distribution. Such magnetic fields can persist in the far downstream region without significant decay. We also confirmed there are the entropy and sound waves in the downstream region during the simulation time. The temperature anisotropy is generated by the diffusion of high-energy particles, being maintained in the far downstream region. The field amplification by the Weibel instability is however, less significant in the downstream region (Section 3). It is finally claimed in Section 4 that the measured temperature anisotropy is large enough to explain the observed GRB afterglows. We also discuss the wavelength of the density fluctuation around the progenitor of GRBs (Section 4).

2. simulation setting

We use the two-dimensional electromagnetic PIC code, pCANS. In order to suppress the numerical Cherenkov instability (Godfrey, 1974), Maxwell’s equations are solved by the implicit method with the Courant-Friedrichs-Lewy number of 1.0 (Ikeya & Matsumoto, 2015). We set a two-dimensional simulation box in the $x$-$y$ plane. Periodic boundary condition is imposed in the $y$-direction. The simulation box size is $L_{x}=1.2\times 10^{4}c/\omega_{\rm pe}$ and $L_{y}=86.4c/\omega_{\rm pe}$, where

[TABLE]

is the plasma skin-depth of the upstream electron-positron plasma with a bulk Lorentz factor $\Gamma$ and upstream mean number density $n_{\rm 0}$. Constants $c,e,m_{\rm e}$ are speed of light, particle charge, particle mass, respectively. The cell size and time step are $\Delta x=\Delta y=0.1c/\omega_{\rm pe}$ and $\Delta t=0.1\omega_{\rm pe}^{-1}$, respectively.

Unmagnetized plasma with the bulk Lorentz factor of $\Gamma=10$ and the thermal velocity of $v_{\rm th}=0.1c$ is injected from the left boundary ($x=0$) and reflected at the right boundary ($x=1.2\times 10^{4}c/\omega_{\rm pe}$). The reflected and incoming particles interact, forming a shock which propagates toward the injection wall. The simulation frame is the downstream rest frame, that is, the downstream mean velocity is zero. We use 40 particles per cell per species. We investigated a dependence of a numerical heating on the number of particles by performing simulations for 10, 20, and 40 particles per cell per species. We confirmed that the more particles, the more suppressed are the numerical heating. We show simulation results for 40 particles per cell per species in this paper.

To investigate a collisionless shock propagating to an inhomogeneous plasma, the spatial distribution of the electron-positron plasma in the upstream region is set to be

[TABLE]

Namely, the density fluctuation is given by a monochromatic wave with $\delta n/n_{0}=0.5$. For simplicity, the upstream temperature is set to be constant, which makes the pressure fluctuation but it is negligible because the shock Mach number is very large. We performed two runs for $\lambda_{x}=4\times 10^{2}c/\omega_{\rm pe}$ and $1.2\times 10^{3}c/\omega_{\rm pe}$ in order to study a wavelength dependence. We confirmed that our results did not significantly depend on $\lambda_{x}$. In the following, we show simulation results for $\lambda_{x}=1.2\times 10^{3}c/\omega_{\rm pe}$ which is much larger than the shock thickness of the Weibel mediated shock, $\sim 50c/\omega_{\rm pe}$.

3. result

3.1. shock structure

First, we compare the shock structure for the case of homogeneous and inhomogeneous density distributions. The left column in Figure 1 shows the shock structures for the homogenous case at $t\omega_{\rm pe}=5000$. From top to bottom, the four panels show the two-dimensional structures of normalized electron-positron density $n/n_{0}$ and the energy fraction of the magnetic field $\epsilon_{\rm B}=\rm{B^{2}}/(32\pi\Gamma n_{0}m_{\rm e}c^{2})$, and one-dimensional profiles averaged over the $y$-direction of $n/n_{0}$ and $\epsilon_{\rm B}$, respectively. The shock front is located at $x\approx 9700~{}c/\omega_{\rm pe}$ and moving leftward with the velocity of $v_{\rm sh}\approx 0.37~{}c$. The left side of the shock front is the upstream region where magnetic filaments are generated by the Weibel instability. The transverse width of the Weibel filaments is about $8~{}c/\omega_{\rm pe}$ in the upstream region. The magnetic field generated by the Weibel instability decays from the shock front downstream. These features are consistent with early studies for relativistic Weibel mediated shocks (Spitkovsky, 2008b; Keshet et al., 2009). Compared with our previous periodic simulation (Tomita & Ohira, 2016) in which the Weibel instability is excited by counter-streaming plasmas, the Weibel-mediated field is maintained for a long time in the shock simulation. This is because the magnetic field with a longer wavelength can be excited in the shock simulation. According to Keshet et al. (2009), the large-scale magnetic field is generated by high-energy particles accelerated by the shock.

The middle and right panels in Figure 1 display the time sequence of the shock structures for the inhomogeneous case at $t\omega_{\rm pe}=$4800 and 5200, respectively. At $t\omega_{\rm pe}=4800$ (middle panels in Figure 1), the shock front is passing through the highest-density region. Although the upstream structures are quite similar to those for the homogeneous case, one can see significant differences in the downstream region. At $x\approx 10000\ c/\omega_{\rm pe}$, there are large-scale filaments and the magnetic-field strength is comparable to that in the shock front. They are generated as the shock front passes through the lower-density region. The right column in Figure 1 shows the shock structures at $t\omega_{\rm pe}=5200$ when the shock front is passing through the lowest-density region. The size of the upstream filaments is significantly larger than that for the homogeneous case. In addition, the magnetic field is widely generated in the broad upstream region (see panels (f) and (l)). The Weibel instability is driven by the cold upstream plasma and the hot plasma leaking from the downstream region. Since the density ratio of the hot leaking plasma to the cold upstream plasma becomes larger in the low-density region, the large momentum dispersion of the hot leaking plasma makes the transverse width of the filament wider (Ruyer et al., 2017). Since the length scale of the Weibel magnetic field is larger when the shock propagates in lower-density regions, it can be advected far downstream without dissipation. This is the reason why there is large-scale magnetic field in the downstream region at $t\omega_{\rm pe}=4800$ (panels (e) and (k)).

In Figure 2, we compare the spatial profile of the transversely averaged $\epsilon_{\rm B}$ in larger downstream region at $t\omega_{\rm pe}=6100$ for the homogeneous (a) and inhomogeneous (b) density distributions, respectively. The shock front is located at $x\approx 9200~{}c/\omega_{\rm pe}$ and passing through the lower-density region for the inhomogeneous case. For the homogeneous case, the magnetic field generated in the shock transition region simply decays downstream. For the inhomogeneous case, we can see two peaks of $\epsilon_{\rm B}$ downstream, which decay much slower than that for the homogeneous case. The number of the peaks is the same as that of the waves which transmitted through the shock. Therefore, the inhomogeneity of the upstream density has an important role on the magnetic field generation far downstream.

Next, we concentrate on the large-scale downstream structures of density and velocity for the inhomogeneous case. The blue and red curves in Figure 3 show the one-dimensional structures of $x$ component of the average velocity, $\langle v_{x}\rangle/c$, and density, $n/n_{0}$, at $t\omega_{\rm pe}=6100$ when the shock front is located at $x\approx 9200~{}c/\omega_{\rm pe}$. As seen in the blue curve, the average velocity oscillates with an amplitude of $\sim 0.2c$. The velocity profile (the blue curve) almost correlates with the large-scale density fluctuation. Therefore, the sinusoidal variation with the wavelength $\lambda_{\rm down}\approx 1100~{}c/\omega_{\rm pe}$ can be interpreted as a sound wave propagating downstream. The propagation velocity in the downstream rest frame is about ${\rm 0.62c}$ that is consistent with the sound velocity of relativistic plasmas, $c/\sqrt{3}$.

There is another intermediate-scale density fluctuation with a wavelength of about $400~{}c/\omega_{\rm pe}$ at $x\approx$ 9500, 10000, 10500, 10900 $c/\omega_{\rm pe}$ in Figure 3, which is not associated with the velocity displacement. They are non-propagating. Therefore, this density fluctuation can be interpreted as an entropy wave. The observed wavelengths of the entropy and sound waves are consistent with results from hydrodynamic analysis of interaction between a shock wave and a monochromatic density wave. This should never be taken for granted because we consider a collisionless system at finite temperature. If no particles change their momentum, density structures have to disappear in a sound crossing time. However, the simulation indicates that entropy and sound waves do not decay in the downstream region. The Weibel magnetic field changes the particle’s momentum and plays an important role in the momentum exchange of particles. The Larmor radius of thermal particles $r_{\rm g}$ is comparable to or slightly smaller than the coherence length of the Weibel magnetic field fluctuation, so that the diffusion coefficient of the thermal particles can be approximately given by

[TABLE]

Then, the decay time scales due to diffusion of the entropy and sound waves in the downstream region are given by

[TABLE]

The decay time scales of the sound and entropy waves are much longer than the simulation time, $t\omega_{\rm pe}=6100$ and sound crossing time, $t\omega_{\rm pe}=1935$. Therefore, the both waves can propagate in the downstream region during the simulation time.

3.2. temperature anisotropy

Figure 4 shows the spatial profile of temperature anisotropy averaged over the $y$-direction in the downstream region for the homogeneous (black solid) and inhomogeneous (red solid) cases at $t\omega_{\rm pe}=6100$. The blue dashed curve shows the $x$ component of the average velocity, $\langle v_{x}\rangle/c$. The temperature anisotropy is defined by $T_{y}/T_{x}-1$, where according to Yoon (2007) the temperatures $T_{x}$ and $T_{y}$ are given by

[TABLE]

Here, $\gamma,f,p_{x},p_{y},$ and $p_{z}$ are the particle Lorentz factor, the distribution function and momentum in the center-of-mass rest frame which is calculated by all particles in region of width $\sim 20c/\omega_{\rm pe}$. The temperature anisotropy is almost zero in the downstream region for the homogeneous case, while for the inhomogeneous upstream case it oscillates synchronously with the sound wave with values ranging from $-0.1$ to 0.05. It has been confirmed that the amplitude of the temperature anisotropy does not change even when the wavelength of upstream fluctuations becomes three times smaller than that in this simulation. As seen in Figure 4, the absolute value of the temperature anisotropy becomes maximum at $x\approx 10100\ c/\omega_{\rm pe}$ and $\approx 11000\ c/\omega_{\rm pe}$ where the divergence of the mean velocity is negative. The origin of the temperature anisotropy will be discussed later.

The temperature anisotropy can generate magnetic field by the Weibel instability in the downstream region. However, the observed temperature anisotropy oscillates downstream. We need to compare the growth time scale of the magnetic field and oscillation period of the temperature anisotropy in order to investigate whether the magnetic field is generated by the Weibel instability or not. The maximum growth rate of the relativistic Weibel instability is given by,

[TABLE]

where $T_{y}>T_{x}$ and $\gamma_{\rm max}\ll\omega_{\rm pe}$ are assumed (Yoon, 2007). For the temperature anisotropy of 0.05 observed in our simulation, the growth time scale is

[TABLE]

The oscillation period of the temperature anisotropy is that of the sound wave, $\lambda_{\rm down}/c_{\rm s}\approx 1.7\times 10^{3}\omega_{\rm pe}^{-1}$. Since $t_{\rm grow}<\lambda_{\rm down}/c_{\rm s}$, the Weibel-unstable mode should grow in the downstream anisotropic temperature plasma. However, we could not identify the magnetic field generated by the Weibel instability in the downstream region. This might be because equation (10) is given by the linear analysis of the Weibel instability in unmagnetized plasma, so that the growth rate is suppressed by the downstream magnetic field generated in the shock front. In addition, although the linear analysis of the Weibel instability also assumed the Maxwell-J$\ddot{u}$ttner distribution, the observed temperature anisotropy is produced by the non-thermal particles. This could also be the reason for the discrepancy in the growth rate. Even if the free energy of the temperature anisotropy is fully converted to the magnetic-field energy during $t_{\rm grow}\approx 20\omega_{\rm pe}^{-1}$, we cannot see it significantly in the downstream field of the strength keeping a few percent. If the magnetic field generated in the shock front sufficiently decays in the far downstream region, the observed temperature anisotropy could play an important role in the generation of the magnetic field in the far-downstream region.

We explain the origin of the downstream temperature anisotropy. As mentioned above, the particle diffusion in the sound wave is negligible for thermal particles. On the other hand, for high-energy particles accelerated by the shock, their gyroradii become larger than the coherent length scale of the magnetic-field fluctuations, so that the diffusion of high-energy particles becomes important. Therefore, the mean motion of the high-energy particles deviates from the sound wave motion. Figure 5 shows the cartoon illustrating the momentum distribution in $p_{x}-p_{y}$ plane at the regions A, B, and C, where mean velocities are $\langle v_{x}\rangle>0,=0$, and $<0$, respectively. The blue dashed circles are isotropic momentum distribution in the rest frame of regions A and C. If they are measured in the rest frame of region B, the momentum distributions are represented as the red circles by the Lorentz transformation. The color contrast (dark or light) represents the density contrast (high or low). Since the sound wave motion concentrates on particles in the region B, high-energy particles in the regions A and C are more concentrated to the region B. As a result, the momentum distribution in the region B has larger momentum dispersion along the $p_{x}$-direction, so that the temperature in the $x$-direction becomes larger than that in the $y$-direction. On the other hand, the regions A and C lose particles with large $p_{x}$, so that, the temperature in the $x$-direction becomes smaller than that in the $y$-direction. Therefore, the temperature anisotropy in the downstream region is generated by the sound wave motion and diffusion of high-energy particles.

4. Discussion

In this paper, we have investigated effects of upstream inhomogeneity on unmagnetized relativistic collisionless shocks by means of PIC simulations. In our simulation, an upstream density sinusoidaly changes with a wavelength of $1.2\times 10^{3}\ c/\omega_{\rm pe}$ in the downstream rest frame, which is much larger than the thickness of the shock transition region, the characteristic length scale of the Weibel instability, and the kinetic scale of plasmas. Although it is still much smaller than the injection scale of the interstellar turbulence ($\sim 1-100\ {\rm pc}$), there are some mechanisms that generate the small-scale density fluctuations in a precursor region of collisionless shocks (Drury & Falle, 1986; Ohira, 2013b, 2014, 2016). In addition, stellar winds contain small-scale density variation compared with the ISM turbulence. These small-scale density fluctuations could play an important role in collisionless shocks.

Our simulation shows that, compared with results for shocks into a uniform plasma, a larger-scale magnetic field is generated in the shock precursor region when the shock propagates into a low-density region in the inhomogeneous plasma. The density ratio of the hot leaking plasma to the upstream cold plasma becomes larger in the low-density region. Then, a finite temperature effect on the Weibel instability becomes significant in the low-density region, so that the characteristic length of the Weibel instability is thought to be larger. Hence, the larger-scale magnetic field hardly decays after they were advected downstream. The long-lived large-scale field is anticipated in GRB afterglows and could make a particle acceleration more efficient although this effect was not observed in this simulation because of a short simulation time. If we consider a longer length scale of upstream density fluctuations, the coherence length and the decay time in the downstream region of generated magnetic fields would be longer.

We found that the upstream density variation generates sound waves in the downstream region. The downstream sound waves are responsible for the temperature anisotropy which can generate the magnetic field by the Weibel instability. Therefore, the energy of the downstream sound waves is expected to be converted to the downstream magnetic-field energy. Our simulations have shown that the ratio of the sound-wave energy density $4n_{0}\Gamma m_{\rm e}\langle v\rangle^{2}/2$ to the kinetic energy density is

[TABLE]

where the downstream density is assumed to be $4n_{0}$ and $n_{0}$ is the upstream number density in the downstream rest frame. Furthermore, our simulation have shown that temperature anisotropy on the order of $10^{-2}$ is generated in the downstream region by the sound wave, which means that the ratio of the free energy of the temperature anisotropy to the downstream kinetic energy is $\epsilon_{\rm ani}\sim 10^{-2}$. Since observations of GRB afterglows suggest the energy fraction of magnetic field with $\epsilon_{\rm B}\sim 10^{-5}$ (Santana et al., 2014) which is smaller than $\epsilon_{\rm sw}$ and $\epsilon_{\rm ani}$, the energies of the downstream sound wave and temperature anisotropy are sufficient to generate the magnetic field in the afterglow emission region. Therefore, for shocks in inhomogeneous plasmas, the magnetic field generation takes place not only at the shock precursor region but also far downstream of the shock.

In order for the upstream density fluctuations to affect the GRB afterglow emission, the wavelength, $\lambda^{\rm u}$ measured in the upstream frame, must be smaller than the deceleration radius of GRB jet, $R_{\rm dec}^{\rm u}$, that is,

[TABLE]

where $E_{\rm iso,53},n_{0}^{\rm u}$ and $\Gamma_{\rm sh,2}^{\rm u}$ are the isotropic equivalent kinetic energy of the blast wave in unit of erg, upstream density in unit of ${\rm cm}^{-3}$ and the initial Lorentz factor of the shock measured in the upstream rest frame, respectively (e.g. Sironi & Goodman (2007)). We have adopted the usual notation $\mathcal{Q}_{n}\equiv\mathcal{Q}/10^{n}$. For a later phase of GRB afterglows, the upstream density fluctuations with longer wavelength generate downstream magnetic field.

In the hydrodynamics, a finite-amplitude sound wave eventually steepens to a shock-like structure, and the wave motion is thermalized by dissipation. The lifetime of a downstream sound wave with the wavelength of $\lambda_{\rm down}$ in the downstream rest frame is given by (Stein & Schwartz, 1972)

[TABLE]

The first term represents a steepening time in which the sound wave becomes a shock-like structure. The second term represents a dissipation time. If the dissipation is mainly due to the shock wave, the dissipation time is given by $t_{\rm dis}=\lambda_{\rm down}/\langle{v\rangle}(M-1)$, where $M$ is a Mach number (Stein & Schwartz, 1972). However, at present we have not understood the main dissipation mechanism and cannot estimate the dissipation time because we consider a collisionless system. In this discussion, we consider only the steepening time as a lower limit of the lifetime of the downstream sound wave. Observed GRB afterglows are already bright until the deceleration time, $t^{\rm u}_{\rm dec}=R^{\rm u}_{\rm dec}/c$, at which observers see the emission peak in the optical band. Hence, the field generation occurs until this epoch, so that $\Gamma^{\rm u}_{\rm sh}t_{\rm lifetime}/\sqrt{2}>t^{\rm u}_{\rm dec}$. Thus, the lower limit of the wavelength of the upstream density fluctuation is given as,

[TABLE]

where we use $\lambda^{\rm u}\sim\Gamma_{\rm sh}^{\rm u}\lambda_{\rm down}$ (Lemoine et al., 2016). Since we ignore the dissipation time in this discussion, the above lower limit actually should be smaller.

In our previous local PIC simulation with periodic boundary conditions (Tomita & Ohira, 2016), we investigated effects of particle diffusion from an entropy mode in the shock downstream region. We showed that if the spacial structure of the entropy mode is anisotropic, particles anisotropically escape from high-density regions. Then, the temperature anisotropy is generated, so that the magnetic field is generated by the Weibel instability. However, the anisotropic escape of particles from the entropy mode was not observed in this shock simulation. This is because the strength of the large-scale long-lived magnetic field is large enough to suppress the particle diffusion from the entropy mode. Instead of that, high-energy accelerated particles escape from high-density regions of the sound wave, which makes temperature anisotropy in the shock downstream region. The large-scale magnetic field generated in the shock precursor region would decay more rapidly in a realistic three-dimensional simulation. In that case, the downstream entropy and sound waves would also rapidly dissipate due to the particle diffusion, so that the temperature anisotropy would be generated by the anisotropic escape from high-density regions. Therefore, a realistic three-dimensional PIC simulation should be performed in the future even though it is very challenging.

In this paper, we study the case where the wave vector of the upstream density fluctuation is parallel to the shock normal. In reality, there are clumpy structures in the ISM or CSM. When a high-density clump passes through the shock front, not only compressible modes but also vortexes (incompressible modes) are excited in the downstream region. If the particle diffusion works, the escape of particles from the vortex can generate the temperature anisotropy as the sound wave makes it. Magnetic-field lines generated by the Weibel instability might be stretched to a large scale by the vorticity. Then, the magnetic field would be further amplified by the turbulent dynamo. In addition, if the particle diffusion works efficiently, the vortex rapidly decays and the turbulent dynamo becomes inefficient compered with ideal MHD case. However, if the shock-compressed background magnetic field is rapidly amplified by the turbulent dynamo before particles escape from the high-density clump, the diffusion is strongly suppressed, so that the small-scale magnetic field cannot be generated by the Weibel instability. A larger simulation box and a longer simulation time are required to understand the effect of vortexes on collisionless shocks, which will be addressed in future work.

The decay of turbulence means heating or acceleration of particles by the turbulence. Although we cannot see in this study a significant difference of the number of high-energy particles between the homogeneous and inhomogeneous cases, the second-order Fermi-acceleration by the downstream large-scale turbulence is expected in larger and longer simulations (Ohira, 2013a; Asano & Terasawa, 2015).

We considered only electron-positron plasma in this study. For the standard model of the GRB afterglow, the external shock wave propagates to an electron-ion plasma. When the shock is ultra-relativistic, the downstream electron and ion temperatures become almost the same (Spitkovsky, 2008a; Kumar et al., 2015). Then, the results for electron-positron plasmas would be applicable to GRB afterglows. However, the shock eventually becomes non-relativistic at a later phase of the GRB afterglow. Therefore, in reality, three-dimensional effects, the upstream magnetic fields, multidimensional density structures, nonrelativistic physics, and ion inertia could play important roles on the GRB afterglow and collisionless shocks propagating to inhomogeneous plasmas. More realistic simulations should be performed to understand these effects. Although it is very challenging for numerical simulations, we may be able to reveal it in the laboratory astrophysics (Takabe et al., 2008). Furthermore, we can apply their knowledge to supernova remnant shocks which are believed to be the cosmic-ray accelerator in our galaxy.

5. Summary

We have studied unmagnetiezed, relativistic collisionless shocks propagating into inhomogeneous plasmas. For the case of inhomogeneous upstream density distribution, large-scale magnetic field is generated by the Weibel instability in the shock precursor region. As the large-scale field is advected downstream, in contrast to the homogeneous case, they did not decay significantly. In this case, the sound and entropy waves can propagate downstream during our simulation time, which produces the downstream temperature anisotropy. We found that the upstream inhomogeneity has free energies large enough to explain the field strength inferred from observed GRB afterglows.

We thank M. Hoshino, Y. Matsumoto, T. Amano for useful comments. The software used in this work was in part developed in pCANS at Chiba University. Numerical computations were carried out on Cray X50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work is supported by Grant-in-Aid for JSPS Research Fellow from Japan Society for the Promotion of Science (JSPS), No. 17J03893 (ST), and in part by JSPS KAKENHI grant numbers 16K17702 (YO), 19H01893 (YO) and 18H01232(RY).