© The Authors 2024

1. Introduction

Fast radio bursts (FRBs) are bright extragalactic transients of millisecond duration (Cordes & Chatterjee 2019; Petroff et al. 2019, 2022). Observationally, it is now clear that FRBs can be produced by magnetars (Bochenek et al. 2020; CHIME/FRB Collaboration 2020). However, despite significant theoretical effort, the physical processes powering FRBs have not been conclusively identified (Lyubarsky 2021; Zhang 2023).

Any emission model must satisfy the basic constraint that FRBs can escape from the magnetar. Due to the extraordinary properties of FRBs, this is not a trivial requirement. One can define the FRB strength parameter as

$$\begin{array}{c}\hfill a_0=\frac{e{E}_0}{{\omega }_{0}mc},\end{array}$$(1)

where E₀ is the peak electric field of the radio wave, ω₀ is the wave frequency, m and −e respectively are the mass and the charge of the electron, and c is the speed of light. When the background magnetic field is smaller than the field of the radio wave, the strength parameter is equal to the peak transverse component of the electron four-velocity in units of the speed of light. Due to the large radio luminosity of FRBs, the strength parameter is a₀ > 1 at distances R < 10¹³ cm from the magnetar (Luan & Goldreich 2014). It is an open question whether strong radio waves with a₀ > 1 can escape from the magnetar or whether FRBs should be produced at large radii where a₀ < 1.

The magnetosphere and the wind of magnetars are strongly magnetized, and they are composed primarily of electron-positron pairs (Kaspi & Beloborodov 2017). The propagation of strong radio waves with a₀ > 1 in magnetized pair plasmas has recently been studied by several authors (Beloborodov 2021, 2022, 2023; Chen et al. 2022; Qu et al. 2022; Golbraikh & Lyubarsky 2023; Lyutikov 2024). These studies focused on the regime where the wave frequency, ω₀, is smaller than the Larmor frequency in the background magnetic field, ω_L. Two interesting cases have been identified: (i) When ω₀ < ω_L/a₀, the particles oscillate in the field of the wave at nonrelativistic speeds because the background magnetic field is larger than the field of the radio wave. In this case, the plasma current is a linear function of the wave electric field, and radio waves can propagate. (ii) When ω_L/a₀ < ω₀ < ω_L, the particle velocity is relativistic, and the plasma current is nonlinear. It has been suggested that in this case, radio waves are damped, which would constrain emission models of FRBs (Beloborodov 2021, 2022, 2023; Chen et al. 2022). However, the latter issue is still debated (Qu et al. 2022; Lyutikov 2024).

In this paper, we study the propagation of strong waves in the regime ω₀ > ω_L, which has received limited attention so far¹. In this regime, the particle velocity is relativistic. We show that (i) when ω₀ > a₀ω_L, the background magnetic field does not affect the particle motion, and the plasma current is a linear function of the wave electric field. Radio waves can then propagate. Interestingly, the dispersion relation is independent of the wave strength parameter. (ii) When ω_L < ω₀ < a₀ω_L, the plasma current becomes nonlinear, and radio waves could be damped. The propagation regimes are summarized in Table 1.

The paper is organized as follows. In Sect. 2, we present the fundamental equations that govern the propagation of strong waves. In Sect. 3, we consider waves that propagate with superluminal phase velocities, as appropriate for ω₀ > ω_L (in Appendix A, we give an intuitive explanation of the results of Sects. 2 and 3). In Sect. 4, we discuss the implications of our results for the modelling of FRBs and summarize our conclusions.

2. Fundamental equations

We considered a linearly polarized electromagnetic wave that propagates in a cold magnetized pair plasma. In the frame where the plasma ahead of the wave packet is at rest (hereafter, the “lab frame,” or equivalently the proper frame of the magnetar wind), the wave frequency is ω₀ and the wave vector is k₀. The wave propagates along the z direction. We considered the extraordinary mode, where the electric field of the wave is directed along y and the magnetic field of the wave is directed along x. The background magnetic field is also directed along x. In Sect. 4, we argue that our key results are valid for any direction of propagation with respect to the background magnetic field and for any direction of the wave electric field. We worked in units where the speed of light is c = 1. We assumed that all physical quantities depend only on z and t.

Electrons and positrons move in the yz plane. They have the same number densities, the same velocities in the z direction, and opposite velocities in the y direction. The equation of motion of the positrons is

$$\begin{array}{cc}\hfill \frac{\partial u_y}{\partial t}+\frac{u_z}{\gamma }\frac{\partial u_y}{\partial z}& =\frac{e}{m}(E+\frac{u_z}{\gamma }B)\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial u_z}{\partial t}+\frac{u_z}{\gamma }\frac{\partial u_z}{\partial z}& =-\frac{e}{m}\frac{u_y}{\gamma }B,\hfill \end{array}$$(3)

where u = u_ye_y + u_ze_z and $\gamma =\sqrt{1+{\mathit{u}}^2}$ respectively are the four-velocity and the Lorentz factor. The electromagnetic fields, E and B, are the sum of the wave field and the background field. The continuity equation is

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\gamma n\right)+\frac{\partial }{\partial z}\left(n{u}_z\right)=0,\end{array}$$(4)

where n is the proper number density. The evolution of the electromagnetic fields is governed by Faraday’s and Ampère’s laws, which are

$$\begin{array}{cc}\hfill \frac{\partial E}{\partial z}& =\frac{\partial B}{\partial t}\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \frac{\partial B}{\partial z}& =8\pi en{u}_y+\frac{\partial E}{\partial t}.\hfill \end{array}$$(6)

We defined the plasma frequency as ${\omega }_{\mathrm{P}}=\sqrt{8\pi e^2{n}_0/m}$, where n₀ is the proper positron number density ahead of the wave packet, and the Larmor frequency as ω_L = eB_bg/m, where B_bg is the background magnetic field in the lab frame. Since magnetar winds are strongly magnetized, we assumed ω_L > ω_P.

It is convenient to apply the operator ∂/∂z to both sides of Eq. (2). Using Eq. (5) to eliminate ∂E/∂z, one finds

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\frac{\partial u_y}{\partial z}-\frac{e}{m}B)+\frac{\partial }{\partial z}\left[\frac{u_z}{\gamma }(\frac{\partial u_y}{\partial z}-\frac{e}{m}B)\right]=0.\end{array}$$(7)

At this point, it is convenient to introduce the new variable:

$$\begin{array}{c}\hfill Q=\frac{\partial u_y}{\partial z}-\frac{e}{m}B+\frac{\gamma n}{n_0}{\omega }_{\mathrm{L}}.\end{array}$$(8)

Ahead of the wave packet, one has u_y = 0 and n = n₀. In the frame where the particles ahead of the wave packet are at rest, one has γ = 1 and (e/m)B = ω_L. An arbitrary Lorentz boost along the direction of propagation gives (e/m)B = γω_L. Then, one finds Q = 0 in any reference frame. From Eqs. (4) and (7), one can see that the variable Q is governed by

$$\begin{array}{c}\hfill \frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\frac{u_z}{\gamma }Q\right)=0.\end{array}$$(9)

Then, one finds Q = 0 also inside the wave packet. Substituting Q = 0 into Eq. (8), one finds

$$\begin{array}{c}\hfill \frac{e}{m}B=\frac{\partial u_y}{\partial z}+\frac{\gamma n}{n_0}{\omega }_{\mathrm{L}}.\end{array}$$(10)

Substituting Eq. (10) into Eq. (2), one finds

$$\begin{array}{c}\hfill \frac{e}{m}E=\frac{\partial u_y}{\partial t}-\frac{n{u}_z}{n_0}{\omega }_{\mathrm{L}}.\end{array}$$(11)

Substituting Eq. (10) into Eq. (3), one finds

$$\begin{array}{c}\hfill \frac{\partial u_z}{\partial t}+\frac{u_y}{\gamma }\frac{\partial u_y}{\partial z}+\frac{u_z}{\gamma }\frac{\partial u_z}{\partial z}=-\frac{n{u}_y}{n_0}{\omega }_{\mathrm{L}},\end{array}$$(12)

or equivalently

$$\begin{array}{c}\hfill \frac{\partial u_z}{\partial t}+\frac{\partial \gamma }{\partial z}=-\frac{n{u}_y}{n_0}{\omega }_{\mathrm{L}}.\end{array}$$(13)

Substituting Eqs. (10) and (11) into Eq. (6), one finds

$$\begin{array}{c}\hfill \frac{{\partial }^2{u}_y}{\partial t^2}-\frac{{\partial }^2{u}_y}{\partial z^2}+{\omega }_{\mathrm{P}}^2\frac{n}{n_0}{u}_y={\omega }_{\mathrm{L}}[\frac{\partial }{\partial t}\left(\frac{n{u}_z}{n_0}\right)+\frac{\partial }{\partial z}\left(\frac{\gamma n}{n_0}\right)].\end{array}$$(14)

Our final set of equations consists of Eqs. (4), (13), and (14).

3. Superluminal waves

We studied waves that propagate with superluminal phase velocities (i.e., ω₀/k₀ > 1). We worked in the reference frame that moves with the group velocity, k₀/ω₀, along the z direction (hereafter, the “wave frame”). As discussed by Clemmow (1974), in this frame the wave vector vanishes, and the wave frequency is

$$\begin{array}{c}\hfill \omega =\sqrt{{\omega }_0^2-k_0^2}.\end{array}$$(15)

We found an approximate solution for Eqs. (4), (13), and (14) by separating the fast oscillations of the physical quantities, which occur on the temporal scale 1/ω, from their secular evolution, which occurs on a much longer temporal scale.

Particles ahead of the wave packet move along the negative z direction with a large four-velocity, u₀ = k₀/ω ≫ 1. The four-velocity can be presented as

$$\begin{array}{cc}\hfill u_y& =\delta u_{0y}+\frac{1}{2}{u}_{\mathrm{f}y}{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill u_z& =-u_0+\delta u_{0z}+\frac{1}{2}{u}_{\mathrm{f}z}{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.\hfill \end{array}$$(17)

The proper density can be presented as

$$\begin{array}{c}\hfill n=n_0+\delta n_0+\frac{1}{2}{n}_{\mathrm{f}}{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.\end{array}$$(18)

We introduced the real functions δu₀ = δu_0ye_y + δu_0ze_z and δn₀, and the complex functions u_f = u_fye_y + u_fze_z and n_f. For example, δu₀(z, t) describes the secular evolution of the four-velocity inside the wave packet, and u_f(z, t) describes the amplitude of the fast oscillations of the four-velocity. When the envelope of the wave packet is much longer than a few wavelengths, these functions vary on spatial and temporal scales ≫1/ω. The dots in Eqs. (16)–(18) indicate higher-order harmonics, and c.c. indicates the complex conjugate of the oscillating terms. The expansion in harmonics is justified when the plasma current is a linear function of the wave electric field (i.e., for ω₀ ≫ a₀ω_L, as we show below). Second-order harmonics would be important to calculate nonlinear corrections to the dispersion relation on the order of u_f²/u₀², which is out of the scope of the paper.

We considered the regime where the four-velocity inside the wave packet is nearly constant and assumed δu₀/u₀ ∼ u_f²/u₀² ≪ 1. Below we show that the condition δu₀/u₀ ∼ u_f²/u₀² ≪ 1 requires ω₀ ≫ a₀ω_L. In this regime, the plasma current is a linear function of the wave electric field.

Retaining terms of the order of δu₀/u₀ and u_f²/u₀², the Lorentz factor can be presented as

$$\begin{array}{c}\hfill \frac{\gamma }{{\gamma }_0}=1-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{{\left|u_{\mathrm{f}y}\right|}^2}{4{\gamma }_0^2}-\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z}}{2{\gamma }_0}{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\end{array}$$(19)

where ${\gamma }_0=\sqrt{1+u_0^2}={\omega }_0/\omega$. The number density is

$$\begin{array}{cc}\hfill \frac{\gamma n}{{\gamma }_0{n}_0}=1+\frac{\delta n_0}{n_0}& -\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{{\left|u_{\mathrm{f}y}\right|}^2}{4{\gamma }_0^2}-\frac{u_0}{{\gamma }_0}\frac{n_{\mathrm{f}}{u}_{\mathrm{f}z}^{\ast }+n_{\mathrm{f}}^{\ast }{u}_{\mathrm{f}z}}{4{\gamma }_0{n}_0}+\hfill \\ \hfill & +\frac{1}{2}(\frac{n_{\mathrm{f}}}{n_0}-\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z}}{{\gamma }_0}){\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.\hfill \end{array}$$(20)

The y component of the current density is

$$\begin{array}{c}\hfill \frac{n{u}_y}{{\gamma }_0{n}_0}=\frac{\delta u_{0y}}{{\gamma }_0}+\frac{n_{\mathrm{f}}{u}_{\mathrm{f}y}^{\ast }+n_{\mathrm{f}}^{\ast }{u}_{\mathrm{f}y}}{4{\gamma }_0{n}_0}+\frac{u_{\mathrm{f}y}}{2{\gamma }_0}{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\end{array}$$(21)

and the z component is

$$\begin{array}{cc}\hfill \frac{n{u}_z}{{\gamma }_0{n}_0}=-\frac{u_0}{{\gamma }_0}& -\frac{u_0}{{\gamma }_0}\frac{\delta n_0}{n_0}+\frac{\delta u_{0z}}{{\gamma }_0}+\frac{n_{\mathrm{f}}{u}_{\mathrm{f}z}^{\ast }+n_{\mathrm{f}}^{\ast }{u}_{\mathrm{f}z}}{4{\gamma }_0{n}_0}+\hfill \\ \hfill & +\frac{1}{2}(\frac{u_{\mathrm{f}z}}{{\gamma }_0}-\frac{u_0}{{\gamma }_0}\frac{n_{\mathrm{f}}}{n_0}){\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\hfill \end{array}$$(22)

where u_f^* and n_f^* respectively indicate the complex conjugates of u_f and n_f.

3.1. Fast oscillations

When the envelope of the wave packet is much longer than a few wavelengths, the functions δu₀, δn₀, u_f, n_f vary on spatial and temporal scales ≫1/ω. The evolution of the physical quantities on the temporal scale 1/ω can be easily determined. Since in the wave frame the wave vector vanishes, one can neglect the spatial derivatives of the physical quantities in Eqs. (4), (13), (14). One can also approximate the time derivatives as ∂/∂t ≃ −iω. From Eq. (4) one sees that γn should be constant on the temporal scale 1/ω. Then, Eq. (20) gives

$$\begin{array}{c}\hfill \frac{n_{\mathrm{f}}}{n_0}=\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z}}{{\gamma }_0}.\end{array}$$(23)

Substituting Eqs. (17) and (21) into Eq. (13), one finds iωu_fz = ω_Lu_fy, or equivalently

$$\begin{array}{c}\hfill u_{\mathrm{f}z}=-\mathrm{i}\frac{{\omega }_{\mathrm{L}}}{\omega }{u}_{\mathrm{f}y}.\end{array}$$(24)

Substituting Eqs. (16), (21), (22) into Eq. (14), one finds

$$\begin{array}{c}\hfill -{\omega }^2{u}_{\mathrm{f}y}+{\omega }_{\mathrm{P}}^2{u}_{\mathrm{f}y}=-\mathrm{i}\omega {\omega }_{\mathrm{L}}(u_{\mathrm{f}z}-u_0\frac{n_{\mathrm{f}}}{n_0}).\end{array}$$(25)

Substituting Eqs. (23) and (24) into Eq. (25), one finds

$$\begin{array}{c}\hfill {\omega }^2(1-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{\omega }^2})={\omega }_{\mathrm{P}}^2.\end{array}$$(26)

Taking into account that ω² = ω₀² − k₀² and γ₀²ω² = ω₀², one sees that Eq. (26) is equivalent to the standard dispersion relation in the lab frame, which is (ω₀² − k₀²)(1 − ω_L²/ω₀²) = ω_P². The wave is superluminal for ω₀ ≫ ω_L.

Next, we discuss the range of validity of our results. We derived the dispersion relation assuming that the particle four-velocity inside the wave packet is nearly constant (i.e., |u_fy|≪γ₀ and |u_fz|≪γ₀). For superluminal waves with ω₀ ≫ ω_L, one has γ₀ = ω₀/ω_P, |u_fy|=a₀, and |u_fz|=(ω_L/ω_P)|u_fy| (the latter two relations follow from Eqs. (11) and (24)). Since the plasma is strongly magnetized, one has |u_fz|≫|u_fy|. The condition |u_fz|≪γ₀ requires²

$$\begin{array}{c}\hfill {\omega }_0\gg a_0{\omega }_{\mathrm{L}}.\end{array}$$(27)

Equation (27) has a simple physical interpretation. As we show in Sect. 3.2, the plasma inside the wave packet is compressed, and consequently the background magnetic field is amplified. In the lab frame, the amplified magnetic field is on the order of a₀²B_bg. Then, Eq. (27) implies that the amplified magnetic field should be smaller than the wave field. We emphasize that Eq. (27) can be satisfied by strong waves with a₀ > 1 when the wave frequency is much larger than the Larmor frequency.

3.2. Secular evolution

In order to study the secular evolution of the physical quantities, one should average Eqs. (4), (13), (14) over the fast oscillations. In this case, the spatial derivatives of the physical quantities should be retained. From Eq. (4), one finds

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\frac{⟨\gamma n⟩}{{\gamma }_0{n}_0}\right)+\frac{\partial }{\partial z}\left(\frac{⟨n{u}_z⟩}{{\gamma }_0{n}_0}\right)=0.\end{array}$$(28)

Eq. (13) gives

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\frac{⟨u_z⟩}{{\gamma }_0}\right)+\frac{\partial }{\partial z}\left(\frac{⟨\gamma ⟩}{{\gamma }_0}\right)=-\frac{⟨n{u}_y⟩}{{\gamma }_0{n}_0}{\omega }_{\mathrm{L}}.\end{array}$$(29)

Eq. (14) gives

$$\begin{array}{cc}\hfill \frac{{\partial }^2}{\partial t^2}\left(\frac{⟨u_y⟩}{{\gamma }_0}\right)& -\frac{{\partial }^2}{\partial z^2}\left(\frac{⟨u_y⟩}{{\gamma }_0}\right)+{\omega }_{\mathrm{P}}^2\frac{⟨u_y⟩}{{\gamma }_0}=\hfill \\ \hfill & ={\omega }_{\mathrm{L}}[\frac{\partial }{\partial t}\left(\frac{⟨n{u}_z⟩}{{\gamma }_0{n}_0}\right)+\frac{\partial }{\partial z}\left(\frac{⟨\gamma n⟩}{{\gamma }_0{n}_0}\right)].\hfill \end{array}$$(30)

In Appendix B, we present the expressions for the average physical quantities that appear in Eqs. (28)–(30).

We found an approximate solution of Eqs. (28)–(30) by neglecting the time derivatives of the average physical quantities. This “quasi-static” approximation has been extensively used to study laser-plasma interaction (e.g., Sprangle et al. 1990). The approximation is justified because we worked in the frame that moves with the group velocity, where the time derivative of the wave envelope can be neglected. The average physical quantities reach a steady state because (as we show below) they can be expressed as a function of a₀². For example, since ∂a₀²/∂t ≪ ∂a₀²/∂z and ⟨u_y⟩ is a function of a₀², one has ∂⟨u_y⟩/∂t ≪ ∂⟨u_y⟩/∂z.

The average quantities vary on a scale comparable to the length of the wave packet, which in the wave frame is much longer than 1/ω_P. Then, for example, one has ∂²⟨u_y⟩/∂z² ≪ ω_P²⟨u_y⟩. In this case, the third term on the left-hand side of Eq. (30) is much larger than the second term, and Eq. (30) can be approximated as ω_P²⟨u_y⟩/γ₀ = ω_L∂(⟨γn⟩/γ₀n₀)/∂z. Using Eqs. (B.1) and (B.8), the latter equation can be presented as

$$\begin{array}{c}\hfill \frac{\delta u_{0y}}{{\gamma }_0}=\frac{{\omega }_{\mathrm{L}}}{{\omega }_{\mathrm{P}}^2}\frac{\partial }{\partial z}(\frac{\delta n_0}{n_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{a_0^2}{4{\gamma }_0^2}-\frac{u_0^2}{{\gamma }_0^2}\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{a_0^2}{2{\gamma }_0^2}).\end{array}$$(31)

Equation (29) can be approximated as ∂(⟨γ⟩/γ₀)/∂z = −ω_L⟨nu_y⟩/γ₀n₀. When substituting Eqs. (B.4) and (B.6) into the latter equation, and using Eq. (31) to eliminate δu_0y, one finds

$$\begin{array}{cc}\hfill \frac{\partial }{\partial z}& (\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}-\frac{a_0^2}{4{\gamma }_0^2})=\hfill \\ \hfill & =\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{\partial }{\partial z}(\frac{\delta n_0}{n_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{a_0^2}{4{\gamma }_0^2}-\frac{u_0^2}{{\gamma }_0^2}\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{a_0^2}{2{\gamma }_0^2}),\hfill \end{array}$$(32)

which implies

$$\begin{array}{c}\hfill \frac{\delta n_0}{n_0}=\frac{u_0}{{\gamma }_0}(1+\frac{{\omega }_{\mathrm{P}}^2}{{\omega }_{\mathrm{L}}^2})\frac{\delta u_{0z}}{{\gamma }_0}+(\frac{2u_0^2}{{\gamma }_0^2}\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}-1-\frac{{\omega }_{\mathrm{P}}^2}{{\omega }_{\mathrm{L}}^2})\frac{a_0^2}{4{\gamma }_0^2}.\end{array}$$(33)

Eq. (28) can be approximated as ∂⟨nu_z⟩/∂z = 0, which implies

$$\begin{array}{c}\hfill ⟨n{u}_z⟩=-n_0{u}_0.\end{array}$$(34)

Substituting Eq. (34) into Eq. (B.9), one finds

$$\begin{array}{c}\hfill \frac{\delta n_0}{n_0}=\frac{{\gamma }_0}{u_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{a_0^2}{2{\gamma }_0^2}.\end{array}$$(35)

Now one can solve Eqs. (33) and (35), and express δu_0z and δn₀ as a function of a₀². Then, one can substitute δu_0z and δn₀ into Eqs. (B.2), (B.4), and (B.8), and also express ⟨u_z⟩, ⟨γ⟩, and ⟨γn⟩ as a function of a₀². This procedure gives

$$\begin{array}{cc}\hfill ⟨u_z⟩& =-u_0+(1+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2})\frac{a_0^2}{4u_0}\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill ⟨\gamma ⟩& ={\gamma }_0-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{a_0^2}{4{\gamma }_0}\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill ⟨\gamma n⟩& =({\gamma }_0+\frac{a_0^2}{4{\gamma }_0})n_0.\hfill \end{array}$$(38)

The average Lorentz factor and the average four-velocity in the lab frame can be determined from Eqs. (36) and (37). Taking into account that ω₀ ≫ a₀ω_L, one finds

$$\begin{array}{cc}\hfill {⟨u_z⟩}_{\mathrm{lab}}& =\frac{a_0^2}{4}\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill {⟨\gamma ⟩}_{\mathrm{lab}}& =1+\frac{a_0^2}{4}.\hfill \end{array}$$(40)

The background magnetic field does not affect the particle motion in the lab frame, as Eqs. (39) and (40) are identical to the case of a test particle that moves in the field of a vacuum wave (e.g., Gunn & Ostriker 1971).

The average number density and the average current density in the lab frame can be determined from Eqs. (34) and (38). One finds

$$\begin{array}{cc}\hfill {⟨n{u}_z⟩}_{\mathrm{lab}}& =\frac{a_0^2}{4}{n}_0\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill {⟨\gamma n⟩}_{\mathrm{lab}}& =(1+\frac{a_0^2}{4})n_0.\hfill \end{array}$$(42)

Substituting Eqs. (41) and (42) into Eqs. (10) and (11), one can determine the average electromagnetic fields in the lab frame, which are given by

$$\begin{array}{cc}\hfill {⟨E⟩}_{\mathrm{lab}}& =\frac{a_0^2}{4}{B}_{\mathrm{bg}}\hfill \end{array}$$(43)

$$\begin{array}{cc}\hfill {⟨B⟩}_{\mathrm{lab}}& =(1+\frac{a_0^2}{4})B_{\mathrm{bg}}.\hfill \end{array}$$(44)

In the lab frame, the number density and the background magnetic field inside the wave packet are amplified by a large factor on the order of a₀². The amplified magnetic field is smaller than the wave field when ω₀ ≫ a₀ω_L. When the latter condition is violated, particles can be trapped inside the wave packet because ⟨u_z⟩ can vanish in the wave frame (see Eq. (36)).

4. Discussion

We have studied the propagation of superluminal electromagnetic waves with strength parameters a₀ > 1 in magnetized pair plasmas where ω_L > ω_P. Our key result is that strong waves are able to propagate when ω₀ > a₀ω_L, where ω₀ is the wave frequency and ω_L is the Larmor frequency in the background magnetic field. When ω₀ > a₀ω_L, the background magnetic field does not affect the particle motion, and the plasma current is a linear function of the wave electric field. The dispersion relation is independent of the wave strength parameter.

Strong superluminal waves could be damped when ω_L < ω₀ < a₀ω_L for the following reasons. (i) The plasma inside the wave packet is compressed, and its density is amplified by a large factor on the order of a₀². Although we demonstrated this for the extraordinary mode, we expect that the density is amplified also by the ordinary mode, as the particle motion is not affected by the background field and is nearly identical to the case of a test particle in the field of a vacuum wave. (ii) Since the background magnetic field is frozen into the plasma, its strength is also amplified by a factor on the order of a₀². When ω₀ < a₀ω_L, the amplified background field would be larger than the wave field. Then, the wave would be damped, as a significant fraction of its energy would be used to compress the plasma and amplify the background field. The only exception may occur when the wave propagates nearly along the background field, as the background field would not be significantly amplified. This exception does not apply to magnetar winds, where the magnetic field is nearly toroidal and the wave propagates in the radial direction. Fully kinetic simulations can be used to study the damping of strong superluminal waves when ω_L < ω₀ < a₀ω_L.

Our results have important implications for the modelling of FRBs. First we summarize the typical parameters of the problem (see, e.g., Sobacchi et al. 2022). The FRB strength parameter is

$$\begin{array}{c}\hfill a_0=20L_{42}^{1/2}{\nu }_9^{-1}{R}_{12}^{-1},\end{array}$$(45)

where L = 10⁴²L₄₂ erg s⁻¹ is the isotropic equivalent FRB luminosity, ν = 1 ν₉ GHz is the observed FRB frequency, and R = 10¹²R₁₂ cm is the distance from the magnetar. The magnetar light cylinder is located at the radius R_lc = cP/2π, where P = 1 P₀ s is the magnetar period. In the magnetosphere (i.e., for R < R_lc), the magnetar has a dipole field, B_bg = μ/R³, where μ = 10³³μ₃₃ G cm³ is the magnetic moment. In the wind (i.e., for R > R_lc), the magnetic field is B_bg = μ/R_lc²R. In the proper frame of the wind (which we called “lab frame” so far), the ratio of the Larmor frequency and the FRB frequency is

$$\begin{array}{c}\hfill \frac{{\omega }_{\mathrm{L}}}{{\omega }_0}=0.2{\mu }_{33}{\nu }_9^{-1}{P}_0^{-2}{R}_{12}^{-1}.\end{array}$$(46)

Equations (45) and (46) are independent of the wind Lorentz factor because a₀ and ω_L/ω₀ are Lorentz invariant.

Next, we show that FRBs can propagate through the magnetar wind only at relatively large radii³. The condition ω₀ > a₀ω_L requires R > R_crit, where

$$\begin{array}{c}\hfill R_{\mathrm{crit}}=2\times {10}^{12}{L}_{42}^{1/4}{\mu }_{33}^{1/2}{\nu }_9^{-1}{P}_0^{-1}\mathrm{cm}.\end{array}$$(47)

Two consistency checks are needed. Since we considered the propagation of the wave in the wind, R_crit should be larger than the light cylinder radius. The condition R_crit > R_lc is satisfied for

$$\begin{array}{c}\hfill P<20L_{42}^{1/8}{\mu }_{33}^{1/4}{\nu }_9^{-1/2}\mathrm{s}.\end{array}$$(48)

Since we assumed that the strength parameter is large, one should have a₀ > 1 for R = R_crit. This condition is satisfied for

$$\begin{array}{c}\hfill P>0.1L_{42}^{-1/4}{\mu }_{33}^{1/2}\mathrm{s}.\end{array}$$(49)

Eqs. (48) and (49) are satisfied by most magnetars, which have a typical period of a few seconds (Kaspi & Beloborodov 2017).

Our results suggest that FRBs produced at radii R < R_crit cannot propagate through the unperturbed magnetar wind, which challenges several magnetospheric models (e.g., Cordes & Wasserman 2016; Dai et al. 2016; Kumar et al. 2017; Yang & Zhang 2018; Lu & Kumar 2018; Kumar & Bošnjak 2020; Lu et al. 2020). There are two alternative possibilities: (i) FRBs could be produced at radii R > R_crit, as in the shock maser model (Lyubarsky 2014; Beloborodov 2017, 2020; Metzger et al. 2019; Sironi et al. 2021; Iwamoto et al. 2024). (ii) Alternatively, the structure of the wind could be modified during a flare, as in the reconnection model (Lyubarsky 2020; Mahlmann et al. 2022). In this model, FRBs produced at radii R < R_crit propagate on the top of a large amplitude fast-magnetosonic pulse launched from the magnetosphere, and FRBs can escape because their electromagnetic field is much smaller than the field of the pulse.

Acknowledgments

This work was supported by the JSPS KAKENHI Grants 20J00280, 20KK0064, and 22H00130 [M.I.], by the Simons Foundation Grant 00001470 to the Simons Collaboration on Extreme Electrodynamics of Compact Sources (SCEECS) [L.S., T.P.], by the DoE Early Career Award DE-SC0023015 [L.S.], by the Multimessenger Plasma Physics Center (MPPC) NSF Grant PHY-2206609 [L.S.], by the NASA Grant 23-ATP23-0074 [L.S.], by the ISF Grant 2126/22 [T.P.], and by the ERC Advanced Grant Multijets [T.P.]. We acknowledge insightful discussions with Andrei Beloborodov and Yuri Lyubarsky.

References

Appendix A: Test particles in a vacuum wave

The motion of test particles in the electromagnetic fields of a strong vacuum wave has been studied for a long time. Here we summarize the main results of previous studies (e.g., Gunn & Ostriker 1971). We work in the lab frame where the particles ahead of the wave packet are at rest. We set ω_L = 0 in Eqs. (4), (13), (14). We also neglect the feedback of plasma current onto the wave and set ω_P = 0 in Eq. (14). In this case, the y component of the four-velocity can be presented as

$$\begin{array}{c}\hfill u_y=a_{0}cos\left({\omega }_0\xi \right),\end{array}$$(A.1)

where ξ = z − t, and a₀ > 1 is the wave strength parameter.

Eq. (13) can be presented as d(γ − u_z)/dξ = 0. In the lab frame where the particles ahead of the wave packet are at rest, the latter equation has the solution γ − u_z = 1, which gives

$$\begin{array}{c}\hfill \gamma =1+\frac{u_y^2}{2},u_z=\frac{u_y^2}{2}.\end{array}$$(A.2)

Taking into account that ⟨u_y²⟩=a₀²/2, one recovers Eqs. (39)-(40). Eq. (4) can be presented as d[(γ − u_z)n]/dξ = 0, where n is the proper number density. Taking into account that γ − u_z = 1, one finds

$$\begin{array}{c}\hfill n=n_0,\end{array}$$(A.3)

where n₀ is the proper number density ahead of the wave packet. Eqs. (A.1) and (A.3) show that in the test particle approximation the plasma current, nu_y, is a linear function of the wave field.

Now we discuss the validity of the test particle approximation. We focus on pair plasmas, where an electrostatic field is not generated. As we demonstrate formally in Sect. 3, two key conditions should be satisfied. First, the particle velocity along the direction of propagation in the lab frame, v_z = 1 − 2/(2 + u_y²), should be smaller than the group velocity of the wave, v_g ≃ 1 − ω_P²/2ω₀², otherwise particles would be trapped inside the wave packet. This condition is satisfied for

$$\begin{array}{c}\hfill {\omega }_0\gg a_0{\omega }_{\mathrm{P}}.\end{array}$$(A.4)

Second, the average magnetic field inside the wave packet in the lab frame, ⟨B⟩, should be smaller than the peak electric field of the wave, E₀ = (m/e)a₀ω₀. The particles number density inside the wave packet is equal to γn₀, whereas the particles number density ahead of it is equal to n₀. Since the average magnetic field is frozen into the plasma, one has ⟨B⟩=⟨γ⟩B_bg = (1 + a₀²/4)B_bg, where B_bg is the strength of the background magnetic field ahead of the wave packet. The average electric field is ⟨E⟩=⟨v_zB⟩=(a₀²/4)B_bg (where v_z = u_z/γ). The condition ⟨B⟩≪E₀ is satisfied for

$$\begin{array}{c}\hfill {\omega }_0\gg a_0{\omega }_{\mathrm{L}}.\end{array}$$(A.5)

In strongly magnetized plasmas where ω_L ≫ ω_P, the current is a linear function of the wave field when Eq. (A.5) is satisfied.

Appendix B: Average physical quantities

The average four-velocity and the average Lorentz factor can be calculated from Eqs. (16), (17), (19), which respectively give

$$\begin{array}{cc}\hfill \frac{⟨u_y⟩}{{\gamma }_0}& =\frac{\delta u_{0y}}{{\gamma }_0}\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill \frac{⟨u_z⟩}{{\gamma }_0}& =-\frac{u_0}{{\gamma }_0}+\frac{\delta u_{0z}}{{\gamma }_0}\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill \frac{⟨\gamma ⟩}{{\gamma }_0}& =1-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{{\left|u_{\mathrm{f}y}\right|}^2}{4{\gamma }_0^2}.\hfill \end{array}$$(B.3)

Taking into account that |u_fy|=a₀, which follows from Eq. (11), one finds

$$\begin{array}{c}\hfill \frac{⟨\gamma ⟩}{{\gamma }_0}=1-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{a_0^2}{4{\gamma }_0^2}.\end{array}$$(B.4)

The average number density and the average current density can be calculated from Eqs. (20)-(22). Taking into account that n_fu_fy^* + n_f^*u_fy = 0 and n_fu_fz^* + n_f^*u_fz = 2n₀u₀|u_fz|²/γ₀², which follow from Eqs. (23)-(24), one finds

$$\begin{array}{cc}\hfill \frac{⟨\gamma n⟩}{{\gamma }_0{n}_0}& =1+\frac{\delta n_0}{n_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{{\left|u_{\mathrm{f}y}\right|}^2}{4{\gamma }_0^2}-\frac{u_0^2}{{\gamma }_0^2}\frac{{\left|u_{\mathrm{f}z}\right|}^2}{2{\gamma }_0^2}\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill \frac{⟨n{u}_y⟩}{{\gamma }_0{n}_0}& =\frac{\delta u_{0y}}{{\gamma }_0}\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill \frac{⟨n{u}_z⟩}{{\gamma }_0{n}_0}& =-\frac{u_0}{{\gamma }_0}-\frac{u_0}{{\gamma }_0}\frac{\delta n_0}{n_0}+\frac{\delta u_{0z}}{{\gamma }_0}+\frac{u_0}{{\gamma }_0}\frac{{\left|u_{\mathrm{f}z}\right|}^2}{2{\gamma }_0^2}.\hfill \end{array}$$(B.7)

Taking into account that |u_fy|=a₀ and |u_fz|=(ω_L/ω_P)|u_fy|=a₀ω_L/ω_P, which follow from Eqs. (11) and (24), one finds

$$\begin{array}{cc}\hfill \frac{⟨\gamma n⟩}{{\gamma }_0{n}_0}& =1+\frac{\delta n_0}{n_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z}}{{\gamma }_0}+\frac{a_0^2}{4{\gamma }_0^2}-\frac{u_0^2}{{\gamma }_0^2}\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{a_0^2}{2{\gamma }_0^2}\hfill \end{array}$$(B.8)

$$\begin{array}{cc}\hfill \frac{⟨n{u}_z⟩}{{\gamma }_0{n}_0}& =-\frac{u_0}{{\gamma }_0}-\frac{u_0}{{\gamma }_0}\frac{\delta n_0}{n_0}+\frac{\delta u_{0z}}{{\gamma }_0}+\frac{u_0}{{\gamma }_0}\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{a_0^2}{2{\gamma }_0^2}.\hfill \end{array}$$(B.9)

All Tables