INTRODUCTION

Merging compact binaries are prime targets for gravitational-wave observations. In the case where one of the binary constituents is a neutron star the merger can be accompanied by electromagnetic counterparts, including afterglows and gamma-ray bursts (Cowperthwaite et al. 2017;Chornock et al. 2017;Villar et al. 2017;Nicholl et al. 2017;Troja et al. 2018;Tanvir et al. 2017;Drout et al. 2017;Abbott et al. 2017;Savchenko et al. 2017;Troja et al. 2017;Margutti et al. 2017Margutti et al. , 2018;;Hajela et al. 2019;Hallinan et al. 2017;Alexander et al. 2017;Ghirlanda et al. 2019;Mooley et al. 2018a,b). In addition, the presence of strong magnetic fields prior to the collision could also give rise to yet unseen precursor transients (Hansen & Lyutikov 2001;Palenzuela et al. 2013b;Most & Philippov 2023a;Cooper et al. 2022), emost@caltech.edu kchatziioannou@caltech.edu which present exciting opportunities for next-generation observatories (Corsi et al. 2024). For example, shortduration gamma ray bursts are associated with neutron star mergers (Kumar & Zhang 2015;Gottlieb et al. 2023;Meszaros 2006) as confirmed by the GRB170817A burst coincident with a gravitational-wave event (Abbott et al. 2017). While the main burst is likely associated with the merger itself, ∼ 10% of the observed long and short bursts are accompanied by precursor emission (Troja et al. 2010;Coppin et al. 2020;Xiao et al. 2022;Dichiara et al. 2023;Coppin et al. 2020). Preceding the main burst by ∼ 100 s to a few seconds for long bursts (Coppin et al. 2020) and ∼ 2s for short bursts (Wang et al. 2020), including a discovery of quasiperiodicity in a burst precursor accompanying a kilonova afterglow and a long gamma-ray burst (Xiao et al. 2022), precursor emission is associated with the binary inspiral. One explanation links precursor emission to resonant excitation of modes in the neutron star crust (Tsang et al. 2012;Penner et al. 2012;Tsang 2013;Neill et al. 2022;Zhang et al. 2022) or in neutron star oceans in the inspiral (Sullivan et al. 2022(Sullivan et al. , 2023)).

In these models, binary tidal interactions excite oscillation modes in the crust-core interface of the neutron star(s) in a compact binary. These excitations cause the neutron star crust to fracture and release large amounts of energy, > 10 46-47 erg/s (Tsang et al. 2012;Penner et al. 2012). For example, torsional modes could be consistent with the 20 Hz quasi-periodic oscillations claimed in GRB211211A (Suvorov et al. 2022). Observation of such modes can put constraints on the equation of state of the crust, and in turn on the nuclear symmetry energy (Neill et al. 2023;Sotani et al. 2023;Neill et al. 2024) or the presence of deconfined quark matter in the core (Pereira et al. 2023). Such constraints would add to ground-based nuclear scattering experiments (Horowitz et al. 2014;Reed et al. 2021Reed et al. , 2024) ) and yield multi-messenger constraints on the equation of state complementary to those obtained through gravitational-wave (Raithel 2019;Chatziioannou 2020) or X-ray observations ( Özel & Freire 2016;Watts et al. 2016).

While estimates of the crustal energy budget and the tidal interactions are straightforward (Blaes et al. 1989;Baiko & Chugunov 2018) (see also Tsao et al. (2021); Sagert et al. (2023) for simulations of crustal oscillations), it remains unclear how waves launched from star quakes could convert into observable radiation. Nonradial shearing modes will likely launch Alfvén waves. Such waves have been investigated in the context of gamma-ray bursts from neutron star quakes (Blaes et al. 1989;Begelman et al. 1993), as well as recently in the context of coincident Fast-Radio and X-ray bursts (Yuan et al. 2020(Yuan et al. , 2022) ) from a galactic magnetar (Andersen et al. 2020;Bochenek et al. 2020). In this picture, Alfvén waves propagate from the neutron star surface into the magnetosphere, become nonlinear, and cause the emission of an ultra-relativistic blast wave (flare), as well as direct dissipation through magnetic reconnection (Yuan et al. 2020).

Whether and how these dynamics can happen in a compact binary magnetosphere is the subject of this work. Numerical simulations of magnetospheric transients have matured considerably over the past years (Palenzuela et al. 2010;Palenzuela 2013;Palenzuela et al. 2013a;Alic et al. 2012;Ponce et al. 2015;Nathanail 2020;Carrasco & Shibata 2020;Carrasco et al. 2021;East et al. 2021;Most & Philippov 2020, 2023a, 2022, 2023b), including studies of the flaring dynamics in single (Parfrey et al. 2013;Carrasco et al. 2019;Mahlmann et al. 2023;Sharma et al. 2023) and binary neutron star magnetospheres (Most & Philippov 2020, 2023a,b, 2022). Leveraging this progress, in this work we perform magnetospheric simulations of crustal shattering flares and their nonlinear Alfvén-wave dynamics.

The aim of this paper is to clarify the magnetospheric dynamics of the crustal shattering scenario. In particular, we target the fraction of energy that can be converted into the final flare/blast wave state, and what the feedback of the flare on the magnetosphere might be, e.g., in the context of shock powered radio emission mechanisms (Beloborodov 2020;Metzger et al. 2019).

To this end, we present global general-relativistic forcefree electrodynamics simulations of Alfvén-wave dynamics in binary neutron star and black hole -neutron star systems prior to merger. We demonstrate the nonlinear evolution of Alfvén waves, show the formation of blast waves at large distances, and place constraints on the energy conversion efficiency in flare launching and reconnection-mediated dissipation. Our results provide insight into the crustal shattering transient picture.

This work is structured as follows: We summarize the basic motivation of our study in Sec. 2. We then describe the numerical setup and the binary configuration we study in Sec. 3. Our main results and analyses are presented in Sec. 4, before concluding in Sec. 5.

BASIC PICTURE

In the resonant crust shattering picture, tidal interactions between the neutron star and its companion excite modes at the crust-core interface prior to merger (Tsang et al. 2012;Penner et al. 2012;Sullivan et al. 2022;Zhang et al. 2022). The crust-core interface (i-)mode frequency is (Neill et al. 2021)

depending on the nuclear symmetry energy. Other modes have also been considered, including g-modes (Passamonti et al. 2021;Kuan et al. 2021a,b). For the i-mode, the excitation timescale is (Tsang et al. 2012;Neill et al. 2021)

where M is the binary chirp mass. The excitation injects energy into the crust, which gets released once the crust fractures (Tsang 2013;Neill et al. 2022) leading to estimated luminosities of (Tsang 2013),

where B NS is the neutron star surface magnetic field strength. One caveat is that the energy stored in the entire crust could be as low as (Baiko & Chugunov 2018)

where r NS is the neutron star radius. Such energies would likely result in luminosities lower than those of gamma-ray burst precursors (Coppin et al. 2020).

Launching Alfvén waves

Crustal shattering will inject perturbations into the surrounding magnetosphere. These perturbations have been mainly associated with non-radial modes (McDermott et al. 1988), including either torsional (Suvorov et al. 2022), or interface (i-)modes at the crust-core boundary (Tsang et al. 2012). Such non-radial modes will predominantly excite Alfvén waves along the magnetic field lines threading the crust. Modeling of a crustal quake has shown that the efficiency factor for converting the crustal energy release into Alfvén waves is f ≃ 1% (Bransgrove et al. 2020), due to large energy absorption by the superfluid core. This reduces the overall luminosity injected into the magnetosphere,

Alfvén waves propagate as an initial perturbation along the magnetic field lines. Following Blaes et al. (1989), the resulting perturbation, δB/B, of the magnetic field, B, at the surface of the neutron star is

As a result of flux freezing and a magnetic dipole background, B ∝ r -3 , decaying with radius, r, the perturbation steepens as

The initial perturbation, (δB/B) surface , will then grow as it propagates outwards, until it becomes nonlinear, approximately corresponding to δB ≃ B (Blaes et al. 1989). In case of multi-polar fields, B ∼ r -n (n > 3), non-linearity can be achieved earlier, and closer to the binary. At the same time, closed field lines can be shorter (compared to a dipole) placing tighter constraints on the required initial perturbation to reach non-linearity. The modified topology of the field would also affect the geometry of wave propagation. We will mainly discuss dipole fields in the following estimates.

Irrespective of the precise structure of the magnetic field, steepening will only efficiently happen on closed magnetic field lines which decay faster with radius than open ones, implying that nonlinear steepening needs to happen before the waves reach the orbital light cylinder,

where M = m 1 + m 2 is the binary total mass, m i are the component masses, a is the orbital separation, and G the gravitational constant. We further define the mass ratio Q = m 1 /m 2 ≥ 1. The largest radial distance an Alfvén wave can propagate inside the light cylinder is from the lighter companion m 2 in the direction of m 1 .foot_0

The lighter binary component is offset from the center of mass by aQ/(1 + Q), implying from Eq. ( 7) that the smallest surface perturbation leading to nonlinear steeping is

To leading order, the time to merger is (Peters 1964)

leading to (δB/B) min ∝ t -9/16 merger . Here c is the speed of light. Equations ( 9) and (10) together yield the minimum perturbation required for the system to feature nonlinear flaring dynamics and are depicted in Fig. 1. Alfvén waves will non-linearly steepen for initial perturbations as low as δB/B ∼ 10 -3 -10 -4 if launched seconds before merger, which is less than the perturbation estimated in Eq. ( 6). This back-of-the-envelop estimate implies that nonlinear Alfvén-wave dynamics is possible for crustal shattering induced by the i-mode resonance of Eq. (1) as indicated by the shaded regions in Fig. 1, see also Tsang et al. (2012); Neill et al. (2022)).

Observational signatures

Assuming for the moment that any Alfvén wave will nonlinearly steepen before leaving the orbital light cylinder, the main question concerns its subsequent evolution. Will it produce a flare/relativistic blast wave that could potentially power electromagnetic transients? How are the flares formed, what energy do they carry? 10 -4 10 -3 10 -2 10 -1 10 0 10

1 10 2 t launcht merger [s] 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 (δB/B) min NS -NS BH -NS Figure 1. Minimum required perturbation that becomes nonlinear within the orbital light cylinder. The launching time, t launch , is given for black hole (BH) -neutron star (NS) and binary neutron star (NS-NS) systems relative to the time of merger tmerger. We adopt m2 = 1.4 M⊙, rNS = 12 km and vary Q = 1 (NS-NS) to Q = 4 (BH-NS). Shaded regions correspond to the expected time window of the crust-core interface mode resonance from Eq. (1).

While nonlinearity is commonly associated with the formation of charge-starved electric zones, numerical work suggests that the dissipation in charge-starved Alfvén waves is likely small (Chen et al. 2022a).

Insights into potential emission mechanisms come from related simulations of star quakes on isolated neutron stars (Yuan et al. 2020(Yuan et al. , 2022)). In these simulations, the Alfvén wave non linearly stretches, opens, and pinches off part of the closed magnetosphere. This leads to the ejection of a large plasmoid, which can be interpreted as a flare, similar to magnetar giant flares (Parfrey et al. 2013) or binary neutron star flaring dynamics (Most & Philippov 2020). Consistent with Yuan et al. (2020), our simulations confirm this picture, Sec. 4. This is similar in spirit to flare collision models proposed for gamma-ray precursors from crustal shattering (Tsang et al. 2012;Neill et al. 2022). We further show that the conversion efficiency, η < 1, from the Alfvén wave into the flare implies a final luminosity of

Furthermore, the flare has a trailing current sheet leading to direct dissipation with luminosities

in X-rays (Yuan et al. 2020(Yuan et al. , 2022)), see also Most & Philippov (2020, 2023a,b) for similar estimates in compact binary magnetospheres.

Collisionless shocks can produce coherent radio emission when a blast wave interacts with a surrounding wind (Metzger et al. 2019;Beloborodov 2020). This process critically depends on the wind properties including the amount of baryon loading, as quantified by the wind magnetization σ w . Following Beloborodov (2020), we estimate the emission peak frequency as

where Γ w is the wind Lorentz factor, and r s the distance from the stellar surface to shock formation. Depending on the formation distance of the collisionless shock such a mechanism can produce radio emission.

One of the challenges of this model is that the wind density cannot be too low, or stated differently, the wind luminosity needs to be a fraction of the flare luminosity. For an orbiting dipole (Hansen & Lyutikov 2001;Ioka & Taniguchi 2000)

≈ 2.8 • 10 39 erg s

where Ω is the orbital angular speed, and µ the magnetic dipole moment. In the context of multiple flares, the first flare could enhance the background in its upstream leading to L wind ∼ 10 -2 L flare (Yuan et al. 2020). As we will demonstrate in Sec. 4.4, this result also holds for binary magnetospheres. The energy, E FRB , available to power a fast-radio burst is

where t burst is the observed burst duration, and ε s is the efficiency of the shock maser (Plotnikov & Sironi 2019;Sironi et al. 2021). In the scenario where the previous flare enhances the wind (see discussion above), the burst duration is (Yuan et al. 2020;Beloborodov 2020)

where ∆t flare ∼ t i-mode is the delay between subsequent flares. Here we have assumed that the crust can fracture multiple times on the resonant pumping timescale t i-mode .

METHODS

We numerically study the propagation and nonlinear dynamics of Alfvén waves in a compact binary magnetosphere, with an application to resonant shattering flares in neutron star mergers. This requires us to model the dynamics of the compact binary magnetosphere from the neutron star radius, r NS , where the perturbation is seeded, to ≃ 50 -100 r NS ≃ 2 -4 r LC , where the flaring dynamics and propagation happens. Bridging this large range of scales is not computationally feasible with kinetic models of magnetospheric plasmas. In addition, we are mainly interested in the global dynamics of the flaring process, as well as bulk estimates of energy conversion efficiencies in the system.

The simplest framework which includes Alfvén waves and can model the dynamics of the pair-plasma filled force-free magnetosphere (Goldreich & Julian 1969;Hansen & Lyutikov 2001;Lyutikov 2019) is force-free electrodynamics (Gruzinov 1999;Palenzuela et al. 2010;Parfrey et al. 2013;Carrasco et al. 2019) in a full generalrelativistic setting, e.g., (Alic et al. 2012;Palenzuela 2013;Paschalidis et al. 2013;Mahlmann et al. 2021;Most & Philippov 2023b;Kim et al. 2024). Specifically, we numerically solve the Maxwell equations,

where J ν is the electric current, and we have defined the Maxwell field strength tensor,

as well as its dual, * F µν , and electric, E µ , and magnetic, B µ , fields in the frame defined by a normal observer n µ . The force-free conditions are enforced by adopting a resistive current (Alic et al. 2012;Palenzuela 2013),

where q is the electric charge density and σ is the parallel conductivity, see, e. 2013), we model the binary spacetime in general relativity using a fixed-orbit numerical solution of Columns mark the fixed orbital frequency, Ω orb , light cylinder, rLC, separation, a, the dimensionless spin of the primary component, χ1, and the inclination angle, θB, of the magnetic dipole field of the secondary neutron star. For the BH-NS system we choose m 1/BH = 5 M⊙ and m 2/NS = 1.4 M⊙, whereas the binary NS system is modelled as an equal mass system with total mass M = 2.7 M⊙. All systems have χ2 = 0, and in the case of NS-NS systems, we set the primary magnetic field to be always anti-aligned with the orbital axis.

the extended conformally thin sandwich system (Grandclement 2006;Taniguchi et al. 2007Taniguchi et al. , 2008;;Foucart et al. 2008;Tacik et al. 2016). The system is solved numerically using the Kadath/FUKA (Papenfort et al. 2021;Grandclement 2010) code. We carry out simulations of both binary neutron star and black hole -neutron star systems. We assume that the stellar light cylinder lies outside the orbital light cylinder, i.e., the stars rotate slower than the orbital period (Bildsten & Cutler 1992;Zhu et al. 2018), in which case stellar spins will not affect the Alfvén wave dynamics. We, therefore, model the neutron star(s) as irrotational. This assumption is valid as long as the orbital frequency is larger than the stellar rotation frequency at the time of launching the initial perturbation. In the case of an i-mode resonance, Eq. ( 1), this would hold except for millisecond pulsars. This implies that the only relevant intrinsic scales in the system (apart from the initial perturbation δB/B) are the orbital light cylinder, r LC , and to a lesser extent the neutron star radius, r NS , and the mass ratio q. System parameters are given in Tab. 1. Simulations use the same numerical framework as in previous works simulating the dynamics of compact binary magnetospheres (Most & Philippov 2020, 2022, 2023a,b). We solve the general-relativistic force-free electrodynamics system using a relaxation approach with an effective parallel conductivity, which enforces the main force-free condition, E i B i = 0 (Alic et al. 2012;Palenzuela et al. 2013a). Different from previous studies using this code, we implement a version of the ECHO scheme (Del Zanna et al. 2007) (see also Most et al. (2019)), combined with WENO-Z reconstruction (Borges et al. 2008), and a Rusanov Riemann solver (Rusanov 1961). The time integration for the stiff relaxation term is handled using the IMEX-SSP433 scheme (Pareschi & Russo 2005), making the scheme formally third-order accurate in time, and fourth-order accurate in space. In addition, we manually enforce the second force-free constraint (E 2 < B 2 ) by rescaling E 2 < 0.999B 2 where necessary. This approach is consistent with the physically expected behavior at magnetospheric shock formation (Beloborodov 2023). The numerical code is implemented on top of the adaptive mesh-refinement (AMR) infrastructure of the AMReX framework (Zhang et al. 2019).

Computational challenges & GPU computing

The large numerical grid resolutions required for tracking the waves as they propagate from the neutron star surface to the orbital light cylinder and become nonlinear poses a computational challenge. Previous work on Alfvén-wave propagation in single-star magnetospheres has employed resolutions of (8192 × 4096) in axisymmetry (Yuan et al. 2020), and (2560 3 ) in three dimensions (Yuan et al. 2022). In this work, we have found it necessary to employ a total of 8 levels of mesh refinement, with a total number of more than 10 9 grid points. The outer boundary extends to about 2, 400 km in order to decouple reflections and spurious artefacts arising near the corners of the outer domain.

Given these computational challenges, we were forced to port our numerical code entirely to GPUs using the functionality of the AMReX infrastructure. With a performance of about 58 Mio cell-updates/s on a single A100 GPU, we have been able to carry out these simulations on 192 A100 GPUs on the NERSC Perlmutter compute cluster. The total cost of the simulations shown here is about 10, 000 GPU node-hours. To check the validity of our results, we have further run one simulation with substantially larger refinement regions on 2,400 V100 GPUs on OLCFs Summit system, which resulted in consistent dynamics.

Wave launching

We model the feedback of crustal shattering on the magnetosphere as a wavetrain of monochromatic Alfvén waves with a fixed period and number of wave cycles. For realistic crust shattering, likely complicated fracturing patterns (Thompson et al. 2017) and resonances in the crust (Bransgrove et al. 2020) will inject a variety of different wavelengths. Since the tidal interaction will mainly excite non-radial modes (Suvorov et al. 2022), we model the surface motion following recent three-dimensional simulations of self-consistent neutron star oscillations. In particular, Sagert et al. (2023) reported that axisymmetric surface motion may resemble a bottle cap (polar) twist, see also Parfrey et al. (2013). Inspired by these results, we adopt

where θ is the azimuthal angle in the coordinate of the neutron star, and we use (σ = 1, θ 0 = 0) for the upper and (σ = -1, θ 0 = π) for the lower hemisphere of the star. The resulting velocity perturbation δv ϕ ∼ δΩ r NS sin θ is shown in Fig. 2. We inject this surface motion via an angular velocity at the outer layer of the neutron star surface, Ω NS = -Ω orb + δΩ, for a total of four periods, after an initial time of t launch , chosen such that the magnetosphere can fully develop after the initial numerical transient has propagated away. In the case of an equal-mass neutron star binary the perturbation is identical and sets in at the same time for both neutron stars. Since the wave will grow regardless of its initial value at the surface (see Sec. 2), we fix a value of δΩ that comfortably leads to steepening of the wave within the computational domain we simulate. Specifically, we adopt a fixed perturbation of θ ¯= π/3, ωM ⊙ = 2π/20, and A = 3πΩ 0 , which corresponds to δB/B ∼ 0.05 at the surface of the neutron star. We show the initial wave launching phase in Fig. 2.

RESULTS

We investigate the nonlinear dynamics of Alfvén waves in the magnetosphere of a compact binary in its final orbits prior to merger. Specifically, we are interested in understanding the formation and dynamics of flares caused by resonant shattering of the neutron star crust due to gravitational tides in the inspiral (Tsang et al. 2012;Penner et al. 2012). We present results as follows: first, we summarize the state of the magnetosphere prior to launching Alfvén waves, Sec. 4.1; second, we describe the magnetospheric dynamics of Alfvén waves and the formation of flares, Sec. 4.2, before presenting a parameter study on the neutron star magnetic field, Sec. 4.3; finally, we summarize several properties of the flares, Sec. 4.4.

Background magnetosphere

We provide a summary of the background compact binary magnetosphere prior to launching of Alfvén waves. The magnetosphere consists of different regions. Similar to a pulsar magnetosphere, the rotational motion of the orbit creates a light cylinder at r LC = c/Ω orb , dividing the magnetosphere into zones of closed, r < r LC , and open, r > r LC , field lines. As discussed in Sec. 2, it is the closed zone where Alfvén waves will primarily steepen. Unlike a pulsar magnetosphere (Spitkovsky 2006;Kalapotharakos & Contopoulos 2009), the binary magnetosphere will in most cases not reach a quasistationary state (Palenzuela et al. 2013a;Most & Philippov 2020;Carrasco et al. 2021).

For neutron star magnetic moments not aligned with the orbital axis, the magnetosphere will feature periodic transient phases. 2 In the case of black hole -neutron star systems, magnetic field lines from the neutron star will thread the black hole (Paschalidis et al. 2013;Most & Philippov 2023b). These connected flux tubes will lead to a built-up of a net twist due to orbital motion, ultimately dissipating energy in a unipolar inductor-like scenario (Goldreich & Lynden-Bell 1969;Piro 2012;Lai 2012). For small magnetic field inclinations, θ B ≲ 45 • , the energy in the twist will largely be dissipated in a trailing current sheet of the black hole (see also magnetic draping (Lyutikov 2023)). For inclinations larger than this threshold, the twisted flux tubes will explode periodically leading to flaring outbursts (Most & Philippov 2023b). The precise threshold for flaring may depend critically on the reconnection rate in the black hole current sheet, which force-free simulations cannot model correctly (Parfrey et al. 2019). From a topological point of view this should happen about twice per orbit (Cherkis & Lyutikov 2021).

For binary neutron star systems the situation is similar. For systems where both stars carry a magnetic field with magnetic moments not fully aligned in the same direction, connected flux tubes can be formed. Once a critical twist is reached, these flux tubes will always reconnect and produce flares (Most & Philippov 2022). Since both neutron stars will have crusts that can shatter, this could lead to counter-propagating Alfvén-waves, which could interact and, in principle, drive turbulence (Goldreich & Sridhar 1995;TenBarge et al. 2021;Ripperda et al. 2021).

Alfvén wave launching and dynamics

Having established the background state of the binary magnetosphere, we now focus on the dynamics of Alfvén waves from their launch through their full nonlinear evolution. In principle, this process will not only produce Alfvén waves and especially should the fields contain any initial twist, the nonradial shearing motion from crust shattering will also produce fast waves (Mahlmann et al. 2023;Sharma et al. 2023), which can steepen into monster shocks and also power transients (Chen et al. 2022b;Beloborodov 2023;Most et al. 2024). However, given that the energetics will be a fraction of that of the Alfvén waves, we restrict to the latter.

The initial monochromatic wave package, see Sec. 3.2, is inject after a simulation time of t launch ≃ 1, 500 km/c. At this time, the initial transient from the start of the simulation has fully propagated away from the inner regions, providing a clean background magnetosphere. As the wave package is propagating outward in the dipole magnetosphere of the neutron star it was launched in, the waves begin to steepen with radius following Eq. ( 7).

Before discussing the simulation outcomes, we first briefly review general propagation properties of Alfvén waves. Mathematically, Alfvén wave propagation is described by the relativistic Elsässer equations (Chandran et al. 2018;Elsasser 1950) in the force-free limit,

where (u

which is an advection equation with the Alfvén vector, v ν A∓ , with an interaction term depending on the orientation and propagation direction of the waves. In particular, the interaction term will vanish if there are no counter-propagating Alfvén waves, which holds for any magnetization. Unlike for magnetosonic fast waves, this means that multiple Alfvén waves propagating in the same direction will not interact, even if they propagate at different speeds 3 and catch up with one another.

3 In force-free electrodynamics, the Alfvén speed v A = √︁ σ/(1 + σ) → 1, as the magnetization σ → ∞. While this is appropriate for magnetospheric dynamics, it also means that individual wave sections can never intersect.

Equipped with this theoretical framework, we can now present and interpret our numerical results. We show the propagation of Alfvén waves in Fig. 3 for both black hole -neutron star (mixed) binaries and binary neutron star systems. In the background dipole magnetosphere (shown in the co-orbiting frame) Alfvén waves will propagate as B ϕ components of the magnetic field.

In the case of a mixed binary (top row), the waves can propagate in the direction of the black hole companion or away from it, which adds an additional redshift of the wavelength of the wave package as seen at infinity. Following discussion of Eq. ( 25), we find no interaction of the Alfvén waves as expected from the absence of counter-propagating waves in the initial setup. As the waves propagate outward, an E ϕ component of the electric field forms, which in this magnetic topology indicates the production of an outgoing flare, see the formation of a plasmoid (closed magnetic field line region) on the equator (Yuan et al. 2021;Bernardi et al. 2024;Mahlmann et al. 2024). Concurrently, the closed zone is ripped open by the nonlinear perturbation of the Alfvén waves, pinching off the outer part of the closed zone and ejecting a large plasmoid on the orbital plane (since the magnetic moment is aligned with the orbital axis). This demonstrates the launching of a flare due to Alfvén-wave nonlinearity. Trailing the flare is a current sheet, which will cause substantial additional dissipation of energy, potentially powering X-ray transients (Yuan et al. 2020), see Sec. 4.4 for discussion.

The situation shown here for the mixed binary is very similar to the case of a star quake on a magnetar (Yuan et al. 2020(Yuan et al. , 2022)). There, the quake launched an Alfvén wave localized to a specific flux bundle, which caused the emission of a flare within the hemisphere of the initial quake. The crustal shattering case examined here with its large-scale perturbation is an extreme limit of that scenario, in which all field lines are excited. That said, the field lines on which fast waves can be produced are strongly constrained by the latitude of the perturbation relative to the magnetic moment. Open field lines cannot convert Alfvén waves into flares/fast waves (Yuan et al. 2021), meaning that Alfvén waves will continue to propagate outward without dissipating. We can spot this as an E ϕ ∼ 0 region in the wave shell corresponding to the north and south pole of the magnetic moment, where field lines are open.

While the overall situation is very similar in the case of a neutron star binary (bottom row), there are subtle differences. Since the anti-aligned configuration we show in Fig. 3 features connected flux tubes, counterpropagating Alfvén waves can be launched along those flux tubes, where they will intersect. This requires the presence of a secondary neutron star and does not happen for the mixed binary case. According to Eq. ( 25), such waves will interact, as shown in the polar region above the binary. In principle, interacting Alfvén waves can trigger a turbulent cascade via three wave interactions (Goldreich & Sridhar 1995). Numerical studies indicate, however, that many interactions (> 100 crossings) are required before a turbulent cascade sets in (Ripperda et al. 2021), which we do not observe here (our simulations feature at most 1-2 crossings). Instead, conversion of Alfvén waves into fast waves (Yuan et al. 2021;Chen et al. 2022bChen et al. , 2024) ) is observed in the inner region of both systems, as denoted by the appearance of large amplitude wave-like patterns in E ϕ . These could provide additional shock-powered dissipation through monster shocks (Chen et al. 2022b;Beloborodov 2023;Most et al. 2024).

Secondary fast waves along the equator are also present (just as in the mixed binary case), but are suppressed compared to the large scale background twist dynamics of the emerging connected flux tube (large positive E ϕ region in Fig. 3. This is a result of orbital flaring, which for realistic orbital separations will however be strongly suppressed compared to the Alfvén wave dynamics (Most & Philippov 2020, 2022). In either case, there is clear production of a flare with trailing current sheets in the equatorial plane, just as the mixed binary. This different angular shape of the flares leads us to conclude that at least for misaligned binary neutron star systems, there will be an angular dependence on the emission of resonant shattering flares.

Dependence on the magnetic field geometry

The alignment of the magnetic field prior to merger is not known, although dipole-dipole interactions may drive the system to near (anti-) alignment (Aykroyd et al. 2023;Lander & Jones 2018). Since tidal excitations have a preferred direction relative to the orbital axis, nonradial shearing motion induced by tidal forces will likely be parallel to the orbital plane. Misalignment between the shearing motion and the magnetic axis (akin to a misalignment between the spin and the magnetic moment in a pulsar) will affect flaring emission, since only Alfvén waves launched on closed field lines will cause a flaring transient. One would then expect a reduction in the flare power as a function of the magnetic inclination angle, θ B . Accordingly, we have performed several simulations investigating the impact of the inclination angle θ B , shown in Fig. 4.

Starting out with mixed binaries (top panels), Alfvén waves in partially inclined configurations, θ B = 30 • , 60 • (left and middle panels), cause similar flaring to aligned configurations, θ B = 0 • (top right panel of Fig. 3). The main differences are in the orientation of the transient current sheet (visible via plasmoids formed in the tearing unstable sheet) and the ejection angle of the main plasmoid, which is approximately perpendicular to the magnetic field axis. Due to the increased twist (likely of shorter field lines), the emission of secondary fast waves behind the main Alfvén wave train is enhanced. For θ B = 90 • , the situation is different. Because now the toroidal perturbation axis and the magnetic moment are misaligned, most Alfvén waves are launched on open field lines. These Alfvén waves steepen fully to nonlinearity and can be seen as large amplitude waves distorting the field lines but they do not dissipate.

For binary neutron star systems (bottom panels) and motivated by the results of previous inclination surveys for compact binary magnetospheres (Most & Philippov 2022;Ponce et al. 2015), we consider the case where one of the neutron stars has a magnetic field anti-aligned with the orbital axis, whereas for the other neutron star we vary the magnetic moment inclination. The main 3 but the final times t -t launch = 2.2 ms for different initial magnetic field inclinations, θB. We do not repeat the θB = 0 • case, as it is overall similar to θB = 30 • . In all cases, Alfvén waves drive flares detaching from the compact binary magnetosphere.

difference in the background field is that rather than having long connected flux tubes extending to high latitudes, misaligned systems feature complicated twist geometries, where the twisted flux tubes extend largely in the angular rather than the radial direction (Cherkis & Lyutikov 2021). In this sense, a small fraction of Alfvén waves can become trapped, never reaching sufficiently large radii to become non-linear. These trapped Alfvén waves may ultimately convert to fast waves, albeit at a lower efficiency than the main dynamics discussed here (Yuan et al. 2021;Bernardi et al. 2024).

For larger inclination angles, θ B = 60 • , similar dynamics sets in. Fast waves and flares are being formed in both hemispheres, but with a stronger inclination angle. In the upper hemisphere, the complicated dynamics of the interacting field lead to an inflated fluxtube connecting both stars, ultimately detaching and causing an orbital-motion powered flare, which is almost orthogonal to the one powered by nonlinear Alfvén waves.

Overall, the inclination angle of the dipole magnetic filed will introduce anisotropies in subsequent emission. In the following, we quantify the corresponding energetics of the flares.

Flare properties

We finally provide a quantitative analysis of the flare properties. To this end, we focus on a representative case for a neutron star binary systems with θ B = 60 • , for which we show the final flaring state in Fig. 5. We are mainly interested in the available dissipation channels. In that regard, the properties of the flare as well as the energy conversion is very similar to results for single star magnetospheres (Yuan et al. 2020;Bernardi et al. 2024;Mahlmann et al. 2024).

First, we consider reconnection. Reconnection in pulsar current sheets can power high energy transients (Uzdensky & Spitkovsky 2014). In binary magnetospheres this dissipation channel contributes about 10% of the overall Poynting flux (Most & Philippov 2020;Carrasco & Shibata 2020;Carrasco et al. 2021;Most & Philippov 2023b), which is in line with Alfvén-wave dynamics in single neutron star magnetospheres (Yuan et al. 2020(Yuan et al. , 2022)).

We show the local dissipation rate J i E i in the top panel. Nonlinearity of the Alfvén waves has created a strong perturbation of the orbital current sheet, which fragments into plasmoids (Uzdensky & Spitkovsky

2014

). In principle, mergers of these plasmoids can produce radio nano-shot emission (Philippov et al. 2019;Lyubarsky 2019) provided that the magnetic field strength is on the order of 10 8 G in the current sheet. A similar behavior has been found for flares caused by orbital motion (Most & Philippov 2020). Energetically, this emission channel will likely be subdominant.

We similarly find that in the core of the large plasmoid of the Alfvén-wave-driven flare substantial dissipation is present. This has been argued by Yuan et al. (2020) to be a source of X-ray emission. Since our prescription for dissipation is ad-hoc (using an artificial parallel conductivity), we cannot meaningfully quantify this dissipation beyond the upper percentage limit quoted above.

Second, and beyond dissipation in current sheets, we also quantify the properties of the converted flare directly. To do so, we compute the drift Lorentz factor, γ d , which is very large γ d ≫ 10 (middle panel). Indeed, the flare will continue to expand, accelerate and become a flat pancake resembling a blast wave with γ d ≫ 100 at distances of 10 13 cm (Yuan et al. 2020). Portions of the initial Alfvén wave on open field lines become nonlinear but have low drift Lorentz factor (upper right corner of the center panel in Fig. 5). This feature is common among all of our simulations.

Finally, we analyze the energetics of the flare. Since our simulations do not model the crust itself, we cannot correlate the energy contained in the flare with that in the crust directly. However, we can compute how much of the energy we initially inject in the Alfvén wave ends up in the flare by computing the energy as a volume integral,

contained in a spherical shell S.

Here e EM = 1 2

(︁ E 2 + B 2 )︁ is the electromagnetic energy density. Choosing S to contain the initial Alfvén wave/final flare we obtain the energies E Alfve ´n and E Flare , respectively. The latter time is chosen to correspond to a propagation distance of at least 500 km from the binary when the Alfvén wave has become fully nonlinear. The resulting energy conversion efficiency, η, is given in Fig. 6. Regardless of dipole inclination the efficiency is always ≃ 20% for black hole -neutron star and ≃ 40% for binary neutron star systems when considering only the part of the outflow that attains high Lorentz factor (γ d > 5), and ≃ 40% -60% in total. Since disentangling the energy content of Alfvén and fast waves is already complicated in axisymmetric spacetimes (Bernardi et al. 2024;Yuan et al. 2021;Mahlmann et al. 2024), we use this distinction as a simple proxy for the energy that ends up in the flare. Overall, these numbers are consistent with findings for single star magnetospheres (Yuan et al. 2020).

In addition to the total luminosity, we also comment on the spatial distribution of the electromagnetic energy (Poynting) flux (Gourgoulhon 2012),

where β i is the coordinate shift vector. This is shown in the lower panel of Fig. 5. Here, most of the Poynting flux is carried by the steepened blast wave. However, in the upstream of this blast wave the Poynting flux is only 1-2 orders smaller. This implies that the local background wind is enhanced substantially over its background state, which in general will be even lower, see Eq. ( 14). While this does not have implications for a single Alfvén wave event as we show here, it may have implications for multiple staggered Alfvén waves as are expected to be launched in a crustal shattering event (Tsang 2013;Neill et al. 2021). Indeed, Yuan et al. (2020) found that this enhancement may strongly increase the feasibility for shock-maser powered radio emission (Metzger et al. 2019;Beloborodov 2020). This finding critically motivates discussion on potential radio transients in Sec. 2.

CONCLUSIONS

Tidal interactions in coalescing binaries involving neutron stars can resonantly excite modes at the crust-core interface (Penner et al. 2012;Tsang et al. 2012) that shatter the neutron star crust (Tsang 2013) and launch waves into the compact binary magnetosphere prior to merger. Since relevant modes will largely be nonradial (Suvorov et al. 2022;Sagert et al. 2023), such waves should predominantly be Alfvén waves that can nonlinearly steepen in the background of the closed zone (Blaes et al. 1989). In the nonlinear phase, the Alfvén waves can convert into a blast wave/flare (Yuan et al. 2020(Yuan et al. , 2022)), which can then power electromagnetic transients.

In this work, we provided an investigation of nonlinear Alfvén-wave dynamics in compact binary magnetospheres using general-relativistic force-free electrodynamics simulations. The computational need to capture the nonlinear phase of Alfvén-wave dynamics required simulations on hundreds of GPUs. Since our simulations did not model the dynamics of the crust, we have injected a monochromatic wave package and tracked its evolution through the fully nonlinear phase. In reality, resonances in the crust (Bransgrove et al. 2020) and fracturing patterns on the surface (Thompson et al. 2017), will likely cause a variety of different oscillation frequencies to be injected into the magnetosphere. This could affect the potential intersection of waves (when realistic magnetizations are included), the presence or onset of Alfvenic turbulence (Goldreich & Sridhar 1995;Ten-Barge et al. 2021;Ripperda et al. 2021), and the power contained in the resulting flares. Consistent models for the shattering of the neutron star crust will be needed to answer these questions.

Our results demonstrate the formation of a high velocity flare in binary neutron star and black hole -neutron star systems irrespective of the inclination of the magnetic field.

Overall, at least 20% -40% of the energy initially injected gets converted into the final flare, slightly less than for single star magnetospheres (Yuan et al. 2020). Additionally, the flare enhanced the wind in its upstream, thus increasing the prospect for shockpowered radio emission in a multi-flare scenario, see Sec. 2. In line with previous work on compact binary magnetospheres (Carrasco & Shibata 2020;Most & Philippov 2020;Carrasco et al. 2021;Most & Philippov 2022, 2023b), the system features copious reconnectionmediated dissipation channels which can power secondary X-ray emission.

Our results represent an important step toward understanding the feasibility of the crustal shattering scenario (Blaes et al. 1989;Tsang et al. 2012;Penner et al. 2012). We demonstrate numerically that the magnetospheric dynamics indeed give rise to flare launching based on local shearing motion on the neutron star surface. Equipped with such a framework, it is possible to study the propagation and dynamics of realistic crustal shattering or star quake events (Bransgrove et al. 2020;Sagert et al. 2023), potentially enabling more direct connections between astrophysical observables and the nuclear physics of the crust (Neill et al. 2021(Neill et al. , 2023)).