Particle-in-cell simulations of electron–positron cyclotron maser forming pulsar radio zebras

Context. The microwave radio dynamic spectra of the Crab pulsar interpulse contain fine structures represented via narrow-band quasiharmonic stripes. This pattern significantly constrains any potential emission mechanism. Similarly to the zebra patterns observed in, for example, type IV solar radio bursts or decameter and kilometer Jupiter radio emission, the double plasma resonance (DPR) e ff ect of the cyclotron maser instability may interpret observations. Aims. We provide insight at kinetic microscales of the zebra structures in pulsar radio emissions originating close to or beyond the light cylinder. Methods. We present the first electromagnetic relativistic particle-in-cell (PIC) simulations of the electron-positron cyclotron maser for cyclotron frequency smaller than the plasma frequency. In four distinct simulation cycles, we focused on the e ff ects of varying plasma parameters on the instability growth rate and saturation energy. The physical parameters were the ratio between the plasma and cyclotron frequency, the density ratio of the ‘hot’ loss-cone to the ‘cold’ background plasma, and the loss-cone characteristic velocity. Results. In contrast to the results obtained from electron–proton plasma simulations (for example, in solar system plasmas), we found that the pulsar electron–positron maser instability does not generate distinguishable X and Z modes. On the contrary, a singular electromagnetic XZ mode is generated in all studied configurations close to or above the plasma frequency. Highest instability growth rates were obtained for the simulations with integer plasma-to-cyclotron frequency ratios. The instability is most e ffi cient for plasma with characteristic loss-cone velocity in the range v th = 0.2–0.3 c . For low density ratios, the highest peak of the XZ mode is at the double frequency of the highest peak of the Bernstein modes, indicating that the radio emission is produced by a coalescence of two Bernstein modes with the same frequency and opposite wave numbers. Our estimate of the radiative flux generated from the simulation is up to ∼ 30mJy from an area of 100km 2 for an observer at 1kpc distance without the inclusion of relativistic beaming e ff ects, which may account for multiple orders of magnitude.


Introduction
Since their first discovery nearly 50 years ago, pulsars remain a widely studied category of astronomical objects.They possess one of the strongest magnetic fields in the known universe, enabling fascinating effects, such as the pair plasma production, which are considerably difficult to reproduce in laboratory.Although the radio emission of pulsars has been subjected to intense research (Sturrock 1971;Ruderman & Sutherland 1975;Usov 1987;Petrova 2009;Philippov et al. 2020;Melrose et al. 2021), there is still no general consensus on the mechanisms behind the origin of pulsar complexly structured radio emission.

Zebra pattern observed in Crab pulsar emission
The dynamic spectra of the Crab pulsar radio emission, obtained via the radio instruments like Karl G. Jansky Very Large Array, Arecibo Observatory, and Robert C. Byrd Green Bank Telescope (Hankins & Eilek 2007;Hankins et al. 2015), contain valuable information about the emission mechanism at plasma kinetic microscales.Remarkably, the dynamic spectra of the Crab pulsar radio interpulse are characterized by fine structures, represented via relatively narrow-band quasiharmonic stripes.The stripes are similar to Zebra patterns observed in Type IV solar radio bursts (Elgarøy 1961;Kuijpers 1975;Zheleznyakov & Zlotnik 1975a) or Jovian decametric to kilometric radio emissions (Litvinenko et al. 2016).Specifically, the dynamic spectra from the Crab pulsar show many similarities to the solar radio bursts (for example Chernov et al. 2005).The spacing between these quasiharmonic stripes is not constant.Instead, the spacing increases with the observed frequency (Zlotnik et al. 2003;Hankins & Eilek 2007).From the observations, summarized in Eilek & Hankins (2016), the relation between the emission frequency spacing between stripes ∆ f obs and the stripe center frequency f obs was found to be ∆ f obs = 6 × 10 −2 f obs . (1) The observed number of stripes varies, with both lower (up to ten, Hankins & Eilek 2007) and higher (more than ten, Hankins & Rankin 2010) amounts of stripes detected.The zebra pattern in the Crab pulsar spectra is detected between 6 and 30 GHz; however, the emission becomes rare at ∼ 10 GHz and higher (Hankins et al. 2016).This discovery significantly constrains any possible emission mechanism.
Beyond the light cylinder distance of the neutron star, the magnetic field strength reduces below ∼ 10 6 G, and the electron-cyclotron maser instability can be relevant.The required magnetic field for the observed frequency range f = 6-10 GHz (Hankins & Eilek 2007) corresponds to relatively low lepton gyro-harmonics and is obtained for weaker magnetic fields B ∼ 10 2 G.The observed frequency f obs is related to the source emission frequency f by the Doppler relation where θ is the angle between the velocity and the direction of the emission in the observer's frame, and β = V/c represents the emission source velocity V related to the speed of light c.
For cos θ ≈ 1 and β ≈ 0.82 (Zheleznyakov 1971), the relation reduces to Assuming the emission is generated near the plasma frequency f p , then the plasma frequency is in the range 1.8 − 3 GHz (Zheleznyakov et al. 2012), and the estimated particle number density at the source is N ≈ 0.4 − 1.1 × 10 11 cm −3 .The requirement of a relatively weak magnetic field and a higher particle number density could be considered as the largest caveat of the proposed model; however, the intensity of the magnetic field can significantly decrease close to the magnetospheric current sheets and at distances of several light cylinder radii where the magnetic field is partially dissipated (Cerutti et al. 2020).Proposed configurations supporting the realization of the cyclotron maser include a neutral current sheet with a transverse magnetic field or a highly elongated magnetic trap.The geometry of the source is discussed in further detail in Zheleznyakov et al. (2012) or Zheleznyakov & Shaposhnikov (2020), showing that the position of the emission bands and their frequency spacing is determined by variations in plasma density and magnetic field intensity along the current sheet, and the emission bands' fine structures depend on the plasma and magnetic field variations orthogonal to the current sheet.

Electron-cyclotron maser and the double plasma resonance
The electron-cyclotron maser is a process capable of producing coherent radiation from plasma.Originally proposed as a somewhat exotic concept, it gained a lot of traction during the last 30 years as a dominant mechanism for producing strong coherent emissions in collisionless magnetized plasmas (Wu & Lee 1979;Melrose & Dulk 1982).The concept of an electron-cyclotron maser (or positron-cyclotron maser) presents a mechanism that emits radiation at frequencies near the electron cyclotron frequency and may emit at its harmonics (Twiss 1958;Schneider 1959).
The maser instability is driven by a positive gradient of velocity perpendicular to the magnetic field in the velocity distribution function (VDF).The instability growth rate in the perpendicular direction to the magnetic field depends on the positive gradient of the VDF f as (Winglee & Dulk 1986) where u ⊥ = p ⊥ /m e is the perpendicular momentum per mass unit.The maser amplification occurs due to a resonant interaction between the plasma waves and energetic, charged particles.Specifically, the resonance condition is where k ∥ is the parallel wave number component of the wave vector k = (k ∥ , k ⊥ ) with respect to the magnetic field B, v ∥ is the parallel component of the particle velocity v = (v ∥ , v ⊥ ), ω c is the particle cyclotron frequency, s is the gyro-harmonic number, and γ = (1 − v 2 c 2 ) − 1 2 is the particle Lorentz factor.Due to the double plasma resonance, the maser instability produces a sharp, exponential increase of the resonant waves in a magnetized plasma.Under the condition when the upper hybrid frequency coincides with the harmonics of cyclotron frequency, the radiation mechanism explains the observed zebra patterns in the radiograms of the Sun, Jupiter, and particularly the Crab pulsar.The analogy between the solar radio zebra patterns and the radiation of the Crab pulsar could be interpreted via the DPR effect, suggesting that the pulsar magnetospheric plasma contains local non-relativistically hot regions with relatively higher particle number density, weaker magnetic fields to satisfy Eq. ( 6), and a presence of unstable loss-cone particle distribution.Owing to the DPR effect, the electrostatic waves are enhanced at frequencies and wave numbers, which fulfill the plasma dispersion relation and the resonance condition (5) in such nonequilibrium plasma.Conversion of the electrostatic waves into electromagnetic waves then leads to strong emission of radiation.According to Pearlstein et al. (1966) and Zheleznyakov & Zlotnik (1975a,b,c), the instability is caused by an unstable particle velocity distribution, steeply increasing at the hybrid frequencies or wave frequencies that fulfill the resonance condition (5), such us Bernstein waves.Dispersion relations of the Bernstein modes (Bernstein 1958) can be obtained by solving equation where I n (λ) is the modified Bessel function with argument λ = k 2 ⊥ v 2 tb /2ω 2 c .In the limit k ⊥ → 0, n = 1, and ω p ≫ ω c , the Bernstein dispersion branch approaches the upper hybrid branch defined as where v tb is the plasma thermal velocity.
Given that the upper hybrid frequency is close to harmonics of the electron cyclotron frequency with s being the gyroharmonic number, we obtain the frequency condition from Eq. 5 for instability growth as Because in the emission model, individual zebra stripes are generated at various ω p /ω c ratios, the resonance condition can be satisfied for several sequential harmonics only in weakly magnetized plasmas where the cyclotron frequency is much lower than the plasma frequency reducing to the approximate equality In this limit, the maser instability supports the growth of the electrostatic upper hybrid waves.However, a coalescence process is required to convert the electrostatic waves into electromagnetic radiation capable of escaping the source region (Winglee & Dulk 1986).
The solutions of plasma electromagnetic dispersion perpendicular to the magnetic field in electron-proton plasma are the "fast" X and O mode waves propagating above the X mode cutoff frequency ω X and ω p respectively, and the "slow" Z mode which is confined below ω UH and thus cannot leave the plasma The X and O modes are polarized perpendicular and parallel, respectively, to the plane of the wave and magnetic field vectors.
Although the Z mode cannot escape the plasma, Li et al. (2021) found a coalescence of the Z mode with the electrostatic mode into the escaping X mode for the solar coronal case.
As the proposed interpulse emission source of pulsars is small in area, it is difficult to measure the properties of the source observationally.Numerical simulations that resolve electromagnetic and relativistic effects on kinetic scales can grasp the system as a whole and remain self-consistent.Given the effects studied in the pulsar magnetosphere under the constraints of the electron-cyclotron maser theory (Zheleznyakov et al. 2016) occur on the kinetic level, the use of the particle-in-cell method (PIC) is advantageous due to its kinetic-scale resolution.The method was successfully used to study the solar radio emission via the electron-cyclotron maser instability and double plasma resonance effects (Benáček & Karlický 2018;Li et al. 2021;Ning et al. 2021).However, the nonlinear evolution, mode conversion, and the resulting electromagnetic waves of the pulsar electron-positron cyclotron maser were not studied to the best of our knowledge yet.

Aims and structure of the paper
In this paper, we investigate the electron-positron plasma for the first time using particle-in-cell simulations that resolve the plasma kinetic microscales at which the radio emission originates.The aims are to investigate how the instability evolves, how it differs from the electron-proton version, and what are instability radiation properties from the perspective of a potential zebra stripe radiation-producing mechanism.
In Section 2, we discuss the used particle-in-cell code, its configurations, and the motivation behind the choice of the simulated plasma parameters.We analyze the results in Section 3, discussing the impacts of varying parameters of the studied plasma (Sect.3.1), the dispersion properties of the system and a comparison between the electron-proton and electron-positron plasma under the same conditions relevant for the studied environment (Sect.3.2), and the evolution of the velocity distribution and its relation to the instability growth rate (Sect.3.3).Section 4 follows, putting our results into a broader physical context, drawing conclusions, and presenting possible future considerations.

Methods
Modified TRISTAN PIC code (Buneman 1993;Matsumoto & Omura 1993) is used for the simulations.The code is threedimensional, fully electromagnetic, and includes effects of special relativity.Particle motion is resolved via the Newton-Lorentz system of equations using the Boris push algorithm (Boris et al. 1970); electromagnetic fields are solved using the Maxwell equations via the Yee lattice (Yee 1966); the current deposition is done using current-conserving scheme by Villasenor & Buneman (1992).The system is evolving self-consistently.TRISTAN uses relative scales for the sake of simplifying the calculations.The computational domain is a rectangular grid with cell dimensions set to ∆x = ∆y = ∆z = ∆ = 1.To study the given system, periodic boundary conditions are used in all three dimensions.Time discretization is ∆t = 1.The size of the time step directly relates to the plasma frequency ω p as ∆t = 0.025 ω −1 p .The value of the speed of light is set to c = 0.5, ensuring the CFL condition holds.
We assume two species of particles, electrons and positrons, which only differ in charge sign, which is either Q = −1 or Q = +1.Initially, both particle species are loaded into the simulation in two distinct populations: a population representing the cold, background thermal plasma with the thermal (background) velocity v tb and a superimposed population of the hot loss-cone plasma component with characteristic loss-cone velocity v th .The background particles are characterized by Maxwell-Juttner distribution.Following Zheleznyakov & Zlotnik (1975a,b) and Winglee & Dulk (1986), hot particles are characterized via the Dory-Guest-Harris (DGH) distribution (Dory et al. 1965) in the form where u ∥ = p ∥ /m e and u ⊥ = p ⊥ /m e are the longitudinal and transverse components of the particle velocity relative to the magnetic field represented in terms of relativistic momentum p = (p ∥ , p ⊥ ) and electron mass m e , and v th is the characteristic loss-cone velocity, at which the loss-cone distribution reaches its maximum.The DGH velocity distribution function proved to be a good approximation for particles trapped in a magnetic field mirror (Dory et al. 1965), characterized by a deficit of particles with small perpendicular velocities and zero mean velocity parallel and perpendicular to the magnetic field.The particles are initially distributed uniformly in space.
The simulations are carried out in four cycles, each focusing on an independent parameter or configuration of multiple variables, and how the parameter values influence the studied instability.The detailed parameter setup of the used simulations is shown in Table 1.
Table 1: Simulation parameters of each simulation cycle.(A) nine simulations covering the frequency ratio range ω p /ω c = 10 -12 equidistantly with increments of ω p /ω c = 0.25; (B) four simulations with gradually increasing characteristic thermal velocity of the hot plasma component v th from 0.15 c to 0.5 c; (C) simulations of four different density ratios between hot loss-cone and background particles for frequency ratio values of ω p /ω c = 3, 5, 11 (twelve simulations in total).The simulation of the electron-proton plasma used the values in column (D).First, the impact of the varying ratio between the plasma and the cyclotron frequency ω p /ω c (Table 1, column A).The simulations were carried out to test whether the instability has the highest growth rate at integer values of the frequency ratio, and how increasing gyro-harmonic number s affects the resulting saturation energy of the formed plasma waves.We studied the range of ω p /ω c = 10-12, which is similar to the environment found in solar plasma (Yasnov et al. 2017;Li et al. 2021).The choice of the frequency ratio range is based on the assumption of reduced magnetic field intensity in the studied region in a light-cylinder distance due to magnetic reconnection and local magnetic traps (Zheleznyakov et al. 2012;Zheleznyakov & Shaposhnikov 2020), although it is worth noting that the choice is somewhat arbitrary.
Second, the behavior of the instability in configurations with a varying loss-cone characteristic velocity v th (Table 1, column B).This way, we investigate the behavior of the more energetic hot population of the plasma by surveying the range v th = 0.15c-0.5cand evaluate the influence of the characteristic velocity on the studied instability by assessing its growth rate and saturation value.Thermal velocity used for the background is fixed to v tb = 0.03c, corresponding to particles with a temperature of T ≈ 5.3 MK.Increasing the background plasma temperature has very little effect on the growth rate of the instability (Benáček et al. 2017;Benáček & Karlický 2018); therefore, we do not study this effect.
Third, the impacts of the varying ratio between the hot and background particle number density n th /n tb (Table 1, column C).For the condition n th ≪ n tb , the dispersion properties are governed by the background plasma population, while the instability is governed by the hot, nonequilibrium plasma component.However, by selecting different values of the number density ratio between the populations, we address the possibility of the hot component also influencing the dispersion properties.In our simulations, the hot component is increasingly more prevalent, from a ratio n th /n tb = 1/10 up to a maximum of n th /n tb = 1/1 (total number of macroparticles remains the same).The former constraints the plasma dominated by its background component, still keeping the growth rate large enough to evolve the instability in a few thousand plasma periods.The latter can occur through magnetic reconnection as was shown, e.g., by Yao et al. (2022).We analyze the behavior of increasing density for frequency ratios ω p /ω c = 11, 5, 3. Similar values were also estimated in the model of the Crab zebra (Zheleznyakov & Shaposhnikov 2020).The intensity of the magnetic field in the simulation under the assumption of f pe = 1.8 − 3 GHz (Zheleznyakov et al. 2012) is 58-97 G for ω p /ω c = 11; 129-214 G for ω p /ω c = 5; and 214-357 G for ω p /ω c = 3.
Fourth, one referential simulation of the electron-proton plasma (Table 1, column D).Configuration of the simulation is the same as the simulation generating the strongest XZ mode from the previous cycle to compare the electron-positron and electron-proton systems, their wave regimes, and how different masses influence the dispersion properties of the studied system.The thermal velocity of the background particles listed in Table 1 applies for electrons and is calculated accordingly for protons.
The choice of the studied parameters is to assess the behavior of the electron-positron cyclotron maser instability in various scenarios of physical constraints required to generate the sought radio emission.The size of the simulations is the same for all the studied cases.The used grid dimensions are 6144 ∆, 12 ∆, and 8 ∆ in the x, y, and z directions respectively.The magnetic field B = (0, 0, B) has a nonzero component in the z direction, as we aim to study the effects occurring in the direction orthogonal to the magnetic field (Zheleznyakov & Shaposhnikov 2020).Each grid cell is initiated with 1100 macroparticles at the start of the computation (summing up to the total of 3.2 × 10 8 macroparticles), sufficient to properly distinguish physical effects from the particle noise.All simulations run for 60 000 time steps, equal to ω p t = 1500.This time is sufficient for the instability to reach saturation for all of the studied scenarios, which either enter or gradually decrease towards relaxation.
From the simulations, the temporal evolution of the electrostatic energy U Ex in the dominant direction and VDF f (v ⊥ , v ∥ ) were computed.For the plasma wave analysis, dispersion diagrams are obtained through the time-space Fourier transform of the electric field components.The obtained frequency resolution is ∆ω/ω p = 0.004 and the wave number resolution is ∆k ⊥ ω p /c = 0.02.

Effects of varying frequency ratio, characteristic velocity, and number density ratio
To test the growth conditions of the electron-positron cyclotron maser with the DPR effect, we computed nine simulations (see Table 1, column (A) with varying plasma-to-cyclotron frequency ratio ω p /ω c from 10 to 12 with increments of 0.25 ω p /ω c for each simulation.We evaluated the instability growth rate Γ/ω p , taken as the exponential coefficient of the time-dependent electrostatic energy U Ex (t) ∝ e 2Γt computed from the x component of the electric field E in each grid point and integrated over the whole domain.The dependency of the growth rate Γ/ω p on the frequency ratio ω p /ω c with a cubic interpolation of data in the studied interval is shown in Fig. 1.Owing to the DPR effect (11), the maxima are located at the integer values of the frequency ratio, and minima are in-between.
Taking a closer look at the configurations with ω p /ω c from 10.5 to 11.5 (interval between two subsequent minima), the temporal evolution of the electrostatic energy scaled to the initial kinetic energy U Ex /E k0 is shown in Fig. 2a.Along with the configuration with the integer value of ω p /ω c having the highest growth rate, it also exerts the highest saturation energy.
To further build on the foundation of the configuration with the highest saturation energy (ω p /ω c = 11, Table 1, column B), five subsequent simulations with increasing loss-cone characteristic velocity v th were computed to investigate whether increased v th has a positive impact on the instability growth rate and to what extent the saturation value can be increased.The temporal evolution of electrostatic energy scaled to the corresponding initial kinetic energy U Ex /E k0 is in Fig. 2b.The configuration with v th = 0.2c shows the highest growth rate Γ relative to the total kinetic energy of the simulation; however, the highest energy ratio U Ex /E k0 is reached in configuration with v th = 0.3c.Therefore, the highest growth rate is not necessarily accompanied by the highest saturation energy.What is apparent is that the instability grows strongly with v th up to 0.3c.Further increasing the characteristic velocity does not increase both the instability growth rate and saturation value.The decrease of the growth rate at sufficiently high characteristic velocities implies that the efficiency of the instability is low for highly relativistic energies (Verdon & Melrose 2011).
Along with the study of increasing characteristic velocity v th , the investigation of an increasing ratio between the hot and background plasma component n th /n tb for four different cases was conducted (Table 1, column C).Fig. 2c shows the temporal evolution of electrostatic energy scaled to the corresponding initial kinetic energy U Ex /E k0 for density ratios n th /n tb = 1/10, 1/3, 1/2, and 1/1.We only consider the x component (with other dimension's components averaged), because we primarily focus on effects perpendicular to the magnetic field.Higher density ratios represent the possibility of a more energetic plasma environment, resulting in the instability forming faster and with a higher growth rate Γ; however, the saturation energy is comparable between all studied cases., where E k0 is the total initial kinetic energy, for different values of the ratio between plasma and cyclotron frequency ω p /ω c .b) Evolution of U Ex /E k0 , for increasing value of the characteristic thermal velocity of the hot plasma component v th .c) Evolution of U Ex /E k0 for increasing value of the ratio between the hot and background plasma density n th /n tb .

Wave analysis and dispersion of the plasma
The wave dispersion is obtained by analyzing the electric field E in the Fourier space.First and foremost, we compare an electron-proton plasma simulation and an electron-positron one (Table 1, column D).The configuration of these referential simulations is identical except for the mass ratio between positive and negative charge particles, where we use m + /m − = 1836 (mass ratio between the proton and electron) for the electronproton plasma and m + /m − = 1 to represent the electron-positron plasma.The only parameters worth noting are the frequency ratio ω p /ω c = 5 and the number density ratio n th /n tb = 1/2.Even though the electron-positron and the electron-proton plasma have different total plasma frequencies, we assume that the contribution of the ion plasma frequency to the total plasma fre- quency is negligible where m e and m i are the electron and ion rest mass, respectively.We scale the presented quantities for both types of plasma to the total plasma frequency, as the emission occurs close to integer multiples of ω p .In the case of the electron-positron plasma, scaling to ω p is justified since we do not see the separate effects of either the electron's or positron's plasma frequency in the dispersion.The comparison between the dispersion of the electron-proton and the electron-positron plasma is shown in Fig. 3. Here, the dispersion in x, y, and z components of the electric field E(k ⊥ , ω) is obtained through the Fourier transform of E(x, t).In the case of the electron-proton plasma, O, X, Z, and Bernstein wave modes are generated, yielding results similar to those obtained in solar plasma simulations (Benáček & Karlický 2019;Ni et al. 2020;Li et al. 2021).
-20 -10 0 10 20 0 For our case, the most important is the dispersion in E y , as it is the transverse field component in which the X and Z modes are generated.The largest difference to the electron-proton plasma is that the plasma does not generate X and Z modes separated in frequency for k ⊥ = 0. Instead, an XZ mode is generated.This mode is generated directly at the plasma frequency ω p as opposed to the X mode being generated above and Z mode below ω p for k ⊥ = 0 in the electron-proton plasma.Another noted difference is the presence of a whistler mode (W) in the dispersion of E y of the electron-positron plasma.Apart from these remarks, further differences between both plasma types are essentially indistinguishable, with both plasmas generating Bernstein waves in E x and O mode waves in E z .The growth rate of the instability is nearly identical between both cases.
Based on the first two simulation cycles, we calculated 12 simulations across a range of different density ratios n th /n tb and plasma to cyclotron frequency ratios ω p /ω c listed in Table 1(C).The dispersion in E x and E y components of E(x) (perpendicular to the magnetic field) are shown in Figs.4-6.Each figure shows a dispersion with a fixed frequency ratio ω p /ω c = 11 (Fig. 4), 5 (Fig. 5), and 3 (Fig. 6) for four different configurations with de- creasing density ratio n th /n tb .The dispersion diagrams are overlaid with frequency ω integrated over wave number k ⊥ , noted as frequency profile F(ω) of the electrostatic Bernstein waves in the longitudinal E x component and the electromagnetic XZ mode in the transverse E y component.For the integration of the wave modes, a of these modes from the background is required.The integration was done accurately in the case of the XZ mode, where the curve could be fitted with a hyperbole analytically and a thin region around the curve of the XZ mode with frequency half-width 0.2 ω p was selected; however, in the case of Bernstein waves there is no straightforward analytic solution for the electron-positron plasma with arbitrary density ra-Article number, page 7 of 12 tio; therefore, a rectangular window encompassing the maxima of the Bernstein wave intensity in the range k ⊥ c/ω p = -15 to 15, ω = (0 -1.5)ω p was used, resulting in increased noise inclusion compared to the integration of the XZ waves.The integrated profiles F(ω) are normalized to the maximum of the top-most case F max , therefore F(ω)/F max is equal to 1 for the top-most case and can be higher or lower for the rest.
Lastly, we present a calculation of the estimated mode energy density for the Bernstein modes and the XZ mode.The estimated energy is obtained through the inverse Fourier transform of the mentioned modes (Bernstein modes in E x and XZ mode in E y ).The temporal evolution of the mode energy scaled to the initial kinetic energy of the simulations Table 1(C) is shown in Fig. 7. Most of the total electrostatic energy is carried by the Bernstein waves (noted as E B ), as can be seen by comparing Fig. 7 (bottom-left) and Fig. 2b.There is only about a 10% difference between the maxima of the integrated Bernstein waves and the electrostatic energy U Ex of the simulation.In other words, Bernstein waves account for 90% of the U Ex .On the contrary, the electromagnetic energy carried by the XZ modes (noted as E XZ ) is reaching maxima that are four orders of magnitude lower than those of the Bernstein waves.However, it should be noted that the integration of the Bernstein waves was less accurate than the analytic solution of the XZ mode because the integration includes more field noise.

Velocity distribution function and its evolution
To comprehend and constrain the effects observed in the investigation of how the increasing characteristic thermal velocity v th of the hot plasma component impacts the formation of the instability, we consider the VDF temporal evolution.The existence of a hot population described via the DGH distribution results in a positive gradient in the distribution of the perpendicular velocity u ⊥ , producing unstable plasma via the DPR effect.The VDF f (v ⊥ , v ∥ ) (a function of perpendicular and parallel velocity component) obtained from simulations Table 1(B) (also Fig. 2b) are shown in Fig. 8. Distributions are computed at the start of the simulation (ω p t = 0) and after 1200 plasma periods (ω p t = 1200).Combining the results from Fig. 2b and Fig. 8, a higher growth rate is accompanied by larger shifts in the velocity distribution.For plasma with lower characteristic velocity (v th ≤ 0.3c), the instability has a higher growth due to its positive gradient in f being located at lower u ⊥ .
The evolution of the VDF is further visualized in Fig. 9, showing the difference between the initial and evolved VDF , overlaid with resonance curves corresponding to the solution of Eq. ( 11) for gyro-harmonic numbers s = 12, 13, 14, and 15.The scenarios with lower values of v th (up to 0.3c), for which the instability grows exponentially, also show a substantial change between the final and initial state of f (v ⊥ , v ∥ ).There is a decrease of particles delimited by the resonance curve of s = 12 towards zero perpendicular velocity; however, the particles are shifted to higher velocities periodically between the resonance curves of higher (s > 12) gyro-harmonic numbers.

Discussion and conclusions
The distinct zebra pattern in the Crab pulsar interpulse emission (Hankins & Eilek 2007) could be explained by electron-positron cyclotron maser emission driven via the double plasma resonance effect (Zheleznyakov et al. 2016).The cyclotron maser requires a relatively weak magnetic field that can be reached close to current sheets of the magnetic reconnection or in the pulsar wind beyond the light cylinder.Using the modified 3D relativistic PIC code TRISTAN, we computed a few series of simulations, surveying impacts of the plasma to cyclotron frequency ratio in the range ω p /ω c = 10-12, characteristic velocity of the hot loss-cone plasma v th = 0.15c -0.5c, and the number density ratio between the background and hot plasma components n th /n tb = 1/1 -1/10 for frequency ratios ω p /ω c = 3, 5, and 11.
Comparing the simulation of the electron-proton plasma with the electron-positron plasma (Fig. 3) for similar plasma parameters, distinguishable X (cutoff above ω p ) and Z (cutoff below ω p ) modes are generated in the electron-proton plasma.In the case of the electron-positron plasma, only one XZ mode is generated.This XZ mode approaches ω p for k ⊥ = 0 with a slight dependence of the cutoff frequency on the density ratio n th /n tb between the hot and background plasma components.Similar dependency of the decreasing cutoff frequency wave mode with increasing electron temperature is also found in the analytic study by Rafat et al. (2019) for relativistically hot plasmas.Although the X mode moves to higher frequencies and Z mode to lower frequencies depending on the gyrofrequency in electron-proton plasmas (Li et al. 2021), the XZ mode is generated exclusively near the plasma frequency of the electronpositron plasma.We found that the growth rate of the instability is the highest at integer values of plasma to cyclotron frequency ratio.Our further evaluation of the instability shows that the highest saturation energy, normalized to initial kinetic energy, is reached for losscone characteristic velocities v th in the range 0.2-0.3c.The highest saturation energy was obtained in simulation with v th = 0.3c while the highest growth rate was obtained for v th = 0.2c.This discrepancy suggests that the highest saturation the electrostatic energy is not necessarily accompanied by the highest instability growth.Furthermore, increasing the characteristic velocity of the loss-cone component diminishes the instability growth, showing the mechanism is less efficient for relativistically hot plasma where higher energy densities are obtained than for lower characteristic velocities.
We observe a correlation between the growth rate of the instability and the evolution of the VDF f (v ⊥ , v ∥ ).The highest growth rates are accompanied by the most changing distributions over the course of the simulation.While some particles shift to lower perpendicualr velocities, others shift to higher velocities.These changes in the VDF are encompassed by DPR resonance curves.
Even though the instability growth rate peaks are found at integer values of the frequency ratio ω p /ω c , this does not necessarily imply that the strongest intensity of the electrostatic and electromagnetic waves is located at frequencies of the cyclotron harmonics.The intensity maxima of the electromagnetic waves, determined in the dispersion of the electric field components (Fig. 4, 5, and 6), are located mainly at frequencies between the gyro-harmonics.Furthermore, these maxima of wave intensity in the ω-k ⊥ domain are contained within a frequency range ω = 1-2 ω p .
The particle number density ratio shifts in the phase space (see Fig. 9) with the resonance curves of the gyro-harmonics slightly above the plasma frequency, and their intensity gradually decreases with increasing the gyro-harmonic number s, effectively constraining the emission within the resulting frequency range.Thus, the electron-positron maser instability driven via the DPR effect could be the mechanism responsible for the generation of narrowband emission observed in the dynamic spectra of the Crab pulsar (Hankins & Eilek 2007).

Constraints on radio emission region
With the currently considered model of the zebra emission, we can constrain the density ratio in the emission region.Because the frequency distances between observed zebra stripes are not constant, the individual stripes are generated in various plasma regions as the ratio between the plasma density and cyclotron frequency changes.That implies that each emission region emits at only one specific frequency, and the observed zebra pattern is produced as a compilation of several emission regions.Therefore, the maser mechanism also has to emit at one specific frequency.Assuming the XZ mode emission, the density ratio n th /n tb ≲ 1/10 is required to fulfill the condition.On the other hand, if a zebra pattern with an equidistant frequency distance between stripes is observed, the emission can be produced by one plasma region with a high density ratio n th /n tb ≳ 1/1.Moreover, the density ratio increases with the number of detected stripes.
The emission generation is a two step process.First, the kinetic instability excites arbitrary electrostatic wave modes.Then, the electromagnetic waves are generated through a coalescence of the electrostatic modes.In the solar plasma, a coalescence of Z mode and Bernstein (UH) mode accounts for the cyclotron maser emission (Ni et al. 2020).In the electron-positron plasma, a coalescence of two Bernstein modes (Kuznetsov 2005) into the XZ mode may account for the observed quasi-harmonic emission of the pulsar radio interpulse.This can be seen in Fig. 5 b, where the highest peak of the XZ mode is at the double frequency of the highest peak of the Bernstein modes.It indicates that the radio emission is produced by a coalescence of two Bernstein modes with the same frequency and opposite wave numbers.The coalescence process, however, is ineffective, as the XZ mode reaches energies approximately four orders of magnitude lower than those of Bernstein modes (Fig. 7).Moreover, since our PIC simulation is grid based, any coalescence based process is inherently less effective due to a limited set of wave vectors k included in the simulation.Because the emission region generates the emission at one frequency, other zebra stripes need to be generated at different locations, which explains the non-constant spacing between zebra stripes.

Estimation of flux density
To present an estimate of the generated flux density of the electromagnetic emission, we base our calculation on the simulation with the strongest XZ mode (ω p /ω c = 5, n th /n tb = 1/2).We estimate the radiative flux under the following assumptions: 1) scaling the energy of the XZ mode to the total kinetic energy of the simulation, 2) an area of the source as 100 km 2 , 3) the distance to the source is 1 kpc (distance of the Crab pulsar is approximately 2 kpc), and 4) neglecting the radiative transfer effects.The kinetic energy of a particle is given by (15) Considering the upper limit of the particle number density (Zheleznyakov et al. 2016) n = 1.2 × 10 11 cm −3 , the estimated total kinetic energy density is The highest energy of the XZ mode obtained from the simulations is ≈ 1.5 × 10 −7 E k .Under the above-mentioned assumptions, we a radiative flux P ∼ 1 mJy.Taking into account the emission source velocity (β = 0.82), the received flux is P ∼ 30 mJy.This value is considerably lower than the observed radiative flux density (Hankins & Eilek 2007) by approximately 4 to 5 orders of magnitude; however, one should bear in mind that the domain size of the presented simulations was small (≈ 350 cm 3 ) and effectively 1D.On the other hand, 2D or even 3D simulations could simulate the effects of an arbitrary propagation angle instead of limiting us to the perpendicular propagation to the magnetic field.Furthermore, due to the grid-based nature of our simulation, any coalescence process is less effective.Relatively low particle number density can be responsible for the widening of the frequency peaks due to the effects of the numerical noise.The Doppler effect could increase the power of the emitted waves if the emission source moves toward the observer with a higher velocity than assumed above.The Lorentz factor for such a Doppler boost can be obtained in the order of γ ∼ 10 2 − 10 4 in the pulsar wind (Lin et al. 2023).That could lead to an increase of the emitted power by factor of ∝ γ 3 for γ ≫ 1 (Rybicki & Lightman 2004, Chapter 4.8).This way, for γ ∼ 10 3 we can obtain the emission power of ∼10 6 Jy.However, we must note that the Doppler effect does not only increase the power, but it also increases the wave frequency.Hence, for the same observed frequency, the Doppler frequency shift would require wave emission at frequencies much lower than the plasma frequency, much lower plasma frequencies, or kinetic and 'XZ' mode energy densities higher than estimated for our case.
We conclude that the detected emission power from our PIC simulations is too low to explain the observed intensities of the zebra pattern.To obtain the observed intensities, the area of the emission region or the conversion efficiency from particles to the electromagnetic waves is required to be ∼ 5 orders of magnitude higher.Therefore, future considerations of the electron-cyclotron maser instability should include twodimensional analysis of the (k ⊥ , k ∥ ) space to determine whether the radiation is exclusively generated in a wide or narrow angles, how the radiative power depends on the emission angle, and what the contribution by the relativistic beaming effect is.Furthermore, larger domain size would allow much less limited set of k vectors, which might increase the efficiency of the coalescence process.

Fig. 1 :
Fig.1: Dependency of the instability growth rate Γ/ω p on the frequency ratio ω p /ω c .Consistent with the theoretical conclusions, the maxima are located at integer values and minima between them.

Fig. 2 :
Fig.2: Evolution of the electrostatic energy of the x field component U Ex scaled to the initial kinetic energy E k0 .a) Evolution of U Ex /E k0 , where E k0 is the total initial kinetic energy, for different values of the ratio between plasma and cyclotron frequency ω p /ω c .b) Evolution of U Ex /E k0 , for increasing value of the characteristic thermal velocity of the hot plasma component v th .c) Evolution of U Ex /E k0 for increasing value of the ratio between the hot and background plasma density n th /n tb .

Fig. 3 :
Fig.3: Dispersion diagrams for electron-proton plasma (left column) and electron-positron plasma (right column).Frequency ratio is ω p /ω c = 5, density ratio between the hot and background plasma is n th /n tb = 1/2.Rows show dispersion in the E x (longitudinal waves), E y (transverse waves), and E z components.Bernstein (B), X, Z, XZ, whistler (W), and ordinary (O) modes are denoted together with plasma frequency ω p (black dotted line).

Fig. 4 :
Fig. 4: Dispersion diagrams for various values of n th /n tb (rows) with frequency ratio ω p /ω c = 11 for the electron-positron plasma.Left column: E x .Right column: E y .Dispersion diagrams are overlaid with an integrated frequency profile scaled to the maximum of the top-most case (n th /n tb = 1).

Fig. 7 :
Fig. 7: Temporal evolution of energy scaled to the initial kinetic energy of the electrostatic Bernstein waves (left column) and the electromagnetic XZ mode (right column) with frequency ratio ω p /ω c = 3, 5, and 11 (rows) for different values of hot/background particle density ratio n th /n tb .

Fig. 8 :
Fig.8: Velocity distribution functions at times ω p t = 0 (left column) and ω p t = 1200 (right column) for increasing value of the loss-cone characteristic velocity.The distributions are normalized to the total number of particles N. Combined with Fig.2b, it is evident that the configurations with higher growth also exert more sizable changes in the velocity distribution function.

Fig. 9 :
Fig. 9: Differences in the velocity distribution functions ∆ f (v ⊥ , v ∥ ) between the final (ω p t = 1200) and the initial state (ω p t = 0) for increasing value of the characteristic thermal velocity v th (rows).Dashed lines represent resonance curves for gyro-harmonic numbers s = 12, 13, 14, 15 with s = 12 being closest to zero velocity.Arrows denote shifts in the VDF contained within the resonance curves.