© The Authors 2024

1. Introduction

Several coherent radio emission models have been proposed to explain the radiation of pulsars, magnetars, and fast radio bursts (FRBs; Ruderman & Sutherland 1975; Blaskiewicz et al. 1991; Kramer et al. 2002; Petrova 2009; Beskin 2018; Xiao et al. 2021). However, no agreement on a specific mechanism has been reached due to a lack of knowledge on the strongly nonlinear evolution of the relativistic plasma instabilities at kinetic microscales and disparaging theoretical predictions for their coherent emissions (Platts et al. 2019; Melrose et al. 2020; Zhang 2020).

The relativistic streaming instability is considered to be a likely cause of coherent radio emission in the magnetospheres of radio pulsars (Weatherall 1994; Zhang 2020) and has also been proposed as the source of radio emissions from magnetar bursts and in magnetar-like models of FRBs (Lu & Kumar 2018; Yang & Zhang 2018; Kumar & Bošnjak 2020; Yang et al. 2020; Lyutikov 2021a). At distances closer than the light cylinder radius, the streaming instabilities can be produced by pair creation cascades in relativistically hot and highly magnetized open magnetic field line regions above polar caps (Sturrock 1971; Ruderman & Sutherland 1975). There, the plasma cannot carry the magnetospheric Goldreich–Julian currents parallel to the magnetic field (Goldreich & Julian 1969), causing the convective electric field to grow to extreme values of ∼10¹² V m⁻¹. This field then accelerates electric charges up to relativistic energies, allowing pair creation avalanches, loosely termed “spark discharges”, to occur. The created dense pair plasma outflows from the spark region are modulated in time, forming trains of plasma bunches or clouds. The emitted bunches can expand as they propagate through the magnetosphere, and there are several conceivable reasons why they may induce streaming instabilities: (1) adjacent bunches can overlap in the phase space, (2) the electrons and positrons can have relative drifts, and (3) the bunches can create instabilities when they propagate through a background plasma (Usov 1987; Ursov & Usov 1988; Asseo & Melikidze 1998; Rahaman et al. 2020a; Manthei et al. 2021). The streaming instability may also be a source of single pulses from the outer magnetosphere, where the magnetospheric reconnection occurs and accelerates the particles to relativistic energies (Cerutti et al. 2012; Werner et al. 2018; Cerutti & Giacinti 2020; Hoshino 2020).

The formation of solitons on kinetic scales in relativistic streaming instabilities has been proposed by Novikov et al. (1984), Sircombe et al. (2005), and Henri et al. (2011). Several types of solitons have been considered in the neutron star magnetospheres (Karpman et al. 1975; Asseo & Melikidze 1998). First, soliton-like wave packets called cavitons (Zakharov 1972; Henri et al. 2011; Romero et al. 2016) can be formed in strong plasma turbulence created by the streaming instability (Eilek & Hankins 2016). In the turbulence, there is an unstable equilibrium between the plasma pressure gradient and the ponderomotive force (Mofiz et al. 1988; Mofiz 1989). This nonequilibrium leads to the creation of plasma cavities (cavitons) filled with a strong electric field and a charge separation inside the so-formed solitons. Second, force-free Alfvén soliton-like waves have been proposed (Mofiz 1993; Lyutikov 2021b). Two types of Alfvén solitons can be formed in the magnetospheres of neutron stars: large-scale waves of the magnetospheric size and waves generated by a fire hose type of current-driven instabilities (Pfeiffer & MacFadyen 2013; Gralla & Jacobson 2014, 2015; Lyutikov 2020). Third, a propagating beam can lead to the Weibel instability in the outer magnetosphere near a current sheet (Doǧan & Ekşi 2020). When the instability saturates, its nonlinear dispersion properties can induce the production of breathers, a localized periodic type of soliton, via self-modulation.

Plasma particles can, in some circumstances, be confined in solitons sufficiently long enough to behave collectively, thereby enabling the emission of coherent radiation. Such cavity type solitons are assumed to contain a large amount of electric charges. The emission power of the coherent emission mechanisms depends on a number of charged particles as N²; hence, the large number of plasma particles associated with the cavitons can be sufficient to explain the observed intensities (Melikidze & Pataraya 1980, 1984; Melikidze et al. 2000). It has been proposed that the radiation from solitons could form the strong single radio pulses in the outer magnetosphere of neutron stars (Main et al. 2021; Lin et al. 2023), or even super-giant pulses (Lyutikov 2007), that can be detected over large distances. The exact evolution, saturation, and the radio emission properties of cavitons in the neutron star magnetospheres are still being investigated (Eilek & Hankins 2016; Melrose & Rafat 2017).

Solitons can radiate via several conceivable emission mechanisms because of their strong nonlinearity and long time stability. The coherent curvature radiation of cavitons moving along curved magnetic field lines of pulsars is a widely considered radio emission mechanism of neutron stars (Weatherall & Benford 1991; Melikidze et al. 2000, 2014; Gil et al. 2004; Mitra et al. 2009; Mitra 2017; Rahaman et al. 2020b; Horvath et al. 2022).

The solitons can also radiate via linear acceleration emission in which the particles undergo acceleration parallel to the magnetic field (Melrose et al. 2009; Melrose & Luo 2009). However, the power of soliton linear acceleration emission was found to be too low to explain the observed radio power of pulsars (Benáček et al. 2023).

The Alfvén force-free soliton-like waves can be associated with free-electron maser or laser radiation as a model of FRBs (Lyutikov 2021a,b). The Alfvén waves can form wigglers for relativistic particles, and when a particle beam propagates through the wiggler, the plasma particles jiggle in the synchronism direction and radiate coherently. A related model has been proposed by Doǧan & Ekşi (2020) in which stimulated radio emission is amplified in klystron-like structures via relativistic bunches created by magnetic reconnection.

The soliton production in neutron star magnetospheres was recently seen in numerical simulations that solve the nonlinear Schrödinger equation in one dimension (Lakoba et al. 2018). In that approach, the solitons can develop if the Lighthill (1967) condition for the particle distribution is fulfilled, Gq > 0, where G is the group velocity parameter and q is the cubic nonlinearity parameter. The solitons were found to be formed with their amplitude growing to infinity because the approach did not consider the feedback of the soliton electric field on the particle distribution. Using that approach, Rahaman et al. (2022) showed that long-tailed particle distributions, such as the Kappa distribution, are more likely to produce the solitons for a wide range of plasma parameters.

We have studied the relativistic streaming instability using kinetic particle-in-cell (PIC) simulations that in previous publications were shown to evolve the particles and waves self-consistently (Manthei et al. 2021; Benáček et al. 2021a). These simulations, however, were restricted to one dimension (1D). They showed that cavitons can be formed if the plasma inverse temperature is ρ = mc²/k_BT ≥ 1.66 (where m is the electron mass, c is the light speed, k_B is the Boltzmann constant, and T is the thermodynamic temperature) or if the Lorentz factor of the beam particles is γ_b > 40. Nonetheless, the required conditions are more complex: the Lorentz factor of the beam can even be γ_b < 40 for cavitons to be created when inverse temperatures are higher than ρ = 1.66. We also studied how the streaming instability can be formed in 1D simulations by the overlapping plasma bunches created by spark events (Benáček et al. 2021b) and found cavitons with similar properties.

In this paper we present the first 2D PIC simulations of the relativistic streaming instability in the inner magnetosphere of a neutron star under the gyro-motion assumption typical for strongly magnetized pair plasma. We investigate the conditions under which the cavitons can be self-consistently formed or collapse in 1D and 2D plasmas. Although “soliton” is a more general term, we use it exclusively from now on with the understanding that it includes the caviton type of solitary waves.

This paper is structured as follows. Section 2 describes the numerical model and its initial 1D and 2D setups. In Sect. 3 we present a detailed analysis and comparison of soliton formation in 1D and 2D simulations and the dispersion properties of the resulting electrostatic waves. We discuss the results and implications for the interpretation of observations by solitons in Sect. 4 and summarize them in Sect. 5.

2. Methods

We selected three example cases, S1–S3, that we consider to be representative of the range of possible scenarios where initial plasma parameters lead to the formation of the streaming instability and Langmuir turbulence. The cases follow the three main scenarios of the bunch evolution that are described in the introduction: In S1, expanding plasma bunches start to overlap in the phase space and form the streaming instability. The streaming instability is characterized by relatively low temperatures and high beam Lorentz factors. In S2, a streaming instability formed between electrons and positrons inside the bunch, Here, the streaming instability is characterized by higher temperatures and lower beam Lorentz factors than in case S1. S3 represents a plasma bunch propagating through a background plasma. The streaming instability is characterized by relativistic temperatures and high beam Lorentz factors.

We studied the formation of nonlinear electrostatic solitary waves by the relativistic streaming instabilities in pair pulsar plasma at kinetic microscales using fully kinetic PIC numerical code ACRONYM (Kilian et al. 2012) with the same initial parameters but in separate runs for the 1D and 2D versions of the code.

In all cases, the code is implemented as 1D3V and 2D3V (one or two spatial dimensions and three-dimensional velocity distribution). Nonetheless, because the particles move strictly along the magnetic field direction, the code is effectively 1D1V and 2D1V (one-dimensional velocity distribution) with the velocity in the direction of the uniform external magnetic field. We utilized the Esirkepov (2001) current-conserving deposition scheme with piecewise quadratic shape function of macroparticles. The piecewise quadratic shape function describes the plasma dispersion properties of the relativistic plasma well according to our tests. The code also contains a higher-order shape function “Weighting with Time-step dependency” (Lu et al. 2020) that we previously used (Manthei et al. 2021; Benáček et al. 2021a,b), but according to our tests, it does not provide any better suppression of the numerical Cerenkov radiation artifacts but consumes more computing time. The simulation setup is summarized in Table 1.

2.1. Numerical methods

As the particles move in strongly magnetized plasmas, they lose their perpendicular momenta by synchrotron radiation in less than one simulation time step. Therefore, we modified the particle pusher to assume that the particle can move only along the magnetic field in a gyro-motion approach. In the approach, the electric field vector is separated into components that are parallel and perpendicular to the magnetic field:

$$\begin{array}{c}\hfill \mathit{E}={\mathit{E}}_{\Vert }+{\mathit{E}}_{\perp }.\end{array}$$(1)

We used the Vay et al. (2011) particle pusher in the gyro-motion (1D velocity space) approximation, which is an energy-conserving ultra-relativistic generalization of the Boris push (Boris 1970). The pusher considers only the parallel component E_∥ to push the particle, setting the perpendicular velocity of every particle to zero, v_⊥ = 0.

We assumed an initially uniform magnetic field of B = 10¹² G directed along the x-axis in the simulation domain. The gyro-motion approximation in the simulation can be considered valid if most of the perpendicular kinetic energy of a particle is radiated faster than in one simulation time step. Our checks using computed radiative energy losses found that the condition is fulfilled for all particles in the simulation if the magnetic fields are ≳10⁸ G. Assuming a dipolar field with an intensity of 10¹² G at the star surface, the approximation is valid at least up to a distance of ∼22 R_NS from the star, where R_NS is the neutron star radius. We note that the applicability of the gyro-motion approximation depends on particle velocity, and the transition between the gyro-motion approach and the requirement for a full 3D velocity calculations is spread over a few orders of magnetic field intensity ≲10⁸ G.

To ensure the numerical precision of a perturbation of the magnetic field produced in the simulation, we separated the magnetic field into a constant external component, B₀, and a perturbed component, δB:

$$\begin{array}{c}\hfill \mathit{B}(x,y,t)={\mathit{B}}_0+\delta \mathit{B}(x,y,t),{\mathit{B}}_0\gg \delta \mathit{B}(x,y,t).\end{array}$$(2)

While the total magnetic field B is exploited in the particle pusher, the Maxwell solver uses only the component δB(x, t) as ∇ × B₀ = 0 and ∂B₀/∂t = 0. The external magnetic field is along the x-axis.

The main numerical challenge in the 2D electromagnetic simulations is the suppression of the numerical Cerenkov radiation appearing as a numerical enhancement of high-frequency waves. We used the modified second-order time-domain and fourth-order finite-difference scheme (M24) proposed by Hadi & Piket-May (1997) and Greenwood et al. (2004) for computation on the Yee lattice to suppress these numerical artifacts. In addition, we combined the field solver with the Friedman (1990) low-pass filtering method with a Friedman parameter θ = 0.1. This value assures that only the high-frequency waves are suppressed while the waves with frequencies around and at a few multiples of plasma frequency are maintained without significant modification. We also tested Cole–Kärkkäinen (CK; Cole 1997; Kärkkäinen & Gjonaj 2006) field solvers with various parameters, which were stable in 1D electrostatic simulations (Benáček et al. 2021a,b). The CK solvers, however, failed to suppress the numerical Cerenkov radiation in electromagnetic 2D simulations while keeping the waves around and above plasma frequency unchanged, even with the applied Friedman filter.

2.2. Simulation setup of the streaming instabilities

The plasma is composed of background and beam particles, initially having the same number of particles per cell n₀ = n₁ = 100. We assumed both components consist of electrons and positrons, initially in a charge equilibrium in each grid cell. The particles have an initial Maxwell–Jüttner velocity distribution function (Jüttner 1911),

$$\begin{array}{cc}\hfill g(u)& =g_0(u)+g_1(u),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill g_0(u)& =\frac{n_0}{2K_1(\rho )}{\mathrm{e}}^{-\rho \gamma },\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill g_1(u)& =\frac{1}{{\gamma }_{\mathrm{b}}}\frac{n_1}{2K_1(\rho )}{\mathrm{e}}^{-\rho {\gamma }_{\mathrm{b}}\gamma (1-\beta {\beta }_{\mathrm{b}})},\hfill \end{array}$$(5)

where u/c = γβ = γv/c and u_b/c = γ_bβ_b = γ_bv_b/c are the four-velocity and velocity of the beam, respectively, both along the x-axis and normalized to the speed of light c; $\gamma ={(1-{\beta }^2)}^{-\frac{1}{2}}$ and ${\gamma }_{\mathrm{b}}={(1-{\beta }_{\mathrm{b}}^2)}^{-\frac{1}{2}}$ are their Lorentz factors, K₁ is the MacDonald function of the first order (modified Bessel function of the second kind), and ρ is the inverse temperature. We assumed that there is no net electric charge and no net electric current in the simulations, or in other words, the magnetospheric Goldreich–Julian currents and charges are zero or negligibly small compared to the bunch density.

To the best of our knowledge, there are no specific predictions regarding which plasma parameters produce most of the radio emission. We therefore followed the estimates of (Arendt & Eilek 2002) for the sparking mechanism and chose the bunches to have an inverse temperature around ρ ≈ 1 before they evolve. This choice appears reasonable in the light of our earlier results for case S1. When the bunches overlapped in the phase space, the distribution tended to narrow and evolve to higher inverse temperatures of ρ ≈ 3 (Benáček et al. 2021b).

The inverse temperature ρ ≈ 2 in the case S2 was estimated in broad parametric studies, for example, Rahaman et al. (2020a, 2022) to fulfill the Lighthill condition. The instability then has a high nonlinearity factor in the Lighthill condition and consequently produces solitons. Similar inverse temperatures can, in addition, be found for highly overlapped plasma bunches. For S3, we assumed a plasma bunch with an inverse temperature ρ ≈ 1 again but that the bunch propagates through a background plasma of the same temperature.

In simulations S1 and S3, we used the beam Lorentz factors for which the largest linear growth rates are obtained and, thereby, expected to produce the highest amount of plasma turbulence for the given temperature. That provides the “best” conditions for soliton formation. Although not rigorously proven, we consider it natural to expect that the largest growth rates of the electrostatic waves will create the strongest instability, the strongest turbulence level, the largest soliton charges, and the largest electromagnetic emission in comparison with the instability with lower linear growth rates. Benáček et al. (2021a) showed that the largest (1D) growth rates occur for beam Lorentz factors γ_b ≈ 80 for ρ = 3 and γ_b ≈ 60 for ρ = 1. The beam Lorentz factor of γ_b = 6.24 for simulations S2 was selected to fulfill the Lighthill condition, as found by Rahaman et al. (2022). That was shown to have a high nonlinearity parameter and to produce solitons in 1D numerical simulations. For simplicity, we assumed the same inverse temperature for the background and beam plasma in all cases.

The grid length of L = 8192Δ_x along the magnetic field is the same for all simulations, where Δ_x = Δ_y is the grid cell size. The 2D simulations have 8192Δ_y in the perpendicular direction. The grid length L allows the highest wave number resolution in each dimension up to Δkd_e ≈ 7.7 × 10⁻³ in S1 and S3 and Δkd_e ≈ 3.8 × 10⁻³ in S2.

The simulation time step is chosen as ω_pΔt = 0.05 for simulations S1 and S3 and ω_pΔt = 0.045 for S2. The highest frequency resolution in the simulation is Δω/ω_p ≈ 1.0 × 10⁻⁴ in all simulations. The boundary conditions are periodic in both directions. Grid cell size and the time step length differ between simulations because we selected a larger grid cell size and slightly smaller simulation time step in S2 to describe better the produced electrostatic waves in the ω − k space for smaller γ_b.

3. Results

3.1. Evolution of the electric field

The evolution of the electric energy W_E = W_Ex + W_Ey in our six simulations is compared in Fig. 1. The energy is always normalized to the initial kinetic energy W_k of the plasma in the given simulation. The energy components $W_{\text{E}x}=\int E_x^2/8\pi \text{d}V=\int w_x\text{d}V$ and $W_y=\int E_y^2/8\pi \text{d}V=\int w_y\text{d}V$ of the electric field components E_x and E_y are then shown for comparison. The energy ratio evolves in the typical manner of a streaming instability: the instability saturates due to nonlinear effects after an initial exponential growth phase.

Fig. 1. Comparison of the energy evolution of electric field components parallel, Ex, and perpendicular, Ey, to the magnetic field between 1D (black lines) and 2D (blue lines) simulations. All energies are normalized to the initial kinetic energy. The parallel components of the energy ratio Wx/Wk are always smaller in 2D than in 1D and decreasing. The 2D simulations grow later because they contain a lower initial noise level than the 1D ones at the simulation start. (a) S1–1D + S1–2D. (b) S2–1D + S2–2D. (c) S3–1D + S3–2D.

The instabilities in the 2D simulations, however, start to grow and saturate later than for 1D simulations because they contain a lower initial noise level, due to a larger number of macro-particles in the Debye volume, from which the instability grows. This effect is not caused by the different size of the simulations as the electric fields in Fig. 1 are always normalized to the initial kinetic energy of the plasma. Though the 2D simulations have 8192 times more grid cells, they contain a higher amount of initial kinetic energy by the same factor. Therefore, the ratio W_E/W_k is comparable between the 1D and 2D cases.

We find that during the whole simulated time interval, the W_Ex/W_k energy ratio of the 2D simulations is always smaller than the W_Ex/W_k energy ratio for 1D simulations. Moreover, while the 1D energy evolution stabilizes after the saturation, the 2D simulations show an exponential decay toward a lower equilibrium value. In the 1D simulations, the W_Ey energy is always zero, as by definition, there is no electric current nor electric field contribution to this component in the 1D1V approximation. The W_Ey energy is plotted only for clarity. However, the W_Ey energy is nonzero in 2D simulations. In simulation S1–2D, that component saturates at lower energy, but it later exceeds the W_Ex energy, which decreases faster. In simulations S2–2D and S3–2D, the W_Ey energy is the same or greater than the W_Ex one. Both components grow and saturate together in 2D simulations; hence, their evolution can be considered to be strongly connected.

We estimated the growth rates, Γ, of all simulations and the decrease rates, δ, of the 2D simulations separately for both amplitudes of the electric field components (see Table 2). The rates were obtained using an exponential fit of the energy ratio evolution profile, and they are normalized to the plasma frequency, ω_p. The time windows of the highest and lowest rates are selected for the growth and decrease, respectively, in each simulation. Because the E_y component equals zero for all 1D simulations, only the parameters for the E_x component are estimated. The growth rates of the energy ratios are assumed in the form of W_E/W_k ≈ W₀e^2Γt + W₁, where W₀, W₁, and Γ are the parameters for the fit. The dissipation rates are assumed similarly W_E/W_k ≈ W₀e^−2δt + W_asymp, where W_asymp is an asymptotic value to which the energy ratio approach after the saturation stage.

The growth rate values of W_Ex/W_k ratio are comparable between 1D and 2D simulations. The simulation S2–2D has the highest growth rate in both W_Ex and W_Ey. While the simulations S2–2D has highest dissipation rate in W_Ex/W_k ratio, S1–2S features, in contrast, the highest dissipation rate in W_Ey/ε_k ratio. The asymptotic energies of both components of 2D simulations are significantly lower than their saturation energies.

3.2. Identification of wave modes

Figures 2 and 3 show dispersion diagrams of the electrostatic waves. Specifically, the E_x component along the magnetic field line is chosen because this component is common for the 1D and 2D simulations.

Fig. 2. Dispersion diagrams of the electrostatic waves along the magnetic field around the saturation time. Comparison between 1D (top row) and 2D (bottom row) simulations with the same plasma parameters. The solitons dispersion branches are denoted by black arrows and the superluminal L-mode waves by red arrows. The superluminal waves appear in 1D and 2D, while the soliton branch appears only in 1D. (a) S1–1D : ωpt = 1500 − 2500. (b) S2–1D : ωpt = 500 − 1500. (c) S3–1D : ωpt = 4000 − 5000. (d) S1–2D : ωpt = 2000 − 3000. (e) S2–2D : ωpt = 1000 − 2000. (f) S3–2D : ωpt = 3000 − 4000.

Fig. 3. Same as Fig. 2 but at the end of simulations in a time interval ωpt = 9000 − 10 000. Solitons do not appear in simulations (b) and (d)–(f). While the most intensive waves in 1D are the solitons and superluminal L-mode waves close to k∥ ≈ 0, in 2D, the most intensive waves are close to the light line, and solitons are not present. (a) S1–1D. (b) S2–1D. (c) S3–1D. (d) S1–2D. (e) S2–2D. (f) S3–2D.

The limited storage capacity of the SuperMUC-NG cluster did not allow us to store the field data for the whole simulation time. We therefore chose a time interval of ω_pδt = 1000 based on the saturation time for the dispersion diagrams, and the similar time intervals were used in Fig. 2 for the saturation regime and in Fig. 3 for the end of the simulation at ω_pt = 9000 − 10 000.

All subplots of Figs. 2 and 3 show longitudinal L-mode electrostatic waves indicated by red arrows. In contrast, we associated the solitons with horizontal branches of the superluminal electrostatic waves. The superluminal L-mode waves can be obtained as solutions of the linearized relativistic dispersion permittivity tensor in the gyro-motion approximation, Λ₃₃ = 0 (Rafat et al. 2019). The soliton branch is, however, a nonlinear solution of the dispersion tensor. The soliton waves form straight, almost horizontal lines in the ω − k frequency–wave number space for kc/ω_p ≲ 1 and ω ≲ ω_p (where ω_p is the plasma frequency), thereby not following the Langmuir dispersion branch. A more detailed analysis of the visible horizontal dispersion features based on the inverse Fourier transformation to real space of the horizontal branch area in the ω − k domain showed that they correspond to stable, soliton-like waves (Benáček et al. 2021a). Furthermore, we therefore identify the visible horizontal structures as the soliton dispersion branch in the rest of the paper. Because the dispersion branch is horizontal in the soliton reference frame, all the waves associated with the soliton oscillate with the same frequency but with different k in that frame. Their group velocities are also the same as their dispersion branch, which has a constant slope. Hence, the soliton branch can be described by a dispersion relation ω(k) = ω_s, where ω_s is the soliton oscillation frequency. The soliton dispersion branches are indicated by black arrows in Figs. 2 and 3.

3.3. Analysis of the simulation results and the evolution of plasma waves in different scenarios

In simulations S1–1D and S3–1D, almost horizontal straight branches can be identified closely below the cutoff frequency of the superluminal L-mode waves. That, and the broad horizontal extent in the ω − k diagram, indicates that there have to be correspondingly narrow spatial features with very low group velocity that are repetitive and spiky in time, which is typical for plasma solitons.

The simulation S2–1D shows the soliton branches during the saturation in Fig. 2 but not at the simulation ends in Fig. 3. The soliton branches in S2–1D are different from S1–1D and S3–1D; mostly subluminal and located more below the L-mode cutoff frequency, typically at lower frequencies around ω = (0 − 0.3)ω_p. We find that the main factor influencing the difference is the Lorentz factor of the beam, which is significantly lower in S2 than in the other two cases. Moreover, compared with the 1D simulations, the 2D simulations show no evidence of a soliton branch neither at saturation time nor at the simulation end. Most energy at the end of simulations S1–2D and S3–2D is present at the L-mode branches in a region of the branch close to the light line ω = ck. In the simulation S2–2D, approximately the same amount of energy is present in the superluminal and the subluminal L-mode branches. However, we find that the superluminal waves continue to dissipate slowly at the simulation end.

In Figs. 4 and 5 we analyze the spatial profiles of the electrostatic energy density w_E = |E|²/8π normalized to the initial kinetic energy density w_k separately for 1D and 2D simulations. The 2D simulations are selected as a cut along x at y = L/2 axis, where L is the size of the simulation domain, and bicubic interpolation is applied in all plots. In 1D simulations, electric fields were stored on disks every 20th time step, but in 2D simulations, we stored the fields only every 200th time step in S2–2D and S3–2D and every 1000th time step in S1–2D because of storage limitations. Therefore, the structures in Fig. 4d are more poorly resolved than in the other panels. The soliton waves appear during the instability saturation as the most intensive red regions in the S1–1D and S3–1D simulations in Fig. 4, and we associated these waves with the soliton branches in Figs. 2 and 3. Although the solitons survive until the end of the simulations in S1–1D and S3–1D, they are not very visible in Fig. 4b because their intensity is smaller than the intensity of the superluminal waves, which can be seen from Fig. 2b. The superluminal waves decrease quickly in S2–1D and S2–2D after their saturation approximately in ω_pt ∼ 1000. In the S2–1D simulation, only the L-mode waves remain, seen as almost horizontal orange-red lines moving with group velocities close to the light speed. The simulations S1–2D and S3–2D also show the rapid decrease in electrostatic energy.

Fig. 4. Evolution of the electrostatic energy density in simulations normalized to the initial kinetic energy density. Although in S1–1D and S3–1D, the solitons survive until the simulations end, in S2–1D, the solitons disappear before ωpt = 2000. No solitons appear in 2D simulations. (a) S1–1D. (b) S2–1D. (c) S3–1D. (d) S1–2D. (e) S2–2D. (f) S3–2D.

Fig. 5. Electrostatic energy density of the 2D simulations normalized to the initial kinetic energy density in the x − y plane during the saturation times (top row) and at the simulation ends (bottom row). The evolution is available as an online movie. Only a simulation sub-domain is selected to show the structures formed by the waves. (a) S1–1D : ωpt = 2500. (b) S2–1D : ωpt = 900. (c) S3–1D : ωpt = 3000. (d) S1–2D : ωpt = 10 000. (e) S2–2D : ωpt = 10 000. (f) S3–2D : ωpt = 10 000.

We plot the electrostatic energy density in 2D simulations in Fig. 5 during the saturation times and at the simulation ends. One quarter of the simulation domain is shown to better resolve the electric field structures in the plot. During the saturation time, more irregularly red-like shaped regions are present in Fig. 2. We find the regions associated with the superluminal L-mode waves in dispersion diagrams of Figs. 2d and e. The vertical green–yellow stripes along y-axis that are narrow in x in Figs. 5c, d and f as well as the small green-yellow regions in Fig. 5e are associated with the L-mode waves close to the light line.

Figure 6 compares the wave vector space of longitudinal electrostatic waves in the 2D simulations at saturation time and simulation end. The longitudinal electric component for a grid cell (i, j) in the k_x − k_y space is obtained as E_long(i, j) = k_xE_x(i, j)+k_yE_y(i, j) for each wave vector k = (k_x, k_y). The wave energy density is calculated as ${\mathit{w}}_{\mathrm{long}}(k_x,k_{\mathit{y}},\omega )=\frac{1}{8\pi }{|E_{\mathrm{long}}(k_x,k_{\mathit{y}},\omega )|}^2$ where E_long is in Gaussian CGS units of the simulation. The waves are obtained in ω_pt = 1000 long intervals and summed over frequencies in a range ω = (0 − 2)ω_p. While in S1–2D and S2–2D simulations, the wave profiles significantly change after the instability saturation, only a small narrowing of the waves appears along the x-axis in S3–2D. Moreover, in S1–2D and S3–2D, the waves have directions mostly along the magnetic field at the simulation end. In S2–2D, the most intense longitudinal waves are directed perpendicular to the magnetic field. These longitudinal waves along the y axis are associated with the intensive E_y component (blue dotted line) in Fig. 1b. We note that the longitudinal waves propagating in the perpendicular direction to the magnetic field are described by the E_y component of the electric field in our configuration. A deeper analysis in the ω − k_x − k_y space, revealed that these perpendicular waves have frequencies ω < 0.01ω_p for |k_yc/ω_p|≳0.5 and ω < 0.05ω_p for |k_yc/ω_p|≲0.5.

Fig. 6. Energy density of longitudinal electrostatic waves of 2D simulations in the wave vector space during the saturation time intervals (top row) and at the simulation ends (bottom row). All analyzed intervals are ωpta = 1000 long. The comparison between rows shows that the longitudinal electrostatic waves close to |k| = 0 decrease in time. (a) S1–2D : ωpt = 2000 − 3000. (b) S2–2D : ωpt = 1000 − 2000. (c) S3–2D : ωpt = 3000 − 4000. (d) S1–2D : ωpt = 9000 − 10 000. (e) S2–2D : ωpt = 9000 − 10 000. (f) S3–2D : ωpt = 9000 − 10 000.

3.4. Evolution of particle distributions

Figure 7 shows the particle distributions γf(p_∥), where f(p_∥) is the particle momentum distribution, at the start of the simulation, at instability saturation time, and at the end of the simulation for 1D and 2D simulations. All simulation macro-particles (electrons and positrons) are included, and local variations in the distributions are averaged. The distributions are normalized as ${\int }_{-\infty }^{\infty }\gamma (p_{\Vert })f(p_{\Vert })\mathrm{d}{p}_{\Vert }=1$. Only the momentum intervals with significant evolution of the distributions are shown. The beams in 1D simulations have a relatively small number of particles, and their distributions contain a high level of noise; the beam parts of these distributions were, therefore, smoothed using the Gaussian filter with a width of four data points.

Fig. 7. Comparison of particle distributions as functions of the particle momenta. All electrons and positrons are used from the simulations. The distributions for 1D and 2D simulations at the simulation start, instability saturation, and the simulation end (top row). Comparison of the distribution difference between the simulation end and saturation time for 1D and 2D simulations (middle row). Comparison of the distribution difference between 1D and 2D simulations for the saturation time and simulation end (bottom row). The vertical dotted gray lines show selected positions of momenta where the distributions change significantly. The dashed gray lines show the initial distributions (top row) and a zero distribution difference (middle and bottom rows). (a) S1–1D + S1–2D. (b) S2–1D + S2–2D. (c) S3–1D + S3–2D.

In all simulations, the initial distributions relax during the linear evolution to the saturation time, generating electrostatic waves (Fig. 1) and thereby closing the gap between the beam and background populations (Figs. 7a–c). Later on, during the saturated stage of the instability, the distribution functions evolve slowly to an asymptotic state with profiles in the velocity space between the initial background and beam populations. We find that the saturation and asymptotic states of the 2D distributions have steeper slopes in comparison with the 1D cases. This implies that more kinetic energy is converted to waves in the 1D simulations than in the 2D simulations, which agrees with the behavior seen in Fig. 1, except for the W_Ey component of S2–2D, which may not be directly related to the average velocity distribution along the magnetic field lines but rather to spatial density fluctuations perpendicular to the magnetic field. Though the 1D and 2D distributions are similar at saturation times, their small differences are significant enough to explain the differences in the evolution of the electrostatic wave energy. The change is caused by a redistribution of the kinetic energy over the distribution function toward smaller momentum gradients of the distributions in regions with positive, p_∥∂f/∂p_∥, resulting in higher plasma stability.

Figures 7d–f show the differences in distributions between the saturation time and the simulation end. There are regions in momenta where the distributions decrease and regions where they increase, respectively. The largest values of the differences are for the background plasma momenta |p_∥/(m_ec)| ≲ 1; however, these differences are not associated with large energy changes because the differences are restricted to very narrow momentum space intervals. Between saturation time and simulation end, the distributions evolve mainly in regions where the particle distributions (1) increase for momenta between the background and beam distributions, (2) decrease where the distribution has a positive gradient, and (3) increase in the beam tail where the distribution has a negative gradient. The momentum distribution gradients become flatter in time in all cases.

The distributions differ between the 1D and 2D simulations as shown in Figs. 7g–i. The largest differences are for beam momenta where the distributions have positive gradients. The distribution functions of 2D simulations have systematically higher values than in 1D simulations in this region. For higher momenta where the distributions have negative distribution gradients for beam particles, the distribution functions of 2D simulations have smaller or similar values. The differences between the background distributions of 1D and 2D simulations are negligible in S1 and S3 simulations because the distributions are very narrow in momentum space. In the S2 simulations, the negative momentum distribution tail p_∥/(m_ec)≲ − 3 has higher values for 1D simulation, for momenta p_∥/(m_ec)≈ − 3 to 1 is higher for 2D simulations, and for momenta p_∥/(m_ec)≈1 − 4. As a result, the 1D distributions have systematically smaller momentum gradients than the 2D ones.

4. Discussion

We studied the formation and evolution of cavitons, soliton-like waves, as possible sources of coherent pulsar radio emissions by utilizing 1D and 2D PIC simulations of the relativistic streaming instability in strongly magnetized pair plasma of neutron star magnetospheres.

Focusing on the instability evolution during and after saturation, we find that the solitons, known to form a horizontal dispersion branch in the 1D simulations, do not appear in 2D simulations because of the additional degree of freedom perpendicular to the magnetic field in the plasma. Only strong superluminal L-mode waves on the relativistic Langmuir-like dispersion branches appear in 2D simulations during the saturation of the instability. These waves decrease on timescales of a few hundred or thousand plasma periods. Moreover, in one of the 1D simulations, the solitons amplitudes were also found to decrease, and the solitons collapse after some time. Solitons already having a somewhat precarious existence in 1D simulations, they cannot be considered a realistic contender for radio emission in strongly magnetized pulsar plasma based solely on the 1D approach. Soliton-like waves do not appear at all in the 2D simulations, a finding that agrees with the already described soliton collapse for unmagnetized or weakly magnetized kinetic plasmas, where the solitons may initially appear but collapse quickly when their wavelengths increase during their evolution (Zakharov & Shabat 1972; Robinson 1997).

In contrast to the 1D soliton stability, there is no stable solution of the nonlinear Schrödinger equation in more than one dimension for space and laboratory plasmas (Robinson 1997, and references therein). The Langmuir solitons in electron–ion plasma were shown to undergo collapse and dissipation when two-dimensional (2D) plasmas are considered, which was shown by solving the set of Zakharov equations and also confirmed by numerical PIC simulations for weakly magnetized plasmas (Newman et al. 1990). The main difference between our and the classical approach based on solving of the Schrödinger equation and Zakharov equations is that the PIC simulations describe the particle–wave interactions self-consistently. Specifically, the PIC simulations consider the evolution of the particle distribution as a function of space and time. And exactly these changes of the distribution function may be the nonlinear effects that suppress the formation of the solitons in the 2D approach. We show, in addition, that the nonlinear soliton dispersion forms a straight horizontal branch in the frequency–wave number space in the 1D approach. However, the analytical description is based on a soliton dispersion that follows the quasi-linear dispersion relation of Langmuir- or L-modes. However, that dispersion branch does not contain the superluminal or subluminal L-modes that are also present.

A soliton, being a plasma structure that affects the electric field and the particle distribution function, can only be influenced by electric fields along the magnetic field in 1D simulations. The physical reason for the absence of solitons in the 2D simulations is connected with the additional plasma effects in directions perpendicular to the magnetic field. In the soliton reference frame, the soliton dispersion branch is horizontal, which agrees with expected soliton stability as all waves forming the soliton oscillate with the same frequency and have the same group velocity. The solitons are assumed infinitely large in the y direction in 1D, but their true perpendicular size in 2D or 3D remained as an open question. In the 2D model, particles located perpendicular to the locally enhanced electric field of the soliton progenitor are free to also interact with their surroundings in the perpendicular direction to the magnetic field and the electric field of the wave can spread out and decrease. These superluminal L-mode waves have been observed to decrease faster in S1–2D and S3–2D than in S2–2D, where the perpendicular size of the waves is smaller. Figure 5 must be seen as supporting evidence that the effects of large perpendicular size contribute significantly to the wave decay.

The particle distributions simultaneously stabilize when the electrostatic wave energy in the 2D simulations decreases after saturation (see Figs. 1 and 7). The momentum gradients of averaged particle distributions, p_∥∂f/∂p_∥, decrease in time, and the instability growth rates approach zero in all 1D and 2D simulations. The 1D distributions have systematically smaller gradients of the distributions at the saturation time than the 2D ones, implying higher energy conversion in 1D simulations than in 2D ones.

The conversion of the electrostatic energy after the instability saturation occurs during the nonlinear evolution phase of the instability. There are two reasons why it is difficult to find direct evidence of the electrostatic energy conversion in the average particle distribution profiles. First, although the parallel (to the magnetic field) subluminal component of electric waves can be converted into particle kinetic energy, the energy is small compared to the total energy that is simultaneously being redistributed between various parts of the particle distribution. Second, the presented particle distributions are spatially averaged, but part of the electrostatic energy is associated with the local density fluctuations that are suppressed by such averaging. The amplitudes of these density fluctuations also decrease after the instability saturation, and that also decreases the amplitudes of the parallel superluminal electrostatic and perpendicular electrostatic waves. There cannot be direct energy conversion of the parallel superluminal waves because wave–particle interaction can only occur for subluminal wave speeds. Direct conversion of the perpendicular wave energy to particle motions is also ruled out because the particles move strictly along the magnetic field, and the electric field is perpendicular to their velocity vector. Both types of waves are instead converted nonlinearly by the fluctuations distributed along and across the magnetic field.

The Lorentz factor of the beam influences the directivity of the longitudinal waves at the end of the 2D simulations (Fig. 6). The longitudinal waves are directed closely along the magnetic field for high values γ_b ≳ 60, but low frequency waves are enhanced mainly in the y direction for lower γ_b ≈ 6 (S2–2D). These waves have higher intensities than the waves along the magnetic field in the same simulation, and the wave amplitudes slowly decrease.

We find that the growth rates for L-waves do not differ strongly between 1D and 2D simulations. As was shown by Manthei et al. (2021), the growth rates obtained from the 1D PIC simulations are comparable with those obtained from a linearized analytical 1D approach. The higher growth rates of W_Ex compared to W_Ey are probably caused by nonlinear coupling between both components. The decay rates are, however, lower than the growth rates in all three 2D simulations.

The 1D and 2D comparison indicates that the gradients of ∂f(x, y, v)/∂y ≠ 0 and ∂E(x, y)/∂y ≠ 0 may be important for the suppression of soliton formation. These gradients are equal to zero in the 1D approach, ∂f(x, y)/∂y = 0 and ∂E/∂y = 0, but are allowed to vary in 2D.

Some laboratory results appear to show the formation of soliton-like waves in (3D) experiments of weakly magnetized plasmas, in contrast to our findings of suppression of solitons in 2D (Borghesi et al. 2002; Sircombe et al. 2005; Sarri et al. 2010; Henri et al. 2011). There have also been in situ space measurements of solitary structures in the magnetospheres of solar planets exhibiting shapes similar to Langmuir solitons (Pickett et al. 2015; Píša et al. 2021). However, after a detailed consideration of these findings, we cannot decisively say that these experimentally and observationally detected waves are soliton-like cavitons or rather a special kind of “linear” Langmuir wave with similar properties. We show in Figs. 2 and 3 that the solitons and Langmuir L-mode waves can be identified by their specific features in the ω − k_∥ dispersion diagram. These were, however, not provided for the mentioned laboratory and planetary observations.

4.1. Implications for the pulsar radio emission

Although the solitons are known to collapse and disappear, they are nevertheless often considered as sources of pulsar radiation because the plasma dispersion properties of extremely magnetized neutron star magnetospheres are conjectured to be well described by the 1D approach.

As the solitons do not appear in 2D because of the existence of additional degrees of freedom for the particles, their formation therefore also is not expected in 3D. When superluminal solitons cannot form for any reasonable set of plasma parameters, then solitons cannot be considered as an essential element of coherent curvature or any other emission mechanism in the inner neutron star magnetospheres. Nevertheless, our findings cannot absolutely exclude the possibility of soliton formation via other plasma kinetic instabilities, and as we have discussed before, the 3D consideration may, however, further aggravate the described formation difficulties.

The absence of solitons does not constrain other possible radiation mechanisms of the streaming instability. The streaming instability still forms strong local plasma charge variations that we found to be associated with the electrostatic waves. Because the plasma charge density follows the positions of the electric field profile in Figs. 5d–f, other strong plasma waves might also lead to coherent curvature radiation.

Only point-like plasma charge concentrations have been considered in the past for predictions of the radio emission properties of solitons. Realistic electric charge structures are, however, far from the point-source approximation usually considered in analytical calculations (Melikidze et al. 2000). The shapes and sizes of the charge structures are important as they significantly influence the observable radio emission patterns Benáček et al. (2023). The structures that formed along the perpendicular direction to the magnetic field, as shown in Fig. 5, can in principle be expected to form strong field aligned curvature radiation patterns as individual radiating elements of the charge distribution may add up positively in the field direction. That alignment will then also enhance the observable radiative flux density. A more detailed knowledge of the spatial structure of solitons would allow a more self-consistent description of their distinctive radio emission properties, especially the directivity, frequency spectrum, and coherency of the produced waves.

The relativistic plasma emission, linear acceleration emission, or the free electron maser were also suggested as alternatives to coherent curvature emission in pulsars (Melrose et al. 2020). The linear acceleration emission, however, is not effective enough because its emission power was found to be too low (≈10¹⁶ erg s⁻¹) to explain the total pulsar radio power (10²⁵ − 10²⁹ erg s⁻¹; Benáček et al. 2023). The specific radio emission properties of the other mechanisms are not well known because they are strongly nonlinear wave–particle processes. Strongly overlapped bunches have, for example, a strong E_y component with a group velocity close to zero, oscillating at frequencies close to zero. If relativistic plasma particles travel through this wave along a magnetic field, the particles will oscillate synchronously and produce radiation by the free electron maser/laser effect (Lyutikov 2020, 2021b). The amplification of free electron maser waves is, nonetheless, constrained by the curvature of the trajectory and the dissipation timescale, and it is still far from certain that it may be responsible for pulsar radio emission.

4.2. Extrapolation of results to 3D and weaker magnetic fields

It is an open question how the instability evolves in 3D space and in a weaker magnetic field; this question can only be resolved using fully fledged 3D simulations. In 2D, the wave energy decrease can be understood as energy flowing from the superluminal wave “volume” (a region where most of the soliton electrostatic energy is located) to a surrounding plasma through a surface of the volume. However, in 3D, the ratio between the surface and volume is even larger than in 2D, and the decrease rate can be expected to be higher.

Our simulations are consistent with magnetic fields ≳10⁸ G, which is valid up to ∼22 R_ns for neutron stars with a field intensity of 10¹² G. Nevertheless, for millisecond pulsars, the condition may only be fulfilled at heights smaller than one star radius above the stellar surface. The plasma particles in weaker magnetic fields, however, have higher pitch angles with respect to the magnetic field, which acts against the soliton formation.

Bret et al. (2010a) found that the streaming instability evolves in a significantly different manner in 2D compared to 1D for a weakly relativistic, cold, electron–ion plasma with weak magnetic fields. They found that oblique modes appear to grow and filamentation instabilities also grow and produce micro-currents in the perpendicular direction, and particles can be accelerated by the Weibel instability (Nishikawa et al. 2003; Bret et al. 2006, 2008, 2010b). In addition, they estimated a typical perpendicular size of the modulated plasma in 3D simulations, finding a value ∼10 d_e. In our 2D simulations, the typical perpendicular size of the L-mode waves is on the order of ∼100 d_e for high Lorentz factors.

5. Conclusions

Streaming instabilities in the magnetospheres of neutron stars were theoretically proposed and later seen to form nonlinear caviton soliton-like waves in various numerical 1D simulations. The instability can be formed because of plasma bunches overlapping in phase space, because of a drift between electrons and positrons, or by bunches propagating through a background plasma. Therefore, we assumed that the 1D approximation is valid and useful since the approach can describe an instability evolution in a strongly magnetized plasma. As solitons can host a large electric charge, they can be exploited to explain the intense pulsar radiation due to, for example, the coherent curvature mechanism.

In this paper we have compared the 1D and 2D PIC simulations of streaming instabilities using simulations for three typical parameter scenarios of the streaming instability based on analytical and numerical investigations that were known to form solitons in 1D. We did not find any indication that solitons could be formed in 2D for those or any of the parameter sets under consideration. It appears that us finding no solitons to form in 2D is because plasma effects perpendicular to the magnetic field suppress the soliton formation.

The absence of the solitons generated by the analyzed relativistic streaming instability in a more realistic dimensionality (2D vs. 1D) in pulsar star magnetospheres implies that the coherent radiation of pulsars, magnetars, and FRBs is not generated by this kind of soliton (which are known to exist only as a consequence of the necessary simplifications in 1D simulations) unless a proven alternative mechanism for their creation can be established. Whether other plasma modes can radiate via the coherent curvature emission mechanism or whether the streaming instability can radiate via other coherent radiation mechanisms, however, remain open questions. Detailed coherent radiation properties of the relativistic streaming instability in 2D will be studied in a follow-up paper.

Movie

Movie 1 associated with Fig. 5. Access here

Acknowledgments

We acknowledge George I. Melikidze for discussions of analytical solutions of nonlinear plasma wave equations. We also acknowledge the support by the German Science Foundation (DFG) projects MU 4255, BU 777-16-1, BU 777-17-1, and BE 7886/2-1. We acknowledge the developers of the ACRONYM code (Verein zur Förderung kinetischer Plasmasimulationen e.V.). The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gauss-centre.eu/) for partially funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (www.lrz.de), projects pn73ne.

References

All Tables

All Figures

Fig. 1. Comparison of the energy evolution of electric field components parallel, Ex, and perpendicular, Ey, to the magnetic field between 1D (black lines) and 2D (blue lines) simulations. All energies are normalized to the initial kinetic energy. The parallel components of the energy ratio Wx/Wk are always smaller in 2D than in 1D and decreasing. The 2D simulations grow later because they contain a lower initial noise level than the 1D ones at the simulation start. (a) S1–1D + S1–2D. (b) S2–1D + S2–2D. (c) S3–1D + S3–2D.

In the text

Fig. 2. Dispersion diagrams of the electrostatic waves along the magnetic field around the saturation time. Comparison between 1D (top row) and 2D (bottom row) simulations with the same plasma parameters. The solitons dispersion branches are denoted by black arrows and the superluminal L-mode waves by red arrows. The superluminal waves appear in 1D and 2D, while the soliton branch appears only in 1D. (a) S1–1D : ωpt = 1500 − 2500. (b) S2–1D : ωpt = 500 − 1500. (c) S3–1D : ωpt = 4000 − 5000. (d) S1–2D : ωpt = 2000 − 3000. (e) S2–2D : ωpt = 1000 − 2000. (f) S3–2D : ωpt = 3000 − 4000.

In the text

Fig. 3. Same as Fig. 2 but at the end of simulations in a time interval ωpt = 9000 − 10 000. Solitons do not appear in simulations (b) and (d)–(f). While the most intensive waves in 1D are the solitons and superluminal L-mode waves close to k∥ ≈ 0, in 2D, the most intensive waves are close to the light line, and solitons are not present. (a) S1–1D. (b) S2–1D. (c) S3–1D. (d) S1–2D. (e) S2–2D. (f) S3–2D.

In the text

	Fig. 4. Evolution of the electrostatic energy density in simulations normalized to the initial kinetic energy density. Although in S1–1D and S3–1D, the solitons survive until the simulations end, in S2–1D, the solitons disappear before ωpt = 2000. No solitons appear in 2D simulations. (a) S1–1D. (b) S2–1D. (c) S3–1D. (d) S1–2D. (e) S2–2D. (f) S3–2D.
In the text

Fig. 5. Electrostatic energy density of the 2D simulations normalized to the initial kinetic energy density in the x − y plane during the saturation times (top row) and at the simulation ends (bottom row). The evolution is available as an online movie. Only a simulation sub-domain is selected to show the structures formed by the waves. (a) S1–1D : ωpt = 2500. (b) S2–1D : ωpt = 900. (c) S3–1D : ωpt = 3000. (d) S1–2D : ωpt = 10 000. (e) S2–2D : ωpt = 10 000. (f) S3–2D : ωpt = 10 000.

In the text

Fig. 6. Energy density of longitudinal electrostatic waves of 2D simulations in the wave vector space during the saturation time intervals (top row) and at the simulation ends (bottom row). All analyzed intervals are ωpta = 1000 long. The comparison between rows shows that the longitudinal electrostatic waves close to |k| = 0 decrease in time. (a) S1–2D : ωpt = 2000 − 3000. (b) S2–2D : ωpt = 1000 − 2000. (c) S3–2D : ωpt = 3000 − 4000. (d) S1–2D : ωpt = 9000 − 10 000. (e) S2–2D : ωpt = 9000 − 10 000. (f) S3–2D : ωpt = 9000 − 10 000.

In the text

Fig. 7. Comparison of particle distributions as functions of the particle momenta. All electrons and positrons are used from the simulations. The distributions for 1D and 2D simulations at the simulation start, instability saturation, and the simulation end (top row). Comparison of the distribution difference between the simulation end and saturation time for 1D and 2D simulations (middle row). Comparison of the distribution difference between 1D and 2D simulations for the saturation time and simulation end (bottom row). The vertical dotted gray lines show selected positions of momenta where the distributions change significantly. The dashed gray lines show the initial distributions (top row) and a zero distribution difference (middle and bottom rows). (a) S1–1D + S1–2D. (b) S2–1D + S2–2D. (c) S3–1D + S3–2D.

In the text