<figure><img src="image_1.png" /></figure>

HAL
open science

# Intra-pulse variability induced by plasmoid formation in pulsar magnetospheres

I. Ceyhun Andaç, Benoît Cerutti, Guillaume Dubus, K. Yavuz Ekşi

► To cite this version:

I. Ceyhun Andaç, Benoît Cerutti, Guillaume Dubus, K. Yavuz Ekşi. Intra-pulse variability induced by plasmoid formation in pulsar magnetospheres. Astronomy and Astrophysics - A&A, 2022, 661, pp.A130. 10.1051/0004-6361/202243152. hal-03672732

HAL Id: hal-03672732

https://hal.science/hal-03672732

Submitted on 10 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Intra-pulse variability induced by plasmoid formation in pulsar magnetospheres

I. C. Andaç<fnref n="1" />,<fnref n="2" />✉, B. Cerutti<fnref n="1" />✉, G. Dubus<fnref n="1" />, and K. Y. Ekşi<fnref n="2" />✉

Received 19 January 2022 / Accepted 21 April 2022

## ABSTRACT

*Context.* Pulsars show irregularities in their pulsed radio emission that originate from propagation effects and the intrinsic activity of the source.

*Aims.* In this work, we investigate the role played by magnetic reconnection and the formation of plasmons in the pulsar wind current sheet as a possible source of intrinsic pulse-to-pulse variability in the incoherent, high-energy emission pattern.

*Methods.* We used a two-dimensional particle-in-cell simulation of an orthogonal pulsar magnetosphere restricted to the plane perpendicular to the star spin axis. We evolved the solution for several tens of pulsar periods to gather a statistically significant sample of synthetic pulse profiles.

*Results.* The formation of plasmons leads to strong pulse-to-pulse variability in the form of multiple short, bright subpulses, which appear only on the leading edge of each main pulse. These secondary peaks of emission are dominated by the dozen plasmons that can grow up to macroscopic scales. They emerge from the high end of the hierarchical merging process occurring along the wind current layer. The flux of the subpulses is correlated with their width in phase. Although the full-scale separation is not realistic, we argue that the simulation correctly captures the demographics and the properties of the largest plasmons, and therefore of the brightest subpulses.

*Conclusions.* The prediction of subpulses at specific pulse phases provides a new observational test of the magnetic reconnection scenario as the origin of the pulsed incoherent emission. High-time-resolution observations of the Crab pulsar in the optical range may be the most promising source to target for this purpose.

**Key words.** pulsars: general – acceleration of particles – magnetic reconnection – radiation mechanisms: non-thermal – methods: numerical

## 1. Introduction

Rotation-powered pulsars are observed in a wide range of the electromagnetic spectrum from radio waves up to gamma-rays (Harding 2022). They exhibit significant fluctuations in their pulse profiles as well as substantial variations in flux (Rickett 1969; Sieber 1982). Intra-pulse variability has been reported in the radio band only (Rankin et al. 1970; Lyne et al. 1971) because of a lack of sensitivity at other wavelengths. The pulse-to-pulse variability in flux is due to both propagation effects and the intrinsic activity of the source of emission itself. While propagation effects are well known and characterized (Gochelashvili & Shishov 1975; Rickett et al. 1984; Blandford & Narayan 1985; Goodman et al. 1987; Narayan 1992; Kumamoto et al. 2021), the origin of the magnetospheric flux variability is still poorly understood. However, it may be tightly connected to the pulsar emission mechanism (Melrose 1995), and therefore the pulse-to-pulse variability carries precious information on how pulsars work and understanding its origin is of prime importance.

Pulsars are well-modeled as rotating, conducting spheres surrounded by a magnetized plasma (Goldreich & Julian 1969). An electron-positron plasma is created near the neutron star surface and along the open field lines by magnetic conversion (Sturrock 1971), and possibly in the outer magnetospheric regions by photon-photon annihilation (Lyubarskii 1996).

A large-scale current sheet forms where both magnetic polarities meet beyond the light-cylinder radius (Coroniti 1990; Bogovalov 1999), $R_{LC}$, where the corotation velocity is equal to the speed of light. A fraction of the pulsar spin-down is dissipated in the current layer via magnetic reconnection (Philippov & Spitkovsky 2014; Philippov et al. 2015; Chen & Beloborodov 2014; Belyaev 2015; Cerutti et al. 2015, 2020; Cerutti & Philippov 2017; Hu & Beloborodov 2021). Being long (light-cylinder scale) and thin (plasma skin-depth scale), this layer rapidly breaks up into a chain of plasmons (Zelenyi & Krasnoselskikh 1979; Uzdensky et al. 2010). This fragmentation leads to fast magnetic reconnection and efficient particle acceleration (Kagan et al. 2015). This may lead in return to intense nonthermal radiation and explain the high-energy pulsed emission (Lyubarskii 1996; Pétri 2012; Arka & Dubus 2013; Cerutti et al. 2016; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018).

In this work, we investigate whether plasmoid formation and their evolution in pulsar magnetospheres could be a possible source of intrinsic noise in the emission pattern. The intermittent nature of magnetic reconnection as well as the expected broad size distribution of plasmons (Uzdensky et al. 2010; Sironi et al. 2016) make this process a good candidate for a significant source of pulse-to-pulse variability in the observed incoherent magnetospheric regions by photon-photon annihilation (Lyubarskii 1996).

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
e-mail: benoit.cerutti@univ-grenoble-alpes.fr</footnote>

<footnote n="2">² Istanbul Technical University, Faculty of Science and Letters, Physics Engineering Department, 34469 Istanbul, Turkey</footnote>
emission. We perform a two-dimensional (2D) particle-in-cell (PIC) simulation of an orthogonal pulsar magnetosphere restricted to the equatorial plane as first introduced in Cerutti & Philippov (2017). The simulation is evolved over many pulsar rotation periods and with a refined time-resolution to well characterize the statistical properties of plasmoids and their imprint on the light-curve noise.

The paper is structured as follows: the numerical methods and the simulation setup are introduced in Sect. 2. Our detailed analysis of the simulation with an emphasis on the formation of plasmoids and their demographics is described in Sect. 3, and the pulse-to-pulse variability and its connection to plasmoid evolution is presented in Sect. 4. In Sect. 5, we discuss the implications of this work and propose future observational tests to validate the magnetic reconnection scenario as the origin of the incoherent pulsar emission.

## 2. The simulation setup

The simulation is performed with the Zeltron PIC code (Cerutti et al. 2013; Cerutti & Werner 2019). We use nearly the same numerical setup as the one presented in Cerutti & Philippov (2017), which we briefly outline below for completeness; the main differences being a smaller radial extent of the box size but a longer integration time.

The numerical box is initialized with a split magnetic monopole configuration (Michel 1973; Bogovalov 1999). The magnetic moment, $\mu$, is perpendicular to the neutron star spin axis, $\mu \cdot \Omega = 0$, where $\Omega$ is the pulsar angular velocity vector. We use spherical coordinates; however, the simulation is restricted in the equatorial plane ($r, \theta = \pi/2, \phi$). The grid spacing is logarithmically stretched in the radial direction and uniform in the azimuthal direction. The domain inner boundary is set at the stellar surface, $r_{\min} = r_{\star}$. The light-cylinder radius is set to $R_{LC} = c/\Omega = 3r_{\star}$, which corresponds to a millisecond pulsar. The outer boundary is located at $r_{\max} = 30r_{\star} = 10R_{LC}$. A damping layer beginning at $r = 9R_{LC}$ absorbs outgoing particles and electromagnetic waves to mimic an open boundary.

The simulation begins in a vacuum state around the neutron star. Field lines at the stellar surface are set in solid rotation by enforcing the co-rotation ideal electric field, $E = -(\Omega \times r) \times B/c$. The plasma particles – electrons and positrons in our case – are injected from the surface at the $E \times B$ drift velocity. An equal number of electrons and positrons is injected with a number density $n_0 = 10n_{GJ}$, where $n_{GJ} = \Omega B_*/2\pi\epsilon c$ is the fiducial critical density (Goldreich & Julian 1969), and $B_*$ is the surface magnetic field at the poles. Pair production is neglected in the simulation. Its effects are mimicked by the injection of plasma at the surface, which is sufficient to establish a filled, quasi-force-free magnetosphere. The fiducial magnetization parameter on the stellar surface is $\sigma_* = B_*^2/4\pi n_* m_e c^2 = 250$ where $m_e$ is the electron rest mass. The radiation-reaction force is activated in the simulation to capture curvature and synchrotron cooling.

The total number of cells is 2048 $\times$ 2048. Two particles are injected per time-step in every cell located at the surface of the star. The simulation time-step is set to half the Courant–Friedrich–Lewy critical step. There are about $2.58 \times 10^4$ time-steps per spin period, $T_{spin}$. To characterize the intrinsic noise due to plasmoid formation in the wind current sheet, the simulation is evolved up to $t_{\max} \approx 70T_{spin}$.

## 3. Formation and evolution of plasmoids

After a few pulsar spin periods, the simulation box is filled with plasma and a nearly steady-state force-free magnetosphere is established. The wind current sheet appears as two thin Archimedean spirals separated by $\Delta\phi = \pi$, such that two consecutive arms are equally spaced along the radial direction by $\pi R_{LC}$. Figure 1 shows a representative snapshot of the plasma density at the late stages of the simulation ($t = 60T_{spin}$). The striking feature outlined in Cerutti & Philippov (2017) is the fragmentation of the current layers into multiple plasma overdensities trapped in magnetic islands, or plasmoids in the following.

Nearly all of the plasmoids are created in the vicinity of the light cylinder where the wind current sheet forms and fragments under the action of the tearing instability, mediating fast magnetic reconnection (Uzdensky et al. 2010; Loureiro et al. 2012), which in return leads to efficient nonthermal particle acceleration. After their formation, plasmoids dynamically grow by merging with each other in a hierarchical manner along the spiral. Only a few macroscopic structures remain along each spiral at the end of the merging process. For example, in Fig. 1, a giant plasmoid of size on the order of the neutron star itself (full width $w \sim r_{\star}$) lies in the upper spiral at $r \approx 4R_{LC}$. Given that plasmoids move away from the star at a significant fraction of the speed of light, the chain of plasmoids must reform continuously at the base of the layer near the light cylinder. Therefore, in steady-state, a given radius statistically probes a different phase of the merging process (see the zoomed-in view on the current layer, the bottom panel in Fig. 1). This process is frozen by relativistic expansion of the wind beyond $r \ge 5R_{LC}$ (Cerutti & Philippov 2017; Cerutti & Giacinti 2021).

In order to better quantify plasmoid formation and evolution, we performed a systematic measurement of the width of all plasmoids being produced in the simulation. Figure 2 shows the size distribution of plasmoids, characterized by their full width $w$ as a function of radius within a constant $1R_{LC}$ bin size outside the light cylinder. Plasmoids are identified along each spiral arm as local maxima of the plasma density. The plasmoid width is then inferred from the full width at half maximum of the density profile with respect to the center of each plasmoid. Their width is measured in the simulation frame along the azimuthal direction to avoid length contraction along the radial direction due to the relativistic motion of the pulsar wind, which becomes significant beyond a few light-cylinder radii (see, e.g., the oblate shape of plasmoids at large radii in the top panel of Fig. 1). The plasmoid size distribution is reconstructed from the entire dataset, cumulating a total of more than $10^4$ islands. Island sizes range from about $2 \times 10^{-2}R_{LC}$ up to $0.3R_{LC} \approx r_{\star}$. Close to the light cylinder, the distribution peaks at $w_{\min} \approx 2 \times 10^{-2}R_{LC}$. With a layer thickness of $\delta \approx 3 \times 10^{-3}R_{LC}$ at the light cylinder, it is consistent with the most unstable tearing mode $w_{\min} \sim \lambda_{TI} = 2\pi\sqrt{3}\delta$. The size distribution within $2-3R_{LC}$ is the broadest and shows hints of a $N \propto w^{-1}$ dependence as predicted by Uzdensky et al. (2010) for a hierarchical merging process. The small-scale separation achieved in the simulation does not allow us to probe a sufficiently broad distribution to constrain the power-law index accurately. At larger radii, the distribution becomes narrower (logarithmically speaking) and freezes. The monotonic shift of the distribution with radius (i.e., $w \propto r$) is entirely explained by the expansion of the wind. The number of plasmoids decreases to almost one-quarter of the initial number with respect to the initial distribution, meaning that plasmoids undergo two merger events on average by the time they escape from the simulation box.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Plasma density map ($r-\phi$ plane) normalized to the fiducial injected plasma density, $n_0$. Top panel: black disk is the neutron star and the dashed circle shows the light cylinder located at $R_{LC} = 3r_*$, Radii are normalized to $R_{LC}$. The white solid lines and arrows show the magnetic field lines. Bottom panel: zoomed-in view on one current layer delimited by the red dashed line shown in the top panel. The layer has been unfolded to better observe the chain of plasmoids forming along the spiral.</figcaption></figure>

# 4. Pulse and subpulse variability

In this section, we evaluate the impact of plasmoid formation on the synchrotron emission pattern. Following Cerutti et al. (2016), we compute the synthetic synchrotron pulse profile from the simulation particles and fields as seen by an ensemble of observers located at infinity and along the equatorial plane. We recover the main expected features: two bright, nearly identical and symmetric pulses of light, each of about $\Delta\Phi \sim 0.2$ width in pulsar phase, and separated by 0.5 in phase. This pattern stems from the radiation emitted in the current layer. Taking into account the finite photon time of flight, the region of the magnetosphere that is being probed by an observer at a given phase is an Archimedean spiral. Whenever this spiral overlaps with the current sheet spiral (also close to being perfectly Archimedean), there is a pulse of light, which therefore means there are two

pulses per period. This is the caustic effect in the wind current layer identified in Cerutti et al. (2016). The average pulse profile is shown by the black solid line in Fig. 3 (bottom panel).

The novel feature highlighted here is the formation of bright, short subpulses on top of the main pulse profile. These structures are intermittent and appear only on the leading edge of each pulse. The orange solid line in Fig. 3 shows the light curve seen by an observer during one rotation period. It is composed of about ten easily identifiable bright subpulses per pulse, each of about $\delta\Phi \sim 0.01$ width in pulsar phase (i.e., about 10% of the main pulse leading edge width). A quantitative analysis of the full dataset shows that the average number of major subpulses per pulse is $11 \pm 3$. Major subpulses are defined as those with a flux exceeding the mean plus one standard deviation. They appear in Fig. 3 (top panel) as oblique filaments, meaning that these subpulses are actually regions of enhanced emission that
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Plasmoids size distribution at different radii in units of $R_{LC}$ (color-coded, ranging from the 1–2$R_{LC}$ interval up to 8–9$R_{LC}$ interval). The dashed line is $N$ proportional to the inverse of the island width $w$ predicted by Uzdensky et al. (2010).</figcaption></figure>

are drifting in phase. Their drifting motion is too fast to be detected by a single, fixed observer. For this reason, subpulses appear randomly from one period to another. Figure 4 shows the flux as a function of time, measured at the fixed pulsar phase $\Phi = 0.23$ (phase “1” in Fig. 3) where the flux is most variable. At this phase, the standard deviation is similar to the mean flux. It is also a phase at which subpulses seem to form quasi-periodically on a timescale on the order of the pulsar period (see the top panel in Fig. 3).

The origin of subpulses is intimately related to plasmoid formation in the pulsar wind. Figure 5 maps the emitting regions at three phases selected in the leading edge of the first pulse. At the base of the pulse (phase 1), the observer sees plasmoids shining in the innermost parts of the current sheet only. This region is very active in terms of the formation and merger of plasmoids, explaining the high level of variability reported above. At phases 2 and 3, more plasmoids located further away begin to contribute to the full flux. At the pulse maximum (phase 3) and beyond in the decaying part of the pulse, all plasmoids contribute within the limits of our box size. Given that the synchrotron flux is strongly peaked near the light cylinder, where the field is strongest, a larger box size would make a negligible difference. The pulse profile shows little variability in the decaying part because the total flux is averaged out by the superposition of emission from all plasmoids along the spiral, and therefore is not very sensitive to large individual plasmoids. This explains the asymmetric nature of the variability pattern in each pulse. At phase 1, the flux can sometimes be below average for an appreciable amount of time; see for example the flux minima at $t \sim 53T_{spin}$ and $t \sim 66T_{spin}$ (Figs. 3 and 4). This phenomenon is related to an episodic thickening of the current layer near the light cylinder, which temporarily prevents the layer from fragmenting and forming new plasmoids. These quiet episodes are usually followed by a recovery phase where the flux is well above average at and around phase 3 (e.g., $t \sim 54T_{spin}$), because the layer eventually breaks up at it stretches outwards, leading to a delay in the phase of the emission. This episodic thickening of the cur-

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Synthetic synchrotron pulse profiles. Top: phase–time diagram of the synchrotron flux. Vertical red dashed lines numbered 1, 2, and 3 are different phases in the left caustic, $\Phi = 0.23, 0.29$, and $0.36$, respectively, and are used in the following figures. Bottom: pulse profile seen by a single observer (orange solid line) generated by crossing the phase–time diagram along the white dashed line labeled “a” in the upper panel. The black solid line shows the pulse profile averaged over times $t \ge 50T_{spin}$. The blue shaded region is bounded by the mean flux which is augmented or diminished by the standard deviation to highlight the variable parts of the pulse profile.</figcaption></figure>

rent sheet may be a numerical artefact due to the unrealistically small scale-separation achieved in the simulation.

The top panel of Fig. 6 reports on the distribution of maximum flux, $F$, and width in pulsar phase, $\delta\Phi$, for all subpulses measured in the full dataset. The subpulse flux is measured with
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Synchrotron flux as a function of time after 50 periods at the pulsar phase $\Phi = 0.23$ labeled “1” in Fig. 3. The average flux is given by the solid black line. The shaded area shows the standard deviation with respect to the mean.</figcaption></figure>

respect to the mean profile, and is normalized to the maximum of the mean flux. Subpulses can be nearly as prominent as the average pulse profile itself. The distribution is broad, but a clear trend emerges: the brighter the subpulse is, the broader it is in phase. The proportionality relation between these two quantities is consistent with expectations for emission from plasmoids. The total synchrotron power in the subpulse is on the order of $P_{\text{sub}} \sim F \times \delta\Phi$. Considering a plasmoid located at a given radius $r$, the synchrotron power it emits will be proportional to the total number of particles it contains, and therefore to its volume. In full 3D, plasmoids become flux ropes, extending at most within the full latitudinal ranges of the stripes. For an orthogonal rotator, as modeled here, the volume of the flux rope can be as large as $V \sim \pi rw^2/4$, such that $P_{\text{sub}} \propto w^2$; thus it is proportional to the surface of the plasmoid. Using the angular size of the plasmoid as a proxy for the pulse width, $\delta\Phi \sim w/r$, yields $P_{\text{sub}} \propto \delta\Phi^2$. Therefore, we obtain $F \propto \delta\Phi$, which is in qualitative agreement with the observed correlation. The bottom panel in Fig. 6 reveals that the subpulse flux is distributed according to the power law $N_{\text{sub}} \propto F^{-1}$ up to about a few times the standard deviation, beyond which the distribution cuts off. A similar distribution is observed for the subpulse width, $N_{\text{sub}} \propto \delta\Phi^{-1}$. This distribution is also consistent with the size distribution of plasmoids, $N \propto w^{-1}$, reported in Sect. 3 (Uzdensky et al. 2010; Sironi et al. 2016).

**5. Discussion and Conclusion**

In this work, we make a case for the existence of bright subpulses in the incoherent synchrotron radiation emitted in the pulsar wind current layer (as opposed to the coherent radio emission). The origin of subpulses in the simulation is undeniably related to plasmoid formation in the layer. We can infer several clear falsifiable predictions from this work. A first prediction is that the subpulse emission should appear mostly on the leading edge of each pulse. This phenomenon stems from the caustic effect in the current sheet. Although subpulses are drifting in phase, their motion is too fast to be observed over multiple pulsar periods by the same observer, in contrast to drifting subpulses observed in the radio band (e.g., Drake & Craft 1968;

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Shining plasmoids and origin of subpulses. Top: map of the synchrotron flux from the current sheet received by the observer at phases 1 (left), 2 (middle), and 3 (right). The oblique white dotted lines are Archimedean spirals encompassing the wind current layer to guide the eye. Bottom: observed pulse profile. The red dots highlight the phases for which the flux maps are displayed.</figcaption></figure>

Weltevrede et al. 2007), suggesting there is a different mechanism at work for the coherent emission (however, see Lyubarsky 2019; Philippov et al. 2019). New subpulses should therefore appear randomly after each period. The drift is due to the relativistic motion of plasmoids traveling away from the pulsar.

Another prediction is that the subpulse flux should be proportional to its width in phase: the brighter they are, the wider they should be in phase. A broad range of pulse flux is reported, but the variability is dominated by the largest plasmoids, which can sometimes be nearly as bright as the mean pulse profile. Although our limited numerical resources do not allow us to reproduce a realistic scale separation between the smallest and the biggest islands (e.g., meter-size to hundreds of kilometers in size in the Crab pulsar), our simulation probably captures the correct demographics of the largest plasmoids, and thus the largest subpulse variability reported here. Indeed, the hierarchical merging model of plasmoids combined with the effect of a relativistic expansion of the wind predicts that the number of large plasmoids is governed by the inverse of the dimensionless
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Demographics of all synchrotron subpulses in the simulation dataset. Top: subpulse width in phase ($\delta\Phi$) and maximum flux with respect to the mean profile ($F$) distributions. The dashed oblique line is $F \propto \delta\Phi$. Bottom: subpulse flux distribution. The black dashed line shows $N_{\text{sub}} \propto F^{-1}$.</figcaption></figure>

reconnection rate, $\beta_{\text{rec}}^{-1} \ge 10$ (Cerutti & Giacinti 2021), which is consistent with the number of subpulses per pulse reported in this work. The duration of the brightest subpulses should then be $\sim\beta_{\text{rec}}^{-1}$ times shorter than the main pulses. As the reconnection rate is nearly independent of plasma parameters in the ultra-relativistic regime ($\sigma \gg 1$, Werner et al. 2018), bright subpulses should be present in all rotation-powered pulsars. We also anticipate that our results will hold in full 3D where the overall picture of flux rope formation is qualitatively similar to 2D equatorial simulations (Cerutti et al. 2020). A possible caveat is that the contribution from other latitudes in 3D may reduce the level of predicted noise.

This work provides yet another observational test for the magnetic reconnection scenario and specifically for the origin of the pulsed incoherent emission. With the intrapulse variability being dominated by at most a dozen large subpulses, we speculate that these may be detectable. The predicted emission does not fall in the radio band, but instead at higher frequencies, from the infrared to the gamma-ray bands. At these frequencies, most extrinsic sources of noise such as propagation effects through the interstellar medium are absent, meaning that any intrinsic source of noise could in principle be more easily identified than in radio.

However, a major limitation of this type of measurement will be the lack of photons and time resolution. Perhaps the most promising avenue would be optical observations of the Crab pulsar with a high time resolution (Germanà et al. 2012). X-ray observations by the NICER instruments may also provide important constraints on the intrinsic variability (Enoto et al. 2021; Hazboun et al. 2022).

**Acknowledgements.** This work received funding and support from Campus France, the French Embassy in Ankara, the COST Action PHAROS (CA16214), and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 863412). Simulations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.

## References

Arka, I., & Dubus, G. 2013, A&A, 550, A101

Belyaev, M. A. 2015, MNRAS, 449, 2759

Blandford, R., & Narayan, R. 1985, MNRAS, 213, 591

Bogovalov, S. V. 1999, A&A, 349, 1017

Cerutti, B., & Giacinti, G. 2021, A&A, 656, A91

Cerutti, B., & Philippov, A. A. 2017, A&A, 607, A134

Cerutti, B., & Werner, G. 2019, Astrophysics Source Code Library [record ascl:1911.012]

Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147

Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606

Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401

Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204

Chen, A. Y., & Beloborodov, A. M. 2014, ApJ, 795, L22

Coroniti, F. V. 1990, ApJ, 349, 538

Drake, F. D., & Craft, H. D. 1968, Nature, 220, 231

Enoto, T., Terasawa, T., Kisaka, S., et al. 2021, Science, 372, 187

Germanà, C., Zampieri, L., Barbieri, C., et al. 2012, A&A, 548, A47

Gochelashvili, K. S., & Shishov, V. I. 1975, Opt. Quant. Electr., 7, 524

Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869

Goodman, J. J., Romani, R. W., Blandford, R. D., & Narayan, R. 1987, MNRAS, 229, 73

Harding, A. K. 2022, in Astrophys. Space Sci. Lib., eds. S. Bhattacharyya, A. Papitto, & D. Bhattacharya, 465, 57

Hazboun, J. S., Crump, J., Lommen, A. N., et al. 2022, ApJ, 928, 67

Hu, R., & Beloborodov, A. M. 2021, ApJ, submitted [arXiv:2109.03935]

Kagan, D., Sironi, L., Cerutti, B., & Giannios, D. 2015, Space Sci. Rev., 191, 545

Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44

Kumamoto, H., Dai, S., Johnston, S., et al. 2021, MNRAS, 501, 4490

Loureiro, N. F., Samtaney, R., Schekochihin, A. A., & Uzdensky, D. A. 2012, Phys. Plasmas, 19, 042303

Lyne, A. G., Smith, F. G., & Graham, D. A. 1971, MNRAS, 153, 337

Lyubarskii, Y. E. 1996, A&A, 311, 172

Lyubarsky, Y. 2019, MNRAS, 483, 1731

Melrose, D. B. 1995, J. Astrophys. Astron., 16, 137

Michel, F. C. 1973, ApJ, 180, 207

Narayan, R. 1992, Phil. Trans. R. Soc. London Ser. A, 341, 151

Pétri, J. 2012, MNRAS, 424, 2023

Philippov, A. A., & Spitkovsky, A. 2014, ApJ, 785, L33

Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94

Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015, ApJ, 801, L19

Philippov, A., Uzdensky, D. A., Spitkovsky, A., & Cerutti, B. 2019, ApJ, 876, L6

Rankin, J. M., Comella, J. M., Craft, H. D., Jr, et al. 1970, ApJ, 162, 707

Rickett, B. J. 1969, Nature, 221, 158

Rickett, B. J., Coles, W. A., & Bourgois, G. 1984, A&A, 134, 390

Sieber, W. 1982, A&A, 113, 311

Sironi, L., Giannios, D., & Petropoulou, M. 2016, MNRAS, 462, 48

Sturrock, P. A. 1971, ApJ, 164, 529

Uzdensky, D. A., Loureiro, N. F., & Schekochihin, A. A. 2010, Phys. Rev. Lett., 105, 235002

Weltevrede, P., Stappers, B. W., & Edwards, R. T. 2007, A&A, 469, 607

Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840

Zelenyi, L. M., & Krasnoselskikh, V. V. 1979, Sov. Astron., 23, 460