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

I N T RO D U C T I O N

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 (Treumann 2006;Melrose 2017). 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 semicoherent 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, '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, Lu & Bhattacharya 2017;Ghisellini & Locatelli 2018;Katz 2018;Lu & Kumar 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 ultrarelativistic 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 1D simulations (more precisely, 1D3V, i.e. we employ one spatial dimension, but all three components of velocities and electromagnetic fields are retained), while multidimensional runs will be presented in a future work (Sironi et al., in preparation; 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 modelling of the electromagnetic precursor wave in relativistic perpendicular shocks (e.g. Langdon, Arons & Max 1988;Gallant et al. 1992;Hoshino et al. 1992;Amato & Arons 2006;Hoshino 2008;Sironi & Spitkovsky 2009;Iwamoto et al. 2017Iwamoto et al. , 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 σ 1, where σ is the ratio of upstream Poynting flux to kinetic energy flux. For σ 1, 1D simulations are adequate, since we find that they agree well with multidimensional results (see Appendix B, for the precursor energetics in 2D and 3D; also Sironi et al., in preparation). We provide an extensive investigation of the dependence on the flow magnetization, from σ = 0.1 to σ = 30, with much longer simulations than previously reported, especially in the σ 1 regime relevant for FRB sources. Previous works arguably never reached a steady state in simulations with σ 1. We employ several PIC codes (TRISTAN-MP, SMILEI and SHOCKAPIC) to check for consistency, and thus confirm the robustness of our results.

At σ 1 the shock converts a fraction f ξ ≈ 2 × 10 -3 /σ of the total incoming energy into the precursor wave, as measured in the post-shock (downstream) frame. In the shock rest frame (SRF), the efficiency is f ξ ≈ 7 × 10 -4 /σ 2 . The spectrum of the precursor wave is narrow-band, ω/ω peak 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 ω peak 3 ω p max[1,

√ σ ], where ω 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 ω peak ≈ 3γ s|u ω p , where γ s|u is the shock Lorentz factor in the upstream frame. At σ 1, where our estimated ω peak differs from earlier works [that quoted ω peak ∝ σ ω p , see Gallant et al. (1992), rather than ω peak ∝ √ σ ω 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 set-up. 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.

S I M U L AT I O N M E T H O D S A N D S E T-U P

We use the PIC codes TRISTAN-MP (Spitkovsky 2005;Sironi & Spitkovsky 2009), SMILEI (Derouillat et al. 2018), and SHOCKAPIC (Plotnikov, Grassi & Grech 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 σ 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 preparation). For lower magnetizations than explored here, i.e. σ 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 set-up described in, e.g. Spitkovsky (2008), which we summarize here for completeness. The upstream flow, composed of electrons and positrons, drifts along thex direction with a speed -β 0 c x. The corresponding bulk Lorentz factor is γ 0 = (1β 2 0 ) -1/2 = 10, but we have also explored higher values of γ 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 = -β 0 B z, 0 . The flow is reflected at a wall located at x = 0. After some time (at least several cyclotron periods ω -1 c ), the shock front forms by magnetic reflection and steadily propagates along the + 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 σ , which we define as the ratio of Poynting to kinetic energy flux

and we vary from σ = 0.1 up to σ = 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 gyrofrequency as ω c = |q|B 0 /(γ 0 m e c) and the plasma frequency as ω p = [8 N 0 q 2 /(γ 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 (i) The skin depth is well resolved with c/ω p = 100 , where is the grid size. This ensures that the typical particle gyroradius σ -1/2 c/ω p is well resolved even for the largest magnetization σ = 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).

(ii) The simulation time-step is defined as c t = 0.5 corresponding to a time resolution of 5 × 10 -3 ω -1 p . (iii) 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).

(iv) The simulation is evolved up to T sim 1.5 × 10 3 ω -1 p for σ < 1 and for longer times (up to 2 × 10 4 ω -1 p ) for σ 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.

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. 2017Iwamoto et al. , 2018;;Plotnikov et al. 2018), and 3D simulations (Spitkovsky 2005;Sironi, Spitkovsky & Arons 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, σ 0.01), the transverse Weibel filamentation instability, which dominates for σ 10 -3 , plays no significant role in shaping the shock structure. This partly justifies the 1D approach adopted here.

In Fig. 1 we present the structure of the shock transition region for two representative magnetizations. The left-hand column (panels ae) presents a shock with σ = 0.3 and the right-hand column (panels f-j) corresponds to σ = 3. The timespan of the simulations is long enough to reach a stationary state: we show results at ω p t = 540 for σ = 0.3 and at ω p t = 1800 for σ = 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 xu x (d and i), and the transverse positron phase space xu y (e and j). Here, we define u α = γ β α 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-hand side of panels (d) and (i) show the particle distribution in momentum space u xu y at the location of the shock front.

The shock front is located at xx front = 0 in Fig. 1. The upstream flow is on the positive side (xx front > 0) and the downstream plasma is on the negative side (xx 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 semicoherent 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 xx front > 0 region). The precursor wave amplitude is larger for σ = 0.3 than for σ = 3. This wave is steadily emitted from the shock front and is linearly polarized. The wave vector 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 B 0 = B 0 ẑ), and the fluctuating electric field is perpendicular to both k and 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 δB 2 z = δB z δE y . 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, ad = 3/2, instead of 4/3 if the downstream plasma was isotropic in all momentum directions. The lack of isotropization is due to the fact that in a σ 1 flow (with downstream plasma magnetically dominated) it will be harder for the plasma to exceed the threshold for velocityspace 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 ∝ σ , 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;Sironi, Keshet & Lemoine 2015;Iwamoto et al. 2017;Pelletier et al. 2017;Plotnikov et al. 2018).

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/ω p ahead of the shock front: 1 5 c/ω p < xx front < 25 c/ω 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

1 This region is delimited by vertical black dashed lines in panels (b) and (g) of Fig. 1.

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

where δE y = E y -E y, 0 = E y + β 0 B 0 . Different solid lines correspond to different values of the magnetization, from σ = 0.1 (blue line) to σ = 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. Fig. 2 shows that:

(i) With increasing σ , 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 β s|d = 0.476 for σ = 0.1 to β s|d = 0.963 for σ = 10, so at higher magnetizations it takes more time for the precursor wave, propagating at β wave 1, to detach from the shock front. Second, there is some interaction z /Bfoot_1 0 , for different values of the upstream magnetization σ . Lines of different colour correspond to a given σ going from 0.1 (blue line) to 10 (black line). The precursor wave energy was extracted from a 20 c/ω p wide slab located at 5 c/ω p < xx front < 25 c/ω p .

occurring between the wave and the upstream plasma, which initially causes a drop in wave efficiency (e.g. at ω p t ∼ 1100 in the black line of Fig. 2). The time required for the wave to selfregulate and settle to a steady state, following this drop, is longer for higher σ (e.g. compare green and black lines in Fig. 2).

(ii) The steady-state value of the normalized wave energy ξ B decreases with increasing magnetization and for σ 1 (cyan, brown, and black lines) it approaches a constant value ξ B 0.01.

(iii) For σ = 0.3 (orange) and σ = 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.

(iv) All the simulations have been evolved for long enough to reach a quasi-stationary state. For the largest explored magnetization σ = 30, the simulation was advanced beyond 2 × 10 4 ω -1 p (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, and spectral width), and the wave strength parameter, as we now describe.

Dependence on the upstream magnetization

In the previous section we have defined the energy fraction in upstream field fluctuations ξ B (see equation 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):

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, β s|d is the velocity of the shock front as measured in the downstream frame of the simulations. The fluctuating components δB u and δ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 (δE x ≈ 0). The latter approximation is fully supported by the simulations and, more fundamentally, by the fact that spacecharge 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 = γ 0 b 0, u . Hence, the magnetization parameter can be rewritten as σ = b 2 0,u /(4πw u ) and we note that δB 2 u /b 2 0,u = γ 2 0 δB 2 u /B 2 z,0 = γ 2 0 ξ 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). 2 We find that the fraction of total incoming energy (including both particle and electromagnetic contributions) that is channelled into the precursor wave can be expressed as

as seen from the DRF. In the following, f ξ 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

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

Concerning the shock front speed and its corresponding bulk Lorentz factor, β s|d and γ 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 β s|d = 0.476c for σ = 0.1 up to β s|d = 0.987c for σ = 30. The Lorentz factor of the shock front tends asymptotically to γ s|d = √ σ , for σ 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 σ 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, σ = 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 σ of the energy fraction parameter f ξ is presented in panel (d) of Fig. 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 per cent for σ = 0.1 down to 0.0065 per cent for σ = 30. The dashed orange line follows the empirical scaling f ξ = 2 × 10 -3 /σ that satisfactorily fits our measured values in the σ > 1 range. The most noticeable result of this panel is that we observe a well-defined scaling f ξ ∝ σ -1 . This result arises from the fact that for σ 1, the normalized wave intensity ξ B is roughly constant and β s|d 1 -1/(2σ ). It follows that in the limit σ 1 the precursor wave carries a constant fraction of the incoming particle kinetic energy, i.e. g ξ 2 × 10 -3 .

Let us emphasize, however, that this σ -dependence of f ξ and g ξ is derived in the DRF (simulation frame). This dependence will be different in the shock front rest frame (SRF), since the front moves with ultrarelativistic speeds for σ 1. This point will be further discussed in Section 4.

Dependence on the upstream bulk Lorentz factor

So far, we have investigated the dependence of the precursor intensity on σ , for a fixed choice of the upstream flow Lorentz factor γ 0 = 10. Here, we demonstrate that ξ B is essentially independent from γ 0 for any value of σ . Let us first consider the dependence at a fixed σ . In Fig. 4 we show the time evolution of the precursor wave energy for σ = 3, when varying γ 0 from 5 to 80. Lines of different colour correspond to different values of γ 0 . Despite large oscillations in time, it appears that ξ B converges to the same value regardless of γ 0 . The time-asymptotic values of ξ 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 γ 0 .

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

In the Appendix A, Fig. A2 shows the values of ξ B obtained for σ ∈ [10 -3 , 1] (horizontal axis) and for γ 0 ranging from 10 to 10 6 (different data sets). This figure shows that in the low magnetization regime σ ∈ [10 -3 , 0.3], the normalized wave intensity ξ B is nearly independent from γ 0 . In the range σ ∈ [0.3, 1] there is a larger scatter among different data sets (which employ different γ 0 ). This range of magnetizations corresponds to the transition cases (see Fig. 2). The most plausible reason for the discrepancy among different data sets 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 ξ B was not yet stabilized (see Fig. 2). In fact, Fig. 4 shows that even at σ > 1 the time-asymptotic value of ξ B is insensitive to the flow Lorentz factor.

Basic considerations

As discussed above, the precursor wave possesses X-mode (extraordinary-mode) polarization, such that its wave vector is perpendicular to B 0 , its fluctuating magnetic field is parallel to B 0 , and its fluctuating electric field is perpendicular to both k and 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)

where double primed quantities are measured in the upstream rest frame (URF). Using Lorentz transformations for ω and k and in the limit γ 2 0 σ , the dispersion relation in the DRF becomes

Interestingly, as long as γ 2 0 σ , 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, dω/dk, with the shock front velocity as:

This relation leads to the cutoff frequency and wavelength

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):

where primed quantities are measured in the SRF. More detailed considerations on the soliton structure of the shock, as presented by Alsop & Arons (1988), give a similar expression for B sh . 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, ω c,sol /ω c = B sh /B 0 .foot_3 Lorentz transforming to the DRF (ω c,sol → ω c,sol ) and using the dispersion relation in equation ( 8) leads to

Based on these arguments, we expect the precursor spectrum to exhibit a low-frequency cutoff at ω cutoff and prominent line-like features at ω c, sol and its harmonics. As we show below, where we compare these scalings with our simulation results, for σ > 1 the predicted ω c, sol systematically overestimates the observed peak frequency ω 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.

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 (kspectrum) and the frequency spectrum (ω-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/ω p < xx front < 105 c/ω 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 ω -1 p , and then calculating its Fourier transform. The spatial window for the k-spectrum and the time interval for the ω-spectrum are chosen so that roughly the same segment of the precursor wave was analysed in the two cases.

In Fig. 5 we present the spectrum of the precursor wave for different σ . Five representative cases are shown from top to bottom, σ = 0.1, 0.3, 1, 3, and 10, respectively. Each panel contains the k-spectrum, plotted using blue solid lines, and the ω-spectrum, plotted using red dashed lines. For the ω-spectrum, the horizontal axis shows ω/ω p , whereas for the k-spectrum we take (k 2 c 2 /ω 2 p + 1) 1/2 . Due to this choice, and given the dispersion relation in equation ( 8), each wavenumber spectrum should nearly overlap with the corresponding frequency spectrum, as Figure 5. Spectrum of the precursor wave for different σ . Five representative cases are shown from top to bottom, σ = 0.1, 0.3, 1, 3, and 10, respectively. The blue solid lines in each panel show the spectrum in k-space and red dashed lines show the same in ω-space. The horizontal axis shows ω/ω p for the ω-spectrum. For the k-spectrum, the choice for the horizontal axis is motivated by the dispersion relation in equation ( 8), so that the k-spectrum should nearly overlap with the corresponding ω-spectrum. The method to compute the spectra is described in the main text. All the spectra are normalized as

The orange vertical lines mark the position of the expected low-frequency cutoff, as given by equation (10).

it is indeed the case. The spectra were normalized such that

In each spectrum, the power drops rapidly below the cutoff frequency given by equation (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 ultrarelativistic energies. The emission is then controlled not by the non-relativistic cyclotron maser, but rather by the ultrarelativistic 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 σ < 1, with the fundamental at ω = ω c, sol or the second harmonic dominating the spectrum at low magnetizations (see the peak at ω 4 ω p for σ = 0.1). In the transition cases with 0.1 < σ < 1, we observe the generation of very strong harmonics up to N = 5, where N = ω/ω c, sol , with high-order harmonics producing stronger lines than the fundamental (see the case with σ = 0.3). For σ > 1 the spectrum shows much less prominent lines. As we will argue later, supplementary amplification mechanisms operate in this regime, and the characteristic frequency ω c, sol given by equation ( 13) no longer controls the location of the spectral peak.

The dependence of the relevant wavelengths and frequencies on the magnetization is presented in Fig. 6. The top row refers to wavelengths, the bottom row to frequencies. The left-hand column shows the variation with σ of the cutoff wavelength λ cutoff (panel a) and cutoff frequency ω 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 equations ( 10) and ( 11), indicated by the black dashed lines. The only exception is the transition case σ 0.3, where the low-frequency cutoff is non-stationary.

The central column (panels c and d) presents the variation with σ 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). The black dashed lines indicate the expectation for soliton emission (equation 13). It is apparent that the peak values obtained in the simulations do not agree with the analytical estimate given by equation ( 13) for any σ > 0.1. 4To understand the disagreement we define two regimes: (i) the transition cases (0.1 < σ < 1) and (ii) the magnetically dominated cases (σ > 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 ω c, sol (see top two panels in Fig. 5). In case (ii), we do not find evidence of any strong line at the expected ω c, sol or its harmonics, but rather we observe less prominent lines at frequencies that have no clear connection with ω c, sol . The measured peak frequency scales as ω peak ≈ 3σ 1/2 ω p . In contrast, from equation ( 13) we would expect a stronger scaling with σ , since ω c, sol → σ ω p in the limit σ 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 σ > 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-hand column (panels e and f, respectively) presents the dependence on σ of the fractional spectral width in wavelength and frequency space ( λ/λ peak and ω/ω peak , respectively). The width ω is the difference between the two frequencies (one above the peak frequency ω peak and one below) where the power drops by a factor of 30 below the peak. The width λ is defined in Figure 6. Characteristic wavelengths and frequencies as a function of σ . Cutoff wavelengths λ cutoff and frequencies ω cutoff are presented in the left-hand column (blue circles in panels a and b, respectively). The dashed black lines show the analytical predictions in equations ( 10) and ( 11). The wavelengths and frequencies at the peak of the spectrum (λ peak and ω peak ) are plotted in the central column (panels c and d) using red squares. The black dashed lines in panels (c) and (d) show the expected ω c, sol from equation ( 13). The green dot-dashed line in panel (c) is the width of the density cavity at the shock front divided by three: L cav /3 (see the text, subsection 3.3.3). The right-hand column presents the dependence on σ of the fractional spectral width λ/λ peak and ω/ω peak (panels e and f, respectively).

an analogous way. This shows quantitatively that the spectrum is narrow, with ω/ω peak 3 nearly independently of σ . The spectra of the cases with σ < 1, that show pronounced line-like features, are even narrower, with line widths of ω/ω peak 1 (see the top two panels in Fig. 5).

Resonating cavity in the shock structure at σ > 1

In the previous subsection we have found that in the magnetically dominated regime σ > 1, the peak frequency in our simulations does not scale as the expected gyration frequency in the soliton, ω c, sol . The physical picture that led to the estimate of ω c, sol must then be revised, since the shock structure for σ > 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 σ 1 regime, we observe for σ 1 the build-up of two density peaks separated by a cavity. 5 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 Fig. 7 we illustrate the structure of the shock transition region for σ = 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 xx shock = 0 and it propagates in the + x direction. The two density peaks near the shock are separated by a cavity of width L cav 1.6 c/ω 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. 5 We 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 σ > 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.

For this particular snapshot, only a mode with wavelength λ = L 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 Fig. 7 we demonstrate the role of the cavity in shaping the precursor spectrum by showing the wavenumber spectrum as a function of λ -1 = k/(2 ). Some characteristic emission wavelengths are easily identified. For instance, the cutoff wavelength at λ cutoff 1.6 c/ω p seems to be closely related to the width of the density cavity L 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 cav /2 and L cav /3, respectively. The latter matches well the position of the strongest emission line. As discussed below, this holds for all σ 1.

To assess the connection between the cavity size and the precursor efficiency we show in panel (c) the time evolution of L cav (blue line) and of the precursor wave energy ξ B multiplied by a factor of 100 (red line). The value of ξ B was computed in a region closer to the shock front than we have done before (here, between 1 and 5 c/ω 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 anticorrelated 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 (ω p t ∼ 1000) as the cavity width. We interpret this behaviour 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 cav has to be roughly equal to λ 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 (ω p t ∼ 1000): we see that large B z fluctuations are first amplified in the cavity, and the emission of 2). The three vertical dashed lines correspond to the cavity width L cav (blue), to L cav /2 (orange), and to L cav /3 (red). Panel (c): time evolution of L cav (blue line) and of the precursor wave energy ξ B multiplied by 100 (red line). The value of ξ B was derived in the region closer to the shock front than previously (between 1 and 5c/ω p ahead of the front). The efficiency settles to a steady state at the same time as the cavity length does. 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 Fig. 6 (panels c and d), where we show that the peak wavelength the precursor emission (red squares) is consistent L cav /3 (green dot-dashed lines in panel c), for σ 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 σ 1 scales as ω peak 3 ω cutoff 3 √ σ ω p 3 ω c in the simulation frame, where we have used that γ s|d √ σ for σ 1. This should be contrasted with equation ( 13), whose scaling (∝ σ ω p in the σ 1 limit) is not supported by our simulations.

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):

where a = e δE y m e cω (16)

is the strength parameter of the wave. δE y is the electric field of the wave and ω 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 on to the wave. We have employed two measures for the wave strength parameter: either from the maximum excursion in u y , a max = max (u y ); or from the root mean square value, a 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/ω p ahead of the front at the final time of the simulations. Fig. 8 presents the dependence on σ of the wave strength parameter derived from our simulations. The maximal value a max decreases from 5 for σ = 0.1 down to 1 for σ = 10, while the root mean square value has the same dependence on σ but it is three times smaller, a std a max /3. The subpanel of this figure shows the dependence on γ 0 . Supplementary simulations were performed for this purpose, where we fixed σ = 3. There is a clear linear dependence of a on γ 0 , as already suggested by Iwamoto et al. (2017).

The linear dependence on γ 0 arises naturally from the fact that ξ B does not depend on γ 0 , combined with the fact that the typical frequency of the precursor wave is 3 ω p for σ 0.1 and 3 ω c for σ > 1. It follows from equation ( 16) that a ≈ γ 0 √ ξ B σ for σ 0.1 and a ≈ γ 0 √ ξ B /3 for σ > 1, which justifies the linear scaling with γ 0 shown in the inset of Fig. 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 Fig. 8 will then be the same in the SRF and in the URF. This is in contrast to the precursor normalized energy ξ B and the precursor spectrum, which are frame-dependent.

From the simulation frame to the shock rest frame

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

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

The ξ B parameter then transforms as (Gallant et al. 1992)

We compute directly ξ B with the following procedure. The values of γ s|d and β s|d obtained from our simulations (see Fig. 3, panels b and c) are used to Lorentz transform the electromagnetic fields into the SRF at a given snapshot of the simulation. Then, ξ B is computed directly, by averaging between 5 and 105 c/ω p ahead of the front (the distance is still measured in the simulation frame).

In Fig. 9 (panel a) we present the dependence of ξ B = ξ B|sh on σ , 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. Secondly, it shows that, beyond the transition cases with 0.1 < σ < 1, where ξ B attains the largest values, the normalized wave energy in the SRF scales as ξ B 7 × 10 -4 σ -2 for σ > 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 σ 1 the wave amplitude in the DRF converges to a constant (i.e. σ -independent) value, ξ B 10 -2 . In addition, the asymptotic shock velocity in the DRF is β s|d 1 -1/(2σ ) for σ 1. Plugging these two scalings into equation ( 19) leads to ξ B = 6.3 × 10 -4 σ -2 , which is very close to the measured scaling.

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 channelled into the precursor). In the SRF, the energy conservation equation including wave contributions can be written as

where 4π) 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

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 subdominant in pair plasmas.

Reminding that the upstream magnetization is Lorentz invariant, we use σ = σ = b 2 0,u /(4πw u ). The fraction of total incoming energy channelled into the precursor wave is then

where β u is the upstream flow velocity measured in SRF. Equivalently, it is the shock front speed in the URF. One can also give the fraction of incoming particle kinetic energy channelled into the precursor wave:

In the latter equation n u is the upstream plasma proper density, including both species (so, n u = 2N 0 /γ 0 ). Getting back to the simulation results, in Fig. 9 (panel b) we present the dependence on σ of the energy fraction f ξ |sh = f ξ , as measured in the SRF. The maximum value is reached at σ ∼ 0.1, where the precursor carries to 5 per cent of the incoming energy. For σ > 0.3 the energy content in the wave rapidly drops. Similarly to the ξ B scaling, there is a clear dependence as ∝ σ -2 for σ > 1 (more precisely, f ξ 7 × 10 -4 σ -2 ). The similarity comes from the fact that the upstream velocity is β u → 1 and the upstream energy content is dominated by the magnetic field (i.e. σ 1). This implies from equation ( 23) that f ξ ξ B .

In the limit σ 1, the conversion efficiency of incoming particle kinetic energy into wave energy scales as g ξ ξ B σ 7 × 10 -4 σ -1 . This should be contrasted with what we have obtained in the DRF, where this quantity became constant in the σ 1 limit. We remark that the scalings reported so far have been obtained from 1D runs. While we expect that the dependence on σ will remain unchanged in 2D and 3D, we speculate that the normalizations of f ξ and g ξ 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. (2017Iwamoto et al. ( , 2018)), performed in the low magnetization regime σ < 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 (σ > 1) of interest for our work, given the rapid decrease of the wave strength parameter with magnetization (see Fig. 8), and so of the wave feedback on to the upstream plasma. This point will be addressed in a forthcoming study (Sironi et al, in preparation).

A P P L I C AT I O N S TO F R B S

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 ultrarelativistic 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σ 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 σ 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 channelled into the precursor waves. If considered in the ejecta frame (post-shock frame), our simulations show that for σ > 1 the emitted wave carries a fraction f ξ = 2 × 10 -3 /σ of the total energy. This corresponds to a fraction g ξ 2 × 10 -3 of the incoming particle kinetic energy regardless of σ . If one considers the energy budget in the SRF, the previous scalings become f ξ 7 × 10 -4 /σ 2 and g ξ 7 × 10 -4 /σ , respectively.

Secondly, 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 (Gajjar et al. 2018;Michilli 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 subdominant in strongly magnetized plasmas.

Thirdly, 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 ω peak 3 ω p for σ 0.1 and as ω peak 3 ω c for σ > 1. Several highorder harmonics characterize the transition region with 0.1 < σ < 1. Joining the two regimes, and neglecting for simplicity the transition cases, we can cast the peak frequency as ω peak 3 ω p max[1, √ σ ]. This can be recast in a simpler form in the pre-shock frame as

as long as the shock is moving with an ultrarelativistic bulk Lorentz factor γ s|u into the upstream medium. The emission frequency for an upstream observer is then

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 γ s|u n e /1 cm -3 ≈ 4 × 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, Margalit & Sironi (2019), this frequency is also consistent with ∼GHz emission from decelerating blast waves produced by flare ejecta in young magnetars.

Finally, the spectrum is narrow-band, ω/ω peak 1 -3 (see Fig. 6), which is again consistent with the observations (e.g. Law et al. 2017;Macquart et al. 2019).

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 γ -1 ) and the particles' pitch-angle distribution is narrowly peaked with spread γ -1 . Here, γ 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).

Secondly, 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 γ ∼ γ 0 , while the inner region of the ring (i.e. at lower γ ) 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 ultrarelativistic, unlike internal shocks. Secondly, the properties of the shell and of the pre-burst wind (as regard to magnetization, temperature, and 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.

S U M M A RY A N D C O N C L U S I O N S

In this work we have investigated by means of 1D PIC simulations the physics of synchrotron maser emission from perpendicular relativistic shocks that propagate in highly magnetized electronpositron plasmas (with magnetization 0.1 ≤ σ ≤ 30). For strongly magnetized shocks, we expect that multidimensional 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 σ and γ 0 . We have found that:

(i) The shock front emits efficiently and steadily a train of highamplitude electromagnetic precursor waves for any σ and γ 0 , in the range 0.1 ≤ σ ≤ 30 and γ 0 ≥ 5 that we have explored. The emission is linearly polarized, with fluctuating magnetic field along the same direction as the upstream mean field.

(ii) 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, ξ B , decreases rapidly from 3.5 at σ = 0.1 down to 0.04 at σ = 1, as measured in the postshock frame of the simulations. For σ 1, this ratio converges to a constant value ξ B 0.01. In the SRF, the asymptotic scaling in the limit σ 1 becomes ξ B ∝ σ -2 . (iii) For σ > 1, the energy output in precursor waves normalized to the total incoming energy scales as f ξ 2 × 10 -3 σ -1 in the post-shock frame and as f ξ 7 × 10 -4 σ -2 in the SRF. The former implies that in the downstream frame, σ > 1 shocks convert a constant fraction of the incoming particle kinetic energy into precursor waves (equal to g ξ 2 × 10 -3 ).

(iv) Magnetically dominated shocks with σ > 1 exhibit a resonating cavity in the shock front structure in between two solitons, instead of the single soliton loop that is observed for σ 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 ω peak 3 ω c = 3 √ σ ω p for σ > 1, whereas earlier works (Gallant et al. 1992) quote a stronger scaling with magnetization, ω peak σ ω p .

(v) The characteristic frequency of the emission, as measured in the post-shock frame, is ω 3ω p for weakly magnetized shocks σ ≤ 0.1, and ω 3ω c for σ 1, as we have just discussed. In the transition region 0.1 < σ < 1, prominent high-order harmonics of ω c, sol (given in equation 13) were observed along with the fundametal at ω 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 ω peak 3 ω p max[1, √ σ ], 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 ω peak ≈ 3γ s|u ω p , where γ s|u is the shock Lorentz factor in the upstream frame.

(vi) The spectrum of the precursor is narrow-band, ω/ω peak 1 -3, with a low-frequency cutoff at ω cutoff = γ s|d ω p (here, γ 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.

(vii) We did not observe any dependence on γ 0 of the energy fraction, ξ B , and of the characteristic emission frequency, ω peak /ω 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. Secondly, we have mostly focused on strongly magnetized (σ > 1) plasmas, a regime that so far has received little attention. Even though this work only presents 1D simulations, we anticipate that the multidimensional physics of σ > 1 shocks (Sironi et al., in preparation) will not depart significantly from what we report here. In contrast, for weaker magnetizations (σ 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 multidimensional effects at low σ 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 ∼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 us 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.

PIC code

t ω p T sim ω p x/(c/ω p ) N ppc k B T 0 /(m e c 2 ) σ min σ max γ 0 TRISTAN-MP 1/200 2 × 10 4 1/100 64 10 -4 10 -1 30 10 SMILEI 1/224 & 1/113 6.7 × 10 3 1/112 20 10 -4 10 -1 30 10 SHOCKAPIC 1/90 1.2 × 10 3 1/44.7 20 10 -4 10 -3 1 1 0 SMILEI (2) 1/90 1 × 10 3 1/44.7 20 10 -6 10 -3 2 1 0 SMILEI (3) 1/90 1.5 × 10 3 1/44.7 20 10 -4 10 -3 2 160

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 A1. 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 γ 0 , and the explored range of σ . 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 the results for any σ negatively. The results become stable for any resolution higher than 40 cells per c/ω p , as long as σ ≤ 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 σ > 1 regime was not explored with this code (as we have discussed, at high σ it takes longer to reach a steady state). With TRISTAN-MP and SMILEI it was possible to reach the stationary state for σ up to 30. Concerning the number of particles per cell, the results are very weakly dependent on N 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.

A1 Precursor energy

In Fig. A1 we present the normalized wave energy ξ B as a function of σ 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 σ , we observe good agreement among different codes. For instance, in the σ > 1 regime SMILEI and TRISTAN-MP give the same values of ξ B . In the range 0.1 < σ ≤ 1, where all codes overlap, the scatter among codes is slightly larger, although the rapid drop in ξ B is common to all codes, and it happens around the same σ . We note that the transition is more abrupt in SHOCKAPIC than in SMILEI and TRISTAN-MP, but differences remain minor.

In Fig. A2 we extend the comparison to different studies in the literature and to different values of γ 0 from 10 to 10 6 . The range of σ 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 σ 1). The blue circles and green diamonds report the values obtained with SMILEI using γ 0 = 10 and 160, respectively. Both give nearly the same values for any explored σ confirming that ξ B does not depend on the flow Lorentz factor. The red squares report the values from TRISTAN- MP using γ 0 = 10 (same as in Fig. A1). The data from the 1D Gallant et al. (1992) using γ 0 = 40 and 10 6 , respectively. We notice that all codes provide the same results in the range of σ ∈ [10 -3 , 0.3], regardless of γ 0 . This demonstrates that the precursor wave normalized energy ξ B is not dependent on γ 0 , and that our study is in very good agreement with earlier results.

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

A2 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 Fig. A3 we compare the precursor spectrum extracted from the three codes for a few representative values of magnetization. From top to bottom, the value of σ is 0.1, 0.5, and 1, respectively. We cannot perform any comparison for σ > 1 as this range was not Figure A4. Comparison of the precursor spectrum in k-space obtained from TRISTAN-MP and SMILEI for σ = 30 (the largest magnetization that we have explored, where differences among codes are most dramatic). The upper panel presents the spectrum from TRISTAN-MP using a fourth-order scheme to solve Maxwell's equations (blue line) and from SMILEI using a Yee-type scheme with CFL number = c t/ x = 0.5 (red line). The latter presents a sharp cutoff at high-k (i.e. for λ < 0.14 c/ω p ) and irregular linelike emission features. The lower panel presents the same comparison but with CFL = 0.99 for SMILEI (red line). The high-k cutoff disappears and a very good agreement with TRISTAN-MP is obtained. We note that the spectra presented in this figure are normalized to unity, instead of the previously adopted normalization |δ Bz

explored with SHOCKAPIC (but see below for a comparison between SMILEI and TRISTAN-MP at σ = 30). The spectrum extracted from TRISTAN-MP is plotted using a solid blue line. The red and orange lines are used for SMILEI and SHOCKAPIC, respectively. There is generally a good agreement among the codes for all values of σ 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 σ = 0.1 and the high-order harmonic line at λ = 0.24 c/ω p for σ = 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 λ -1 c/ω p ∼ 1, for σ = 0.5 and σ = 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. σ > 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 Fig. A3 for σ = 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 t/ 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 t/ 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ω) 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 lightwave 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 σ .

A P P E N D I X B : P R E C U R S O R E N E R G E T I C S : 1 D V E R S U S M U LT I D I M E N S I O NA L S I M U L AT I O N S

In order to support our claim that the precursor wave energy does not significantly decrease due to multidimensional effects (in the σ ≥ 1 regime of interest for this work), here we present a preliminary analysis of 2D and 3D simulations performed with TRISTAN-MP. 6 We explore a range of σ ∈ [0.1, 10] in 2D, and σ ∈ [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/ω p . We find that a transverse width of more than 2 -3 c/ω p is sufficient to capture multidimensional effects. In particular, the effects of wave filamentation 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/ω p (in both y and z directions). The spatial resolution in 3D runs is set to 25 cells per c/ω 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 σ = 3 run, since the spectrum extends to higher frequencies, which are not captured properly if the resolution is insufficient. in 2D and 3D simulations is only a factor of two smaller than in 1D, while for σ < 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ω p = 500 and of the orange line after tω p = 1000). So, we can confidently state that the decrease in precursor efficiency due to multidimensional effects is much less severe in the highmagnetization case σ = 1 than for σ < 0.5.

Let us note that in this appendix we have redefined the ξ B parameter. Here, ξ B corresponds to the normalized Poynting flux in the x-direction, ξ B = δE y δB z -δE z δB y /B 2 0 . The average is done over the region between 5 and 25 c/ω 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 δB 2 z = δE y δB z -δE z δB y = δE y δB z , but this equality is not obviously satisfied in multidimensional simulations with σ ≤ 0.6. The choice of defining ξ 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 Fig. B2, using a suite of 1D, 2D, and 3D simulations, we show the dependence on σ of the normalized Poynting flux of the precursor wave ξ B (panel a), of the energy fraction parameter as measured in the simulation frame f ξ (panel b), and of the energy fraction parameter as measured in the SRF f ξ = f ξ |sh (panel c). The definition of the latter two is given in main body of the article: equation ( 5) and equation ( 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 ξ B in 2D and 3D simulations was done by considering the asymptotic values in the time evolution for each σ , as shown in Fig. B1 for the particular case of σ = 1. Error bars quantify uncertainties due to temporal oscillations of the timeevolution curves. Knowing ξ B and measuring directly the shock front velocities from simulations, the values in panels (b) and (c) were produced using equations ( 5) and ( 23), respectively. Fig. B2 shows that:

(i) In 2D and 3D (red and green symbols), for σ = 0.1 the Poynting flux of the precursor wave ξ B , the energy fractions f ξ and f ξ |sh are reduced by a factor of ≈10-20 as compared to 1D (blue circles). This is in agreement with Iwamoto et al. (2017).

(ii) The suppression in efficiency becomes gradually smaller when σ increases from 0.1 to 3. For σ 1, the difference between 1D and multidimensional results becomes negligible.

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

(iv) If the precursor energy fraction is cast in the SRF, panel (c) shows that f ξ 10 -3 for σ ∼ 0.1 -0.4, instead of ∼0.01 in 1D. For σ > 1, multidimensional simulations converge towards 1D values and follow the scaling f ξ ≈ 5 × 10 -4 /σ 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 (δB y component, since δB ⊥ B 0 for this mode) versus X-mode (δB z component, since δB B 0 for this mode) and beaming of the emitted precursor wave. By extracting systematically the values of δB 2 y and δB 2 z in 3D simulations, we find that the O-mode is subdominant for all magnetizations explored here, i.e. δB 2 y / δB 2 z ∼ 0.2 -0.5. This implies that the precursor wave retains (at the 99 per cent 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 shocknormal direction. The ratio is | y | x ∼ 5 × 10 -2 for any σ ∈ [0.1, 3], where the Poynting vector of the wave is defined as π = δE × δB/B 2 0 . 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 multidimensional runs demonstrates that 1D simulations provide accurate numbers in the σ 1 regime, in agreement with 2D out-of-plane and 3D simulations.