INTRODUCTION

Neutron star mergers are not only sources of gravitational waves, but are also accompanied by electromagnetic counterparts. Such counterparts provide insights into properties of dense matter in neutron stars and their magnetic fields. Canonical counterparts include a prompt gamma-ray burst (GRB (Meszaros 2006)) and a kilonova (Li & Paczynski 1998;Metzger 2020). Both were observed from the neutron star merger GW170817 (see, e.g., Abbott et al. (2017) for a summary). The GRB followed the merger with a 2-s delay and was likely emitted by a blast wave breaking out at the photosphere of the merger ejecta (Murguia-Berthier et al. 2017;Gottlieb et al. 2018;Xie et al. 2018;Beloborodov et al. 2020). The kilonova was emitted on a day timescale and powered by nuclear decay in the expanding ejecta.

The explosion picture usually includes collimated relativistic jets from the black hole promptly formed after the merger. The jets are extremely bright sources of gamma-rays if observed face-on, and also emit a broad-band afterglow when they decelerate in an am-bient medium (Gottlieb et al. 2018). Similar jets are expected if the merger forms a short-lived neutron star with a debris accretion disk, see Rezzolla et al. (2011); Ruiz et al. (2016); Kawamura et al. (2016); Kiuchi et al. (2023) for magnetohydrodynamic (MHD) simulations of such systems. Such models predict the collapse of the neutron star in milliseconds to a second after the merger, triggered by accretion or the loss of differential rotation.

However, intermediate mass mergers may not promptly collapse to black holes and form jets (Bauswein et al. 2013;Köppel et al. 2019;Bauswein et al. 2020;Kashyap et al. 2022;Tootle et al. 2021;Kölsch et al. 2022;Schianchi et al. 2024). Instead, they may form a long-lived neutron-star remnant (Baiotti & Rezzolla 2017;Radice et al. 2020). Its initial differential rotation is expected to generate magnetic fields up to B = 10 16 G, which buoyantly emerge from the remnant and form its external magnetosphere (Kluzniak & Ruderman 1998;Most & Quataert 2023;Combi & Siegel 2023;Giacomazzo & Perna 2013;Giacomazzo et al. 2015;Kiuchi et al. 2015;Mösta et al. 2020;Aguilera-Miret et al. 2023). Then, the remnant cools and forms a young magnetar. It likely has a twisted magnetosphere, filled with e ± pairs of a small mass density ρ. Its expected magnetization parameter σ = B 2 / (︁ 4πρc 2 )︁ is huge, exceeding 10 10 (Beloborodov 2023).

Eventually, the remnant emits its angular momentum in a magnetized wind, loses its rotational support, and can collapse into a black hole (Lasky et al. 2014;Ravi & Lasky 2014;Dall'Osso et al. 2015). The spindown occurs on the timescale t sd ≈ c 3 I/µ 2 Ω 2 ≈ 10 4 µ -2 33 (ν/300 Hz) -2 s, where I ≈ 10 45 g cm 2 is remnant's moment of inertia, ν = Ω/2π is its rotation rate and µ is the dipole moment of its magnetosphere (we normalized it to 10 33 G cm 3 ). The delayed collapse after t ∼ t sd occurs suddenly, on a ms timescale. Proposed electromagnetic signals from the delayed collapse include low-frequency waves (Falcke & Rezzolla 2014) and GRBs emitted by collimated jets (Ciolfi & Siegel 2015).

Previous simulations of neutron-star collapse followed the dynamics of the external magnetosphere either in vacuum or using force-free electrodynamics (FFE) (Baumgarte & Shapiro 2003;Lehner et al. 2012;Dionysopoulou et al. 2013;Palenzuela 2013;Most et al. 2018). Both frameworks do not allow plasma heating. Furthermore, both neglect plasma inertia and so are unable to track magnetospheric shock formation. By contrast, a full magnetohydrodynamic analysis predicts monster shocks, which can generate gamma-rays (Beloborodov 2023).

This Letter reports the first magnetohydrodynamic simulations of the magnetosphere evolution in the dynamic spacetime of the collapsing magnetar. It demonstrates shocks and ejection of a hot outflow that will emit a GRB. We also find that the outflow carries information about quasi-normal modes of the nascent black hole, which may be observed in the GRB time profile.

METHODS

We set up a simple initial state: a dipole magnetic field is attached to a rotating star with µ ∥ Ω, using the vector potential given by Shibata et al. (2011) with a surface magnetic field of B * ∼ 10 16 G.

We artificially reduce the rotation rate from a realistic ν ≳ 300 Hz to 56 Hz, so that the light cylinder R LC = c/Ω is beyond our computational box. This simplification avoids the challenging preparation of an equilibrium rotating magnetosphere with the equatorial current sheet at r > R RL . Rotation is still essential, as it enables the ring-down effect when the star collapses into a black hole (e.g. Kokkotas & Schmidt (1999)). The amplitude of ring-down oscillations scales with an-gular momentum. Therefore, we also simulate collapse with a high initial ν = 900 Hz; this additional simulation misses details of the outer magnetospheric dynamics but demonstrates the enhanced ring-down effect on the magnetized outflow from the nascent black hole.

We use the RNS code (Stergioulas & Friedman 1995) to initialize the star as an unstable general relativistic polytrope (with polytropic coefficient K=100) of mass M = 1.7M ⊙ and equatorial coordinate radius R ⋆ ≈ 12 km; the details of its internal structure are unimportant, as we focus on the external magnetosphere. The magnetosphere has a low mass density ρ and a large magnetization parameter σ ≫ 1. At time t = 0, we add a small pressure perturbation and the star begins to collapse.

The simulation tracks the spacetime of the collapsing star and evolves its magnetosphere according to general relativistic magnetohydrodynamics (GRMHD) equations (Duez et al. 2005). Since the plasma magnetosphere has a very small effective resistivity on the dynamical timescales of interests t dyn ∼ R/c, it can be treated as an ideal conductor everywhere except the sites of magnetic reconnection, which develops later in the simulation (and is mediated by numerical resistivity).

We use the Frankfurt/IllinoisGRMHD (FIL) code (Most et al. 2019;Etienne et al. 2015) which is built on top of the Einstein Toolkit infrastructure (Loffler et al. 2012). The spacetime dynamics is tracked using the Z4c formulation of the Einstein equations (Hilditch et al. 2013). We use moving puncture coordinates (Alcubierre et al. 2003) in the simulation and the presentation of results below. The GRMHD equations assume an ideal fluid with thermal pressure proportional to eρc 2 , where e is the fluid energy density including rest mass. The equations are evolved using the ECHO scheme (Del Zanna et al. 2007) with vector potentialbased constraint transport (Etienne et al. 2010(Etienne et al. , 2012)).

No spatial symmetries are imposed during the simulation. We employ a fixed three-dimensional Cartesian grid with 6 levels of mesh refinement; the highest resolution has 78 grid points per R ⋆ . The grid extends to 3.25 × R ⋆ (750 km) in each direction. The boundaries of the computational domain are far from the shock and do not affect the results presented below.

Unlike most MHD simulations in dynamical spacetimes (e.g. Liu et al. 2008;Kiuchi et al. 2014;Palenzuela et al. 2015;Ciolfi et al. 2017;Most et al. 2019), our simulation follows the magnetosphere with a high magnetization parameter (see also Paschalidis et al. (2015); Ruiz et al. (2016)). In particular, the initial "background" magnetosphere has σ bg ∼ 25, and its perturbation during the collapse leads to ultra-relativistic motions with Lorentz factors γ ∼ 10. Such simulations are challenging in terms of accuracy requirements and numerical stability of the algorithm, especially when MHD is coupled to a dynamically evolved spacetime.

We have changed the code to improve its robustness (see also Most (2023)), including different primitive recovery schemes (Kastaun et al. 2021;Kalinani et al. 2022), drift floors (Ressler et al. 2017), and bounds on σ and β = 8πP/B 2 where P is the fluid pressure and B is the magnetic field in fluid rest frame. Specifically, before black hole formation, we enforce bounds, 10 < σ < 50, and β > 0.02, outside of the star. We distinguish between magnetospheric and stellar matter using a passive scalar (see also Ref. Parfrey & Tchekhovskoy (2017)). Even with all these improvements, we have found that the fourth-order derivative corrector performs poorly in the shock region. Since the constraint-transport algorithm prevents us from switching it off selectively, we have disabled it everywhere. Thus, the simulation maintains second-order accuracy, different from all previous GRMHD simulations carried out with FIL (Most et al. 2019(Most et al. , 2021;;Chabanov et al. 2023;Most & Quataert 2023;Most 2023).

RESULTS

The evolution observed in the simulation may be summarized as follows.

Collapse and wave launching

The collapsing star quickly forms a black hole with the apparent horizon radius R h ≈ 2.5 km and the ergosphere around it. Effectively, the magnetospheric footprints on the star are quickly pulled in from R ⋆ to R h , and the star's magnetic flux Ψ is now threading the smaller sphere of radius R h . As a result, a strong quasimonopolar magnetic field B ∼ Ψ/r 2 is created in the radial zone R h < r < R ⋆ . This inner zone with the amplified magnetic pressure launches a compressive wave into the surrounding magnetosphere, which propagates with nearly speed of light, v wave /c ≈ 1σ -1 bg . The compressive MHD waves (known as "fast magnetosonic modes") have electric field E ∥ k × B bg where k is the (approximately radial) wavevector and B bg is the initial background dipole magnetic field. So, the wave has a toroidal electric field E ϕ . The launched wave of E ϕ continues to propagate outward for the rest of the simulation; its snapshot is shown in Figure 1 at a late time t = 0.67 ms (as measured by a distant observer), near the end of the simulation. The magnetic field ahead of the wave is the dipole B bg , and the magnetic field behind the wave is close to the split monopole configuration. The wave dynamics at radii r ≫ R h can be approximately described neglecting general-relativistic corrections.

Monster shock formation

Magnetosonic waves were recently shown to accelerate plasma to a huge radial 4-velocity u r = γv r /c (Beloborodov 2023). This effect occurs when the wave reaches radius R × = (cµ 2 /8L) 1/4 in the equatorial plane, where µ = r 3 B bg and L is the wave power. In the FFE limit (σ bg → ∞) fluid expansion in the rarefaction phase of the compressive wave would diverge at R × as E 2 -B 2 touches zero and u r → -∞. For a finite σ bg , fluid develops a finite u r ∝ σ bg . The ultrarelativistic fluid motion is directed toward the star and the wave immediately develops a monster shock.

Shock formation in magnetosonic waves has been demonstrated by kinetic plasma simulations (Chen et al. 2022) and by MHD calculations using characteristics (Beloborodov 2023). For waves with frequency ω ≫ c/R × a simple analytical MHD solution has been obtained. It demonstrates that at r > R × the wave profile E(tr/c) develops a plateau of width W p ∼ c/ω where E 2 ≈ B 2 . The plateau forms a linear accelerator, so the wave pushes the fluid 4-velocity to a huge value (Be- loborodov 2023)

The accelerated flow dissipates its energy in a monster shock.

The monster shock appears in our simulation where E 2 -B 2 approaches zero, as predicted. We also observe the development of an E-plateau where u r develops a steep linear profile, reaching values consistent with Eq. (1) (see Fig. 2). The shock is the sudden jump of u r from a large negative value u r ∼ -6 back to a moderate u r . These unique features of monster shocks are clear in the simulation despite numerical inaccuracies accumulated in the shock region. We also observe the expected strong heating localized at the shock (see the temperature panel in Fig. 2).

The magnetization parameter σ bg ∼ 25 used in the simulation is far below its real value in a magnetar, and the shock strength should be scaled to a larger σ bg according to Eq. ( 1). The corresponding large u r ∝ σ bg will make the shock highly radiative, i.e. the accelerated flow will radiate its energy before crossing the shock and joining the downstream flow (Beloborodov 2023). Future simulations could attempt to track radiative transfer with self-consistent creation of e ± pairs (Beloborodov 2021), which immediately make the flow optically thick.

The dissipated energy is inevitably thermalized behind the shock, creating an opaque, radiation-dominated outflow. Our simulation assumes that the released energy remains trapped in the fluid. The observed relativistic outflow trails the shock, which expands with speed c.

Black hole ring-down

The ring-down of the nascent spinning black hole lasts ∼ 100R h /c. It involves quasi-periodic oscillations of the horizon with frequencies characteristic of black hole quasi-normal modes, whose amplitude decays exponentially with time (e.g. Kokkotas & Schmidt (1999)). For a stationary black hole, the quasi-monopolar magnetic flux Ψ threading the horizon would be stuck for a significant time -its decay would be controlled by magnetic reconnection in the equatorial plane on a timescale ≳ 100R h /c (Bransgrove et al. 2021). By contrast, the oscillating black hole quickly and quasi-periodically sheds magnetic flux, losing significant δΨ each oscillation period. The discharged δΨ forms a quasi-periodic MHD outflow with the characteristic cusps of magnetic field lines in the equatorial plane. The cusps are inherited from the earlier split-monopole shape of field lines at the time of their decoupling from the oscillating black hole. This effect is seen in the simulation with the reduced rotation rate ν = 56 Hz (Fig. 1) and becomes stronger in the simulation with fast rotation ν = 900 Hz (Fig. 3).

Previous collapse simulations with a vacuum or FFE magnetosphere showed electromagnetic waves from ringdown (Baumgarte & Shapiro 2003;Lehner et al. 2012;Palenzuela 2013;Most et al. 2018). Vacuum ring-down is usually described as a coupling of black-hole quasinormal modes to outgoing gravitational and electromagnetic waves (Teukolsky 1972(Teukolsky , 1973) ) using Newman-Penrose scalars

Here C µνκλ is the Weyl curvature tensor and F µν is the electromagnetic tensor; vectors m and n are conveniently chosen as m = (θ ˆ+ iϕ ˆ)/ √ 2 and n = (t ˆr ˆ)/ √ 2 (so that m, n, and l = (t ˆ+ r ˆ)/ √ 2 form an orthonormal null tetrad). Then, ϕ 2 ∝ E θ + iE ϕ represents two polarization states of the outgoing electromagnetic waves (in MHD, E θ and E ϕ correspond to the fast magnetosonic and Alfvén waves, respectively). We have verified that the dominant (quadrupole) component of ψ 4 observed in our simulation is consistent with the quasinormal mode computed using qnm code (Stein 2019). We have also calculated the dominant (dipole) component of ϕ 2 , which also approximately matches the corresponding quasi-normal mode frequency (Fig. 4 shows the evolution of dominant spherical harmonics in Im ϕ 2 and ψ 4 ). Note however that ϕ 2 was designed for vacuum electromagnetic fields, and so the oscillation of ϕ 2 may not accurately represent the modulation of MHD outflow. The frequency of compressive modulations observed in Fig. 3 may be directly estimated as ν mod = v/λ, where v ≲ c is the outflow speed and λ ∼ 7 km is the spatial modulation period. The subsequent evolution of the balding black hole will crucially depend on the resistivity in the current sheet (Bransgrove et al. 2021;Selvi et al. 2024). In the absence of a controlled resistivity in our simulations, we defer the study of the late phase to future work.

3.4. Gamma-ray burst.

Our simulation demonstrates that the magnetospheric destruction during the neutron-star collapse involves strong dissipation and creates a powerful magnetized outflow with a characteristic peak duration ∼ 1 ms. The hot outflow is launched behind the leading monster shock and has a quasi-periodic tail. The modulated tail is generated by the nascent black hole, as it rings down and quickly sheds most of its magnetic flux initially inherited from the neutron star. The modulation frequency lies naturally in the kilohertz band, suggesting a connection with recently reported kilohertz quasi-periodic oscillations in some GRBs (Chirenti et al. 2023).

Note that the magnetosphere around the neutron star prior to collapse has a minute plasma density. Therefore, the explosion triggered by collapse is practically clean from baryons. The created e ± plasma in the hot outflow is initially opaque to scattering. Most of the heat density U is contained in trapped blackbody radiation U ≈ aT 4 (a is the radiation constant), as the photons far outnumber the e ± pairs. The outflow expansion to large radii is not followed by our simulation, however its basic features can be predicted in analogy with the well-known "fireball" model for cosmological GRBs (Paczynski 1986;Goodman 1986). The outflow Lorentz factor γ will grow and its temperature T will drop due to adiabatic cooling. Eventually, most of e ± pairs annihilate, and the trapped photons are released, producing a GRB. The burst duration is set by the outflow duration, which lasts only a few ms after the collapse.

The burst spectrum will peak at photon energies ∼ 3k B T γ, where k B is the Boltzmann constant. Simplest GRB models assume adiabatic outflows with no magnetic fields, and then γT ≈ const = T 0 (Paczynski 1986). In that case, the observed temperature is weakly changed during the outflow expansion and remains close to the initial temperature. Outflows predicted by our simulation are magnetically dominated, and their temperature may be affected by additional dissipation of magnetic energy at large Regulation of the GRB spectral peak in dissipative outflows is discussed in Beloborodov (2013).

The energy of the GRB outflow is set by the precollapse energy of the neutron-star magnetosphere E m ∼ 10 47 µ 2 33 erg. During the collapse the magnetospheric energy is amplified by compression, and then most of it is ejected in the outflow. Assuming that ∼ 0.1 of this energy is eventually emitted in the GRB, we can roughly estimate the GRB luminosity as L ∼ 0.1E m /1 ms ∼ 10 49 B 2 15 erg/s, where B 15 = B/10 15 G is the magnetic field near the neutron star prior to collapse. Note that the outflow (and hence the GRB) is anisotropic, but not strongly collimated, unlike jet-powered GRB models.

During its lifetime prior to collapse, the rotating magnetar produced a magnetic wind outside the dipole magnetosphere, with a termination shock of a large radius behind the merger ejecta (e.g., Metzger et al. 2007;Dessart et al. 2009) 1 . The collapse launches a relativistic shock, which will continue to expand into the cold wind. The shock is expected to emit a fast radio burst (FRB) at radii r ∼ 10 13 cm (Beloborodov 2020). This suggests a mechanism for a delayed FRB from neutron star mergers, emitted together with the delayed GRB. However, the FRB can hardly escape through the surrounding shell of mass M ej ∼ 10 -2 M ⊙ ejected earlier by the merger (Bhardwaj et al. 2023;Radice et al. 2023). A similar problem is faced by the recently proposed FRB-GRB connection (Rowlinson et al. 2023).

The massive shell is also a threat to the GRB predicted by our simulation. The GRB will not be blocked by M ej if the magnetar collapse occurs with a sufficient delay ∆t ≳ 10 hr after the merger, so that the massive ejecta expands to radius R ∼ 3 × 10 14 (∆t/10 hr) cm and its optical depth to gamma-rays τ ∼ κM ej /4πR 2 drops below unity (κ is the ejecta opacity to Compton scattering). ris et al. 2020), qnm (Stein 2019), scipy (Virtanen et al. 2020)