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

HAL
open science

# Enhanced particle acceleration in a pulsar wind interacting with a companion

Valentina Richard-Romei, Benoît Cerutti

► To cite this version:

Valentina Richard-Romei, Benoît Cerutti. Enhanced particle acceleration in a pulsar wind interacting with a companion. Astronomy & Astrophysics - A&A, 2024, 689, pp.A251. 10.1051/0004-6361/202450605. hal-04700973

HAL Id: hal-04700973

https://hal.science/hal-04700973v1

Submitted on 18 Sep 2024

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
# Enhanced particle acceleration in a pulsar wind interacting with a companion

Valentina Richard-Romei* and Benoît Cerutti*

Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

Received 3 May 2024 / Accepted 16 June 2024

## ABSTRACT

**Context.** Pulsar winds have been shown to be preferred sites of particle acceleration and high-energy radiation. Numerous studies have been conducted to better characterize the general structure of such relativistic plasmas in isolated systems. However, many pulsars are found in binary systems and there are currently no ab initio models available that would include both the pulsar magnetosphere and the wind of the pulsar in interaction with a spherical companion.

**Aims.** We investigate the interaction between a pulsar wind and a companion to probe the rearrangement of the pulsar wind, assess whether it leads to an enhancement of particle acceleration, and predict the high-energy radiative signature that stems from this interaction. We consider the regime where the companion is small enough to hold between two successive stripes of the wind.

**Methods.** We performed two-dimensional (2D) equatorial particle-in-cell simulations of an inclined pulsar surrounded by a spherical, unmagnetized, perfectly conducting companion settled in its wind. Different runs correspond to different distances and sizes of the companion.

**Results.** We find that the presence of the companion significantly alters the structure of the wind. When the companion lies beyond the fast magnetosonic point, a shock is established and the perturbations are advected in a cone behind the companion. We observe an enhancement of particle acceleration due to forced reconnection as the current sheet reaches the companion surface. Hence, high-energy synchrotron radiation is also amplified. The orbital light curves display two broad peaks reaching up to 14 times the high-energy pulsed flux emitted by an isolated pulsar magnetosphere. These effects increase with the growth of the companion size and with the decrease of the pulsar-companion separation.

**Conclusions.** The present study suggests that a pulsar wind interacting with a companion induces a significant enhancement of high-energy radiation that takes the form of an orbital-modulated hollow cone of emission, which should be detectable by galactic-plane surveys, possibly with long-period radio transient counterparts.

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

## 1. Introduction

A few percent of galactic pulsars are found in binary systems (Lorimer 2008; Breton 2009). While tight binary systems involving magnetically coupled neutron stars have been studied in the past (e.g., Hansen & Lyutikov 2001; Palenzuela et al. 2013; Crinquand et al. 2019; Most & Philippov 2020), fewer studies have drawn attention to the interaction between a pulsar and a companion lying in its wind. However, it is a generic problem that has many astrophysical applications starting from binary pulsar systems, such as: (i) pulsar-neutron star such as PSR B1913+16 (Hulse & Taylor 1975) or the double pulsar system PSR J0737−3039 (Burgay et al. 2003; Lyne et al. 2004); (ii) pulsar-black hole (no confirmed detection so far; however, see Barr et al. 2024 and Chen & Liu 2024); (iii) pulsar-main sequence star, including spider pulsars such as PSR B1957+20 (Fruchter et al. 1988) or PSR J2051−0827 (Stappers et al. 1996); and (iv) pulsar-white dwarf, such as PSR J1141−6545 (Kaspi et al. 2000), PSR J1909−3744 (Jacoby et al. 2003), PSR J1738+0333 (Jacoby et al. 2007), and PSR J0348+0432 (Antoniadis et al. 2013). Planets orbiting pulsars (see Niu et al. 2022 and references therein) form another class of relevant appli-

cations from which we expect characteristic signatures, as suggested by Mishra et al. (2023), analogously to the Solar system planets and moons (e.g., Neubauer 1980 for the Jupiter-Io interaction). We should also consider the case of asteroids interacting with a pulsar wind, which have been proposed to trigger repeating fast radio bursts (Dai et al. 2016; Mottez et al. 2020; Decoene et al. 2021). As of today, the basics of such interactions remain uncertain and we are left with opened questions regarding the possible rearrangement of the magnetosphere, the strength and location of particle acceleration, and the high-energy emission originating from such systems. Ultimately, we consider whether a pulsar wind interacting with a companion may lead to a new class of long-period high-energy transients.

Modelling the interaction between a pulsar wind and a companion requires us to capture both the magnetosphere's global dynamics and the microphysical processes of the relativistic plasma responsible for particle acceleration and non-thermal radiation. To do so, we resorted to global particle-in-cell (PIC) simulations. Many studies have been conducted in order to characterize the general structure of isolated pulsar magnetospheres, whether aligned (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Cerutti et al. 2015; Belyaev 2015a; Hu & Beloborodov 2022) or inclined (Philippov et al. 2015a; Cerutti et al. 2016; Kalapotharakos et al. 2018). These models

<footnote>* Corresponding authors; valentina.richard-romei@univ-grenoble-alpes.fr, benoit.cerutti@univ-grenoble-alpes.fr</footnote>
have played a major role in the understanding of the magnetospheric structure and the kinetic processes at play. They have demonstrated that a sizeable fraction of the pulsar spindown power can be dissipated into kinetic energy. Magnetic reconnection has been identified as the predominant cause of particle acceleration and localized in the wind current sheet beyond the pulsar light cylinder. Accelerated particles escaping the reconnection layers have been shown to emit synchrotron radiation, mainly from the inner parts of the wind, thereby reproducing pulsars' magnetospheric radiative signatures (Cerutti et al. 2016; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018). To our knowledge, no simulations of a pulsar magnetosphere and a spherical companion lying in the pulsar wind have been realized, neither in the framework of magnetohydrodynamics nor in that of PIC. However, the regime where a millisecond pulsar interacts with a low-mass companion star (i.e., the spider pulsar regime) has been explored with 2D PIC simulations by focusing on the outer parts of the wind, where the pulsar wind was modeled by a plane parallel striped wind (see Cortés & Sironi 2022, 2024).

In this work, we aim to probe intermediate scales, namely, scales enabling us to simulate an inclined pulsar magnetosphere and its wind. We focus our study on a spherical unmagnetized companion settled in the pulsar wind, small enough so that its diameter is shorter than the wind stripes width, which excludes the regime of spider pulsars from our field of study. We investigate the structure of the magnetosphere, the enhancement of particle acceleration, and the high-energy emission arising from this interaction. Furthermore, we probe the impacts of the companion size and separation on these diagnostics.

In the following, we first introduce the numerical model employed in this study (Section 2). Section 3 describes our reference case, the isolated pulsar magnetosphere, recalling the main characteristics of the striped wind and the main kinetic and observational features. Our results on the pulsar-companion interaction are reported in Section 4, where we first describe the wind global dynamics before discussing particle acceleration and high-energy radiation. In Section 5, we outline the main results and their implications, and we discuss future prospects.

## 2. Numerical model

We resort to the ab-initio PIC method, which allows us to capture both the global pulsar magnetosphere as well as the kinetic processes at play. We employ the 2.5D version of the relativistic electromagnetic Zeltron code (Cerutti et al. 2013; Cerutti & Werner 2019), in spherical coordinates (Cerutti et al. 2015, 2016), restricted to the equatorial plane. We used nearly the same setup as the one presented in Cerutti & Philippov (2017), to which we added a companion in the pulsar wind. For the sake of completeness, the full setup is described below.

### 2.1. Initial setup

The computational domain is a disk made of (4096 × 4096) cells (see Fig. 1). Particles are therefore confined in the $r\phi$-plane, but we keep the three spatial components of the electromagnetic field and of the particle velocities. The spherical grid is linear in $\phi$ and logarithmic in $r$, so as to keep the cell aspect ratio constant with radius. It is well suited for pulsar winds since the particles density and the field amplitudes are decreasing functions of the radius, meaning that the relativistic plasma skin depth, $c/\omega_p$ (where $\omega_p$ is the plasma frequency) and the Larmor radius, $r_L$, increase with radius. The surface of the star, $r = r_\star$, determines the inner boundary and the box extends radially until 24

<figure><img src="image_3.png" /><figcaption>Fig. 1. Initial numerical setup. The pulsar, of a radius, $r_\star$ and angular velocity, $\Omega$, lies at the center of the spherical box of coordinates $(r, \phi)$ and radius $r_{\text{max}}$. The companion, of radius $r_{\text{comp}}$, is settled at a fixed distance $d_{\text{comp}}$ along the $\phi = 0$ direction. Initially, the box is empty of particles. Magnetic field lines (blue and red arrows) are radial, directed outwards on one hemisphere and inwards on the other hemisphere, according to the split-monopole configuration. The plane inclination in between both hemispheres is perpendicular to the magnetic moment $\mu$ shown in green. An absorbing layer is implemented at $r_{\text{absorb}} = 0.9 r_{\text{max}}$. The light-cylinder radius $R_{\text{LC}}$ is represented by the black dashed circle and the fast magnetosonic radius $r_{\text{rms}}$ is given by the red dashed circle.</figcaption></figure>

$R_{\text{LC}}$, where $R_{\text{LC}}$ is the light-cylinder radius set at $R_{\text{LC}} = 3r_\star$. The inner boundary absorbs all particles and an absorbing layer is implemented at $r_{\text{absorb}} = 0.9 r_{\text{max}}$ for both particles and fields in order to mimic an open boundary after which no information can go back inwards (Cerutti et al. 2015; Belyaev 2015b). This assumption is reasonable given that the fast magnetosonic point lies well inside the outer boundary (see Section 2.2).

The neutron star is modeled at the center of the box as a spherical perfect conductor with a constant angular velocity of $\Omega_{\text{spin}} = (c/R_{\text{LC}}) u_z$, where the unity vector, $u_z$, points in the out-of-plane direction. Initially, the star is in vacuum with magnetic field lines anchored at its surface. The implementation of a dipolar magnetic field for such equatorial configuration was proven to be inadequate (see Cerutti & Philippov 2017). We therefore replace it by a split-monopole configuration (Michel 1973; Bogovalov 1999), for which the initial magnetic field is purely radial and reverses across the plane perpendicular to the magnetic moment $\mu$, so as to ensure $\nabla \cdot \mathbf{B} = 0$ (see Fig. 1). By construction of the equatorial setup, the magnetic axis must be inclined at an angle $\chi = \pi/2$. The perfect conductor condition applied to the pulsar surface in the corotating frame implies, by Lorentz transformation, a non-zero electric field in the simulation frame:

$$E_{\star} = -\frac{(\Omega_{\text{spin}} \times r_{\star}) \times B_{\star}}{c}, \qquad (1)$$

where $E_\star$ and $B_\star$ are the fields at the surface of the star. This constraint starts the rotation of the magnetic field lines at $t = 0$.
The magnetosphere is initially empty and becomes progressively filled with electron-positron pairs that are evenly injected from the surface of the star at a rate of one macroparticle per cell per timestep for each species. They account for polar cap discharge and pair creation processes (Sturrock 1971; Ruderman & Sutherland 1975; Timokhin & Arons 2013; Philippov et al. 2020). Indeed, we neglected any other pair creation process away from the star surface, given that pairs are mainly produced in the inner magnetosphere. The presence of ions (largely dominated by pairs in terms of number density) would not modify the magnetospheric structure and their radiative signature would be largely exceeded by pair losses. This justifies our choice not to model the ion extraction from the surface. We set a high surface multiplicity of $\kappa_*=n_*/n_{GJ}=10$, where $n_*$ is the density of the plasma injected at the surface of the star and $n_{GJ} = \Omega B_*/2\pi ec$ is the Goldreich-Julian density (Goldreich & Julian 1969) at the surface of the star. This ensures that the plasma is sufficiently dense to screen the parallel electric field ($E \cdot B=0$). We also set a high plasma magnetization at the surface of the star: $\sigma_*=B_*^2/4\pi\Gamma_*n_*m_e c^2$, where $\Gamma_*$ is the surface plasma bulk Lorentz factor. Having very high multiplicity and plasma magnetization (i.e., $\kappa \gg 1$ and $\sigma \gg 1$) ensures that the quasi force-free limit is achieved in the magnetosphere and in the wind (except deviations in the current sheet), meaning that the Lorentz force dominates over any other force: the conservation of momentum equation therefore reads $\rho E + J \times B/c = 0$, where $\rho$ and $J$ are respectively the charge and current densities. Pair plasma is injected along the magnetic field lines with an initial velocity given by the drift velocity expected for the monopole solution in the force-free limit: $V_D = cE \times B/B^2$.

The scale separation in the computing box is reduced by several orders of magnitude (~$10^{-4}$) compared to a realistic pulsar to resolve the plasma kinetic scales. Our simulations are therefore best suited for millisecond pulsars. The plasma skin depth $d_e$ is resolved by 1.14 cells ($\Delta r$) at the surface of the star (i.e., at its minimum value $d_e^*$), ($d_e/\Delta r$)$_{LC} \sim 10$ at $r=R_{LC}$, and $d_e$ globally increases with radius up to a resolution of ~24 cells at the outer boundary. The current sheet that forms in the magnetosphere (see Section 3) has a width of the order of $d_e$. At $r=R_{LC}$, the Larmor radius resolution goes from ~1 cell in the wind to ~70 cells inside the current sheet. The Larmor radius resolution then gradually increases with radius up to ~10 cells in the wind at the edge of the box.

## 2.2. Companion implementation

We added a companion in the wind of the central pulsar (see Fig. 1). For simplification purposes, the companion is chosen to be an unmagnetized perfect conductor, with no wind and no intrinsic spin. Particles hitting its surface are absorbed. We aim to study the impact of the binary separation and of the companion radius on the pulsar magnetosphere and wind. Therefore, we ran six simulations for different choices of separations ($d_{comp}$) and radii ($r_{comp}$). We refer to Table 2 for the full set of parameters. The pulsar wind is globally made of two nested Archimedean-shaped stripes of alternating magnetic polarity, with a wavelength of $\lambda = 2\pi R_{LC}$, separated by current sheets (see Section 3). We restricted our study to companion sizes ensuring

$$ \delta_{cs} < r_{comp} < \pi R_{LC}, \qquad (2) $$

where $\delta_{cs}$ is the current sheet width. Indeed, we want the companion radius to be smaller than the semi-stripe wavelength so that the whole diameter holds between two successive reconnection layers. Such choice excludes the spider pulsar regime. Given

our choice of companion separations and radii (see Table 2), the $\delta_{cs}/r_{comp}$ ratio varies from 0.15 (run D2R1) to 0.9 (run D9R05).

We probe different zones of the wind by changing the separation of the companion $d_{comp}$. In particular, we wish to investigate the impact of the companion position with respect to the fast magnetosonic point, $r_{fms}$. This point is defined as the radius for which the wind velocity exceeds the Alfvén speed (Michel 1969; Kirk et al. 2009):

$$ \Gamma_{fms} = \left( \frac{B^2}{4\pi nm_e c^2} \right)^{1/3} = (\Gamma \sigma)^{1/3}, \qquad (3) $$

where $\sigma$ is the plasma magnetization. We need the fast magnetosonic point to be far enough from the light cylinder radius as well as from the outer boundary of the box. This constraint leads us to fix a magnetization of $\sigma_* = 250$ at the neutron star surface, meaning $\sigma_{LC} \sim 60$ at the light-cylinder radius outside the current sheet.

The pulsar-companion separation is assumed to be constant in time. Indeed, even considering a compact object binary made of millisecond pulsars and for our closest separation (i.e., $d_{comp} = 2R_{LC}$), the merger time due to the emission of gravitational waves is 300 times longer than the pulsars spin period ($P_{spin}$) and 30 times longer than their orbital period ($P_{orb}$), according to the approximate analytical expression of merger time derived in Iorio et al. (2023). The companion is settled at rest in the simulation. Even for the closest binary separation ($d_{comp} = 2R_{LC}$), the orbital period computed considering the Keplerian orbit of a millisecond pulsar binary is 10 times longer than the pulsar spin period ($P_{spin}$). This assumption, added to the perfect conductor condition, implies that $E=0$ at the surface and everywhere inside the companion.

## 2.3. Fields and particle evolution

Starting from the initial conditions given above, the PIC method allows us to couple the particles evolution along with time-dependent electromagnetic fields in a self-consistent way. Firstly, the Boris push (Birdsall & Langdon 1991) is used to evolve particles positions and velocities according to the Abraham-Lorentz-Dirac equation. This equation accounts for radiative energy losses by adding a radiation-reaction force ($f_{rad}$) to the Lorentz force:

$$ \frac{d(\gamma m_e v)}{dt} = q(E + \beta \times B) + f_{rad}, \qquad (4) $$

where, for each particle, $v = \beta c$ is the 3-velocity, $q$ is the electric charge, and $\gamma$ is the Lorentz factor. The radiation-reaction force is implemented according to the Landau-Lifshitz formula (Landau & Lifshitz 1971), in the Tamburini et al. (2010) approximation (see Cerutti et al. 2016 for full details) as:

$$ f_{rad} = \frac{2}{3} r_e^2 [(E + \beta \times B) \times B + (\beta \cdot E)E] - \frac{2}{3} r_e^2 \gamma^2 [(E + \beta \times B)^2 - (\beta \cdot E)^2] \beta, \qquad (5) $$

where $r_e = e^2/m_e c^2$ is the classical electron radius. Due to numerical costs, each simulated particle represents a large number of real particles, given by the weight, $w_k$, with the same $q/m$ ratio and therefore following the same trajectory in phase space. Secondly, charges and currents are deposited on the grid, based on an area-weighting deposition scheme. Thirdly, electromagnetic fields are evolved on the grid by solving Maxwell-Faraday
<table><caption>**Table 1.** Simulation parameters. Values evaluated at $R_{LC}$ or $r_{fms}$ are computed for the isolated pulsar magnetosphere.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Number of cells</td><td>4096 (r) × 4096 (φ)</td></tr><tr><td>Inner boundary</td><td>r<sub>★</sub></td></tr><tr><td>R<sub>LC</sub></td><td>3 r<sub>★</sub></td></tr><tr><td>r<sub>absorb</sub></td><td>24 R<sub>LC</sub></td></tr><tr><td>(d<sub>e</sub>/Δr)<sub>LC</sub></td><td>16.2</td></tr><tr><td>σ<sub>LC</sub></td><td>60</td></tr><tr><td>P<sub>spin</sub>/Δt</td><td>4.3 × 10<sup>4</sup></td></tr><tr><td>r<sub>fms</sub></td><td>5.1 R<sub>LC</sub></td></tr><tr><td>Γ<sub>fms</sub></td><td>3.9</td></tr><tr><td>d<sub>e</sub><sup>★</sup>/r<sub>★</sub></td><td>1.8 × 10<sup>-3</sup></td></tr><tr><td>Plasma composition</td><td>electrons and positrons</td></tr><tr><td>Injection model</td><td>from the star surface</td></tr></tbody></table>

and Maxwell-Ampère equations through the finite-difference-time-domain algorithm (Yee 1966). While $\nabla \cdot \mathbf{B} = 0$ is automatically verified under these conditions, the conservation of charge is not ensured, due to machine truncation errors. We then solve the Poisson equation with the iterative Gauss-Seidel method every 25 timesteps ($\Delta t$). The timestep is given by half the Courant-Friedrich-Lewy critical condition. We note that general relativistic corrections on the electrodynamics of the system are neglected in this work (see however Philippov et al. 2015a; Philippov & Spitkovsky 2018; Chen et al. 2020; Torres et al. 2024). The full list of numerical and physical parameters chosen for the simulations is given in Table 1.

## 2.4. Radiation modeling

Synchrotron radiation is the main source of high-energy emission beyond the light cylinder (Cerutti et al. 2016). Each macroparticle of the simulation emits a macrophoton, representing a set of physical photons and emitting the following power (Blumenthal & Gould 1970):

$$P_{\text{rad}} = \frac{2}{3} r_{e}^{2} c \gamma^{2} \tilde{B}_{\perp}^{2}, \qquad (6)$$

with power spectrum

$$\frac{dP_{\text{rad}}}{d\nu} = \frac{\sqrt{3}e^3 \tilde{B}_{\perp}}{m_e c^2} \left(\frac{\nu}{v_c}\right) \int_{v/v_c}^{+\infty} K_{5/3}(x)dx, \qquad (7)$$

where $K_{5/3}$ is the modified Bessel function of order $5/3$, $\nu$ is the radiation frequency, $v_c = 3e\tilde{B}_{\perp}\gamma^2/4\pi m_e c$ is the critical frequency, and $\tilde{B}_{\perp}$ is defined in Eq. (8). Around pulsars, the electric and magnetic field strengths are comparable. Instead of referring to the classical expression $B_{\perp} = B \sin \alpha$ (where $\alpha$ is the angle between **B** and the direction of the particle), valid for synchrotron radiation in a pure magnetic field, we use the effective perpendicular magnetic field strength, $\tilde{B}_{\perp}$, that takes into account an arbitrary electromagnetic field (see Cerutti et al. 2016 for further details):

$$\tilde{B}_{\perp} = \sqrt{(E + \beta \times \mathbf{B})^2 - (\beta \cdot E)^2}. \qquad (8)$$

Quantum electrodynamical effects reached above the critical magnetic field, $B_{QED} = m_e^2 c^3 / \hbar e$, are neglected since we globally have $\gamma \tilde{B}_{\perp} \ll B_{QED}$ in the magnetosphere. Photons are emitted along the particles momentum. This is a good approximation

in the presence of strong relativistic beaming since the emission cone has a semi-aperture angle of $\sim 1/\gamma \ll 1$. Photons then propagate freely at the speed of light, without interacting between themselves nor with the magnetic field. However, they are absorbed if they hit the star or the companion. General relativistic effects on the photon trajectories are neglected. We collect the photons on a screen at infinity.

In this problem, we focus on the orbital modulation of the light curves. As previously mentioned, the companion is at rest in the simulation, under the assumption that $P_{\text{orb}} \gg P_{\text{spin}}$. However, assuming a circular orbit, we can model orbital-modulated light curves by placing observers all around the box. An isolated pulsar emits a pulsed radiation over its spin period. To compute such pulsed light curves, we need to take into account the time delay between photons emitted at different locations of the box when they reach the screen (Cerutti et al. 2016). Nevertheless, when we focus on orbital-modulated light curves, the delay times are way shorter than the orbital timescales and we do not need to consider them.

## 3. Reference case: Isolated pulsar magnetosphere

The initial split-monopole condition represented in Figure 1 quickly evolves as the magnetic field lines start rotating at $t > 0$ to ensure the perfect conductor condition of the star (see Section 2.1). Pairs injected from the surface gradually fill the box. The magnetosphere converges to a stationary solution after $\sim 4.5 P_{\text{spin}}$ (Fig. 2). Within the light cylinder, the magnetic equatorial regions are characterized by closed field lines that co-rotate with the star and trap the plasma, whereas open field lines of opposite polarities originating from the polar caps escape the light cylinder due to the rotation of the star and reconnect along the current sheet (Kirk et al. 2009; Arons 2012; Pétri 2016; Cerutti & Beloborodov 2017). When the magnetic axis is inclined with respect to the rotation axis (recall $\chi = \pi/2$ here), the current sheet takes the shape of an oscillatory structure in $\theta$, with wavelength $2\pi R_{LC}$ and angular aperture $2\chi$, referred to as the “striped wind” (Coroniti 1990; Bogovalov 1999). A 2D equatorial cut of the current sheet results in two Archimedean spirals separated by stripes $\pi R_{LC}$ wide of alternating magnetic field polarities, conveying the outflowing cold relativistic wind (see Fig. 2a).

Figure 2c shows a snapshot of the bulk Lorentz factor of the wind at $t = 12.4 P_{\text{spin}}$. On average, the bulk Lorentz factor of the wind increases almost linearly until it reaches $\Gamma = 3.9$ at the fast magnetosonic point, $r_{\text{FMS}} = 5.1 R_{LC}$. After this point, the wind continues to accelerate at a slower rate up to $\Gamma \sim 7.5$ at the box edge. Due to magnetic reconnection locally attracting particles, the wind is slowed down on the leading edge of the current sheet whereas its acceleration increases on the trailing edge of the spiral. Right after its formation near the light cylinder, the current sheet fragments due to the relativistic tearing instability (Zelenyi & Krasnoselskikh 1979), giving rise to a chain of plasma overdensities confined in magnetic islands, called plasmoids (Cerutti & Philippov 2017; Cerutti et al. 2020; Hakobyan et al. 2023). Plasmoids then gradually grow by merging with each other while flowing outwards along the spirals. In between plasmoids, short current layers (referred to as X-points) allow for relativistic magnetic reconnection (Uzdensky et al. 2010), which is the main physical process responsible for magnetic energy dissipation into particles kinetic energy (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Philippov et al. 2015b; Cerutti et al. 2015, 2020; Cerutti & Philippov 2017; Belyaev 2015).
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Snapshot of the isolated pulsar simulation after 12.4 pulsar spin periods. The central pulsar is represented by a black disk on left panels and a white disk on right panels, for better visualization purposes. *From top-left to bottom-right:* toroidal magnetic field ($B_\phi/B_\star$)($R_{LC}r/r_\star^2$), where solid lines are the magnetic field lines; the mean Lorentz factor of particles ($\langle\gamma\rangle$); the bulk Lorentz factor of the plasma $\Gamma$; and the isotropic synchrotron power emitted by the particles per unit of volume ($dP_{rad}/dV(r^2/r_\star^2)$). Radii are given in units of the light cylinder radius, represented by the red circle.</figcaption></figure>

2015a; Hu & Beloborodov 2022). About 24% of the Poynting flux reservoir (i.e., the spindown power extracted from the star) is consumed via magnetic reconnection and converted into kinetic energy between the light cylinder and the outer part of the box. Reconnection is particularly efficient at accelerating particles for $r < 2R_{LC}$. Accelerated particles escaping the X-points are trapped by plasmons. Figure 2b shows the mean particle acceleration ($\langle\langle\gamma\rangle\rangle$).

A significant fraction of kinetic energy is radiated away through synchrotron emission (see Fig. 2d). As expected,

non-thermal radiation is emitted from the current sheet and decays with radius given that the density and $\tilde{B}_\perp^2$ decrease with radius as $1/r^2$. In the end, 0.7% of the spindown power is converted into high-energy radiation. Due to time delay effects, the synchrotron radiation is pulsed, with two short bright pulses per stellar spin period (Lyubarskii 1996; Pétri 2012; Arka & Dubus 2013; Cerutti et al. 2016; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018, 2023). We defined a “reference flux,” computed by averaging this pulsed radiation over one pulsar spin period, to later compare it with the orbital-modulated
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Snapshot of the D9R1 run after 11.7 pulsar spin periods. The central pulsar and the companion (of distance $d_{\text{comp}} = 9R_{\text{LC}}$ and size $r_{\text{comp}} = r_\star$) are represented by black disks on left panels and white disks on right panels, for better visualization purposes. *From top left to bottom right:* toroidal magnetic field ($B_\phi/B_\star$($R_{\text{LC}}r/r_\star^2$), where solid lines are the magnetic field lines; the mean Lorentz factor of particles, $\langle\gamma\rangle$; the bulk velocity of the wind, $\Gamma$; the isotropic synchrotron power emitted by the particles per volume unit ($dP_{\text{rad}}/dV$(r$^2/r_\star^2$). Radii are given in units of the light cylinder radius, represented by the red circle. The white dashed lines indicate the aperture angle of the shocked zone appearing on panel (c); we report them in panel (d) for comparison with the emitting hollow cone.</figcaption></figure>

light curves computed in the presence of the companion (see Section 4.3).

**4. Pulsar-companion interaction**

We carry out a parametric study in order to probe the impact of the companion on the pulsar wind, depending on the binary separation ($d_{\text{comp}}$) and on the companion size ($r_{\text{comp}}$). Table 2 details the parameters chosen for the six simulations we ran.

**4.1. Wind dynamics**

Figure 3 presents the same diagnostics as the ones we show for the isolated pulsar (Section 3), but in the presence of a companion at $d_{\text{comp}} = 9R_{\text{LC}}$ and $r_{\text{comp}} = r_\star$ (run D9R1). The conductor is settled beyond the fast magnetosonic point ($r_{\text{fms}} \sim 5R_{\text{LC}}$) for this run. We observe an alteration of the magnetosphere (Fig. 3, panel a), but (as expected) the perturbations are advected by the pulsar wind and, therefore, they stay in a cone behind the
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Bulk Lorentz factor of the wind ($\Gamma_{wind}$) averaged over several spin periods, for each run. Radii are given in units of the light-cylinder radius, represented by the red circle. The fast magnetosonic surface is represented by the dashed green curve. When $d_{comp} > r_{fms}$, a shock of angular aperture $\zeta$ (see upper right panel) is established. Otherwise, no shock appears and the whole magnetosphere is slowed down.</figcaption></figure>

<table><caption>**Table 2.** Parameters chosen for each run.</caption><thead><tr><td>Run</td><td>d<sub>comp</sub> (R<sub>LC</sub>)</td><td>r<sub>comp</sub> (r<sub>★</sub>)</td></tr></thead><tbody><tr><td>D2R1</td><td>2</td><td>1</td></tr><tr><td>D5R1</td><td>5</td><td>1</td></tr><tr><td>D9R05</td><td>9</td><td>0.5</td></tr><tr><td>D9R1</td><td>9</td><td>1</td></tr><tr><td>D9R2</td><td>9</td><td>2</td></tr><tr><td>D15R1</td><td>15</td><td>1</td></tr></tbody></table>

On the contrary, when $d_{comp} \le r_{fms}$ (runs D2R1 and D5R1), perturbations propagate faster than the outflowing wind and the whole magnetosphere is affected. As the wind reaches the companion, its bulk Lorentz factor sharply drops along a wider cone, indicating the presence of a shock (see Fig. 3c). Given the strong magnetization of the magnetosphere, we did not expect large discontinuities in the magnetic field, nor in the density, that would provide evidence of a significant compression (Kennel & Coroniti 1984).

Figure 4 presents, for each of our simulations, the bulk Lorentz factor of the wind averaged over several spin periods of the pulsar, to which we added the fast magnetosonic surfaces. Averaging over several spin periods enables us to get rid of the wind striped structure and of local inhomogeneities, so as to keep only the overall bulk motion. We see that for pulsar-

companion separations larger than the radius of the fast magnetosonic surface (lying at $r_{fms} \sim 5R_{LC}$ according to the pulsar parameters set in this study), a shock is systematically established. However, when the companion is settled before the fast magnetosonic surface $d_{comp} \le 5r_{fms}$ (i.e., for runs D2R1 and D5R1), no shock appears. The wind is slowed down over the whole box, more isotropically and more intensely as the orbital separation decreases. We notice that, in this case, the fast magnetosonic surface gets disrupted but keeps encompassing the companion. The local enhancement of $\Gamma$ that appears in the wake of the companion is due to artefacts of the very low plasma density.

When a shock is formed in a plane-parallel and uniform flow, the cone aperture angle ($\zeta$) can be related to the relativistic Alfvénic Mach number ($M$) computed at the apex of the cone via the following relation:

$$ \sin\left(\frac{\zeta}{2}\right) \sim \frac{1}{M} \qquad (9) $$

where $M = \beta\Gamma/\beta_A\Gamma_A$, given that $\beta_A = \sqrt{\sigma/(1+\sigma)}$ (Kirk et al. 2009) and $\Gamma_A$ are the Alfvén velocity and the Alfvén Lorentz factor, respectively. The measured values from the simulation (see $\zeta$ on the top right panel of Fig. 4) are in reasonably good agreement with the law presented in Eq. (9), with relative errors of less than 28%. The concordance is particularly remarkable given the high inhomogeneity of the striped wind and the tilt in the cone direction, which lead to different wind velocities reached for a same distance from the apex of the cone. Increasing the
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Synchrotron power per unit of volume ($dP_{\text{rad}}/dV$) on successive snapshots of run D9R1, zoomed-in around the companion, during a current sheet-companion interaction. The white dashed line is an Archimedean spiral tracing the geometry that the current sheet (the one we focus on) would adopt in the absence of the companion. The current sheet flows almost radially at $d_{\text{comp}} = 9R_{\text{LC}}$. It bends and lights up when reaching the companion surface and finally breaks apart.</figcaption></figure>

companion radius ($r_{\text{comp}}$) leads to an increase of the cone aperture angle ($\zeta$). The altered part of the magnetosphere (Fig. 3 panel a) displaying enhanced particle acceleration and radiation lies inside the shocked zone, so that the aperture angle of the shocked cone also gives an upper boundary to the aperture of the emitting cone (see Fig. 3d).

## 4.2. Particle kinetics

The presence of the companion in the magnetosphere of the pulsar enhances particle acceleration regardless of the binary separation. This comes from the compression of the magnetic field lines reaching the companion, which results in a forced reconnection. Indeed, as can be seen in Figure 5, the forward part of the current sheet is first slowed down while approaching the companion (starting from a distance of about $1r_{\text{comp}}$). It then bends backwards but ends up reaching the companion surface as the current sheet continues to progress radially. After the companion surface is reached, the current sheet eventually breaks apart around the obstacle and the two branches of enhanced particle acceleration flow outwards, creating two radial lines of enhanced particle acceleration. In fact, the whole cone between these two lines represents a favorable zone for particle acceleration (see Fig. 3 panel b). However, the density is extremely low in the wake of the companion, so that most of the highly accelerated particles flow along the borders of the cone. This can be seen in Fig. 3d representing the emitted radiation power $P_{\text{rad}}$, which is a good tracer of particle acceleration since the mean Lorentz factor ($\langle\gamma\rangle$) is weighted by the local density. Particle acceleration continues to increase with radius until the end of the box, since magnetic reconnection continues to operate and radiative losses diminish with radius.

Figure 6 shows the particles energy spectra for all runs. The presence of the companion induces a localized bump in the spectrum, centered on different energies depending on the binary separation. While the separation of $2R_{\text{LC}}$ leads to a bump centered on $\gamma \sim 17$, the separation of $5R_{\text{LC}}$ seems to be the most favourable to reach the highest energies, with the bump being centered on $\gamma \sim 200$. The radius of the companion does not seem to play a significant role on the particles energy spectra (see orange lines on Fig. 6). A domain decomposition based on the shocked zone (see Fig. 7) confirms that the observed hardening of the spectra entirely comes from the shocked part of the wind while, outside the shocked surface, the spectral shape is the same as for an isolated pulsar.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Particles energy spectra for all runs, averaged over several spin periods. Spectrum of the isolated pulsar magnetosphere (black solid line) and spectra for all binary configurations (coloured solid lines). Orange dashed and dotted lines stand for different companion sizes.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Domain decomposition of the total particles energy spectrum (blue solid line) for the D9R1 run. Spectrum computed inside the shocked zone (blue dashed line) and spectrum computed outside of the shocked zone (blue dotted line). Comparison with the isolated pulsar case (black solid line). The spectrum alteration in the companion presence exclusively comes from the shocked zone.</figcaption></figure>

## 4.3. Electromagnetic signature

As previously mentioned (Section 2.4), synchrotron losses account for most of the high-energy radiation in the studied system. As can be seen in Figure 2 for the isolated pulsar configuration, synchrotron losses originate from the particles
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Synchrotron spectra for all runs, averaged over several spin periods. Spectrum of the isolated pulsar magnetosphere (black solid line) and spectra for all binary configurations (coloured solid lines). Orange dashed and dotted lines stand for different companion sizes.</figcaption></figure>

previously accelerated in the reconnection layers and confined inside the plasmons along the current sheet. In particular, synchrotron emission is predominant in the inner parts of the magnetosphere where the magnetic field is stronger ($B \propto 1/r$) and the density higher ($n \propto 1/r^2$).

However, the presence of a companion adds a predominant contribution to the usual synchrotron losses. Indeed, when the current sheet hits the conductor surface, particles are further accelerated (see previous section) and the current sheet lights up again (see Fig. 5). The current sheet is torn apart by the companion and the two separated branches resulting from this interaction flow radially beyond the companion, creating a hollow cone of light (as shown in Figure 3, panel d).

Figure 8 compares the synchrotron spectra for all runs. As expected, for all binary separations, the presence of the companion significantly enhances synchrotron losses and shifts the peaks of the spectra to higher frequencies. Here, again, the separation $d_{\text{comp}} = 5R_{\text{LC}}$ seems to optimise the emission of synchrotron radiation peaking at a higher frequency compared to the other runs. The excess of synchrotron emission decreases with the binary separation. This mainly comes from the decrease of the density and the magnetic field strength with radius, despite the slower growth of $\gamma$ with radius. We notice a slight upward shift of the spectrum when doubling the companion size ($r_{\text{comp}} = 2r_{\star}$), and a slight downward shift when decreasing the companion radius ($r_{\text{comp}} = 0.5r_{\star}$). Figure 9 makes a comparison, for the run D9R1, between the synchrotron spectrum emitted by the shocked part of the magnetosphere and the one emitted by the unshocked part. The unshocked part of the magnetosphere displays the same spectrum as the isolated pulsar magnetosphere, indicating that the additional contribution exclusively comes from the shocked region of the magnetosphere.

Orbital light curves for all runs are shown in Figure 10. We represent the high-energy synchrotron emission taken above the fiducial synchrotron frequency $v_0 = 3eB_*/4\pi mc$, which was proven in previous studies to be an appropriate threshold for studying synchrotron emission from the current sheet (see Cerutti et al. 2016). Light curves are integrated over the full polar angle range and are normalized by the reference flux, computed by averaging over $P_{\text{spin}}$ the pulsed radiation emitted by the isolated pulsar magnetosphere (see Section 3). All light curves present a significant enhancement of the radiation flux, with two broad peaks corresponding to the passing of the observer's line of sight through the edges of the hollow cone of emission (see Fig. 3d). The peak height varies from 3.5 to 14 times the ref-

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Domain decomposition of the total synchrotron spectrum (blue solid line) for the D9R1 run. Spectrum computed inside the shocked zone (blue dashed line) and spectrum computed outside of the shocked zone (blue dotted line). Comparison with the isolated pulsar case (black solid line). The spectrum alteration in the companion presence exclusively comes from the shocked zone.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** High-energy ($v > v_0 = 3eB_*/4\pi mc$) synchrotron light curves for all runs, as a function of the orbital phase, $\Phi$. Fluxes are integrated over the full polar angle range. Fluxes are normalized by the reference flux computed for the isolated pulsar magnetosphere (black solid line). Different colours stand for different binary separations. Orange dashed and dotted lines stand for different companion sizes. The pulsar-companion direction is given by the black dashed line. Two intense and broad peaks (corresponding to a hollow cone of light in 3D) are induced by the pulsar-companion interaction, which properties depend on the companion distance and size. The angle between peaks is referred to as $\Delta\Phi$ and the light cone tilt with respect to the pulsar-companion direction is designated by $\delta_{\text{offset}}$.</figcaption></figure>

ference flux, depending on the companion separation and size. Higher companion sizes (orange lines on Fig. 10) induce higher peak intensities and higher orbital phases between peaks. The binary separation has an impact on the intensity of the peaks, their widths, the separation between them ($\Delta\Phi$), and the light cone orientation with respect to the pulsar-companion direction ($\delta_{\text{offset}}$).

The exact dependency between these parameters and the pulsar-companion separation is presented in Figure 11. The upper panel shows that the total synchrotron power emitted over an orbit scales as $d_{\text{comp}}^{-0.56}$ (red line). The orbital phases between both peaks (light blue points on the middle panel of Fig. 11) decrease with radius as $d_{\text{comp}}^{-0.68}$ (light blue line). We note that they scale similarly to the corresponding $1/\Gamma$ factors, computed for each run at the apex of the shocked cone (dark blue points on the middle panel of Fig. 11), suggesting that the peak
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Light curves properties as a function of the pulsar-companion separation. Upper panel: integrated synchrotron power over the whole orbit (red points), normalized by the power computed for the isolated pulsar. The total synchrotron power scales as $d_{\text{comp}}^{-0.56}$ (red line). Middle panel: the orbital phase difference between peaks (light blue points) scales as $d_{\text{comp}}^{-0.68}$ (light blue line). We plot the corresponding $1/\Gamma$ factors normalized by a factor 0.55, computed before the companion for each run (dark blue crosses). Lower panel: angular offsets of the emitting cone with respect to the pulsar-companion direction (green points) and ratio of the wind drift velocity components $V_{\phi}/V_{\Gamma} = R_{\text{LC}}/r$ (green line) given by the monopole analytical solution (Michel 1973).</figcaption></figure>

separation is shaped by relativistic beaming effects. In Figure 11, the points corresponding to the $1/\Gamma$ scaling are normalized by a factor of 0.55. This discrepancy can be explained by the values of $\Gamma$ inferred from the averaged maps (shown in Fig. 4), which erase any inhomogeneities. Indeed, as can be seen in Figure 3c, $\Gamma$ can locally be amplified by up to a factor of 2 inside the current sheet near the companion. The lower panel demonstrates that the angular offsets of the light hollow cones with respect to the pulsar-companion direction are determined by the plasma drift velocity at the companion radius. Indeed, the values of $\tan(\delta_{\text{offset}})$ for each companion separation (green crosses) are given by the corresponding drift velocity ratios $V_{\phi}/V_{\Gamma} = R_{\text{LC}}/r$ (green line) taken at the companion location and computed according to the monopole analytical solution (derived in Michel 1973). We note that the shocked cones represented in Figure 4 are tilted, for each run, by exactly the corresponding angle $\delta_{\text{offset}}$.

**5. Discussion and conclusion**

In this work, we study the interaction of a pulsar wind with a similarly sized companion to assess the extent to which the pulsar magnetosphere reshapes itself. We probe the enhancement and location of particle acceleration to predict the non-thermal radiation originating from such systems. We show that the interaction between a pulsar wind and a companion alters significantly the dynamical and energetic properties of the wind. The pulsar wind is slowed down rather isotropically if the companion lies within the fast magnetosonic point. Otherwise, when the companion is settled beyond the fast magnetosonic point, a shock is established and all the perturbations are advected in a cone behind the companion. Each time the outflowing wind current sheet impacts the companion surface, we observe a forced reconnection that leads to a significant enhancement of particle acceleration. Reaccelerated particles form a hollow cone behind the

companion. By doing so, they induce an orbital-modulated hollow cone of high-energy synchrotron radiation, whose intensity and aperture depend on the orbital separation and the size of the companion. Hence, non-thermal radiation from such systems is significantly enhanced as well, compared to the reference pulsed flux of an isolated pulsar magnetosphere.

The 2D framework used here may introduce some artifacts. First, the magnetic field is confined in the equatorial plane, which may have an impact on the rearrangement of the field lines. It is worth noting that our results give an upper limit in terms of particle acceleration and high-energy emission. Indeed, the equatorial configuration necessarily implies that all the magnetic field lines accumulating ahead of the companion reconnect, possibly overestimating the amount of dissipation into kinetic energy. In addition, in 2D, particles accelerated by magnetic reconnection escape the X-points and are all trapped by the neighbouring plasmons, however, a third dimension would allow for the leakage of the accelerated particles outside of plasmons (Zhang et al. 2021; Chernoglazov et al. 2023). This could slightly alter the shape of the light hollow cone. As mentioned at the beginning of this study, the scale separation that we have been able to reach in these PIC simulations is reduced by several orders of magnitude compared to a realistic pulsar magnetosphere. This numerical limitation constrains the frequency range of the high-energy radiation, leaving many uncertainties on the spectral limits and the spectral shape – but not on the total energy flux that is conserved. The presence of a companion wind would imply a bigger effective surface, but should not radically alter the nature of the shock. A slight asymmetry between the two peaks appears in the light curves. Taking into account the orbital motion of the companion could potentially exacerbate this asymmetry.

We predict an orbital-modulated emission that could originate from the orbital motion of asteroids, planets, or neutron stars around the pulsar, with a higher chance of observing such systems if seen close to edge-on. The characteristic frequency of the non-thermal radiation is expected to fall in the soft gamma-ray band. We therefore expect these transients to be observable only on galactic distances. While not investigated in this work, radio counterparts are expected in addition to the high-energy non-thermal radiation. Indeed, in the wake of the companion where the wind is altered, we observe the propagation of small-scale fast magnetosonic modes. According to Lyubarsky (2019) and Philippov et al. (2019), the collision of plasmons with each other and with the magnetic field perturbates the magnetic field and induces short fast magnetosonic pulses that eventually escape the plasma as radio waves. Other coherent mechanisms of radio emission could also be at play, such as the cyclotron maser instability that was initially retained to describe fast radio bursts due to an outflowing pulsar wind crossing Alfvén wings (Neubauer 1980; Mottez & Heyvaerts 2011) formed behind a companion (Heyvaerts et al. 2012; Mottez & Heyvaerts 2012; Mottez & Zarka 2014; Mottez et al. 2020; Decoene et al. 2021). Galactic radio counterparts would be of significant interest, especially in the light of the recently discovered galactic long-period radio transients (Hurley-Walker et al. 2022, 2023). These periodic radio transients are characterized by very long periods, ranging from 10s to 1000s. Their sources have not been uniquely identified yet; however, binary systems including a pulsar have been noted as possible source candidates (Rea et al. 2022, 2024).

**Acknowledgements.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 863412). Computing resources were provided by TGCC under the allocation A0150407669 made by GENCI.
References

Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, *Science*, **340**, 448

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

Arons, J. 2012, *Space Sci. Rev.*, **173**, 341

Barr, E. D., Dutta, A., Freire, P. C. C., et al. 2024, *Science*, **383**, 275

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

Belyaev, M. A. 2015b, *Nature Astron.*, **3**, 37

Birdsall, C. K., & Langdon, A. B. 1991, *Plasma Physics via Computer Simulation* (Bristol: Adam Hilger)

Blumenthal, G. R., & Gould, R. J. 1970, *Rev. Mod. Phys.*, **42**, 237

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

Breton, R. P. 2009, *PhD Thesis*, McGill University, Canada

Burgay, M., D'Amico, N., Possenti, A., et al. 2003, *Nature*, **426**, 531

Cerutti, B., & Beloborodov, A. M. 2017, *Space Sci. Rev.*, **207**, 111

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

Chen, A. Y., Cruz, F., & Spitkovsky, A. 2020, ApJ, **889**, 69

Chen, Z.-C., & Liu, L. 2024, arXiv e-prints [arXiv:2401.12889]

Chernoglazov, A., Hakobyan, H., & Philippov, A. 2023, ApJ, **959**, 122

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

Cortés, J., & Sironi, L. 2022, ApJ, **933**, 140

Cortés, J., & Sironi, L. 2024, arXiv e-prints [arXiv:2404.03700]

Crinquand, B., Cerutti, B., & Dubus, G. 2019, *A&A*, **622**, A161

Dai, Z. G., Wang, J. S., Wu, X. F., & Huang, Y. F. 2016, ApJ, **829**, 27

Decoene, V., Kotera, K., & Silk, J. 2021, *A&A*, **645**, A122

Fruchter, A. S., Stinebring, D. R., & Taylor, J. H. 1988, Nature, **333**, 237

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

Hakobyan, H., Philippov, A., & Spitkovsky, A. 2023, ApJ, **943**, 105

Hansen, B. M. S., & Lyutikov, M. 2001, *MNRAS*, **322**, 695

Heyvaerts, J., Lehner, T., & Mottez, F. 2012, *A&A*, **542**, A128

Hu, R., & Beloborodov, A. M. 2022, ApJ, **939**, 42

Hulse, R. A., & Taylor, J. H. 1975, ApJ, **195**, L51

Hurley-Walker, N., Zhang, X., Bahramian, A., et al. 2022, *Nature*, **601**, 526

Hurley-Walker, N., Rea, N., McSweeney, S. J., et al. 2023, *Naure*, **619**, 487

Iorio, G., Mapelli, M., Costa, G., et al. 2023, *MNRAS*, **524**, 426

Jacoby, B. A., Bailes, M., Ord, S. M., Knight, H. S., & Hotan, A. W. 2007, ApJ, **656**, 408

Jacoby, B. A., Bailes, M., van Kerkwijk, M. H., et al. 2003, ApJ, **599**, L99

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

Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., & Kazanas, D. 2023, ApJ, **954**, 204

Kaspi, V. M., Lyne, A. G., Manchester, R. N., et al. 2000, ApJ, **543**, 321

Kennel, C. F., & Coroniti, F. V. 1984, ApJ, **283**, 694

Kirk, J. G., Lyubarsky, Y., & Petri, J. 2009, Astrophys. Space Sci. Lib., **357**, 421

Landau, L. D., & Lifshitz, E. M. 1971, *The Classical Theory of Fields* (Pergamon Press)

Lorimer, D. R. 2008, *Liv. Rev. Relat.*, **11**, 8

Lyne, A. G., Burgay, M., Kramer, M., et al. 2004, *Science*, **303**, 1153

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

Lyubarsky, Y. 2019, *MNRAS*, **483**, 1731

Michel, F. C. 1969, ApJ, **158**, 727

Michel, F. C. 1973, ApJ, **180**, L133

Mishra, R., Čemeljić, M., Varela, J. & Falanga, M. 2023, ApJ, **959**, L13

Most, E. R., & Philippov, A. A. 2020, ApJ, **893**, L6

Mottez, F., & Heyvaerts, J. 2011, *A&A*, **532**, A21

Mottez, F., & Heyvaerts, J. 2012, in SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. S. Boissier, P. de Laverny, N. Nardetto, et al., 597

Mottez, F., & Zarka, P. 2014, *A&A*, **569**, A86

Mottez, F., Zarka, P., & Voisin, G. 2020, *A&A*, **644**, A145

Neubauer, F. M. 1980, *J. Geophys. Res.*, **85**, 1171

Nitu, I. C., Keith, M. J., Stappers, B. W., Lyne, A. G., & Mickaliger, M. B. 2022, *MNRAS*, **512**, 2446

Palenzuela, C., Lehner, L., Ponce, M., et al. 2013, *Phys. Rev. Lett.*, **111**, 061105

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

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

Philippov, A. A., Cerutti, B., Tchekhovskoy, A., & Spitkovsky, A. 2015a, ApJ, **815**, L19

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

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

Philippov, A., Timokhin, A., & Spitkovsky, A. 2020, *Phys. Rev. Lett.*, **124**, 245101

Pétri, J. 2012, *MNRAS*, **424**, 605

Pétri, J. 2016, *J. Plasma Phys.*, **82**, 635820502

Rea, N., Coti Zelati, F., Dehman, C., et al. 2022, ApJ, **940**, 72

Rea, N., Hurley-Walker, N., Pardo-Araujo, C., et al. 2024, ApJ, **961**, 214

Ruderman, M. A., & Sutherland, P. G. 1975, ApJ, **196**, 51

Stappers, B. W., Bailes, M., Lyne, A. G., et al. 1996, ApJ, **465**, L119

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

Tamburini, M., Pegoraro, F., Di Piazza, A., Keitel, C. H., & Macchi, A. 2010, *New J. Phys.*, **12**, 123005

Timokhin, A. N., & Arons, J. 2013, *MNRAS*, **429**, 20

Torres, R., Grismayer, T., Cruz, F., Fonseca, R., & Silva, L. 2024, *New Astron.*, **112**, 102261

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

Yee, K. 1966, *IEEE Trans. Anten. Propag.*, **14**, 302

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

Zhang, H., Sironi, L., & Giannios, D. 2021, ApJ, **922**, 261