Neutron stars have solid crusts threaded by strong magnetic fields. Perturbations in the crust can excite nonradial oscillations which can in turn launch Alfvn waves into the magnetosphere. In the case of a compact binary close to merger involving at least one neutron star this can happen through tidal interactions causing resonant excitations that shatter the neutron star crust. We present the first numerical study that elucidates the dynamics of Alfvn waves launched in a compact binary magnetosphere. We seed a magnetic field perturbation on the neutron star crust which we then evolve in fully generalrelativistic forcefree electrodynamics using a GPUbased implementation. We show that Alfvn waves steepen nonlinearly before reaching the orbital light cylinder form flares and dissipate energy in a transient current sheet. Our results predict radio and Xray precursor emission from this process.
['Neutron stars', 'Magnetospheric radio emissions', 'Plasma astrophysics', 'Radio transient sources', 'X-ray transient sources', 'Astronomical simulations', 'General relativity', 'Gravitational wave sources', '1108', '998', '1261', '2008', '1852', '1857', '641', '677', 'Astrophysics - High Energy Astrophysical Phenomena', 'Physics - Plasma Physics']
Nonlinear Alfvnwave Dynamics and Premerger Emission from Crustal Oscillations in Neutron Star Mergers
{"Nonlinear Alfv'en-Wave Dynamics and Premerger Emission from Crustal Oscillations in Neutron Star Mergers": '1, 2  \n1, 3  \n3, 4  \n<!-- image -->  \nElias R. Most  \n,  \nYoonsoo Kim  \n,  \nKaterina Chatziioannou  \n,  \nand Isaac Legred  \n1 TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA 2 Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA 3 Department of Physics, California Institute of Technology, Pasadena, CA 91125, USA 4 LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA', 'ABSTRACT': "Neutron stars have solid crusts threaded by strong magnetic fields. Perturbations in the crust can excite non-radial oscillations, which can in turn launch Alfv'en waves into the magnetosphere. In the case of a compact binary close to merger involving at least one neutron star, this can happen through tidal interactions causing resonant excitations that shatter the neutron star crust. We present the first numerical study that elucidates the dynamics of Alfv'en waves launched in a compact binary magnetosphere. We seed a magnetic field perturbation on the neutron star crust, which we then evolve in fully general-relativistic force-free electrodynamics using a GPU-based implementation. We show that Alfv'en waves steepen nonlinearly before reaching the orbital light cylinder, form flares, and dissipate energy in a transient current sheet. Our results predict radio and X-ray precursor emission from this process.  \nKeywords: Black holes(162), General relativity(641), Gravitational wave sources (677), Magnetospheric radio emissions (998), Neutron stars (1108), Plasma astrophysics(1261), Radio bursts(1339), Radio transient sources (2008), X-ray transient sources(1852)", '1. 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. 2017, 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),  \nemost@caltech.edu  \nkchatziioannou@caltech.edu  \nwhich 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;  \nZhang et al. 2022) or in neutron star oceans in the inspiral (Sullivan et al. 2022, 2023).  \nIn 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. 2021, 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 ( Ozel & Freire 2016; Watts et al. 2016).  \nWhile 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'en 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, 2022) from a galactic magnetar (Andersen et al. 2020; Bochenek et al. 2020). In this picture, Alfv'en 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).  \nWhether 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  \n2020, 2023a,b, 2022). Leveraging this progress, in this work we perform magnetospheric simulations of crustal shattering flares and their nonlinear Alfv'en-wave dynamics.  \nThe 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'en-wave dynamics in binary neutron star and black hole - neutron star systems prior to merger. We demonstrate the nonlinear evolution of Alfv'en 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.  \nThis 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.", '2. 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)  \nf i -mode ≈ 130 -170 Hz , (1)  \ndepending 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)  \nt i -mode ≈ 0 . 04 s ( 1 . 2 M ⊙ M ) 5 / 6 ( 150 Hz f i -mode ) 11 / 6 , (2)  \nwhere 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),  \nL crust ≃ 10 47 ( B NS 10 13 G ) 2 erg / s , (3)  \nwhere B NS is the neutron star surface magnetic field strength. One caveat is that the energy stored in the  \nentire crust could be as low as (Baiko & Chugunov 2018)  \nE crust ≳ 2 . 5 × 10 41 ( r NS 12 km ) 2 erg , (4)  \nwhere 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).', "2.1. Launching Alfv'en 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'en 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'en 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,  \nL Alfv'en = f L crust ≃ 10 45 ( f 0 . 01 )( B NS 10 13 G ) 2 erg / s . (5)  \nAlfv'en 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  \nδB B ∣ ∣ ∣ ∣ surface ≃ 0 . 02 ( L Alfv'en 10 45 erg / s ) 1 / 2 ( 10 13 G B NS ) . (6)  \nAs a result of flux freezing and a magnetic dipole background, B ∝ r -3 , decaying with radius, r , the perturbation steepens as  \nδB B = ( δB B ) surface ( r r NS ) 3 / 2 . (7)  \nThe 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,  \nsteepening 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,  \nr LC ≃ a 3 / 2 ( GM ) 1 / 2 , (8)  \nwhere 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'en wave can propagate inside the light cylinder is from the lighter companion m 2 in the direction of m 1 . 1 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  \n( δB B ) min = r NS 3 / 2 [ a ( √ a GM + Q 1+ Q )] 3 / 2 ≈ ( GM ) 3 / 4 r NS 3 / 2 a 9 / 4 . (9)  \nTo leading order, the time to merger is (Peters 1964)  \nt merger = 5 256 c 5 G 3 a 4 m 3 2 Q (1 + Q ) , (10)  \nleading 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'en 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'en-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)).", '2.2. Observational signatures': "Assuming for the moment that any Alfv'en 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?  \nFigure 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 t merger . We adopt m 2 = 1 . 4 M ⊙ , r NS = 12km 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).  \n<!-- image -->  \n-  \nWhile nonlinearity is commonly associated with the formation of charge-starved electric zones, numerical work suggests that the dissipation in charge-starved Alfv'en waves is likely small (Chen et al. 2022a).  \nInsights into potential emission mechanisms come from related simulations of star quakes on isolated neutron stars (Yuan et al. 2020, 2022). In these simulations, the Alfv'en 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'en wave into the flare implies a final luminosity of  \nL flare ≃ 2 · 10 44 ( η 0 . 2 ) ( f 0 . 01 )( B NS 10 13 G ) 2 erg / s . (11)  \nFurthermore, the flare has a trailing current sheet leading to direct dissipation with luminosities  \nL x -ray ≤ 0 . 1 L Alfv'en , (12)  \nin X-rays (Yuan et al. 2020, 2022), see also Most & Philippov (2020, 2023a,b) for similar estimates in compact binary magnetospheres.  \nCollisionless 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  \nν peak ≃ 5 GHz ( 10 13 cm r s ) 3 / 2 ( Γ w 20 )( E flare 10 41 erg ) 1 / 2 , (13)  \nwhere Γ 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.  \nOne 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)  \nL wind = 4 15 c 5 µ 2 a 2 Ω 6 , (14) ≈ 2 . 8 · 10 39 erg s ( B NS 10 13 G ) 2 ( M 2 . 7 M ⊙ ) 3 ( 120 km a ) 7 ,  \nwhere Ω 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.  \nThe energy, E FRB , available to power a fast-radio burst is  \nE FRB ≃ 10 38 erg ( ε s 10 -3 ) ( L flare 10 44 erg / s )( t burst 1 ms ) , (15)  \nwhere 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)  \nt burst = ∆ t flare 4 σ w ( L w L flare ) 1 / 2 (16)  \n≃ 1 ms 1 σ w ( t i -mode 40 ms )( L wind 10 -2 L flare ) 1 / 2 , (17)  \nwhere ∆ 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 .", '3. METHODS': "We numerically study the propagation and nonlinear dynamics of Alfv'en 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.  \nThe simplest framework which includes Alfv'en 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,  \n∇ µ ∗ F µν = 0 , (18)  \n∇ µ F µν = -J ν , (19)  \nwhere J ν is the electric current, and we have defined the Maxwell field strength tensor,  \nF µν = n µ E ν -n ν E µ + ε µνκ B κ , (20)  \nas 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),  \nJ i = qε ijk E j B k B 2 + σ E j B j B 2 B i , (21)  \nwhere q is the electric charge density and σ is the parallel conductivity, see, e.g., Palenzuela (2013); Dionysopoulou et al. (2013); Paschalidis et al. (2013); Palenzuela et al. (2013b,a); Ponce et al. (2015); East et al. (2021); Carrasco et al. (2021) for applications in compact binary magnetospheres.  \nFollowing Most & Philippov (2023b) and Paschalidis et al. (2013), we model the binary spacetime in general relativity using a fixed-orbit numerical solution of  \nTable 1. System parameters of the black hole (BH) - neutron star (NS) and NS-NS binaries considered in this work. Columns mark the fixed orbital frequency, Ω orb , light cylinder, r LC , 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 BHNS 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.  \nthe extended conformally thin sandwich system (Grandclement 2006; Taniguchi et al. 2007, 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'en 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.  \nSimulations 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  \nthe 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 . 999 B 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).", '3.1. 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'en-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.  \nGiven 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.", '3.2. Wave launching': "We model the feedback of crustal shattering on the magnetosphere as a wavetrain of monochromatic Alfv'en 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.  \nSince the tidal interaction will mainly excite non-radial  \nFigure 2. Initial wave packet of Alfv'en waves launched from the neutron star (gray disk) for the aligned configuration, θ B = 0. The orbital angular momentum points along the positive z -axis. Shown is the out-of-plane magnetic field component B ϕ , rescaled to the surface magnetic field strength B NS , as well as the relative surface perturbation velocity δv ϕ .  \n<!-- image -->  \nmodes (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  \nδ Ω( θ ) = σA exp ( -( θ -θ 0 ) 2 2 ¯ θ 2 ) cos ( ω ( t -t 0 )) , (22)  \nwhere θ 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.", '4. RESULTS': "Weinvestigate the nonlinear dynamics of Alfv'en 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'en waves, Sec. 4.1; second, we describe the magnetospheric dynamics of Alfv'en 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.", '4.1. Background magnetosphere': "We provide a summary of the background compact binary magnetosphere prior to launching of Alfv'en 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'en 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).  \nFor 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  \nhole 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).  \nFor 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'en-waves, which could interact and, in principle, drive turbulence (Goldreich & Sridhar 1995; TenBarge et al. 2021; Ripperda et al. 2021).", "4.2. Alfv'en wave launching and dynamics": "Having established the background state of the binary magnetosphere, we now focus on the dynamics of Alfv'en waves from their launch through their full nonlinear evolution. In principle, this process will not only produce Alfv'en 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'en waves, we restrict to the latter.  \nThe 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).  \nBefore discussing the simulation outcomes, we first briefly review general propagation properties of Alfv'en waves. Mathematically, Alfv'en wave propagation is described by the relativistic Elsasser equations (Chandran et al. 2018; Elsasser 1950) in the force-free limit,  \n∇ ν ( z µ ± z ν ∓ ) = -( 3 4 z µ ± z ν ∓ + 1 4 z µ ∓ z ν ± + 1 2 g µν ) ∂ ν b 2 b 2 , (23)  \nz µ ± = u µ d ∓ b µ b , (24)  \nwhere ( u d ) i = γ d v i with γ d = √ B 2 / ( B 2 -E 2 ) being the drift Lorentz factor, v i = ε ijk E j B k /B 2 being the drift velocity, and b µ = B µ /γ d the magnetic field in the drift frame. Alfv'en waves propagate as perturbations,  \nFigure 3. Alfv'en wave propagation in a compact binary spacetime with aligned magnetic fields, θ B = 0 · , in a frame co-rotating with the binary. (Top) Neutron star (NS) - black hole (BH), (Bottom) NS-NS. In the NS-NS case, the companion NS has an anti-aligned magnetic field. Shown are three different times, t , stated relative to the time of launch t launch . Colors denote the out-of-plane magnetic, B ϕ , and electric, E ϕ , fields normalized to the poloidal magnetic field B P and the field strength at the surface B NS respectively.  \n<!-- image -->  \nδz µ ± = z µ ± -〈 z µ ± 〉 , relative to the background 〈 z µ ± 〉 . Here, the sign in z µ ± indicates the propagation of the wave relative to the magnetic field direction. Most importantly, the relativistic Elsasser equations imply that for a magnetically dominated plasma (TenBarge et al. 2021),  \nv ν A ∓ ∂ ν δz µ ± = -δz ν ∓ ∂ ν δz µ ± , (25)  \nwhich is an advection equation with the Alfv'en 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'en waves, which holds for any magnetization. Unlike for magnetosonic fast waves, this means that multiple Alfv'en waves propagating in the same direction will not interact, even if they propagate at different speeds 3 and catch up with one another.  \nEquipped with this theoretical framework, we can now present and interpret our numerical results. We show the propagation of Alfv'en 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'en waves will propagate as B ϕ components of the magnetic field.  \nIn 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'en 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  \nAlfv'en 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'en-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.  \nThe situation shown here for the mixed binary is very similar to the case of a star quake on a magnetar (Yuan et al. 2020, 2022). There, the quake launched an Alfv'en 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'en waves into flares/fast waves (Yuan et al. 2021), meaning that Alfv'en 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.  \nWhile 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'en 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'en 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'en waves into fast waves (Yuan et al. 2021; Chen et al. 2022b, 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).  \nSecondary 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  \npositive 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'en 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.", '4.3. 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'en 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.  \nStarting out with mixed binaries (top panels), Alfv'en 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'en 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'en waves are launched on open field lines. These Alfv'en waves steepen fully to nonlinearity and can be seen as large amplitude waves distorting the field lines but they do not dissipate.  \nFor 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  \nFigure 4. Same as Fig. 3 but showing 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'en waves drive flares detaching from the compact binary magnetosphere.  \n<!-- image -->  \ndifference 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'en waves can become trapped, never reaching sufficiently large radii to become non-linear. These trapped Alfv'en 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).  \nFor 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'en waves.  \nOverall, 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.", '4.4. 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).  \nFirst, 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'en-wave dynamics in single neutron star magnetospheres (Yuan et al. 2020, 2022).  \nWe show the local dissipation rate J i E i in the top panel. Nonlinearity of the Alfv'en waves has created a strong perturbation of the orbital current sheet, which fragments into plasmoids (Uzdensky & Spitkovsky  \nFigure 5. Late stage compact magnetosphere evolution for a binary neutron star system with magnetic field inclination angle θ B = 60 · in the meridional plane. Shown in color are the current sheets, in terms of the dissipative current J i , the drift Lorentz factor γ d , as well as the radial Poynting flux, S r EM normalized to its maximum value, S EM , flare , in the outgoing flare.  \n<!-- image -->  \n∣∣∣  \n∣∣∣  \nη  \nFigure 6. Conversion efficiency η of the initial Alfv'en wave into the final flare. Shown are results for black hole (BH) - neutron star (NS) and NS-NS systems. Solid lines correspond to waves with drift Lorentz factor γ d > 5, which approximately filters out the fast wave component, whereas dashed lines correspond to the entire energy in the final outflow. Error bars quantify the uncertainty in choosing the cut-off on γ d .  \n<!-- image -->  \n2014). 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.  \nWe similarly find that in the core of the large plasmoid of the Alfv'en-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.  \nSecond, 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'en 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.  \nFinally, 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'en wave ends  \nup in the flare by computing the energy as a volume integral,  \nE S = ∫ S d V e EM = ∫ S d V 1 2 ( E 2 + B 2 ) , (26)  \ncontained 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'en wave/final flare we obtain the energies E Alfv'en 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'en 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'en 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).  \nIn addition to the total luminosity, we also comment on the spatial distribution of the electromagnetic energy (Poynting) flux (Gourgoulhon 2012),  \nS i EM = ε ijk E j B k -β i e EM , (27)  \nwhere β 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'en wave event as we show here, it may have implications for multiple staggered Alfv'en 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.", '5. 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  \nshatter 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'en waves that can nonlinearly steepen in the background of the closed zone (Blaes et al. 1989). In the nonlinear phase, the Alfv'en waves can convert into a blast wave/flare (Yuan et al. 2020, 2022), which can then power electromagnetic transients.  \nIn this work, we provided an investigation of nonlinear Alfv'en-wave dynamics in compact binary magnetospheres using general-relativistic force-free electrodynamics simulations. The computational need to capture the nonlinear phase of Alfv'en-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; TenBarge 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.  \nOur 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.  \nOverall, 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.  \nOur 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  \nstudy 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, 2023).", 'ACKNOWLEDGMENTS': 'The authors are grateful for insightful discussions with A. Beloborodov, O. Blaes, A. Bransgrove, D. Neill, A. Philippov, E.S. Phinney, J. Read, B. Ripperda, L. Sironi, A. Suvorov, C. Thompson, Y. Yuan, and B. Zhang. ERM acknowledges support from the National Science Foundation under Grant No. AST-2307394. IL and KC were supported by a grant from the Simons Foundation (MP-SCMPS-00001470). The simulations were performed on DOE NERSC supercomputer  \nPerlmutter under grant m4575. Additional simulations were performed on DOE OLCF Summit under allocation AST198. Preliminary work was also performed on NSF Frontera supercomputer under grant AST21006. ERM further acknowledges the use of Delta at the National Center for Supercomputing Applications (NCSA) through allocation PHY210074 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.  \nSoftware: AMReX (Zhang et al. 2019), FUKA (Papenfort et al. 2021), Kadath (Grandclement 2010), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020)', 'REFERENCES': "Abbott, B. P., et al. 2017, Astrophys. J. Lett., 848, L13,  \n```\ndoi: 10.3847/2041-8213/aa920c Alexander, K. D., et al. 2017, Astrophys. J. Lett., 848, L21, doi: 10.3847/2041-8213/aa905d Alic, D., Mosta, P., Rezzolla, L., Zanotti, O., & Jaramillo, J. L. 2012, Astrophys. J., 754, 36, doi: 10.1088/0004-637X/754/1/36 Andersen, B. C., et al. 2020, Nature, 587, 54, doi: 10.1038/s41586-020-2863-y Aykroyd, C., Bourgoin, A., Le Poncin-Lafitte, C., Mathis, S., & Angonin, M.-C. 2023, Astron. Astrophys., 675, A32, doi: 10.1051/0004-6361/202346171 Baiko, D. A., & Chugunov, A. I. 2018, Mon. Not. Roy. Astron. Soc., 480, 5511, doi: 10.1093/mnras/sty2259 Begelman, M. C., Meszaros, P., & Rees, M. J. 1993, MNRAS, 265, L13, doi: 10.1093/mnras/265.1.L13 Beloborodov, A. M. 2020, Astrophys. J., 896, 142, doi: 10.3847/1538-4357/ab83eb -. 2023. https://arxiv.org/abs/2307.12182 Bernardi, D., Yuan, Y., & Chen, A. Y. 2024. https://arxiv.org/abs/2405.02199 Bildsten, L., & Cutler, C. 1992, Astrophys. J., 400, 175, doi: 10.1086/171983 Blaes, O., Blandford, R., Goldreich, P., & Madau, P. 1989, ApJ, 343, 839, doi: 10.1086/167754 Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59, doi: 10.1038/s41586-020-2872-x Borges, R., Carmona, M., Costa, B., & Don, W. 2008, Journal of Computational Physics, 227, 3191, doi: 10.1016/j.jcp.2007.11.038\n```  \nBransgrove, A., Beloborodov, A. M., & Levin, Y. 2020, ApJ, 897, 173, doi: 10.3847/1538-4357/ab93b7  \nBransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, Astrophys. J. Lett., 958, L9,  \n```\ndoi: 10.3847/2041-8213/ad0556 Carrasco, F., & Shibata, M. 2020, Phys. Rev. D, 101, 063017, doi: 10.1103/PhysRevD.101.063017 Carrasco, F., Shibata, M., & Reula, O. 2021, Phys. Rev. D, 104, 063004, doi: 10.1103/PhysRevD.104.063004 Carrasco, F., Vigan'o, D., Palenzuela, C., & Pons, J. A. 2019, Mon. Not. Roy. Astron. Soc., 484, L124, doi: 10.1093/mnrasl/slz016 Chandran, B. D. G., Foucart, F., & Tchekhovskoy, A. 2018, J. Plasma Phys., 84, 905840310, doi: 10.1017/S0022377818000387 Chatziioannou, K. 2020, Gen. Rel. Grav., 52, 109, doi: 10.1007/s10714-020-02754-3 Chen, A. Y., Yuan, Y., Beloborodov, A. M., & Li, X. 2022a, Astrophys. J., 929, 31, doi: 10.3847/1538-4357/ac59b1 Chen, A. Y., Yuan, Y., & Bernardi, D. 2024. https://arxiv.org/abs/2404.06431 Chen, A. Y., Yuan, Y., Li, X., & Mahlmann, J. F. 2022b. https://arxiv.org/abs/2210.13506 Cherkis, S. A., & Lyutikov, M. 2021, Astrophys. J., 923, 13, doi: 10.3847/1538-4357/ac29b8 Chornock, R., et al. 2017, Astrophys. J. Lett., 848, L19, doi: 10.3847/2041-8213/aa905c Cooper, A. J., Gupta, O., Wadiasingh, Z., et al. 2022, doi: 10.1093/mnras/stac3580 Coppin, P., de Vries, K. D., & van Eijndhoven, N. 2020, PhRvD, 102, 103014, doi: 10.1103/PhysRevD.102.103014 Corsi, A., et al. 2024. https://arxiv.org/abs/2402.13445 Cowperthwaite, P. S., et al. 2017, Astrophys. J. Lett., 848, L17, doi: 10.3847/2041-8213/aa8fc7\n```  \nWada, T., Shibata, M., & Ioka, K. 2020, PTEP, 2020, 103E01, doi: 10.1093/ptep/ptaa126  \nWang, J.-S., Peng, Z.-K., Zou, J.-H., Zhang, B.-B., & Zhang, B. 2020, Astrophys. J. Lett., 902, L42,  \ndoi: 10.3847/2041-8213/abbfb8  \n- Watts, A. L., et al. 2016, Rev. Mod. Phys., 88, 021001, doi: 10.1103/RevModPhys.88.021001  \nXiao, S., et al. 2022. https://arxiv.org/abs/2205.02186  \n- Yuan, Y., Beloborodov, A. M., Chen, A. Y., & Levin, Y. 2020, Astrophys. J. Lett., 900, L21,  \ndoi: 10.3847/2041-8213/abafa8  \n- Yuan, Y., Beloborodov, A. M., Chen, A. Y., et al. 2022, Astrophys. J., 933, 174, doi: 10.3847/1538-4357/ac7529\n- Yuan, Y., Levin, Y., Bransgrove, A., & Philippov, A. 2021, Astrophys. J., 908, 176, doi: 10.3847/1538-4357/abd405\n- Zhang, W., Almgren, A., Beckner, V., et al. 2019, Journal of Open Source Software, 4, 1370,  \ndoi: 10.21105/joss.01370  \n- Zhang, Z., Yi, S.-X., Zhang, S.-N., Xiong, S.-L., & Xiao, S. 2022, Astrophys. J. Lett., 939, L25,  \ndoi: 10.3847/2041-8213/ac9b55  \n- Zhu, X., Thrane, E., Oslowski, S., Levin, Y., & Lasky, P. D. 2018, Phys. Rev. D, 98, 043002,  \ndoi: 10.1103/PhysRevD.98.043002"}