© The Authors 2025

1. Introduction

Theoretical models of relativistic quasi-perpendicular electron-positron shocks suggest that diffusive shock acceleration is inhibited because particles sliding along the magnetic field lines cannot cross the shock multiple times. Then, the downstream particles should have a quasi-thermal distribution (Begelman & Kirk 1990; Gallant et al. 1992; Sironi & Spitkovsky 2009). Pulsar wind termination shocks are prototypical examples of such shocks. Multiwavelength observations show that the pulsar wind nebulae are filled with non-thermal particles (Gaensler & Slane 2006; Hester 2008; Bühler & Blandford 2014; Kargaltsev et al. 2015; Reynolds et al. 2017). The hardness of the radio spectrum implies that the acceleration mechanism transfers most of the available energy to a small fraction of the particles, which is not what one would expect in diffusive shock acceleration. However, the morphology of the nebulae indicates that non-thermal particles are accelerated at the termination shock (Komissarov & Lyubarsky 2003, 2004; Del Zanna et al. 2004, 2006; Porth et al. 2014; Olmi et al. 2015). It is clear that theoretical models of shocks lack some key element that allows the acceleration of non-thermal particles.

Most studies attempt to revive particle acceleration by modifying the properties of the upstream plasma. Some possibilities are as follows. (i) The plasma carries alternating magnetic field lines separated by current sheets, which reconnect at the termination shock (Lyubarsky 2003; Pétri & Lyubarsky 2007; Lyubarsky & Liverts 2008; Sironi & Spitkovsky 2011). However, reconnection-driven acceleration of non-thermal particles requires extremely large particle multiplicities. (ii) The plasma is contaminated with protons (Hoshino et al. 1992; Amato & Arons 2006). However, in order to efficiently accelerate the electron-positron pairs, protons should carry a significant fraction of the energy flux, which seems inconsistent with the modeling of pulsar wind nebulae (Bucciantini et al. 2011). (iii) The magnetic field intensity has a large scale gradient in the latitudinal direction (Cerutti & Giacinti 2020). In this scenario, particles are accelerated by shear flow in an elongated cavity that forms in the equatorial region of the shock and by turbulence that develops in the shock downstream.

Another line of research emphasizes the role of the electromagnetic precursor in the shock dynamics. It is now well established that relativistic quasi-perpendicular electron-positron shocks produce strong electromagnetic waves that propagate into the upstream plasma (Langdon et al. 1988; Gallant et al. 1992; Iwamoto et al. 2017, 2018; Plotnikov & Sironi 2019; Babul & Sironi 2020; Margalit et al. 2020; Sironi et al. 2021). The strength parameter of the electromagnetic precursor can be large, namely a₀ > 1 (the strength parameter, a₀, is defined as the peak transverse component of the particles four-velocity in units of the speed of light). Lyubarsky (2018, 2019) argues that absorption of the strong electromagnetic precursor could accelerate particles. These studies consider synchrotron absorption and induced Compton scattering as absorption mechanisms. However, the applicability of these absorption mechanisms to the termination shock of pulsar wind nebulae is not discussed in detail.

The filamentation of the electromagnetic precursor has received limited attention so far, although Lyubarsky (2018, 2019) correctly argues that filamentation can be important for shock dynamics. The main difficulty is that one should consider the regime of large wave strength parameters, a₀ > 1. The filamentation instability has been extensively studied in different fields, including laser-plasma interaction (Kaw et al. 1973; Drake et al. 1974; Max et al. 1974) and fast radio bursts (Sobacchi et al. 2021, 2022, 2023; Ghosh et al. 2022; Iwamoto et al. 2023). These studies focus on the non-relativistic limit a₀ ≪ 1. In this limit, the physics of the instability is understood as follows. Modulations of the radiation intensity in the direction transverse to the wave propagation are expected to be unstable because the ponderomotive force pushes particles out of regions where the intensity is enhanced. The resulting modulation of the plasma refractive index creates a converging lens that further enhances the radiation intensity, thus forming a positive feedback loop. In our previous work (Sobacchi et al. 2024b), we studied the filamentation of electromagnetic waves in unmagnetized electron-positron plasmas in the regime a₀ > 1. However, our results cannot be applied to the precursor of relativistic shocks, where the upstream plasma is magnetized.

In this paper, we study the filamentation instability in magnetized electron-positron plasmas. We focus on the regime a₀ > 1, as appropriate for relativistic shocks. We show that the electromagnetic precursor at the termination shock of pulsar wind nebulae is almost certainly filamented. The instability has important implications for shock dynamics. In the shock frame, the bulk Lorentz factor of the upstream plasma is significantly reduced when the plasma is illuminated by the electromagnetic precursor. Then, the flow becomes sheared when the precursor is filamented. The shear flow distorts the background magnetic field lines and generates a field component parallel to the shock normal that reverses across each radiation filament, potentially triggering magnetic reconnection. Relativistic shear flow and magnetic reconnection may accelerate particles before the plasma enters the shock.

We show that the physics of the filamentation instability can be modified by the presence of a strong background magnetic field. As we explained above, the instability develops because the density of the plasma increases outside radiation filaments and decreases inside the filaments. When the magnetization of the plasma is large, the ponderomotive force is not the main reason why the density is modified. Instead, the dominant effect is the following. Outside radiation filaments the density increases because the tension of the distorted magnetic field lines pushes the plasma in the direction of wave propagation, whereas inside the filaments the density decreases because the direction of the tension force is reversed.

The paper is organized as follows. In Sect. 2, we review the properties of the electromagnetic precursor. In Sect. 3, we derive the equations that govern the evolution of the precursor. In Sect. 4, we calculate the growth rate and the most unstable wave number of the filamentation instability. In Sect. 5, we discuss the implications of the instability for relativistic shocks. In Sect. 6, we summarize our conclusions. Readers not interested in technical details may skip Sects. 3 and 4.

2. Properties of the electromagnetic precursor

The properties of the electromagnetic precursor are extensively discussed in the literature (Langdon et al. 1988; Gallant et al. 1992; Iwamoto et al. 2017, 2018; Plotnikov & Sironi 2019; Babul & Sironi 2020; Margalit et al. 2020; Sironi et al. 2021). In the following, we summarize the available estimates for the frequency and strength parameter of the precursor.

In the pulsar frame, the termination shock is stationary. The pulsar wind has a large Lorentz factor, Γ_w, with respect to the shock. In the proper frame of the pulsar wind (hereafter, the “wind frame”), the particles ahead of the precursor are at rest. We define the plasma frequency as ${\omega }_{\mathrm{P}}=\sqrt{8\pi n_0{e}^2/m}$, where n₀ is the proper positron number density ahead of the precursor, and e and m, respectively, are the charge and the mass of the positron. We define the Larmor frequency as ω_L = eB_bg/mc, where B_bg is the background magnetic field in the wind frame ahead of the precursor, and c is the speed of light. The bulk Lorentz factor of the downstream plasma with respect to the shock is of the order of $\sqrt{1+\sigma }$, where the magnetization is defined as σ = ω_L²/ω_P².

In the wind frame, the frequency of the precursor is (Iwamoto et al. 2017; Plotnikov & Sironi 2019; Sironi et al. 2021)

$$\begin{array}{c}\hfill {\omega }_0\simeq 3{\mathrm{\Gamma }}_{\mathrm{w}}{\omega }_{\mathrm{P}}=3\frac{{\mathrm{\Gamma }}_{\mathrm{w}}}{\sqrt{\sigma }}{\omega }_{\mathrm{L}}.\end{array}$$(1)

When σ > 1, the strength parameter can be presented as a₀ ≃ 0.03 Γ_w, d, where ${\mathrm{\Gamma }}_{\mathrm{w},\mathrm{d}}={\mathrm{\Gamma }}_{\mathrm{w}}/\sqrt{1+\sigma }$ is the Lorentz factor of the wind with respect to the downstream plasma (Plotnikov & Sironi 2019; Sironi et al. 2021). For a given Γ_w, d, the strength parameter increases for decreasing magnetization when 0.1 < σ < 1 and then decreases when 10⁻³ < σ < 0.1 (see Fig. 11 of Iwamoto et al. 2017, and Fig. 8 of Plotnikov & Sironi 2019). In this paper, we assume

$$\begin{array}{c}\hfill a_0\simeq 0.03\frac{{\mathrm{\Gamma }}_{\mathrm{w}}}{\sqrt{1+\sigma }},\end{array}$$(2)

which is the exact expression when σ > 1 and is accurate within a factor of a few when 10⁻³ < σ < 1. Eqs. (1) and (2) are also valid in the regime a₀ > 1 (Margalit et al. 2020). In the wind frame, the magnetic field of the precursor, B₀ = a₀mcω₀/e, can be much larger than the background magnetic field. Using Eqs. (1) and (2), one finds $B_0/B_{\mathrm{bg}}=a_0{\omega }_0/{\omega }_{\mathrm{L}}\simeq 0.1{\mathrm{\Gamma }}_{\mathrm{w}}^2/\sqrt{\sigma (1+\sigma )}$.

In the following, we estimate the frequency and strength parameter of the electromagnetic precursor at the termination shock of the Crab pulsar wind. The Crab pulsar has a period of P = 33 ms and a magnetic moment of μ = 5 × 10³⁰ G cm³. At the light cylinder radius, R_LC = cP/2π, the magnetic field in the pulsar frame is B_LC = μ/R_LC³, and the particle number density is N_LC = κN_GJ, where N_GJ = B_LC/ecP is the Goldreich-Julian density (Goldreich & Julian 1969), and κ is the particle multiplicity. Conservation of energy implies κΓ_w(1 + σ) = eB_LCP/4πmc (Kargaltsev et al. 2015). The latter relation can be presented as Γ_w = 6 × 10⁵(1 + σ)⁻¹κ₅⁻¹, where κ₅ = κ/10⁵. Then, Eq. (2) implies a₀ ≃ 2 × 10⁴(1 + σ)^−3/2κ₅⁻¹.

The particle multiplicity, κ, and the magnetization at the termination shock, σ, are difficult to calculate from first principles. The modeling of the Crab pulsar wind nebula suggests κ > 10⁵ (Bucciantini et al. 2011), which is consistent with the multiplicity that one would expect from pair production in the pulsar polar cap (Timokhin & Harding 2015, 2019). The magnetization of the upstream plasma at the termination shock, σ, depends on the amount of magnetic energy that is dissipated by reconnection in the wind (Lyubarsky & Kirk 2001; Kirk & Skjæraasen 2003). Recent fully kinetic simulations show that alternating magnetic field lines could reconnect relatively close to the light cylinder (Cerutti et al. 2020). After reconnection occurs, the plasma is threaded by a uniform magnetic field. In the equatorial region, where stripes of alternating magnetic field are roughly symmetric before reconnection occurs, one would expect σ < 1. In any case, the wave strength parameter is large for any reasonable value of κ and σ, since condition a₀ > 1 requires κ < 2 × 10⁹(1 + σ)^−3/2.

The Crab pulsar wind terminates at the radius R_TS = 4 × 10¹⁷ cm (Weisskopf et al. 2000). At the termination shock, the proper number density of the positrons is n₀ = N_LCR_LC²/2Γ_wR_TS². Taking into account that N_LC = κB_LC/ecP = 3 × 10¹¹κ₅ cm⁻³, one can calculate the plasma frequency, which is ${\omega }_{\mathrm{P}}=0.02{\kappa }_5\sqrt{1+\sigma }\mathrm{Hz}$. Then, from Eq. (1) one finds ${\omega }_0\simeq 30\mathrm{kHz}/\sqrt{1+\sigma }$.

3. Evolution of the electromagnetic precursor

We consider a linearly polarized quasi-monochromatic wave packet that propagates in a cold magnetized electron-positron plasma. In the wind frame, the wave frequency is ω₀, and the wave vector is $k_0{\mathbf{e}}_z$. The electric field of the wave is directed along ${\mathbf{e}}_y$, and the background magnetic field is directed along ${\mathbf{e}}_x$. The wave polarization corresponds to the extraordinary mode, which is the dominant component of the electromagnetic precursor for large magnetization (Sironi et al. 2021).

We work in the frame that moves with the group velocity, k₀/ω₀, along ${\mathbf{e}}_z$ (hereafter, the “wave frame”). In our units, the speed of light is c = 1. 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}$$(3)

In the wave frame, particles ahead of the wave packet have a large four velocity, ${\mathbf{u}}_0=-u_0{\mathbf{e}}_z$, where u₀ = k₀/ω (the corresponding Lorentz factor is γ₀ = ω₀/ω). The magnetic field is $\mathbf{B}=(m/e){\gamma }_0{\omega }_{\mathrm{L}}{\mathbf{e}}_x$. Since in the wave frame the particles ahead of the wave packet are moving, there is an electric field $\mathbf{E}=(m/e)u_0{\omega }_{\mathrm{L}}{\mathbf{e}}_y$. Inside the wave packet, the magnetic field can be presented as $\mathbf{B}=(m/e)[{\gamma }_0{\omega }_{\mathrm{L}}{\mathbf{e}}_x+\mathrm{\nabla }\times \stackrel{\sim }{\mathbf{a}}]$, and the electric field can be presented as $\mathbf{E}=(m/e)[u_0{\omega }_{\mathrm{L}}{\mathbf{e}}_y-\partial \stackrel{\sim }{\mathbf{a}}/\partial t]$, where $(m/e)\stackrel{\sim }{\mathbf{a}}$ is the vector potential of the wave. We work in the Coulomb gauge, which is $\mathrm{\nabla }·\stackrel{\sim }{\mathbf{a}}=0$. As we demonstrate in Appendix A, the electrostatic potential vanishes because in electron-positron plasmas the wave does not produce any charge separation. As discussed by Bourdier & Fortin (1979), the equation of motion of the particles can be presented as

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}({\mathit{u}}_{\pm }\pm \stackrel{\sim }{\mathit{a}})-\frac{{\mathit{u}}_{\pm }}{{\gamma }_{\pm }}\times [\mathrm{\nabla }\times ({\mathit{u}}_{\pm }\pm \stackrel{\sim }{\mathit{a}})]=-\mathrm{\nabla }{\gamma }_{\pm }\pm (\frac{{\mathit{u}}_{\pm }}{{\gamma }_{\pm }}+\frac{u_0}{{\gamma }_0}{\mathit{e}}_z)\times {\gamma }_0{\omega }_{\mathrm{L}}{\mathit{e}}_x,\end{array}$$(4)

where ${\mathbf{u}}_{\pm }$ is the four-velocity, and γ_± is the Lorentz factor (the plus and minus signs respectively refer to positrons and electrons). The evolution of the proper number density n_± is governed by the continuity equation, which is

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left({\gamma }_{\pm }{n}_{\pm }\right)+\mathrm{\nabla }·\left(n_{\pm }{\mathit{u}}_{\pm }\right)=0.\end{array}$$(5)

The evolution of the wave vector potential is governed by Ampère’s law, which is

$$\begin{array}{c}\hfill \frac{{\partial }^2\stackrel{\sim }{\mathit{a}}}{\partial t^2}-{\mathrm{\nabla }}^2\stackrel{\sim }{\mathit{a}}=\frac{4\pi e^2}{m}(n_{+}{\mathit{u}}_{+}-n_{-}{\mathit{u}}_{-}).\end{array}$$(6)

In Sects. 3.1–3.4, we manipulate Eqs. (4)–(6) to derive the equations that govern the evolution of the wave envelope.

3.1. Weakly non-linear regime

When the Lorentz factor of the particles is nearly constant inside the wave packet (i.e., γ_± ≃ γ₀ in the wave frame), the non-linear terms of Eqs. (4)–(6) can be treated as a small perturbation, because the current $n_{\pm }{\mathbf{u}}_{\pm }$ is roughly proportional to $\stackrel{\sim }{\mathbf{a}}$ (Sobacchi et al. 2024a,b). The key condition γ_± ≃ γ₀ is satisfied when a₀ ≪ ω₀/max[ω_P, ω_L] (Sobacchi et al. 2024a). In the weakly non-linear regime, high-order harmonics of the fields can be neglected. Second-order harmonics would be important to calculate non-linear corrections to the dispersion relation. However, these corrections do not affect the filamentation of electromagnetic waves (Sobacchi et al. 2024b). Then, the vector potential can be presented as

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathit{a}}={\mathit{\psi }}_0+\frac{1}{2}a{\mathit{e}}_y{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\end{array}$$(7)

where c.c. indicates the complex conjugate of the oscillating term proportional to e^−iωt, and the dots indicate high-order harmonics. The complex function $a(\mathbf{x},t)$ describes the envelope of the wave packet. The strength parameter of the wave is a₀ = max[|a|]. We emphasize that we do not assume a₀ ≪ 1. The real function ${\psi }_0(\mathbf{x},t)$ describes the secular evolution of the vector potential. We will show that ${\psi }_0$ does not vanish because the electromagnetic wave compresses the plasma inside the wave packet and consequently amplifies the background magnetic field. The functions a and ${\psi }_0$ vary on spatial and temporal scales ≫1/ω.

As discussed earlier, the particles ahead of the wave packet have a large four-velocity, u₀ = k₀/ω, in the wave frame. Since we are interested in the weakly non-linear regime where the Lorentz factor of the particles inside the wave packet is nearly constant, the four-velocity can be presented as

$$\begin{array}{c}\hfill {\mathit{u}}_{\pm }=-u_0{\mathit{e}}_z+\delta {\mathit{u}}_{0\pm }+\frac{1}{2}{\mathit{u}}_{\mathrm{f}\pm }{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\end{array}$$(8)

where the real function $\delta {\mathbf{u}}_{0\pm }$ describes the secular evolution of the four-velocity, and the complex function ${\mathbf{u}}_{\mathrm{f}\pm }$ describes the envelope of the fast oscillations.

In the following, we study separately the fast oscillations of the physical quantities, which occur on the time scale 1/ω, and their secular evolution. We retain secular terms of the order of $|\delta {\mathbf{u}}_{0\pm }|/{\gamma }_0$ and $|{\mathbf{u}}_{\mathrm{f}\pm }{|}^2/{\gamma }_0^2$, and oscillating terms of the order of $|{\mathbf{u}}_{\mathrm{f}\pm }||\delta {\mathbf{u}}_{0\pm }|/{\gamma }_0^2$. We neglect oscillating terms of the order of $|{\mathbf{u}}_{\mathrm{f}\pm }{|}^3/{\gamma }_0^3$, which do not affect the filamentation of electromagnetic waves (Sobacchi et al. 2024b). The Lorentz factor can be approximated as

$$\begin{array}{c}\hfill \frac{{\gamma }_{\pm }}{{\gamma }_0}=1-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0}+\frac{{\left|{\mathit{u}}_{\mathrm{f}\pm }\right|}^2}{4{\gamma }_0^2}-\frac{u_0^2}{{\gamma }_0^2}\frac{{\left|u_{\mathrm{f}z\pm }\right|}^2}{4{\gamma }_0^2}-[\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z\pm }}{2{\gamma }_0}(1+\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0})-\frac{{\mathit{u}}_{\mathrm{f}\pm }·\delta {\mathit{u}}_{0\pm }}{2{\gamma }_0^2}]{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.\end{array}$$(9)

From Eqs. (8)–(9), one finds

$$\begin{array}{cc}\hfill \frac{{\mathit{u}}_{\pm }}{{\gamma }_{\pm }}=& -\frac{u_0}{{\gamma }_0}(1+\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0}-\frac{{\left|{\mathit{u}}_{\mathrm{f}\pm }\right|}^2}{4{\gamma }_0^2}+\frac{u_0^2}{{\gamma }_0^2}\frac{3{\left|u_{\mathrm{f}z\pm }\right|}^2}{4{\gamma }_0^2}){\mathit{e}}_z+\frac{\delta {\mathit{u}}_{0\pm }}{{\gamma }_0}+\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z\pm }{\mathit{u}}_{\mathrm{f}\pm }^{\ast }+u_{\mathrm{f}z\pm }^{\ast }{\mathit{u}}_{\mathrm{f}\pm }}{4{\gamma }_0^2}+\hfill \\ \hfill & +[\frac{{\mathit{u}}_{\mathrm{f}\pm }}{2{\gamma }_0}(1+\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0})+\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z\pm }}{2{\gamma }_0}\frac{\delta {\mathit{u}}_{0\pm }}{{\gamma }_0}]{\mathrm{e}}^{-\mathrm{i}\omega t}-[\frac{u_0^2}{{\gamma }_0^2}\frac{u_{\mathrm{f}z\pm }}{2{\gamma }_0}(1+\frac{u_0}{{\gamma }_0}\frac{3\delta u_{0z\pm }}{{\gamma }_0})-\frac{u_0}{{\gamma }_0}\frac{{\mathit{u}}_{\mathrm{f}\pm }·\delta {\mathit{u}}_{0\pm }}{2{\gamma }_0^2}]{\mathit{e}}_z{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\hfill \end{array}$$(10)

where ${\mathbf{u}}_{\mathrm{f}\pm }^{\ast }$ denotes the complex conjugate of ${\mathbf{u}}_{\mathrm{f}\pm }$.

3.2. Fast oscillations

In the following, we study the fast oscillations of the physical quantities on the time scale 1/ω. We exploit the fact that the spatial derivatives do not affect the fast oscillations because the wave vector vanishes in the wave frame. The proper number density can be presented as

$$\begin{array}{c}\hfill \frac{n_{\pm }}{n_0}=1+\frac{\delta n_{0\pm }}{n_0}+\frac{n_{\mathrm{f}\pm }}{2n_0}{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.,\end{array}$$(11)

where the real function δn_0± describes the secular evolution of the proper number density, and the complex function n_f± describes the envelope of the fast oscillations. The function n_f± can be determined from Eq. (5). Since spatial derivatives do not affect fast oscillations, the oscillating terms of γ_±n_± should vanish. Using Eq. (9) to eliminate γ_±, one finds

$$\begin{array}{c}\hfill \frac{n_{\pm }}{n_0}=1+\frac{\delta n_{0\pm }}{n_0}+[\frac{u_0}{{\gamma }_0}\frac{u_{\mathrm{f}z\pm }}{2{\gamma }_0}(1+\frac{\delta n_{0\pm }}{n_0}+\frac{u_0}{{\gamma }_0}\frac{2\delta u_{0z\pm }}{{\gamma }_0})-\frac{{\mathit{u}}_{\mathrm{f}\pm }·\delta {\mathit{u}}_{0\pm }}{2{\gamma }_0^2}]{\mathrm{e}}^{-\mathrm{i}\omega t}+\dots +\mathrm{c}.\mathrm{c}.\end{array}$$(12)

From Eqs. (9) and (12), one finds

$$\begin{array}{c}\hfill \frac{{\gamma }_{\pm }{n}_{\pm }}{{\gamma }_0{n}_0}=1+\frac{\delta n_{0\pm }}{n_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0}+\frac{{\left|{\mathit{u}}_{\mathrm{f}\pm }\right|}^2}{4{\gamma }_0^2}-\frac{u_0^2}{{\gamma }_0^2}\frac{3{\left|u_{\mathrm{f}z\pm }\right|}^2}{4{\gamma }_0^2}+\dots +\mathrm{c}.\mathrm{c}.\end{array}$$(13)

From Eqs. (8) and (12), one finds

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

In Eqs. (12) and (14), we retained oscillating terms of the order of $|{\mathbf{u}}_{\mathrm{f}\pm }||\delta {\mathbf{u}}_{0\pm }|/{\gamma }_0^2$ and $|{\mathbf{u}}_{\mathrm{f}\pm }||\delta n_{0\pm }|/{\gamma }_0{n}_0$.

3.2.1. Dispersion relation

In the following, we determine the dispersion relation of the wave at leading order (i.e., we neglect the non-linear terms). These terms are important for the evolution of the wave envelope, which is studied in Sect. 3.2.2. First, one should express the fast oscillations of the four-velocity as a function of a. Substituting Eqs. (7), (8), (10) into Eq. (4), neglecting the spatial derivatives, one finds

$$\begin{array}{c}\hfill (1-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{\omega }^2})u_{\mathrm{f}y\pm }\simeq \mp a,(1-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{\omega }^2})u_{\mathrm{f}z\pm }\simeq \mathrm{i}\frac{{\omega }_{\mathrm{L}}}{\omega }a.\end{array}$$(15)

Then, one should substitute Eqs. (7) and (14) into Eq. (6). On the left-hand side and on the right-hand side, neglecting the non-linear terms, one can approximate respectively ${\partial }^2\stackrel{\sim }{\mathbf{a}}/\partial t^2-{\mathrm{\nabla }}^2\stackrel{\sim }{\mathbf{a}}\simeq -{\omega }^{2}a{\mathbf{e}}_y$ and $n_{+}{\mathbf{u}}_{+}-n_{-}{\mathbf{u}}_{-}\simeq 2n_0{u}_{\mathrm{f}y+}{\mathbf{e}}_y$, where u_fy+ is given by Eq. (15). Then, the dispersion relation is given by

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

Taking into account that ω² = ω₀² − k₀² and γ₀²ω² = ω₀², one sees that Eq. (16) is equivalent to the standard dispersion relation in the wind frame, which is (ω₀² − k₀²)(1 − ω_L²/ω₀²) = ω_P². For ω₀ ≫ ω_L, Eq. (16) implies ω ≃ ω_P, and therefore γ₀ ≃ ω₀/ω_P.

To derive the dispersion relation, we assumed that the Lorentz factor of the particles inside the wave packet is nearly constant (i.e., γ_± ≃ γ₀). Eq. (9) shows that the Lorentz factor is nearly constant for u_fy± ≪ γ₀ and u_fz± ≪ γ₀. Since for ω₀ ≫ ω_L Eq. (15) gives u_f ± y ≃ ∓a and u_f ± z ≃ i(ω_L/ω_P)a, the condition γ_± ≃ γ₀ is satisfied for

$$\begin{array}{c}\hfill a_0\ll \frac{{\omega }_0}{max[{\omega }_{\mathrm{P}},{\omega }_{\mathrm{L}}]},\end{array}$$(17)

where a₀ = max[|a|]. Eq. (17) has a simple physical interpretation (Sobacchi et al. 2024a). The condition a₀ ≪ ω₀/ω_P implies that in the wind frame the longitudinal velocity of the particles inside the wave packet should be smaller than the group velocity. The condition a₀ ≪ ω₀/ω_L implies that the background magnetic field inside the wave packet, which is amplified by a factor of order a₀² in the wind frame, should be smaller than the wave field. We emphasize that Eq. (17) is always satisfied by electromagnetic precursors in relativistic shocks. Taking into account that ω₀ and a₀ are respectively given by Eqs. (1) and (2), one finds $a_{0}max[{\omega }_{\mathrm{P}},{\omega }_{\mathrm{L}}]/{\omega }_0=a_0{\omega }_{\mathrm{P}}max[1,\sqrt{\sigma }]/{\omega }_0\simeq 0.01$.

3.2.2. Evolution of the wave envelope

The equation that governs the evolution of the wave envelope on time scales ≫1/ω can be determined by substituting Eqs. (7) and (14) into the y component of Eq. (6). Considering the resonant terms (i.e., the terms proportional to e^−iωt, or e^iωt), one finds

$$\begin{array}{c}\hfill \frac{{\partial }^{2}a}{\partial t^2}-2\mathrm{i}\omega \frac{\partial a}{\partial t}-{\omega }^{2}a-{\mathrm{\nabla }}^{2}a=\frac{{\omega }_{\mathrm{P}}^2}{2}[\frac{\delta n_{0+}}{n_0}{u}_{\mathrm{f}y+}-\frac{\delta n_{0-}}{n_0}{u}_{\mathrm{f}y-}+\frac{u_0}{{\gamma }_0}(\frac{\delta u_{0y+}}{{\gamma }_0}{u}_{\mathrm{f}z+}-\frac{\delta u_{0y-}}{{\gamma }_0}{u}_{\mathrm{f}z-})]+\frac{{\omega }_{\mathrm{P}}^2}{2}(u_{\mathrm{f}y+}-u_{\mathrm{f}y-}).\end{array}$$(18)

The first term on the left-hand side of Eq. (18) can be neglected because ∂²a/∂t² ≪ ω∂a/∂t. In the first term on the right-hand side, one can approximate u_fy± ≃ ∓a and u_fz± ≃ i(ω_L/ω_P)a, which follow from Eq. (15) since max[ω_P, ω_L]≪ω₀. In the second term, one should retain the non-linear corrections to u_fy+ − u_fy−. Substituting Eqs. (7), (8), (10) into Eq. (4), one finds

$$\begin{array}{c}\hfill (1-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{\omega }^2})(u_{\mathrm{f}y+}-u_{\mathrm{f}y-})=-2a-\mathrm{i}\frac{{\omega }_{\mathrm{L}}}{\omega }(\frac{\delta u_{0y+}}{{\gamma }_0}-\frac{\delta u_{0y-}}{{\gamma }_0})a.\end{array}$$(19)

On the right-hand side of Eq. (19), we approximated u₀ ≃ γ₀, as appropriate because γ₀ ≫ 1. Substituting Eq. (19) into Eq. (18), and using Eq. (16) to eliminate ω, one finds

$$\begin{array}{c}\hfill \frac{\mathrm{i}}{{\omega }_{\mathrm{P}}}\frac{\partial a}{\partial t}=-\frac{1}{2{\omega }_{\mathrm{P}}^2}{\mathrm{\nabla }}^{2}a+\frac{1}{2}\delta \stackrel{\sim }{\rho }a,\end{array}$$(20)

where we defined

$$\begin{array}{c}\hfill \delta \stackrel{\sim }{\rho }=\frac{1}{2}(\frac{\delta n_{0+}}{n_0}+\frac{\delta n_{0-}}{n_0}).\end{array}$$(21)

The evolution of the wave envelope is governed by Eq. (20).

3.3. Secular evolution

Studying the secular evolution of the physical quantities is key for the wave envelope, as Eq. (20) depends on the secular variable $\delta \stackrel{\sim }{\rho }$. One should substitute Eqs. (7)–(14) into Eqs. (4)–(6) and average over the fast oscillations. We start by considering Ampère’s law. Substituting Eqs. (7) and (14) into Eq. (6), and taking into account that u_fy± ≃ ∓a and u_fz± ≃ i(ω_L/ω_P)a, which follow from Eq. (15), one finds

$$\begin{array}{c}\hfill \frac{{\partial }^2}{\partial t^2}\frac{{\mathit{\psi }}_0}{{\gamma }_0}-{\mathrm{\nabla }}^2\frac{{\mathit{\psi }}_0}{{\gamma }_0}=\frac{{\omega }_{\mathrm{P}}^2}{2}[\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}-\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0}-\frac{u_0}{{\gamma }_0}(\frac{\delta n_{0+}}{n_0}-\frac{\delta n_{0-}}{n_0}){\mathit{e}}_z].\end{array}$$(22)

The transverse component of Eq. (22) can be presented as

$$\begin{array}{c}\hfill \frac{\delta {\mathit{u}}_{0\perp +}}{{\gamma }_0}-\frac{\delta {\mathit{u}}_{0\perp -}}{{\gamma }_0}=\frac{2}{{\omega }_{\mathrm{P}}^2}(\frac{{\partial }^2}{\partial t^2}\frac{{\mathit{\psi }}_{0\perp }}{{\gamma }_0}-{\mathrm{\nabla }}^2\frac{{\mathit{\psi }}_{0\perp }}{{\gamma }_0}),\end{array}$$(23)

where we defined $\delta {\mathbf{u}}_{0\perp \pm }=\delta u_{0\pm x}{\mathbf{e}}_x+\delta u_{0\pm y}{\mathbf{e}}_y$ and ${\psi }_{0\perp }={\psi }_{0x}{\mathbf{e}}_x+{\psi }_{0y}{\mathbf{e}}_y$. Now we consider the particles equation of motion. Substituting Eqs. (7)-(10) into Eq. (4), and taking into account that u_fy± ≃ ∓a and u_fz± ≃ i(ω_L/ω_P)a, one finds

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}(\frac{\delta {\mathit{u}}_{0\pm }}{{\gamma }_0}\pm \frac{{\mathit{\psi }}_0}{{\gamma }_0})+\frac{u_0}{{\gamma }_0}{\mathit{e}}_z\times [\mathrm{\nabla }\times (\frac{\delta {\mathit{u}}_{0\pm }}{{\gamma }_0}\pm \frac{{\mathit{\psi }}_0}{{\gamma }_0})]& \pm \mathrm{i}\frac{{\omega }_{\mathrm{L}}}{{\omega }_{\mathrm{P}}}[\frac{a}{2{\gamma }_0}\frac{\partial }{\partial y}\left(\frac{a^{\ast }}{2{\gamma }_0}\right)-\frac{a^{\ast }}{2{\gamma }_0}\frac{\partial }{\partial y}\left(\frac{a}{2{\gamma }_0}\right)]{\mathit{e}}_z=\hfill \\ \hfill & =-\mathrm{\nabla }(\frac{{\left|a\right|}^2}{4{\gamma }_0^2}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0})\pm (\frac{\delta {\mathit{u}}_{0\pm }}{{\gamma }_0}-\frac{u_0^2}{{\gamma }_0^2}\frac{\delta u_{0z\pm }}{{\gamma }_0}{\mathit{e}}_z+\frac{{\left|a\right|}^2}{4{\gamma }_0^2}{\mathit{e}}_z)\times {\omega }_{\mathrm{L}}{\mathit{e}}_x,\hfill \end{array}$$(24)

where a^* denotes the complex conjugate of a. In the derivation of Eq. (24), we neglected terms of order a₀²/γ₀⁴. It is convenient to present two separate equations for the variables $\delta {\mathbf{u}}_{0+}+\delta {\mathbf{u}}_{0-}$ and $\delta {\mathbf{u}}_{0+}-\delta {\mathbf{u}}_{0-}$. As we demonstrate in Appendix B, the equation for the evolution of $\delta {\mathbf{u}}_{0+}+\delta {\mathbf{u}}_{0-}$ is

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}+\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0})-\frac{u_0}{{\gamma }_0}\frac{\partial }{\partial z}(\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}+\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0})=-\mathrm{\nabla }\frac{{\left|a\right|}^2}{2{\gamma }_0^2}+\frac{1}{{\gamma }_0^2}(\frac{\delta u_{0z+}}{{\gamma }_0}-\frac{\delta u_{0z+}}{{\gamma }_0}){\omega }_{\mathrm{L}}{\mathit{e}}_y+\frac{2{\omega }_{\mathrm{L}}}{{\omega }_{\mathrm{P}}^2}({\mathrm{\nabla }}^2\frac{{\psi }_{0y}}{{\gamma }_0}-\frac{{\partial }^2}{\partial t^2}\frac{{\psi }_{0y}}{{\gamma }_0}){\mathit{e}}_z.\end{array}$$(25)

The first term on the right-hand side of Eq. (25) is the ponderomotive force. The equation for the evolution of $\delta {\mathbf{u}}_{0+}-\delta {\mathbf{u}}_{0-}$ is

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}(\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}-\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0}+\frac{2{\mathit{\psi }}_0}{{\gamma }_0})& -\frac{u_0}{{\gamma }_0}\frac{\partial }{\partial z}(\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}-\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0}+\frac{2{\mathit{\psi }}_0}{{\gamma }_0})+\frac{u_0}{{\gamma }_0}\mathrm{\nabla }\frac{2{\psi }_{0z}}{{\gamma }_0}+\mathrm{i}\frac{2{\omega }_{\mathrm{L}}}{{\omega }_{\mathrm{P}}}[\frac{a}{2{\gamma }_0}\frac{\partial }{\partial y}\left(\frac{a^{\ast }}{2{\gamma }_0}\right)-\frac{a^{\ast }}{2{\gamma }_0}\frac{\partial }{\partial y}\left(\frac{a}{2{\gamma }_0}\right)]{\mathit{e}}_z=\hfill \\ \hfill & =-(\frac{\delta u_{0y+}}{{\gamma }_0}+\frac{\delta u_{0y-}}{{\gamma }_0}){\omega }_{\mathrm{L}}{\mathit{e}}_z+\frac{1}{{\gamma }_0^2}(\frac{\delta u_{0z+}}{{\gamma }_0}+\frac{\delta u_{0z-}}{{\gamma }_0}){\omega }_{\mathrm{L}}{\mathit{e}}_y+\frac{{\left|a\right|}^2}{2{\gamma }_0^2}{\omega }_{\mathrm{L}}{\mathit{e}}_y.\hfill \end{array}$$(26)

Finally, we consider the continuity equation. Substituting Eqs. (13)-(14) into Eq. (5), and taking into account that u_fy± ≃ ∓a and u_fz± ≃ i(ω_L/ω_P)a, one finds

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\frac{\delta n_{0\pm }}{n_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z\pm }}{{\gamma }_0}+\frac{{\left|a\right|}^2}{4{\gamma }_0^2}-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{2{\gamma }_0^2})+\mathrm{\nabla }·(\frac{\delta {\mathit{u}}_{0\pm }}{{\gamma }_0}-\frac{u_0}{{\gamma }_0}\frac{\delta n_{0\pm }}{n_0}{\mathit{e}}_z+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{2{\gamma }_0^2}{\mathit{e}}_z)=0.\end{array}$$(27)

In the derivation of Eq. (27), we neglected terms of order a₀²/γ₀⁴. It is convenient to present an equation for the variable $\delta \stackrel{\sim }{\rho }$, which appears in Eq. (20). One finds

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(2\delta \stackrel{\sim }{\rho }-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z+}}{{\gamma }_0}-\frac{u_0}{{\gamma }_0}\frac{\delta u_{0z-}}{{\gamma }_0}+\frac{{\left|a\right|}^2}{2{\gamma }_0^2}-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{{\gamma }_0^2})-2\frac{u_0}{{\gamma }_0}\frac{\partial \delta \stackrel{\sim }{\rho }}{\partial z}+\frac{\partial }{\partial z}\left(\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{{\gamma }_0^2}\right)+\mathrm{\nabla }·(\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}+\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0})=0.\end{array}$$(28)

In the following, we make some algebraic manipulations to simplify our set of equations. In Sect. 3.3.1, we study the evolution of the background electromagnetic fields. Then, in Sect. 3.3.2, we study the velocity and density of the fluid.

3.3.1. Background electromagnetic fields

The background electromagnetic fields are described by the function ${\psi }_0$. To study the secular evolution of the magnetic field, it is convenient to introduce the variable $\delta {\omega }_{\mathrm{L}}=\mathrm{\nabla }\times {\psi }_0$, or equivalently

$$\begin{array}{c}\hfill \delta {\omega }_{\mathrm{L}x}=\frac{\partial {\psi }_{0z}}{\partial y}-\frac{\partial {\psi }_{0y}}{\partial z},\delta {\omega }_{\mathrm{L}y}=\frac{\partial {\psi }_{0x}}{\partial z}-\frac{\partial {\psi }_{0z}}{\partial x},\delta {\omega }_{\mathrm{L}z}=\frac{\partial {\psi }_{0y}}{\partial x}-\frac{\partial {\psi }_{0x}}{\partial y}.\end{array}$$(29)

Since $\mathrm{\nabla }·\delta {\omega }_{\mathrm{L}}=0$, one has

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{L}z}}{\partial z}=-\frac{\partial \delta {\omega }_{\mathrm{L}x}}{\partial x}-\frac{\partial \delta {\omega }_{\mathrm{L}y}}{\partial y}.\end{array}$$(30)

To study the secular evolution of the electric field, it is convenient to introduce the variable $\delta {\omega }_{\mathrm{E}}=-\partial {\psi }_0/\partial t$, or equivalently

$$\begin{array}{c}\hfill \delta {\omega }_{\mathrm{E}x}=-\frac{\partial {\psi }_{0x}}{\partial t},\delta {\omega }_{\mathrm{E}y}=-\frac{\partial {\psi }_{0y}}{\partial t},\delta {\omega }_{\mathrm{E}z}=-\frac{\partial {\psi }_{0z}}{\partial t}.\end{array}$$(31)

Since $\mathrm{\nabla }·{\psi }_0=0$ in the Coulomb gauge, one has

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial z}=-\frac{\partial \delta {\omega }_{\mathrm{E}x}}{\partial x}-\frac{\partial \delta {\omega }_{\mathrm{E}y}}{\partial y}.\end{array}$$(32)

From Eqs. (29) and (31), one finds

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{L}x}}{\partial t}=\frac{\partial \delta {\omega }_{\mathrm{E}y}}{\partial z}-\frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial y},\frac{\partial \delta {\omega }_{\mathrm{L}y}}{\partial t}=\frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial x}-\frac{\partial \delta {\omega }_{\mathrm{E}x}}{\partial z}.\end{array}$$(33)

It is also convenient to introduce the variable

$$\begin{array}{c}\hfill \delta \mathbf{v}=\frac{1}{2}(\frac{\delta {\mathit{u}}_{0+}}{{\gamma }_0}+\frac{\delta {\mathit{u}}_{0-}}{{\gamma }_0}).\end{array}$$(34)

First, we consider the x component of Eq. (26). Using Eq. (23) to eliminate δu_0 + x − δu_0 − x and taking into account that ψ_0x varies on spatial and temporal scales ≫1/ω_P, one sees that terms proportional to δu_0 + x − δu_0 − x can be neglected. Then, the x component of Eq. (26) can be approximated as

$$\begin{array}{c}\hfill \delta {\omega }_{\mathrm{E}x}=-\frac{u_0}{{\gamma }_0}\delta {\omega }_{\mathrm{L}y}.\end{array}$$(35)

Substituting Eq. (35) into Eq. (33), one finds

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{L}y}}{\partial t}-\frac{u_0}{{\gamma }_0}\frac{\partial \delta {\omega }_{\mathrm{L}y}}{\partial z}=\frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial x}.\end{array}$$(36)

Next, we consider the y component of Eq. (26). As we discussed in the derivation of Eq. (35), terms proportional to δu_0 + y − δu_0 − y can be neglected. Then, one can approximate

$$\begin{array}{c}\hfill \delta {\omega }_{\mathrm{E}y}=\frac{u_0}{{\gamma }_0}\delta {\omega }_{\mathrm{L}x}-\frac{{\left|a\right|}^2}{4{\gamma }_0}{\omega }_{\mathrm{L}}-\frac{\delta v_z}{{\gamma }_0}{\omega }_{\mathrm{L}}.\end{array}$$(37)

Substituting Eq. (37) into Eq. (33), one finds

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{L}x}}{\partial t}-\frac{\partial }{\partial z}(\frac{u_0}{{\gamma }_0}\delta {\omega }_{\mathrm{L}x}-\frac{{\left|a\right|}^2}{4{\gamma }_0}{\omega }_{\mathrm{L}}-\frac{\delta v_z}{{\gamma }_0}{\omega }_{\mathrm{L}})=-\frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial y}.\end{array}$$(38)

Substituting Eqs. (35) and (37) into Eq. (32), one finds

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial z}=\frac{u_0}{{\gamma }_0}\frac{\partial \delta {\omega }_{\mathrm{L}y}}{\partial x}-\frac{\partial }{\partial y}(\frac{u_0}{{\gamma }_0}\delta {\omega }_{\mathrm{L}x}-\frac{{\left|a\right|}^2}{4{\gamma }_0}{\omega }_{\mathrm{L}}-\frac{\delta v_z}{{\gamma }_0}{\omega }_{\mathrm{L}}).\end{array}$$(39)

Eqs. (30) and (39) can be presented in a more convenient form. Applying the operator ∂/∂t − (u₀/γ₀)∂/∂z to both sides of Eq. (30), and using Eqs. (36) and (38) to eliminate δω_Lx and δω_Ly, one finds

$$\begin{array}{c}\hfill \frac{\partial }{\partial z}(\frac{\partial \delta {\omega }_{\mathrm{L}z}}{\partial t}-\frac{u_0}{{\gamma }_0}\frac{\partial \delta {\omega }_{\mathrm{L}z}}{\partial z})=\frac{{\partial }^2}{\partial x\partial z}(\frac{{\left|a\right|}^2}{4{\gamma }_0}{\omega }_{\mathrm{L}}+\frac{\delta v_z}{{\gamma }_0}{\omega }_{\mathrm{L}}),\end{array}$$(40)

which implies

$$\begin{array}{c}\hfill \frac{\partial \delta {\omega }_{\mathrm{L}z}}{\partial t}-\frac{u_0}{{\gamma }_0}\frac{\partial \delta {\omega }_{\mathrm{L}z}}{\partial z}=\frac{\partial }{\partial x}(\frac{{\left|a\right|}^2}{4{\gamma }_0}{\omega }_{\mathrm{L}}+\frac{\delta v_z}{{\gamma }_0}{\omega }_{\mathrm{L}}).\end{array}$$(41)

Similarly, applying the operator ∂/∂t − (u₀/γ₀)∂/∂z to both sides of Eq. (39), one finds

$$\begin{array}{c}\hfill \frac{u_0}{{\gamma }_0}(\frac{{\partial }^2\delta {\omega }_{\mathrm{E}z}}{\partial x^2}+\frac{{\partial }^2\delta {\omega }_{\mathrm{E}z}}{\partial y^2}+\frac{{\partial }^2\delta {\omega }_{\mathrm{E}z}}{\partial z^2})=\frac{\partial }{\partial t}[\frac{\partial \delta {\omega }_{\mathrm{E}z}}{\partial z}-\frac{\partial }{\partial y}(\frac{{\left|a\right|}^2}{4{\gamma }_0}{\omega }_{\mathrm{L}}+\frac{\delta v_z}{{\gamma }_0}{\omega }_{\mathrm{L}})].\end{array}$$(42)

The evolution of δω_Lx, δω_Ly, δω_Lz, δω_Ez is governed by Eqs. (36), (38), (41), (42). After determining the solution of these equations, one can calculate δω_Ex from Eq. (35) and δω_Ey from Eq. (37).

3.3.2. Velocity and density of the fluid

The velocity and density of the fluid are described by the functions $\delta \stackrel{\sim }{\rho }$ and $\delta \mathbf{v}$. The transverse component of Eq. (25) gives

$$\begin{array}{c}\hfill \frac{\partial \delta {\mathit{v}}_{\perp }}{\partial t}-\frac{u_0}{{\gamma }_0}\frac{\partial \delta {\mathit{v}}_{\perp }}{\partial z}=-{\mathrm{\nabla }}_{\perp }\frac{{\left|a\right|}^2}{4{\gamma }_0^2}+\frac{{\omega }_{\mathrm{L}}}{{\gamma }_0^2}\delta j_z{\mathit{e}}_y,\end{array}$$(43)

where we defined

$$\begin{array}{c}\hfill \delta j_z=\frac{1}{2}(\frac{\delta u_{0z+}}{{\gamma }_0}-\frac{\delta u_{0z-}}{{\gamma }_0}).\end{array}$$(44)

Using Eq. (C.2) to simplify the right-hand side, the z component of Eq. (25) gives

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\delta v_z+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{4{\gamma }_0^2})-\frac{u_0}{{\gamma }_0}\frac{\partial }{\partial z}(\delta v_z-\frac{{\left|a\right|}^2}{4{\gamma }_0^2}-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{4{\gamma }_0^2})=\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{\partial }{\partial x}\left(\frac{\delta {\omega }_{\mathrm{L}z}}{{\gamma }_0{\omega }_{\mathrm{L}}}\right)-\frac{u_0}{{\gamma }_0}\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{\partial }{\partial y}\left(\frac{\delta {\omega }_{\mathrm{E}z}}{{\gamma }_0{\omega }_{\mathrm{L}}}\right)-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{\omega }_{\mathrm{P}}^2}\frac{\partial }{\partial z}\left(\frac{\delta {\omega }_{\mathrm{L}x}}{{\gamma }_0{\omega }_{\mathrm{L}}}\right).\end{array}$$(45)

The evolution of δj_z is governed by the z component of Eq. (26), which gives

$$\begin{array}{c}\hfill \frac{\partial \delta j_z}{\partial t}-\frac{u_0}{{\gamma }_0}\frac{\partial \delta j_z}{\partial z}=\frac{\delta {\omega }_{\mathrm{E}z}}{{\gamma }_0}-{\omega }_{\mathrm{L}}\delta v_y-\mathrm{i}\frac{{\omega }_{\mathrm{L}}}{{\omega }_{\mathrm{P}}}[\frac{a}{2{\gamma }_0}\frac{\partial }{\partial y}\left(\frac{a^{\ast }}{2{\gamma }_0}\right)-\frac{a^{\ast }}{2{\gamma }_0}\frac{\partial }{\partial y}\left(\frac{a}{2{\gamma }_0}\right)].\end{array}$$(46)

The evolution of $\delta \stackrel{\sim }{\rho }$ is governed by Eq. (28), which gives

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\delta \stackrel{\sim }{\rho }-\frac{u_0}{{\gamma }_0}\delta v_z+\frac{{\left|a\right|}^2}{4{\gamma }_0^2}-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{2{\gamma }_0^2})-\frac{\partial }{\partial z}(\frac{u_0}{{\gamma }_0}\delta \stackrel{\sim }{\rho }-\delta v_z-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{{\left|a\right|}^2}{2{\gamma }_0^2})=-{\mathrm{\nabla }}_{\perp }·\delta {\mathit{v}}_{\perp }.\end{array}$$(47)

The evolution of the fluid is governed by Eqs. (43), (45), (46), (47).

3.4. Final set of equations

To simplify the equations that govern the secular evolution of the system, we neglect the time derivative of |a|. This “quasi-static” approximation is extensively used to study the laser-plasma interaction (e.g., Sprangle et al. 1990). The approximation is justified because we work in the frame that moves with the group velocity, where the time derivative of the wave envelope can be neglected. As we demonstrate in Appendix D, the equations that govern the evolution of the system can be approximated as

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\left(\frac{\delta {\omega }_{\mathrm{L}z}}{{\gamma }_0{\omega }_{\mathrm{L}}}\right)-\frac{\partial }{\partial z}\left(\frac{\delta {\omega }_{\mathrm{L}z}}{{\gamma }_0{\omega }_{\mathrm{L}}}\right)& =\frac{\partial }{\partial x}\frac{{\left|a\right|}^2}{4{\gamma }_0^2}\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill \frac{\partial \delta {\mathit{v}}_{\perp }}{\partial t}-\frac{\partial \delta {\mathit{v}}_{\perp }}{\partial z}& =-{\mathrm{\nabla }}_{\perp }\frac{{\left|a\right|}^2}{4{\gamma }_0^2}+\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2}\left(\frac{\delta j_z}{{\omega }_{\mathrm{L}}}\right){\mathit{e}}_y\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}\left(\frac{\delta j_z}{{\omega }_{\mathrm{L}}}\right)-\frac{\partial }{\partial z}\left(\frac{\delta j_z}{{\omega }_{\mathrm{L}}}\right)& =-\delta v_y\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill \frac{\partial \delta \rho }{\partial t}-\frac{\partial \delta \rho }{\partial z}& =-{\mathrm{\nabla }}_{\perp }·\delta {\mathit{v}}_{\perp }+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{\partial }{\partial x}\left(\frac{\delta {\omega }_{\mathrm{L}z}}{{\gamma }_0{\omega }_{\mathrm{L}}}\right),\hfill \end{array}$$(51)

where we defined

$$\begin{array}{c}\hfill \delta \rho =\delta \stackrel{\sim }{\rho }-\frac{{\left|a\right|}^2}{4{\gamma }_0^2}-\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2}\frac{3{\left|a\right|}^2}{4{\gamma }_0^2}.\end{array}$$(52)

Since we work in the frame that moves with the group velocity of the wave packet, it is interesting to consider the case where the variables δω_Lz, $\delta {\mathbf{v}}_{\perp }$, δj_z and δρ depend only on z (and are independent of x, y, t). In this approximation, one would recover the results of Sobacchi et al. (2024a).

The evolution of the wave envelope is governed by Eq. (20). Since non-linear terms of order a₀²/γ₀² do not affect the filamentation of electromagnetic waves (Sobacchi et al. 2024b), Eq. (20) is equivalent to

$$\begin{array}{c}\hfill \frac{\mathrm{i}}{{\omega }_{\mathrm{P}}}\frac{\partial a}{\partial t}=-\frac{1}{2{\omega }_{\mathrm{P}}^2}{\mathrm{\nabla }}^{2}a+\frac{1}{2}\delta \rho a.\end{array}$$(53)

The evolution of the system is governed by Eqs. (48), (49), (50), (51), (53). When the plasma is unmagnetized (i.e., ω_L = 0), one would recover the evolution model of Sobacchi et al. (2024b).

4. Filamentation instability

Our final set of equations (i.e., Eqs. 48, 49, 50, 51, 53) have an exact solution, which is $\delta {\omega }_{\mathrm{L}z}=\delta {\mathbf{v}}_{\perp }=\delta j_z=\delta \rho =0$ and a = a₀ (where a₀ is a constant, which we assumed to be real). In the following, we show that this solution is unstable, and the wave breaks into filaments. We study the evolution of a small perturbation of the wave intensity by defining a = a₀(1 + δa). It is convenient to derive from Eq. (53) two equations for δa + δa^* and δa − δa^*, where δa^* denotes the complex conjugate of δa. Assuming that δa + δa^* and δa − δa^* are proportional to $exp[\mathrm{i}(\mathbf{K}·\mathbf{x}-\mathrm{\Omega }t)]$, one finds

$$\begin{array}{c}\hfill \frac{\mathrm{\Omega }}{{\omega }_{\mathrm{P}}}(\delta a+\delta a^{\ast })=\frac{K^2}{2{\omega }_{\mathrm{P}}^2}(\delta a-\delta a^{\ast }),\frac{\mathrm{\Omega }}{{\omega }_{\mathrm{P}}}(\delta a-\delta a^{\ast })=\frac{K^2}{2{\omega }_{\mathrm{P}}^2}(\delta a+\delta a^{\ast })+\delta \rho .\end{array}$$(54)

To determine the dispersion relation of the filamentation instability, one should use Eqs. (48)–(51) to express δρ as a function of δa + δa^*. Taking into account that |a|² = a₀²(1 + δa + δa^*) at leading order, one finds

$$\begin{array}{c}\hfill {(\mathrm{\Omega }+K_z)}^2\delta \rho =[K_x^2(1+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2})+\frac{K_y^2}{1-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{(\mathrm{\Omega }+K_z)}^2}}]\frac{a_0^2}{4{\gamma }_0^2}(\delta a+\delta a^{\ast }).\end{array}$$(55)

The dispersion relation can be determined from Eqs. (54)–(55). One finds

$$\begin{array}{c}\hfill {(\mathrm{\Omega }+K_z)}^2(\frac{4{\mathrm{\Omega }}^2}{K^2}-\frac{K^2}{{\omega }_{\mathrm{P}}^2})=\frac{a_0^2}{2{\gamma }_0^2}[K_x^2(1+\frac{{\omega }_{\mathrm{L}}^2}{{\omega }_{\mathrm{P}}^2})+\frac{K_y^2}{1-\frac{{\omega }_{\mathrm{L}}^2}{{\gamma }_0^2{(\mathrm{\Omega }+K_z)}^2}}].\end{array}$$(56)