The synchrotron maser emission from relativistic shocks in Fast Radio Bursts: 1D PIC simulations of cold pair plasmas

TL;DR

This study uses 1D PIC simulations to explore synchrotron maser emission from relativistic shocks in cold pair plasmas, offering insights into the FRB emission mechanism within the magnetar scenario.

Contribution

It demonstrates that relativistic shocks can generate narrow-band, line-like precursor waves consistent with FRB observations, advancing understanding of the emission process.

Findings

Shock converts ~7e-4/σ^2 of energy into precursor waves.

Wave spectrum is narrow-band with fractional width 1-3%.

Peak frequency depends on shock Lorentz factor and plasma frequency.

Abstract

The emission process of Fast Radio Bursts (FRBs) remains unknown. We investigate whether the synchrotron maser emission from relativistic shocks in a magnetar wind can explain the observed FRB properties. We perform particle-in-cell (PIC) simulations of perpendicular shocks in cold pair plasmas, checking our results for consistency among three PIC codes. We confirm that a linearly polarized X-mode wave is self-consistently generated by the shock and propagates back upstream as a precursor wave. We find that at magnetizations $\sigma\gtrsim 1$ (i.e., ratio of Poynting flux to particle energy flux of the pre-shock flow) the shock converts a fraction $f_\xi' \approx 7 \times 10^{-4}/\sigma^2$ of the total incoming energy into the precursor wave, as measured in the shock frame. The wave spectrum is narrow-band (fractional width $\lesssim 1-3$), with apparent but not dominant line-like features as many resonances concurrently contribute. The peak frequency in the pre-shock (observer) frame is $\omega^{\prime \prime}_{\rm peak} \approx 3 \gamma_{\rm s | u} \omega_{\rm p}$, where $\gamma_{\rm s|u}$ is the shock Lorentz factor in the upstream frame and $\omega_{\rm p}$ the plasma frequency. At $\sigma\gtrsim1$, where our estimated $\omega''_{\rm peak}$ differs from previous works, the shock structure presents two solitons separated by a cavity, and the peak frequency corresponds to an eigenmode of the cavity. Our results provide physically-grounded inputs for FRB emission models within the magnetar scenario.

Figures

Tables

Table 1. Table 1: Typical parameters of the PIC simulations presented in this study: Δ ​ t Δ 𝑡 \Delta t is the time-step in units of the inverse plasma frequency ω p − 1 superscript subscript 𝜔 p 1 \omega_{\rm p}^{-1} (defined with both species), T sim subscript 𝑇 sim T_{\rm sim} is the simulation timespan, Δ ​ x Δ 𝑥 \Delta x is the cell size in units of c / ω p 𝑐 subscript 𝜔 p c/\omega_{\rm p} , k B ​ T 0 subscript 𝑘 𝐵 subscript 𝑇 0 k_{B}T_{0} is the upstream thermal energy in units of m e ​ c 2 subscript 𝑚 𝑒 superscript 𝑐 2 m_{e}c^{2} , N ppc subscript 𝑁 ppc N_{\rm ppc} is the number of particles-per-cell for each species, σ min subscript 𝜎 min \sigma_{\rm min} and σ max subscript 𝜎 max \sigma_{\rm max} are the minimal and maximal values of the magnetization explored with a given code.

PIC code	$\Delta t{\omega }_p$	$T_{\mathrm{sim}}{\omega }_p$	$\Delta x/\left(c/{\omega }_p\right)$	$N_{\mathrm{ppc}}$	$k_B{T}_0/\left(m_e{c}^2\right)$	${\sigma }_{\min }$	${\sigma }_{\max }$	${\gamma }_0$
Tristan-MP	1/200	$\lesssim 2\times {10}^4$	1/100	64	${10}^{-4}$	${10}^{-1}$	30	10
Smilei	$1/224$ & $1/113$	$6.7\times {10}^3$	1/112	20	${10}^{-4}$	${10}^{-1}$	30	10
Shockapic	1/90	$1.2\times {10}^3$	1/44.7	20	${10}^{-4}$	${10}^{-3}$	1	10
Smilei (2)	$1/90$	$1\times {10}^3$	1/44.7	20	${10}^{-6}$	${10}^{-3}$	2	10
Smilei (3)	$1/90$	$1.5\times {10}^3$	1/44.7	20	${10}^{-4}$	${10}^{-3}$	2	160

Equations

$$\sigma =\frac{B_{z,0}^2}{8\pi {\gamma }_0{N}_0{m}_e{c}^2}={\left(\frac{{\omega }_{\mathrm{c}}}{{\omega }_{\mathrm{p}}}\right)}^2,$$

$$⟨\delta B_z^2⟩=⟨(B_z-B_0{)}^2⟩.$$

$${\xi }_B=\frac{⟨\delta B_z^2⟩}{B_0^2}=\frac{⟨\delta B_z\delta E_y⟩}{B_0^2},$$

$${\gamma }_0^2\left({\beta }_0+{\beta }_{s|d}\right)\left(w_u+\frac{b_{0,u}^2}{4\pi }\right)-\left(1-{\beta }_{s|d}\right)\frac{\delta B_u^2}{4\pi }=$$

$${\beta }_{s|d}\left(w_d-p_d+\frac{b_{0,d}^2}{8\pi }+\frac{\delta B_d^2}{4\pi }\right),$$

$$f_{\xi }={\xi }_B\left(\frac{\sigma }{1+\sigma }\right)\left(\frac{1-{\beta }_{s|d}}{{\beta }_0+{\beta }_{s|d}}\right),$$

$$g_{\xi }=f_{\xi }(1+\sigma ).$$

$$\frac{k^{\prime \prime 2}{c}^2}{{\omega }^{\prime \prime 2}}$$

$$k^2{c}^2\simeq {\omega }^2-{\omega }_{\mathrm{p}}^2.$$

$$c\sqrt{1-\frac{{\omega }_{\mathrm{p}}^2}{{\omega }^2}}={\beta }_{s|d}c.$$

$${\omega }_{\mathrm{cutoff}}$$

$${\lambda }_{\mathrm{cutoff}}$$

$$\frac{B_{\mathrm{sh}}^{\prime }}{B_0^{\prime }}\approx \sqrt{1+\frac{2}{\sigma }},$$

$${\omega }_{\mathrm{c},\mathrm{sol}}\approx \left(\sqrt{\sigma +2}+\sqrt{\sigma +2-{\beta }_{s|d}^2}\right){\gamma }_{s|d}{\omega }_{\mathrm{p}}.$$

$$u_x$$

$$u_y$$

$$a=\frac{e\delta E_y}{m_{e}c\omega }$$

$$B_0^{\prime }={\gamma }_{s|d}\left(B_0+{\beta }_{s|d}{E}_0\right)={\gamma }_{s|d}{B}_0\left(1+{\beta }_{s|d}{\beta }_0\right),$$

$$\delta B_z^{\prime }={\gamma }_{s|d}\left(\delta B_z-{\beta }_{s|d}\delta E_y\right)\simeq {\gamma }_{s|d}\delta B_z\left(1-{\beta }_{s|d}\right)$$

$${\xi }_B^{\prime }={\xi }_{B|\mathrm{sh}}=\frac{⟨\delta B_z^{\prime 2}⟩}{B_0^{\prime 2}}={\left(\frac{1-{\beta }_{s|d}}{1+{\beta }_0{\beta }_{s|d}}\right)}^2\frac{⟨\delta B_z^2⟩}{B_0^2}={\left(\frac{1-{\beta }_{s|d}}{1+{\beta }_0{\beta }_{s|d}}\right)}^2{\xi }_B.$$

$${\gamma }_u^{\prime 2}{\beta }_u^{\prime }{W}_u-\frac{\delta B_u^{\prime 2}}{4\pi }={\gamma }_d^{\prime 2}{\beta }_d^{\prime }{W}_d+\frac{\delta B_d^{\prime 2}}{4\pi }.$$

$$\frac{\delta B_u^{\prime 2}}{4\pi }$$

$$\frac{\delta B_d^{\prime 2}}{4\pi }$$

$$f_{\xi }^{\prime }=f_{\xi |\mathrm{sh}}=\frac{\mathrm{EM}\mathrm{wave}\mathrm{energy}}{\mathrm{Total}\mathrm{incoming}\mathrm{energy}}=\frac{{\xi }_B^{\prime }}{{\beta }_u^{\prime }}\frac{\sigma }{1+\sigma },$$

$$g_{\xi }^{\prime }=\frac{\mathrm{EM}\mathrm{wave}\mathrm{energy}}{\mathrm{Incoming}\mathrm{kinetic}\mathrm{energy}}=\frac{\delta B_u^{\prime 2}}{4\pi {\gamma }_u^{\prime 2}{\beta }_u^{\prime }{n}_{u}m{c}^2}=\frac{{\xi }_B^{\prime }}{{\beta }_u^{\prime }}\sigma .$$

$${\omega }_{\mathrm{peak}}^{\prime \prime }\approx 3{\gamma }_{s|u}{\omega }_{\mathrm{p}},$$

$${\nu }^{\prime \prime }=\frac{{\omega }^{\prime \prime }}{2\pi }\approx 2.7\times 10^4{\gamma }_{s|u}{\left(\frac{n_e}{1{\mathrm{cm}}^{-3}}\right)}^{1/2}\mathrm{Hz},$$

Full text

††pagerange: The synchrotron maser emission from relativistic shocks in Fast Radio Bursts: 1D PIC simulations of cold pair plasmas–LABEL:lastpage††pubyear: 2018

1 Introduction

Synchrotron masers are known to produce strong decametric radio emission in the Jovian magnetosphere and kilometric emission in the terrestrial magnetosphere (Auroral Kilometric Emission). They are driven by a “population inversion” of energetic electrons gyrating in an intense magnetic field. The driver of the emission is either a loss-cone or ring-like electron distribution function (Melrose, 2017; Treumann, 2006). Such a population inversion occurs also in strongly magnetized relativistic perpendicular shocks, where a coherent cold-ring distribution of particles is self-consistently produced as part of the shock evolution (Gallant et al., 1992; Amato & Arons, 2006). This distribution is unstable to the synchrotron maser instability (Hoshino & Arons, 1991) that causes the emission of a train of high amplitude semi-coherent electromagnetic waves propagating from the shock front into the unshocked (upstream) medium (Gallant et al., 1992; Hoshino et al., 1992). The possible importance of the synchrotron maser in astrophysical sources was anticipated long ago (Sazonov, 1973). Yet, to our present knowledge, there is no firm demonstration that shock-powered synchrotron maser emission is dominant in any non-heliospheric astrophysical environment, even though a few potential scenarios have been proposed (e.g., Sagiv & Waxman, 2002).

Recently, however, the discovery of Fast Radio Bursts (FRBs; Lorimer et al., 2007; Keane et al., 2012; Thornton et al., 2013; Spitler et al., 2014; Marcote et al., 2017) has revived the interest in this mechanism. These events are bright ($\sim 1\,$Jy) pulses of millisecond duration detected in the $\sim$ GHz band. Their extremely high brightness temperature, $T_{\rm B}\sim 10^{37}$K, requires a coherent emission mechanism (Katz, 2016; Popov et al., 2018). Young magnetars have emerged as one of the most likely progenitors of FRBs, at least for the repeating class (e.g., Popov & Postnov, 2013; Murase et al., 2016; Lyutikov, 2017; Metzger et al., 2017; Kashiyama & Murase, 2017; Long & Pe’er, 2018; Margalit et al., 2018; Margalit & Metzger, 2018). Magnetars can naturally explain the short FRB durations, large energy requirements, and ordered magnetic fields needed for coherent emission. Thus far, most works on the FRB emission mechanism are based on mere considerations of energetics and timescales, in which it is assumed that a fraction of the free energy of the system is radiated away at GHz frequencies by coherent charge “bunches” via, e.g., curvature or synchrotron maser processes. In the case of curvature radiation, the emission is postulated to be a product of magnetic reconnection close to the magnetar surface (e.g., Lyutikov, 2002; Kumar et al., 2017; Lu & Kumar, 2018; Ghisellini & Locatelli, 2018; Katz, 2018). In the case of the synchrotron maser, the emission is thought to occur at relativistic shocks propagating in the magnetar wind or nebula (Lyubarsky, 2014; Murase et al., 2016; Beloborodov, 2017) or inside the ultra-relativistic shell ejected from the central compact object (Waxman, 2017). Yet, in either case the conditions for coherent emission and the very existence of charge bunches with the required properties are often postulated ad hoc, resulting in models with little predictive power.

The purpose of this work is to demonstrate from first principles that the synchrotron maser at relativistic shocks in the magnetar wind can naturally explain the observed FRB properties. By means of particle-in-cell (PIC) simulations, we investigate how the efficiency and spectrum of the electromagnetic wave emitted by the shock into the pre-shock medium (which we shall call “precursor wave”) depend on the physical conditions in the magnetar wind. In this work, the first of a series, we present results from one-dimensional (1D) simulations (more precisely, 1D3V, i.e., we employ one spatial dimension, but all three components of velocities and electromagnetic fields are retained), while multi-dimensional runs will be presented in a future work (Sironi et al., in prep.; see also Appendix B, for the precursor energetics in 2D and 3D). We focus on the case of a cold pair-dominated plasma.

There exists extensive literature on PIC modeling of the electromagnetic precursor wave in relativistic perpendicular shocks (e.g., Langdon et al., 1988; Gallant et al., 1992; Hoshino et al., 1992; Amato & Arons, 2006; Hoshino, 2008; Sironi & Spitkovsky, 2009; Iwamoto et al., 2017, 2018). Our work is motivated by the poor exploration of the extreme regime where the energy content of the plasma is dominated by magnetic fields, as it is supposedly the case in magnetar winds. In other words, we focus on magnetizations $\sigma\gtrsim 1$, where $\sigma$ is the ratio of upstream Poynting flux to kinetic energy flux. For $\sigma\gtrsim 1$, 1D simulations are adequate, since we find that they agree well with multi-dimensional results (see Appendix B, for the precursor energetics in 2D and 3D; also Sironi et al., in prep.). We provide an extensive investigation of the dependence on the flow magnetization, from $\sigma=0.1$ to $\sigma=30$, with much longer simulations than previously reported, especially in the $\sigma\gtrsim 1$ regime relevant for FRB sources. Previous works arguably never reached a steady state in simulations with $\sigma\gtrsim 1$. We employ several PIC codes (Tristan-MP, Smilei and Shockapic) to check for consistency, and thus confirm the robustness of our results.

At $\sigma\gtrsim 1$ the shock converts a fraction $f_{\xi}\approx 2\times 10^{-3}/\sigma$ of the total incoming energy into the precursor wave, as measured in the post-shock (downstream) frame. In the shock rest frame, the efficiency is $f_{\xi}^{\prime}\approx 7\times 10^{-4}/\sigma^{2}$. The spectrum of the precursor wave is narrow-band, $\Delta\omega/\omega_{\rm peak}\lesssim 1-3$, with apparent but not dominant line-like features as many resonances concurrently contribute. The peak frequency scales in the post-shock frame as $\omega_{\rm peak}\simeq 3\,\omega_{\rm p}\max[1,\sqrt{\sigma}]$, where $\omega_{\rm p}$ is the plasma frequency. In the pre-shock frame (which coincides with the observer frame, if the magnetar wind is non-relativistic), this can be recast in a simpler form as $\omega^{\prime\prime}_{\rm peak}\approx 3\gamma_{\rm s|u}\omega_{\rm p}$, where $\gamma_{\rm s|u}$ is the shock Lorentz factor in the upstream frame. At $\sigma\gtrsim 1$, where our estimated $\omega_{\rm peak}$ differs from earlier works (that quoted $\omega_{\rm peak}\propto\sigma\omega_{\rm p}$, see Gallant et al. (1992), rather than $\omega_{\rm peak}\propto\sqrt{\sigma}\omega_{\rm p}$ as we find) we see that the shock structure displays two solitons separated by a cavity, and the peak frequency of the spectrum corresponds to an eigenmode of the cavity. The efficiency and spectrum of the precursor wave do not depend on the bulk Lorentz factor of the pre-shock flow.

The paper is organized as follows. In section 2 we present the methods and the numerical setup. In section 3 we discuss the main results, in the post-shock (downstream) rest frame. Section 4 discusses the energy content of precursor waves in the frame of the shock front. In section 5 we present the implications of our results for FRB emission models, and we conclude in section 6.

2 Simulation methods and setup

We use the particle-in-cell (PIC) codes Tristan-MP (Spitkovsky, 2005; Sironi & Spitkovsky, 2009), Smilei (Derouillat et al., 2018), and Shockapic (Plotnikov et al., 2018) to perform 1D3V simulations, where we retain one spatial direction, but all three components of velocities and electromagnetic fields. Mainly Tristan-MP and Smilei are employed for large simulations. The pseudo-spectral code Shockapic is used to check for consistency in shorter simulations. In the main body of the paper, we present only results obtained with Tristan-MP, unless stated otherwise. In the Appendix A, we demonstrate the agreement between different codes across the whole range of $\sigma$ explored in this work.

The use of a reduced 1D spatial geometry is justified in the limit of magnetically-dominated plasmas, as we will demonstrate with 2D and 3D simulations in a forthcoming study (Sironi et al, in prep.). For lower magnetizations than explored here, i.e., $\sigma\lesssim 0.5$, Iwamoto et al. (2017) found that the precursor wave energy is reduced by at most an order of magnitude in 2D simulations as compared to 1D. This difference is smaller or even negligible in the high magnetization limit explored here (see Appendix B).

The shock is initialized using the common setup described in, e.g., Spitkovsky (2008), which we summarize here for completeness. The upstream flow, composed of electrons and positrons, drifts along the $-\hat{x}$ direction with a speed $-\beta_{0}c\hat{x}$. The corresponding bulk Lorentz factor is $\gamma_{0}=(1-\beta_{0}^{2})^{-1/2}=10$, but we have also explored higher values of $\gamma_{0}$, up to $10^{5}$. The upstream pair plasma is cold, with thermal spread $k_{B}T_{0}/(m_{e}c^{2})=10^{-4}$. The pre-shock plasma carries a frozen-in magnetic field $B_{0}$ oriented along $z$ (so, $B_{z,0}=B_{0}$), i.e., perpendicular to the flow propagation, and a motional electric field $E_{y,0}=-\beta_{0}B_{z,0}$. The flow is reflected at a wall located at $x=0$. After some time (at least several cyclotron periods $\omega_{\rm c}^{-1}$), the shock front forms by magnetic reflection and steadily propagates along the $+\hat{x}$ direction with a speed that is in good agreement with the Rankine-Hugoniot conditions. The resulting simulation frame coincides with the frame where the downstream plasma is at rest (downstream rest frame; DRF).

The shock physics is sensitive to the upstream magnetization $\sigma$, which we define as the ratio of Poynting to kinetic energy flux

[TABLE]

and we vary from $\sigma=0.1$ up to $\sigma=30$. Here, $N_{0}$ is the number density of upstream electrons (the overall particle number density is then $2N_{0}$), $m_{e}$ is the electron (or positron) mass, and $c$ is the speed of light in vacuum. We also define the typical gyro-frequency as $\omega_{\rm c}=|q|B_{0}/(\gamma_{0}m_{e}c)$ and the plasma frequency as $\omega_{\rm p}=[8\pi N_{0}q^{2}/(\gamma_{0}m_{e})]^{1/2}$, where $q$ is the elementary electric charge. Both quantities are based on the upstream values of magnetic field and plasma density measured in the simulation frame.

Typical numerical parameters used with Tristan-MP are:

• •

The skin depth is well resolved with $c/\omega_{\rm p}=100\,\Delta$, where $\Delta$ is the grid size. This ensures that the typical particle gyro-radius $\simeq\sigma^{-1/2}\,c/\omega_{\rm p}$ is well resolved even for the largest magnetization $\sigma=30$ explored in this study. A good spatial resolution is also essential to capture the high-frequency part of the spectrum of precursor waves (see also Iwamoto et al., 2017, for a discussion on the required spatial resolution).

• •

The simulation time-step is defined as $c\Delta t=0.5\,\Delta$, corresponding to a time resolution of $5\times 10^{-3}\omega_{\rm p}^{-1}$.

• •

The number of particles per cell initialized in the upstream plasma is $N_{ppc}=64$ per species (values between 20 and 200 were tested with no appreciable differences).

• •

The simulation is evolved up to $T_{\rm sim}\simeq 1.5\times 10^{3}\omega_{\rm p}^{-1}$ for $\sigma<1$ and for longer times (up to $2\times 10^{4}\omega_{\rm p}^{-1}$) for $\sigma\gg 1$. This is required in order to reach a steady state in which the precursor emission maintains a constant amplitude (see section 3.2).

In the Appendix A we also report the typical simulation parameters for the other two PIC codes used in this study. They are not presented here since in the following sections we mainly discuss the results obtained with Tristan-MP.

3 Results

In this section we explore the physics of relativistic highly magnetized electron-positron shocks, focusing on the properties of the electromagnetic precursor. In section 3.1 we discuss the typical shock structure and show the presence of electromagnetic precursor waves. In sections 3.2 and 3.3 we show the dependence on $\sigma$ of the precursor wave intensity and spectrum, respectively. The wave strength parameter is discussed in section 3.4.

3.1 Shock layer structure

As shown in several works using 1D simulations (Langdon et al., 1988; Gallant et al., 1992; Amato & Arons, 2006; Lyubarsky, 2006; Hoshino, 2008), 2D simulations (Sironi & Spitkovsky, 2009, 2011; Iwamoto et al., 2017, 2018; Plotnikov et al., 2018) and 3D simulations (Spitkovsky, 2005; Sironi et al., 2013), highly magnetized perpendicular relativistic shocks form by magnetic reflection and generate a strong electromagnetic wave propagating from the shock into the upstream region. For high magnetizations (typically, $\sigma\gtrsim 0.01$), the transverse Weibel filamentation instability, which dominates for $\sigma\lesssim 10^{-3}$, plays no significant role in shaping the shock structure. This partly justifies the 1D approach adopted here.

In Figure 1 we present the structure of the shock transition region for two representative magnetizations. The left column (panels a-e) presents a shock with $\sigma=0.3$ and the right column (panels f-j) corresponds to $\sigma=3$. The timespan of the simulations is long enough to reach a stationary state: we show results at $\omega_{\rm p}t=540$ for $\sigma=0.3$ and at $\omega_{\rm p}t=1800$ for $\sigma=3$. From top to bottom we show the electron number density $N_{e}/N_{0}$ (panels a and f), the transverse magnetic field $B_{z}/B_{0}$ (b and g), the transverse electric field $E_{y}/B_{0}$ (c and h), the longitudinal positron phase space $x-u_{x}$ (d and i), and the transverse positron phase space $x-u_{y}$ (e and j). Here, we define $u_{\alpha}=\gamma\beta_{\alpha}$ as the dimensionless four-velocity. The phase space of electrons is identical to the one of positrons in virtue of mass symmetry, except for the opposite sign in variations of $u_{y}$. The electrostatic field $E_{x}$ is not plotted since it is completely negligible in pair plasmas (we have systematically checked this conclusion). The vertical dashed lines in panels (b) and (g) delimit the region where we have extracted the wave properties, such as amplitude and spectrum, that will be discussed in the sections below. Small insets in the upper right side of panels (d) and (i) show the particle distribution in momentum space $u_{x}-u_{y}$ at the location of the shock front.

The shock front is located at $x-x_{\rm front}=0$ in Figure 1. The upstream flow is on the positive side $(x-x_{\rm front}>0)$ and the downstream plasma is on the negative side $(x-x_{\rm front}<0)$. The existence of a well developed shock is confirmed by the jump in the electron number density and in the $B_{z}$ field at the front location. The shock front itself exhibits a soliton-like structure (see, e.g., Alsop & Arons, 1988), where the particle distribution forms a semi-coherent cold ring in momentum space (see insets in panels d and i). The presence of a large amplitude electromagnetic precursor wave is evidenced in the upstream region of the $B_{z}/B_{0}$ and $E_{y}/B_{0}$ panels, for both magnetizations (see the ripples in the $x-x_{\rm front}>0$ region). The precursor wave amplitude is larger for $\sigma=0.3$ than for $\sigma=3$. This wave is steadily emitted from the shock front and is linearly polarized. The wave vector $\mn@boldsymbol{k}$ lies along the shock direction of propagation (i.e., along $x$), the fluctuating magnetic field is along $z$ (i.e., along the same direction as the upstream field $\mn@boldsymbol{B}_{0}=B_{0}\hat{z}$), and the fluctuating electric field is perpendicular to both $\mn@boldsymbol{k}$ and $\mn@boldsymbol{B}_{0}$. The wave is then identified with the extraordinary mode (X-mode). The phase velocity of the wave is slightly superluminal, as expected for X-mode propagation in a plasma, while its electromagnetic nature is confirmed by the fact that the space-averaged $\langle\delta B_{z}^{2}\rangle=\langle\delta B_{z}\delta E_{y}\rangle$.

We note that the field-aligned component of the particle momentum $u_{z}$ is not affected by the shock. The incoming particles are efficiently isotropized in the $xy$ plane perpendicular to the field, but the post-shock particle distribution remains largely confined to this plane. It follows that the downstream effective adiabatic index corresponds to a 2D relativistically hot gas, $\Gamma_{\rm ad}=3/2$, instead of $4/3$ if the downstream plasma were isotropic in all momentum directions. The lack of isotropization is due to the fact that in a $\sigma\gg 1$ flow (with downstream plasma magnetically dominated) it will be harder for the plasma to exceed the threshold for velocity-space instabilities that feed off the particle temperature anisotropy. For example, the plasma will go unstable via the mirror mode if the temperature anisotropy is above a threshold that scales as $\propto\sigma$, which is harder to exceed at higher magnetizations. In addition, the 1D spatial geometry employed here will further suppress the growth of field-aligned modes leading to momentum isotropization.

The downstream particle energy spectrum (not shown) resembles a 2D Maxwell-Jüttner distribution whose temperature is slightly lower than the one expected from the Rankine-Hugoniot jump conditions (the difference is due to the energy transferred to the precusor waves). No non-thermal tail is observed for the runs presented in this work. This is in agreement with the inefficiency of particle acceleration expected at relativistic strongly magnetized perpendicular shocks (Sironi & Spitkovsky, 2009; Lemoine & Pelletier, 2010; Sironi et al., 2013, 2015; Pelletier et al., 2017; Iwamoto et al., 2017; Plotnikov et al., 2018).

3.2 Precursor wave energy

We now focus on the dependence of the precursor properties, and specifically of its amplitude, on the upstream magnetization $\sigma$. The dependence on the upstream bulk Lorentz factor $\gamma_{0}$ will be discussed in the last part of this subsection.

3.2.1 Temporal evolution of the precursor wave

After an initial transient — whose duration depends on the upstream magnetization, as we show below — the intensity of the precursor wave settles to its asymptotic value. We measure the wave intensity in a region between 5 and 25 $c/\omega_{\rm p}$ ahead of the shock front:111This region is delimited by vertical black dashed lines in panels (b) and (g) of Figure  1. $5\,c/\omega_{\rm p}<x-x_{\rm front}<25\,c/\omega_{\rm p}$. This region is far enough from the shock not to be affected by the front structure itself, and it contains a large number of precursor wavelengths so that we can obtain a solid measure of the precursor average properties.

The wave intensity is then calculated as the spatial average

[TABLE]

In Figure 2 we show for different magnetizations the time evolution of the normalized wave intensity, defined as

[TABLE]

where $\delta E_{y}=E_{y}-E_{y,0}=E_{y}+\beta_{0}B_{0}$. Different solid lines correspond to different values of the magnetization, from $\sigma=0.1$ (blue line) to $\sigma=10$ (black line), as indicated in the legend. By computing the temporal variation of the wave intensity, we can assess when the precursor wave has reached a steady state. Figure 2 shows that:

• •

With increasing $\sigma$, a longer time is required for the precursor to settle at its time-asymptotic state. This is due to the combination of two effects. First, the shock velocity increases from $\beta_{\rm s|d}=0.476$ for $\sigma=0.1$ to $\beta_{\rm s|d}=0.963$ for $\sigma=10$, so at higher magnetizations it takes more time for the precursor wave, propagating at $\beta_{\rm wave}\simeq 1$, to detach from the shock front. Second, there is some interaction occurring between the wave and the upstream plasma, which initially causes a drop in wave efficiency (e.g., at $\omega_{\rm p}t\sim 1100$ in the black line of Figure 2). The time required for the wave to self-regulate and settle to a steady state, following this drop, is longer for higher $\sigma$ (e.g., compare green and black lines in Figure 2).

• •

The steady-state value of the normalized wave energy $\xi_{B}$ decreases with increasing magnetization and for $\sigma\gg 1$ (cyan, brown and black lines) it approaches a constant value $\xi_{B}\simeq 0.01$.

• •

For $\sigma=0.3$ (orange) and $\sigma=0.4$ (yellow), which will be called “transition cases” in the following, the wave intensity varies between periods of high efficiency and phases of low efficiency.

• •

All the simulations have been evolved for long enough to reach a quasi-stationary state. For the largest explored magnetization $\sigma=30$, the simulation was advanced beyond $2\times 10^{4}~{}\omega_{\rm p}^{-1}$ (this case is not shown in the figure but reported in subsequent figures).

Once the wave intensity has settled to a steady state, we have extracted a number of wave properties, such as the energy, the spectrum (peak frequency, low-frequency cutoff, spectral width), and the wave strength parameter, as we now describe.

3.2.2 Dependence on the upstream magnetization

In the previous section we have defined the energy fraction in upstream field fluctuations $\xi_{B}$ (see Eq. 3), as the ratio of the precursor wave magnetic energy to the background field energy. In order to get a global idea of the energetics, we also need to complement it with a parameter that quantifies the energy fraction of the incoming plasma (including both kinetic and electromagnetic content) that is radiated from the shock front in the form of precursor waves. When the electromagnetic field fluctuations induced by the precursor are taken into account in the jump conditions across the shock, the energy conservation equation expressed in the simulation frame is (Gallant et al., 1992; Plotnikov et al., 2018):

[TABLE]

where the subscripts ‘$u$’ and ‘$d$’ refer to the upstream and downstream regions, respectively. Here, $w_{i}$ and $b_{0,i}$ are respectively the fluid enthalpy density and mean magnetic field, both measured in the fluid rest-frame. As above, $\beta_{\rm s|d}$ is the velocity of the shock front as measured in the downstream frame of the simulations. The fluctuating components $\delta B_{u}$ and $\delta B_{d}$ are measured in the simulation frame. We have made the approximation of negligible thermal pressure upstream (strong shock limit) and we have assumed that electrostatic effects are negligible both upstream and downstream ($\delta E_{x}\approx 0$). The latter approximation is fully supported by the simulations and, more fundamentally, by the fact that space-charge effects are expected to be negligible in pair plasmas. The strong shock limit means that the upstream plasma pressure can be neglected and the upstream fluid enthalpy density (in the fluid rest-frame) is then $w_{u}=n_{u}m_{e}c^{2}$, where $n_{u}$ is the upstream plasma proper density. The mean upstream magnetic field $b_{0,u}$ in the upstream frame is related to the pre-shock magnetic field $B_{z,0}$ in the simulation frame via a Lorentz boost: $B_{z,0}=\gamma_{0}b_{0,u}$. Hence, the magnetization parameter can be rewritten as $\sigma=b_{0,u}^{2}/(4\pi w_{u})$ and we note that $\delta B_{u}^{2}/b_{0,u}^{2}=\gamma_{0}^{2}\delta B_{u}^{2}/B_{z,0}^{2}=\gamma_{0}^{2}\xi_{B}$, because the fluctuating part is measured directly in the simulation frame.

As we focus on the precursor wave propagating upstream, here we only consider the left hand side of equation 4.222Since we forego the discussion of the downstream part of the energy conservation equation, an interested reader will find details in the aforementioned works (Gallant et al., 1992; Plotnikov et al., 2018). We find that the fraction of total incoming energy (including both particle and electromagnetic contributions) that is channeled into the precursor wave can be expressed as

[TABLE]

as seen from the DRF. In the following, $f_{\xi}$ will be identified as the “energy fraction parameter”. It is also convenient to define the fraction of incoming particle kinetic energy that is converted into precursor emission

[TABLE]

The time-asymptotic values of $\xi_{B}$, $\beta_{\rm s|d}$, and $f_{\xi}$ measured in our simulations are presented in panels (a), (b) and (d) of Figure 3, as a function of magnetization. Error bars indicate the standard deviation of our time measurements. As regard to $\xi_{B}$, we observe a rapid decrease from $\xi_{B}\simeq 3.53$ at $\sigma=0.1$ down to $\xi_{B}\simeq 0.04$ at $\sigma=1$. The inflection point of the transition occurs at $\sigma\simeq 0.35$. This decrease accompanies a change in the shock front structure that for $\sigma>1$ presents a coherent soliton-like shape (compare left and right columns in figure 1 at the shock). For $\sigma>1$, $\xi_{B}$ slowly decreases and eventually approaches a constant value $\xi_{B}\simeq 10^{-2}$.

Concerning the shock front speed and its corresponding bulk Lorentz factor, $\beta_{\rm s|d}$ and $\gamma_{\rm s|d}$, panels (b) and (c) demonstrate an excellent agreement between our measured values, plotted with blue symbols, and the predictions of ideal MHD jump conditions (e.g., Appendix B of Plotnikov et al., 2018), as indicated by the red solid line. The front speed increases from $\beta_{\rm s|d}=0.476c$ for $\sigma=0.1$ up to $\beta_{\rm s|d}=0.987c$ for $\sigma=30$. The Lorentz factor of the shock front tends asymptotically to $\gamma_{\rm s|d}=\sqrt{\sigma}$, for $\sigma\gg 1$. We remark that the MHD equations used here to derive the jump conditions do not incorporate modifications due to the precursor wave. The accurate agreement of our results with ideal MHD jump conditions for $\sigma\gg 0.1$ is then due to the fact that at high magnetizations the precursor wave is relatively weak, and it does not have an appreciable dynamical effect on the shock. In contrast, in the case when the precursor wave is the strongest, $\sigma=0.1$, the agreement is the worst, because the emission of the large amplitude wave can slow down the shock front, as compared to the ideal MHD prediction.

The dependence on $\sigma$ of the energy fraction parameter $f_{\xi}$ is presented in panel (d) of figure 3. It was calculated by plugging the values from panels (a) and (b) into equation 5. It shows that the energy fraction in the precursor wave decreases from $10$% for $\sigma=0.1$ down to $0.0065$% for $\sigma=30$. The dashed orange line follows the empirical scaling $f_{\xi}=2\times 10^{-3}/\sigma$ that satisfactorily fits our measured values in the $\sigma>1$ range. The most noticeable result of this panel is that we observe a well-defined scaling $f_{\xi}\propto\sigma^{-1}$. This result arises from the fact that for $\sigma\gg 1$, the normalized wave intensity $\xi_{B}$ is roughly constant and $\beta_{\rm s|d}\simeq 1-1/(2\sigma)$. It follows that in the limit $\sigma\gg 1$ the precursor wave carries a constant fraction of the incoming particle kinetic energy, i.e., $g_{\xi}\simeq 2\times 10^{-3}$.

Let us emphasize, however, that this $\sigma$-dependence of $f_{\xi}$ and $g_{\xi}$ is derived in the DRF (simulation frame). This dependence will be different in the shock front rest frame, since the front moves with ultra-relativistic speeds for $\sigma\gg 1$. This point will be further discussed in section 4.

3.2.3 Dependence on the upstream bulk Lorentz factor

So far, we have investigated the dependence of the precursor intensity on $\sigma$, for a fixed choice of the upstream flow Lorentz factor $\gamma_{0}=10$. Here, we demonstrate that $\xi_{B}$ is essentially independent from $\gamma_{0}$, for any value of $\sigma$. Let us first consider the dependence at a fixed $\sigma$. In Figure 4 we show the time evolution of the precursor wave energy for $\sigma=3$, when varying $\gamma_{0}$ from $5$ to $80$. Lines of different color correspond to different values of $\gamma_{0}$. Despite large oscillations in time, it appears that $\xi_{B}$ converges to the same value, regardless of $\gamma_{0}$. The time-asymptotic values of $\xi_{B}$, with corresponding error bars, are plotted in the figure inset. Within the error bars, we can assert that there is no obvious dependence on $\gamma_{0}$.

In order to generalize this conclusion to any $\sigma$, it is worth noting that in the seminal study of Gallant et al. (1992), two very different values of the bulk Lorentz factor ($\gamma_{0}=40$ and $10^{6}$) were used, for a range of $\sigma\in[10^{-3},5]$. The authors did not notice any dependence on $\gamma_{0}$. Also, Iwamoto et al. (2017) performed 1D simulations with $\gamma_{0}=40$ and explored $\sigma$ values between $10^{-3}$ and $0.5$, finding similar values as in Gallant et al. (1992).

In the Appendix A, Figure 11 shows the values of $\xi_{B}$ obtained for $\sigma\in[10^{-3},1]$ (horizontal axis) and for $\gamma_{0}$ ranging from $10$ to $10^{6}$ (different datasets). This figure shows that in the low magnetization regime $\sigma\in[10^{-3},0.3]$, the normalized wave intensity $\xi_{B}$ is nearly independent from $\gamma_{0}$. In the range $\sigma\in[0.3,1]$ there is a larger scatter among different datasets (which employ different $\gamma_{0}$). This range of magnetizations corresponds to the transition cases (see Figure 2). The most plausible reason for the discrepancy among different datasets is that the simulations from earlier studies were not evolved long enough in order to reach the asymptotic state of the transition cases, so the value of $\xi_{B}$ was not yet stabilized (see figure 2). In fact, Figure 4 shows that even at $\sigma>1$ the time-asymptotic value of $\xi_{B}$ is insensitive to the flow Lorentz factor.

3.3 Precursor spectrum

After discussing the wave energy, we now address the dependence on $\sigma$ of the precursor spectrum and of the typical wavelength of the emission. Our results will be presented in the downstream frame of the simulations. It is important to note, however, that the wave propagates in the upstream plasma and that its emitter is the shock front. Both move with respect to the simulation frame. Hence, when comparing simulation results with the expected scalings, we need to consider the wave dispersion relation first in the upstream frame, and then transform it to the DRF. Also, the typical emission frequency is most naturally estimated in the shock rest frame, and then it should be transformed to the DRF in order to compare with simulation results.

In this section, we first present basic analytical considerations and then we compare them with our simulation results. A special feature of $\sigma>1$ shocks, where a density and magnetic field cavity is observed in the front structure, is discussed at the end of this section. As we argue, the cavity is instrumental in setting the precursor power and determining its dominant frequency.

3.3.1 Basic considerations

As discussed above, the precursor wave possesses X-mode (extraordinary-mode) polarization, such that its wave vector is perpendicular to $\mn@boldsymbol{B}_{0}$, its fluctuating magnetic field is parallel to $\mn@boldsymbol{B}_{0}$, and its fluctuating electric field is perpendicular to both $\mn@boldsymbol{k}$ and $\mn@boldsymbol{B}_{0}$. Some basic properties of the extraordinary mode in the context of the shock emission were derived by Gallant et al. (1992) and Iwamoto et al. (2017). We reproduce here their estimations for completeness.

The dispersion relation of the extraordinary mode in the frame where the background plasma is at rest reads (see, e.g., Hoshino & Arons, 1991)

[TABLE]

where double primed quantities are measured in the upstream rest frame (URF). Using Lorentz transformations for $\omega$ and $k$ and in the limit $\gamma_{0}^{2}\gg\sigma$, the dispersion relation in the DRF becomes

[TABLE]

Interestingly, as long as $\gamma_{0}^{2}\gg\sigma$, this is identical to the dispersion relation of a simple electromagnetic wave propagating in an unmagnetized plasma.

The motion of the shock front imposes a cutoff frequency below which the wave cannot escape into the upstream medium. It follows that little or no power should be observed in the upstream precursor spectrum below the cutoff frequency. This cutoff frequency is obtained by equating the group velocity of the wave, ${\rm d}\omega/{\rm d}k$, with the shock front velocity as:

[TABLE]

This relation leads to the cutoff frequency and wavelength

[TABLE]

As regard to the characteristic frequency of the precursor wave, the most natural assumption is that it corresponds to the collective cyclotron motion of the bunching particles at the shock front, which we now evaluate. First, the magnetic field at the shock can be roughly estimated by assuming that, in the shock frame, all the momentum of the incoming particles is stored in the magnetic field at that point (Alsop & Arons, 1988):

[TABLE]

where primed quantities are measured in the shock rest frame (SRF). More detailed considerations on the soliton structure of the shock as presented by Alsop & Arons (1988), give a similar expression for $B_{\rm sh}^{\prime}$. For particles with Lorentz factors comparable to the upstream bulk Lorentz factor, the ratio of the expected emission frequency (which we label “sol” since it is emitted by the soliton at the shock) to the upstream cyclotron frequency is then equal to the magnetic field enhancement ratio, $\omega_{\rm c,sol}^{\prime}/\omega_{\rm c}=B_{\rm sh}^{\prime}/B_{0}^{\prime}$.333There is no prime on the upstream cyclotron frequency as it is Lorentz-invariant for perpendicular shocks. Lorentz transforming to the DRF ($\omega_{\rm c,sol}^{\prime}\to\omega_{\rm c,sol}$) and using the dispersion relation in Eq. 8 leads to

[TABLE]

Based on these arguments, we expect the precursor spectrum to exhibit a low-frequency cutoff at $\omega_{\rm cutoff}$ and prominent line-like features at $\omega_{\rm c,sol}$ and its harmonics. As we show below, where we compare these scalings with our simulation results, for $\sigma>1$ the predicted $\omega_{\rm c,sol}$ systematically over-estimates the observed peak frequency $\omega_{\rm peak}$. In section 3.3.3, we propose a new model for the precursor peak frequency in the high-magnetization regime, and we show that it is in good agreement with our simulation results.

3.3.2 Spectrum dependence on the upstream magnetization

To characterize the spectrum of the precursor wave, we have employed two complementary diagnostics, one spatial and one temporal. They were used to construct the wavenumber spectrum ($k$-spectrum) and the frequency spectrum ($\omega$-spectrum), respectively.

The wavenumber spectrum was calculated by extracting the spatial profile of $B_{z}(x)-B_{z,0}$ in the region located at $5\,c/\omega_{\rm p}<x-x_{\rm front}<105\,c/\omega_{\rm p}$, at a time when the precursor has reached the steady state, and then computing its Fourier transform. The frequency spectrum was constructed by recording the temporal variation of $B_{z}(t)-B_{z,0}$ at one selected grid point in the upstream region, during a time interval of $100\,\omega_{\rm p}^{-1}$, and then calculating its Fourier transform. The spatial window for the $k$-spectrum and the time interval for the $\omega$-spectrum are chosen so that roughly the same segment of the precursor wave was analyzed in the two cases.

In Figure 5 we present the spectrum of the precursor wave for different $\sigma$. Five representative cases are shown from top to bottom, $\sigma=0.1,0.3,1,3,$ and $10$, respectively. Each panel contains the $k$-spectrum, plotted using blue solid lines, and the $\omega$-spectrum, plotted using red dashed lines. For the $\omega$-spectrum, the horizontal axis shows $\omega/\omega_{\rm p}$, whereas for the $k$-spectrum we take $(k^{2}c^{2}/\omega_{\rm p}^{2}+1)^{1/2}$. Due to this choice, and given the dispersion relation in Eq. 8, each wavenumber spectrum should nearly overlap with the corresponding frequency spectrum, as it is indeed the case. The spectra were normalized such that $\int|\delta\tilde{B}_{z}(k)|^{2}/B_{0}^{2}\,{\rm d}k=\int|\delta\tilde{B}_{z}(\omega)|^{2}/B_{0}^{2}\,{\rm d}\omega=\xi_{B}$.

In each spectrum, the power drops rapidly below the cutoff frequency given by Eq. 10, which is indicated by an orange vertical line in each panel. This is expected, since for lower frequencies (or wavenumbers) the group velocity is smaller than the shock speed, so the wave cannot propagate ahead of the shock.

The spectra are narrow-band, but they are not consistent with a unique line, as it would be expected for cyclotron emission. This is due to the fact that the ring-like particle distribution at the shock front possesses ultra-relativistic energies. The emission is then controlled not by the non-relativistic cyclotron maser, but rather by the ultra-relativistic synchrotron maser instability, that generates a large number of harmonics with comparable growth rate to the fundamental (Hoshino & Arons, 1991).

Prominent line-like features are observed at $\sigma<1$, with the fundamental at $\omega=\omega_{\rm c,sol}$ or the second harmonic dominating the spectrum at low magnetizations (see the peak at $\omega\simeq 4\,\omega_{\rm p}$ for $\sigma=0.1$). In the transition cases with $0.1<\sigma<1$, we observe the generation of very strong harmonics up to $N=5$, where $N=\omega/\omega_{\rm c,sol}$, with high-order harmonics producing stronger lines than the fundamental (see the case with $\sigma=0.3$). For $\sigma>1$ the spectrum shows much less prominent lines. As we will argue later, supplementary amplification mechanisms operate in this regime, and the characteristic frequency $\omega_{\rm c,sol}$ given by Eq. 13 no longer controls the location of the spectral peak.

The dependence of the relevant wavelengths and frequencies on the magnetization is presented in Figure 6. The top row refers to wavelengths, the bottom row to frequencies. The left column shows the variation with $\sigma$ of the cutoff wavelength $\lambda_{\rm cutoff}$ (panel a) and cutoff frequency $\omega_{\rm cutoff}$ (panel b). The values derived from our simulations are plotted using blue circles, and they are in very good agreement with the analytical predictions of Eqs. 10 and 11, indicated by the black dashed lines. The only exception is the transition case $\sigma\simeq 0.3$, where the low-frequency cutoff is non-stationary.

The central column (panels c and d) presents the variation with $\sigma$ of the peak wavelength and frequency (red squares are the results of our simulations), defined as the location where the precursor spectrum peaks (see Fig. 5). Black dashed lines indicate the expectation for soliton emission, Eq. 13. It is apparent that the peak values obtained in the simulations do not agree with the analytical estimate given by Eq. 13 for any $\sigma>0.1$.444We note, however, that in simulations with $\sigma<0.1$, not presented here, we have obtained a very good agreement between the measured $\omega_{\rm peak}$ and $\omega_{\rm c,sol}$ given in Eq. 13 (see also Gallant et al., 1992). To understand the disagreement we define two regimes: (i) the transition cases ($0.1<\sigma<1$) and (ii) the magnetically dominated cases ($\sigma>1$).

In case (i), high-order harmonics in the precursor spectrum are stronger than the fundamental, and the spectral peak is not at the fundamental frequency. If we artificially select the lowest frequency corresponding to a local maximum in the spectrum, we find that its location is in reasonable agreement with the expected fundamental frequency $\omega_{\rm c,sol}$ (see top two panels in Fig. 5). In case (ii), we do not find evidence of any strong line at the expected $\omega_{\rm c,sol}$ or its harmonics, but rather we observe less prominent lines at frequencies that have no clear connection with $\omega_{\rm c,sol}$. The measured peak frequency scales as $\omega_{\rm peak}\approx 3\sigma^{1/2}\omega_{\rm p}$. In contrast, from Eq. 13 we would expect a stronger scaling with $\sigma$, since $\omega_{\rm c,sol}\to\sigma\omega_{\rm p}$ in the limit $\sigma\gg 1$. As discussed in the next subsection, we attribute the observed scaling to the presence of a resonant cavity in the shock structure, that builds up only for $\sigma>1$. We show below that the peak wavelength in case (ii) corresponds to an eigenmode of the cavity, and it is roughly three times shorter than the cavity width (see the green dot-dashed lines in panels c and d).

The right column (panels e and f, respectively) presents the dependence on $\sigma$ of the fractional spectral width in wavelength and frequency space ($\Delta\lambda/\lambda_{\rm peak}$ and $\Delta\omega/\omega_{\rm peak}$, respectively). The width $\Delta\omega$ is the difference between the two frequencies (one above the peak frequency $\omega_{\rm peak}$ and one below) where the power drops by a factor of 30 below the peak. The width $\Delta\lambda$ is defined in an analogous way. This shows quantitatively that the spectrum is narrow, with $\Delta\omega/\omega_{\rm peak}\lesssim 3$ nearly independently of $\sigma$. The spectra of the cases with $\sigma<1$, that show pronounced line-like features, are even narrower, with line widths of $\Delta\omega/\omega_{\rm peak}\lesssim 1$ (see the top two panels in Figure 5).

3.3.3 Resonating cavity in the shock structure at $\sigma>1$

In the previous subsection we have found that in the magnetically dominated regime $\sigma>1$, the peak frequency in our simulations does not scale as the expected gyration frequency in the soliton, $\omega_{\rm c,sol}$. The physical picture that led to the estimate of $\omega_{\rm c,sol}$ must then be revised, since the shock structure for $\sigma>1$ appears to be different than for lower magnetizations. In fact, instead of one density peak defining the shock front, as it is the case in the $\sigma\lesssim 1$ regime, we observe for $\sigma\gtrsim 1$ the build-up of two density peaks separated by a cavity.555We believe that the structure of the shocks studied here is controlled by wave dispersion (rather than dissipation), given the importance of the precursor emission from the shock. For $\sigma>1$, the amount of dispersion provided by the leading soliton becomes insufficient to sustain the shock structure, and a secondary soliton forms to provide additional dispersion. As we now argue, it appears that the density cavity plays an essential role in amplifying the precursor emission and in selecting a well-defined wavelength for the precursor waves that corresponds to an eigenmode of the cavity.

In Figure 7 we illustrate the structure of the shock transition region for $\sigma=10$ at a well-advanced stage of the simulation when the precursor power has reached a steady state. Panel (a) of this figure shows the profile of the electron density (blue line) and of the transverse magnetic field $B_{z}/B_{0}$ (red line). The shock front is located at $x-x_{\rm shock}=0$ and it propagates in the $+\hat{x}$ direction. The two density peaks near the shock are separated by a cavity of width $L_{\rm cav}\simeq 1.6\,c/\omega_{\rm p}$, just behind the shock front. The magnetic field profile peaks at the positions of the two density spikes, but in addition it exhibits a wave-like pattern within the density cavity. For this particular snapshot, only a mode with wavelength $\lambda=L_{\rm cav}/2$ is clearly seen in the cavity. However, the cavity is dynamic in nature, and different eigenmodes are distinctly seen at different times.

In panel (b) of Figure 7 we demonstrate the role of the cavity in shaping the precursor spectrum, by showing the wavenumber spectrum as a function of $\lambda^{-1}=k/(2\pi)$. Some characteristic emission wavelengths are easily identified. For instance, the cutoff wavelength at $\lambda_{\rm cutoff}\simeq 1.6\,c/\omega_{\rm p}$ seems to be closely related to the width of the density cavity $L_{\rm cav}$, which is indicated by a vertical dashed blue line. The other two vertical lines (red and orange, respectively) correspond to wavelengths equal to $L_{\rm cav}/2$ and $L_{\rm cav}/3$, respectively. The latter matches well the position of the strongest emission line. As discussed below, this holds for all $\sigma\gtrsim 1$.

To assess the connection between the cavity size and the precursor efficiency we show in panel (c) the time evolution of $L_{\rm cav}$ (blue line) and of the precursor wave energy $\xi_{B}$ multiplied by a factor of 100 (red line). The value of $\xi_{B}$ was computed in a region closer to the shock front than we have done before (here, between $1$ and $5\,c/\omega_{\rm p}$ ahead of the front), which allows to probe more directly the causal connection between the precursor efficiency and the instantaneous shock structure. This panel shows that the cavity width (blue line) initially increases, then it decreases and finally settles to a steady state. The time evolution of the precursor efficiency appears to be anti-correlated to the cavity width: when the cavity size is larger the emitted precursor is weaker (no amplification), and the wave intensity settles to a steady state at the same time ($\omega_{\rm p}t\sim 1000$) as the cavity width. We interpret this behavior as a self-regulation in the shock structure, such that the cavity width self-tunes to the value where it can efficiently channel the precursor emission into the upstream, i.e., $L_{\rm cav}$ has to be roughly equal to $\lambda_{\rm cutoff}$ (see also panel b). When this condition is met, the wave is amplified and its efficiency settles to the steady state. The critical role of the cavity for efficient wave emission is also revealed by inspecting the shock profile at the time when the precursor intensity sharply increases, right before settling to a steady state ($\omega_{\rm p}t\sim 1000$): we see that large $B_{z}$ fluctuations are first amplified in the cavity, and the emission of a strong precursor propagating upstream is then the consequence of partial transmission of these waves from the cavity through the leading soliton.

The validity of our “resonating cavity” interpretation is tested in Figure 6 (panels c and d), where we show that the peak wavelength of the precursor emission (red squares) is consistent with $L_{\rm cav}/3$ (green dot-dashed lines in panel c), for all $\sigma\gtrsim 1$. In other words, for magnetically dominated plasmas the wave amplification inside the cavity plays an important role in selecting the dominant wavelength of the emitted precursor, as an eigenmode of the cavity. It follows that the peak frequency for $\sigma\gtrsim 1$ scales as $\omega_{\rm peak}\simeq 3\,\omega_{\rm cutoff}\simeq 3\sqrt{\sigma}\omega_{\rm p}\simeq 3\,\omega_{\rm c}$ in the the simulation frame, where we have used that $\gamma_{\rm s|d}\simeq\sqrt{\sigma}$ for $\sigma\gg 1$. This should be contrasted with equation 13, whose scaling ($\propto\sigma\omega_{\rm p}$ in the $\sigma\gg 1$ limit) is not supported by our simulations.

3.4 Wave strength parameter

The wave strength parameter (also known as “wiggler”) measures the dynamical effect of the propagating wave on the background plasma. It is defined through the equation of motion of particles in a high-amplitude wave (Lyubarsky, 2006; Iwamoto et al., 2017):

[TABLE]

where

[TABLE]

is the strength parameter of the wave. $\delta E_{y}$ is the electric field of the wave and $\omega$ is the wave frequency; $e$ and $m_{e}$ are the particle charge and mass, respectively. When $a>1$, the particle quiver motion becomes relativistic and the plasma back-reacts strongly onto the wave. We have employed two measures for the wave strength parameter: either from the maximum excursion in $u_{y}$, $a_{\rm max}=\max(u_{y})$; or from the root mean square value, $a_{\rm std}$. These choices are motivated by the form of equation 15, where the wiggler parameter controls the $y$-oscillations of the particle 4-velocity. In either case, we have extracted the measurement from the region between $5$ and $105\,c/\omega_{\rm p}$ ahead of the front at the final time of the simulations.

Figure 8 presents the dependence on $\sigma$ of the wave strength parameter, derived from our simulations. The maximal value $a_{\rm max}$ decreases from $\simeq 5$ for $\sigma=0.1$ down to $\simeq 1$ for $\sigma=10$, while the root mean square value has the same dependence on $\sigma$ but it is three times smaller, $a_{\rm std}\simeq a_{\rm max}/3$. The sub-panel of this figure shows the dependence on $\gamma_{0}$. Supplementary simulations were performed for this purpose, where we fixed $\sigma=3$. There is a clear linear dependence of $a$ on $\gamma_{0}$, as already suggested by Iwamoto et al. (2017).

The linear dependence on $\gamma_{0}$ arises naturally from the fact that $\xi_{B}$ does not depend on $\gamma_{0}$, combined with the fact that the typical frequency of the precursor wave is $\simeq 3\,\omega_{p}$ for $\sigma\ll 0.1$ and $\simeq 3\,\omega_{\rm c}$ for $\sigma>1$. It follows from Eq. 16 that $a\approx\gamma_{0}\sqrt{\xi_{B}\sigma}$ for $\sigma\ll 0.1$ and $a\approx\gamma_{0}\sqrt{\xi_{B}}/3$ for $\sigma>1$, which justifies the linear scaling with $\gamma_{0}$ shown in the inset of Figure 8.

The wiggler parameter is Lorentz-invariant under transformations along the shock propagation direction, in virtue of equation 15. Alternatively, one can note that the electric field of the wave transforms in the same way as its frequency. Values presented in figure 8 will then be the same in the shock rest frame and in the upstream rest frame. This is in contrast to the precursor normalized energy $\xi_{B}$ and the precursor spectrum, which are frame-dependent.

4 Energetics in the shock front rest frame

The shock front rest frame (SRF) is, by definition, the frame where the shock is stationary. In this frame the upstream plasma flows along the shock normal with a negative velocity in the $x$ direction (whose magnitude is larger than in the DRF). The downstream plasma recedes from the front along the negative $x$ direction. This frame can be naturally employed to quantify the incoming (and outgoing) momentum and energy, and so to derive the energy fraction channeled into the precursor wave.

4.1 From the simulation frame to the shock rest frame

So far, all the quantities related to the precursor waves have been given in the DRF, so we need to Lorentz transform them to the SRF. We will employ primed variables for the SRF. The amplitude of the mean magnetic field transforms as

[TABLE]

where we have used the shortcut notations $B_{0}=B_{z,0}$ and $E_{0}=E_{y,0}$. Since by transforming into the shock frame we are “catching up” with the precursor wave, the precursor amplitude will decrease as

[TABLE]

The $\xi_{B}$ parameter then transforms as (Gallant et al., 1992)

[TABLE]

We compute directly $\xi_{B}^{\prime}$ with the following procedure. The values of $\gamma_{\rm s|d}$ and $\beta_{\rm s|d}$ obtained from our simulations (see figure 3, panels b and c) are used to Lorentz transform the electromagnetic fields into the SRF at a given snapshot of the simulation. Then, $\xi_{B}^{\prime}$ is computed directly, by averaging between $5$ and $105\,c/\omega_{\rm p}$ ahead of the front (the distance is still measured in the simulation frame).

In Figure 9 (panel a) we present the dependence of $\xi_{B}^{\prime}=\xi_{B|\rm sh}$ on $\sigma$, obtained independently with the three PIC codes used in this study: orange squares for Shockapic, red circles for Smilei and blue diamonds for Tristan-MP. First, the figure demonstrates excellent agreement between the three codes. Second, it shows that, beyond the transition cases with $0.1<\sigma<1$, where $\xi_{B}^{\prime}$ attaints the largest values, the normalized wave energy in the SRF scales as $\xi_{B}^{\prime}\simeq 7\times 10^{-4}\sigma^{-2}$ for $\sigma>1$. This scaling is plotted with a dashed black line, and it can be easily justified. In fact, in section 3 we have shown that for $\sigma\gg 1$ the wave amplitude in the DRF converges to a constant (i.e., $\sigma$-independent) value, $\xi_{B}\simeq 10^{-2}$. In addition, the asymptotic shock velocity in the DRF is $\beta_{\rm s|d}\simeq 1-1/(2\sigma)$ for $\sigma\gg 1$. Plugging these two scalings into Eq. 19 leads to $\xi_{B}^{\prime}=6.3\times 10^{-4}\sigma^{-2}$, which is very close to the measured scaling.

4.2 Energy budget in the precursor

Let us also discuss the global energy budget as seen from the SRF (i.e., the fraction of total incoming energy channeled into the precursor). In the SRF, the energy conservation equation including wave contributions can be written as:

[TABLE]

where $W_{i}=w_{i}+b_{0,i}^{2}/(4\pi)$ is the generalized enthalpy, expressed in the proper frame of the fluid. The left-hand side corresponds to the upstream total energy content and the right-hand side to the downstream energy content. We have used that the upstream and downstream electromagnetic wave energies can be expressed as:

[TABLE]

respectively. The brackets represent either space averages at a given time or equivalently time averages at one spatial position. We have neglected the contribution from electrostatic waves, since it is largely sub-dominant in pair plasmas.

Reminding that the upstream magnetization is Lorentz invariant, we use $\sigma=\sigma^{\prime}=b_{0,u}^{2}/(4\pi w_{u})$. The fraction of total incoming energy channeled into the precursor wave is then

[TABLE]

where $\beta_{u}^{\prime}$ is the upstream flow velocity measured in SRF. Equivalently, it is the shock front speed in the upstream rest frame. One can also give the fraction of incoming particle kinetic energy channeled into the precursor wave:

[TABLE]

In the latter equation $n_{u}$ is the upstream plasma proper density, including both species (so, $n_{u}=2N_{0}/\gamma_{0}$).

Getting back to the simulation results, in Figure 9 (panel b) we present the dependence on $\sigma$ of the energy fraction $f_{\xi|\rm sh}=f_{\xi}^{\prime}$, as measured in the SRF. The maximum value is reached at $\sigma\sim 0.1$, where the precursor carries up to 5% of the incoming energy. For $\sigma>0.3$ the energy content in the wave rapidly drops. Similarly to the $\xi_{B}^{\prime}$ scaling, there is a clear dependence as $\propto\sigma^{-2}$ for $\sigma>1$ (more precisely, $f_{\xi}^{\prime}\simeq 7\times 10^{-4}\sigma^{-2}$). The similarity comes from the fact that the upstream velocity is $\beta_{u}^{\prime}\to 1$ and the upstream energy content is dominated by the magnetic field (i.e., $\sigma\gg 1$). This implies from equation 23 that $f_{\xi}^{\prime}\simeq\xi^{\prime}_{B}$.

In the limit $\sigma\gg 1$, the conversion efficiency of incoming particle kinetic energy into wave energy scales as $g_{\xi}^{\prime}\simeq\xi^{\prime}_{B}\sigma\simeq 7\times 10^{-4}\sigma^{-1}$. This should be contrasted with what we have obtained in the DRF, where this quantity became constant in the $\sigma\gg 1$ limit.

We remark that the scalings reported so far have been obtained from 1D runs. While we expect that the dependence on $\sigma$ will remain unchanged in 2D and 3D, we speculate that the normalizations of $f_{\xi}$ and $g_{\xi}$ will decrease due to transverse effects that cannot be captured in 1D, e.g., wave filamentation and self-focusing through interaction with the upstream plasma. In fact, the 2D simulations of Iwamoto et al. (2017, 2018), performed in the low magnetization regime $\sigma<0.5$, demonstrated that the wave energy is reduced typically by a factor of 3 (and up to 10), when going from 1D to 2D. However, we expect that the efficiency drop from 1D to 2D (and 3D) will be much less severe in the magnetically-dominated regime ($\sigma>1$) of interest for our work, given the rapid decrease of the wave strength parameter with magnetization (see Figure 8), and so of the wave feedback onto the upstream plasma. This point will be addressed in a forthcoming study (Sironi et al, in prep.).

5 Applications to FRBs

During a magnetar flare, in response to the motions of the neutron star crust, the above-lying magnetosphere is violently twisted and a strongly magnetized pulse is formed, which propagates away through the magnetar wind. The FRB can be potentially generated at ultra-relativistic shocks resulting from the collision of the magnetized pulse with the steady magnetar wind produced by its spin-down luminosity or by the cumulative effect of earlier flares (Lyubarsky, 2014; Beloborodov, 2017; Waxman, 2017). The train of electromagnetic waves emitted by the shock front via the synchrotron maser is the candidate FRB. Most works up to now assumed empirical values for the conversion efficiency of the shock kinetic energy into the precursor waves. These values were primarily motivated by the work of Gallant et al. (1992) where, however, high-$\sigma$ simulations were not evolved long enough to reach a stationary state. Here, we use long-term simulations to quantify the steady-state energetics and spectrum of the precursor waves, for a wide range of magnetizations (up to $\sigma\gg 1$). As we now argue, our work can provide a physically-grounded model for the origin of coherent emission in FRBs.

First, the synchrotron maser at shocks is a coherent process, which helps explaining the extremely high brightness temperatures of FRBs. In this work, we have derived the fraction of incoming flow energy channeled into the precursor waves. If considered in the ejecta frame (post-shock frame), our simulations show that for $\sigma>1$ the emitted wave carries a fraction $f_{\xi}=2\times 10^{-3}/\sigma$ of the total energy. This corresponds to a fraction $g_{\xi}\simeq 2\times 10^{-3}$ of the incoming particle kinetic energy, regardless of $\sigma$. If one considers the energy budget in the shock rest frame, the previous scalings become $f_{\xi}^{\prime}\simeq 7\times 10^{-4}/\sigma^{2}$ and $g_{\xi}^{\prime}\simeq 7\times 10^{-4}/\sigma$, respectively.

Second, the precursor emission is linearly polarized, in agreement with the observations of several non-repeating FRBs (Ravi et al., 2016; Petroff et al., 2017; Caleb et al., 2018) and of the repeating FRB 121102 (Michilli et al., 2018; Gajjar et al., 2018). Linear polarization is a natural consequence of the resonance of bunching particles with the extraordinary mode (X-mode). This mode can escape out of the plasma and become a vacuum electromagnetic wave. A contribution from the ordinary mode (O-mode) was also observed in the 2D simulations of Iwamoto et al. (2018), but it was found to be largely sub-dominant in strongly magnetized plasmas.

Third, the spectral peak can fall in the GHz range for a reasonable choice of parameters. In particular, we have found that in the post-shock frame the emission peak frequency scales as $\omega_{\rm peak}\simeq 3\,\omega_{\rm p}$ for $\sigma\lesssim 0.1$ and as $\omega_{\rm peak}\simeq 3\,\omega_{\rm c}$ for $\sigma>1$. Several high-order harmonics characterize the transition region with $0.1<\sigma<1$. Joining the two regimes, and neglecting for simplicity the transition cases, we can cast the peak frequency as $\omega_{\rm peak}\simeq 3\,\omega_{\rm p}\max[1,\sqrt{\sigma}]$. This can be recast in a simpler form in the pre-shock frame as

[TABLE]

as long as the shock is moving with an ultra-relativistic bulk Lorentz factor $\gamma_{\rm s|u}$ into the upstream medium. The emission frequency for an upstream observer is then

[TABLE]

where $n_{e}$ is the pre-shock electron density. If we assume that the upstream frame corresponds to the observer frame (which is true if the pre-burst wind expands with a non-relativistic velocity), then the combination $\gamma_{\rm s|u}\sqrt{n_{e}/1~{}{\rm cm}^{-3}}\approx 4\times 10^{4}$ is required for the shock to emit in the GHz band, in rather good agreement with the estimates of Beloborodov (2017). As recently found by Metzger et al. (2019) this frequency is also consistent with $\sim$GHz emission from decelerating blast waves produced by flare ejecta in young magnetars.

Finally, the spectrum is narrow-band, $\Delta\omega/\omega_{\rm peak}\lesssim 1-3$ (see figure 6), which is again consistent with the observations (e.g., Law et al., 2017; Macquart et al., 2018).

5.1 Comment on criticisms to the synchrotron maser

A number of criticisms have recently been moved against the synchrotron maser emission as a source of the coherent FRB radiation. Lu & Kumar (2018) looked into a wide variety of maser mechanisms operating in either vacuum or plasma and found that none of them can explain the high luminosity of FRBs without invoking unrealistic or fine-tuned plasma conditions. Here, we argue that the synchrotron maser at relativistic shocks — due to its unique properties — still remains a viable candidate for powering FRBs.

First, it was argued that the synchrotron maser in vacuum requires fine-tuned plasma conditions where the magnetic field is nearly uniform (to within an angle $\gamma^{-1}$) and the particles’ pitch-angle distribution is narrowly peaked with spread $\lesssim\gamma^{-1}$. Here, $\gamma$ is the typical Lorentz factor of the emitting particles. This is indeed the natural configuration expected at a relativistic magnetized shock, if the pre-shock particles have non-relativistic temperatures (which is anyway a requirement for efficient synchrotron maser emission). In the shock transition region, the magnetic field is nearly uniform, and the particles coherently rotate in a plane perpendicular to the field (with negligible pitch angle spread).

Second, it was argued that it is unclear how the mechanism for the population inversion required by the maser is achieved. Once again, this is naturally realized in the shock transition of a magnetized relativistic shock, where the particles form a ring in momentum space at fixed Lorentz factor $\gamma\sim\gamma_{0}$, while the inner region of the ring (i.e., at lower $\gamma$) is devoid of particles, as indeed required for the existence of a population inversion.

Also, it was argued that during the maser amplification process, high-energy electrons radiate faster than low-energy ones, so the population inversion condition may be quickly destroyed. This is indeed true for each generation of particles passing through the shock, since the synchrotron maser instability relaxes by “filling up” the hollow ring in momentum space, thus destroying the population inversion. However, while this happens, a new generation of particles is entering into the shock. They establish a new ring in momentum space, and keep sustaining the radiated train of precursor waves. In other words, the continuous passage of plasma through the shock ensures that the population inversion is steadily maintained (yet, at each time by different particles).

Finally, Lu & Kumar (2018) considered more specifically the maser synchrotron emission at shocks, which they named as “bunching in the gyration phase.” In order to minimize the effect of induced Compton scattering, they estimated that the radiative efficiency of the shock must be extremely small. However, they considered only internal shocks occurring in between two identical consecutive density shells propagating inside the pre-burst wind, and not the leading shock moving directly into the wind. Aside from the limitations of induced Compton scattering, it is anyway hard for internal shocks to be efficient emitters of maser synchrotron radiation, since they propagate into a relativistically hot shocked plasma (the downstream region of the leading shock). The arguments by Lu & Kumar (2018) will not apply to the leading shock. First, this shock is likely to be ultra-relativistic, unlike internal shocks. Second, the properties of the shell and of the pre-burst wind (as regard to magnetization, temperature, composition) are generally different, in contrast to what Lu & Kumar (2018) implicitly assumed. We believe that the quantitative results on precursor energetics and spectrum that we provide in this work will help revisit the estimates provided by Lu & Kumar (2018), for the case of the leading shock.

6 Summary and conclusions

In this work we have investigated by means of 1D Particle-In-Cell simulations the physics of synchrotron maser emission from perpendicular relativistic shocks that propagate in highly magnetized electron-positron plasmas (with magnetization $0.1\leq\sigma\leq 30$). For strongly magnetized shocks, we expect that multi-dimensional simulations (to be discussed in a forthcoming work) will not yield very different results than what we present here. We have explored the efficiency and spectrum of the electromagnetic precursor emission as a function of $\sigma$ and $\gamma_{0}$. We have found that:

The shock front emits efficiently and steadily a train of high-amplitude electromagnetic precursor waves for any $\sigma$ and $\gamma_{0}$, in the range $0.1\leq\sigma\leq 30$ and $\gamma_{0}\geq 5$ that we have explored. The emission is linearly polarized, with fluctuating magnetic field along the same direction as the upstream mean field. 2. 2.

Thanks to unprecedentedly long simulations, we have been able to reach the stage when the precursor emission settles to a steady state, which allows to systematically extract the wave properties (energetics and spectrum). We find that the ratio of the wave energy to the upstream magnetic energy, $\xi_{B}$, decreases rapidly from $3.5$ at $\sigma=0.1$ down to $0.04$ at $\sigma=1$, as measured in the post-shock frame of the simulations. For $\sigma\gg 1$, this ratio converges to a constant value $\xi_{B}\simeq 0.01$. In the shock rest frame, the asymptotic scaling in the limit $\sigma\gg 1$ becomes $\xi_{B}^{\prime}\propto\sigma^{-2}$. 3. 3.

For $\sigma>1$, the energy output in precursor waves normalized to the total incoming energy scales as $f_{\xi}\simeq 2\times 10^{-3}\sigma^{-1}$ in the post-shock frame and as $f_{\xi}^{\prime}\simeq 7\times 10^{-4}\sigma^{-2}$ in the shock rest frame. The former implies that in the downstream frame, $\sigma>1$ shocks convert a constant fraction of the incoming particle kinetic energy into precursor waves (equal to $g_{\xi}\simeq 2\times 10^{-3}$). 4. 4.

Magnetically dominated shocks with $\sigma>1$ exhibit a resonating cavity in the shock front structure in between two solitons, instead of the single soliton loop that is observed for $\sigma\ll 1$ shocks. This cavity plays an essential role in amplifying the radiation and selecting the dominant emission frequency as an eigenmode of the cavity. This effect causes the peak emission frequency, as measured in the downstream frame, to scale as $\omega_{\rm peak}\simeq 3\,\omega_{\rm c}=3\sqrt{\sigma}\omega_{\rm p}$ for $\sigma>1$, whereas earlier works (Gallant et al., 1992) quote a stronger scaling with magnetization, $\omega_{\rm peak}\simeq\sigma\omega_{\rm p}$. 5. 5.

The characteristic frequency of the emission, as measured in the post-shock frame, is $\omega\simeq 3\omega_{\rm p}$ for weakly magnetized shocks $\sigma\leq 0.1$, and $\omega\simeq 3\omega_{\rm c}$ for $\sigma\gg 1$, as we have just discussed. In the transition region $0.1<\sigma<1$, prominent high-order harmonics of $\omega_{\rm c,sol}$ (given in Eq. 13) were observed along with the fundametal at $\omega_{\rm c,sol}$. Aside from the transition cases, we can interpolate between the low- and high-magnetization results and state that the peak emission occurs at $\omega_{\rm peak}\simeq 3\,\omega_{\rm p}\max[1,\sqrt{\sigma}]$, as measured in the downstream frame. In the pre-shock frame (which coincides with the observer frame, if the magnetar wind is non-relativistic), this can be recast in a simpler form as $\omega^{\prime\prime}_{\rm peak}\approx 3\gamma_{\rm s|u}\omega_{\rm p}$, where $\gamma_{\rm s|u}$ is the shock Lorentz factor in the upstream frame. 6. 6.

The spectrum of the precursor is narrow-band, $\Delta\omega/\omega_{\rm peak}\lesssim 1-3$, with a low-frequency cutoff at $\omega_{\rm cutoff}=\gamma_{\rm s|d}\omega_{\rm p}$ (here, $\gamma_{\rm s|d}$ is the shock Lorentz factor in the downstream frame) set by the requirement that the group velocity be faster than the shock speed. 7. 7.

We did not observe any dependence on $\gamma_{0}$ of the energy fraction, $\xi_{B}$, and of the characteristic emission frequency, $\omega_{\rm peak}/\omega_{\rm p}$, in the post-shock frame.

We conclude with a few caveats. First, we have assumed that the upstream plasma has negligible thermal spread, $k_{B}T_{0}/m_{e}c^{2}=10^{-4}$. Higher temperatures are likely to suppress high-order harmonics and reduce the global energy of the wave. Second, we have mostly focused on strongly magnetized ($\sigma>1$) plasmas, a regime that so far has received little attention. Even though this work only presents 1D simulations, we anticipate that the multi-dimensional physics of $\sigma>1$ shocks (Sironi et al., in prep.) will not depart significantly from what we report here. In contrast, for weaker magnetizations ($\sigma\ll 10^{-2}$), transverse effects (e.g., Weibel-driven filamentation) will significantly reduce the energy carried by the precursor waves (Sironi et al., 2013; Iwamoto et al., 2017). In summary, both higher pre-shock temperatures and multi-dimensional effects at low $\sigma$ are expected to degrade the precursor efficiency, which might become too low to explain the FRB emission.

Finally, we have only considered electron-positron shocks. Recently, a very large Faraday Rotation Measure (RM) of $\sim 10^{5}$ rad m*-2* was reported from the repeating FRB 121102 (Michilli et al., 2018). This challenges the pure electron-positron composition assumed in this study, since the presence of an appreciable fraction of ions is required to produce non-zero RM (Margalit & Metzger, 2018). This urges to explore the shock physics for electron-proton and electron-positron-proton compositions. Yet, it is still possible that the FRB pulse is produced in localized regions with pristine electron-positron composition, even though most of the magnetar wind (which inflates the surrounding nebula, where the RM accumulates) is proton-dominated.

Acknowledgments

IP acknowledges discussions with Anatoly Spitkovsky, Patrick Crumley and Yuri Cavecchi. LS is grateful to Brian Metzger for many inspiring discussions. IP was supported by NSF grants PHY-1804048 and PHY-1523261. This work was facilitated by the Max-Planck/Princeton Center for Plasma Physics. LS acknowledges support from NASA ATP 80NSSC18K1104. The simulations were performed on Habanero cluster at Columbia University, NERSC (Edison) and NASA (Pleiades) resources, PICSciE-OIT High Performance Computing Center and Visualization Laboratory at Princeton University, and on CALMIP supercomputing resources at Université de Toulouse (France) under the allocation 2016-p1504.

Appendix A Codes comparison

In this appendix we show how the results from the different codes compare. The synchrotron maser emission occurs through the resonance of the cyclotron harmonics with the X-mode branch. It is not guaranteed that a typical PIC code can capture accurately a large number of these resonances, especially at the high-frequency end of the branch. For instance, in typical Yee-type second order solvers of Maxwell’s equations the numerical speed of electromagnetic waves is known to be artificially suppressed at high $\omega$ if the CFL number is smaller than unity (Birdsall & Langdon, 1991).

For this reason we undertook an extensive comparison of three different PIC codes: two Finite-Difference Time-Domain (FDTD) codes, Tristan-MP and Smilei, and the pseudo-spectral code Shockapic. In principle, Shockapic is the best suited to capture the dispersion relation of waves in a plasma, but it is the least optimized among the three codes, making it challenging to perform long-term simulations. Concerning Smilei, it is a well-optimized code, but in 1D setups it currently has only a standard Yee solver. Tristan-MP is the most optimized for shock setups. Also, it allows to use a fourth-order scheme to solve Maxwell’s equations (Greenwood et al., 2004) that reproduces accurately the dispersion relation of electromagnetic waves even at low CFL numbers. This is the reason why the runs presented in the main body of the paper were performed with Tristan-MP. In general, each code employs different algorithms and implementations. The agreement between the three codes will then be a strong indication of the physical robustness of our results.

The simulation parameters for each code are presented in Table 1. The table reports the space and time resolution, the simulation timespan, the number of particles per cell, the values of the upstream temperature $T_{0}$ and bulk Lorentz factor $\gamma_{0}$, and the explored range of $\sigma$. For better comparison we used comparable space and time resolutions: the skin depth was resolved with 100 cells in Tristan-MP simulations, with 112 cells in Smilei simulations, and with 44.7 cells in Shockapic simulations. The latter has a twice smaller resolution due to code performance limitations (not parallelized). We noticed that a resolution lower than 20 cells per skin depth affected negatively the results for any $\sigma$. The results become stable for any resolution higher than 40 cells per $c/\omega_{\rm p}$, as long as $\sigma\leq 10$. Similar conclusions were reached by Iwamoto et al. (2017). For this reason a high spatial and time resolution was employed in the simulations presented in the main body of the paper. Only short simulations were affordable with Shockapic. For this reason, the $\sigma>1$ regime was not explored with this code (as we have discussed, at high $\sigma$ it takes longer to reach a steady state). With Tristan-MP and Smilei it was possible to reach the stationary state for $\sigma$ up to 30. Concerning the number of particles per cell, the results are very weakly dependent on $N_{\rm ppc}$, as long as at least a dozen of particles per cell are initialized.

In the following we present in more detail the comparison of precursor energy and spectrum as derived from different codes.

A.1 Precursor energy

In figure 10 we present the normalized wave energy $\xi_{B}$ as a function of $\sigma$, obtained with the three codes. Values obtained with Shockapic, Smilei, and Tristan-MP are plotted using orange squares, red circles, and blue diamonds, respectively. In the overlapping range of $\sigma$, we observe good agreement among different codes. For instance, in the $\sigma>1$ regime Smilei and Tristan-MP give the same values of $\xi_{B}$. In the range $0.1<\sigma\leq 1$, where all codes overlap, the scatter among codes is slightly larger, although the rapid drop in $\xi_{B}$ is common to all codes, and it happens around the same $\sigma$. We note that the transition is more abrupt in Shockapic than in Smilei and Tristan-MP, but differences remain minor.

In figure 11 we extend the comparison to different studies in the literature and to different values of $\gamma_{0}$, from $10$ to $10^{6}$. The range of $\sigma$ in this figure is from $10^{-3}$ to $1$, since other studies did not explore highly magnetized cases with sufficiently long simulations (i.e., they did not reach a steady state in the regime $\sigma\gg 1$). Blue circles and green diamonds report the values obtained with Smilei using $\gamma_{0}=10$ and $160$, respectively. Both give nearly the same values for any explored $\sigma$, confirming that $\xi_{B}$ does not depend on the flow Lorentz factor. Red squares report the values from Tristan-MP using $\gamma_{0}=10$ (same as in figure 10). The data from the 1D simulations of Iwamoto et al. (2017) using $\gamma_{0}=40$ are plotted with orange stars. Their values are slightly smaller than what is found in this study, though generally in good agreement. Violet stars and light-blue triangles report the values from Gallant et al. (1992) using $\gamma_{0}=40$ and $10^{6}$, respectively. We notice that all codes provide the same results in the range of $\sigma\in[10^{-3},0.3]$, regardless of $\gamma_{0}$. This demonstrates that the precursor wave normalized energy $\xi_{B}$ is not dependent on $\gamma_{0}$, and that our study is in very good agreement with earlier results.

For $\sigma>0.3$ there is a noticeable scatter between different simulations. The most plausible reason for the discrepancy among different datasets is that the high-$\sigma$ simulations from earlier studies were not evolved long enough to reach the asymptotic state, so the value of $\xi_{B}$ was not yet stabilized (see, Fig. 2 for the time convergence of the efficiency).

A.2 Precursor spectrum

We now compare the precursor $k$-spectrum among the three codes. Some differences are expected, since the numerical schemes for the integration of Maxwell’s equations differ among the codes.

In figure 13 we compare the precursor spectrum extracted from the three codes for a few representative values of magnetization. From top to bottom, the value of $\sigma$ is $0.1$, $0.5$ and $1$, respectively. We cannot perform any comparison for $\sigma>1$ as this range was not explored with Shockapic (but see below, for a comparison between Smilei and Tristan-MP at $\sigma=30$). The spectrum extracted from Tristan-MP is plotted using a solid blue line. Red and orange lines are used for Smilei and Shockapic, respectively. There is generally a good agreement among the codes for all values of $\sigma$ as regard to the low-$k$ cutoff wavenumber, the high-$k$ slope, and the main peaks in the spectrum. For example, the dominant emission line for $\sigma=0.1$ and the high-order harmonic line at $\lambda=0.24\,c/\omega_{\rm p}$ for $\sigma=0.5$ are exactly at the same wavelength for the three codes. One difference can be noted: the spectral energy density is slightly smaller in Shockapic than in the two FDTD codes around $\lambda^{-1}c/\omega_{\rm p}\sim 1$, for $\sigma=0.5$ and $\sigma=1$. Yet, this difference is not systematic and the overall energy in the precursor is very close among the three codes.

As an exception and a word of caution, we noticed that the use of a small CFL number with a Yee-type solver of Maxwell’s equations (as used in the Smilei code) has a negative impact on the results for the largest magnetizations explored here, i.e., $\sigma>10$. In fact, the emission peaks at high frequencies where the light-wave branch is affected by the artificial reduction of the phase speed. The spectrum of the precursor is then sharply cut at high frequencies, affecting the overall energy output in the precursor. This effect is evidenced in figure 13 for $\sigma=30$ (the largest value explored in this work). The upper panel of the figure compares the spectrum from Tristan-MP (blue), where a fourth-order scheme was used, with the spectrum from Smilei (red), which employs a Yee-type scheme with $c\Delta t/\Delta x=0.5$. There is an artificial suppression in the high-$k$ region in the Smilei simulation. The bottom panel shows the same comparison, but with $c\Delta t/\Delta x=0.99$ being used with Smilei. In this case, the spectra agree very well, up to details in line-like features. This conveys that the high-$k$ (and so, high-$\omega$) part of the precursor spectrum can be properly captured only when the numerical integrator is capable of reproducing correctly the dispersion relation of electromagnetic waves. This problem does not arise in Tristan-MP (with high-order spatial solver) and Shockapic, since for them the numerical dispersion of the light-wave branch is much closer to the realistic one even for small CFL numbers.

All Smilei simulations that use a CFL number as close as possible to unity (CFL=0.99) display spectra that are in very good agreement with the other two codes for any $\sigma$.

A.3 Concluding remark

We find that our results do not depend on the code that we employ if these three conditions are realized: (i) a high spatial resolution (i.e., large $c/\omega_{\rm p}$) is employed; (ii) in a Yee-type based code, the CFL number is as close as possible to unity; (iii) the simulations are sufficiently long to reach the steady state.

Appendix B Precursor energetics: 1D vs multi-dimensional simulations

In order to support our claim that the precursor wave energy does not significantly decrease due to multi-dimensional effects (in the $\sigma\geq 1$ regime of interest for this work), here we present a preliminary analysis of 2D and 3D simulations performed with Tristan-MP. 666The code accuracy and stability in multi-dimensional simulations of relativistic shocks was assessed in several studies (Spitkovsky, 2005, 2008; Sironi & Spitkovsky, 2009, 2011; Sironi et al., 2013). We explore a range of $\sigma\in[0.1,10]$ in 2D, and $\sigma\in[0.1,3]$ in 3D. In 2D simulations we focus on the out-of-plane configuration: the simulation plane is the xy plane, the shock front propagates in the $x$-direction, and the upstream magnetic field is along the $z$-direction. We do not present any in-plane 2D simulation results here because we find that 3D simulations are in excellent agreement with 2D out-of-plane results.

In 2D simulations we keep all parameters the same as in 1D, except that the number of particles per cell per species is set to 8 (values between 2 and 32 have been tested with no significant differences). The transverse dimension of the simulation box is set to $14\,c/\omega_{\rm p}$.We find that a transverse width of more than $2-3\,c/\omega_{\rm p}$ is sufficient to capture multi-dimensional effects. In particular, the effects of wave filametation and self-focusing that lead to efficient pre-heating of the upstream plasma in the longitudinal momentum are properly captured with a box width of a few skin depths.

In 3D simulations we reduce the transverse dimension to $4\,c/\omega_{\rm p}$ (in both $y$ and $z$ directions). The spatial resolution in 3D runs is set to 25 cells per $c/\omega_{\rm p}$ (four times lower than in 1D and 2D) and the number of particles per cell per species is varied between 3 and 18 (again, with little differences). This was necessary to produce sufficiently long runs while still capturing the relevant physics. The effect of a lower spatial resolution was only apparent in the $\sigma=3$ run, since the spectrum extends to higher frequencies, which are not captured properly if the resolution is insufficient.

Figure 14 presents the time evolution of the precursor energy for $\sigma=1$ as measured in 1D (blue line), 2D (red line) and 3D (orange line) simulations. The direct comparison between 1D and multi-dimensional simulations shows that the asymptotic value of $\xi_{B}$ in 2D and 3D simulations is only a factor of two smaller than in 1D, while for $\sigma<0.5$ — as we will show below, and see also Iwamoto et al. (2017) — the energy of the wave decreases by a factor of about $3-10$ when going from 1D to 2D and 3D configurations. It also shows that the 2D and 3D energetics are in very good agreement. The only difference between 2D and 3D is that it takes more time in 3D to settle into the steady state (see the rise of the red line after $t\omega_{\rm p}=500$ and of the orange line after $t\omega_{\rm p}=1000$). So, we can confidently state that the decrease in precursor efficiency due to multi-dimensional effects is much less severe in the high-magnetization case $\sigma=1$ than for $\sigma<0.5$.

Let us note that in this appendix we have redefined the $\xi_{B}$ parameter. Here, $\xi_{B}$ corresponds to the normalized Poynting flux in the $x$-direction, $\xi_{B}=\langle\delta E_{y}\delta B_{z}-\delta E_{z}\delta B_{y}\rangle/B_{0}^{2}$. The average is done over the region between $5$ and $25\,c/\omega_{\rm p}$ ahead of the shock front, for consistency with our 1D results, and over all the transverse directions ($y$ in 2D; $y$ and $z$ in 3D). In 1D we have systematically verified that $\langle\delta B_{z}^{2}\rangle=\langle\delta E_{y}\delta B_{z}-\delta E_{z}\delta B_{y}\rangle=\langle\delta E_{y}\delta B_{z}\rangle$, but this equality is not obviously satisfied in multi-dimensional simulations with $\sigma\leq 0.6$. The choice of defining $\xi_{B}$ as the precursor Poynting flux is due to the fact that the most relevant measure of the electromagnetic energy output of the shock is the Poynting flux of the escaping wave in the shock-normal direction.

In Figure 15, using a suite of 1D, 2D and 3D simulations, we show the dependence on $\sigma$ of the normalized Poynting flux of the precursor wave $\xi_{B}$ (panel a), of the energy fraction parameter as measured in the simulation frame $f_{\xi}$ (panel b), and of the energy fraction parameter as measured in the shock rest frame $f_{\xi}^{\prime}=f_{\rm\xi|sh}$ (panel c). The definition of the latter two is given in the main body of the article: Eq. 5 and Eq. 23, respectively. Values from 1D, 2D and 3D simulations are plotted using blue circles, red squares and green stars, respectively. The results of 2D out-of-plane simulations of Iwamoto et al. (2017) are plotted using orange triangles in panel (a). The measurement of $\xi_{B}$ in 2D and 3D simulations was done by considering the asymptotic values in the time evolution for each $\sigma$, as shown in Figure 14 for the particular case of $\sigma=1$. Error bars quantify uncertainties due to temporal oscillations of the time-evolution curves. Knowing $\xi_{B}$ and measuring directly the shock front velocities from simulations, the values in panels (b) and (c) were produced using Eq. 5 and Eq. 23, respectively.

Figure 15 shows that:

• •

In 2D and 3D (red and green symbols), for $\sigma=0.1$ the Poynting flux of the precursor wave $\xi_{B}$, the energy fractions $f_{\xi}$ and $f_{\xi|\rm sh}$ are reduced by a factor of $\approx 10-20$ as compared to 1D (blue circles). This is in agreement with Iwamoto et al. (2017).

• •

The suppression in efficiency becomes gradually smaller when $\sigma$ increases from $0.1$ to $3$. For $\sigma\gtrsim 1$, the difference between 1D and multi-dimensional results becomes negligible.

• •

Values from 2D out-of-plane and 3D simulations are generally in very good agreement, except for $\sigma=0.3$ and $0.4$ (which we have called “transition cases” in the main body of the text).

• •

If the precursor energy fraction is cast in the shock rest frame, panel (c) shows that $f_{\xi}^{\prime}\simeq 10^{-3}$ for $\sigma\sim 0.1-0.4$, instead of $\sim 0.01$ in 1D. For $\sigma>1$, multi-dimensional simulations converge towards 1D values and follow the scaling $f_{\xi}^{\prime}\approx 5\times 10^{-4}/\sigma^{2}$, only slightly lower than reported in the main text for 1D simulations only.

Using 3D simulations we can address other aspects of the precursor physics, such as the importance of the O-mode ($\delta B_{y}$ component, since $\mathbf{\delta B\perp B_{0}}$ for this mode) versus X-mode ($\delta B_{z}$ component, since $\mathbf{\delta B\parallel B_{0}}$ for this mode) and beaming of the emitted precursor wave. By extracting systematically the values of $\langle\delta B_{y}^{2}\rangle$ and $\langle\delta B_{z}^{2}\rangle$ in 3D simulations, we find that the O-mode is subdominant for all magnetizations explored here, i.e., $\langle\delta B_{y}^{2}\rangle/\langle\delta B_{z}^{2}\rangle\sim 10^{-3}$. This implies that the precursor wave retains (at the $99\%$ level, or more) the linear polarization of the X-mode, with magnetic field of the wave lying in the same direction as the upstream background field.

Concerning the beaming of the precursor wave in 3D, we considered the components of the Poynting vector in different directions. We find that the Poynting flux along the $y$-direction (and $z$-direction) is largely subdominant as compared to the shock-normal direction. The ratio is $|\Pi_{y}|/\Pi_{x}\sim 5\times 10^{-4}$ for any $\sigma\in[0.1,3]$, where the Poynting vector of the wave is defined as $\mathbf{\Pi}=\mathbf{\delta E}\times\mathbf{\delta B}/B_{0}^{2}$. This shows that the emitted wave is strongly beamed in the shock-normal direction. For an external observer the beaming will be further enhanced by Lorentz transformation from the simulation frame to the observer frame (in the case of shocks in magnetar winds, from the post-shock frame to the pre-shock frame).

This preliminary analysis of multi-dimensional runs demonstrates that 1D simulations provide accurate numbers in the $\sigma\gg 1$ regime, in agreement with 2D out-of-plane and 3D simulations.