∎

11institutetext: N.Sarin and P.D. Lasky22institutetext: School of Physics and
Astronomy, Monash University, Vic 3800, Australia OzGrav: The ARC Centre of
Excellence for Gravitational Wave Discovery, Clayton VIC 3800, Australia
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# The evolution of binary neutron star post-merger remnants: a review

Nikhil Sarin and Paul D. Lasky

(Received: date / Accepted: date)

###### Abstract

Two neutron stars merge somewhere in the Universe approximately every 10
seconds, creating violent explosions observable in gravitational waves and
across the electromagnetic spectrum. The transformative coincident
gravitational-wave and electromagnetic observations of the binary neutron star
merger GW170817 gave invaluable insights into these cataclysmic collisions,
probing bulk nuclear matter at supranuclear densities, the jet structure of
gamma-ray bursts, the speed of gravity, and the cosmological evolution of the
local Universe, among other things. Despite the wealth of information, it is
still unclear when the remnant of GW170817 collapsed to form a black hole.
Evidence from other short gamma-ray bursts indicates a large fraction of
mergers may form long-lived neutron stars. We review what is known
observationally and theoretically about binary neutron star post-merger
remnants. From a theoretical perspective, we review our understanding of the
evolution of short- and long-lived merger remnants, including fluid, magnetic-
field, and temperature evolution. These considerations impact prospects of
detection of gravitational waves from either short- or long-lived neutron star
remnants which potentially allows for new probes into the hot nuclear equation
of state in conditions inaccessible in terrestrial experiments. We also review
prospects for determining post-merger physics from current and future
electromagnetic observations, including kilonovae and late-time x-ray and
radio afterglow observations.

###### Keywords:

Binary neutron star mergers remnants gravitational-waves gamma-ray bursts
kilonovae

##  1 Introduction

The coincident gravitational-wave and electromagnetic observations of binary
neutron star merger GW170817/GRB170817A Abbott  _et al._ (2017a, b, c) was a
watershed moment, signaling the beginning of a new field of multi-messenger
gravitational-wave astronomy. Gravitational-wave emission from the inspiral
phase was detected by the Advanced LIGO Aasi  _et al._ (2015) and Advanced
Virgo Acernese  _et al._ (2015) interferometers Abbott  _et al._ (2017a). No
gravitational-wave signal from the merger or post-merger phase was observed
Abbott  _et al._ (2017d, 2019a). Approximately 1.7 s after the inferred merger
time, GRB 170817A was observed by the Fermi Gamma-ray Burst Monitor von
Kienlin  _et al._ (2017); Meegan  _et al._ (2009) and Integral Savchenko  _et
al._ (2017), with an optical/UV counterpart detected by an array of
instruments less than eleven hours later Abbott  _et al._ (2017b); Arcavi  _et
al._ (2017); Coulter  _et al._ (2017); Lipunov  _et al._ (2017); Soares-Santos
_et al._ (2017); Tanvir  _et al._ (2017); Valenti  _et al._ (2017). Subsequent
observations across a majority of the electromagnetic spectrum have continued
for more than 1000 days (e.g., Fong  _et al._ , 2019; Hajela  _et al._ , 2019,
2020; Troja  _et al._ , 2020).

Gravitational-wave observations of the inspiral phase of binary neutron star
mergers such as GW170817 and the more-recent GW190425 Abbott  _et al._ (2020)
provide valuable insight into the progenitor neutron stars, including their
masses and spins, as well as the cold equation of state of nuclear matter
Abbott  _et al._ (2017a, c, 2019b, 2020). The lack of gravitational-wave
observations of the merger and post-merger phase limits our inference of the
remnant’s evolution. We rely instead on indirect observations of the post-
merger remnant derived from electromagnetic observations of ejected and
stripped material. Understanding the post-merger evolution has the potential
to provide valuable, complementary insights into the hot nuclear equation of
state, as well as details about short gamma-ray bursts and kilonovae hitherto
unknown.

This review is dedicated to understanding what we know about binary neutron
star post-merger remnants from both an observational and theoretical
perspective, and what we hope to learn in the near future as both
gravitational-wave and electromagnetic observations increase in both number
and detail.

---
Figure 1: The fate of binary neutron star merger remnants. Two neutron stars
coalesce, losing orbital angular momentum through the emission of
gravitational waves until they eventually merge (panels A$\rightarrow$B).
Depending on the mass of the remnant, it will either promptly collapse to form
a black hole with an accretion torus and jet (panels B$\rightarrow$C), or form
a rapidly differentially-rotating neutron star (panels $B\rightarrow$D).
Depending on the mass of this neutron star, it will either be hypermassive, in
which case it will collapse to form a black hole in
$\mathcal{O}(1\,\mathrm{s})$ (panels D$\rightarrow$E), it will be
supramassive, collapsing to form a black hole in $\lesssim 10^{5}\,\mathrm{s}$
(panels F$\rightarrow$G), or it will form an infinitely stable neutron star
(panels F$\rightarrow$H).

There are four possible evolutionary pathways for a neutron star post-merger
remnant. These depend primarily on the remnant mass and the unknown neutron
star equation of state. The latter dictates the Tolman-Oppenheimer-Volkoff
mass $M_{\rm TOV}$, which is the maximum mass a non-rotating neutron star can
sustain Tolman (1939); Oppenheimer and Volkoff (1939). Observations of pulsars
in binary systems indicate $M_{\rm TOV}\gtrsim 2.0\,\mathrm{M_{\odot}}$
Demorest  _et al._ (2010); Antoniadis  _et al._ (2013); Cromartie  _et al._
(2019). Determining the evolutionary pathway of both individual binary neutron
star mergers and populations will therefore provide insights into the nuclear
equation of state.

Given a remnant mass $M$, the four evolutionary pathways (shown schematically
in Fig. 1) are:

  * •

$M\gtrsim 1.5\,M_{\rm TOV}$: the system promptly collapses to a black hole.
See path A$\rightarrow$B$\rightarrow$C of Fig. 1 and Sec. 3.

  * •

$1.2\,M_{\rm TOV}\lesssim M\lesssim 1.5\,M_{\rm TOV}$: a hypermassive neutron
star survives the collision, but collapses to form a black hole on dynamical
timescales. See path A$\rightarrow$B$\rightarrow$D$\rightarrow$E of Fig. 1 and
Sec. 4.

  * •

$M_{\rm TOV}<M\lesssim 1.2\,M_{\rm TOV}$: a supramassive neutron star will
survive the collision and will collapse to form a black hole on secular
timescales. See path A$\rightarrow$B$\rightarrow$D$\rightarrow$F$\rightarrow$G
of Fig. 1 and Sec. 5.

  * •

$M\leq M_{\rm TOV}$: a stable neutron star will survive the merger. See path
A$\rightarrow$B$\rightarrow$D$\rightarrow$F$\rightarrow$H of Fig. 1 and Sec. 5

Neutron stars can sustain more mass than their non-rotating limit $M_{\rm
TOV}$ only when rapidly rotating or extremely hot, in which case either
centrifugal support or thermal gradients provide an extra term in the force-
balance equation of hydrodynamic equilibrium. As the star spins down and/or
cools, this extra support is lost and the star eventually reaches a point
where it can no longer support its own mass and collapses to form a black
hole. In the case of hypermassive stars where $M\gtrsim 1.2M_{\rm TOV}$,
uniform rotation cannot provide enough centrifugal support to prevent
collapse, implying the star collapses as soon as enough differential rotation
is quenched. While this necessarily happens on the system’s dynamical
timescale, the exact timescale is unknown. We discuss this in detail in Sec.
4.

In the supramassive case, even once differential rotation ceases and the star
is uniformly rotating, it can still have enough centrifugal support to prevent
gravitational collapse. Secular timescales associated with magnetic dipole
radiation and gravitational-wave emission become relevant to establish the
timescale for collapse in this case. It was previously believed that collapse
would necessarily happen on a timescale of $\lesssim 10^{5}$ s Ravi and Lasky
(2014), but this is dependent on the strength of the external dipole magnetic
field. Recent afterglow observations of GW170817 may controversially shed new
light on this topic (e.g Yu  _et al._ , 2018; Li  _et al._ , 2018; Piro  _et
al._ , 2019; Ai  _et al._ , 2020; Troja  _et al._ , 2020). We discuss this in
detail in Sec. 2.1.

Each of the pathways shown in Fig. 1 has different gravitational-wave and
electromagnetic signatures, providing hope that one will be able to use such
observations to make measurements of the nuclear equation of state. For
example, the specific kilonova color depends on the survival time of merger
remnants (e.g., Li and Paczyński, 1998; Metzger  _et al._ , 2010; Metzger and
Fernández, 2014). In this article, we review theoretical and observational
aspects of each of the pathways shown in Fig. 1.

The article is set out as follows. In Sec. 2, we review observational features
of GW170817/GRB170817A that potentially allow us to discriminate the post-
merger evolutionary pathway. We detail conflicting reports that independently
suggest either a short- or long-lived neutron star survived the merger. There
is compelling evidence that we have observed electromagnetic emission from
numerous other binary neutron star mergers seen as short gamma-ray bursts. In
Sec. 2.2, we review observational features of short gamma-ray bursts that
potentially hint at long-lived neutron star remnants. Following the
observational review, we discuss more theoretical aspects of post-merger
remnant dynamics, separating the discussion into the different pathways
outlined in Fig. 1. In Sec. 3 we discuss the prompt formation of black holes,
in Sec. 4 we discuss dynamics and evolution of short-lived hypermassive
remnants, and in Sec. 5 we discuss the evolution of longer-lived supramassive
and stable neutron star remnants.

##  2 Observational evidence for post-merger remnants

###  2.1 The fate of GW170817

The smoking-gun observation to determine the nature of a post-merger remnant
are gravitational waves from the hot, differentially rotating nascent neutron
star. Searches for gravitational waves from possible post-merger remnants of
GW170817 or GW190425 have not detected a signal (Abbott  _et al._ , 2017d,
2019a, 2020)111van Putten and Della Valle (2019) contentiously claim a
detection of gravitational waves following GW170817, although see Oliver  _et
al._ (2019a) for a thorough rebuttal of this work.. This lack of detection was
expected given current detector sensitivities and theoretical models (e.g.,
Clark  _et al._ , 2016a; Sarin  _et al._ , 2018, and references therein). In
the following, we concentrate on GW170817 for three primary reasons: first,
given GW170817’s relative proximity and loudness compared to GW190425; second
the fact that the former has a plethora of electromagnetic observations,
whereas the latter had no counterpart detections; and third, because of the
larger total mass for GW190425, implying it likely promptly collapsed to form
a black hole (i.e., path A$\rightarrow$B$\rightarrow$C in Fig. 1) Abbott  _et
al._ (2020).

Although the non-detection of gravitational waves means we are unable to
definitively confirm the fate of the post-merger remnant of GW170817, much can
be inferred through the various electromagnetic observations, albeit with
somewhat conflicting conclusions. Here we elaborate on the possible fates of
the post-merger remnant of GW170817, the observations that support and
contradict each scenario.

Parameter estimation of the gravitational-wave inspiral signal constrained the
total mass of the system to $2.74^{+0.04}_{-0.01}M_{\odot}$ Abbott  _et al._
(2019b). A small fraction of this total mass $\approx 0.07\,M_{\odot}$ was
ejected and powered the optical kilonova AT2017gfo Cowperthwaite  _et al._
(2017); Smartt  _et al._ (2017); Coulter  _et al._ (2017); Soares-Santos  _et
al._ (2017); Arcavi  _et al._ (2017); Smartt  _et al._ (2017); Chornock  _et
al._ (2017); Abbott  _et al._ (2017b). AT2017gfo was first detected as a
luminous blue source with a thermal spectrum peaking in optical and
ultraviolet frequencies (e.g., Evans  _et al._ , 2017), evolving over the
course of a few days to become dominated by emission in the near-infrared
(e.g., Tanvir  _et al._ , 2017; Pian  _et al._ , 2017). The late-time near-
infrared observations agree well with “red” kilonova models predicted by the
radioactive decay of r-process nuclei (Li and Paczyński, 1998; Metzger  _et
al._ , 2010). Similarly, the early-time “blue” observations are well explained
by lower-opacity radioactive material (Metzger  _et al._ , 2010) as would be
expected if the outer layers of the ejecta are composed exclusively of light
r-process nuclei formed from matter with high electron fractions (Metzger and
Fernández, 2014). We point the interested reader to Refs. Fernández and
Metzger (2016); Metzger (2017) for detailed reviews of kilonovae and Sec. 4.2,
where we discuss the theoretical impact of merger remnants on kilonovae.

Different neutron-richness of the material implies different ejecta sources.
Two possible sources are the dynamical ejecta launched by tidal forces
(Rosswog  _et al._ , 1999; Radice  _et al._ , 2016), and matter launched from
shock heating at the contact boundary of the merger (Bauswein  _et al._ ,
2013; Hotokezaka  _et al._ , 2013; Margalit and Metzger, 2017). The former is
ejected along the equatorial plane with typically lower electron fractions
than the latter, which is typically launched in the polar direction (e.g.,
Sekiguchi  _et al._ , 2016). Another source of ejecta are outflows from the
accretion torus of a central engine (Metzger  _et al._ , 2008; Siegel and
Ciolfi, 2016; Siegel and Metzger, 2018), typically with broad electron-
fraction distributions that can increase with time due to constant neutrino
irradiation from a central neutron-star engine.

As mentioned, observations of AT2017gfo showed evidence for at least two
distinctive components: an early-time “blue” component and a late-time “red”
component. There are also hints at a third “purple” component Villar  _et al._
(2017), a point we discuss in more detail below. Common interpretations of
these observations suggest that the two dominant components of the kilonova
were powered by two distinct ejecta components. The early blue component by a
lanthanide-poor $\approx 0.02M_{\odot}$ of material with electron fraction
$Y_{e}\gtrsim 0.25$ from accretion-disc outflows along the binary polar axis
(e.g., Smartt  _et al._ , 2017; Evans  _et al._ , 2017; Tanvir  _et al._ ,
2017; Pian  _et al._ , 2017), or alternatively from dynamical ejecta launched
through shock heating. On the other hand, the late-time red component was
likely powered by lanthanide-rich $\approx 0.05M_{\odot}$ ejecta with electron
fraction $Y_{e}\lesssim 0.25$ (e.g., Cowperthwaite  _et al._ , 2017). The
total amount of ejecta and the blue component in particular, offers the first
clue into the nature of the post-merger remnant. The amount of ejecta required
to produce these observations is incompatible for a remnant that promptly
collapsed into a black hole (e.g., Radice  _et al._ , 2018a). Prompt collapse
would have resulted in a primarily red and dimmer kilonova (e.g., Margalit and
Metzger, 2017; Piro  _et al._ , 2019).

The merger was accompanied by a gamma-ray burst jet that was most-likely
structured and off axis (Alexander  _et al._ , 2018; Mooley  _et al._ , 2018a,
b). In classical gamma-ray burst models, the jet is launched through accretion
onto a black hole, leading to interpretations of the $1.7\,\mathrm{s}$ delay
between the gravitational-wave signal and gamma-ray burst (Goldstein  _et al._
, 2017; Savchenko  _et al._ , 2017) as the maximum lifetime of a putative
neutron star before it collapses into a black hole (e.g., Metzger  _et al._ ,
2018). This is further necessitated by claims that the region around the poles
needs to be relatively free of ejecta to efficiently launch an ultra-
relativistic jet (e.g., Ciolfi  _et al._ , 2019), a point we return to in Sec.
5.1. This interpretation is contentious. There are numerous short and long
gamma-ray burst observations with evidence of neutron-star central engines
(e.g., Rowlinson  _et al._ , 2010, 2013; Lü  _et al._ , 2015; Sarin  _et al._
, 2020a), providing observational evidence jets are not only launched from
accretion tori around black holes. In addition, numerical-relativity
simulations show jets could potentially be launched given sufficiently large,
but realistic, magnetic field strengths of the remnant neutron star ($B\gtrsim
10^{14}\,\mathrm{G}$) (e.g., Giacomazzo and Perna, 2013; Ciolfi  _et al._ ,
2017; Mösta  _et al._ , 2020). We discuss these points in greater detail in
Sec. 5.1. These observations and simulations disfavour the hypothesis that the
remnant of GW170817 must have collapsed into a black hole in order to launch
the jet, opening the possibility for the remnant to be a long-lived
supramassive or an infinitely stable neutron star.

With the above caveats in mind, if the remnant collapsed to a black hole
before launching the jet, strong constraints on the non-rotating neutron star
maximum mass can be derived, indicating $M_{\rm TOV}\lesssim
2.2\,\mathrm{M_{\odot}}$ (e.g., Margalit and Metzger, 2017; Rezzolla  _et al._
, 2018; Ai  _et al._ , 2020). Although see Shibata  _et al._ (2019) who derive
a more conservative estimate $M_{\rm TOV}\lesssim 2.3\,\mathrm{M_{\odot}}$ by
relaxing some of the assumptions on energy emitted in gravitational-wave and
neutrinos indirectly imposed by other analyses.

While the blue color and total ejecta mass are helpful in ruling out prompt
collapse, the exact source of the ejecta mass is unclear. Observations suggest
an ejecta mass of $\approx 0.02M_{\odot}$ with a mean velocity and electron
fraction of $\approx 0.25c$ and $Y_{e}\gtrsim 0.25$, respectively (Smartt  _et
al._ , 2017; Evans  _et al._ , 2017; Tanvir  _et al._ , 2017; Pian  _et al._ ,
2017). Although the velocity and high electron fraction agree with a shock-
heated dynamical ejecta source, the quantity of material is difficult to
explain. Relativistic hydrodynamics simulations of equal-mass progenitor
systems typically only show ejecta mass $\gtrsim 0.01M_{\odot}$ for soft
equations of state that have small neutron star radii $R\lesssim
11\,\mathrm{km}$ (Bauswein  _et al._ , 2013; Nicholl  _et al._ , 2017a;
Metzger  _et al._ , 2018). However, such compact progenitors have less ejecta
in the tidal tails, and hence less mass in the resultant torus, which is
inconsistent with the observations of the red component of the kilonova
discussed above Radice  _et al._ (2018a); Metzger  _et al._ (2018). This
raises doubt about the dynamical source of the blue component of the kilonova.
An alternate explanation for the origin of this ejecta posits that it is a
magnetized neutrino-irradiated wind from a hypermassive neutron star that
survived $\sim 1\,\mathrm{s}$ before collapsing into a black hole and had a
strong poloidal magnetic field $B_{p}\approx 1-3\times 10^{14}\,\mathrm{G}$
Metzger  _et al._ (2018). Kilonova observations a few hours after the merger,
had they existed, could have provided observational support for this
hypothesis.

There are other interpretations of the kilonova observations that imply a
different fate of the post-merger remnant. For example, there is speculation
that observations are best fit by a three-component model Villar  _et al._
(2017). Here, the early-time blue kilonova is of similar mass as inferred by
other groups, but the late-time observations are dominated by an intermediate-
opacity purple component with a significantly weaker red component. These
purple/red components may be sourced by the accretion disk around a central
engine, however it is difficult to interpret why there is a large variation in
the opacity of these components, especially if they are coming from the same
source. The purple component may be naturally supported by a long-lived
remnant neutron star Yu  _et al._ (2018), where the high-energy emission from
the remnant’s wind ionizes the surrounding material (Metzger and Fernández,
2014). Li  _et al._ (2018) compared models that could account for both the
peak luminosity and time of the kilonova observation, concluding that the
observations are best fit with a long-lived neutron star.

A long-lived remnant supports the low-significance x-ray flare
$155\,\mathrm{d}$ following the merger (Piro  _et al._ , 2019), potentially
due to untwisting toroidal magnetic field similar to x-ray flares from older
magnetars (Thompson and Duncan, 1995, 1996; Piro  _et al._ , 2019) (although,
see Ref. Lin  _et al._ (2019) for an alternate explanation). Energetic
arguments imply the toroidal component of the magnetic field must be
$B_{t}\gtrsim 10^{14}\,\mathrm{G}$, albeit with a relatively low $\approx
10^{12}\,\mathrm{G}$ poloidal field such that gravitational-wave emission
dominates the spin down for long times Piro  _et al._ (2019). Such a low
poloidal but high toroidal magnetic field structure is perhaps concerning but
not dissimilar to the soft gamma-ray repeater SGR 0418+5729 (Rea  _et al._ ,
2013) (although one should treat such measurements cautiously; see Ref.
(Mastrano  _et al._ , 2013)). Prolonged x-ray excess almost 1,000 d after the
merger is still consistent with a long-lived central engine driving the
emission Troja  _et al._ (2020). If indeed this excess is due to a long-lived
neutron star, it is almost certainly not supramassive, but likely has a mass
below the TOV mass, i.e., it follows path
A$\rightarrow$B$\rightarrow$D$\rightarrow$F$\rightarrow$H in Fig. 1. Although
the above interpretations suggest a long-lived neutron star may have formed in
the aftermath of GW170817, there are also problems with this interpretation.
For example, if long-lived, the surface poloidal component of the magnetic
field must be $B_{p}\lesssim 10^{11}-10^{12}\,\mathrm{G}$ (Yu  _et al._ ,
2018; Piro  _et al._ , 2019; Ai  _et al._ , 2018). This constraint is
problematic as it is energetically difficult to launch a Poynting-flux-
dominated jet with $B_{p}\lesssim 10^{14}\,\mathrm{G}$ (e.g., Ciolfi  _et al._
, 2017, and references therein). Furthermore, the lack of signature of the
rotational energy of the stable neutron star on the kilonova Margalit and
Metzger (2017) and radio (e.g., Ricci  _et al._ , 2020; Schroeder  _et al._ ,
2020) implies a considerable amount of energy emitted in gravitational waves,
which requires a large ellipticity (Ai  _et al._ , 2020). If the remnant of
GW170817 was long-lived, it perhaps provides the most interesting constraints
on the maximum mass of neutron stars $M_{\rm TOV}\gtrsim
2.4\,\mathrm{M_{\odot}}$ (Ai  _et al._ , 2020).

If, dear reader, you are not yet confused enough about the nature of the
remnant of GW170817, there is one final scenario consistent with all
observations. Namely, the remnant was a supramassive neutron star that spun-
down primarily through gravitational waves, collapsing into a black hole after
losing centrifugal support approximately 300 s after the merger (e.g., Ai  _et
al._ , 2020). Although one can impose longer collapse times for different
magnetic field configurations. Such a scenario supports the kilonova
observations and the potential lack of observational signature of a rapidly
spinning neutron star in other electromagnetic bands at later times. In
general, one may expect to see the signature of such a remnant on the x-ray
afterglow of the short gamma-ray burst (e.g., Rowlinson  _et al._ , 2010,
2013, and Sec. 2.2). However, Swift did not observe the region until $\sim
0.039\,\mathrm{d}$ after the gamma-ray burst trigger Evans  _et al._ (2017),
placing an upper limit on the collapse time $t_{\rm col}\lesssim
0.039\,\mathrm{d}$ of such a supramassive neutron star. If supramassive and
collapsing in less than $\sim 0.039\,\mathrm{d}$, the maximum neutron star
mass $2.1\lesssim M_{\rm TOV}/M_{\odot}\lesssim 2.4$ Ai  _et al._ (2020).

The first multimessenger binary neutron star merger GW170817 offered an
unprecedented opportunity for a detailed study into the the aftermath of such
a collision. Unfortunately, while electromagnetic observations were plentiful
and extensive, they remain inconclusive. The lack of a smoking-gun observation
of gravitational waves from a post-merger remnant make inferring the nature of
the remnant difficult. In practice, the only scenario everyone seems to agree
can be ruled out is of prompt collapse into a black hole.

###  2.2 The fate of other binary neutron star merger remnants

The coincident detection of gravitational waves from a binary neutron star
merger GW170817 and the short gamma-ray burst GRB170817A confirmed that at
least some of the latter are caused by the former. While sensitivity
improvements in gravitational-wave detectors will see increased numbers and
regularity of binary neutron star merger detections, it will remain true that
most observed gravitational-wave signals will not be accompanied with
electromagnetic signatures. Likewise, the foreseeable future will see most
short gamma-ray burst observations not accompanied by gravitational-wave
detections. But there is already a wealth of short gamma-ray burst
observational data at our disposal that can, and is, used to understand the
remnants of binary neutron star mergers.

The x-ray afterglows of some short gamma-ray bursts exhibit extended plateaus
that indicate energy injection from rapidly rotating neutron star central
engines (e.g., Rowlinson  _et al._ , 2013; Lü  _et al._ , 2015). This even
includes the observations of a putative off-axis gamma-ray burst seen only as
an x-ray transient CDF-S XT2 that is consistent with these other x-ray
afterglows Xue  _et al._ (2019); Xiao  _et al._ (2019). Such observations are
difficult to interpret in the standard afterglow model of synchrotron
radiation from shocks produced by a jet interacting with the surrounding
interstellar medium. The first such observation consistent with energy
injection from a short gamma-ray burst was GRB051221A Fan and Xu (2006), which
was followed by a catalogue of bursts Rowlinson  _et al._ (2013). These were
shown to be consistent with models of energy injection where the spin down of
the nascent neutron star is driven by magnetic dipole radiation Dai and Lu
(1998); Zhang and Mészáros (2001), a model that has further been extended to
include spin down with arbitrary braking indices Lasky  _et al._ (2017), akin
to what is seen in the spin down of radio pulsars.

For the majority of short gamma-ray bursts with extended x-ray plateaus,
debate continues to rage about the origin of the x-ray flux. Case in point,
evidence for an achromatic jet break in the broadband observations of
GRB140903A. This achromatic break has been used to argue that the long-lived
emission is due to a combination of jet geometry and dynamics of the fireball
Troja  _et al._ (2016), in stark contrast to other works that showed the
afterglow is consistent with the spin down of a long-lived central engine
Lasky  _et al._ (2017); Zhang  _et al._ (2017). Systematic Bayesian model
comparison between the two scenarios overwhelmingly favours the latter
explanation Sarin  _et al._ (2019).

Model comparison between fireball dynamics and a long-lived central engine for
another gamma ray-burst GRB130603B yields intriguingly different results Sarin
_et al._ (2019). The discriminator between the two models is a quantity called
the odds ratio, but this itself depends on the unknown equation of state. The
punch line is that, if the maximum neutron star mass is $M_{\rm TOV}\lesssim
2.3\,\mathrm{M_{\odot}}$, the data favours the fireball-shock model.
Conversely, if $M_{\rm TOV}\gtrsim 2.3\,\mathrm{M_{\odot}}$, the data favours
the existence of a long-lived neutron star central engine. It is worth
mentioning the above conclusion relies on knowing the underlying binary
neutron star mass distribution, which in light of GW190425, we do not.
Moreover, more detailed modeling for each scenario is required, which may
further discriminate and potentially yield different conclusions.

The evolution of the x-ray luminosity in these afterglows does allow us to
understand somewhat the dynamical evolution of the central engine. For
example, gamma-ray bursts GRB130603B and GRB140903A, are spinning down with
braking index $n=2.9\pm 0.1$ and $n=2.6\pm 0.1$, respectively Lasky  _et al._
(2017), where a braking index of $n=3$ is dipole magnetic spin down in vacuum.
It is worth stressing that all but one radio pulsar with accurately measured
braking index falls below $n=3$ where magnetic torques are believed to
dominate spindown Archibald  _et al._ (2016); Clark  _et al._ (2016b);
Marshall  _et al._ (2016), although see Ref. Parthasarathy  _et al._ (2019)
for a census of the highly uncertain nature of braking indices in young
pulsars. Moreover, it is also worth emphasising that calculations of realistic
braking indices for pulsars predicts they fall below 3 (e.g., Melatos,
1997)—we discuss theoretical expectations for the dynamics of neutron-star
spin down on these relatively long timescales in Sec. 5.

The highly dynamic nature of the newly-born neutron star implies the assumed
constant braking index and smooth spin-down evolution are almost certainly
naive assumptions. For example, the magnetic field’s inclination angle will
likely evolve as a function of time, leading to a changing inferred braking
index (see Şasmaz Mus  _et al._ (2019) for interpretations of long gamma-ray
bursts in this context and Sec. 5 for details of the relevant physics), or the
triaxial nature of the remnant may cause precession and a flux-modulated
amplitude of the light curve Melatos (2000); Suvorov and Kokkotas (2020). In
addition, the radiative efficiency is likely not constant throughout the spin
down and may, for example, depend on the luminosity of the central engine
itself Xiao and Dai (2019). This radiative efficiency likely changes as the
shock front decelerates as it ploughs into the interstellar medium Cohen and
Piran (1999); Dall’Osso  _et al._ (2011); Stratta  _et al._ (2018); Sarin  _et
al._ (2020b), and may evolve through plerionic-like emission as electrons fill
the cavity within the gamma-ray burst blast wave Strang and Melatos (2019). In
reality, the dynamical evolution of the remnant and the resultant emission
that eventually reaches the observer will be affected by all of these physical
processes and more, although which are truly the most dominant is still an
open question.

In all, approximately 70% of short gamma-ray bursts exhibit behaviour
consistent with long-lived remnants, split into $\approx 30\%$ with
supramassive, and $\approx 30\%$ stable neutron stars Gao  _et al._ (2016),
although these numbers are highly uncertain (e.g., Margalit and Metzger,
2019). In principle, understanding these fractions together with the
progenitor mass distribution has the ability to provide strong constraints on
the neutron star equation of state through the maximum mass. In practice, our
understanding of the numbers are not yet mature enough to make quantifiably-
reliable estimates for three reasons. First, determining the nature of the
remnant from the x-ray data alone is fraught with difficulties, as highlighted
by the various analyses of GRB130603B and GRB140903A mentioned above. Second,
while we previously thought we understood the mass distribution of binary
neutron stars from galactic radio observations of double neutron star systems
(e.g., Kiziltan  _et al._ , 2013; Alsing  _et al._ , 2018; Farrow  _et al._ ,
2019), only one of the two extragalactic neutron star mergers has progenitor
masses consistent with that distribution Abbott  _et al._ (2020). Third, it is
possible that some short gamma-ray bursts may be misidentified as long gamma-
ray bursts caused by the collapse of massive stars, or that they are caused
not by merging neutron stars, but by a neutron star-black hole merger, or that
they represent a biased sample of neutron star mergers only including ones
that produced a black hole which could launch an ultra-relativistic jet. All
of these effects would cause systematic problems with measuring the maximum
mass. Understanding both these systematic effects by collecting more
gravitational-wave observations to ameliorate the former issue, and more x-ray
plateau observations the latter, will eventually provide interesting and
stringent constraints on $M_{\rm TOV}$.

Perhaps our best understanding of post-merger behaviour comes from a subset of
eighteen short gamma-ray bursts that not only exhibit x-ray plateaus, but also
show sudden drops in the x-ray luminosity tens to thousands of seconds after
the prompt emission (e.g., Rowlinson  _et al._ , 2013; Sarin  _et al._ ,
2020a). Such dramatic changes in flux are particularly difficult to explain in
the standard fireball-shock scenario, but fit well the premise that a
supramassive neutron star was born in a neutron star merger and collapses to
form a black hole simultaneously shutting of the energy injection.

The supramassive neutron star observations again provide a tantalising way of
developing our understanding of the dynamics of the nascent neutron star and
the equation of state of nuclear matter (e.g., Fan  _et al._ , 2013; Lasky
_et al._ , 2014; Ravi and Lasky, 2014; Li  _et al._ , 2016; Gao  _et al._ ,
2016; Drago  _et al._ , 2016; Drago and Pagliara, 2018). The procedure is
straight forward: if we understand the progenitor mass distribution (which we
do not), as well as the dominant spin down mechanism (we do not understand
that either), and the spin-down rate/braking index (not really), then we can
rearrange the set of equations governing the system’s evolution to find that
the time of collapse is a function of the unknown maximum neutron star mass,
which we can therefore infer. This procedure has been performed a number of
times in different works, each arriving at different answers depending on the
underlying assumptions at each of the step. The vanilla assumptions of dipole
vacuum spin down of hadronic stars does not well fit the data Fan  _et al._
(2013); Ravi and Lasky (2014), leading some authors to infer that quark stars,
rather than hadronic stars, best explain the data (e.g., Li  _et al._ , 2016;
Drago  _et al._ , 2016), while others infer that gravitational radiation
dominates the star’s angular momentum loss rather than magnetic dipole
radiation (e.g Fan  _et al._ , 2013; Gao  _et al._ , 2016).

The correct way to do the above procedure rigorously is to include all sources
of uncertainty and marginalise over the unknown parameters such as those
describing the progenitor mass distribution, the braking index, the equation
of state, etc (for details, see Sarin  _et al._ , 2020a). Hierarchical
Bayesian inference then allows posterior probability distributions to be
calculated for each of these parameters, including those that describe the
population as a whole, rather than individual gamma-ray burst afterglows. As
one would expect, including all uncertainties implies less-well constrained
parameters. Ultimately, the eighteen-known short gamma-ray bursts allow us to
constrain $M_{\rm TOV}=2.31^{+0.36}_{-0.21}\,\mathrm{M_{\odot}}$, with 68%
uncertainties. Perhaps more interestingly, $69^{+21}_{-39}$% of remnants are
inferred to be spinning down predominantly through gravitational-wave
emission, potentially providing interesting consequences for gravitational-
wave detection of a post-merger remnant, or indirectly through a stochastic
gravitational-wave background. Furthermore, the observations show tentative
evidence for these neutron stars to be composed of deconfined quark matter,
suggesting a phase transition in the course of merger that may be visible
through gravitational-wave measurements of the inspiral (e.g., Chatziioannou
_et al._ , 2017; Bauswein  _et al._ , 2019).

Having reviewed the observational aspects of short gamma-ray bursts and what
can be learned about the remnants of the binary neutron star mergers that
power them, we now move onto more theoretical aspects, following the
evolutionary scenarios discussed alongside Fig. 1d.

##  3 Prompt formation of black holes

---
Figure 2: A post-merger remnant of mass $M\gtrsim 1.5\,M_{\rm TOV}$ will
immediately collapse to form a black hole with an accretion torus and jet.

Perhaps the least interesting of outcomes of a binary neutron star merger is
the prompt collapse to a black hole (Fig 2; panels B$\rightarrow$C). Upon
collapse, a thick accretion torus forms in the black hole’s equatorial plane,
potentially driving a Blandford-Znajek jet Blandford and Znajek (1977) that is
seen as the short gamma-ray burst. Such a scenario likely leaves an
undetectable gravitational-wave signal, and little direct electromagnetic
signal that can be used to infer properties of the central engine.

Prompt black hole formation implies that, following coalescence, the
gravitational-wave signal simply shuts down, with nearly spherical collapse
generating comparatively minimal gravitational-wave emission. We can
approximate the lowest quasi-normal mode ringdown frequency for a remnant
black hole as Echeverria (1989)

$\displaystyle f_{\rm gw}^{\rm qnm}\approx
11\,\mathrm{kHz}\left(\frac{M}{3M_{\odot}}\right)^{-1}\left[1-0.63\left(1-a\right)^{3/10}\right],$
(1)

where $M$ is the remnant mass and $a$ the dimensionless spin. Faster rotating
black holes emit higher frequency gravitational-wave signals; a non-rotating
($a=0$) black hole of $M\approx 3\,M_{\odot}$ emits its lowest quasinormal
mode signal at $f_{\rm gw}^{\rm qnm}\approx 4\,\mathrm{kHz}\,$. Prompt
collapse will more typically result in a remnant with $a\sim 0.7$–$0.8$ for
which $f_{\rm gw}^{\rm qnm}\gtrsim 6\,\mathrm{kHz}$. In these regimes, the
sensitivity of current interferometers Abbott  _et al._ (2016); Aasi  _et al._
(2015); Acernese  _et al._ (2015) and even proposed future detectors (e.g.,
Miller  _et al._ , 2015; Punturo  _et al._ , 2010; Abbott  _et al._ , 2017e;
Martynov  _et al._ , 2019; Ackley  _et al._ , 2020) is not sufficient to
detect such a signal at relevant distances.

Potentially, the prompt formation of a black hole can have implications for
electromagnetic emission, in particular through a lack of sustained energy
injection into the x-ray and optical afterglow signal and in the amount of
ejecta both dynamical and from an accretion disc. For example, kilonova
observations would likely differ from those of GW170817/AT2017gfo in that they
would be primarily red due to the lack of neutrino irradiation of the ejected
material and inferred to have less mass Margalit and Metzger (2017); Piro  _et
al._ (2019). If observations can concretely say no remnant survived the
collisions, tight constraints could be placed on the equation of state of
nuclear matter. We return to this in subsequent sections.

##  4 Short-lived, hypermassive neutron stars

---
Figure 3: A post-merger remnant of mass $1.2\,M_{\rm TOV}\gtrsim M\gtrsim
1.5\,M_{\rm TOV}$ will form a hypermassive neutron star which will collapse to
a black hole on a timescale $\mathcal{O}(1\,\mathrm{s})$.

A dominant fraction of binary neutron star mergers likely form post-merger
neutron star remnants that are either hypermassive, supramassive, or long-
lived (e.g., Gao  _et al._ , 2016; Margalit and Metzger, 2019). In all these
cases, the remnant will undergo a short period of highly-dynamic activity
(Fig. 3; panels B$\rightarrow$D) before either settling down into rigid-body
rotation (Fig. 1; panels D$\rightarrow$F), or collapsing to form a black hole
(Fig. 3; panels D$\rightarrow$E). The latter scenario describes that of a
hypermassive neutron star. We review those dynamics in Sec. 4.1, and
consequences for electromagnetic and gravitational-wave observations in Secs.
4.2 and 4.3, respectively.

We note that the dynamics and emission scenarios outlined in this section are
also relevant for the early evolution of a supramassive (Fig. 1: panels
D$\rightarrow$F$\rightarrow$G) and an eternally-stable neutron star (panels
D$\rightarrow$F$\rightarrow$H).

###  4.1 Dynamics and the Collapse Time

The short-term dynamics of binary neutron star merger remnants in the first
tens of milliseconds and up to a few seconds depends on complicated physics.
Immediately following the merger, the rotational profile of the remnant and
the amount of mass in the disk depend principally on the mass ratio of the
progenitor Bauswein  _et al._ (2013); Lehner  _et al._ (2016). The hot remnant
is rotating differentially, has an incredibly strong magnetic field, several
large-amplitude oscillation modes, and cools primarily through neutrino
emission. If the mass of the remnant is above the supramassive limit, then it
will collapse within a few seconds of formation. General relativistic
magnetohydrodynamic simulations are currently inadequate to robustly and
reliably determine that lifetime for several reasons. First, not all of the
aforementioned physics is adequately resolved (Kiuchi  _et al._ , 2018, and
see below), and second, the simulations generally only last $\lesssim
50\,\mathrm{ms}$ Ciolfi  _et al._ (2017) (although see Ref. Ciolfi  _et al._
(2019) for recent, long-lived simulations lasting up to 95 ms post merger).

Differential rotation of the remnant and thermal pressure are expected to
provide extra centrifugal support to sustain the remnant above the maximum
rigidly-rotating mass limit. When either of these is sufficiently quenched,
the remnant will collapse rapidly to form a black hole. Conventional wisdom
(e.g., Baumgarte  _et al._ , 2000; Shapiro, 2000) states that differential
rotation is suppressed on an Alfvén timescale, which can be approximately
written as

$\displaystyle\tau_{\rm A}\approx
0.3\,\mathrm{s}\left(\frac{\left<B\right>}{10^{15}\,\mathrm{G}}\right)\left(\frac{M}{2M_{\odot}}\right)^{1/2}\left(\frac{R}{10\,\mathrm{km}}\right)^{-1/2},$
(2)

where $\left<B\right>$ is the volume-averaged magnetic field inside the star,
and $R$ is the star’s radius. Below we argue that magnetic-field quenching
should occur on longer a longer timescale than expressed in Eq. 2. Clearly,
not all of the differential rotation must be suppressed for the remnant to
collapse, but only the critical threshold must be reached which will depend
ultimately on a number of factors, including how much more massive the remnant
is than the supramassive-mass threshold.

The speed with which the internal magnetic field grows, its saturation
strength, and the role of magnetic-field instabilities are potentially the
most significant unknowns when considering the suppression of differential
rotation (e.g., Ferrario  _et al._ , 2015). General relativistic, three-
dimensional magnetohydrodynamic simulations of mergers show the magnetic field
at the shock interface between the two stars grows rapidly at initial times
(e.g., Giacomazzo and Perna, 2013; Kiuchi  _et al._ , 2014, 2015; Ciolfi  _et
al._ , 2017; Kiuchi  _et al._ , 2018; Ciolfi  _et al._ , 2019). This field
growth is primarily due to the Kelvin-Helmholtz instability that develops at
the shear boundary between the two stars. The highest-resolution simulations
show amplification of the average field up to $\sim 10^{3}$ times the seed
field Kiuchi  _et al._ (2018), although this is still under-resolved and
considered a lower limit on the potential of the Kelvin-Helmholtz instability
Kiuchi  _et al._ (2018); Duez and Zlochower (2019). Magnetic winding, the
magnetorotational instability, and turbulence can subsequently amplify the
field to $\gtrsim 10^{16}\,\mathrm{G}$ on short timescales Duez  _et al._
(2006a, b); Zrake and MacFadyen (2013); Siegel  _et al._ (2013); Ciolfi  _et
al._ (2019). It is worth noting that these latter effects are well-beyond the
resolution limit of current numerical simulations Baiotti and Rezzolla (2017);
Kiuchi  _et al._ (2018), implying quantitative evolution of the star’s
magnetic-field growth should not be trusted.

Most numerical-relativity simulations are performed under the assumption of
ideal magnetohydrodynamics, where infinite conductivity of the fluid implies
the magnetic field is frozen into the dynamic fluid. However, the hot remnant
should almost certainly have regions of finite conductivity, implying flux
freezing is not the correct assumption. The scarce numerical-relativity
simulations of post-merger remnants that include the effects of resistive
magnetohydrodynamics in the core indeed show the expected result that the
magnetic field lags behind the fluid Dionysopoulou  _et al._ (2015). This
implies that magnetic winding is less efficient than previously believed, and
the Alfvén timescale of Eq. 2 should be taken as a rough lower bound for the
quenching of differential rotation through magnetic-field winding.

In addition to differential rotation, the nascent star is supported through
thermal pressures. Merger simulations ubiquitously show temperatures at the
shock interface in excess of $3\times 10^{11}\,\mathrm{K}$, and the majority
of the star above $10^{11}\,\mathrm{K}$ (e.g., Sekiguchi  _et al._ , 2011a, b;
Foucart  _et al._ , 2016; Perego  _et al._ , 2019). Simulations that take
cooling into account show the remnant may not collapse for a few cooling
timescales Paschalidis  _et al._ (2012), which may be as long as a few seconds
Duez and Zlochower (2019). We return to the cooling of the nascent neutron
star in Sec. 5.3, but the above argument suggests the dominant physical
processes dictating the collapse time of the remnant is magnetic winding
rather than thermal dissipation.

Quite clearly there is complicated physics dictating the collapse time of a
hypermassive remnant. If one believes that the collapse of the hypermassive
star is required to launch the gamma-ray burst, then 1.74 s delay between the
gravitational-wave inferred merger time of GW170817 and the gamma-ray
detection could be partially explained by the collapse time of the
hypermassive star. This is complicated by the additional unknown time it takes
for the jet to launch from the black hole, and for the jet to break out from
the merger ejecta. It would therefore be nice to be able to directly measure
the collapse time of the hypermassive star, which may be done with future
gravitational-wave detections—see Sec. 4.3.

Ultimately, when the remnant collapses to form a black hole, it will do so
approximately on a dynamical timescale Stark and Piran (1985). The dynamics
and timescale of the collapse depend heavily on the system’s angular momentum
Baiotti  _et al._ (2005, 2007) and the degree of differential rotation
Giacomazzo  _et al._ (2011). Although the free-fall timescale of a massive
neutron star is $\lesssim 0.01\,\mathrm{ms}$, the collapse timescale for a
rapidly, differentially rotating neutron star is $\sim 1\,\mathrm{ms}$
Giacomazzo  _et al._ (2011). Whether the collapse time can be measured upon a
successful detection of post-merger gravitational waves is an open question
(see Sec. 4.3), however the lifetime of the remnant is believed to indirectly
impact the spectral evolution of the electromagnetic signature; a point to
which we now turn.

###  4.2 Electromagnetic consequences

Kilonova emission is a direct result of radioactive decay of heavy elements
produced by the merger ejecta. As this review deals specifically with the
merger remnant, rather than the merger itself, we do not review the physics of
kilonovae emission directly. Instead, we focus on the effects the hypermassive
remnant’s lifetime, dynamics, and evolution have on the kilonova and other
electromagnetic emission channels. For reviews of kilonova itself, see e.g.,
Refs. Fernández and Metzger (2016); Metzger (2017).

The ejecta in a neutron star merger can be broadly split into two categories:
dynamical ejecta produced in the merger itself, and the (secular) outflow from
the accretion disc formed around the remnant object. We note that there may be
additional sub-dominant channels that contribute to the total ejecta. Both the
amount of ejecta and its properties (velocity and electron fraction) are
intrinsically connected to binary parameters and the fate of the merger
remnant, with the electron fraction being perhaps the most critical as it
directly impacts what elements can be synthesised and therefore the color of
the kilonova.

Dynamical ejecta usually constitutes two sources: shock-heated ejecta from the
contact interface between the two merging neutron stars distributed
isotropically, and spiral arms from the tidal interactions in the merger which
is predominantly in the equatorial plane. The total quantity of dynamical
ejecta depends sensitively on the fate of the merger remnant and the binary
mass ratio, if the remnant promptly collapses into a black hole, there will be
little shock-heated ejecta, as the region is promptly swallowed up Bauswein
_et al._ (2013); Ciolfi  _et al._ (2017); Radice  _et al._ (2018b). In
general, asymmetric binaries tend to produce more ejecta Rezzolla  _et al._
(2010); Bauswein  _et al._ (2013). However, we emphasise that this
relationship is not wel understood quantitatively. Numerical simulations
suggest that the total dynamical ejecta in a merger is in the range of
$10^{-4}-10^{-2}M_{\odot}$ with velocities in the range $0.1-0.3$c (see e.g.,
Hotokezaka  _et al._ , 2013) and a broad electron fraction distribution,
$Y_{e}\sim 0.1-0.4$ Radice  _et al._ (2016) which dictates what elements can
be synthesised from this ejecta, potentially up to an atomic mass number,
$A\sim 195$ Wanajo  _et al._ (2014); Martin  _et al._ (2015). We note that
this is an area of active research, with significant uncertainties in many
critical nuclear reaction quantities Zhu  _et al._ (2020).

The other source of ejecta in a neutron star merger is the outflow from the
accretion disk that forms around the remnant object. The quantity of mass in
the accretion disk range from $\sim 0.01-0.3M_{\odot}$ depending on the binary
parameters (e.g., Oechslin and Janka, 2006) and the fate of the merger
remnant, with prompt formation likely resulting in less mass around the
remnant object (e.g., Perego  _et al._ , 2014; Metzger and Fernández, 2014;
Martin  _et al._ , 2015; Metzger, 2017). Depending on the lifetime of the
remnant neutron star (as we discuss below), outflows from this disc likely
contribute more mass to the ejecta than the dynamical ejecta launched in the
merger itself (Perego  _et al._ , 2014; Fernández and Metzger, 2016; Siegel
and Metzger, 2018), a statement seemingly verified by the inferred properties
of the kilonova, AT2017gfo (e.g., Smartt  _et al._ , 2017; Kasen  _et al._ ,
2017; Metzger  _et al._ , 2018). It is the properties of this outflow that are
most affected by the nature and lifetime of the remnant and make the biggest
impact on the kilonova. In particular, the cooling of the nascent neutron star
through neutrino losses. We discuss the thermal evolution of neutron stars in
detail in Sec 5.3.

The accretion disk itself evolves viscously and cools through neutrino
emission, driving a wind similar to proto-neutron stars born in supernovae
(e.g., Beloborodov, 2008; Metzger  _et al._ , 2008). The mass loss through
this channel is dependent on the neutrino flux, which as we elaborate below is
connected to the fate of the remnant. Depending on the merger outcome, this
process results in a mass loss of up to $10^{-3}M_{\odot}$ either from the
disc, remnant neutron star or a combination of the two. Further evolution of
the disc is dictated by angular momentum transport, either through turbulence
generated by the magneto-rotational instability (e.g., Metzger and Fernández,
2014; Siegel and Metzger, 2017, 2018) or by spiral density waves excited by
oscillations of a neutron star remnant which expand the disc outwards (e.g.,
Nedora  _et al._ , 2019; Metzger, 2017). Initially, the disk accretes matter
at a relatively high rate, but once this rate drops below a critical
threshold, cooling through neutrinos is ineffective and the disk thermally
expands. In this process, free nucleons recombine into $\alpha-$particles
which releases enough energy to unbind a significant fraction of the disc
(e.g., Beloborodov, 2008; Metzger  _et al._ , 2008, 2010; Fernández and
Metzger, 2013; Perego  _et al._ , 2014; Martin  _et al._ , 2015; Metzger and
Fernández, 2014; Siegel and Metzger, 2017; Fernández  _et al._ , 2019). The
amount of ejecta this process unbinds is again connected to the fate of the
remnant, with numerical simulations suggesting prompt black hole formation
unbinds up to $\sim 40\%$ of the disc (e.g., Fernández and Metzger, 2016)
while a neutron star remnant unbinds potentially up to $\sim 90\%$ (e.g.,
Siegel and Metzger, 2017) due in large part to the additional neutrino flux
from the remnant neutron star.

In principle, the quantity of ejecta driven by the outflow of the accretion
disc is closely linked to the fate of the remnant and can principle be used to
infer the fate of the remnant. However, there are substantial quantitative
uncertainties (e.g., Bernuzzi, 2020) due to simplified neutrino treatments,
numerical artifacts from limited resolution, and additional unmodelled
processes such as a magnetised neutrino driven wind (e.g., Metzger  _et al._ ,
2018). Moreover, there are significant systematic uncertainties associated
with nuclear reaction networks, opacities, etc., that can led to substantial
bias in inferring properties of the kilonova from observations Zhu  _et al._
(2020); Barnes  _et al._ (2020).

Ignoring the aforementioned complications, a more reliable discriminator of
the fate of the remnant is to consider the impact of neutrino radiation on the
electron fraction of the ejecta. As we discuss in more detail in Sec 5.3,
nascent neutron stars cool by emitting neutrinos. The additional neutrino flux
increases the electron fraction with time through $\nu_{e}+n\rightarrow
p+e^{-}$. The electron fraction continues to increase with longer remnant
lifetimes making it increasingly difficult to synthesize heavier $r$-process
elements (e.g., Metzger and Fernández, 2014; Kasen  _et al._ , 2015; Lippuner
_et al._ , 2017; Kawaguchi  _et al._ , 2020), which directly affects the
colour of the kilonova. For $Y_{e}\lesssim 0.25$, a predominantly red kilonova
is expected with elements greater than atomic mass $A\gtrsim 140$, while for
electron fractions $Y_{e}\gtrsim 0.25$, smaller elements are expected and the
colour of the kilonova is predominantly blue (e.g., Metzger, 2017). Numerical
calculations suggest that a remnant lifetime longer than $\sim 300$ ms will
make $Y_{e}\gtrsim 0.25$ in the ejecta from the disk outflow Lippuner  _et
al._ (2017), although other calculations suggest a larger lifetime of $\sim
1$s for a similar electron fraction Sekiguchi  _et al._ (2016); Kawaguchi  _et
al._ (2020).

###  4.3 Gravitational-wave emission and detection

Gravitational-wave emission from the first second post-merger is expected to
have a relatively large strain amplitude, possibly comparable to the peak
amplitude of the inspiral phase, albeit with frequencies in the kHz range.
This makes them an interesting target for current and future ground-based
gravitational-wave observatories. In the following two subsections we review
state-of-the-art predictions for gravitational-wave emission and detection,
respectively.

####  4.3.1 Gravitational-wave emission

Numerical-relativity simulations of binary neutron star mergers and their
subsequent post-merger evolution show that gravitational-wave emission from a
hypermassive remnant is dominated by the quadrupolar $f$-mode (e.g., Xing  _et
al._ , 1994; Ruffert  _et al._ , 1996). Depending on the equation of state,
this occurs anywhere from $\sim 2$ to $4\,\mathrm{kHz}$ (e.g., Takami  _et
al._ , 2015), and is strongly correlated with the star’s compactness and tidal
deformability (e.g., Bauswein  _et al._ , 2012; Bauswein and Janka, 2012;
Hotokezaka  _et al._ , 2013; Read  _et al._ , 2013; Takami  _et al._ , 2014;
Bauswein and Stergioulas, 2019). It is actually somewhat surprising that the
frequency of the dominant $f$-mode post-merger correlates so well with the
tidal deformability and compactness—these are quantities calculated for cold,
non-rotating neutron stars, whereas the hypermassive post-merger remnant is
rapidly rotating and has a temperature $\gtrsim{\rm few}\times
10^{10}\,\mathrm{K}$ (e.g., Sekiguchi  _et al._ , 2011a, b; Foucart  _et al._
, 2016). This suggests rotational and temperature effects play a minor role in
the dominant properties of the gravitational-wave signal, and broadly implies
that a successful measurement of the dominant post-merger gravitational-wave
frequency is a robust measurement of the nuclear equation of state.
Importantly, if the gravitational-wave frequency of the post-merger
oscillations do not match up to the tidal deformation measured from the
inspiral phase, this could be the signature of a first-order hadron-quark
phase transition occurring in the core of neutron stars at high temperatures
and pressures Bauswein  _et al._ (2019).

In Fig. 4 we plot an example gravitational waveform from a hypermassive post-
merger remnant Bernuzzi  _et al._ (2015)222This waveform is publicly available
through the CoRe database of binary neutron star merger waveforms (waveform ID
BAM:0035; Dietrich  _et al._ , 2018). This shows the merger of two
$1.35\,\mathrm{M_{\odot}}$ neutron stars with the H4 equation of state Lackey
_et al._ (2006) at a distance of $40\,\mathrm{Mpc}$. The strain amplitude is
maximal at time $t=0$, which we take as a proxy for the merger time (other
definitions are often used, including the first contact between the two
stars). For this simulation, the neutron star collapsed after time $t\approx
14\,\mathrm{ms}$, although we reiterate this collapse time is not a reliable
prediction for all of the reasons given in Sec. 4.1.

---
Figure 4: Gravitational-wave strain from a numerical relativity simulation of
a binary neutron star merger. Gravitational-waves from the hypermassive post-
merger remnant can have amplitudes comparable to that of the peak of the
inspiral. In this example, the gravitational-wave signal shuts off rapidly
approximately $14\,\mathrm{ms}$ after the merger, signifying the collapse of
the hypermassive neutron star to a black hole. This simulation Bernuzzi  _et
al._ (2015); Dietrich  _et al._ (2018) is of an equal-mass binary with
$M=1.35\,\mathrm{M_{\odot}}$ component masses using the H4 equation of state
at a distance of $40\,\mathrm{Mpc}$.

In Fig. 5 we plot the gravitational-wave amplitude spectral density for the
waveform shown in Fig. 4. The amplitude spectra of the full waveform, which
includes almost ten full orbits of the binary prior to merger, is shown as the
solid black curve, whereas the dashed black curve includes only the post-
merger component. For comparison, we also plot the design amplitude noise
spectral densities for three instruments; Advanced LIGO (solid blue curve;
Aasi  _et al._ , 2015), the Einstein Telescope (red dashed curve; Punturo  _et
al._ , 2010) and Cosmic Explorer (dot-dashed green curve; Abbott  _et al._ ,
2017e).

For the three sensitivity curves shown alongside the predicted gravitational-
wave spectrum in Fig. 5 we can calculate the expectation value of the single-
detector, optimal matched filter

$\displaystyle\left<{\rm S/N}\right>^{2}=4{\rm Re}\int
df\frac{\left|\tilde{h}(f)^{2}\right|}{S_{h}(f)},$ (3)

where $\tilde{h}(f)$ is the Fourier transform of the gravitational-wave time
series, and $S_{h}(f)$ is the noise power spectral density. We calculate this
for the post-merger signal only (i.e., for time $t>0$) and find $\left<{\rm
S/N}\right>=1.0$, 9.0, and 12.7, for Advanced LIGO, Einstein Telescope, and
Cosmic Explorer, respectively.

Together with the latest prediction for the merger rates derived from the
first gravitational-wave observation of a binary neutron star merger Abbott
_et al._ (2017a), one can estimate the expected event rates for post-merger
gravitational-wave detection. Marginalising over a range of equations of
state, Ref. Martynov  _et al._ (2019) predicted $\lesssim 2$, ${\rm S/N}>5$
detections of a post-merger remnant per year with Einstein Telescope, and
$\lesssim 10$ such detections per year with Cosmic Explorer (see also Clark
_et al._ , 2016a). According to these calculations, the pay-off is likely to
only come with third-generation detectors where the high-frequency
($\gtrsim$kHz) sensitivity is increased by a factor of at least ten over
Advanced LIGO design sensitivity. This has prompted many to also think about
the potential for stacking multiple sub-threshold post-merger signals with
second- or third-generation inteferometers Yang  _et al._ (2018), or to build
dedicated high-frequency gravitational-wave instruments with the primary
science goal to detect tidal effects during the inspiral of binary neutron
star mergers as well as their post-merger remnants (e.g., Miao  _et al._ ,
2018; Martynov  _et al._ , 2019; Ackley  _et al._ , 2020).

---
Figure 5: Gravitational-wave amplitude spectrum from the hypermassive post-
merger remnant shown in Fig. 4. The solid black curve shows the full spectral
density including almost $30\,\mathrm{ms}$ if inspiral, whereas the dashed
black curve shows the spectral density of just the post-merger remnant (i.e.,
truncating the time series at the merger). The solid blue, dashed red, and
dot-dashed green curves are the projected amplitude spectral densities of
Advanced LIGO at design sensitivity, the Einstein Telescope, and Cosmic
Explorer, respectively. This post-merger signal has a single-detector signal-
to-noise ratio of 1.0, 9.0 and 12.7 for the three detectors, respectively.

When the hypermassive star collapses to form a black hole, it does so on a
relatively short timescale ($\sim\mathrm{ms}$; see Sec. 4.1 and Fig. 4),
implying the signal is potentially upward of $\sim 1\,\mathrm{kHz}$. However,
the amplitude of the signal is weak; optimistic estimates from simulations of
collapsing, differentially-rotating stars suggest they may be detectable at a
distance of up to $10\,\mathrm{Mpc}$ with third-generation detectors such as
the Einstein Telescope or Cosmic Explorer Giacomazzo  _et al._ (2011).
Numerical relativity simulations of binary neutron star mergers with
subsequent hypermassive star formation and collapse seem to show no
discernible burst of radiation above that seen from the star’s oscillations
immediately prior to collapse.

It is currently an open question as to whether gravitational-wave parameter
estimation methods targeting hypermassive neutron star signals—see next
section—can infer the collapse time with a successful gravitational-wave
detection. If they can, this would allow us to probe the complex physics that
governs the quenching of differential rotation discussed in Sec. 4.1.

####  4.3.2 Gravitational-wave detection methods

The LIGO/Virgo collaborations have a number of methods in place to search for
and characterise gravitational waves in the immediate aftermath of a neutron-
star merger. Two algorithms were used to search for short-lived gravitational-
wave signals following GW170817 Abbott  _et al._ (2017d, 2019b): Coherent Wave
Burst (cWB; Klimenko  _et al._ , 2016) and BayesWave Cornish and Littenberg
(2015); Littenberg and Cornish (2015); Chatziioannou  _et al._ (2017). The cWB
algorithm searches for coherent excess power in multi-resolution wavelet
transformations, while BayesWave uses Bayesian inference, modelling the
gravitational-wave signal itself as a linear superposition of wavelets; see
Refs. Abbott  _et al._ (2017d) and Abbott  _et al._ (2019b), respectively, for
details of the specific implementation and setup of the two algorithms
searching for gravitational waves following GW170817. Neither methods found
any hint of a signal, but placed upper limits on the total energy emitted in
gravitational waves.

The two methods described above do not rely on waveform models to search for
the gravitational-wave signal. While they are therefore more robust than
modelled searches that use template waveforms, they are also less sensitive
Tsang  _et al._ (2019); Easter  _et al._ (2020). However, modelled searches
are in their infancy due to a paucity of enough reliable gravitational-wave
templates to perform matched-filter searches. This is rapidly changing, with
analytic approximations Bauswein  _et al._ (2016); Bose  _et al._ (2018),
principal component decompositions Clark  _et al._ (2016a), and machine-
learning algorithms Easter  _et al._ (2019) showing promising results fitting
to numerical-relativity waveforms.

Two recent Bayesian methods have been independently developed that use
analytic waveforms. Using only a single oscillation mode modelled as a damped
sinusoid (or Lorentzian function in the frequency domain), Tsang  _et al._
(2019) showed an average mismatch between numerical-relativity injections and
recovered signals of 0.15. A single damped sinusoid allows the main $f_{2}$
peak to be measured, which is the key peak for determining the equation of
state. This method is therefore capable of distinguishing inspiral and post-
merger inferences of the equation of state Tsang  _et al._ (2019) to, for
example, determine potential quark deconfinement in the stellar core that only
occurs at high temperatures Bauswein  _et al._ (2019).

Instead of a single damped sinusoid, Easter  _et al._ (2020) modelled the full
waveform as a linear sum of three damped sinusoids (inspired by Refs. Bauswein
_et al._ , 2016; Bose  _et al._ , 2018), also allowing all three frequencies
to drift linearly in time. They found an average mismatch of only 0.03,
implying the method is $\approx 15\%$ more sensitive than Ref. Tsang  _et al._
(2019). The addition of the extra mode oscillations in the analytic waveform
approximations is unlikely to improve equation of state estimates, however
realistic modelling of the frequency drift of the fundamental $f_{2}$ mode may
have some, as yet undetermined effect.

In reality, given the potential inaccuracies of gravitational waveforms from
numerical-relativity simulations (see previous section), both modelled and
unmodelled searches will be needed as this field hopefully moves from the
development to the observational stage.

Of course, many binary neutron star coalescences will likely be detected
before a bona fide post-merger detection. This opens the possibility to
effectively stack sub-threshold events by either multiplying Bayes factors
from individual events or by coherent summation of signals relying on pre-
merger phase information Yang  _et al._ (2018).

##  5 Long-lived neutron star remnants

---
Figure 6: A post-merger remnant of mass $1.0\,M_{\rm TOV}\leq M\gtrsim
1.2\,M_{\rm TOV}$ will form a supramassive neutron star, which will collapse
to a black hole on a timescale $\lesssim 10^{5}\,\mathrm{s}$
---
Figure 7: A post-merger remnant of mass $M\leq 1\,M_{\rm TOV}$ will form an
infinitely-stable neutron star.

Post-merger remnants less massive than $\sim 1.2M_{\rm TOV}$ will survive for
more than one second. They are deemed supramassive if their mass is greater
than the non-rotating limit $M_{\rm TOV}$—Fig. 6—or infinitely stable if their
mass is less than $M_{\rm TOV}$—Fig. 7. Observational evidence from the x-ray
afterglow of short gamma-ray bursts suggest a non-negligible fraction of
binary neutron star mergers result in these outcomes (e.g., Rowlinson  _et
al._ , 2013; Lü  _et al._ , 2015).

Figures 6 and 7 emphasise an important point: even supramassive and stable
remnants will still undergo a period of strong differential rotation
immediately after the merger (panel B), implying the gravitational-wave
emission and detection discussion from Secs. 4.1 and 4.3 are still relevant.
However, electromagnetic emission due to r-process nucleosynthesis discussed
in Sec. 4.2 may be different depending on the lifetime of the remnant; a point
we return to in Sec. 5.5.

###  5.1 Can a long-lived neutron star launch a jet?

The question whether a neutron star can launch an ultra-relativistic jet that
can produce a short gamma-ray burst continues to engage theorists and
observers. Before we review our more recent understanding of this question, it
is intriguing to start near the beginning and briefly review our understanding
of short gamma-ray bursts themselves. We refer the interested reader to more
detailed reviews of gamma-ray bursts (e.g., Nakar, 2007; Kumar and Zhang,
2015; Nakar, 2019).

Several lines of evidence of short gamma-ray bursts suggest that the central
engine must be able to launch an ultra-relativistic outflow (e.g., Ackermann
_et al._ , 2010), a constraint demanded by the compactness problem Ruderman
(1975). Constraints on the total energy requires that the outflow is
collimated into a narrow jet (e.g., Nakar, 2007), while rapid variability of
the prompt emission implies the engine must be a compact object (e.g., Nakar,
2007; Berger, 2014). Immediately, these constraints imply two potential
central engines, a black hole or a neutron star. Moreover, one of the implicit
requirements imposed by an ultra-relativistic outflow is having a relatively
baryon-free environment where the jet is launched (e.g., Nakar, 2007; Ciolfi,
2018), such that the ultra-relativistic jet can launch efficiently, break out
of the environment and stay relativistic. It is this constraint of having a
clean environment, sometimes refereed to as the baryon-loading problem (Shemi
and Piran, 1990), that has led to debate about whether long-lived neutron
stars can launch jets and power short gamma-ray bursts.

Assuming the central engine is a black hole, there are two main mechanisms for
generating an ultra-relativistic jet: neutrino-antineutrino annihilation along
the black hole spin axis (Ruffert and Janka, 1998), and the Blandford-Znajek
mechanism (Blandford and Znajek, 1977). Although both processes may be
responsible for launching jets, the Blandford-Znajek mechanism is more
favoured due to energetic constraints, with the neutrino-antineutrino
annihilation unable to produce high energy gamma-ray bursts (e.g., Ciolfi,
2018). In the case of a long-lived neutron star the mechanism for launching
the jet is far from certain, although it is clear that if a jet is launched,
it is likely through magnetic processes that tap into the large rotational-
energy reservoir (e.g., Bucciantini  _et al._ , 2012; Ciolfi  _et al._ ,
2019). However, the problem as we alluded to above is the baryon-load of the
environment. We refer the interested reader to more detailed reviews of gamma-
ray burst jet launching (e.g., Ciolfi, 2018).

As we discussed earlier in Sec 4.2, neutrino radiation from the nascent
neutron star unbinds a non-negligible (potentially $\gtrsim 10^{-2}M_{\odot}$)
amount of matter (e.g., Fernández and Metzger, 2016; Metzger, 2017). This
additional ejecta pollutes the surrounding environment, particularly along the
spin axis (e.g., Ciolfi  _et al._ , 2017). This pollution means that even
tapping into the entire $\approx 10^{52}$ erg of rotational energy, the
maximum Lorentz factor of a jet that forms is $\mathcal{O}(10)$, an order of
magnitude smaller than required to alleviate the compactness problem (Nakar,
2007; Murguia-Berthier  _et al._ , 2014, 2017), and two orders of magnitude
smaller than measurements from short gamma-ray bursts (e.g., Ackermann  _et
al._ , 2010). This led Murguia-Berthier  _et al._ (2017) to argue that the
formation of an ultra-relativistic jet has to wait until the formation of a
black hole. One may think this opens up the possibility of a supramassive
neutron star producing a short gamma-ray burst, however detailed analysis by
Margalit  _et al._ (2015) showed that for a range of equations of state,
accretion disc and therefore jet formation is unlikely after the collapse of a
supramassive neutron star. Given our current constraints on the maximum mass
of neutron stars, this suggests up to 80% of binary neutron star mergers will
not produce a short gamma-ray burst.

Beniamini  _et al._ (2017) showed that to efficiently launch an ultra-
relativistic jet, the energy per baryon ($\eta\propto\dot{E}/\dot{M}$) at the
base of the jet must exceed a critical threshold ($\eta\gtrsim 100$).
Following this result, one can derive a critical timescale for when this
threshold is exceeded Beniamini  _et al._ (2020). For $B\sim 3\times
10^{16}\,\mathrm{G}$ this timescale is $\approx 0.2\,\mathrm{s}$, short enough
to be consistent with the delay seen in GRB170817A. Alternatively, for more
typical parameters, this timescale is $\gtrsim 3\,\mathrm{s}$, inconsistent
with the delay of GRB170817A. This suggests that typical magnetars are unable
to produce short gamma-ray bursts with a delay consistent with GRB170817A,
although it is important to emphasise that it is not necessary for all short
gamma-ray bursts to have the same delay from merger, implying a timescale
$\gtrsim 3\,\mathrm{s}$ is not necessarily problematic. Furthermore, given the
complicated physics at play for launching a jet, it is worth questioning how
reliably these scalings can distinguish between timescales within an order of
magnitude. More importantly however is the dependency of this timescale on the
magnetic field. As discussed in earlier sections, numerical simulations do not
resolve the magneto-rotational and Kelvin-Helmholtz instabilities,
significantly underestimating the amplification of the magnetic field. If the
magnetic field amplification is higher this will naturally alleviate this
problem with the timescale for efficiently launching ultra-relativistic jets.

While this work in the past painted a pessimistic view on a long-lived neutron
star launching an ultra-relativistic jet. Recent numerical simulations have
changed the state of this research area Mösta  _et al._ (2020). Including
detailed physics of neutrinos in their simulation, Mösta  _et al._ (2020)
showed that neutrinos emitted by the nascent neutron star predominantly around
the polar region end up preventing baryon pollution and aiding formation of a
relativistic jet, in contrast to simulations which did not include neutrino
effects and had polluted polar regions (e.g., Ciolfi  _et al._ , 2017, 2019).
To confirm this result, Mösta  _et al._ (2020) performed identical simulations
without including neutrino effects resulting in a polluted polar region. These
latest simulations strongly suggest that long-lived neutron stars are viable
engines for short gamma-ray bursts.

While theoretical support is slowly arriving, there is considerable, although
perhaps subjective, observational evidence for long-lived neutron stars being
viable engines of short gamma-ray bursts. The $1.74$ s delay between the
gravitational-wave signal, GW170817 and gamma-ray burst, GRB170817A has been
calculated to be dominated by the timescale for the relativistic jet to reach
the $\gamma$-ray emitting radius (Ren  _et al._ , 2020; Beniamini  _et al._ ,
2020). Independently, kilonova observations have been suggested to require a
magnetised wind from a hypermassive neutron star that survived for up to $1$ s
(Metzger  _et al._ , 2018). In combination, these two independent arguments
suggest that the relativistic jet in GRB170817A must have been launched by the
differentially rotating neutron star before it collapsed into a black hole. If
a differentially rotating neutron star can launch such a jet then long-lived
neutron stars which also go through this evolutionary phase must be able to
launch a jet as well.

Another more tentative observational evidence for a neutron star engine is the
recent detection of a luminous kilonova from short GRB200522A (Fong  _et al._
, 2020). The kilonova is significantly brighter than expected from r-process
nucleosynthesis alone and the characteristic features of the emission suggest
additional energy from a neutron star. If true, this neutron star is most
definitely stable and given the kilonova is detected as a counterpart to a
short gamma-ray burst, the neutron star engine must have launched the jet that
powered this gamma-ray burst. We note that given there are large systematic
uncertainties with kilonova modelling (e.g., Zhu  _et al._ , 2020) there may
be alternative explanations for this exceptionally bright kilonova.

Perhaps the best observational evidence for long-lived neutron stars being
able to launch relativistic jets and power short gamma-ray bursts are x-ray
afterglow observations. The observations of a plateau (e.g., Fan and Xu, 2006;
Dall’Osso  _et al._ , 2011; Rowlinson  _et al._ , 2013; Lü  _et al._ , 2015),
a sharp drop in luminosity (e.g., Rowlinson  _et al._ , 2010, 2013; Sarin  _et
al._ , 2020a) and late-time x-ray flares (e.g., Fan and Xu, 2006) are all are
best interpreted by invoking the spin-down energy of a long-lived neutron star
(Zhang and Mészáros, 2001; Lasky  _et al._ , 2017). If long-lived neutron
stars are powering the x-ray afterglows they must by definition also be
responsible for launching the ultra-relativistic jet and producing prompt
emission.

Although these observational lines of evidence are not definitive, they do
strongly suggest long-lived neutron stars can launch jets that produce short
gamma-ray bursts, with more recent numerical simulations (e.g., Mösta  _et
al._ , 2020) backing up this claim.

###  5.2 Supramassive or stable

Back-of-the-envelope calculations suggest supramassive stars born from
neutron-star mergers are expected to survive for anywhere between $\sim
10\,\mathrm{s}$ to $\sim 4\times 10^{4}\,\mathrm{s}$ Ravi and Lasky (2014),
which is broadly consistent with observationally-inferred collapse times from
short gamma-ray burst x-ray afterglows Rowlinson  _et al._ (2013). However, a
detailed look at the expected and observed distributions show the stars tend
to collapse on shorter timescales than one might expect Fan  _et al._ (2013);
Ravi and Lasky (2014); Sarin  _et al._ (2020a). Suggestions for the resolution
of this discrepancy include excess gravitational-wave emission at early times
Fan  _et al._ (2013); Gao  _et al._ (2016); Lü  _et al._ (2017); Lin and Lu
(2019), or deconfined quark matter changing the star’s moment of inertia from
that with normal matter Li  _et al._ (2016, 2017); Drago and Pagliara (2018).
Hierarchical Bayesian inference of the population of collapse times for 18
supramassive neutron stars identified by the $72$ Swift x-ray afterglow
observations suggest the story may be some combination of these two affects
Sarin  _et al._ (2020a).

These constraints have important consequences. Firstly, if these nascent
neutron stars are really composed of deconfined quarks, then it suggests there
is a temperature dependent hadron-quark phase transition. Understanding where
in the nuclear phase diagram such a transition occurs is critical to our
understanding of the behaviour of nuclear matter (e.g., Bauswein  _et al._ ,
2019; Chatziioannou and Han, 2020). An unexplored consequence of this phase
transition is on the kilonova; additional neutrino flux from the cooling of a
hybrid star cf. a purely hadronic star may mean that the electron fraction in
the accretion disk outflow is higher than expected. In turn, this could make
it difficult to infer the time the remnant collapsed based on the r-process
elements synthesised (e.g., see Lippuner  _et al._ , 2017; Kawaguchi  _et al._
, 2020). We discuss the cooling of nascent neutron stars in more detail in
Sec. 5.3.

The implication that supramassive neutron stars are spinning down
predominantly through gravitational waves has important consequences for the
dynamics of these neutron stars, a point we discuss in more detail in later
sections. More immediately relevant is the implication significant
gravitational-wave spin down has on the electromagnetic observations. As
discussed in Sec 2.1, energetic constraints on GW170817 from kilonova (e.g.,
Margalit and Metzger, 2017) and radio (e.g., Ricci  _et al._ , 2020; Schroeder
_et al._ , 2020) observations put strong limits on the total energy of the
ejecta. It has been proposed that this energy constraint indirectly rules out
a long-lived neutron star in GW170817 Margalit and Metzger (2017). However, if
the supramassive star spins down predominantly through gravitational-wave
emission and therefore collapses earlier than expected, then these constraints
are weakened significantly. This also implies that future analyses utilising
energy constraints to infer the nature of the remnant need to carefully
consider energy losses through gravitational waves. However, it is worth
stressing we have an incomplete, quantitative understanding of how much energy
is radiated in gravitational waves.

Finally, it is worth mentioning that the hierarchical analysis of collapse
times Sarin  _et al._ (2020a) predict a large fraction of supramassive stars
collapse in $\lesssim 100\,\mathrm{s}$. This is significant, because this is
also approximately the time it takes Swift to slew, implying a number of post-
merger remnants may be observationally misidentified as being hypermassive
because Swift is not able to slew in time to see the result of the energy
injection from the central engine.

###  5.3 Thermal evolution

The thermal evolution and history of a newly born neutron star has important
consequences on several aspects of its evolution, including, for example,
neutrino emission, the role of viscosity in dynamical and secular
instabilities, the freezing out of the crust, and the transformation of the
core to neutron superfluidity and proton superconductivity.

At very early times the star cools through neutrino emission via modified Urca
cooling Shapiro and Teukolsky (1983). In general, the neutrinos do not free
stream, implying the cooling timescale is set by the time it takes for the
neutrinos to diffuse out of the star (e.g., Rosswog and Liebendörfer, 2003).
It is these neutrinos that affect the electron fraction of the ejecta, and
hence alter the colour of the kilonova (see Sec. 4.2). An approximate cooling
timescale is Paschalidis  _et al._ (2012)

$\displaystyle t_{\rm cool}\approx
1200\left(\frac{M}{2.8\,\mathrm{M_{\odot}}}\right)\left(\frac{R}{15\,\mathrm{km}}\right)\left(\frac{E_{\nu}}{10\,\mathrm{MeV}}\right)^{2}{\rm
ms},$ (4)

where $E_{\nu}=10\,\mathrm{MeV}$ is the root-mean-square value of the neutrino
energy found in simulations Rosswog and Liebendörfer (2003). In general, the
cooling timescale for modified URCA cooling has a strong temperature
dependence, with $T\propto t^{1/6}$.

During the cooling of the remnant, the first major structural change is the
formation of the crystalline lattice crust. This typically takes
$\mathcal{O}({\rm days})$ for the base of the crust to begin forming, with a
complete crust not forming for up to a year following birth (e.g., Krüger  _et
al._ , 2015, and references therein). In the stellar core, neutron
superfluidity and proton superconductivity are expected at temperatures
$\lesssim(5-9)\times 10^{8}\,\mathrm{K}$ and $\lesssim(2-3)\times
10^{8}\,\mathrm{K}$, respectively (see Page  _et al._ , 2011; Shternin  _et
al._ , 2011, and references therein). While the development for superfluidity
and superconductivity are certainly relevant for late-time evolution of
neutron stars, e.g., in understanding glitches in young radio pulsars, it is
not clear what observational impact it has on the evolution of post-merger
remnants.

###  5.4 Dynamical evolution

Nascent neutron stars born in binary neutron star mergers are differentially
rotating with large poloidal magnetic fields in the range of $\sim 10^{14-16}$
G. This differential rotation winds up a toroidal component of the magnetic
field roughly symmetrical to the rotation axis. One expects this toroidal
component $B_{t}$ to dominate over the poloidal component $B_{p}$ deforming
the neutron star into a prolate ellipsoid (e.g., Cutler, 2002; Lander and
Jones, 2018). Misalignment between the rotation and magnetic axes result in
precession (e.g., Dall’Osso  _et al._ , 2018; Lander and Jones, 2018).
Internal dissipative processes drive the magnetic axis towards orthogonality
with the spin axis (Mestel and Takhar, 1972; Jones, 1976; Cutler, 2002),
minimising the energy of the system and also making the system an optimal
emitter of gravitational waves (e.g., Cutler, 2002; Lasky and Glampedakis,
2016; Dall’Osso  _et al._ , 2018). This instability is known as the “spin-flip
instability”.

In reality, the evolution of the misalignment angle $\chi$ between the star’s
magnetic and rotation axes is significantly more complicated. In general, for
a star deformed by a dominantly toroidal field, viscous dissipation increases
$\chi$. Conversely, $\chi$ decreases for a star deformed predominantly through
a poloidal field. Early efforts to model the evolution of $\chi$ focused on
determining the effects of this viscous dissipation. However, the evolution of
$\chi$ is in fact a complex mixture of dissipation, neutrino cooling, and the
effects of the external torque from spin-down which also acts to decrease
$\chi$ (Lander and Jones, 2018).

The coupling between the spin-down and viscous effects proves to be critical
in large parts of the neutron star parameter space and determining whether
$\chi\rightarrow 0^{\circ}$ (i.e., an aligned rotator that will not emit
gravitational waves) or $\chi\rightarrow 90^{\circ}$ (i.e., an orthogonal
rotator, an optimal emitter of gravitational waves) (Lander and Jones, 2018).
This coupling is ignored in Dall’Osso  _et al._ (2018) by assuming that the
timescale for spin flip is much faster than the spin down. This may not be
true, and depends on the dissipation timescale, which is inversely
proportional to the size of the deformation (Ipser and Lindblom, 1991). The
dissipation timescale is closely related to the spin-flip time-scale, an
approximation of which is (e.g., Lasky and Glampedakis, 2016),

$\tau_{\text{sf}}\approx
8\,\mathrm{s}\left(\frac{R}{15\,\mathrm{km}}\right)^{2}\left(\frac{\left<\rho\right>}{10^{15}\,\mathrm{gcm^{-3}}}\right)\left(\frac{\epsilon_{B}}{10^{-3}}\right)^{-1}\left(\frac{\epsilon_{\Omega}}{0.3}\right)^{-2}\left(\frac{\xi}{10^{30}\,\mathrm{g\,cm^{-1}\,s^{-1}}}\right)^{-1}$
(5)

Here, $\left<\rho\right>$ is the volume-averaged density, $R$ is the neutron
star radius, $\epsilon_{B}$ is the magnetic-field induced ellipticity,
$\epsilon_{\Omega}$ is the rotational ellipticity, and $\xi$ is the bulk
viscosity coefficient. For these parameters, it is difficult to imagine the
spin-flip timescale being much faster than the spin-down timescale especially
considering there are short gamma-ray bursts with observations on
significantly shorter timescales where the emission is potentially derived
from spin-down. However, this timescale is extremely sensitive to the
temperature $\xi\sim T^{6}$ and therefore the cooling history of the neutron
star, which is not well understood.

The evolution of $\chi$ is critical for developing an understanding what
mechanism is responsible for radiating away the rotational energy of the long-
lived neutron star (e.g., Margalit  _et al._ , 2018; Lander and Jones, 2020).
In the context of magnetars born in long gamma-ray bursts, Margalit  _et al._
(2018) showed that an aligned rotator will almost exclusively power a gamma-
ray burst, while a misaligned rotator will deposit some fraction of its
rotational energy onto the supernova. Similar behaviour is to be expected in a
long-lived neutron star born in a binary neutron star merger. In the first
several seconds while $\chi\rightarrow 0^{\circ}$ the rotational energy will
be lost in powering the gamma-ray burst (Lander and Jones, 2020). However once
the system begins to orthogonalise i.e., $\chi\rightarrow 90^{\circ}$ the
energy will be lost in gravitational waves and in powering the kilonova.

For neutron stars born in long gamma-ray bursts, the evolution of the braking
index through coupling of the spin and magnetic axes has potentially been
measured directly (Şasmaz Mus  _et al._ , 2019; Çıkıntoğlu  _et al._ , 2020).
Such a model has not yet been fit to short gamma-ray burst observations.
However, measurements of braking indices $n\lesssim 3$ of long-lived neutron
stars from short gamma-ray bursts hint towards the evolution of $\chi$ for a
number of objects (Lasky  _et al._ , 2017; Sarin  _et al._ , 2020b), although
it is worth noting there are other ways to explain such measurements (Lasky
_et al._ , 2017). If $\chi$ is evolving, it will be important to understand
the timescales and the long-term evolution. This will provide important clues
into the evolution of these objects into ordinary magnetars we see in our
galaxy today. While more immediately, it will allow for more informed
inferences about the radiation mechanisms that are tapping into the large
rotational-energy reservoir of long-lived neutron stars.

###  5.5 Electromagnetic observations

In Sec 5.1 we discussed whether a neutron star could launch an ultra-
relativistic jet and produce a short gamma-ray burst. Here, we discuss the
broader electromagnetic imprints of a long-lived neutron star remnant.

The longest numerical simulations of binary neutron stars and their remnants
last approximately 100 ms post merger (e.g., De Pietri  _et al._ , 2018;
Ciolfi  _et al._ , 2019; De Pietri  _et al._ , 2020), significantly shorter
than the time where they could be used to provide insight into the
electromagnetic signature of a long-lived neutron star. Such insights
therefore rely predominantly on analytical and semi-analytic models (e.g.,
Zhang and Mészáros, 2001; Dall’Osso  _et al._ , 2011; Metzger and Piro, 2014).

The diverse predictions of electromagnetic signatures from long-lived post-
merger remnants can all be primarily attributed to the large reservoir of
rotational energy of the long-lived neutron star. Unlike hypermassive neutron
stars which can trap a significant amount of their spin-down energy as they
collapse into black holes (e.g., Metzger, 2017; Shibata  _et al._ , 2019),
supramassive and infinitely stable neutron stars will radiate a large fraction
of this energy away with several electromagnetic consequences.

One of the signatures of spin-down energy is on the kilonova itself (e.g., Yu
_et al._ , 2013). The injection of spin-down energy is either via a Poynting
flux from a collimated jet (Bucciantini  _et al._ , 2012) or photons generated
from dissipation of a magnetar wind (Thompson  _et al._ , 2004; Zhang, 2013).
This increased energy creates a distinct impact on the kilonova lightcurve,
altering the peak time and duration, while drastically increasing the
luminosity of the kilonova compared to a kilonova only powered by radioactive
decay. This latter characteristic implies the colour of the kilonova becomes
bluer, simply due to higher luminosity translating to a higher effective
temperature for a similar photospheric radius. Such effects are identical to
the differences seen between ordinary and magnetar-driven supernovae (Kasen
and Bildsten, 2010; Nicholl  _et al._ , 2020; Margalit  _et al._ , 2018). As
mentioned in Sec. 2.1, this is one of the interpretations of the kilonova
following GW170817 (Yu  _et al._ , 2018). The recent identification of an
exceptionally bright kilonova following short GRB200522A (Fong  _et al._ ,
2020) provides further tantalising evidence for this magnetar-driven kilonova
scenario. An additional imprint of this increased ejecta energy will be seen
on the synchrotron radio signal generated from the interaction of the ejecta
with the interstellar medium (Hotokezaka and Piran, 2015; Horesh  _et al._ ,
2016; Ricci  _et al._ , 2020; Schroeder  _et al._ , 2020). Such constraints
have been used to infer the fate of several gamma-ray bursts, suggesting that
$<50\%$ of binary neutron star mergers make a long-lived neutron star
(Schroeder  _et al._ , 2020). However, we note that this may be overly
constraining as they ignore rotational energy lost through other channels.

Even disregarding the impact of the spin-down luminosity, a long-lived remnant
will show a distinctive characteristic feature on the kilonova. As discussed
in Secs. 4.2 and 5.3, the nascent neutron star cools rapidly through neutrino
cooling. This additional neutrino flux dramatically increases the electron
fraction of the ejecta and suppresses the production of lanthanides (e.g.,
Lippuner  _et al._ , 2017; Metzger, 2017; Kawaguchi  _et al._ , 2020). Such
ejecta would naturally be less opaque than more lanthanide-rich ejecta making
the resulting kilonova dominantly “blue” regardless of how much rotational
energy is deposited into the kilonova ejecta.

Potentially some of the best electromagnetic observations of nascent neutron
stars is the evidence of energy injection in the x-ray afterglows of gamma-ray
bursts. A large fraction of x-ray afterglows have plateaus followed by sharp
drops in luminosity (e.g., Rowlinson  _et al._ , 2013), which are challenging
to explain with the canonical model for afterglows that model the interaction
of a relativistic jet with the surrounding environment. The spin-down energy
from a nascent neutron star can naturally explain both these features (e.g.,
Zhang and Mészáros, 2001; Rowlinson  _et al._ , 2013). However, the exact
mechanism that extracts this rotational energy is still uncertain. Either the
spin-down energy is extracted directly from the remnant itself (e.g., Zhang
and Mészáros, 2001; Rowlinson  _et al._ , 2013; Lasky  _et al._ , 2017; Strang
and Melatos, 2019) with a constant efficiency, or indirectly through energy
injection at the afterglow shock interface (e.g., Cohen and Piran, 1999;
Dall’Osso  _et al._ , 2011; Sarin  _et al._ , 2020b). Systematic model
selection suggests the latter model, with a generalised braking index better
explaining the observations (Sarin  _et al._ , 2020b), at least for the small
subset (eight) of short and long gamma-ray bursts analysed. Such a model also
self-consistently explains x-ray flares that are seen at the onset of the
plateau phase, providing additional support for this model.

An alternative way to generate x-ray emission with a magnetar is through the
interaction of a magnetar wind with the merger ejecta (e.g., Yu  _et al._ ,
2013; Metzger and Piro, 2014; Siegel and Ciolfi, 2016; Strang and Melatos,
2019; Strang  _et al._ , 2021). Assuming the magnetar is completely enshrouded
in the merger ejecta, the magnetar wind energy will be dissipated via shocks
or magnetic re-connection creating a hot nebula behind the ejecta. This hot
nebula will be comprised of photons and electron/positron pairs analogous to a
pulsar wind nebula (e.g., Metzger and Piro, 2014). Initially, most of the
spin-down energy will be lost in expanding the nebula and ejecta along with
it. This expansion will eventually reduce the optical depth allowing photons
at various wavelengths to diffuse out once the diffusion timescale becomes
shorter than the expansion timescale (Metzger and Piro, 2014). While such a
model explains several observational features of x-ray afterglows, it is
problematic as this model cannot explain x-ray emission at early-times
($T\lesssim 10\,\mathrm{hr}$) as the ejecta is still optically thick. In this
model, early-time x-ray observations must be from synchrotron emission from
the interaction of the relativistic jet with the surrounding environment, with
the neutron star component becoming dominant after the ejecta becomes
optically thin. Such a constraint is, however, problematic, given sharp drops
in x-ray plateaus are seen as early as 100 s after the prompt (e.g., Sarin
_et al._ , 2020a). There are two ways to reconcile this issue. First, the
timescale for the ejecta to become optically thin need not be as long as 10
hr, as this timescale is dependent on the opacity of $r$-process elements,
which are far from certain (e.g., Tanaka  _et al._ , 2020; Barnes  _et al._ ,
2020). Second, and perhaps more importantly, the simple assumption that a
magnetar is wholly enshrouded in the merger ejecta may not be valid, with
potential holes due to piercing by the ultra-relativistic jet that produced
the gamma-ray burst or by Rayleigh-Taylor instabilities in the ejecta itself
(e.g., Strang and Melatos, 2019; Strang  _et al._ , 2021).

An important question to consider with the spin-down energy is whether the
radiation is emitted isotropically or beamed in a particular direction,
perhaps along with the jet that produces the gamma-ray burst prompt emission
itself. The answer to this question could shed critical insight into the
emission mechanism. Observations and interpretations of CDF-S XT2 (Xue  _et
al._ , 2019) as a magnetar seen off-axis suggests the magnetar spin-down
energy is emitted isotropically. However, isotropic emission potentially leads
to violations of the total energy budget for certain gamma-ray bursts believed
to be powered by nascent neutron stars (e.g., Beniamini and Mochkovitch,
2017). By contrast, if the spin-down energy is collimated, it must be through
a mechanism that can be sustained for long timescales to explain the late-time
x-ray observations. One explanation for the collimation originates from the
interaction of the magnetar wind and the surrounding environment.
Magnetohydrodynamic simulations suggest that the deceleration of the magnetar
wind due to the dense environment of the merger ejecta may collimate the wind
into a jet for sufficiently high spin-down luminosity (Bucciantini  _et al._ ,
2012). However, such a mechanism seems unfeasible, as the jet will become
unstable and susceptible to magnetohydrodynamic instabilities once the spin-
down luminosity of the long-lived remnant drops (Porth  _et al._ , 2013;
Metzger and Piro, 2014).

Ultimately, testing detailed models with better observations will allow us to
determine what physical processes are relevant. Beyond these more immediate
electromagnetic observations, long-lived neutron stars may also be responsible
for fast radio bursts, to which we now turn our attention.

###  5.6 Fast radio bursts

Fast radio bursts are millisecond duration pulses of coherent radio emission.
In general they come from cosmological distances, with dispersion measures
significantly larger than galactic values. We refer the reader to Refs. Cordes
and Chatterjee (2019); Zhang (2020a) for recent reviews on fast radio bursts.

Ever since the original discovery (Lorimer  _et al._ , 2007), subsequent
identification of fast-radio burst repeaters (Spitler  _et al._ , 2014, 2016)
and localisation of its host galaxy (Chatterjee  _et al._ , 2017; Tendulkar
_et al._ , 2017), nascent magnetars have been invoked to explain these
enigmatic astrophysical phenomena (e.g., Popov and Postnov, 2013; Lyubarsky,
2014; Kulkarni  _et al._ , 2014; Katz, 2016; Beloborodov, 2017; Metzger  _et
al._ , 2017, 2019; Katz, 2018; Lu and Kumar, 2018; Margalit  _et al._ , 2019;
Lu  _et al._ , 2020). While it is possible that there are two distinct
populations of fast radio bursts, characterised as repeaters and isolated
bursts (e.g., Falcke and Rezzolla, 2014; Zhang, 2017, 2020b), a non-repeating
population may be disfavoured on rate-based arguments (e.g., Nicholl  _et al._
, 2017b; Ravi, 2019). We note there are numerous other alternative progenitor
models, in fact, as recently as 2019, theoretical models outnumbered the
number of events themselves (Platts  _et al._ , 2019). However, this all
changed with the recent watershed discovery of a fast radio burst from a
galactic magnetar, SGR1935+2154, which provided the smoking-gun observation
for the magnetar origin for fast radio bursts.

SGR1935+2154 is a galactic soft gamma repeater first identified by Swift as a
potential gamma-ray burst candidate (Stamatikos  _et al._ , 2014), it is
associated with a supernova remnant at a distance $d\approx 6$ kpc (Gaensler,
2014; Zhou  _et al._ , 2020), with an estimated surface magnetic field
$B\approx 10^{14}$ G and age (based on the supernova remnant association)
$\gtrsim 16$ kyr (Zhou  _et al._ , 2020). On April 28, 2020, a millisecond
duration radio pulse was independently detected from this source by CHIME (The
Chime/Frb Collaboration  _et al._ , 2020) and STARE2 (Bochenek  _et al._ ,
2020) in coincidence with a bright X-ray burst (Zhang  _et al._ , 2020).
Various analyses confirmed the analogous nature of the coherent radio emission
with cosmological fast radio bursts albeit with significantly lower energy (at
least $\sim 25$ times) than typical cosmological fast radio bursts. This
observation provided unequivocal evidence that magnetars are the progenitors
of at least some fast radio bursts (e.g., Margalit  _et al._ , 2020), with the
discrepancy in energies attributed to the old age of this magnetar, weaker
magnetic field, and slower rotation (e.g., Lu  _et al._ , 2020; Beloborodov,
2020).

Observations of SGR1935+2154 provide evidence that young, rapidly rotating
magnetars produce some fast radio bursts. However, it is unclear whether these
magnetars are ones born in core-collapse supernovae or in binary neutron star
mergers. The first repeating fast radio burst FRB121102 was localised to a low
metallicity dwarf star-forming galaxy (Chatterjee  _et al._ , 2017; Tendulkar
_et al._ , 2017) and also associated with a persistent radio source (Marcote
_et al._ , 2017). The host galaxy properties and the persistent radio source
are best interpreted as the emission from a young magnetar embedded in the
expanding supernova ejecta (e.g., Omand  _et al._ , 2018; Margalit and
Metzger, 2018) pointing towards a superluminous supernovae/long gamma-ray
burst origin for FRB121102. The Ausrtralian Square Kiklometer Array Pathfinder
(ASKAP) has since localised another fast radio burst FRB180924 (Bannister  _et
al._ , 2019), finding the host galaxy properties and offsets to be comparable
to short gamma-ray bursts. FRB180924, therefore, provided the first possible
evidence for a binary neutron star remnant merger origin for fast radio
bursts. Since then, more fast radio bursts have been localised in host
galaxies strongly suggestive of binary neutron star merger remnant origins
(e.g., Ravi  _et al._ , 2019). Rate-based arguments suggest all fast radio
bursts are repeaters (e.g., Ravi, 2019), implicitly demanding long-lived
sources, i.e., long-lived neutron star remnants. Moreover, as supramassive
neutron stars have been shown to collapse on relatively short timescales
(e.g., Sarin  _et al._ , 2020a) where the kilonova ejecta is still optically
thick, implying fast radio bursts would not escape to an external observer
(e.g., Margalit  _et al._ , 2019). The timescale for this ejecta to become
optically thin ranges from weeks to months, implying any fast radio bursts
associated with binary neutron star mergers must be from infinitely stable
neutron stars (e.g., Margalit  _et al._ , 2019).

While the progenitor model of fast radio bursts can be confidently confirmed
as a nascent magnetar, the exact mechanism that generates the coherent radio
emission is unknown. The various models can be broadly divided into two
categories: those for which the emission is close to magnetar i.e., a pulsar-
like mechanism involving the magnetosphere (e.g., Katz, 2016; Lu  _et al._ ,
2020), and those for which the emission mechanism involves relativistic shocks
similar to gamma-ray bursts, where the coherent radio emission is generated
far from the magnetar (e.g., Beloborodov, 2017; Metzger  _et al._ , 2019).
Future observations of fast radio bursts will be able to shed light into which
mechanism is correct.

###  5.7 Gravitational-wave emission and detection

####  5.7.1 Gravitational-wave emission

Long-lived neutron stars also emit gravitational waves, but we have an
incomplete understanding of which mechanisms are relevant, how long they are
active for, or how much energy is emitted. As we discuss later, this
uncertainty weakens our ability to detect gravitational waves from a long-
lived neutron star while also weakening our ability to infer the fate from
indirect energetic constraints from electromagnetic observations (e.g.,
Shibata  _et al._ , 2019).

Three main instabilities are relevant for producing gravitational waves in
long-lived post-merger remnants. These are the spin-flip, bar-mode and r-mode
instabilities (e.g., Cutler, 2002; Lai and Shapiro, 1995; Shapiro, 2000;
Shibata  _et al._ , 2000; Andersson and Kokkotas, 2001; Andersson, 2003). Note
there are potentially other mechanisms, such as gravitational-wave emission
due to the formation of mountains from fall-back accretion (Sur and Haskell,
2020). Here we briefly discuss these gravitational-wave emission mechanisms
following on from the discussion in Sec 5.4. For a detailed review of
gravitational waves from neutron stars we refer the reader to (e.g., Lasky,
2015, and references therein).

The precessional spin-flip instability drives a nascent neutron star to become
an orthogonal rotator and therefore an optimal emitter of gravitational waves
(Cutler, 2002). Past studies (e.g., Lander and Jones, 2018) have shown that
once a system evolves to being near-aligned or near-orthogonal, it does not
further evolve. However, more recent work (e.g., Lander and Jones, 2020) has
shown that including the effect of a neutrino-driven wind can change the late-
time behaviour, with orthogonal rotators slowly decreasing the angle between
magnetic and rotation axes over hundreds of years to become more aligned.
Lander and Jones (2018) also showed that nascent neutron stars would become
orthogonal rotators for sufficiently large toroidal magnetic fields
$B_{t}\gtrsim 10^{14}$G. Such magnetic fields are naturally expected in binary
neutron star mergers (e.g., Giacomazzo and Perna, 2013; Kiuchi  _et al._ ,
2014; Ciolfi  _et al._ , 2019; Mösta  _et al._ , 2020), making it likely that
long-lived neutron stars born in binary neutron star mergers become orthogonal
rotators and stay that way for a long time.

As discussed in Sec 5.4, the timescale for orthogonalisation is far from
certain, but is a critical component for building gravitational-wave waveform
models for the spin-flip instability. An aligned rotator does not emit
gravitational waves. As the angle between the magnetic and rotation axes
increases, the star will emit gravitational waves dominantly at the spin
frequency, or twice the spin frequency, with a host of other potential values
that depend on the precession timescale (e.g., Jones and Andersson, 2002;
Lasky and Melatos, 2013). Building gravitational-wave waveform models for this
evolution requires careful modelling of the orthogonlisation timescale,
including the evolution of the misalignment angle between the rotation and
magnetic axes. If the orthogonlisation timescale is short compared to the
overall gravitational-wave emission timescale, then the uncertainty associated
with the evolution of the misalignment angle can be largely ignored. Such an
assumption would be detrimental to detection prospects if the
orthogonalisation timescale is comparable to the timescale for gravitational-
wave emission.

Assuming orthogonalisation has already taken place, Sarin  _et al._ (2018)
built a model for the gravitational-wave signature of long-lived post-merger
remnants. Unfortunately, detection prospects are not good with second-
generation gravitational-wave interferometers; for detection at greater than a
few Mpc, the amount of emitted gravitational-wave energy must exceed the
rotational energy budget of the system.

The critical quantities determining the energy produced in gravitational waves
for gravitational-wave dominated spin down (i.e., the most optimistic scenario
cf. when the spin down is dominated by electromagnetic torques) are the
initial spin period and ellipticity. The strong magnetic fields expected in a
binary neutron star merger deform the nascent neutron star, with the size of
this deformation dependent on the magnetic field and the internal geometry of
the field. For simple stellar models, the ellipticity can be approximated as
Cutler (2002)

$\epsilon_{B}\approx 10^{-6}\left(\frac{\left\langle
B_{\mathrm{t}}\right\rangle}{10^{15}\mathrm{G}}\right)^{2},$ (6)

where $\langle B_{\mathrm{t}}\rangle$ is the volume-averaged toroidal magnetic
field strength. Observational inferences based on the collapse times of
supramassive neutron stars from short gamma-ray bursts suggests
$\epsilon_{B}\sim 10^{-3}$ (Gao  _et al._ , 2016). Such large ellipticities
require $\langle B_{\mathrm{t}}\rangle\sim 10^{16}$–$10^{17}\,\mathrm{G}$,
i.e., toroidal fields that are 1-2 orders of magnitude stronger than inferred
values of the poloidal field. It is an open question whether such magnetic
fields are dynamically stable (e.g., Braithwaite, 2009; Lasky  _et al._ ,
2011; Ciolfi and Rezzolla, 2012; Akgün  _et al._ , 2013; Herbrik and Kokkotas,
2017; Sur and Haskell, 2020).

Another mechanism to generate gravitational waves in a newly-born neutron is
through unstable $f$ modes, also known as the bar-mode instability. This comes
in two varieties: dynamical and secular. To activate the dynamical
instability, the ratio of a star’s rotational kinetic energy $T$ to
gravitational binding energy $W$ must be $T/|W|\geq 0.27$ (Lai and Shapiro,
1995; Shapiro, 2000; Corsi and Mészáros, 2009), while the secular instability
is active for $T/|W|\geq 0.14$. The former instability is only possible for
realistic equations of state in differentially rotating neutron stars (e.g.,
Friedman and Schutz, 1978; Corsi and Mészáros, 2009; Ravi and Lasky, 2014;
Lasky and Glampedakis, 2016) with the growth of the dynamical instability
halted when the magnetic field damps differential rotation, governed by the
Alfvén timescale of $\mathcal{O}(1)$s (Shapiro, 2000). Therefore, the
dynamical bar mode is only expected to be active in the first seconds of a
neutron star’s life.

More relevant for long-term gravitational-wave emission is the secular
instability. It is worth noting that several equations of state that support
masses $\gtrsim 2\,\mathrm{M_{\odot}}$ do not have rigidly-rotating
equilibrium solutions with $T/|W|\geq 0.14$, implying the secular bar-mode
instability cannot be active for those equations of state Ravi and Lasky
(2014). However, if active, the secular bar mode grows to large non-linear
amplitudes on timescales of $\mathcal{O}(10$–$100)$ s (Doneva  _et al._ ,
2015). For the instability to grow, the nascent neutron star has to cool to a
temperature below $\sim 10^{10}\,\mathrm{K}$, such that bulk viscosity does
not suppress the instability. As we have discussed in Sec 5.3, a newly-born
neutron star is cooled to this temperature a few seconds after the merger,
making the secular-bar mode a potentially-dominant mechanism for
gravitational-wave emission Corsi and Mészáros (2009).

Another oscillation mode that will likely lead to the gravitational-wave
radiation are the $r$ modes; low-frequency toroidal oscillations for which the
Coriolis force is the restoring force. These oscillations are retrograde in
the co-moving frame and prograde in the inertial frame,(Andersson and
Kokkotas, 2001; Andersson, 2003), making them always unstable to the
Chandrasekhar-Friedman-Schutz (CFS) instability (Chandrasekhar, 1970; Friedman
and Schutz, 1978). Whether this instability is active is dependent on a
delicate balance between gravitational-wave radiation which drives up the size
of the instability, and viscous processes that dampen the oscillations. For
nascent neutron stars, the dominant viscous force is bulk viscosity, which in
turn depends sensitively on neutron star microphysics, cooling history, and
rotation (e.g., Andersson, 2003; Lasky, 2015).

Recently, an interpretation of the x-ray afterglow of GRB090510 suggests that
the observations support gravitational-wave losses due to $r$-mode
oscillations through the measurement of the braking index $n=7$ (Lin and Lu,
2019). This interpretation is contentious, most notably because a braking
index of $n=7$ is not necessarily a reliable indicator of spin-down through
$r$-mode gravitational-wave emission (e.g., Alford and Schwenzer, 2014, 2015)

Long-lived neutron stars could also emit gravitational waves from mountains
formed through fall-back accretion (e.g., Piro and Thrane, 2012; Melatos and
Priymak, 2014; Sur and Haskell, 2020). However, given the relatively little
amount of material ejecta in a binary neutron star merger, this mechanism is
likely only relevant for neutron stars born in supernovae and long gamma-ray
bursts.

All the different mechanisms described above suggest it is quite likely that
long-lived neutron stars spin down, at least in some part, through
gravitational wave emission. This is to an extent observationally verified by
the observed collapse time distribution (Sarin  _et al._ , 2020a), which
suggests $\sim 70\%$ of these neutron stars are spinning down predominantly
through gravitational-wave radiation. However, knowing which mechanism is
active and for how long can only be verified through the direct detection of
gravitational waves from these objects, something to which we now turn.

####  5.7.2 Gravitational-wave detection methods

Searching for gravitational waves from long-lived neutron stars suffers from
some of the same problems as traditional searches for continuous gravitational
waves (see Riles, 2013, and references therein). The long-duration signals
expected implies traditional matched-filtering approaches are computationally
unfeasible (e.g., Brady  _et al._ , 1998) requiring semi-coherent methods that
are less sensitive than fully coherent methods. Similarly, the uncertainty in
the gravitational-wave modelling necessitates the use of unmodelled searches,
which by design are not as sensitive as modelled searches.

The LIGO/Virgo Collaborations searched for gravitational waves from a
potential long-lived post-merger remnant from GW170817 on intermediate
($\lesssim 500$ s; Abbott  _et al._ , 2017d) and long ($\lesssim 8$ d; Abbott
_et al._ , 2019a) timescales (the latter timescale being set by the length of
data available following the merger). No viable candidate was found, however
upper limits on the gravitational-wave strain and energy were derived.

A number of complementary pipelines were used in the LIGO/Virgo searches. The
robustly-named Stochastic Transient Analysis Multi-detector Pipeline (STAMP;
Thrane  _et al._ , 2011) looks for tracks of excess power in cross-correlated
data from the two detectors using seedless clustering algorithms Thrane and
Coughlin (2013, 2015). Coherent Wave Burst (cWB; Klimenko  _et al._ , 2016)
was used with a similar setup as for the short-duration search; see Sec. 4.3.
The Hidden Markov Tracking method using the Viterbi algorithm Suvorova  _et
al._ (2016); Sun and Melatos (2019) searches for quasichromatic signals with
unknown frequency evolution and stochastic timing noise. The Adaptive
Transient Hough Krishnan  _et al._ (2004); Oliver  _et al._ (2019b) assumes
the signal’s frequency evolution can be modelled as a power law, with the
amplitude and frequency described by the generalised millisecond magnetar
model Lasky  _et al._ (2017); Sarin  _et al._ (2018). Finally, the Generalized
FrequencyHough algorithm Antonucci  _et al._ (2008); Astone  _et al._ (2014);
Miller  _et al._ (2018) is a pattern-recognition technique that maps time-
frequency points to lines in frequency-spin down space.

Four of the aforementioned algorithms are unmodelled searches, and one is
modelled. Upper limits for the unmodelled searches were derived using two
theoretically-motivated signal models; the generalised millisecond magnetar
model that describes a spinning-down neutron star with arbitrary braking index
Lasky  _et al._ (2017); Sarin  _et al._ (2018), and a waveform model Corsi and
Mészáros (2009) where the star’s spin evolution is dictated by the
gravitational-wave driven secular Chandrasekhar-Friedman-Schutz instability
Chandrasekhar (1970); Friedman and Schutz (1978).

The derived upper limits from the searches do not necessarily bode well for
future detections of gravitational waves from long-lived neutron star
remnants. Throughout the entire explored parameter space, the distance at
which a source could have been observed from the remnant of GW170817 was, at
best, just $\sim 1\,\mathrm{Mpc}$, cf. the actual distance of $\sim
40\,\mathrm{Mpc}$.

Although aLIGO/Virgo’s sensitivity is predicted to improve by a factor of a
few over the sensitivity at the time GW170817 was detected, the prospect for
detection with second-generation telescopes is still grim. From an
observational perspective, aLIGO/Virgo could eventually be sensitive to
mergers at best at $\sim 10\,\mathrm{Mpc}$. Theoretical estimates are
consistent with this (e.g., see Refs. Corsi and Mészáros (2009); Fan  _et al._
(2013); Dall’Osso  _et al._ (2015); Doneva  _et al._ (2015), although see
Refs. Dall’Osso  _et al._ (2018); Lasky and Glampedakis (2016); Sarin  _et
al._ (2018) for more pessimistic estimates). However, in practice, the most
optimistic of these estimates may require nonphysical quantities of
gravitational-wave energy, for example in excess of the total rotational
energy budget of the system Sarin  _et al._ (2018). A successful detection of
gravitational waves from a long-lived post-merger remnant may therefore have
to wait until A+ Miller  _et al._ (2015), or even third-generation detectors
such as Einstein Telescope Punturo  _et al._ (2010) or Cosmic Explorer Abbott
_et al._ (2017e).

##  6 Conclusions and Outlook

The era of gravitational-wave and electromagnetic multi-messenger astronomy
began spectacularly with GW170817 Abbott  _et al._ (2017a, b, c), a watershed
event that provided invaluable insights across several domains of high energy
astrophysics. Observations of GW190425 transformed our understanding of the
formation of binary neutron stars, highlighting the flaw in the simple
assumption that the mass distribution of cosmological binary neutron stars
follows the same mass distribution as ones observed locally in our galaxy
(e.g., Abbott  _et al._ , 2020). The future of electromagnetic and
gravitational-wave multi-messenger astronomy looks bright, with upgrades to
both gravitational-wave and electromagnetic detectors in progress (e.g.,
Abbott  _et al._ , 2016; Bellm  _et al._ , 2019) and proposed for the future
(e.g., Ivezić  _et al._ , 2019; Punturo  _et al._ , 2010; Abbott  _et al._ ,
2017e; Ackley  _et al._ , 2020).

Notwithstanding our limited understanding of the maximum mass of neutron stars
and the binary neutron star mass distribution, we expect potentially up to
$\sim 80\%$ of neutron star mergers to produce some form of neutron-star
remnant Margalit and Metzger (2019). However, at least for the near future,
the smoking gun observation of gravitational waves from short or long-lived
remnants is unlikely. As such, the fate of the remnant will need to be
inferred through the suite of electromagnetic observations.

In the near future, electromagnetic observations, in particular, the early
x-ray afterglow and kilonova observations provide the best probe into the
nature of the remnant. Late time radio limits on the energy of the ejecta will
also provide clues into the nature of the remnant. The latter analyses will
require a good understanding of the different channels of energy emission, in
particular gravitational waves. If future research concludes that long-lived
remnants cannot launch ultra-relativistic jets capable of producing short
gamma-ray bursts. This would imply that short gamma-ray bursts are a biased
and relatively small fraction of binary neutron star mergers, something that
will become telling as the rate of binary neutron star mergers becomes better
constrained. Furthermore, any discrepancy in binary neutron star merger and
gamma-ray burst rates must therefore be explained by neutron star-black hole
mergers, a statement which has significant implications on the properties of
such binaries.

Although x-ray afterglow and kilonova observations are the best
electromagnetic channels for determining the fate of the post-merger remnant,
they alone may also not be definitive without further development and testing
of models, requiring input from theorists, simulations and observers. In
particular, the idea of using the colour of kilonova to probe the lifetime of
the remnant is fraught with difficulties such as viewing angle dependence
(e.g., Darbha and Kasen, 2020), and the uncertain impact of the jet-ejecta
interaction (e.g., Nativi  _et al._ , 2020). These uncertainties, combined
with our incomplete knowledge of nuclear reactions and opacities of r-process
elements, implies significant systematic uncertainties for inferring the
properties of kilonovae from observations (e.g., Zhu  _et al._ , 2020).

Inferring the fate of the remnant from the early-time x-ray afterglow
observations may be more reliable, particularly if the emission from the long-
lived neutron star is isotropic (as suggested by the observations of CDF-S XT2
(e.g., Xue  _et al._ , 2019)). This is different from the emission expected
from the interaction of the jet with the interstellar medium (the physics
known to be responsible for the afterglow in normal circumstances), which is
strongly affected by relativistic beaming (e.g., Totani and Panaitescu, 2002;
Granot  _et al._ , 2002). However, models for the emission from nascent
neutron stars are in their infancy, with significant development required such
that they accurately reflect all the critical physics (e.g., Metzger and Piro,
2014; Lasky  _et al._ , 2017; Şasmaz Mus  _et al._ , 2019; Strang and Melatos,
2019; Sarin  _et al._ , 2020b). Jet structure could also explain such
observations without requiring a neutron star remnant; in such a scenario,
systematic model selection may provide the answer (see e.g., Sarin  _et al._ ,
2019).

The observations of an x-ray plateau with a sharp drop in luminosity following
a binary neutron star merger may provide the most definitive electromagnetic
evidence for the fate of a binary neutron star merger. Here again there may be
different emission mechanisms such as a reverse shock formed from the
interaction of a relativistic jet with the surrounding interstellar medium
(e.g., van Eerten, 2018; Lamb and Kobayashi, 2019), or radial stratification
of the jet such that it is refreshed at late times (e.g., Lamb  _et al._ ,
2020). High-latitude emission (e.g., Oganesyan  _et al._ , 2020; Ascenzi  _et
al._ , 2020), or fall back accretion Desai  _et al._ (2019) that may also
explain the sharp drop in luminosity without requiring a supramassive neutron
star. Fortunately, if the gamma-ray burst is observed off-axis, several of
these scenarios become less likely to be the cause (e.g., van Eerten, 2018).
However, to ensure this sharp drop in luminosity is observed, the
electromagnetic counterpart of a binary neutron star merger must be identified
quickly, on timescales as short as $100$s, which will be difficult in the near
future. In light of these theoretical and observational issues, it is the
combination of various electromagnetic phenomena and richer data confronted
with more detailed models that can shed light into the nature of the remnant.

In this review, we have discussed the fate of binary neutron star mergers
GW170817 and GW190425 and potential other neutron star mergers seen as short
gamma-ray bursts. We have reviewed all possible outcomes of a binary neutron
star merger from the prompt collapse into a black hole to the formation of an
infinitely stable neutron star, discussing their observational signature,
evolution, and prospects for gravitational-wave detection. As observations of
binary neutron star mergers grow in number, understanding the fate of the
remnant will become increasingly more important due to its far-reaching
implications, such as on the nuclear equation of state, gamma-ray bursts,
kilonovae, fast radio bursts, and beyond. At least in the near future, the
lack of smoking-gun gravitational-wave observations means the fate must be
inferred from electromagnetic observations. The promise of richer and more
frequent observations confronted with better models ensures this will be an
exciting endeavour.

###### Acknowledgements.

We are grateful to Carl Knox for the creation of Fig. 1. We also thank Andrew
Melatos and Toni Font for their helpful comments on the manuscript. This work
is supported through Australian Research Council (ARC) Centre of Excellence
CE170100004, ARC Future Fellowship FT160100112 and ARC Discovery Project
DP180103155. N.S is supported by an Australian Government Research Training
Program (RTP) Scholarship.

## References

  * Abbott  _et al._ (2017a) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , Phys. Rev. Lett. 119, 161101 (2017a).
  * Abbott  _et al._ (2017b) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , Astrophys. J. 848, L12 (2017b).
  * Abbott  _et al._ (2017c) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , Astrophys. J. 848, L13 (2017c).
  * Aasi  _et al._ (2015) J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, _et al._ , Classical and Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, _et al._ , Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ (2017d) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , ApJ 851, L16 (2017d).
  * Abbott  _et al._ (2019a) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , ApJ 875, 160 (2019a), arXiv:1810.02581 [gr-qc] .
  * von Kienlin  _et al._ (2017) A. von Kienlin, C. Meegan, and A. Goldstein, GRB Coordinates Network, Circular Service 21520 (2017).
  * Meegan  _et al._ (2009) C. Meegan, G. Lichti, P. N. Bhat, E. Bissaldi, M. S. Briggs, _et al._ , ApJ 702, 791 (2009), arXiv:0908.0450 [astro-ph.IM] .
  * Savchenko  _et al._ (2017) V. Savchenko, C. Ferrigno, E. Kuulkers, A. Bazzano, E. Bozzo, _et al._ , ApJ 848, L15 (2017), arXiv:1710.05449 [astro-ph.HE] .
  * Arcavi  _et al._ (2017) I. Arcavi, G. Hosseinzadeh, D. Howell, C. McCully, D. Poznanski, _et al._ , Nature 551, 64 (2017).
  * Coulter  _et al._ (2017) D. A. Coulter, R. J. Foley, C. D. Kilpatrick, M. R. Drout, A. L. Piro, _et al._ , Science 358, 1556 (2017), arXiv:1710.05452 [astro-ph.HE] .
  * Lipunov  _et al._ (2017) V. M. Lipunov, E. Gorbovskoy, V. G. Kornilov, N. . Tyurina, P. Balanutsa, A. Kuznetsov, D. Vlasenko, D. Kuvshinov, I. Gorbunov, D. A. H. Buckley, A. V. Krylov, R. Podesta, C. Lopez, F. Podesta, H. Levato, C. Saffe, C. Mallamachi, S. Potter, N. M. Budnev, O. Gress, Y. Ishmuhametova, V. Vladimirov, D. Zimnukhov, V. Yurkov, Y. Sergienko, A. Gabovich, R. Rebolo, M. Serra-Ricart, G. Israelyan, V. Chazov, X. Wang, A. Tlatov, and M. I. Panchenko, ApJ 850, L1 (2017), arXiv:1710.05461 [astro-ph.HE] .
  * Soares-Santos  _et al._ (2017) M. Soares-Santos, D. E. Holz, J. Annis, R. Chornock, K. Herner, _et al._ , ApJ 848, L16 (2017), arXiv:1710.05459 [astro-ph.HE] .
  * Tanvir  _et al._ (2017) N. R. Tanvir, A. J. Levan, C. González-Fernández, O. Korobkin, I. Mandel, _et al._ , ApJ 848, L27 (2017), arXiv:1710.05455 [astro-ph.HE] .