© ESO 2019

1. Introduction

On August 2017, the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo detectors observed a gravitational wave signal that coincided with the measurement of a short-gamma-ray burst (SGRB) by the Fermi/Gamma-ray Burst Monitor (GBM) instrument (Abbott et al. 2017a). This unique joint detection allowed unambiguous association of SGRBs with binary neutron star mergers (Abbott et al. 2017b). It is expected on theoretical grounds that these binary systems release a large amount of energy prior to the merger. First developed by Goldreich & Lynden-Bell to explain decametric emission from the Io-Jupiter system (Goldreich & Lynden-Bell 1969), the Direct Current (DC) model was recently adapted to other types of astrophysical systems, including binary neutron stars (Piro 2012; Lai 2012). In this framework, the energy dissipation rate can go up to ∼10⁴⁴ erg s⁻¹ during the late stage of the inspiral, in the case of neutron stars with high magnetic fields (∼10¹³ G). Some of this energy is then converted into observable radiation; it may also partly be consumed by the creation of a plasma “fireball” (Vietri 1996; Hansen & Lyutikov 2001; Metzger & Zivancev 2016). The overall dissipated power predicted by these models is ∼10^{44 − 46} erg s⁻¹. In this case, the precursor should be observable from Earth.

The detection of an electromagnetic precursor would be very helpful to trigger multi-wavelength observations of upcoming mergers. Aside from this observational motivation, the interaction of plasma-filled magnetospheres has so far been poorly studied. In particular, it is not known whether the presence of another pulsar can enhance particle acceleration and how it affects pulsar spindown. To investigate a possible electromagnetic precursor to the merger, it is necessary to understand how particles can be accelerated, but the mechanisms at play are too complex to be analytically solved. Magnetohydrodynamic simulations were recently performed to model binary neutron star coalescences (Palenzuela et al. 2013; Dionysopoulou et al. 2015). Although this approach is appropriate to model the overall plasma structure, it fails to capture particle acceleration and cannot produce the shape of the electromagnetic outburst.

In recent years, kinetic modeling of pulsar magnetospheres using particle-in-cell (PIC) simulations has lead to significant progress in our understanding of magnetic dissipation, pair production, and particle acceleration (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; Cerutti et al. 2015; Philippov et al. 2015a,b; Belyaev 2015; Cerutti & Philippov 2017; Brambilla et al. 2018). In particular, these models proved fit to account for the two-peaked structure of gamma-ray pulses (Cerutti et al. 2016a; Philippov & Spitkovsky 2018; Kalapotharakos et al. 2018), and optical polarization properties of the Crab pulses (Cerutti et al. 2016b). Spherical simulations were able to reproduce a quasi force-free magnetosphere and identify magnetic reconnection in the equatorial current sheet as the main acceleration mechanism (Cerutti et al. 2015).

In this work we perform PIC simulations of a merging binary neutron star system in a 2D axisymmetrical cylindrical setup. To retain cylindrical symmetry, we keep both pulsars aligned with the symmetry axis, thus neglecting their relative orbital motion. It is likely that pulsar obliquity decreases with age, as can be seen on theoretical (Philippov et al. 2014) and observational (Tauris & Manchester 1998; Young et al. 2010) grounds (see, however, Beskin et al. 1993). Therefore, it is relevant to study pulsars with aligned magnetic moments and rotation axes. Although the problem is intrinsically 3D, performing 2D simulations allows us to obtain a much better resolution by decreasing the numerical cost, and to more easily gain physical insight into this complex problem. We explore the mechanisms of particle acceleration resulting from magnetospheric interactions in several configurations. Furthermore, we compute the high-energy radiation light curve that would be produced during the inspiral phase of the merger, in order to assess the observability of the precursor. In Sect. 2, we describe the numerical methods and the particular setup we simulate. We also report how the consistency of this new code was tested. The results are presented in Sect. 3.

2. Numerical methods

2.1. Particle-in-cell technique

In this study, we use the relativistic PIC code Zeltron (Cerutti et al. 2013, 2015). In PIC methods, the plasma is treated as a set of charged macro-particles. During one time step three operations must be performed. First, the positions and velocities of the particles are evolved using the Boris push (Birdsall & Langdon 2005). The electromagnetic source terms are then collected on a Yee grid (Yee 1966), using an area-weighting procedure. Finally, the fields are updated by solving discretized versions of Maxwell-Ampère’s and Maxwell-Faraday’s equations on the grid. Given the spatial resolution, the time step is enforced by the Courant-Friedrich-Lewy condition (Sironi & Cerutti 2017) so that the numerical scheme remains stable. This scheme does not enforce Maxwell-Gauss’ equation, so the total electric charge deposited on the grid is not conserved. Small errors can accumulate and give rise to nonphysical charge densities. The electric field is therefore corrected every 25 time steps, by solving Poisson’s equation: Δ(δΦ) = − (4πρ − ∇ ⋅ E), where ρ is the charge density, E₀ the uncorrected field, and E = E₀ − ∇(δΦ) the corrected electric field.

Since we intend to model a binary pulsar, the system no longer retains spherical symmetry, which is a natural consideration when studying isolated pulsars. Instead we have implemented a 2D axisymmetric cylindrical version of the code: particles move in 3D space but all quantities are invariant under rotations about the symmetry axis (Oz). Maxwell’s equations are then solved in cylindrical coordinates (R, z), where R is the distance from the symmetry axis. The simulated domain is defined as (R ∈ [R_min, R_max], z ∈ [ − R_max, R_max]), and is centered at the origin. All simulation parameters are given in Table 1. The simulated pulsars must have aligned magnetic moments and spin axes. This amounts to neglecting the orbital motion of the two pulsars, which are constantly aligned (see Fig. 1).

Fig. 1. Numerical setup. The initial magnetic field lines (black solid lines) are dipolar, the magnetic moments are here identical in direction and norm. The angular velocity vectors of the two stars Ω1 and Ω2 are either parallel or anti-parallel, and have the same magnitude. The particles are launched from the shell R⋆ ≤ ri ≤ R⋆ + δr with corotation velocity. The size of the injecting shell is exaggerated in the figure. Positrons and electrons are denoted by blue and orange disks. The dashed line marks the radius of the common light cylinder RLC = c/Ω. The striped zone has a resistivity profile λ(r) (see Eq. (2)). The pair-creation process is pictured on the right, a new pair being created at the expense of the energy of the primary particle.

2.2. Particles

The numerical setup is presented in Fig. 1. The stellar radius is denoted by R_⋆. The angular velocity vector and magnetic moment of the top (resp. the bottom) pulsar will be denoted Ω₁ and ℳ₁ (resp. Ω₂ and ℳ₂). They are aligned with the symmetry axis (Oz). The simulation starts in complete vacuum, and the stars remain empty of particles at all times. In order to fill the magnetospheres with plasma, electron/positron (e^±) pairs are uniformly injected within thin injection shells between R_⋆ and R₀ = R_⋆ + δr surrounding the stars. Charges with the right sign are then extracted from the surface, while the others rush back to the star. The particle distribution in this shell is reinitialized at every time step, in order to avoid numerical over-densities. These particles are initially corotating with the star, with toroidal velocity v_φ, i = RΩ_i. Particles that leave the simulation domain are removed.

The density of the surface plasma is set at the Goldreich-Julian density n_GJ, i = B_iΩ_i/(2πce), where B_i is the polar field strength at the stellar surface (Goldreich & Julian 1969). This is consistent with the actual charge density that is extracted from the star to screen the surface electric fields. The whole electromagnetic cascade, filling the magnetosphere with e^± plasma, is implemented in a simplified way. Following Philippov et al. (2015a), a critical Lorentz factor γ_c is empirically defined as a fraction ζ of the maximal energy available to a particle, which is given by the total voltage drop $\delta V_{\text{max},i}=R_0^2{B}_i{\Omega }_i/c$ from the pole to the equator at the stellar surface. If a particle reaches a Lorentz factor¹γ = γ_c, i = ζeδV_max, i/mc², and provided the local pair plasma multiplicity (defined as the ratio between the plasma density and the Goldreich-Julian density) remains smaller than ∼30, an e^± pair is created at the location of the parent particle, whose energy is consequently reduced. The resulting pair carries a fixed fraction of the energy of the parent particle. The ζ factor was chosen so as to sufficiently fill the magnetosphere and reach a force-free regime (Cerutti & Beloborodov 2017). This was checked by comparing the spindown of a single pulsar to the expected values in the monopole and dipole cases. Some numerical problems arise because of the discretization of the pulsar sphere by Cartesian cells in the (R, z) plane. For example, it can cause important angular inhomogeneities in the plasma injection, which impact the whole magnetosphere. This difficulty is circumvented by moderately increasing the thickness of the injecting shell and the number of macro-particles per cell.

2.3. Electromagnetic fields

In most simulations (see Sect. 3.5 for asymmetric simulations with Ω₁ ≠ Ω₂ and ℳ₁ ≠ ℳ₂), Ω₁ and Ω₂ have the same magnitude, so both pulsars have the same light cylinder R_LC = c/Ω. The two pulsars are initially separated by a distance a₀, and disposed symmetrically with respect to the origin. The initial magnetic field B_tot is the sum of two dipole fields generated by each pulsar. As the simulation proceeds, this field is enforced inside the stars (for r_i ≤ R₀, with r_i the distance from the center of pulsar “i”). Pulsar spin is implemented by imposing a poloidal electric field within the pulsar (for r_i ≤ R₀), given by

$$\begin{array}{c}\hfill {\mathit{E}}_i(\mathit{r})=-{\displaystyle \frac{1}{c}}({\mathbf{\Omega }}_i\times \mathit{r})\times {\mathit{B}}_{\mathrm{tot}},\end{array}$$(1)

arising in the laboratory frame because the electric field is zero in the corotating frame, assuming that the neutron star interior is perfectly conducting. The toroidal electric field in the star is set to zero. In order to mimic an open outer boundary, a numerical resistivity λ(r) (where r is the distance from the origin) is given to the medium for r ≥ R_pml, with R_pml = 0.9 R_max. This amounts to adding a resistive term to Maxwell’s equations: for instance Maxwell-Faraday’s equation becomes ∂B/∂t = −λ(r)B − c ∇ × E. In this part of the simulated domain the fields are critically damped. The chosen resistivity profile is selected to be smooth, so as to avoid unwanted reflections (Cerutti et al. 2015):

$$\begin{array}{c}\hfill \lambda (r)\propto {\left({\displaystyle \frac{r-R_{\mathrm{pml}}}{R_{\max }-R_{\mathrm{pml}}}}\right)}^3,\text{for}r\ge R_{\mathrm{pml}}.\end{array}$$(2)

On the outer boundaries, a zero-gradient condition is enforced. Cylindrical symmetry requires the radial (B_R, E_R) and toroidal (B_φ, E_φ) fields to be zero on the symmetry axis, whereas the vertical (B_z, E_z) components have no gradient on the axis. The correction potential δΦ is set to zero inside the pulsars, with an additional zero-gradient condition for the electric field at the stellar surfaces.

Initially, the electric field is zero everywhere but in the pulsar. A simulation starts with the launching of a spherical wave (E_θ, B_φ) distributing the electric field in the domain. If the toroidal magnetic field B_φ is set to zero for r ≤ R₀, the simulation cannot start properly. On the other hand, if B_φ is simply left free to evolve, some spurious B_φ grows at a slow but constant rate inside the injection shell. This is due to a round-off error in the computation of the φ component of ∇ × E in the shell (R_⋆ ≤ r_i ≤ R₀), which is supposed to be zero inside the star. This problem does not occur in spherical coordinates, because it is sufficient to use a single radial cell as an injection domain. In this cylindrical setup the shell must be several cells wide to enable homogeneous plasma injection. To alleviate this difficulty, a resistivity profile similar to Eq. (2) was given to the injection shell, where B_φ was left free to evolve. The induced damping rate was chosen to be slightly greater than the spurious growth rate.

2.4. Pulsar separation

Since the orbiting binary pulsar loses energy because of gravitational wave emission, the two stars get closer and closer. The orbital motion of the stars was not simulated, but the analytical expression for the separation a(t) between them was implemented ad hoc. For a pair of binary stars with masses m₁, m₂, separated by a distance a, the rate of energy loss is given by (Carroll 2003):

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=-{\displaystyle \frac{32G^4(m_1+m_2)m_1^2{m}_2^2}{5c^5{a}^5}}·\end{array}$$(3)

Injecting the gravitational energy of a two-body system E = −Gm₁m₂/2a and solving the associated differential equation yields

$$\begin{array}{c}\hfill a(t)=a_0{(1-{\displaystyle \frac{t}{\tau }})}^{1/4},\end{array}$$(4)

where a₀ is the initial separation. Numerical evaluation for the inspiral time τ yields

$$\begin{array}{c}\hfill \tau =0.003{\left({\displaystyle \frac{a_0}{30}}\right)}^4\text{s}.\end{array}$$(5)

The inspiral time was rescaled so as to match the code units, with τ ∼ P(a₀/3R_⋆)⁴.

The inspiral begins in the simulation only once both magnetospheres have been filled with plasma and have reached a quasi-steady state (the topology of the magnetosphere and the field strength no longer vary), which usually requires about 2.5 rotation periods. Once this phase has elapsed, the pulsar separation decreases according to Eq. (4). The simulation ends when the pulsars touch. Otherwise, no general relativity effect was included in the simulation.

2.5. Modeling radiation

In order to evaluate the electromagnetic signature of the merger, the radiation of accelerated particles must be computed. The procedure is described in Cerutti et al. (2016a). Particles cool through synchrotron radiation. The radiated power spectrum emitted by a single particle with Lorentz factor γ moving in a perpendicular magnetic field B_⊥ in the ultra-relativistic limit γ ≫ 1 is (Blumenthal & Gould 1970)

$$\begin{array}{c}\hfill F_{\nu }(\nu )={\displaystyle \frac{\sqrt{3}{e}^3{B}_{\perp }}{m{c}^2}}{\displaystyle \frac{\nu }{{\nu }_c}}{\displaystyle {\int }_{\nu /{\nu }_c}^{+\infty }{K}_{5/3}(x)\mathrm{d}x,}\end{array}$$(6)

where K_5/3 is the modified Bessel function of order 5/3. The critical frequency ν_c is given by ν_c = 3eB_⊥γ²/4πmc. In the presence of an electric field, a standard procedure consists in replacing the perpendicular magnetic field B_⊥ by an effective field strength (Kelner et al. 2015; Cerutti et al. 2016a):

$$\begin{array}{c}\hfill {\stackrel{\sim }{B}}_{\perp }=\sqrt{{(\mathit{E}+(\mathit{\nu }/c)\times \mathit{B})}^2-{(\mathit{\nu }·\mathit{E}/c)}^2}.\end{array}$$(7)

This procedure is valid only in the ultra-relativistic limit, with ∥ν∥ ≈ c, and takes into account both synchrotron and curvature radiation. This procedure allows one to deal with the E × B drift. Indeed, in the wind zone, the motion of the plasma is dominated by this drift. Injecting the drift velocity U = cE × B/B² into Eq. (7) yields ${\stackrel{\sim }{B}}_{\perp }\propto \mathit{E}·\mathit{B}\approx 0$. Physically, the plasma does not radiate significantly when the flow is dominated by the E × B drift, because in the drift frame where the electric field is zero the magnetic field strength is reduced by a factor Γ ≫ 1, the Lorentz factor of the bulk flow. Keeping the actual B_⊥ instead of the effective field would lead to highly overestimated synchrotron losses.

Moving charged particles emit “macro-photons” along their direction of motion. This is only valid in the ultra-relativistic limit, as the emission cone has a semi-aperture angle 1/γ ≪ 1. Once emitted, these photons escape the simulation (unless they hit a star). Each macro-photon conveys a power spectrum given by Eq. (6). The output of the code is νF_ν(ν, R, z), where F_ν = hν dN/dν is the locally emitted power spectrum. In order to find the high-energy emission sites, the code outputs ${\displaystyle {\int }_{{\nu }_0}^{+\infty }\nu }{F}_{\nu }\text{d}\ln \nu$, where ν₀ = qB₀/mc is a typical synchrotron frequency. In this 2D axisymmetric setup the spherical surface of a pulsar must be discretized in Cartesian coordinates (see Sects. 2.2 and 2.3). Right at the stellar surface, fields suffer discontinuities that produce spurious particle acceleration. Consequently, we choose to discard macro-photons emitted within 5% of the stellar surfaces. Similarly, particles that would undergo nonphysically high accelerations near the surfaces are also removed from the simulation.

The emitted macro-photons are collected at infinity, making an angle α_obs with the symmetry axis. The propagation time of the photons depends on the location of its emitting particle and its direction of motion. The time delay δt between a photon emitted at polar coordinates (ρ, θ) in the (R, z) plane, propagating with an angle α with respect to the symmetry axis, and one emitted at the origin in the same direction, is given by

$$\begin{array}{c}\hfill \delta t={\displaystyle \frac{\rho }{c}}cos(\alpha -\theta ).\end{array}$$(8)

All photons emitted at a given time step are distributed on a grid according to their time delay. Since the merger is inherently unsteady, we collect photons every ten time steps.

2.6. Consistency checks

Since PIC simulations of pulsar magnetospheres had not yet been conducted in cylindrical coordinates, we compared the output of the code to previous results. Some analytical solutions are known for an isolated pulsar magnetosphere. First, the electromagnetic solver alone was tested in vacuum. For a monopolar field B = B₀(R₀/r)²e_r, the electric field can be readily computed in spherical coordinates (r, θ, φ) in the whole space (Cerutti & Beloborodov 2017; Jackson 1998), using the boundary condition at r = R₀ given by Eq. (1):

$$\begin{array}{c}\hfill \{\begin{array}{c}\hfill E_{\mathrm{r}}=-{\displaystyle \frac{2B_0{R}_0^4}{R_{\mathrm{LC}}}}{\displaystyle \frac{cos\theta }{r^3}},\hfill \\ \hfill E_{\theta }=-{\displaystyle \frac{B_0{R}_0^4}{R_{\mathrm{LC}}}}{\displaystyle \frac{sin\theta }{r^3}}·\hfill \end{array}\end{array}$$(9)

The simulated electric field of an empty magnetosphere was found to be consistent with Eq. (9), the angular and radial dependence agreeing with this theoretical prediction. The toroidal field remained zero in the steady state. The code was then tested on a monopole force-free magnetosphere, for which an analytical solution also exists (Michel 1973a,b; Henriksen & Rayburn 1971). In the presence of plasma a toroidal field develops, given by

$$\begin{array}{c}\hfill B_{\phi }=-B_0{\displaystyle \frac{R_0}{R_{\mathrm{LC}}}}{\displaystyle \frac{R_0}{r}}sin\theta .\end{array}$$(10)

We recovered this behavior with good accuracy. Finally, the single force-free aligned dipole configuration was also simulated. The spindown power of the dipole pulsar L₀ can be analytically estimated (Cerutti & Beloborodov 2017) to scale like

$$\begin{array}{c}\hfill L_0={\displaystyle \frac{c{B}_0^2{R}_{\star }^6}{R_{\mathrm{LC}}^2}},\end{array}$$(11)

where B₀ is the equatorial surface field of the pulsar. This served as another quantitative test for the 2D axisymmetric code. Previous 2D spherical simulations (Cerutti et al. 2015; Belyaev 2015) showed that the outgoing Poynting flux indeed scales like L₀, provided the force-free regime is achieved. If the pair creation parameters are properly adjusted, so as to reach a force-free regime, we confirm this result in the axisymmetric setup. The quantity L₀ provides a useful reference to convert energy flux from code units to physical units.

3. Results

Aside from Sect. 3.5, both pulsars have equal magnetic field strengths and rotation periods in the following. We first present simulations that have reached a quasi-steady state, in order to study the magnetospheric interactions at play. In Sects. 3.4 and 3.5 we then study the inspiral phase and the resulting electromagnetic outburst.

3.1. Magnetospheric structure

Eight configurations can be simulated with this axisymmetric setup, depending on the relative orientations of the magnetic moments and the rotation axes. We focus on the four configurations where the pulsars have their magnetic moments parallel or anti-parallel, and their rotation axes parallel or anti-parallel. The top pulsar has ℳ = ℳ ⋅ e_z >  0 and Ω >  0, whereas the bottom pulsar can have its magnetic moment or rotation axis reversed. Until Sect. 3.5 we only study the two configurations with the magnetic moments aligned for both stars. Indeed, symmetric simulations with opposed magnetic moments are qualitatively similar, regardless of the relative spin orientation (this is not shown in the symmetric case, but see Sect. 3.5). The two opposed dipolar fields exclude one another, so no field line links the pulsars, and no reconnection layer occurs in between.

3.1.1. Parallel configuration

Figure 2a shows the steady-state toroidal magnetic field. Since the two pulsars have no relative motion before the merger is initiated, a large zone between them remains closed, with no outflowing currents and no toroidal fields. The main feature of this configuration is the equatorial sheet, which supports a discontinuity in the toroidal (like in the isolated pulsar case) and poloidal field. A discontinuity in the toroidal field induces a radial current sheet, visible in Fig. 2b. We note that this “inter-pulsar” current sheet qualitatively differs from the “proper” current sheets; it carries a current with the opposite sign, which means electrons undergo greater acceleration than with a single pulsar. It is a third promising site for particle acceleration, apart from the two proper current sheets. The location of the Y-point of this current sheet is determined exclusively by the separation between the pulsars, rather than the light cylinder. This means that unlike the isolated pulsar configuration, here reconnection occurs well within the light cylinder where the fields are much stronger, although this is more obvious at closer separations. The field lines are represented in three dimensions in Fig. 4.

Fig. 2. Panel a: normalized toroidal magnetic field Bφ/B0 in the parallel configuration. The poloidal magnetic field lines are depicted as black solid lines. A large magnetic island is visible in the inter-pulsar current sheet on the right; it can be recognized by its closed magnetic field lines. Panel b: normalized radial current density JR/cρGJ in the parallel configuration. The stellar separation is fixed at a = 3RLC.

Fig. 3. Panel a: normalized toroidal magnetic field Bφ/B0 in the anti-parallel configuration. The poloidal magnetic field lines are depicted as black solid lines. Panel b: normalized vertical current density Jz/cρGJ in the anti-parallel configuration. The stellar separation is fixed at a = 3RLC.

Fig. 4. Magnetic field lines in the parallel configuration. The stellar separation is fixed at a = 3 RLC. The pulsars are depicted as two black spheres.

Densities for the two species are shown in Fig. 6, as well as the bulk velocity field lines. The small discrepancies between the electronic and positronic densities imply that the plasma is globally quasi-neutral, except in the closed zones and within the current sheets. Some density gaps are visible surrounding the current sheets, as in Philippov et al. (2015a). In the wind zone, the solution resembles Michel’s split monopole (Michel 1973a). The bulk velocity field is similar for both species in the wind zone, where it is approximately radial, as prescribed by the E × B drift velocity of the particles. A difference can be found in the inter-pulsar current sheet for instance: electrons flow outward whereas positrons move back towards the pulsars.

3.1.2. Anti-parallel configuration

The initial magnetic field is the same as in the parallel case, but here the two stars spin in opposite directions. Consequently, the field lines linking the two stars are twisted (see Fig. 5), while a strong toroidal field develops in between, as can be seen in Fig. 3a. We can relate this configuration to the DC model (Lai 2012; Piro 2012) which considers a magnetic neutron star orbiting a non-magnetic but perfectly conducting companion. The motion of the non-magnetic star with respect to the magnetic field of the primary induces an electromotive force, which in turn drives currents along the magnetic field lines emerging from the magnetic star. The energy is dissipated because of the plasma and stellar resistivity. This model allows us to shed light on this magnetic configuration: the toroidal field between the pulsars induces a strong current flowing from the down pulsar to the top, as can be seen in Fig. 3b. No instability prevents B_φ from growing in these 2D simulations.

Fig. 5. Magnetic field lines in the anti-parallel configuration. The stellar separation is fixed at a = 3 RLC. The pulsars are depicted as two black spheres.

Fig. 6. Logarithm of the electronic (panel a) and positronic (panel b) densities, normalized by the Goldreich-Julian density and in the parallel configuration. The bulk velocity field lines of the electrons (on the left) and positrons (on the right) are displayed in white. The stellar separation is fixed at a = 3 RLC.

In both configurations, the single pulsar current sheets seem to be repelled from the midway equator, although they both carry opposite electric charges. This is a magnetic pressure effect. The presence of another pulsar generates a stronger toroidal field outside the light cylinder, which induces a magnetic pressure acting perpendicular to the field lines. The equilibrium position of the current sheets results from a balance between magnetic and kinetic plasma pressure. Close to the outer boundary, the current sheets become unstable due to the kink instability. It is also subject to the tearing instability. Magnetic islands are generated at the Y-point before being advected away (see Fig. 2a).

The latitudinal stripes visible in Figs. 2b and 3b likely stem from the discretization of the sphere, as mentioned in Sect. 2.2, and should not be considered physical. However, the radial current oscillations have a physical basis. If some unscreened parallel electric field lies at the stellar surface, particles are accelerated and induce pair creation. Then the local plasma density rises, so the electric field is momentarily screened and particle acceleration is quenched. As the plasma flows away, the same process starts again, giving rise to these oscillations.

3.2. Spindown

The interaction between the two pulsars significantly impacts their spindown. Two steady-state simulations were performed, with a = 3 R_LC and a = 1.25 R_LC. We compute the Poynting energy flux flowing out of each pulsar at r_i = 1.1 R₀ separately, then we sum the results to get the total luminosity L_⋆. This is compared to the Poynting flux flowing out of the simulation L_out, through a sphere of radius R_pml centered on the origin. The dissipation rate is then computed as ϵ = (L_⋆ − L_out)/L_⋆. The results are shown in Table 2. The spindown values lie between 2 L₀ and 4 L₀, L₀ being the spindown of an isolated pulsar Eq. (11). If the separation a tends to infinity, the spindown will be 2L₀. If a tends to 0, the magnetic moment of the star will double and its spindown will quadruple. For intermediate separations, one pulsar can intercept outgoing field lines from the other, thus reducing its Poynting flux.

At fixed separation, the anti-parallel total spindown is higher than the parallel spindown. The spindown does not increase when a decreases in the parallel case because the presence of another pulsar does not significantly affect the magnitude of the toroidal field (since they have opposite polarities), and the Poynting flux approximately scales like $B_{\phi }^2$. More field lines outgoing from one pulsar are captured by the other, diminishing the amount of energy extracted along the open field lines. On the other hand, in the anti-parallel configuration the toroidal field is amplified near the stars, which leads to a higher Poynting flux.

The fraction of dissipated spindown increases as the separation decreases. It is greater in the parallel configuration. This fact seems to support the idea that the inter-pulsar current sheet is a prominent site for dissipation and particle acceleration. The absence of this layer in the anti-parallel configuration could imply that the DC mechanism is a less efficient dissipation mechanism than magnetic reconnection. However, the amplification of the toroidal field in this configuration compensates for this effect.

3.3. Particle acceleration and high-energy radiation

The sites of pair creation can be investigated through the local quantity α(θ) = J_r(R₀, θ)/cρ_GJ, where J_r is the (spherically) radial current flowing out of a pulsar, θ is the polar angle defined with respect to this pulsar, and ρ_GJ = en_GJ is the Goldreich-Julian charge density. If 1/α <  1, the charge-separated current escaping from the star fails to sustain the current required by the current sheet, so a parallel electric field develops which accelerates particles and produces pairs (Beloborodov 2008). Pair creation must occur to supply current. This diagnostic was computed for both the up and down pulsars, following Parfrey et al. (2012).

The anti-parallel configuration is presented in Fig. 7a. The presence of the second pulsar greatly increases α, which can reach values from 3 to 4 at the south pole of the top pulsar and at the north pole of the bottom pulsar (the “inward” poles). This is consistent with the DC model: the electromotive force between the stars drives a large current flowing from the “down” to the “up” pulsar. Pair creation is extremely efficient in between the pulsars, hence the formation of a dense pair plasma. However, this plasma remains highly magnetized and does not radiate much. At closer separation, the locations of the peaks in α are almost unchanged, whereas the magnitude of the inward pole peaks rises a lot. Figure 7b shows the same diagnostic for the parallel configuration. Here, the magnitudes of the inward and outward pole peaks are similar. At closer separation, the inward pole peaks move towards the stellar equators.

Fig. 7. Plots of α(θ) = Jr(R0, θ)/cρGJ in the parallel (panel a) and anti-parallel (panel b) configurations, for the top (top panel) and bottom (bottom panel) pulsars. The stellar distance is fixed at a = 3 RLC. The dashed lines mark the locations of θpc and π − θpc (see Eq. (12)). The grey-shaded zones indicate the inward pole regions of emission, and the color-shaded zones mark the outward pole emission.

In both cases, pair creation seems to be slightly enhanced at the outward poles with respect to the isolated pulsar, α remaining close to 1. Indeed, this is why periodic oscillations in current density can occur at the poles in Figs. 2b and 3b. We also find that the closed zone gets smaller and the polar cap larger. This is proven by the fact that the return current sheets (α <  0 marks the location of the separatrix, separating closed field lines from open ones) are located closer to the equator and further from the theoretical polar cap angles θ_pc and π − θ_pc, where θ_pc is given by

$$\begin{array}{c}\hfill {sin}^2{\theta }_{\mathrm{pc}}\approx {\displaystyle \frac{R_0}{R_{\mathrm{LC}}}}·\end{array}$$(12)

This likely results from the magnetic pressure exerted by the other pulsar. This effect is more significant in asymmetric simulations (see Sect. 3.5).

The high-energy radiation maps for electrons and positrons are shown in Figs. 8a and 8b for the parallel configuration, and Figs. 8c and 8d for the anti-parallel configuration. As mentioned above, particles too close to the stellar surface are spurious and should not contribute to the radiation diagnostic. The anti-parallel maps show that at large separation, high-energy radiation is mainly due to the current sheets belonging to each pulsar. Furthermore, in accordance with the diagnostic of Fig. 7a, particle acceleration is enhanced at the inward poles with respect to the outward polar cap emission. In the parallel configuration, the inter-pulsar current sheet radiates predominantly, as expected since the Y-point lies within the light cylinder where the field is stronger. From the radiation data, one can infer the fraction of the spindown that was channeled to electromagnetic radiation by computing the total radiative efficiency in the total simulation domain 𝒱 (excluding a thin shell of width 0.05 R₀ at the stellar surface):

$$\begin{array}{c}\hfill {\eta }_{\mathrm{rad}}={\displaystyle \frac{1}{L_0}}{\displaystyle {\int }_0^{+\infty }dln\nu {\int }_{\mathcal{V}}\mathrm{d}V\nu (F_{\nu ,\mathrm{el}}+F_{\nu ,\mathrm{ions}}).}\end{array}$$(13)

Fig. 8. High-energy radiation maps, as defined in Sect. 2.5, for the parallel configuration (electrons: panel a, and positrons: panel b) and anti-parallel configuration (electrons: panel c, and positrons: panel d).

In this study, the radiative efficiency of an isolated pulsar is η_rad = 2.7%. The emitted radiation in the steady state two-pulsar setup is much more intense than in the isolated pulsar case. At a = 3 R_LC, we find η_∥ = 22% in the parallel configuration and η_∦ = 23% in the anti-parallel configuration. At the closer separation a = 1.25 R_LC, and the efficiencies are η_∥ = 29% and η_∦ = 81%, respectively.

To sum up, even though the reconnection mechanism is more efficient at dissipating electromagnetic energy than the DC mechanism, the amplification of B_φ in the anti-parallel configuration compensates for its relative inefficiency. The amount of dissipated energy is similar in both cases (around 1.0 L₀ at a = 1.25 R_LC, inferred from Table 2). Eventually, the radiative efficiency increases faster with decreasing separation in the anti-parallel configuration. This can be explained by the low densities achieved in the inter-pulsar reconnection layer, with respect to the dense zone between the stars. Although particles experience greater accelerations in the parallel configuration, this affects a smaller number of particles, so that more energy is converted into particle kinetic energy in the anti-parallel configuration.

3.4. The inspiral phase

The light curves obtained using the procedure described in Sect. 2.5 are presented in Fig. 9. In these simulations, the computation of the light curves started after 1.25 P (marked by a dash-dotted line in Fig. 9), while the inspiral was initialized at 2.5 P (marked by a dashed line). Computing the light curve after a steady state is reached provides a baseline to evaluate the increase in the measured flux, and makes it possible to check that the routine works correctly by comparing the steady-state value to the radiative efficiency.

Fig. 9. Light curves in the parallel (black circles) and anti-parallel (blue triangles) configurations, normalized by the isolated pulsar spindown power L0. The diagnostic starts at the dash-dotted line. The merger starts at the dashed line, until the two stars touch at the dotted line. Photons are still collected by the observer until the solid line is reached. The time is normalized by the spin period of the pulsar.

The two curves are qualitatively different before the inspiral starts. In the parallel simulation, a steady state is reached when the merger is initiated. Conversely, the anti-parallel light curve has not yet reached a steady state when we start approaching the stars. The received flux rises steadily. The regular increase in radiated flux is due to the accumulation of toroidal field between the pulsars. This can be confronted with the DC model mentioned earlier. Lai (2012) noticed that flux tubes twisted beyond B_φ/B_z ≳ 1 break down, because the magnetic pressure exerted by the toroidal field is too large. This occurs if the resistivity of the plasma is too low, and the currents too high. The axisymmetry of the tube is broken by unstable kink modes, which distort the tube. Nonlinear evolution of the tube disrupts it completely. Magnetic energy is then dissipated by reconnection. Thus a quasi-periodic circuit can be expected: after the flux tube breaks, reconnection between the inflated field lines restores the linkage between the two stars and the cycle repeats. However, this phenomenon can only be captured in 3D simulations. The staircase appearance of the anti-parallel light curve possibly results from a discharge phenomenon similar to what was explained in Sect. 3.1. Nonetheless, the presence of such a twist in 2D simulations is promising as to the strength of the outburst, since even more energy should be released by the breaking of the flux tube.

After the merger starts, both curves eventually meet. The peak of the light curve is reached after the end of the merger (marked by the dotted line in Fig. 9). The full width at half maximum of the peak is T_merge ≈ 0.2 P. This can be compared to the final orbital frequency in the GW170817 neutron star merger f_max ≈ 400 HZ (Abbott et al. 2017a,b). This means that the energy of the magnetospheric outburst is released during the last orbit of the binary. The angular distribution of the pulse is shown in Fig. 10. The radiation is mainly emitted in the equatorial plane in both configurations, although the signal is not strongly anisotropic. Consequently, the direction of the line of sight is not critical regarding observability. The anti-parallel angular distribution is broader than the parallel one, because the emission mechanism occurring before the merger is different. In contrast to magnetic reconnection, which mainly accelerates particles near the equatorial plane, in the anti-parallel configuration particles are accelerated when heading from one pulsar to the other, thus pointing towards high latitudes.

Fig. 10. Angular distribution of the light-curve signal during the inspiral phase in the two configurations. Time integration ran from the start of the merger until the end of the pulse. The dotted line marks the equator αobs = π/2.

The similarity between the two light curves in the late phase strongly indicates that close to the merger, the physics that accelerates particles and dissipates energy is the same. In both configurations, there is a toroidal current sheet resulting from the discontinuity in the poloidal magnetic field. At high separations its effect is negligible with respect to the other current sheets (the proper pulsar current sheets, or the inter-pulsar radial current sheet in the parallel configuration). However, as the pulsars close up, it becomes predominant in both cases because the poloidal field dominates over the toroidal component. The snapshots of J_φ/cρ_GJ right before the pulsars touch show that the toroidal sheet is indeed predominant, and that it has a similar shape in both configurations. At the end of the inspiral, dissipation mainly occurs in the toroidal current sheet through magnetic reconnection, which converts a great amount of energy, as the fields so close to the stars are huge. Plus, magnetic reconnection is forced as the pulsars close up, which increases the reconnection rate.

In general, the pulsar inspiral leads to a great increase in radiated power. The instantaneous radiative efficiency (displayed in Fig. 11) shows an increase in the bolometric luminosity by a factor of approximately ten for both configurations between the beginning and the end of the inspiral phase. This amounts to an increase by a factor of approximately 100 with respect to the case of two quasi-isolated pulsars, with a ≫ R_LC. Even so, the dissipated power is well below the previous theoretical expectations from the DC model. Including the orbital motion in 3D simulations would probably increase the simulated dissipated energy, as the relative motion and therefore the electromotive force would be greater. However, these theoretical estimates concern the total output energy, yet only the radiated energy will leave an observable signature.

Fig. 11. Total instantaneous radiative efficiency, ηrad, given by Eq. (13) and plotted as a function of time in both configurations. The beginning of the merger is marked by the dashed line. The time is normalized by the spin period of the pulsars.

3.5. Asymmetric simulations

So far we have focused on symmetric binaries. Actual binary pulsars are likely to be asymmetric (Tauris et al. 2017; Piro 2012), comprising a recycled pulsar spinning rapidly and a younger, slower pulsar with a stronger magnetic field. This was confirmed by the discovery of the double pulsar J0737–3039A/B (Kramer & Stairs 2008), for which Ω₁/Ω₂ ≈ 10² and B₁/B₂ ≈ 5 × 10⁻³. However, such ratios are far too extreme for our setup. The light cylinders of both pulsars must be included in the simulation domain, otherwise the closed magnetic field lines of the slow pulsar interact with the outer damping layer. In this case the closed corotating zone is artificially reduced, which induces spurious particle acceleration. Instead, we performed simulations with Ω₁/Ω₂ = 4 and B₁/B₂ = 0.25, and increased the size of the box to include the whole magnetosphere of the slow pulsar. The bottom pulsar is the old recycled pulsar; it has the same spin and magnetic strength as in the symmetric simulations, whereas the top pulsar has a stronger field and slower spin.

A map of the normalized toroidal field in the configuration with aligned spins and magnetic moments is shown in Fig. 12a. In contrast to the case of symmetric simulations, the relative orientations of the spins are of little consequence. The parallel and anti-parallel configurations are similar, since the toroidal field of the quickly recycled pulsar dominates anyway. Strongly twisted magnetic field lines connect the two pulsars, like in the symmetric anti-parallel setup. The closed zone of the weak pulsar is compressed because of the influence of the strong pulsar above it. The inter-pulsar current sheet present in the symmetric parallel configuration (see Fig. 2a) is missing from Fig. 12a. Snapshots of the toroidal field at early times show that these two current sheets merge as the toroidal field of the fast pulsar takes over. The resulting current sheet is deflected upwards in this parallel configuration, but downwards if they are anti-parallel. This is because the compensation of the top pulsar toroidal field makes the magnetic pressure initially smaller above the down pulsar than below.

Fig. 12. Normalized toroidal magnetic field Bφ/B0, down for asymmetric simulation with Ωdown/Ωup = 4 and B0, down/B0, up = 0.25, in the configuration with identical spins and aligned (panel a) or anti-aligned (panel b) magnetic moments. The simulation has reached a steady state. The young (resp. old) pulsar is at the top (resp. bottom). The poloidal magnetic field lines are depicted as black solid lines. The black dashed (resp. dotted) line denotes the light cylinder of the bottom (resp. top) pulsar. The stellar separation is fixed at a = 1.25 RLC.

In this situation it is interesting to study also the case of two opposed magnetic moments. As shown in Fig. 12b, the situation is quite different. No field line can connect the two stars, so the influence of the strong pulsar on the weak one can only be indirect. In addition, no current can flow from one pulsar to the other. As a result, the relative orientation of the spins has rigorously no impact on the magnetic topology; the polarities of the toroidal field and currents due to the down pulsar are simply reversed. At t = 0, the magnetosphere of the weak pulsar is confined within the strong field of the other pulsar. Its magnetic field lines then open up and extend to infinity as rotation is imposed. The current sheet of the weak pulsar is strongly deflected. Remarkably, at high separations, this configuration radiates no more than the isolated pulsar. However, after the merger starts, the radiated flux strongly increases by a factor of approximately 100, until it reaches usual outburst values (see Table 3).

We identify multiple dissipation mechanisms. Magnetic reconnection occurs at the Y-point of the current sheet, which is pushed closer to the star by the strong pulsar. In the case of aligned magnetic moments, strong pair creation and particle acceleration occur between the pulsars, as in the symmetric anti-parallel configuration. The higher outburst values in asymmetric simulations can be understood in this framework (see Table 3). Asymmetric simulations exhibit an intense magnetic field, which compresses the magnetosphere of the old pulsar. The Y-point at the foot of the proper current sheet of the recycled pulsar lies within its light cylinder. Reconnection is therefore more efficient than in symmetric simulations. We also think that in contrast to symmetric simulations, emission from the north polar cap of the strong pulsar may be essential. High-energy emission in this area has very contrasted stripes, which indicates a violent pair-creation cascade. The presence of the fast pulsar below implies strong currents must flow from the south polar cap of the strong pulsar. Therefore strong currents must also flow from its north polar cap, and pair creation is triggered to supply them. Then the huge surface parallel electric field accelerates particles. This last mechanism was negligible in symmetric simulations. However, the radiation diagnostics in the asymmetric simulations must be taken with care. Since the magnetic field of the top pulsar is increased in this setup, the plasma quantities are less resolved, and the output data are not fully reliable close to the pulsars. The exact acceleration processes are therefore harder to infer.

4. Discussion and observational prospects

We have performed global PIC simulations of a binary pulsar, with various magnetic moments and spins. Our strongest conclusion is that the total radiated power increases by two orders of magnitude during the whole inspiral. The bolometric luminosity can reach up to ten times the spindown power of an isolated pulsar. The light curve presents a peak of roughly 0.2 P in width, approximately corresponding to the last orbit of the double pulsar. The shape of the peak depends very little on the relative orientation of the spins or magnetic moments. The radiative power is concentrated within the equatorial regions but presents a weak anisotropy. We find that magnetic reconnection is a key ingredient in particle acceleration and the emission of high-energy synchrotron radiation in all the configurations explored here. In the symmetric configuration, dissipation occurs mainly because of an inter-pulsar current sheet forming at the interface between the two pulsar magnetospheres. As a result, electromagnetic energy is released even in the case where the two pulsars spin synchronously. Before the last stage of the merger, a significant amount of energy is also dissipated between the two stars in the symmetric anti-parallel configuration. Furthermore, pair creation at the outward poles is amplified with respect to the isolated pulsar configuration. In particular, in asymmetric simulations, pair creation in the vicinity of the slow magnetized pulsar is revived by the presence of the fast-spinning pulsar.

A natural extension of this work would be to perform 3D PIC simulations of a binary aligned pulsar, and more generally with arbitrary inclinations. This would have several benefits. First, we would be able to capture plasma instabilities that can only develop in 3D. For instance, the symmetric anti-parallel configuration displays a highly twisted flux tube, which should be unstable in a 3D setup. This could be another source of energy dissipation. Second, a 3D setup would allow us to simulate more realistic configurations, but more importantly, we would be able to take orbital motion into account. Since the orbital frequency peaks at the merger, it might lead to an increased electromotive force, stronger currents, and more dissipation. The luminosity of the precursor predicted here should therefore be seen as a lower limit. Nevertheless, for the first time we are able to self-consistently solve the magnetospheric interaction between two coalescing pulsars. This simplified setup allows us to elucidate the physical mechanisms at play.

Assuming a binary pulsar with a powerful Crab-like pulsar (Bühler & Blandford 2014), we predict a luminosity for the electromagnetic precursor of its coalescence of ∼10³⁸ erg s⁻¹. In the optimistic case where most of this energy is emitted in the Fermi/GBM range, this sets a distance upper limit of ≲4 Mpc for detection. Therefore, only the local group galaxies could be probed this way. For comparison, the neutron star merger event GRB 170817A, which occurred at a distance of 40 Mpc, released an amount of ∼10⁴⁶ erg s⁻¹ in the γ range (Abbott et al. 2017b). The situation is similarly grim in X-ray, but radio emission could be of better use. Some radio detectors such as MeerKAT, or the SKA array in the near future, map the whole sky in real time with high sensitivity. There are some caveats however: The radio luminosity of most pulsars is only a small faction of the spindown power (usually around 10⁻⁶). Dispersion by the interstellar medium is also likely to delay the arrival of the radio burst, which cannot arrive ahead of the merger. Besides, radio emission is coherent and does not directly originate from particle acceleration, but may rather be related to pair creation (Philippov et al. 2015b). Although we observe a strong pair-creation cascade at the polar caps, especially during the inspiral phase, we are not able to compute the radio efficiency of the binary. Still, radio detection is a promising candidate for the prospects of neutron star merger observability. In particular, fast radio bursts have a duration and energy range consistent with pulsar coalescence, and could be counterparts to nonrepeating events (Pen 2018).

Acknowledgments

We would like to thank the referee Richard Henriksen for valuable comments on the manuscript. We acknowledge the support from CNES and the Université Grenoble Alpes (IDEX-IRS grant). This work was granted access to the HPC resources of CINES on Occigen under the allocation A0030407669 made by GENCI.

References

All Tables

All Figures

Fig. 1. Numerical setup. The initial magnetic field lines (black solid lines) are dipolar, the magnetic moments are here identical in direction and norm. The angular velocity vectors of the two stars Ω1 and Ω2 are either parallel or anti-parallel, and have the same magnitude. The particles are launched from the shell R⋆ ≤ ri ≤ R⋆ + δr with corotation velocity. The size of the injecting shell is exaggerated in the figure. Positrons and electrons are denoted by blue and orange disks. The dashed line marks the radius of the common light cylinder RLC = c/Ω. The striped zone has a resistivity profile λ(r) (see Eq. (2)). The pair-creation process is pictured on the right, a new pair being created at the expense of the energy of the primary particle.

In the text

Fig. 2. Panel a: normalized toroidal magnetic field Bφ/B0 in the parallel configuration. The poloidal magnetic field lines are depicted as black solid lines. A large magnetic island is visible in the inter-pulsar current sheet on the right; it can be recognized by its closed magnetic field lines. Panel b: normalized radial current density JR/cρGJ in the parallel configuration. The stellar separation is fixed at a = 3RLC.

In the text

	Fig. 3. Panel a: normalized toroidal magnetic field Bφ/B0 in the anti-parallel configuration. The poloidal magnetic field lines are depicted as black solid lines. Panel b: normalized vertical current density Jz/cρGJ in the anti-parallel configuration. The stellar separation is fixed at a = 3RLC.
In the text

	Fig. 4. Magnetic field lines in the parallel configuration. The stellar separation is fixed at a = 3 RLC. The pulsars are depicted as two black spheres.
In the text

	Fig. 5. Magnetic field lines in the anti-parallel configuration. The stellar separation is fixed at a = 3 RLC. The pulsars are depicted as two black spheres.
In the text

	Fig. 6. Logarithm of the electronic (panel a) and positronic (panel b) densities, normalized by the Goldreich-Julian density and in the parallel configuration. The bulk velocity field lines of the electrons (on the left) and positrons (on the right) are displayed in white. The stellar separation is fixed at a = 3 RLC.
In the text

Fig. 7. Plots of α(θ) = Jr(R0, θ)/cρGJ in the parallel (panel a) and anti-parallel (panel b) configurations, for the top (top panel) and bottom (bottom panel) pulsars. The stellar distance is fixed at a = 3 RLC. The dashed lines mark the locations of θpc and π − θpc (see Eq. (12)). The grey-shaded zones indicate the inward pole regions of emission, and the color-shaded zones mark the outward pole emission.

In the text

	Fig. 8. High-energy radiation maps, as defined in Sect. 2.5, for the parallel configuration (electrons: panel a, and positrons: panel b) and anti-parallel configuration (electrons: panel c, and positrons: panel d).
In the text

Fig. 9. Light curves in the parallel (black circles) and anti-parallel (blue triangles) configurations, normalized by the isolated pulsar spindown power L0. The diagnostic starts at the dash-dotted line. The merger starts at the dashed line, until the two stars touch at the dotted line. Photons are still collected by the observer until the solid line is reached. The time is normalized by the spin period of the pulsar.

In the text

	Fig. 10. Angular distribution of the light-curve signal during the inspiral phase in the two configurations. Time integration ran from the start of the merger until the end of the pulse. The dotted line marks the equator αobs = π/2.
In the text

	Fig. 11. Total instantaneous radiative efficiency, ηrad, given by Eq. (13) and plotted as a function of time in both configurations. The beginning of the merger is marked by the dashed line. The time is normalized by the spin period of the pulsars.
In the text

Fig. 12. Normalized toroidal magnetic field Bφ/B0, down for asymmetric simulation with Ωdown/Ωup = 4 and B0, down/B0, up = 0.25, in the configuration with identical spins and aligned (panel a) or anti-aligned (panel b) magnetic moments. The simulation has reached a steady state. The young (resp. old) pulsar is at the top (resp. bottom). The poloidal magnetic field lines are depicted as black solid lines. The black dashed (resp. dotted) line denotes the light cylinder of the bottom (resp. top) pulsar. The stellar separation is fixed at a = 1.25 RLC.

In the text