Persistence of precursor waves in two-dimensional relativistic shocks

TL;DR

This study uses two-dimensional particle-in-cell simulations to analyze the persistence and efficiency of precursor electromagnetic waves in relativistic shocks, revealing their potential role in electron acceleration.

Contribution

It demonstrates that precursor wave emission persists in 2D shocks even at low magnetizations, contrasting with 1D results, and explores implications for wakefield acceleration.

Findings

Precursor waves persist at low magnetization in 2D simulations.

Wave amplitudes are smaller than in 1D but still significant.

Electromagnetic waves disturb upstream medium, forming density structures.

Abstract

We investigated the efficiency of coherent upstream large-amplitude electromagnetic wave emission via synchrotron maser instability at relativistic magnetized shocks by using two-dimensional particle-in-cell simulations. We considered the purely perpendicular shock in an electron-positron plasma. The coherent wave emission efficiency was measured as a function of the magnetization parameter {\sigma}, which is defined by the ratio of the Poynting flux to the kinetic energy flux. The wave amplitude was systematically smaller than that observed in one-dimensional simulations. However, it continued to persist, even at a considerably low magnetization rate, where the Weibel instability dominated the shock transition. The emitted electromagnetic waves were sufficiently strong to disturb the upstream medium, and transverse filamentary density structures of substantial amplitude were produced. Based on this result, we discuss the possibility of the wakefield acceleration model for the production of non-thermal electrons in a relativistic magnetized ion-electron shock.

Equations

$${\omega }_{pe}=\sqrt{\frac{4\pi N_1{e}^2}{{\gamma }_1{m}_e}},$$

$${\omega }_{ce}=\frac{e{B}_1}{{\gamma }_1{m}_{e}c}.$$

$${\sigma }_e\equiv \frac{B_1^2}{4\pi {\gamma }_1{N}_1{m}_e{c}^2}=\frac{{\omega }_{ce}^2}{{\omega }_{pe}^2},$$

$$\delta B\frac{\lambda }{c}\lesssim \frac{B_1}{{\omega }_{ce}},$$

$${ϵ}_B\lesssim {\left(\frac{\lambda }{c/{\omega }_{pe}}\right)}^{-2},$$

$${\int }_1^{\infty }{f}_e(\gamma )\mathrm{d}\gamma =1.$$

$$f(\gamma )\mathrm{d}\gamma \propto \gamma \exp (-\frac{\gamma m{c}^2}{kT}).$$

$$\frac{\omega }{{\omega }_{pe}}\sim 3.$$

$$a={\gamma }_1\sqrt{{\sigma }_e}\frac{{\omega }_{pe}}{\omega }\frac{\delta E}{B_1}\simeq {\gamma }_1\sqrt{{\sigma }_e}\frac{{\omega }_{pe}}{\omega }\frac{\delta B}{B_1}\sim \frac{\sqrt{2}}{3}{\gamma }_1\sqrt{{ϵ}_{conv}},$$

$$n^{\prime 2}=1-\frac{2{\omega }_{pe}^2}{{\omega }^{\prime 2}-{\sigma }_e{\omega }_{pe}^2},$$

$$\frac{{\omega }^2}{{\omega }_{pe}^2}=\frac{2}{1-n^2}+\frac{{\sigma }_e}{{\gamma }_1^2(1+{\beta }_{1}n\cos \theta {)}^2},$$

$${\omega }^2\simeq 2{\omega }_{pe}^2+k^2{c}^2.$$

$$v_g\equiv \frac{d\omega }{dk}=\frac{k{c}^2}{\sqrt{2{\omega }_{pe}^2+k^2{c}^2}}.$$

$$k_x={\beta }_{shock}{\gamma }_{shock}\sqrt{k_y^2+\frac{2{\omega }_{pe}^2}{c^2}}.$$

$${\omega }^4-(2{\omega }_{pe}^2+c^2{k}^2){\omega }^2-{\omega }_{pe}^2{c}^2{k}^2=0$$

$${\Gamma }_{max}\sim {\omega }_{pe},$$

Full text

1 Introduction

The origin of cosmic rays (CRs) has been a long-standing problem in astrophysics. CRs with energies up to $10^{15.5}$ eV are commonly believed to be generated by diffusive shock acceleration (DSA) (Bell, 1978; Blandford & Ostriker, 1978; Drury, 1983; Blandford & Eichler, 1987) at supernova remnant shocks in our galaxy. Observations of broadband non-thermal emission from young supernova remnants support this paradigm (e.g., Koyama et al., 1995; Bamba et al., 2006; Helder et al., 2009). In contrast, the origin of high-energy CR population ($>10^{15.5}$ eV), presumably of extragalactic origin, is not well understood. Active galactic nuclei (AGN) and gamma-ray bursts (GRBs) are among the possible sources of such ultra-high-energy CRs (UHECRs) (e.g., Hillas, 1984; Milgrom & Usov, 1995; Vietri, 1995; Waxman, 1995). Observations show that AGN produce high-energy photons via synchrotron and inverse Compton emission (e.g., Abdo et al., 2010a, b; Nolan et al., 2012), which clearly shows the presence of ultrarelativistic electrons. AGN are typically associated with relativistic jets with Lorentz factors of a few tens (e.g., Lister et al., 2016). The energetic electrons are often assumed to be produced at shockwaves associated with the jet (e.g., Marscher et al., 2008). The intense radiation from GRBs is also thought to originate from a relativistic outflow whose Lorentz factor may exceed several hundreds (e.g., Pe’er et al., 2007; Ackermann et al., 2014). The standard model for GRBs assumes the presence of shockwaves wherein charged particles are accelerated and produce photons through synchrotron radiation (see, e.g., Piran, 2005).

DSA has also been one of the most commonly proposed mechanisms for producing UHECRs in these sources. However, it is well known that DSA becomes less efficient at a relativistic shock propagating in a magnetized plasma (Begelman & Kirk, 1990). The Lorentz transformation of the magnetic field from the upstream rest frame into the shock rest frame boosts only the transverse component by the Lorentz factor of the shock. Therefore, a highly-relativistic shock is almost always superluminal, indicating that particles moving along the magnetic field line cannot diffuse back into the upstream. In such a shock, particles must diffuse across the magnetic field to be accelerated by the DSA mechanism. Since perpendicular diffusion is typically much slower than parallel diffusion, this may not be a plausible solution for particle acceleration in highly-relativistic shocks (see, e.g., Takamoto & Kirk, 2015).

Chen et al. (2002) proposed an alternative model for UHECR acceleration at GRBs. The model invokes a large-amplitude Alfvén wave in a relativistic plasma that generates an electrostatic wave behind it. A particle may be accelerated by this electric field to high energies, possibly to the UHECR-energy range. This is essentially an application of wakefield acceleration (WFA), as discussed in the context of laser-plasma interactions in the laboratory (Tajima & Dawson, 1979). WFA in laboratory plasmas occurs when an intense laser pulse (or transverse electromagnetic wave) propagates in a plasma. A Langmuir wave is excited via Raman scattering, in which the ponderomotive force exerted by the laser pulse expels electrons from the region of high laser intensity, whereas ions are nearly unaffected. Consequently, a large-amplitude charge density fluctuation is generated, which is associated with a longitudinal electric field called a wakefield. the phase velocity of the excited Langmuir wave is nearly equal to the group velocity of the laser pulse, and is close to the speed of light . Therefore, the wakefield is able to accelerate particles to highly-relativistic energies via the Landau resonance.

WFA may also be applicable to relativistic magnetized shocks in nature. It is known that large-amplitude electromagnetic waves are excited at relativistic shock fronts by synchrotron maser instability (SMI) driven by particles reflected off the shock-compressed magnetic field (Hoshino & Arons, 1991). The instability results from the resonance between the relativistic particle cyclotron motion and an electromagnetic wave of the extraordinary mode (X-mode). The fluctuating magnetic field of the X-mode wave is parallel to the ambient magnetic field and perpendicular to the wavenumber vector. In a pair plasma, the fluctuating electric field is perpendicular to both the fluctuating magnetic field and the wavenumber vector and thus the X-mode wave is linearly polarized. The electromagnetic waves may be emitted toward both the upstream and downstream directions; however, only the waves with upstream-directed group velocities greater than the shock propagating speed can escape upstream. Consequently, the appearance of such high-frequency precursor waves ahead of the shock is a common feature of relativistic magnetized shocks reproduced by one-dimensional (1D) particle-in-cell (PIC) simulations (e.g., Langdon et al., 1988; Gallant et al., 1992; Hoshino et al., 1992; Amato & Arons, 2006). The precursors emitted in 1D PIC simulations of relativistic magnetized shocks in electron-ion plasmas were sufficiently strong to generate a wakefield of substantial amplitude wherein electrons were accelerated (Lyubarsky, 2006; Hoshino, 2008). Hoshino (2008) demonstrated the generation of non-thermal electrons by the WFA and estimated the maximum attainable energy $\varepsilon_{max}$ as $\varepsilon_{max}/\gamma_{1}m_{e}c^{2}\sim L_{sys}/(c/\omega_{pe})$, where $\gamma_{1}$ is the upstream bulk Lorentz factor and $\omega_{pe}$ is the proper electron plasma frequency. According to this estimate, the maximum energy is proportional to the system size $L_{sys}$ of an astronomical object. Hoshino (2008) found that the maximum energy can exceed the theoretical limit of conventional shock acceleration models such as DSA and that the WFA plays an important role in relativistic magnetized shocks, especially at a low magnetization rate.

WFA in the context of relativistic shocks has so far been discussed solely through 1D simulations. It is not well known whether the same mechanism can operate in more realistic multidimensional systems. In general, inhomogeneity (or waves) may appear in the transverse direction of the shock and the waves emitted from different positions at the shock may not necessarily be coherent in phase. Because the precursor waves are a superposition of such waves, the efficiency of the WFA may deteriorate in an incoherent precursor. Another possible problem is the competition between the SMI and Weibel instability (WI) (Weibel, 1959; Fried, 1959). The WI is driven unstable by effective temperature anisotropy generated by the reflected particles in the shock-transition region (e.g., Kato, 2007; Chang et al., 2008). The mode is unstable for wavenumbers perpendicular to the shock normal, thereby appearing only in multidimensional simulations. The maximum growth rate of the WI including relativistic effects scales as the proper electron plasma frequency $\omega_{pe}$ for sufficiently strong anisotropy (see, e.g., Yoon & Davidson, 1987; Yang et al., 1993; Achterberg et al., 2007; Schaefer-Rolffs & Tautz, 2008). In contrast, the growth rate of the SMI is on the order of the relativistic electron cyclotron frequency $\omega_{ce}$ (Hoshino & Arons, 1991). As both instabilities are excited from the same free-energy source, the precursor wave emission efficiency may be affected or even completely shut off in a low-magnetization regime ($\omega_{ce}\ll\omega_{pe}$), where the WI grows more quickly than the SMI. Previous PIC simulation studies in multiple dimensions indeed showed that the shock transition is dominated by the WI at low magnetization $\sigma\lesssim 10^{-2}$ (e.g., Spitkovsky, 2005; Kato & Takabe, 2010; Sironi et al., 2013), where $\sigma$ is the ratio of the Poynting flux to the kinetic energy flux.

In fact, some earlier studies of two-dimensional (2D) magnetized shocks gave negative results for the WFA. Sironi & Spitkovsky (2009) reported that the precursor waves were seen only in the initial phase and had little effect in their long-term calculations. Later, Sironi & Spitkovsky (2011) found that although a wakefield was induced in their 2D simulations of magnetized shocks in an ion-electron plasma, the amplitude was not sufficiently strong to produce non-thermal particles, which is in clear contrast to the 1D results reported by Hoshino (2008). However, according to our numerical convergence study (see Appendix A), the numerical resolutions used in the earlier studies were insufficient to capture the precursor waves. In this study, we quantify the efficiency of precursor wave emission using high-resolution 2D PIC simulations to investigate the applicability of the WFA model to astrophysical relativistic magnetized shocks. For the purpose of estimating the precursor wave emission efficiency, we consider only shocks in electron-positron plasmas. We note that, although the WFA needs a finite inertial difference between the positive and negative charges, the emission mechanism itself is identical between the pair and ion-electron plasmas. Therefore, the precursor wave emission efficiency measured in a pair plasma shock will also give a good estimate for an ion-electron plasma.

This work is organized as follows. First, we describe our simulation setup in $\S$2. $\S$3 describes global structures of relativistic magnetized shocks in the presence of large-amplitude electromagnetic precursor waves. In $\S$4, the properties of precursor waves are comprehensively analyzed. In $\S$5, we discuss the applicability of the WFA model to relativistic magnetized shocks in an ion-electron plasma. Finally, $\S$6 summarizes this study.

2 Simulation setup

Our simulations of relativistic shocks in magnetized pair plasmas were performed in 2D systems in the $x$-$y$ plane using a fully relativistic electromagnetic PIC code (Matsumoto et al., 2013, 2015). The basic configuration of our simulations is illustrated in Figure 1. A cold pair plasma flow (with zero thermal spread) is continuously injected from the right-hand boundary of the simulation domain toward the $-x$ direction at a bulk Lorentz factor of $\gamma_{1}=40$. The particles are specularly reflected at the left-hand boundary, and the conducting-wall boundary condition is applied for the electromagnetic field. The boundary condition in the $y$ direction is periodic for both the field and the particles. The particles reflected off the wall interact with the incoming plasma, and a shockwave propagating in the $+x$ direction is formed. Consequently, the simulation frame corresponds to the downstream rest frame.

In this study, we solely focus on the purely perpendicular shock. In 2D simulations, there is a degree of freedom in the choice of direction of the ambient magnetic field. In general, due to the neglect of the third coordinate ($z$ in our coordinate system), the rotation of the ambient magnetic field around the shock normal may change the shock dissipation physics. Nevertheless, we only present simulation results with the ambient magnetic field pointing in the out-of-plane direction ($B_{z}$ in our coordinate system) in this study for simplifying our discussion. Note that our conclusion is not limited to this particular configuration. The results with an in-plane magnetic field configuration ($B_{y}$) will be published elsewhere.

The length and time were measured in units of the electron inertial length $c/\omega_{pe}$ and the inverse proper electron plasma frequency $\omega_{pe}^{-1}$, respectively. The proper electron plasma frequency is defined as follows:

[TABLE]

where $N_{1}$ is the upstream electron number density measured in the simulation frame. The number of particles/cell is $N_{1}\Delta x^{2}=64$ for both electrons and positrons upstream, where $\Delta x=\Delta y$ is the spatial grid size. The grid size is fixed to $\Delta x/(c/\omega_{pe})=1/40$ and the number of grids in each direction is $N_{x}\times N_{y}=20,000\times 1,680$ throughout this study. The relativistic electron cyclotron frequency $\omega_{ce}$ is also the characteristic timescale of the relativistic magnetized shocks, and it is defined as follows:

[TABLE]

It is known that in a PIC simulation involving a relativistic plasma flow, there appears a so-called numerical Cherenkov instability (Godfrey, 1974), which will grow to substantial amplitudes and cause non-negligible amounts of plasma heating, even in homogeneous systems with relativistic bulk flows. As the unphysical heating in the cold upstream flow will change the physics of the shock, the simulation box size and integration time are severely limited by the growth timescale of the numerical instability. To minimize numerical artifacts, the method developed by Ikeya & Matsumoto (2015) is applied to our code. This method is a variant of the method originally reported by Godfrey & Vay (2013) and Xu et al. (2013), who found that a particular choice of CFL number can dramatically reduce the numerical instability growth rate. However, the magic CFL number depends on the numerical scheme used for solving Maxwell’s equations. Ikeya & Matsumoto (2015) performed numerical experiments and found that the magic CFL number is equal to unity for an implicit Maxwell solver employed in the code. As we used a fixed grid size of $\Delta x/(c/\omega_{pe})=1/40$, the time step was automatically determined as $\omega_{pe}\Delta t=1/40$. In addition, we used a moving injector for the upstream boundary. The particles are continuously injected from an injector that moves away from the shock (i.e., toward $+x$ direction) with the speed of light. This minimizes the propagation time for the undisturbed upstream uniform plasma, enabling us to follow the long-term evolution of the shock without being affected by the numerical instability.

Note that one has to carefully check numerical artifacts in the code for an accurate estimate of the precursor wave emission efficiency. Since the precursor waves are high-frequency (and thus short-wavelength) electromagnetic waves, they may easily be damped (e.g., in the case of the insufficient resolution). In addition, the application of digital filtering often used to suppress the numerical Cherenkov noise may underestimate the emission efficiency. Therefore, we performed 1D simulations and investigated numerical convergence with respect to both the grid size and the number of particles (see Appendix A). The numerical resolution employed in the 2D simulations is motivated by the numerical convergence study.

We investigated the dependence on the magnetization parameter $\sigma_{e}$ defined as the ratio of the Poynting flux to the kinetic energy flux:

[TABLE]

which is the main controlling parameter of a relativistic magnetized shock. Our simulations were performed for the following seven cases: $\sigma_{e}=1$, $3\times 10^{-1}$, $1\times 10^{-1}$, $3\times 10^{-2}$, $1\times 10^{-2}$, $3\times 10^{-3}$ and $1\times 10^{-3}$.

3 Global Shock Structure

3.1 High-$\sigma_{e}$ Case

First, an overview of the global shock structure for a relatively high-$\sigma_{e}$ case is discussed. Figure 2 illustrates the global shock structure at $\omega_{pe}t=500$ for $\sigma_{e}=3\times 10^{-1}$. Shown from top to bottom are the electron number density $N_{e}$, the transversely averaged electron number density $\langle N_{e}\rangle$, the out-of-plane magnetic field $B_{z}$, the 1D profile for $B_{z}$ along $y/(c/\omega_{pe})=21$, the transversely averaged electrostatic field $\langle E_{x}\rangle$, and the phase-space plots of electrons in the $x-u_{xe}$ and $x-u_{ye}$ planes. All quantities are normalized by the corresponding upstream values. At this time, the shock front is clearly seen at $x/(c/\omega_{pe})\sim 280$. It has already propagated sufficiently far away from the left-hand boundary and a well-developed shock structure is formed. A shock formation may be identified by the density increase by a factor of three, as expected from the MHD Rankine-Hugoniot relation (see Appendix B for finite deviation from the theoretical prediction).

In the upstream region, large-amplitude fluctuations in $B_{z}$ are clearly seen. These are the electromagnetic waves emitted and propagated away from the shock front. The wave magnetic field is polarized in the $z$ direction (i.e., parallel to the ambient magnetic field direction). This is the signature of the X-mode wave and is consistent with the linear theory of the SMI (Hoshino & Arons, 1991) as well as previous 1D PIC simulation results (Langdon et al., 1988; Gallant et al., 1992; Hoshino et al., 1992; Amato & Arons, 2006). At around the tip of the precursor wave region ($450\leq x/(c/\omega_{pe})\leq 480$), the wavefront is roughly uniform in the $y$ direction and is obviously different from the other part of the precursor. The waves in this region are generated in the early phase of simulations and are contaminated by the initial and boundary conditions. Therefore, we excluded this region from our analysis presented below.

In the precursor region, transverse filamentary structures are clearly identified in both the density and magnetic fields. We think that these filaments are generated by a filamentation instability discussed by Kaw et al. (1973) and Drake et al. (1974) for an ion-electron plasma. This is a kind of parametric instability driven by the ponderomotive force of a large-amplitude electromagnetic wave. Strictly speaking, to the authors’ knowledge, this instability has not been analyzed for a pair plasma. However, the similarity with the ion-electron plasma case indicates that it is probably related to the filamentation instability. More comprehensive studies of this instability will be presented in a future publication. In any case, the appearance of filamentary structures in the precursor region is the strong evidence that the precursor waves still retain coherence in 2D systems.

As is seen in the phase-space density plot, the precursor waves cause strong heating in $u_{x}$ in the upstream region. Note that this is not clearly seen in $u_{y}$ because the Lorentz transformation (from the upstream rest frame to the laboratory frame) increases the thermal spread only in the $u_{x}$ direction. Because the heating region coincides with the density filament region, this may be attributed to precursor waves. As discussed in $\S$4.2, this apparent heating is merely a result of coherent oscillations in velocity in the strong wave electromagnetic field and cannot be referred to as a true heating.

The density filaments in the precursor are convected by the upstream flow and eventually hit the shock. Consequently, the shock surface is strongly modulated and large amplitude fluctuations both in density and magnetic field are seen in the downstream region. The downstream turbulence may be explained by the Richtmyer-Meshkov instability (Richtmyer, 1960; Meshkov, 1969) It is worth noting that the precursor wave activity survives even in the presence of such a strong feedback effect of the emission to the shock. As discussed in $\S$4.1, the wave emission efficiency has nearly reached a quasi-stationary state at this stage. This leads to the firm conclusion that coherent precursor waves persist in 2D systems.

3.2 Low-$\sigma_{e}$ Case

Further, we discuss the global shock structure for a relatively low $\sigma_{e}$ case. Figure 3 illustrates the global shock structure at $\omega_{pe}t=500$ for $\sigma_{e}=3\times 10^{-3}$. The format is the same as Figure 2. As we have already mentioned, our main concern in a low magnetization regime is the competition between the SMI and the WI. As expected, we observe filamentary density and magnetic field structures at $x/(c/\omega_{pe})=240-280$ in Figure 3, which are attributed to the WI (see Appendix C). Because the WI results from velocity anisotropy near the shock front, the length of the magnetic field filament corresponds to that of the reflected particle beam (or simply the particle Larmor radius for a perpendicular shock). It is clear that the density filaments are extended well ahead of the Weibel-dominated shock transition region where no back-streaming particles exist. We believe that these filaments are generated by the interaction between coherent radiation and upstream plasmas, just as in the high-$\sigma_{e}$ case discussed in $\S$3.1 .

Surprisingly, precursor waves are observed with appreciable amplitudes in this case as well. Although the amplitude of the density filaments is smaller than that in the high-$\sigma_{e}$ case, presumably due to smaller precursor wave amplitudes (see $\S$4.4), the formation of filaments indicates that the precursors strongly interact with the upstream plasma. Particle heating in the precursor region is not observed in Figure 3 again because of lower wave amplitude. Nevertheless, finite heating due to quiver motion in the wave electromagnetic field is indeed identified (see $\S$4.2). All these results suggest that the precursor waves remain coherent even in this case, and that the SMI and the WI may somehow coexist at a relativistic magnetized shock.

4 Precursor wave properties

4.1 Time Evolution

Here, we discuss the time evolution of precursor wave power. Figure 4 shows the time evolution of the wave energy averaged over the $y$ direction in the range $300\leq\omega_{pe}t\leq 500$. The wave energy is normalized by the upstream bulk kinetic energy. Two different runs with $\sigma_{e}=3\times 10^{-1}$ (left) and $\sigma_{e}=3\times 10^{-3}$ (right) are shown.

This study primarily focuses on a sufficiently long-term evolution wherein the effects of the initial and boundary conditions may be neglected. For this purpose, we focus on precursor waves in the immediate upstream of the shock front, as they are newly generated waves at the shock in the highly-disturbed medium. Figure 4 indicates that the precursor wave amplitude near the shock front gradually decreases over time. However, it remains finite and reaches a quasi-steady state by the end of the simulations.

We quantified the precursor wave power by taking the average power in the range: $r_{L}<x-X_{sh}<r_{L}+50c/\omega_{pe}$, where $X_{sh}$ is the position of the shock and $r_{L}=c/\omega_{ce}$ is the relativistic Larmor radius calculated with the upstream flow Lorentz factor $\gamma_{1}$. The shock position $X_{sh}$ is determined from the $y$-averaged density profile $\langle N_{e}\rangle$, assuming that the shock propagation speed is constant over time. Figure 5 shows the precursor wave energy density as a function of time in the range $300\leq\omega_{pe}t\leq 500$ for $\sigma_{e}=3\times 10^{-1}$ (left) and $3\times 10^{-3}$ (right). Note that the results with low resolution $\Delta x/(c/\omega_{pe})=1/20$ are also shown for comparison. Each data point represents the temporal average during the time interval $20/\omega_{pe}$, and the error bars were estimated from the standard deviation during the time interval. The error bars for $\sigma_{e}=3\times 10^{-1}$ were systematically larger than those for $\sigma_{e}=3\times 10^{-3}$ because of low-frequency fluctuations in $B_{z}$ associated with the density filaments. Nevertheless, the results suggest that the precursor wave power is nearly saturated at around $\omega_{pe}t\simeq 450$ in both cases at our fiducial resolution.

On the other hand, the lower resolution runs showed continuous decrease (for $\sigma_{e}=3\times 10^{-3}$) or saturation at a lower level (for $\sigma_{e}=3\times 10^{-1}$). Therefore, it is natural that even lower resolution simulations reported previously could not accurately reproduce the precursor waves (Sironi & Spitkovsky, 2009, 2011).

4.2 Precursor Heating

As shown in Figure 2, the $x$ component of the electron four-velocity $u_{xe}$ in the precursor wave region for $\sigma_{e}=3\times 10^{-1}$ has a substantial spread. However, we found that particles’ velocities are merely oscillating in the precursor wave electromagnetic field and that the apparent broadening in the velocity distribution does not indicate true heating.

To confirm this idea, Figures 6 and 7 show individual electron four-velocities (blue circle) and the Lorentz force (red solid line) calculated at each grid point along $y/(c/\omega_{pe})=21$ for $\sigma_{e}=3\times 10^{-1}$ and $\sigma_{e}=3\times 10^{-3}$, respectively. Note that $E_{x}$ is weak because of the perfect symmetry in a pair plasma. The top and bottom panels show the $x$ and $y$ components, respectively, at $\omega_{pe}t=500$. In addition, only  1% of the electrons located within one grid cell in the $y$ direction are shown.
The electrons obey these equations quite well, suggesting that the particles in the upstream are not actually thermalized but forced to oscillate in a coherent manner by the wave electromagnetic field.

It is easy to confirm that the electromagnetic waves propagating in a plasma always have superluminal phase speeds. Therefore, they are not susceptible to resonant wave-particle interactions and will not directly thermalize the particles. Heating may occur when a parametric instability produces longitudinal waves with subluminal phase speeds, allowing energy to be absorbed by the particles via wave-particle resonances. The WFA is indeed one of the examples of such heating. However, in the present simulations, although filamentation instability occurs, the electric field associated with the density filaments is weak due to the symmetry in a pair plasma. This is probably the reason why there is no true heating in the precursor even in the presence of filamentation instability.

In any case, the fast particle quiver motion may behave as an effective temperature that potentially reduces the growth rate of the SMI. Note that the SMI results from harmonic resonances between the relativistic particle gyromotion and the X-mode wave $\omega=n\omega_{ce}$, where integer $n$ denotes the harmonic number of the resonance. Amato & Arons (2006) reported that a finite thermal spread in the ring distribution suppresses the growth rate at higher harmonics, and the suppression becomes increasingly significant at lower harmonic numbers as the thermal spread of the particles increases. Therefore, we think that the effective temperature may be responsible for the gradual decrease of the precursor wave amplitude. The quasi-steady state in the final phase of the simulations is probably determined by the balance between wave generation and heating.

4.3 Wavenumber Spectra

Figure 8 shows the precursor wave power spectra in wavenumber space normalized by the upstream ambient magnetic field energy density. The spectra were obtained for a snapshot at $\omega_{pe}t=500$ by taking the Fourier transform in the region: $r_{L}<x-X_{sh}<r_{L}+50c/\omega_{pe}$. We applied the Hanning window to remove edge effects. Note that the Nyquist wavenumber for our simulation is $k_{N}=\pi/\Delta x\sim 120\omega_{pe}/c$, and the precursor waves are well resolved.

Each of these spectra has a lower cutoff wavenumber. This can be explained by the wavenumber above which the upstream-directed group velocity of the precursor wave is greater than the shock propagation velocity (Gallant et al., 1992). We estimated the theoretical cutoff wavenumber (see Appendix B), which is shown by the white solid line in Figure 8. The theoretical lower cutoff wavenumber is roughly consistent with the measured one. The wave power to the right of the white solid line may be attributed to the precursor waves propagating away from the shock front. Note that both spectra have non-negligible wave power around $k_{x}=0$, and in particular the power for $\sigma_{e}=3\times 10^{-1}$ has a substantial fraction of the total wave power. Recall that we excluded the region contaminated by the WI from this analysis. Therefore, we think that waves around $k_{x}=0$ are associated with the filamentation instability. These filamentary structures generated in the upstream region are convected by the upstream plasma flow, which then interact with the shock.

4.4 $\sigma_{e}$ Dependence

Now we discuss the $\sigma_{e}$ dependence of the precursor wave amplitude. The amplitude was calculated by integrating the power spectra as described above over the whole wavenumber space. Figure 9 shows the precursor wave power as a function of $\sigma_{e}$ with two different normalizations: one normalized to the upstream ambient magnetic field (left) and another to the upstream bulk kinetic energy (right). Figure 9 also shows 1D simulation results and the results reported by Gallant et al. (1992) for comparison. The amplitudes of the precursor waves in 2D were systematically smaller than those in the 1D results. This is more or less to be expected. Because the inhomogeneity along the shock surface reduces the coherence of shock-reflected particles, the cold ring distribution assumption may no longer valid and the growth rate and saturation level of the SMI may become smaller than that in 1D simulations. The $\sigma_{e}$ dependence shows a different behavior for $\sigma_{e}\leq 10^{-2}$, which may be attributed to the WI. Although the wave power was smaller than that in 1D results by roughly an order of magnitude in the Weibel-dominated regime, the coherent precursor waves continued to exist with amplitudes sufficiently strong to substantially disturb the upstream medium, confirming that the coherent electromagnetic precursor wave emission is indeed an intrinsic aspect of relativistic shocks.

5 Discussion

5.1 SMI in the Weibel Region

Our simulation results suggest that precursor wave emission via the SMI occurs even in the regime where the WI produces a large-amplitude magnetic turbulence. This sounds rather counterintuitive in the sense that the WI and the SMI both grow from the same unstable distribution function: the dominant WI in a low-$\sigma_{e}$ shock may substantially distort the distribution function before electromagnetic waves are amplified by the SMI. However, in reality, this is not the case and both instabilities may actually coexist to some extent. This may be understood as follows.

We may roughly estimate the effect of the Weibel-generated turbulence on the particle trajectory by comparing the Lorentz force exerted by the ambient $B_{1}$ and the fluctuating magnetic field components $\delta B$. For this, we assume that the fluctuation is random at scale length longer than the coherence length $\lambda$ of the turbulence. By assumption, the turbulent magnetic field effect disappears if the time scale is longer than the particle transit time over the coherence length $\lambda/c$. Therefore, if the Lorentz force arising from the fluctuation during this time scale $\propto\delta B\lambda/c$ is smaller than the average Lorentz force over the unperturbed cyclotron period $B_{1}/\omega_{ce}$:

[TABLE]

the particle performs, on average, a gyromotion with respect to the ambient magnetic field. This condition (Eq. 4) may be written as follows:

[TABLE]

where $\epsilon_{B}=\delta B^{2}/4\pi N_{1}\gamma_{1}m_{e}c^{2}$ is the energy conversion rate from upstream plasma kinetic energy into Weibel-generated magnetic field energy. Equivalently, this may be understood as the condition such that the Larmor radius defined with respect to the fluctuation magnetic field is larger than the coherence length of the turbulence. Therefore, the particle cannot complete a full gyromotion around the turbulent magnetic field.

It is interesting to note that the above condition is solely determined by the properties of the WI, and independent of the shock parameters such as $\sigma_{e}$ or $\gamma_{1}$. Noting that our simulations give $\epsilon_{B}\sim O(0.1)$ and $\lambda/(c/\omega_{pe})\sim O(1)$, we find that the condition (Eq. 5) is always satisfied, and the unperturbed particle gyromotion is sustained. We believe that this is the reason why we observe large-amplitude waves even in the presence of the strong Weibel turbulence.

Nevertheless, the effect of turbulence may not completely be neglected. Since the turbulent magnetic field fluctuation gives random kicks during the gyromotion, the coherence of the particle motion may be broken. In other words, the fluctuation introduces an effective thermal spread in the original otherwise cold ring distribution in momentum. As we already mentioned in $\S$4.2, a finite temperature reduces the growth rate of the SMI at higher harmonics (Amato & Arons, 2006). Therefore, the wave power at higher frequency will be strongly reduced in the Weibel turbulence. This probably leads to the lower wave emission efficiency in the Weibel-dominated regime, but the lower frequency waves being less affected by the turbulence may persist with finite amplitude.

On the other hand, we expect that there will be a lower limit for the magnetization rate below which coherent precursor waves cannot be generated. A lower magnetization rate requires a higher harmonic number $n$ for a generated wave to propagate upstream because the wave frequency must be larger than the cutoff frequency $\omega=n\omega_{ce}\gtrsim\sqrt{2(1+\beta_{shock}\gamma_{shock})}\omega_{pe}$ (see Appendix B). Therefore, the suppression of higher harmonic resonance due to the finite temperature effect becomes progressively important as decreasing $\sigma_{e}$. More quantitative estimate for the lower limit will be addressed in the future.

5.2 Particle Energy Spectra

Figure 10 shows the downstream energy spectra of electrons for $\sigma_{e}=3\times 10^{-1}$ and $3\times 10^{-3}$ in the range $r_{L}<x<X_{sh}-r_{L}$, which are normalized as follows:

[TABLE]

The positron energy spectra are identical to those of electrons. The time evolution from $\omega_{pe}t=50$ to $\omega_{pe}t=500$ is shown. For $\sigma_{e}=3\times 10^{-1}$, the measured distribution reaches a steady state at $\omega_{pe}t=400$ and can be well fitted with the 2D relativistic Maxwellian,

[TABLE]

This indicates that the particles downstream are completely thermalized and that the non-thermal particles are not generated. In contrast, for $\sigma_{e}=3\times 10^{-3}$, a suprathermal tail is clearly visible. The tail gradually approaches a power law with a spectral index of $p\sim 2.6$. The energy spectrum reaches a steady state at $\omega_{pe}t\sim 450$, and the maximum Lorentz factor was approximately at $\gamma_{sat}\sim 600$. Sironi et al. (2013) reported that efficient particle acceleration occurs in a relativistic perpendicular shock propagating in a pair plasma due to the strong turbulence generated by the WI when $\sigma_{e}\lesssim 3\times 10^{-3}$, and that the acceleration efficiency is independent of the bulk Lorentz factor $\gamma_{1}$ if $\gamma_{1}\gtrsim 10$. Note that $\gamma_{1}=40$ in our simulations satisfies this condition. However, the energy spectrum in our simulation saturated faster than that in their simulation possibly because our high-resolution simulation without the application of digital filtering can accurately resolve Weibel-generated turbulence. So far, we have not found any evidence of precursor wave emission contributing positively for the particle acceleration in a pair plasma shock.

5.3 Applicability of WFA

Based on the simulation results, we now discuss the application of the WFA model to astrophysical relativistic shocks. As mentioned earlier, Lyubarsky (2006) and Hoshino (2008) both presented 1D simulation results. In contrast, the numerical resolution of the 2D simulations presented by Sironi & Spitkovsky (2011) was probably insufficient for this purpose. However, another approach was taken by Kuramitsu et al. (2008), who assumed the presence of precursor waves in their 2D PIC simulations. Focusing only on the interaction between the precursor and plasma in the upstream region, they injected large-amplitude electromagnetic waves into a uniform ion-electron plasma and investigated particle acceleration efficiency. Although based on the strong assumption, they found that efficient particle acceleration indeed occurs when the strength parameter of the electromagnetic wave $a=e\delta E/m_{e}c\omega$ is greater than unity, where $\delta E$ is the amplitude of the wave electric field, and $\omega$ is the wave frequency. The accelerated particle exhibited a power-law like spectrum $N(\gamma)\propto\gamma^{-p}$ with a spectral index of roughly $p\simeq 2$.

To discuss the relation with the results by Kuramitsu et al. (2008), we estimated the strength parameter of the precursor waves from our simulation results. As the strength parameter is the amplitude of the particle quiver motion under the wave electromagnetic field, it may be estimated from the transverse particle velocity in the upstream region. As discussed in $\S$4.2, the particles in the upstream region are forced to oscillate in velocity by the strong electromagnetic field of the precursors. In the precursor region defined in $\S$3.1, the oscillation amplitude of the particle velocity was determined by integrating the Fourier spectrum of the transverse velocity (i.e., first-order velocity moment) fluctuations over the wavenumber at $\omega_{pe}=500$. The obtained strength parameter is shown in Figure 11 by the solid line. Another way to estimate the strength parameter is to use the precursor wave amplitude. Typical wavenumbers $kc/\omega_{pe}=2-4$ of the observed precursor wave combined with the dispersion relation give the wave frequency as follows:

[TABLE]

The strength parameter $a$ may then be estimated as

[TABLE]

where $\epsilon_{conv}=\delta B^{2}/8\pi N_{1}\gamma_{1}m_{e}c^{2}$ is the energy conversion rate from the upstream plasma kinetic energy into the precursor wave energy. Figure 11 also shows the strength parameter estimated in this manner by the dashed line. It is readily apparent that the former method (based on the particle quiver velocity) gives a smaller value than the latter method (based on the wave amplitude). This is natural because the latter estimate (Eq. 9) uses $\epsilon_{conv}$ as a proxy of the amplitude of a monochromatic wave, whereas the actual spectra are rather broadband in wavenumber. Nevertheless, it is worth noting that both estimates give the same order of magnitude and essentially the same tendency. The strength parameter $a\gtrsim 1$ implies that the precursor waves are strongly unstable against parametric instability. The formation of the observed density filaments is consistent with this understanding.

We again emphasize that the WFA should not occur in a pair plasma. However, the observation of a strong precursor-plasma interaction in the upstream region indicates that the same precursors will generate large-amplitude Langmuir waves if propagating in an ion-electron plasma.

The latter estimate predicts that the strength parameter of the precursor waves linearly scales with the shock Lorentz factor. Therefore, highly-relativistic shocks associated with, for instance, external shocks of GRBs may be an important site of the application of the WFA model.

6 Summary

In this study, we investigated the efficiency of coherent electromagnetic precursor wave emission by the SMI at relativistic magnetized shocks in pair plasmas via 2D PIC simulations. The precursor wave emission efficiency was measured as a function of the magnetization parameter $\sigma_{e}$, which well characterizes the relativistic shocks. Although the wave power in 2D simulations was systematically smaller than that in the corresponding 1D results, it was sufficiently strong to produce large-amplitude density filaments in the upstream medium. We think that the formation of density filaments is due to a filamentation instability excited by intense electromagnetic waves propagating in the upstream plasma. The precursor wave emission continues even well after the highly-disturbed upstream medium interacts with the shock, and the amplitude reaches a quasi-steady state level. At low $\sigma_{e}$, the power was roughly an order of magnitude smaller than that in 1D simulations, which may be attributed to the presence of the WI dominating the shock transition region. Nevertheless, we found that large-amplitude precursor waves persist even in the Weibel-dominated regime. Therefore, we conclude that the emission of the coherent precursor waves is indeed intrinsic to relativistic shocks, even if the multidimensional effect is considered.

Based on the simulation results, we discussed the applicability of the WFA model to astrophysical relativistic shocks. We concluded that the precursor wave power may be sufficiently strong for the WFA at highly-relativistic shocks in an ion-electron plasma. External shocks in the relativistic jets from GRBs may be important sites for the production of UHECRs via the WFA. However, the actual particle acceleration efficiency must be comprehensively examined by directly performing simulations for relativistic ion-electron shocks. This will be a subject for future research.

Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A Numerical Convergence

We also performed 1D PIC simulations for numerical convergence study. The simulation setup was identical to that described in $\S$2, except for the number of grids and particles per cell. We used the same 2D simulation code with two grid points in the $y$ direction. For technical reasons, this is the minimum number of grid points in the transverse direction.

First, we investigated the dependence on the number of particles per cell $N_{1}\Delta x$ for each species and found that the results did not change as long as $N_{1}\Delta x\gtrsim 64$. However, for a smaller number of particles, the measured amplitude shows non-negligible decrease, probably due to enhanced numerical collisions. In general, the effect of numerical collision becomes smaller as the dimensionality increases (Kato, 2013; May et al., 2014). Therefore, $N_{1}\Delta x=64$ is sufficient for 2D simulations.

Next, we studied the dependence on the grid size $\Delta x$. For this study, we fixed the box size to $N_{x}\Delta x/(c/\omega_{pe})=500$ where $N_{x}$ is the number of gird points. For each $\sigma_{e}$, we performed simulations for the following four different resolutions: $\Delta x/(c/\omega_{pe})=1/10$, $1/20$, $1/40$ and $1/80$.

Figure 12 shows the overall shock structure at $\omega_{pe}t=500$ for $\sigma_{e}=3\times 10^{-1}$ for $\Delta x/(c/\omega_{pe})=1/10$ (left) and $\Delta x/(c/\omega_{pe})=1/40$ (right). Shown from top to bottom are the electron number density $N_{e}$, the $z$ component of magnetic field $B_{z}$, the $x$ component of the electric field $E_{x}$, and the phase-space plots in $x-u_{x}$ and $x-u_{y}$, respectively. All quantities are normalized by the corresponding upstream values. It is clear that the positions of the shock front are different between the two cases: at $x/(c/\omega_{pe})\sim 270$ for $\Delta x/(c/\omega_{pe})=1/10$ and at $x/(c/\omega_{pe})\sim 240$ for $\Delta x/(c/\omega_{pe})=1/40$. We think that the difference in shock propagation velocity arises due to different precursor wave emission efficiencies, as shown in Figure 12. The increased precursor wave emission indicates a larger fraction of the energy flux carried away from the shock, thereby modifying the Rankine-Hugoniot relation. In Appendix B, we show actually that the precursor emission has substantial effect on the shock propagation speed in this parameter range. In contrast, the precursor region ahead of the shock is more extended in higher resolution runs than in lower resolution runs. This is in contrast to the shock propagation speed, and may be explained by analyzing the emitted wave spectrum.

The power spectra of the precursor waves are shown in Figure 13. The spectra were obtained for snapshots at $\omega_{pe}t=500$ by taking the Fourier transform in the region $r_{L}<x-X_{sh}<r_{L}+50c/\omega_{pe}$. The Nyquist wavenumbers are $k_{N}\sim 30\omega_{pe}/c$ and $k_{N}\sim 120\omega_{pe}/c$ for $\Delta x/(c/\omega_{pe})=1/10$ and $\Delta x/(c/\omega_{pe})=1/40$, respectively. Both the maximum wavenumber and the peak power are greater for $\Delta x/(c/\omega_{pe})=1/40$ than for $\Delta x/(c/\omega_{pe})=1/10$. Because the group velocity of electromagnetic waves increases with the wavenumber, the precursor waves in the higher resolution run can propagate farther away from the shock front than at lower resolution. In addition, numerical damping (which is more significant at lower resolution) may also contribute to waves during long-distance propagation.

Figure 14 shows the $\sigma_{e}$ dependence of the wave energy density normalized by the total upstream bulk kinetic energy for $\Delta x/(c/\omega_{pe})=1/10,1/20,1/40$ and $1/80$. The amplitude was calculated by integrating the power spectra over the whole wavenumber space. Error estimates were obtained by taking the standard deviation during the time interval $500\leq\omega_{pe}t\leq 520$. The wave amplitude systematically increases along with the resolution for all $\sigma_{e}$ and more or less converges for $\Delta x/(c/\omega_{pe})\leq 1/40$. This numerical convergence study confirms that the precursor wave emission is indeed very sensitive to numerical resolution. Based on this result, we used $\Delta x/(c/\omega_{pe})=1/40$ for the 2D simulations discussed in the main text.

Appendix B Cutoff Wavenumber

We here estimate a theoretical cutoff wavenumber, above which electromagnetic waves may escape from the shock upstream using the X-mode dispersion relation in a cold pair plasma. In the plasma rest frame, this relation is given by

[TABLE]

where the prime indicates physical quantities measured in the plasma rest frame and $n^{\prime}\equiv ck^{\prime}/\omega^{\prime}$ is the refraction index. By performing Lorentz transformation to the simulation frame, the dispersion relation upstream becomes

[TABLE]

where $\theta$ is the angle between the wave propagation direction and the $x$ axis. As the precursor wave traveling upstream is considered, $\theta$ varies from $-\pi/2$ to $\pi/2$. Therefore, $1+\beta_{1}n\cos\theta$ is always greater than unity, and we can safely neglect the second term in Eq. B2 for $\sigma_{e}/\gamma_{1}^{2}\ll 1$, which is satisfied in our simulations. This leads to the dispersion relation B2, which is identical to that for a weakly magnetized plasma ($\sigma_{e}\ll 1$) in the plasma rest frame:

[TABLE]

Now, the group velocity of the precursor wave can be expressed as

[TABLE]

Equating the $x$ component of the group velocity $v_{g}\cos\theta$ with the shock velocity $c\beta_{shock}$, we obtain the cutoff wavenumber

[TABLE]

Determining the shock velocities is essential to calculate the theoretical cutoff wavenumber. The shock velocities were measured based on the time evolution of the electron number density averaged over the $y$ direction. The $\sigma_{e}$ dependence of the shock velocity is shown in Figure 15. The filled circles are the simulation results, which clearly deviate, particularly at higher $\sigma_{e}$, from the dashed line, showing the ideal MHD Rankine-Hugoniot prediction calculated with an adiabatic index of $\Gamma=3/2$. This is because the electromagnetic precursor waves are not considered in the ideal MHD Rankine-Hugoniot condition. We examined the consistency between the simulations and theory considering the effect of precursor wave emission in the energy and momentum conservation laws. Namely, the shock velocity is calculated by solving the conservation laws using the energy and momentum fluxes of the precursors in the simulations (see Gallant et al., 1992). The open triangles in Figure 15 show the estimate based on the modified Rankine-Hugoniot conditions. This agrees quite well with the observed shock velocity, which we used to estimate the cutoff wavenumber.

Appendix C Weibel Instability

To confirm that the fluctuating magnetic fields near the shock front in the simulation with $\sigma_{e}=3\times 10^{-3}$ result from the WI, we here compare the time evolution of the magnetic field with the linear growth rate of the WI. The maximum value of the fluctuating magnetic field energy averaged over the $y$ axis is determined for each snapshot (solid line in Figure 16). The fluctuating magnetic field energy is normalized by the upstream bulk kinetic energy. For a low $\sigma_{e}$ shock wherein the WI substantially grows in amplitude, we may use a linear growth rate obtained for an unmagnetized plasma. Schaefer-Rolffs & Tautz (2008) discussed the linear theory of the WI for the monochromatic, waterbag, bi-Maxwellian and $\kappa$ distribution functions. In all these distributions, the dispersion relation reduces to

[TABLE]

for sufficiently strong anisotropy. Note that the relativistic effect was appropriately considered when calculating the linear growth rate. The maximum linear growth rate $\Gamma_{max}$ can be estimated from the dispersion relation C1

[TABLE]

which is indicated in Figure 16 by the dashed line. We see that the linear theory is consistent with the simulation result. The energy density of the generated magnetic field saturates at around $10-20\%$ of the upstream bulk kinetic energy. This result is also consistent with those of the previous studies (Kato, 2007; Chang et al., 2008; Sironi & Spitkovsky, 2011). Therefore, we conclude that the filamentary magnetic field structure in the shock transition region results from the WI.