Fast radio bursts as precursor radio emission from monster shocks

The propagation and dissipation of large amplitude waves in strongly magnetized plasma is an issue of considerable interest in high-energy astrophysics. Such waves have long been suspected to be responsible for immense cosmic eruptions, including magnetar flares [1, 2, 3, 4, 5, 6], fast radio bursts (FRBs) [7, 8, 9, 10, 11, 12], delayed gamma-ray emission from a collapsing magnetar [13], X-ray precursors in BNS mergers [5], and conceivably gamma-ray flares from blazars and other sources.

Disturbance of a neutron star magnetosphere, e.g., by star quakes, collapse or collision with a compact companion, generates MHD waves that propagate in the magnetosphere. In general, both Alfvén and magnetosonic modes are expected to be produced during abrupt magnetospheric perturbations, with millisecond periods - a fraction of the stellar radius [5]. The amplitude of such a wave is usually small near the stellar surface, but gradually grows as the wave propagates down the decreasing background dipole field. As the wave enters the nonlinear regime, it is strongly distorted. A periodic fast magnetosonic (FMS) wave, in particular, steepens and eventually forms a shock (termed monster shock) when the wave fields reach a value at which $B^{2}-E^{2}\approx 0$ [3, 5]. Half of the energy carried by the wave is dissipated in the shock and radiated away, producing a bright X-ray burst. The other half can escape the inner magnetosphere without being significantly distorted. This is also true in the case of FMS pulses in which the electric field does not reverse sign. The escaping wave (or pulse) can generate a fast radio burst (FRB), either by compressing the magnetar current sheet [14, 6], or through a maser shock produced far out via collision of the FMS pulse with surrounding matter [7, 9, 15, 10, 16]. As will be shown below, GHz waves are also produced at the precursor of monster shocks and may provide another production mechanism for FRBs.

Steepening and breakdown of FMS waves may also be a viable production mechanism of rapid gamma-ray flares in blazars. In this picture, episodic magnetic reconnection in the inner magnetosphere, at the base of magnetically dominated jet [e.g., 17], can excite large amplitude MHD waves that will propagate along the jet, forming radiative shocks close to the source, that can, conceivably, give rise to large amplitude, short duration flares. Whether detailed calculations support this picture remains to be investigated. But if this mechanism operates effectively in black hole jets, it can alleviate the issue of dissipation of relativistic force-free jets. Striped jets can also produce rapid flares, but the mean power of such jets is likely to be considerably smaller than that of jets produced by ordered fields in MAD states [18, 19]. Dissipation of ordered fields must rely on instabilities, like the current driven kink instability, which are expected to be too slow to account for the rapid variability seen in many blazars, if generated at all.

In this Letter, we study the steepening and breaking of a FMS wave by means of first principles plasma simulations. By comparing simulation results with an analytic model, we are able to derive scaling relations of fundamental properties. As hinted above, we also identify, for the first time, the generation of a high-frequency precursor wave by a synchrotron maser instability at the monster shock. We propose that it might explain the enigmatic fast radio bursts and perhaps other radio transients.

We begin by analyzing the general properties of nonlinear wave propagation and the onset of wave breaking. We consider a FMS wave propagating in a medium having a background magnetic field $\pmb{B}_{\rm bg}=B_{\rm bg}\hat{z}$, and proper density $\rho_{\rm bg}$, where $B_{\rm bg}$ and $\rho_{\rm bg}$ may depend on $x$. The magnetization of the background medium is defined as

with $\sigma_{\rm bg}\gg 1$ for the force-free magnetospheres under consideration here. The wave is injected at time $t=0$ from some source located at $\pmb{x}_{0}$, and propagates in the positive $x$ direction. The electric and magnetic field components are given, respectively, by $\pmb{E}=E(x,t)\hat{y}$ and $\pmb{B}_{\rm w}=\pmb{B}-\pmb{B}_{\rm bg}=B_{\rm w}(x,t)\hat{z}$. In the limit of ideal MHD, $F_{\mu\nu}u^{\nu}$=0, where $F_{\mu\nu}$ is the electromagnetic tensor and $u^{\mu}=(\gamma,\gamma v,0,0)$ the plasma 4-velocity, the plasma 3-velocity equals the drift velocity: $\pmb{v}=\pmb{E}\times\pmb{B}/B^{2}=E/B\,\hat{x}$. The corresponding Lorentz factor is $\gamma=B/\sqrt{B^{2}-E^{2}}$. The simple wave considered here moves along the characteristic $C_{+}$, defined by:

where $a$ is the fast magnetosonic speed measured in the fluid rest frame. In the appendix, it is shown that for a cold plasma $a(\lambda)=\tanh(\lambda/2+c)$ and $v_{+}(\lambda)=\tanh(3\lambda/2+c)$, where $c=\ln\left(\sqrt{\sigma_{\rm bg}}+\sqrt{\sigma_{\rm bg}+1}\right)$ and

For a uniform background field, ${\rm d}B_{\rm bg}/{\rm d}x=0$, $\lambda$ is conserved along the characteristic $C_{+}$ (a Riemann invariant). However, in general, $\lambda$ varies along $C_{+}$ with a dependence on the profile of $B_{\rm bg}$.

For a cold plasma, the wave magnetization relates to $a$ through $\sigma=a^{2}/(1-a^{2})=\sinh^{2}(\lambda/2+c)$. In the regime $\sigma_{\rm bg}\gg 1$ this approximates to:

Since for ideal MHD, the continuity equation implies that $B^{\prime}/\rho$ is conserved, with $B^{\prime}=B/\gamma$ being the magnetic field measured in the fluid rest frame, one finds a compression factor of $B^{\prime}/B_{\rm bg}=\rho/\rho_{\rm bg}=e^{\lambda}$, and $B/B_{\rm bg}=\gamma B^{\prime}/B_{\rm bg}=(1-v)^{-1}$. The relations $E/B=v$ and $B_{w}=B-B_{\rm bg}$ then yield:

Note that $E$ can vary between $E=-B_{\rm bg}/2$ at $v=-1$ ($\lambda\to-\infty$) and $E\to\infty$ at $v\to 1$ ($\lambda\to\infty$). Note also that $(B^{2}-E^{2})/B_{\rm bg}^{2}=e^{2\lambda}$. Thus, $B^{2}-E^{2}\to 0$ as $\lambda\to-\infty$ or $E\to-B_{\rm bg}/2$.

In cases where the background field declines along the characteristic $C_{+}$, viz., ${\rm d}B_{\rm bg}/{\rm d}x<0$, energy conservation implies that $|E|/|B_{\rm bg}|$ increases. For instance, for a wave propagating in the inner magnetosphere of a neutron star, which to a good approximation is a dipole, $|E|/|B_{\rm bg}|\propto r^{2}$. This means that $\lambda$ and, hence, $B^{2}-E^{2}$, decrease along $C_{+}$, ultimately approaching $B^{2}-E^{2}=0$. From the expression for the wave velocity given below Eq. (2), it is seen that $v_{+}=0$ at $\lambda=-2c/3$, or equivalently $(B^{2}-E^{2})/B_{\rm bg}^{2}=(4\sigma_{\rm bg})^{-2/3}$, $\gamma v=\sinh\lambda\approx-(4\sigma_{\rm bg})^{1/3}$. At even smaller $\lambda$ values $v_{+}$ changes sign and the characteristic turns back towards the injection point. Consequently, wave breaking is anticipated around this location. The exact value of $\lambda$ at the moment of shock birth, $\lambda_{\rm s}$, depends on the details. For example, in case of a periodic planar wave with electric field $E(x_{0},t)=E_{0}\sin(\omega t)$ at the boundary $x_{0}$, propagating in a background magnetic field $B_{\rm bg}=B_{0}(x_{0}/x)$, we find $\lambda_{\rm s}=\ln(\omega^{2}x_{0}^{2}/16a\sigma_{\rm bg}^{2})^{1/4}$, where $a=E_{0}/B_{0}$ - see appendix for details. As the wave propagates, the plasma upstream of the shock continues to accelerate backward ($\gamma v=\sinh\lambda<0$), the magnetization declines (Eq. 4), and the shock strengthens. Only wave phases that satisfy $E>-B_{\rm bg}/2$ survive. The other parts are erased through shock dissipation.

We illustrate the analytical model and study the kinetic evolution of the wave after shock formation with 1D PIC simulations performed with the relativistic electromagnetic code Tristan-MP V2 [20]. At $\omega_{\rm p}t\,=\,0$, we launch a periodic FMS wave from $x=0$, letting it propagate along $+\hat{x}$ in a cold pair plasma with a declining static and external background magnetic field, $\pmb{B}_{\rm bg}(x)=B_{0}\,(1+x/\mathcal{R}_{0})^{-1}\,\hat{z}$, constant density and initial background magnetization $\sigma_{0}=1600$, so that $\sigma_{\rm bg}(x)=\sigma_{0}\,(1+x/\mathcal{R}_{0})^{-2}$. To capture MHD scales, the ratio of the wavelength, $\lambda_{w}$, to skin depth, $c/\omega_{\rm p}=\sqrt{m_{\rm e}c^{2}/4\pi n_{\rm bg}e^{2}}$, was taken to be $\lambda_{w}\,\omega_{\rm p}/c\,=\,2\,\omega_{\rm p}/\omega\,=\,1.06\times 10^{4}$ with 60 cells per skin depth, $c\Delta t\,=\,\tfrac{1}{2}\Delta x$, and 50 particles per cell. The amplitude of the wave is set to $0.4\,B_{0}$. The gradient length scale of the background magnetic field, $\mathcal{R}_{0}$, is $\mathcal{R}_{0}\,=\,10\lambda_{w}\,\simeq\,1.06\times 10^{5}\,c/\omega_{\rm p}$. For further details on the derivation and numerics, see the appendix.

Figure 1 illustrates the propagation in the laminar (left), steepening (middle), and shock wave (right) regimes along the gradient. We observe the onset of strong steepening around $x\,=\,0.211\,\mathcal{R}_{0}$. In the steepening region where $B_{z}\gtrsim-E_{y}$ [panel (a.2)], the flow accelerates up to relativistic speeds [panel (c.2)]. The wave then breaks, forming a shock, at $x_{\rm s}\,=\,0.225\,\mathcal{R}_{0}$, corresponding to about 2.25 wavelengths from the injection point. A second shock subsequently forms [Fig. 1(a-b.3)]. The shock forms near the wave trough, as expected from the analytic derivation in the appendix. Panels 1(a-b.3) also reveal a soliton structure, further detailed below, characterized by a sharp peak in density and electromagnetic field at the shock and associated with significant bulk flow heating. A careful examination indicates that the value of the compression factor at wave breaking is $\lambda_{\rm s}=\ln\,(0.1)$, corresponding to a plasma 4-velocity of $u_{x}=-5$. The corresponding value of the invariant at this location is $(B^{2}-E^{2})/B_{\rm bg}^{2}=10^{-2}$. These values are in good agreement with the analytic results derived in the appendix, but note that the density profile adopted there differs from the one employed in the simulations. Following shock formation, the wave develops a plateau at phases where $E=-B_{\rm bg}/2$. This is clearly seen in Fig. 1(a.3). This plateau extends in size over time until complete eradication of the lower part of the wave. The corresponding wave energy is dissipated at the shock. The plasma in the plateau accelerates towards the shock before crossing it, while the magnetization decreases [see Fig. 1(c.3)]. The Lorentz factor just upstream of the shock increases as the wave evolves, reaching a maximum towards the end of the simulations, and then starts declining.

Figure 1 also reveals the generation of high-frequency precursor waves, likely formed by a synchrotron maser [21, 22, 23, 24, 25, 26, 27]. It is seen only in the leading shock, at $x\,=\,0.475\,\mathcal{R}_{0}$ in Fig. 1(c.1-3). The reason why the trailing shock does not generate a precursor wave is the high temperature of the upstream plasma caused by the passage of the leading shock (see panel c.3 in Fig. 1); the synchrotron maser operates efficiently only when the thermal spread of emitting particles is sufficiently small [28]. We emphasize that under realistic conditions, fast cooling of the shocked plasma, which can change the conditions at the trailing shock, is anticipated [5] but ignored in the present analysis. To elucidate the basic features of the shock structure and the precursor wave, we present in Fig. 2 an enlarged view of the immediate shock vicinity. We identify a (double) soliton-like structure [29], where the particle distribution forms a semicoherent cold ring in momentum space (see c.3 in the bottom panel of Fig. 2). A similar structure has been reported earlier for infinite planar shocks [27]. The precursor wave is dominated by a linearly polarized X-mode that propagates superluminaly against the upstream plasma. We also identified a precursor O-mode (Fig. 3), which is considerably weaker than the X-mode and seems unimportant.

The spectrum of the precursor wave (Fig 3), as measured in the Lab frame, exhibits several sharp harmonic peaks around $\omega\simeq 10\,\omega_{p}^{\prime}$, where $\omega^{\prime}_{p}$ denotes the proper plasma frequency, measured in the local rest frame of the fluid in the immediate upstream of the double soliton structure. These peaks correspond to the lowest harmonics of the resonant cavity defined by the double soliton structure seen in Fig. 2. In good agreement with [27], the wavelength of the highest peak and the width of the cavity are proportional, $\lambda_{\rm peak}\sim L_{\rm sh}/3$, with a weak dependence of $L_{\rm sh}$ on the upstream magnetization, $\sigma_{u}$. This implies that the peak frequency, as measured in the Lab frame, should scale as the shock Lorentz factor with respect to the Lab frame, $\gamma_{\rm sh|Lab}\sim\sqrt{\sigma_{u}}$ (see appendix). From our simulations we estimate $\omega_{peak}\approx 1.6\sqrt{\sigma_{u}}\omega_{p}^{\prime}$ for the highest peak. Above the peak, the spectrum extends up to about $70\,\omega_{p}^{\prime}$. Below the peak, it cuts off at a frequency below which the wave is trapped by the shock (that is, the group velocity is smaller than the shock velocity). The latter can be estimated most easily in the rest frame of the upstream plasma. In this frame, the dispersion relation can be written as [30],

where $\omega^{\prime}_{c}=\sigma_{u}\omega^{\prime}_{p}$ is the cyclotron frequency. Upon transforming to the Lab frame, the dispersion relation reduces to $k^{2}=\omega^{2}-\omega_{p}^{{}^{\prime}2}$ in the limit $\gamma^{2}_{\rm sh|u}\gg\sigma_{u}\gg 1$ (see appendix). The cutoff frequency and k-vector are obtained by equating the group velocity, $v_{g}=d\omega/dk$, with the shock velocity with respect to the Lab frame, $v_{\rm sh|Lab}$:

From the analysis of the data presented in Fig 2 we estimate $\sigma_{u}=25$ and $\sigma_{u}/\gamma^{2}_{\rm sh|u}=6\times 10^{-3}$, from which we obtain $u_{\rm sh|Lab}=4.3$ upon transforming to the Lab frame. The cutoff frequency and k-vector given in Eq. (7) are in good agreement with those seen in Fig 3.

Regarding the efficiency of the precursor wave, we find that the X-mode carries a fraction of about $\epsilon_{X}=10^{-3}$ of the total energy dissipated in the shock (the bottom half wave; see Fig. 3). As found in [27], the fraction of incoming kinetic energy converted into precursor waves in the post-shock frame is independent of the magnetization in the limit of $\sigma\gg 1$. The increase in efficiency by a factor of unity ($\lesssim 3$) is consistent with the increase in downstream frame incoming particle kinetic energy across the simulation. We thus anticipate $\epsilon_{X}$ to be weakly sensitive to the conditions prevailing in the medium in which the FMS wave propagates.

We now consider a wave with luminosity $L=cE^{2}r^{2}/2=10^{43}L_{43}$ erg/s, generated near the surface of a magnetar, and assume for simplicity that the wave propagates in the equatorial plane of the dipole background field, where $B_{\rm bg}(r)=10^{15}B_{15}(R/r)^{3}$ G; $R=10^{6}$ cm being the stellar radius. The background density depends on the pair multiplicity $\mathcal{M}$ in the magnetosphere, which is uncertain but expected to be large [31]; we henceforth adopt ${\cal M}=10^{6}{\cal M}_{6}$. With this normalization, the background number density, $n_{\rm bg}=\rho_{\rm bg}/m_{\rm e}$, can be expressed as $n_{\rm bg}=\mathcal{M}n_{\rm GJ}\approx 10^{19}\mathcal{M}_{6}B_{15}\Omega(R/r)^{3}$ cm^-3, where $\Omega$ is the angular velocity of the star, measured in rad/s, and $n_{\rm GJ}=\Omega B_{\rm bg}/2\pi ec$ is the Goldreich–Julian density. The corresponding plasma frequency of the background pair plasma is $\omega_{p}=10^{14}\mathcal{M}_{6}^{1/2}B_{15}^{1/2}(R/r)^{3/2}$ Hz.

The wave propagates nearly undisturbed in the magnetosphere until reaching a radius $r_{\rm s}$ at which the wave breaks. This happens when $|E|=B_{\rm bg}/2$, or $r_{\rm s}\simeq 3\times 10^{8}B_{15}^{1/2}L_{43}^{-1/4}$ cm. As discussed above, a shock forms and continues evolving with the wave. It can be shown [5] that the Lorentz factor just upstream of the shock quickly reaches a maximum value, $\gamma_{u,max}\sim c\sigma_{\rm bg}/\omega r_{\rm s}$ and then gradually declines, roughly as $(r/r_{\rm s})^{-4}$ up to $\sim 3r_{\rm s}$, where the lower half of the wave is completely erased. The majority of the dissipated energy will be released in the $x$-ray and $\gamma$-ray bands, and a fraction of $\sim 10^{-3}$ in the form of a precursor wave. For our choice of parameters, and wave frequency $\nu=\omega/2\pi=10^{4}$ Hz, $\gamma_{u,max}\sim 6\times 10^{5}(B_{15}\Omega\mathcal{M}_{6})^{-1}L_{43}$. The upstream magnetization reaches a minimum, $\sigma_{\rm u}=\sigma_{\rm bg}/2\gamma_{\rm u,max}\sim\omega r_{\rm s}/2c\sim 300B_{15}^{1/2}L_{43}^{-1/4}$, and only slightly increases thereafter [5]. The proper plasma frequency in the shock upstream satisfies $\omega_{p}^{\prime}=\omega_{p}/\gamma_{u}^{1/2}\approx 10^{8}(r/r_{\rm s})^{1/2}\mathcal{M}_{6}B_{15}^{1/4}L_{43}^{-1/8}\Omega^{1/2}$ Hz, and barely changes in the wave dissipation zone, $r_{\rm s}<r<3r_{\rm s}$. Adopting the scaling we find from the simulations, $\omega_{\rm peak}\approx\sqrt{\sigma_{u}}\,\omega_{p}^{\prime}$, we anticipate the observed spectrum of the precursor wave to appear in the GHz band, as seen in many fast radio bursts.

Whether the precursor wave can escape the inner magnetosphere is unclear at present. It has been argued that GHz waves will be strongly damped in the inner magnetosphere by nonlinear decay into Alfvén waves [32], steepening or kinetic effects [33] (but cf. [34]). We note, however, that in the dissipation zone, where the precursor wave is formed, the peak frequency of the precursor wave is smaller than the background frequency by a factor $\sqrt{\gamma_{u}/\sigma_{u}}\sim$ a few, so the wave transitions from the MHD to the kinetic regime. Thus, the damping mechanisms mentioned above need to be reassessed. Moreover, the precursor wave is trapped by the kHz FMS wave, and the effect this has on the damping of the precursor wave needs to be studied. Finally, it has been shown [5] that under magnetar conditions, the cooling rate of the upstream plasma entering the shock is comparable to or even shorter than the Larmor frequency. How this might affect the efficiency of the precursor wave has yet to be determined.

In summary, we have demonstrated the self-consistent steepening, monster shock formation, and precursor wave emission emerging from a fast magnetosonic wave propagating along a declining background magnetic field. The analytical properties of wave steepening are found to be in good agreement with ab initio fully kinetic simulations. The ensuing shock formation leads to efficient dissipation of the bottom part of the wave over the dynamical time of the wave crossing. A fraction of about $10^{-3}$ of the total dissipated energy is imparted to X-modes. The associated spectrum propagating upstream shows pronounced peaks corresponding to harmonics of the cavity forming between two leading solitons at the shock front. Finally, our results open promising avenues to study the fate of electromagnetic pulses propagating in strongly magnetized environments to address self-consistently their dissipation and, if any, the escaping signal.